北航研究生数值分析作业第一题

北航研究⽣数值分析作业第⼀题
北航研究⽣数值分析作业第⼀题：
⼀、算法设计⽅案
1.要求计算矩阵的最⼤最⼩特征值，通过幂法求得模最⼤的特征值，进⾏⼀定
判断即得所求结果；
2.求解与给定数值接近的特征值，可以该数做漂移量，新数组特征值倒数的绝
对值满⾜反幂法的要求，故通过反幂法即可求得；
3.反幂法计算时需要⽅程求解中间过渡向量，需设计Doolite分解求解；
4.|A|=|B||C|，故要求解矩阵的秩，只需将Doolite分解后的U矩阵的对⾓线相
乘即为矩阵的Det。

算法编译环境：vlsual c++6.0
需要编译函数：幂法，反幂法，Doolite分解及⽅程的求解
⼆、源程序如下：
#include
#include
#include
#include
int Max(int value1,int value2);
int Min(int value1,int value2);
void Transform(double A[5][501]);
double mifa(double A[5][501]);
void daizhuangdoolite(double A[5][501],double x[501],double b[501]); double fanmifa(double A[5][501]); double Det(double A[5][501]);
/***定义2个判断⼤⼩的函数，便于以后调⽤***/
int Max(int value1,int value2)
{
return((value1>value2)?value1:value2);
}
int Min(int value1,int value2)
{
return ((value1
}
/*****************************************/
/***将矩阵值转存在⼀个数组⾥，节省空间***/
void Transform(double A[5][501],double b,double c)
{
int i=0,j=0;
A[i][j]=0,A[i][j+1]=0;
for(j=2;j<=500;j++)
A[i][j]=c;
i++;
j=0;
A[i][j]=0;
for(j=1;j<=500;j++)
A[i][j]=b;
i++;
for(j=0;j<=500;j++)
A[i][j]=(1.64-0.024*(j+1))*sin(0.2*(j+1))-0.64*exp(0.1/(j+1)); i++;
for(j=0;j<=499;j++)
A[i][j]=b;
A[i][j]=0;
i++;
for(j=0;j<=498;j++)
A[i][j]=c;
A[i][j]=0,A[i][j+1]=0;
}
/***转存结束***/
//⽤于求解模最⼤的特征值，幂法
double mifa(double A[5][501])
{
int s=2,r=2,m=0,i,j;
double b2,b1=0,sum,u[501],y[501];
for (i=0;i<=500;i++)
{
u[i] = 1.0;
}
do
{
sum=0;
if(m!=0)b1=b2;
m++;
for(i=0;i<=500;i++)
sum+=u[i]*u[i];
for(i=0;i<=500;i++)
y[i]=u[i]/sqrt(sum);
for(i=0;i<=500;i++)
{
u[i]=0;
for(j=Max(i-r,0);j<=Min(i+s,500);j++)
u[i]=u[i]+A[i-j+s][j]*y[j];
}
b2=0;
for(i=0;i<=500;i++)
b2=b2+y[i]*u[i];
}
while(fabs(b2-b1)/fabs(b2)>=exp(-12));
return b2;
}
//带状DOOLITE分解，并且求解出⽅程组的解
void daizhuangdoolite(double A[5][501],double x[501],double b[501]) { int i,j,k,t,s=2,r=2;
double B[5][501],c[501];
for(i=0;i<=4;i++)
{
for(j=0;j<=500;j++)
B[i][j]=A[i][j];
}
for(i=0;i<=500;i++)
c[i]=b[i];
for(k=0;k<=500;k++)
{
for(j=k;j<=Min(k+s,500);j++)
{
for(t=Max(0,Max(k-r,j-s));t<=k-1;t++)
B[k-j+s][j]=B[k-j+s][j]-B[k-t+s][t]*B[t-j+s][j]; }
for(i=k+1;i<=Min(k+r,500);i++)
{
for(t=Max(0,Max(i-r,k-s));t<=k-1;t++)
B[i-k+s][k]=B[i-k+s][k]-B[i-t+s][t]*B[t-k+s][k]; B[i-k+s][k]=B[i-k+s][k]/B[s][k];
}
}
for(i=1;i<=500;i++)
for(t=Max(0,i-r);t<=i-1;t++)
c[i]=c[i]-B[i-t+s][t]*c[t];
x[500]=c[500]/B[s][500];
for(i=499;i>=0;i--)
{
x[i]=c[i];
for(t=i+1;t<=Min(i+s,500);t++)
x[i]=x[i]-B[i-t+s][t]*x[t];
x[i]=x[i]/B[s][i];
}
}
//⽤于求解模最⼤的特征值，反幂法double fanmifa(double A[5][501])
{
int s=2,r=2,m=0,i;
double b2,b1=0,sum=0,u[501],y[501];
for (i=0;i<=500;i++)
{
u[i] = 1.0;
}
do
{
if(m!=0)b1=b2;
m++;
sum=0;
for(i=0;i<=500;i++)
sum+=u[i]*u[i];
for(i=0;i<=500;i++)
y[i]=u[i]/sqrt(sum);
daizhuangdoolite(A,u,y);
b2=0;
for(i=0;i<=500;i++)
b2+=y[i]*u[i];
}
while(fabs(b2-b1)>=fabs(b1)*exp(-12));
return 1/b2;
}
//⾏列式的LU分解，U的主线乘积即位矩阵的DET double Det(double A[5][501]) {
int i,j,k,t,s=2,r=2;
for(k=0;k<=500;k++)
{
for(j=k;j<=Min(k+s,500);j++)
{
for(t=Max(0,Max(k-r,j-s));t<=k-1;t++)
A[k-j+s][j]=A[k-j+s][j]-A[k-t+s][t]*A[t-j+s][j];
}
for(i=k+1;i<=Min(k+r,500);i++)
{
for(t=Max(0,Max(i-r,k-s));t<=k-1;t++)
A[i-k+s][k]=A[i-k+s][k]-A[i-t+s][t]*A[t-k+s][k];
A[i-k+s][k]=A[i-k+s][k]/A[s][k];
}
}
double det=1;
for(i=0;i<=500;i++)
det*=A[s][i];
return det;
}
void main()
{
double b=0.16,c=-0.064,p,q;
int i,j;
double A[5][501];
Transform(A,b,c); //进⾏A的赋值
cout.precision(12); //定义输出精度
double lamda1,lamda501,lamdas;
double k=mifa(A);
if(k>0) //判断求得最⼤以及最⼩的特征值.如果K>0，则它为最⼤特征值值，
//并以它为偏移量再⽤⼀次幂法求得新矩阵最⼤特征值，即为最⼤ //与最⼩的特征值的差{
lamda501=k;
for(i=0;i<=500;i++)
A[2][i]=A[2][i]-k;
lamda1=mifa(A)+lamda501;
for(i=0;i<=500;i++)
A[2][i]=A[2][i]+k;
}
else //如果K<=0，则它为最⼩特征值值，并以它为偏移量再⽤⼀次幂法
//求得新矩阵最⼤特征值，即为最⼤与最⼩的特征值的差
{
lamda1=k;
for(i=0;i<=500;i++)
A[2][i]=A[2][i]-k;
lamda501=mifa(A)+lamda1;
for(i=0;i<=500;i++)
A[2][i]=A[2][i]+k;
}
lamdas=fanmifa(A);
FILE *fp=fopen("result.txt","w");
fprintf(fp,"λ1=%.12e\n",lamda1);
fprintf(fp,"λ501=%.12e\n",lamda501);
fprintf(fp,"λs=%.12e\n\n",lamdas);
fprintf(fp,"\t要求接近的值\t\t\t实际求得的特征值\n");
for(i=1;i<=39;i++) //反幂法求得与给定值接近的特征值{
p=lamda1+(i+1)*(lamda501-lamda1)/40;
for(j=0;j<=500;j++)
A[2][j]=A[2][j]-p;
q=fanmifa(A)+p;
for(j=0;j<=500;j++)
A[2][j]=A[2][j]+p;
fprintf(fp,"µ%d: %.12e λi%d: %.12e\n",i,p,i,q);
}
double cond=fabs(mifa(A)/fanmifa(A));
double det=Det(A);
fprintf(fp,"\ncond(A)=%.12e\n",cond);
fprintf(fp,"\ndetA=%.12e\n",det);
}
三、程序运⾏结果
λ1=-1.069936345952e+001
λ501=9.722283648681e+000
λs=-5.557989086521e-003
要求接近的值实际求得的特征值
µ1: -9.678281104107e+000 λi1: -9.585702058251e+000µ2: -9.167739926402e+000 λi2: -9.172672423948e+000µ3: -8.657198748697e+000 λi3: -8.652284007885e+000µ4: -8.146657570993e+000 λi4: -8.0934********e+000
µ5: -7.636116393288e+000 λi5: -7.659405420574e+000µ6: -7.125575215583e+000 λi6: -7.119684646576e+000µ7: -6.615034037878e+000 λi7: -6.611764337314e+000µ8: -6.104492860173e+000 λi8: -6.0661********e+000
µ9: -5.593951682468e+000 λi9: -5.585101045269e+000µ10: -5.0834********e+000 λi10: -5.114083539196e+000
µ11: -4.572869327058e+000 λi11: -4.578872177367e+000µ12: -4.062328149353e+000 λi12: -4.096473385708e+000µ13: -3.551786971648e+000 λi13: -3.554211216942e+000µ14: -3.0412********e+000 λi14: -3.0410********e+000
µ15: -2.530704616238e+000 λi15: -2.533970334136e+000µ16: -2.020*********e+000 λi16: -2.003230401311e+000
µ17: -1.509622260828e+000 λi17: -1.503557606947e+000
µ18: -9.990810831232e-001 λi18: -9.935585987809e-001
µ19: -4.885399054182e-001 λi19: -4.870426734583e-001
µ20: 2.200127228676e-002 λi20: 2.231736249587e-002
µ21: 5.325424499917e-001 λi21: 5.324174742068e-001
µ22: 1.043083627697e+000 λi22: 1.052898964020e+000
µ23: 1.553624805402e+000 λi23: 1.589445977158e+000
µ24: 2.064165983107e+000 λi24: 2.060330427561e+000
µ25: 2.574707160812e+000 λi25: 2.558075576223e+000
µ26: 3.0852********e+000 λi26: 3.080240508465e+000
µ27: 3.595789516221e+000 λi27: 3.613620874136e+000
µ28: 4.106330693926e+000 λi28: 4.0913********e+000
µ29: 4.616871871631e+000 λi29: 4.603035354280e+000
µ30: 5.127413049336e+000 λi30: 5.132924284378e+000
µ31: 5.637954227041e+000 λi31: 5.594906275501e+000
µ32: 6.148495404746e+000 λi32: 6.080933498348e+000
µ33: 6.659036582451e+000 λi33: 6.680354121496e+000
µ34: 7.169577760156e+000 λi34: 7.293878467852e+000
µ35: 7.680118937861e+000 λi35: 7.717111851857e+000
µ36: 8.190660115566e+000 λi36: 8.225220016407e+000
µ37: 8.701201293271e+000 λi37: 8.648665837870e+000
µ38: 9.211742470976e+000 λi38: 9.254200347303e+000
µ39: 9.722283648681e+000 λi39: 9.724634099672e+000
cond(A)=1.925042185755e+003
detA=2.772786141752e+118
四、分析
如果初始向量选择不当,将导致迭代中X1的系数等于零.但是，由于舍⼊误差的影响,经若⼲步迭代后,.按照基向量展开时,x1的系数可能不等于零。

把这⼀向量看作初始向量，⽤幂法继续求向量序列,仍然会得出预期的结果,不过收敛速度较慢.如果收敛很慢,可改换初始向量,所以初始向量的选取会影响到迭代的次数。

五、流程图。