北航数值分析大作业一.docx

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

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

3、turn 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=

4、0; j+; i+; /第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) /*带

5、原点平移的幂法*/ int 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-1 for(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-

6、p*yi; /*引入平移量*/ /求k for(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; f

7、or(i=0;i=500;i+) ai=0; aj=1;void LU(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=ma

8、x3(0,k-2,j-2);t=(k-1);t+) arrayj-k+2k-=arrayj-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-1 for(i=0;i=50

9、0;i+) a=a+ui*ui; a=sqrt(a); /求yk-1 for(i=0;i=500;i+) yi=ui/a; /回带过程，求解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; /求k for(i=0;i=500;i+) b+=yi*ui; if(fabs(b-c)/b)=E) /*达到精度要求，迭代终止*/ break; c=b; return (p+(1

10、/b); /*直接返回距离原点P最接近的A的特征值*/主函数main() int i; double d1,d501,ds,d,a; double u501; double MatrixC5501; printf( 数值分析计算实习题目第一题n); printf( SY1103120 朱舜杰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; else d501=d;

11、d1=a; printf(1=%.12en,d1); printf(501=%.12en,d501);/用反幂法求s chuzhi(u); ds=fmifa(u,MatrixC,0); printf(s=%.12en,ds);/用带原点平移的反幂法求ik for(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

12、(A的（谱范数）条件数cond2=%.12en,d);/求detA assignment(MatrixC); LU(MatrixC); a=1; for(i=0;i=500;i+) a*=MatrixC2i; 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)=

13、0) printf(n); printf(Press any key to continuen); getchar();三. 程序结果：四分析初始向量选择对计算结果的影响 矩阵的初始向量选择，对结果的影响很大，选择不同的初始向量可能会得到的特征值。以幂法为例（反幂法原理相同），选取初始迭代向量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，则以此类推。