}
//定义输入数据函数
voidsolution:
:
solution_input（）
{
/**********将矩阵A存放到c[5][501]中*******/
for（i=0;i<5;i++）
{
for（j=0;j<501;j++）
{c[i][j]=0;}
}
for（j=0;j<501;j++）
{c[2][j]=（（1.64-0.024*（j+1））*sin（0.2*（j+1））-0.64*exp（0.1/（j+1）））;}
for（j=1;j<501;j++）
{c[1][j]=0.16;}
for（j=2;j<501;j++）
{c[0][j]=-0.064;}
for（j=0;j<500;j++）
{c[3][j]=0.16;}
for（j=0;j<499;j++）
{c[4][j]=-0.064;}
/**********将矩阵c[5][501]备份*******/
for（i=0;i<5;i++）
{
for（j=0;j<501;j++）
{a[i][j]=c[i][j];}
}
}
//定义幂法函数，自变量r是平移量
doublesolution:
:
solution_powermethod（doubler）
{
/**********初始化特征向量*******/
for（i=0;i<501;i++）
{v[i]=1;}
b=0;
/**********幂法迭代过程*******/
for（;;）
{
sum=0;
for（i=0;i<501;i++）
{sum+=pow（v[i],2）;}
p=sqrt（sum）;
for（i=0;i<501;i++）
{y[i]=v[i]/p;}
for（i=0;i<501;i++）
{
sum=0;
for（j=0;j<501;j++）
{
if（fabs（i-j）<3）
{
if（i==j）
sum+=（c[i-j+2][j]-r）*y[j];
else
sum+=c[i-j+2][j]*y[j];
}
}
v[i]=sum;
}
d=0;
for（i=0;i<501;i++）
{
d+=y[i]*v[i];
}
if（fabs（d-b）/fabs（d）<=1e-12）//设置迭代精度
break;
b=d;
}
returnd+r;//返回值即为按模最大特征值
}
//定义反幂法函数，自变量r为平移量
doublesolution:
:
solution_inversepowermethod（doubler）
{
/**********初始化各向量*******/
for（j=0;j<501;j++）
{x[j]=0;}
for（j=0;j<501;j++）
{y[j]=0;}
for（i=0;i<5;i++）
{
for（j=0;j<501;j++）
{c[i][j]=a[i][j];}
}
for（j=0;j<501;j++）
{
c[2][j]=c[2][j]-r;
}
for（i=0;i<501;i++）
{v[i]=1;}
b=0;
/**********先对矩阵A进行Doolittle分解，分解所得的L、U矩阵仍存放到c矩阵中*******/
for（k=0;k<501;k++）
{
for（j=k;j<=solution:
:
s_min（k+2,500）;j++）
{
sum=0;
for（t=solution:
:
s_max（0,solution:
:
s_max（k-2,j-2））;t<=k-1;t++）
{sum+=c[k-t+2][t]*c[t-j+2][j];}
c[k-j+2][j]=c[k-j+2][j]-sum;
}
for（i=k+1;i<=solution:
:
s_min（k+2,500）;i++）
{
if（k>=500）
break;
sum=0;
for（t=solution:
:
s_max（0,solution:
:
s_max（k-2,i-2））;t<=k-1;t++）
{sum+=c[i-t+2][t]*c[t-k+2][k];}
c[i-k+2][k]=（c[i-k+2][k]-sum）/c[2][k];
}
}
/**********求第k代特征向量y，并备份*******/
for（;;）
{
sum=0;
for（i=0;i<501;i++）
{sum+=pow（v[i],2）;}
q=sqrt（sum）;
for（i=0;i<501;i++）
{
y[i]=v[i]/q;
x[i]=y[i];
}
/**********用Doolittle分解法求解第k代迭代向量v以及特征值*******/
for（i=1;i<501;i++）
{
sum=0;
for（t=solution:
:
s_max（0,i-2）;t<=（i-1）;t++）
{sum+=c[i-t+2][t]*y[t];}
y[i]=y[i]-sum;
}
v[500]=y[500]/c[2][500];
for（i=499;i>=0;i--）
{
sum=0;
for（t=（i+1）;t<=solution:
:
s_min（i+2,500）;t++）
{sum+=c[i-t+2][t]*v[t];}
v[i]=（y[i]-sum）/c[2][i];
}
d=0;
for（i=0;i<501;i++）
{
d+=x[i]*v[i];
}
if（fabs（（1/d）-（1/b））/fabs（1/d）<=1e-12）//设置迭代精度
break;
b=d;
}
returnr+1/d;//返回值即为按模最小特征值
}
//定义求行列式函数
doublesolution:
:
solution_detA（）
{
/********将备份的矩阵A的值赋给矩阵c*********/
for（i=0;i<5;i++）
{
for（j=0;j<501;j++）
{c[i][j]=a[i][j];}
}
/**********先对矩阵A进行Doolittle分解，分解所得的L、U矩阵仍存放到c矩阵中*******/
for（k=0;k<501;k++）
{
for（j=k;j<=solution:
:
s_min（k+2,500）;j++）
{
sum=0;
for（t=solution:
:
s_max（0,solution:
:
s_max（k-2,j-2））;t<=k-1;t++）
{sum+=c[k-t+2][t]*c[t-j+2][j];}
c[k-j+2][j]=c[k-j+2][j]-sum;
}
for（i=k+1;i<=solution:
:
s_min（k+2,500）;i++）
{
if（k>=500）
break;
sum=0;
for（t=solution:
:
s_max（0,solution:
:
s_max（k-2,i-2））;t<=k-1;t++）
{sum+=c[i-t+2][t]*c[t-k+2][k];}
c[i-k+2][k]=（c[i-k+2][k]-sum）/c[2][k];
}
}
/**********detA为矩阵u对角线上运算的乘积*******/
sum=1;
for（j=0;j<501;j++）
{sum*=c[2][j];}
returnsum;
}
//定义取最小值函数
doublesolution:
:
s_min（doublea,doubleb）
{
if（a<=b）returna;
elsereturnb;
}
//定义取最大值函数
doublesolution:
:
s_max（doublea,doubleb）
{
if（a>=b）returna;
elsereturnb;
}
三、算法运行结果
将计算结果分别输出到屏幕和txt文件中，如下所示：
************第1小题结果*************
Lamda_1=-1.070011361502e+001
Lamda_501=9.724634098777e+000
Lamda_s=-5.557910794230e-003
************第2小题结果*************
L[1]=-1.018293403315e+001
L[2]=-9.585707425068e+000
L[3]=-9.172672423928e+000
L[4]=-8.652284007898e+000
L[5]=-8.093483808675e+000
L[6]=-7.659405407692e+000
L[7]=-7.119684648691e+000
L[8]=-6.611764339397e+000
L[9]=-6.066103226595e+000
L[10]=-5.585101052628e+000
L[11]=-5.114083529812e+000
L[12]=-4.578872176865e+000
L[13]=-4.096470926260e+000
L[14]=-3.554211215751e+000
L[15]=-3.041090018133e+000
L[16]=-2.533970311130e+000
L[17]=-2.003230769563e+000
L[18]=-1.503557611227e+000
L[19]=-9.935586060075e-001
L[20]=-4.870426738850e-001
L[21]=2.231736249575e-002
L[22]=5.324174742069e-001
L[23]=1.052898962693e+000
L[24]=1.589445881881e+000
L[25]=2.060330460274e+000
L[26]=2.558075597073e+000
L[27]=3.080240509307e+000
L[28]=3.613620867692e+000
L[29]=4.091378510451e+000
L[30]=4.603035378279e+000
L[31]=5.132924283898e+000
L[32]=5.594906348083e+000
L[33]=6.080933857027e+000
L[34]=6.680354092112e+000
L[35]=7.293877448127e+000
L[36]=7.717111714236e+000
L[37]=8.225220014050e+000
L[38]=8.648666065193e+000
L[39]=9.254200344575e+000
************第3小题结果*************
detA=2.772786141752e+118
condA=1.925204273902e+003
***************实习作业第一题完成******************
四、讨论迭代初始向量的选取对计算结果的影响
矩阵的初始向量选择，对结果的影响很大，选择不同的初始向量可能会得到不同阶的特征值。
以幂法为例（反幂法原理相同），常见的初始向量选择有两种：
1、
2、
，此处以
为例讨论。
通过第一种方法构造迭代初始向量，求得第一小题的结果是：
************第1小题结果*************
Lamda_1=-1.070011361502e+001
Lamda_501=9.724634098777e+000
Lamda_s=-5.557910794230e-003
通过第二种方法构造迭代初始向量，求得第二小题的结果是：
************第1小题结果*************
Lamda_1=-2.080981085336e+000
Lamda_501=9.978750038032e-001
Lamda_s=2.668886923785e-002
结果发现两种计算结果相差较大，而第一种初始向量得到的特征值模更大，实际上，此种初始化向量的方法能得到最准确的按模最大特征值，但迭代数次会较第2种方法多很多。
当增加迭代初始向量中非零元素的个数时，计算出的结果会越接近准确值。
实际上，当初始向量
中的0元素较多时，可能
的情况较为普遍，许多
都有可能等于0，此时计算出的结果便与最大特征值差距较大。
实际操作中可以选择不同形式的初始向量，利用求得的特征值，找到其中最大的即为需求的按模最大特征值。
升级会员 