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

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

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

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

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

1-A ，乘幂法对应的模最大特征值λ取倒数为λ

1，变为模最小特征值，由此得到反幂法，可反复迭代：)1()(-=k k Ay y ，)1()(/-k k y y 收敛于A 的主特征值，)(k y 为对应的特征向量。

5.1 取步长2.0=h ，用显示欧拉法求⎩⎨⎧=--='1

)0(2

y xy y y 在6.0=x 处y 的近似值。

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

运行结果为：

5.2 已知矩阵⎥⎥⎥⎦

⎤⎢⎢⎢⎣⎡=361641593642A ，用改进后的幂乘法求A 得主特征值及其对应的特征向量。

运行结果为：

程序代码

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",&a,&b);

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

scanf("%lf",&h);

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

x[0]=a;

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

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",i,x[i],i,y[i]);

getch();

}

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);

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",&n);

input(a,v,n);

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

scanf("%ld",&maxk);

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

scanf("%lf",&epsilon);

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

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",n);

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

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

scanf("%lf",&a[i][j]);

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

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

scanf("%lf",&v[i]);

}

/* 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’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 */