北航数值分析大作业一(2)

for(i=0;i<501;i++) B[i]=A[i]-mu;

beta_k=inv_power(B)+mu;

printf(\\\t= %.12e\\t迭代次数：%d\\n\ }

printf(\ //第三问

printf(\ }

double power(double a[501]) //幂法 {

int i=0,N=5000;

double b=0.16,c=-0.064; double u[501],y[501]; double m=1,beta; for(i=0;i<501;i++) u[i]=1; j=0;

while(j

u[0]=a[0]*y[0]+b*y[1]+c*y[2];

u[1]=b*y[0]+a[1]*y[1]+b*y[2]+c*y[3];

u[499]=c*y[497]+b*y[498]+a[499]*y[499]+b*y[500]; u[500]=c*y[498]+b*y[499]+a[500]*y[500]; for(i=2;i<499;i++)

{u[i]=c*y[i-2]+b*y[i-1]+a[i]*y[i]+b*y[i+1]+c*y[i+2];} beta=0;

for(i=0;i<501;i++) { if(fabs(u[i])>=fabs(beta)) beta=u[i]; }

if(beta<0)

if(fabs(fabs(beta)-fabs(m))/fabs(beta)

if(fabs(beta-m)/fabs(beta)

return beta; }

double inv_power(double a[501]) //反幂法 {

double p[501],r[501],t[501],q[501],u[501],y[501]; double beta,m=1; int i,N=1000;

p[0]=a[0];t[0]=b/p[0];r[1]=b; p[1]=a[1]-r[1]*t[0]; q[0]=c/p[0];q[1]=c/p[1]; t[1]=(b-r[1]*q[0])/p[1];

for(i=2;i<501;i++) { r[i]=b-c*t[i-2]; p[i]=a[i]-c*q[i-2]-r[i]*t[i-1]; q[i]=c/p[i]; t[i]=(b-r[i]*q[i-1])/p[i]; }

for(i=0;i<501;i++) u[i]=1; j=0;

while(j=0;i--) u[i]=u[i]-t[i]*u[i+1]-q[i]*u[i+2];

beta=0;

for(i=0;i<501;i++) { if(fabs(u[i])>=fabs(beta)) beta=u[i]; }

if(beta<0)

if(fabs(fabs(beta)-fabs(m))/fabs(beta)

if(fabs(beta-m)/fabs(beta)

return 1/beta; }

double det(double a[501]) //求det {

double det_A=1;

double p[501],r[501],t[501],q[501]; int i;

p[0]=a[0];t[0]=b/p[0];r[1]=b; p[1]=a[1]-r[1]*t[0]; q[0]=c/p[0];q[1]=c/p[1]; t[1]=(b-r[1]*q[0])/p[1];

for(i=2;i<501;i++) { r[i]=b-c*t[i-2]; p[i]=a[i]-c*q[i-2]-r[i]*t[i-1]; q[i]=c/p[i]; t[i]=(b-r[i]*q[i-1])/p[i]; }

for(i=0;i<501;i++) det_A=det_A*p[i]; return det_A; }

四 程序结果

五 计算过程中的现象

使用 |?k??k?1|/|?k|??作为终止迭代条件时，出现迭代无法终止的情况，通过调试发现按模最大特征值为负时，当k充分大后，迭代向量uk各分量不断变号，使得?k与?k?1异号，判别式|?k??k?1|/|?k|不收敛。

因此将终止迭代条件修改为||?k|?|?k?1||/|?k|??，程序实现如下： if(beta<0)

if(fabs(fabs(beta)-fabs(m))/fabs(beta)

if(fabs(beta-m)/fabs(beta)

从迭代次数可以看出?1与?501收敛较慢，由按模最大特征值与按模次大特征值的比值越小，收敛速度越慢，可知存在与?1和?501的模相近的特征值。