return（t）;
}
voidfuzhi（doubleA[5][501]）/赋值函数，对矩阵进行旋转45度，然后赋值，不储存多余的非0数值/
{
ints=2,r=2,i,j;
doublearray[501],b=0.16,c=-0.064;
A[0][0]=A[0][1]=A[1][0]=A[3][500]=A[4][499]=A[4][500]=0;
for（j=s;j<501;j++）
A[0][j]=c;
for（j=s-1;j<501;j++）
A[1][j]=b;
for（i=0;i<501;i++）
array[i]=（1.64-0.024*i）*sin（0.2*i）-0.64*exp（0.1/i）;
for（j=s-2;j<501;j++）
A[2][j]=array[j];
for（j=0;j<500;j++）
A[3][j]=b;
for（j=0;j<499;j++）
A[4][j]=c;
}
doublemifa（doubleA[5][501]）/幂法，求模最大的特征值/
{
ints=2,r=2,i,j;
doublesum=0,f=1,u[501],y[501],b[501],d,a2,a1=0;
for（i=0;i<501;i++）
u[i]=1;
while（f>1.0e-12）
{
a1=a2;
for（i=0;i<501;i++）
{b[i]=u[i]*u[i];
sum+=b[i];
}
d=sqrt（sum）;
for（i=0;i<501;i++）
y[i]=u[i]/d;
for（i=1;i<501;i++）
{
sum=0;
for（j=max（i-r,1）;j<=min（i+s,500）;j++）
sum+=A[i-j+s][j-1]*y[j-1]
u[i-1]=sum;
}
sum=0;
for（i=0;i<501;i++）
sum+=y[i]*u[i];
a2=sum;
f=fabs（a2-a1）/fabs（a2）;
}
returna2;
}
voidtrianglefj（doubleA[5][501],doublex[501],doubley[501]）/三角分解函数，对A矩阵进行LU分解，然后通过AX=Y求解X/
{
ints=2,r=2,i,j,k,t;
doublec[5][501],b[501],sum;
for（i=0;i<5;i++）
{
for（j=0;j<501;j++）
c[i][j]=A[i][j];
}
for（i=0;i<501;i++）
b[i]=y[i];
for（k=1;k<=501;k++）
{
for（j=k;j<=min（k+s,501）;j++）
{
sum=0;
for（t=max3（1,k-r,j-s）;t<=k-1;t++）
{
sum=c[k-t+s][t-1]*c[t-j+s][j-1];
c[k-j+s][j-1]=c[k-j+s][j-1]-sum;
}
}
for（i=k+1;i<=min（k+r,501）;i++）
{
sum=0;
for（t=max3（1,i-r,k-s）;t<=k-1;t++）
sum=c[i-t+s][t-1]*c[t-k+s][k-1];
c[i-k+s][k-1]=（c[i-k+s][k-1]-sum）/c[s][k-1];
}
}
for（i=2;i<=501;i++）
{sum=0;
for（t=max（1,i-r）;t<=i-1;t++）
sum+=c[i-t+s][t-1]*b[t-1];
b[i]=b[i]-sum;
}
x[500]=b[500]/c[s][500];
for（i=500;i>=1;i--）
{sum=0;
for（t=i+1;t<=min（i+s,501）;t++）
sum=c[i-t+s][t-1]*x[t-1];
x[i-1]=（b[i-1]-sum）/c[s][i-1];
}
}
doublefanmf（doubleA[5][501]）/反幂法函数，求解模最小的特征值/
{
ints=2,r=2,i;
doublea2,a1=0,sum=0,u[501],y[501],b[501],f=1,d;
for（i=0;i<501;i++）
u[i]=1;
while（f>1.0e-12）
{
a1=a2;
for（i=0;i<501;i++）
{
b[i]=u[i]*u[i];
sum+=b[i];
}
d=sqrt（sum）;
for（i=0;i<501;i++）
y[i]=u[i]/d;
trianglefj（A,u,y）;
sum=0;
for（i=0;i<501;i++）
sum+=y[i]*u[i];
a2=sum;
f=fabs（a2-a1）/fabs（a2）;
}
return1/a2;
}
doubledet（doublec[5][501]）/求值函数，求解矩阵C的行列式的值并返回。
/
{
inti,j,k,r=2,s=2,t;
doubledet=1,sum;
for（k=1;k<=501;k++）
{
for（j=k;j<=min（k+s,501）;j++）
{
sum=0;
for（t=max3（1,k-r,j-s）;t<=k-1;t++）
{
sum=c[k-t+s][t-1]*c[t-j+s][j-1];
c[k-j+s][j-1]=c[k-j+s][j-1]-sum;
}
}
for（i=k+1;i<=min（k+r,501）;i++）
{
sum=0;
for（t=max3（1,i-r,k-s）;t<=k-1;t++）
sum=c[i-t+s][t-1]*c[t-k+s][k-1];
c[i-k+s][k-1]=（c[i-k+s][k-1]-sum）/c[s][k-1];
}
}
for（i=0;i<501;i++）
det*=c[2][i];
returndet;
}
main（）
{
doublelanmda1,lanmdas,lanmda501,k,k1,A[5][501],u[39]={0},h,cond,Det;
inti,j;
fuzhi（A）;/对A矩阵进行旋转赋值/
k=mifa（A）;/用幂法求解A的模最大特征值，并且赋值给K/
for（i=0;i<501;i++）
A[2][i]=A[2][i]-k;/把A矩阵第3行所有元素都平移k个单位/
k1=mifa（A）+k;/用幂法再求解平移后的模最大的特征值，赋值给k1，后面进行判断赋值给lanmda1和lanmda501/
for（i=0;i<501;i++）
A[2][i]=A[2][i]+k;/将平移矩阵还原/
if（k>k1）/判断k和k1的值，并相应的赋值给lanmda1和lanmda501/
{
lanmda1=k1;
lanmda501=k;
}
else
{
lanmda1=k;
lanmda501=k1;
}
fuzhi（A）;
lanmdas=fanmf（A）;/用反幂法直接求解模最小的特征值，并赋值给lanmdas/
printf（"最大特征值λ为%.12e最小特征值λ为%.12e\n",lanmda501，lanmda1）;
printf（"λs=%.12e\n",lanmdas）;
fuzhi（A）;
for（i=0;i<=39;i++）/通过循环用反幂法求解离各个u值最近的特征值/
{
u[i]=lanmda1+（i+1）*（lanmda501-lanmda1）/40;
for（j=0;j<501;j++）
A[2][j]=A[2][j]-u[i];
h=fanmf（A）+u[i];
for（j=0;j<501;j++）
A[2][j]=A[2][j]+u[i];
printf（"A的与数u%d=%.12e最接近的特征值λ%d为%.12e\n",i+1,u[i].i+1,h）;
}
cond=fabs（mifa（A）/fanmf（A））;
fuzhi（A）;
Det=det（A）;
printf（"A的条件数cond（A）为%.12e\n",cond）;
printf（"矩阵A的行列式detA为%.12e\n",Det）;
}
三．计算截图
上图为迭代初始值u[]全取1时的值。
上图为迭代初始值全为15时的值。
上图为迭代初始值全为0.1
四．不同迭代初始值的分析
从不同的截图可以看出，无论初始值怎么选取，计算出来的结果都是一样的，也就是说迭代初始值的选取不影响特征值的计算，仅仅影响迭代的次数。
但当迭代初始值取1个1，其余500个都是0的时候，算法无法继续进行，因为计算机需要迭代的次数太多，计算机无法计算从而得不出正确的结果。
升级会员 