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

《数值分析计算实习第一次大作业.doc》由会员分享，可在线阅读，更多相关《数值分析计算实习第一次大作业.doc（9页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、学生：黄小琦 学号：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执行 j=k,k+1,，min(k+2,500) i=k+1,k+2,，min(k+2,500)；kn4）若平移量x=0，即对矩阵A进行LU分解，所得a(5,501)第二行元素的乘积即为矩阵A的行列

2、式的值，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；利用3）矩阵的LU分解结

3、果求解方程组,求解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)存储于原来的存储单元内，覆盖原来的值，执行如下循环：对k=0，1，2，500（500为最大迭代次数，超过这个次数仍不

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

5、其中close(39)用来存放A中与数s(39)最接近的特征值，执行如下循环： 对k=0，1，38 ； 调用反幂法子函数，取平移量； ，其中q为反幂法子函数的带回值；6）输出和cond(A)2，并输出所有close（39）；全部源程序：#include#includedouble min(double x) 编辑反幂法子程序，该子程序中也包含的矩阵的LU分解double a5501=0,u501,y501,z501,b=0,n=0;double m,s,sum,deta,comp,q=0;int i,j,k,t,temp,min,max;for(i=0;i501;i+)if(i+2501) a

6、0i+2=-0.064;if(i+1=0) a3i-1=0.16;if(i-2=0) a4i-2=-0.064;for(k=0;k501;k+)min=(k+2500)?k+2:500;for(j=k;jk-2)?0:k-2;max=(tempj-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;i501;i+) ui=1;if(k500)for(i=k+1;ii-2)?0:i-2;max=(tempk-2)?temp:k-2;for(t=max;t=k-1;t+)sum=sum

7、+ai-t+2t*at-k+2k;ai-k+2k=(ai-k+2k-sum)/a2k;deta=1;if(x=0)for(i=0;i501;i+)deta=a2i*deta;printf(deta=%.12en,deta);if(deta=1e-12) return(x);elsefor(k=0;k500;k+)comp=q;sum=0;for(j=0;j501;j+)sum=uj*uj+sum;n=sqrt(sum);for(i=0;i501;i+)yi=ui/n;z0=y0;for(i=1;ii-2)?0:i-2;sum=0;for(t=max;t=0;i-)min=(i+2500)?i+

8、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;i501;i+)b=yi*ui+b;q=1/b;if(fabs(comp-q)=1e-12) break;return(q);double max(double x) 编辑幂法子程序，带入平移量的数值，输出按模最大特征值double a5501=0,v501,w501=0,m=0,c=0;double temp,sum,d,comp,s;int i,j,k;for(i=0;i501;i+)if(i+2501) a0i+2=-0.064;if

9、(i+1=0) a3i-1=0.16;if(i-2=0) a4i-2=-0.064;for(i=0;i501;i+) vi=1;for(k=0;k500;k+)comp=c;sum=0;for(j=0;j501;j+)sum=vj*vj+sum;m=sqrt(sum);for(i=0;i501;i+)wi=vi/m;for(i=0;i501;i+)d=0;for(j=0;j=0&i-j+25)d=wj*ai-j+2j+d;vi=d;c=0;for(j=0;j501;j+)c=wj*vj+c;temp=fabs(c-comp);if(temp0)root501=root;root1=comp;e

10、lseroot1=root;root501=comp;printf(root1=%.12e,root501=%.12e,roots=%.12en,root1,root501,roots);printf(conda=%.12en,conda);for(k=0;k39;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.925204

11、272920e+003；对于k=1,2,39，依次为如下值：1.018293403315e+001 9.585707425059e+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.554211215751

12、e+000 3.041090018133e+0002.533974279002e+000 2.003230769562e+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+

13、000 4.603035378279e+000 5.132924283898e+000 5.594906348083e+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，会使叠加的结果变小，对应于幂法得

14、到的特征值小于实际的特征值，而对反幂法所得特征值会大于实际的特征值。2. 将幂法和反幂法分别编辑成两个子程序，并且将LU分解放在反幂法的子程序中，每次运行幂法和反幂法之前都对二维数组a（5,501）赋值，在赋值时对a（5,501）第三行元素，即原矩阵A的对角线元素设置平移量，这样在求解时可以直接带入平移量，得到返回值为平移后矩阵的特征值，避免的程序的重复编写。3. 带状矩阵A存储时只存储带内元素，存为a（5,501），Aij=ai-j+2,j，应注意在C语言中存储是下标的变化。4. 对原矩阵A进行原点平移时，应该先判断平移后是否使得矩阵的特征值变为0，若平移后矩阵特征值变为0，则平移量就是A的

15、 特征值；若平移后矩阵的特征值不是0，才能进行反幂法的计算，以免产生错误。5. C语言中利用循环叠加时要注意初始变量清零，否则会在初始变量上一直叠加，使计算结果错误。6. 反幂法中利用LU分解时，Ly=b，Ux=y，为了节省存储空间，会将计算所得y值存储于原来b 中，覆盖原来的b值；但是在反幂法运算求解时，必须设置中间数组存储，不能将计算所得的中间数组直接存储于y中，因为y还在下面的运算中使用，若将y改变，将会得到错误的结果。7. 在幂法和反幂法中，设置一个比较的变量comp（取compare的意思）来存储前一次运行得到的特征值，与后一次运行得到的相比较，若小于给定的精度水平，则停止迭代。8. 有些C程序运行时与电脑配置有关，同一个C程序在不同电脑上有的能运行，有的不能运行，会出现“16位MS-DOS子系统遇到无效的指令”，我所编译的程序在我的电脑运行上就会出现此错误，但换一台电脑可以正常运行。