实验5 常微分方程初值问题和矩阵特征值的计算

实验5 常微分方程初值问题和矩阵特征值的计算

思考问题：欧拉折线法与泰勒级数法的关系是什么？如何把乘幂法变为反幂法来求解按模最小特征值？

答：欧拉折线法就是用一系列折线近似的代替曲线，欧拉折线法思想就是泰勒级数展开，其解得精度是差分步长的一阶精度。泰勒公式展开是一种高阶显式一步法，理论上它具有任意阶精度；

把乘幂法变为反幂法来求解按模最小特征值的方法为该求方阵A的逆矩阵A，乘幂法

11对应的模最大特征值取倒数为，变为模最小特征值，由此得到反幂法，可反复迭代：y(k)Ay(k1)，y(k)/y(k1)收敛于A的主特征值，y(k)为对应的特征向量。

yyxy25.1 取步长h0.2，用显示欧拉法求y(0)1在x0.6处y的近似值。

提示：循环3次，用3段折线来逼近函数y。

运行结果为：

1

246A391541636，用改进后的幂乘法求A得主特征值及其对应的特征向量。 5.2 已知矩阵

运行结果为：

程序代码

2

5.1

#include

#include

#define MAXSIZE 50

double f(double x,double y);

void main(void)

{

double a,b,h,x[MAXSIZE],y[MAXSIZE];

long i,n;

printf(\"\\nPlease input interval a,b:\");

scanf(\"%lf,%lf\

printf(\"\\nPlease input step length h:\");

scanf(\"%lf\

3

n=(long)((b-a)/h);

x[0]=a;

printf(\"\\nPlease input the start point of x[0]=%lf ordinate y[0]:\

scanf(\"%lf\

for(i=0;i{

x[i+1]=x[i]+h;

y[i+1]=y[i]+h * f(x[i],y[i]);

}

printf(\"\\nThe calcuate result is:\");

for(i=0;i<=n;i++)

printf(\"\\nx[%ld]=%lf,y[%ld]=%lf\

getch();

4

}

double f(double x,double y)

{

return(-y-x*y*y); /* Calcuate and return to function value f(x,y) */

}

5.2

#include

#include

#include

#define MAXSIZE 50

void input (double a[][MAXSIZE],double v[],long n);

void matrix_product (double a[][MAXSIZE],double u[],double v[],long n);

void normalization (double u[],double v[],long n,double * pm1);

5

void output (double v[],long n,double m1);

void main (void)

{

double a[MAXSIZE][MAXSIZE],u[MAXSIZE],v[MAXSIZE];

double epsilon,m0,m1;

long n,maxk,k;

printf(\"\\nPlease input the square A's order number: \");

scanf(\"%ld\

input(a,v,n);

printf(\"\\nPlease input the max number of iteratings:\");

scanf(\"%ld\

printf(\"\\nPlease input the main eigenvalue'sprecision:\");

scanf(\"%lf\

6

matrix_product(a,u,v,n);

normalization(u,v,n,&m1);

for (k=1;k<=maxk;k++)

{

m0=m1;

matrix_product(a,u,v,n);

normalization(u,v,n,&m1);

if (fabs(m0-m1)<=epsilon) break;

}

if (k<=maxk)

output (v,n,m1);

else

7

printf(\"\\nThe number of iterations has the ultra limit.\");

}

/* A subroutine 1: Read in square A and the initial vector */

void input (double a[][MAXSIZE],double v[],long n)

{

long i,j;

printf(\"\\nPlease input %ld order square A:\\n\

for (i=0;i<=n-1;i++)

for (j=0;j<=n-1;j++)

scanf(\"%lf\j]);

printf(\"\\nPlease input the initial iterations vector:\");

for (i=0;i<=n-1;i++)

scanf(\"%lf\

8

}

/* A subroutine 2: Calcuate U=AV */

void matrix_product (double a[][MAXSIZE],double u[],double v[],long n)

{

long i,j;

for (i=0;i<=n-1;i++)

{

u[i]=0;

for (j=0;j<=n-1;j++)

u[i]+=a[i][j]*v[j];

}

}

/* A subroutine 3: To normalized the vector U’9

s length */

void normalization (double u[],double v[],long n,double * pm1)

{

long i;

*pm1=u[0];

for (i=1;i<=n-1;i++)

if(fabs(*pm1)*pm1=u[i];

for (i=0;i<=n-1;i++)

v[i]=u[i]/(*pm1);

}

/* A subroutine 4: Output the result of calcuate */

void output (double v[],long n,double m1)

{

10

long i;

printf(\"\\nThe square A's eigenvalue is about: %f\

printf(\"\\nIts corresponding eigenvectors is about: \\n\");

for (i=0;i<=n-1;i++)

printf(\" %lf\getch();

}

11

12