数值分析计算实习第一次大作业

学生：黄小琦 学号：S201150110 联系方式：13261252103

算法：

1.编辑反幂法的子程序

1）设置二维数组a（5，501）来存储矩阵A 的带内元素，Aij=ai-j+2,j；

2）对于主对角线上的元素，即a2j设置平移量n，并存储到原来单元内，使得对于j=0，1，2，?，500，a2j=a2j-n；

3）对矩阵A进行LU分解，相对应于啊（5,501）中的元素执行如下循环： 对k=0，1，2，?，500执行

ak?j?2,j?ak?j?2,j?k?t?2,tt?max(0,k?2,j?2)k?1?ak?1at?j?2,j [j=k,k+1,?，min(k+2,500)]

ai?k?2,k?(ai?k?2,k?i?t?2,tt?max(0,i?2,k?2)?aat?k?2,k)/a2,k [i=k+1,k+2,?，min(k+2,500)；k

4）若平移量x=0，即对矩阵A进行LU分解，所得a(5,501)第二行元素的乘积即为矩阵A的行列式的值，detA=a0,0*a0,1*?a0,500;

5）进行反幂法迭代前，先对平移后的矩阵的行列式做判断，看其是否为0，若行列式为0，则绝对值最小的特征值为0。相对于原来的矩阵A，其特征值为平移量x。若行列式不为0，则开始进行反幂法的迭代求解；

6）取非零向量组u(501)作为初始值，取u的所有元素都为1，每迭代一次就将新求得的u(501)存储于原来的存储单元内，覆盖原来的值，执行如下循环：

对k=0，1，2，?，500（500为最大迭代次数，超过这个次数仍不收敛则停止计算） ①设置一个数comp来存储前一次运行所得的特征值，以便与第二次运行的值相比较，求出误差范围，comp=q； ②n?uTu ③y?u/n

④利用3）矩阵的LU分解结果求解方程组Au?y,求解Lz=y和Uu=z的计算公式是：

z0?y0 zi?yi?i?t?2,tt?max(0,i?2)?ai?1zt (i=1，2，?，500)

u500?y500/a2,500

min(i?2,500)ui?(zi?t?i?1?ai?t?2,tut)/a2,i (i=499，498，?，0)

⑤b?yTu，q?1 b?12⑥判断，若comp?q??，??1?10，结束循环转7），否则继续循环

7）将q值带回主函数。

2.编辑幂法的子程序

1）设置二维数组a（5，501）来存储矩阵A 的带内元素，Aij=ai-j+2,j；

2）对于主对角线上的元素，即a2j设置平移量n，并存储到原来单元内，使得对于j=0，1，2，?，500，a2j=a2j-n；

3）取非零向量组v(501)作为初始值，取u的所有元素都为1，每迭代一次就将新求得的v(501)存储于原来的存储单元内，覆盖原来的值，执行如下循环：

对k=0，1，2，?，500（500为最大迭代次数，超过这个次数仍不收敛则停止计算） ①设置一个数comp来存储前一次运行所得的特征值，以便与第二次运行的值相比较，求出误差范围，comp=c； ②m?vTv

③w?v/m

④v?Aw，但是由于A中的元素已被存为a(5,501)，对于v中每一个分量vi和w中

的wi需对迭代公式做一下变换： 对i=0，1，2，?，500，

(i?j?2??0)?(i?j?2?5) vi? ⑤c?wv

T?waij?0i?j?2,j；

?12⑥判断，若comp?c??，??1?10，结束循环转7），否则继续循环

7）将c值带回主函数。

3．编辑主函数

1）调用反幂法子函数，令平移量为0，得到矩阵A的按模最小特征值?s，并输出A的行列式的值；

2）调用幂法子函数，令其平移量为0，得到按模最大特征值?m；

3）调用幂法子函数，令平移量为???m，得到的??c??m，其中c是幂法子函数的带回值；

4）如果??0，则?501??m，?1??；若??0，则?1??m，?501??；

5）设置数组s(39)，close(39)，其中close(39)用来存放A中与数s(39)最接近的特征值，执行如下循环： 对k=0，1，?，38 ①sk??1?k?501??140；

②调用反幂法子函数，取平移量??sk； ③closek?q?sk，其中q为反幂法子函数的带回值；

6）输出?1,?501,?s和cond(A)2，并输出所有close（39）；

全部源程序：

#include #include

double min(double x) ——编辑反幂法子程序，该子程序中也包含的矩阵的LU分解 {

double a[5][501]={0},u[501],y[501],z[501],b=0,n=0; double m,s,sum,deta,comp,q=0; int i,j,k,t,temp,min,max; for(i=0;i<501;i++) {

if(i+2<501) a[0][i+2]=-0.064; if(i+1<501) a[1][i+1]=0.16; s=sin(0.2*(i+1)); m=exp(0.1/(i+1));

a[2][i]=(1.64-0.024*(i+1))*s-0.64*m-x; if(i-1>=0) a[3][i-1]=0.16; if(i-2>=0) a[4][i-2]=-0.064; }

for(k=0;k<501;k++) {

min=(k+2<500)?k+2:500; for(j=k;j<=min;j++) {

sum=0;

temp=(0>k-2)?0:k-2;

max=(temp>j-2)?temp:j-2; for(t=max;t<=k-1;t++)

sum=sum+a[k-t+2][t]*a[t-j+2][j]; a[k-j+2][j]=a[k-j+2][j]-sum; }

for(i=0;i<501;i++) u[i]=1; if(k<500) {

for(i=k+1;i<=min;i++) {

sum=0;

temp=(0>i-2)?0:i-2;

max=(temp>k-2)?temp:k-2; for(t=max;t<=k-1;t++)

sum=sum+a[i-t+2][t]*a[t-k+2][k]; a[i-k+2][k]=(a[i-k+2][k]-sum)/a[2][k]; }

} }

deta=1; if(x==0) {

for(i=0;i<501;i++) deta=a[2][i]*deta;

printf(\}

if(deta<=1e-12) return(x); else {

for(k=0;k<500;k++) {

comp=q; sum=0;

for(j=0;j<501;j++) sum=u[j]*u[j]+sum; n=sqrt(sum);

for(i=0;i<501;i++) y[i]=u[i]/n; z[0]=y[0];

for(i=1;i<501;i++) {

max=(0>i-2)?0:i-2; sum=0;

for(t=max;t<=i-1;t++) sum=a[i-t+2][t]*z[t]+sum; z[i]=y[i]-sum; }

u[500]=z[500]/a[2][500]; for(i=499;i>=0;i--) {

min=(i+2<500)?i+2:500; sum=0;

for(t=i+1;t<=min;t++) sum=sum+a[i-t+2][t]*u[t]; u[i]=(z[i]-sum)/a[2][i]; } b=0;

for(i=0;i<501;i++) b=y[i]*u[i]+b; q=1/b;

if(fabs(comp-q)<=1e-12) break;

本文来源：https://www.bwwdw.com/article/lhxv.html