北航数值分析课程第一次大作业

1、数值分析A计算实习题目第一题一 算法设计方案：1.矩阵A的存储与检索将带状线性矩阵A501501转存为一个矩阵MatrixC5501 .由于C语言中数组角标都是从0开始的，所以在数组MatrixC5501中检索A的带内元素aij的方法是：A的带内元素aij=C中的元素ci-j+2,j2.求解1，501，s 首先分别使用幂法和反幂法迭代求出矩阵按摸最大和最小的特征值max和min。min即为s；如果max0,则501=max；如果max0,则1=，max+p；如果max0,则501=，max+p。3.求解A的与数k=1+k（501-1）/40的最接近的特征值ik（k=1,2，39）。 使用带原点

2、平移的反幂法，令平移量p=k，即可求出与k最接近的特征值ik。4.求解A的（谱范数）条件数cond（A）2和行列式detA。cond（A）2=|1/n|，其中1和n分别是矩阵A的模最大和最小特征值。矩阵A的行列式可先对矩阵A进行LU分解后，detA等于U所有对角线上元素的乘积。二 源程序(VS2010环境下，C+语言)#include#include#include #include #include#include#include#define E 1.0e-12 /*定义全局变量相对误差限*/int max2(int a,int b) /*求两个整型数最大值的子程序*/if(ab) ret

3、urn a;else return b; int min2(int a,int b) /*求两个整型数最小值的子程序*/if(ab) return b;else return a; int max3(int a,int b,int c) /*求三整型数最大值的子程序*/ int t;if(ab)t=a; else t=b;if(tc) t=c;return(t);void assignment(double array5501) /*将矩阵A转存为数组C5501*/ int i,j,k;/所有元素归零 for(i=0;i=4;) for(j=0;j=500;) arrayij=0;j+; i+

4、; /第0,4行赋值 for(j=2;j=500;) k=500-j; array0j=-0.064; array4k=-0.064; j+; /第1,3行赋值 for(j=1;j=500;) k=500-j; array1j=0.16; array3k=0.16; j+; /第2行赋值 for(j=0;j=500;) k=j;j+; array2k=(1.64-0.024*j)*sin(double)(0.2*j)-0.64*exp(double)(0.1/j); double mifa(double u501,double array5501,double p) /*带原点平移的幂法*/i

5、nt i,j; /* u501为初始迭代向量*/double a,b,c=0; /* array5501为矩阵A的转存矩阵*/double y501; /*p为平移量*/for(;)a=0;b=0; /*选用第一种迭代格式*/求k-1for(i=0;i=500;i+) a=a+ui*ui;a=sqrt(a); /求yk-1for(i=0;i=500;i+)yi=ui/a; /求uk for(i=0;i=500;i+)ui=0;for(j=max2(i-2,0);j=min2(i+2,500);j+)ui+=arrayi-j+2j*yj;ui=ui-p*yi; /*引入平移量*/ /求kfor(

6、i=0;i=500;i+)b+=yi*ui;if(fabs(b-c)/b)=E) /*达到精度水平，迭代终止*/break; c=b;return (b+p); /*直接返回A的特征值*/ void chuzhi(double a) /*用随机数为初始迭代向量赋值*/int i; srand(int)time(0);for(i=0;i=500;i+)ai=(10.0*rand()/RAND_MAX); /*生成010的随机数*/ void chuzhi2(double a,int j) /*令初始迭代向量为ei*/int i;for(i=0;i=500;i+)ai=0;aj=1;void LU

7、(double array5501) /*对矩阵A进行Doolittle分解*/ /*矩阵A转存在C5501中*/int j,k,t; /*分解结果L，U分别存在C5501的上半部与下半部*/for(k=0;k=500;k+)for(j=k;j=min2(k+2),500);j+)for(t=max3(0,k-2,j-2);t=(k-1);t+)arrayk-j+2j-=arrayk-t+2t*arrayt-j+2j;if(k500)for(j=k+1;j=min2(k+2),500);j+)for(t=max3(0,k-2,j-2);t=(k-1);t+)arrayj-k+2k-=array

8、j-t+2t*arrayt-k+2k;arrayj-k+2k=arrayj-k+2k/array2k;double fmifa(double u501,double array5501,double p) /*带原点平移的反幂法*/int i,j;double a,b,c=0;double y501;/引入平移量for(i=0;i=500;i+)array2i-=p; /先将矩阵Doolittle分解LU(array);for(;)a=0;b=0; /求k-1for(i=0;i=500;i+)a=a+ui*ui;a=sqrt(a); /求yk-1for(i=0;i=500;i+)yi=ui/a

9、; /回带过程，求解uk for(i=0;i=500;i+)ui=yi;for(i=1;i=500;i+)for(j=max2(0,(i-2);j=0;i-)for(j=i+1;j=min2(i+2),500);j+)ui-=arrayi-j+2j*uj;ui=ui/array2i;/求kfor(i=0;i=500;i+)b+=yi*ui;if(fabs(b-c)/b)=E) /*达到精度要求，迭代终止*/break; c=b;return (p+(1/b); /*直接返回距离原点P最接近的A的特征值*/主函数int main() int i; double d1,d501,ds,d,a;do

10、uble u501; double MatrixC5501;printf( 数值分析计算实习题目第一题n);printf( sy1405317 梁天骄n);/将矩阵A转存为MatrixCassignment(MatrixC);/用带原点平移的幂法求解1，501 chuzhi(u);d=mifa(u,MatrixC,0);chuzhi(u);a=mifa(u,MatrixC,d);if(d0)d1=d;d501=a;elsed501=d;d1=a; printf(1=%.12en,d1);printf(501=%.12en,d501);/用反幂法求schuzhi(u);ds=fmifa(u,Ma

11、trixC,0);printf(s=%.12en,ds);/用带原点平移的反幂法求ikfor(i=1;i=39;i+)a=d1+(i*(d501-d1)/40;assignment(MatrixC);chuzhi(u);d=fmifa(u,MatrixC,a);printf(与%02d=%+.12e最接近的特征值i%02d=%+.12en,i,a,i,d);/求A的条件数d=fabs(d1/ds);printf(A的（谱范数）条件数cond2=%.12en,d);/求detAassignment(MatrixC);LU(MatrixC);a=1;for(i=0;i=500;i+)a*=Matr

12、ixC2i;printf(行列式detA=%.12en,a);/测试不同迭代初始向量对1计算结果的影响。printf(改变迭代初始向量对max计算结果的测试如下：n); assignment(MatrixC);for(i=0;i=500;i+) chuzhi2(u,i);d1=mifa(u,MatrixC,0); printf(u%03d,max=%+e ,i,d1);if(i+1)%3)=0)printf(n);printf(Press any key to continuen);getchar();return 0;三. 程序结果：四分析初始向量选择对计算结果的影响矩阵的初始向量选择，对结

13、果的影响很大，选择不同的初始向量可能会得到的特征值。以幂法为例（反幂法原理相同），选取初始迭代向量ui=ei(i=0,1，500)；即uj=0,ji；uj=1,j=i。测试结果如下：试验结果发现只有当i取特定的一些值时才能得到正确的结果，即得到按摸最大的特征值。i取不同值时，即取不同的初始向量时，可能得到不同的特征值。这是因为以A的n个线性无关的特征向量为一组基，将初始向量线性表出时，按摸最大的特征值1对应的特征向量x1的系数1如果为0，就无法求出特征值1。如果按摸第二大的特征值2对应的特征向量x2的系数2不为0，则求出该特征值。若为0，则以此类推。五你觉得在本次试验中碰到什么印象深刻的问题？你是怎么解决的？首先拿到题目的时候，我不知道该如何设计算法，后来经过查阅教材以及例题，明白了要求出矩阵按摸最大和最小的特征值max和min，要分别使用幂法和反幂法迭代，其中min即为s，再使用带原点平移的幂法（mifa（）函数），令平移量p=max，求出对应的按摸最大的特征值，max，同时使用带原点平移的反幂法，令平移量p=k，即可求出与k最接近的特征值ik。在求矩阵的行列式时，需要用LU分解，detA就等于U所有对角线上元素的乘积。另外在编程过程中，由于我这一次使用的是C+