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

1、【精品文档】如有侵权，请联系网站删除，仅供学习与交流数值分析计算实习第一次大作业.精品文档.学生：黄小琦 学号：S201150110 联系方式法：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执行 j=k,k+1,，min(k+2,500) i=k+1,k+2,，min(k+2,500)；k<

2、;n4）若平移量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来存储前一

3、次运行所得的特征值，以便与第二次运行的值相比较，求出误差范围，comp=q；利用3）矩阵的LU分解结果求解方程组,求解Lz=y和Uu=z的计算公式是： (i=1，2，500) (i=499，498，0)，判断，若，结束循环转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)存储于原来的存储

4、单元内，覆盖原来的值，执行如下循环：对k=0，1，2，500（500为最大迭代次数，超过这个次数仍不收敛则停止计算）设置一个数comp来存储前一次运行所得的特征值，以便与第二次运行的值相比较，求出误差范围，comp=c； ，但是由于A中的元素已被存为a(5,501)，对于中每一个分量和中的需对迭代公式做一下变换： 对i=0，1，2，500， 判断，若，结束循环转7），否则继续循环7）将c值带回主函数。3编辑主函数1）调用反幂法子函数，令平移量为0，得到矩阵A的按模最小特征值，并输出A的行列式的值；2）调用幂法子函数，令其平移量为0，得到按模最大特征值；3）调用幂法子函数，令平移量为，得到的，其

5、中c是幂法子函数的带回值；4）如果，则，；若，则，；5）设置数组s(39)，close(39)，其中close(39)用来存放A中与数s(39)最接近的特征值，执行如下循环： 对k=0，1，38 ； 调用反幂法子函数，取平移量； ，其中q为反幂法子函数的带回值；6）输出和cond(A)2，并输出所有close（39）；全部源程序：#include<stdio.h>#include<math.h>double min(double x) 编辑反幂法子程序，该子程序中也包含的矩阵的LU分解double a5501=0,u501,y501,z501,b=0,n=0;doubl

6、e 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) a0i+2=-0.064;if(i+1<501) a1i+1=0.16;s=sin(0.2*(i+1);m=exp(0.1/(i+1);a2i=(1.64-0.024*(i+1)*s-0.64*m-x;if(i-1>=0) a3i-1=0.16;if(i-2>=0) a4i-2=-0.064;for(k=0;k<501;k+)min=(k+2<500)?k+2:500;for(j=k;j<=m

7、in;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+ak-t+2t*at-j+2j;ak-j+2j=ak-j+2j-sum;for(i=0;i<501;i+) ui=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+ai-t+2t*at-k+2k;ai-k+2k=(

8、ai-k+2k-sum)/a2k;deta=1;if(x=0)for(i=0;i<501;i+)deta=a2i*deta;printf("deta=%.12en",deta);if(deta<=1e-12) return(x);elsefor(k=0;k<500;k+)comp=q;sum=0;for(j=0;j<501;j+)sum=uj*uj+sum;n=sqrt(sum);for(i=0;i<501;i+)yi=ui/n;z0=y0;for(i=1;i<501;i+)max=(0>i-2)?0:i-2;sum=0;for(t

9、=max;t<=i-1;t+)sum=ai-t+2t*zt+sum;zi=yi-sum;u500=z500/a2500;for(i=499;i>=0;i-)min=(i+2<500)?i+2:500;sum=0;for(t=i+1;t<=min;t+)sum=sum+ai-t+2t*ut;ui=(zi-sum)/a2i;b=0;for(i=0;i<501;i+)b=yi*ui+b;q=1/b;if(fabs(comp-q)<=1e-12) break;return(q);double max(double x) 编辑幂法子程序，带入平移量的数值，输出按模最大

10、特征值double a5501=0,v501,w501=0,m=0,c=0;double temp,sum,d,comp,s;int i,j,k;for(i=0;i<501;i+)if(i+2<501) a0i+2=-0.064;if(i+1<501) a1i+1=0.16;s=sin(0.2*(i+1);m=exp(0.1/(i+1);a2i=(1.64-0.024*(i+1)*s-0.64*m-x;if(i-1>=0) a3i-1=0.16;if(i-2>=0) a4i-2=-0.064;for(i=0;i<501;i+) vi=1;for(k=0;k&

11、lt;500;k+)comp=c;sum=0;for(j=0;j<501;j+)sum=vj*vj+sum;m=sqrt(sum);for(i=0;i<501;i+)wi=vi/m;for(i=0;i<501;i+)d=0;for(j=0;j<501;j+)if(i-j+2>=0&&i-j+2<5)d=wj*ai-j+2j+d;vi=d;c=0;for(j=0;j<501;j+)c=wj*vj+c;temp=fabs(c-comp);if(temp<=1e-12) break;return(c);void main() 编辑主函数

12、，调用反幂法子程序和幂法子程序，求解特征值int k;double s39,close39,conda,root,roots,root1,root501,comp;roots=min(0);root=max(0);conda=fabs(root/roots);comp=max(root)+root;if(root>0)root501=root;root1=comp;elseroot1=root;root501=comp;printf("root1=%.12e,root501=%.12e,roots=%.12en",root1,root501,roots);printf

13、("conda=%.12en",conda);for(k=0;k<39;k+)sk=root1+(k+1)*(root501-root1)/40;closek=min(sk)+sk;if(k%3=0) printf("n");printf("%.12e ",closek);程序的输出结果：detA=2.772786141752e+118；1.070011361514e+001；cond(A)2=1.925204272920e+003；对于k=1,2,39，依次为如下值：1.018293403315e+001 9.58570742

14、5059e+000 9.172672423928e+0008.652284007898e+000 8.093483808672e+000 7.659405407692e+0007.119684648691e+000 6.611764339397e+000 6.066103226595e+0005.585101052628e+000 5.114083529812e+000 4.578872176865e+0004.096470926271e+000 3.554211215751e+000 3.041090018133e+0002.533974279002e+000 2.003230769562e

15、+000 1.503557611227e+0009.935586060075e-001 4.870426738850e-001 2.231736249575e-002 5.324174742069e-001 1.052898962694e+000 1.589445881881e+000 2.060330460171e+000 2.558075597072e+000 3.080240509307e+000 3.613620867692e+000 4.091378510451e+000 4.603035378279e+000 5.132924283898e+000 5.594906348083e+

16、000 6.080933857026e+000 6.680354092112e+000 7.293877448128e+000 7.717111714236e+000 8.225220014050e+000 8.648666065193e+000 9.254200344575e+000讨论：1. 迭代初始向量中所含0的个数对计算结果有很大的影响，所含0越多，利用幂法和反幂法求得的数就越小。因为若初始时将该分量赋值为零，在计算中，该分量仍然为0，会使叠加的结果变小，对应于幂法得到的特征值小于实际的特征值，而对反幂法所得特征值会大于实际的特征值。2. 将幂法和反幂法分别编辑成两个子程序，并且将LU

17、分解放在反幂法的子程序中，每次运行幂法和反幂法之前都对二维数组a（5,501）赋值，在赋值时对a（5,501）第三行元素，即原矩阵A的对角线元素设置平移量，这样在求解时可以直接带入平移量，得到返回值为平移后矩阵的特征值，避免的程序的重复编写。3. 带状矩阵A存储时只存储带内元素，存为a（5,501），Aij=ai-j+2,j，应注意在C语言中存储是下标的变化。4. 对原矩阵A进行原点平移时，应该先判断平移后是否使得矩阵的特征值变为0，若平移后矩阵特征值变为0，则平移量就是A的 特征值；若平移后矩阵的特征值不是0，才能进行反幂法的计算，以免产生错误。5. C语言中利用循环叠加时要注意初始变量清零，否则会在初始变量上一直叠加，使计算结果错误。6. 反幂法中利用LU分解时，Ly=b，Ux=y，为了节省存储空间，会将计算所得y值存储于原来b