北航数值分析大作业一

1、数值分析大作业一1、 算法设计方案由于算法要求所有零元素都不存储，因此采用带状矩阵来存储矩阵,主要子程序为LU分解子程序、求解方程组子程序、带偏移量幂法求特征值子程序、带偏移量反幂法求特征值子程序具体算法设计思路如下：利用幂法求按模最大的特征值使用带原点平移的幂法求的按模最大的特征值，则为另外一特征值比较与的大小，较大值即为，较小者即为使用反幂法求的按模最小的特征值，利用LU分解求取每一次迭代后的 使用带偏移量的反幂法，偏移量分别为，即可求得分别与之相距最近的特征值为非奇异的是对称矩阵，因此，其中和分别为按模最大的特征值和按模最小的特征值本题中LU分解采用doolittle方法，分解后所得矩阵

2、即为对A进行初等行列式变换后所得矩阵，切变换过程中行列式值不变，因此2、 源程序#include <math.h>#include <stdio.h>#include <stdlib.h>#define N 501/矩阵为501*501的矩阵#define s 2/上半带宽#define r 2/下半带宽3个宏定义为方便LU分解及求解方程组过程/*全局变量定义*/double A5501;double u501,y501;double lambda41;/lambda0为1及最小特征值，lambda40为501及最大特征值double lambda_s,la

3、mbda_m;/按模最小(大)的特征值；double LU5501;/*子函数声明*/void init_A(double A501);/初始化A矩阵double module_value_u(double tt);/求u501的模值void init_u(double tt);/初始化迭代初始向量u0double power_method(double offset);/带原点偏移的幂法，返回值为特征值double inverse_power_method(double offset);/带原点偏移的反幂法，返回值为特征值，子函数中打印出偏移量，求得的特征值，误差，迭代次数void LU_d

4、ecomposition(double c501);/参数c为矩阵LU5501首地址，程序进行完后，保存分解后的L和U缩减后的矩阵void solve(double c501,double b,double x);/第1个参数为LU分解完的LU矩阵，第2个参数b为已知的右端值，第3个参数x为求得的解向量存储位置int max_2(int a,int b);int max_3(int a,int b,int c);int min_2(int a,int b);void main()/主程序int i;double uk;/偏移量double det;/A的行列式的值init_A(A);/初始化矩

5、阵Alambda_m=power_method(0);/求按模最大的特征值lambda40=power_method(lambda_m);/求相反方向另一个端点的特征值if(lambda_m>lambda40)/若大小反向，则交换两个元素中的值，得到1和501lambda0=lambda40;lambda40=lambda_m;elselambda0=lambda_m;lambda_s=inverse_power_method(0);det=1;for(i=0;i<N;i+)det=det*LUsi;for(i=1;i<40;i+)uk=lambda0+(lambda40-l

6、ambda0)*i/40;lambdai=inverse_power_method(uk);printf("-The results are as follows-n");printf("1=%1.11e n501=%1.11en",lambda0,lambda40);printf("s=%1.11en",lambda_s);printf("cond(A)=%1.11en",fabs(lambda_m/lambda_s);printf("detA=%1.11e n",det);for (i=1;

7、i<40;i+)printf("i%d=%1.11e ",i,lambdai);if(i% 3=0) printf("n");printf("n");void init_A(double A501)/初始化A矩阵int i,j;for(i=0;i<5;i+)if(abs(i-2)=0)for(j=0;j<501;j+)Aij=(1.64-0.024*(j+1)*sin(0.2*j+0.2)-0.64*exp(0.1/(j+1);else if(abs(i-2)=1)for(j=0;j<501;j+)Aij=0.

8、16;elsefor(j=0;j<501;j+)Aij=-0.064;void init_u(double tt)/初始化迭代初始向量,自变量为初始向量的数组int i;for(i=0;i<501;i+)tti=1;double module_value_u(double tt)/求u501的模值int i=501;double t=0;for(i=0;i<501;i+)t=t+tti*tti;return sqrt(t);double power_method(double offset)/带原点偏移的幂法，返回值为特征值double b=0,b0=0,e;double m

9、;/求得的向量模值int i,j,k=0;/k表示迭代次数init_u(u);dok+;b0=b;m=module_value_u(u);for(i=0;i<501;i+)yi=ui/m;for(i=0;i<501;i+)ui=0;for(j=max_2(0,i-s);j<min_2(i+s+1,N);j+)ui=Ai-j+2j*yj+ui;if(i=j)ui=ui-offset*yj;b=0;for(i=0;i<501;i+)b=b+yi*ui;while(fabs(e=b-b0)>1e-12);b=b+offset;printf("=%f e=%e

10、k=%d n",b,e,k);/分别为特征值，迭代误差，迭代次数return (b);void LU_decomposition(double c501)/参数c为矩阵LU5501首地址，程序进行完后，保存分解后的L和U缩减后的矩阵int j,k,t;for(k=0;k<N;k+)/求U矩阵的第k行以及L矩阵的第k列for(j=k;j<min_2(k+s,N-1)+1;j+)/for(t=max_3(0,k-r,j-s);t<k;t+)/求U矩阵的行中各个元素的循环ck-j+sj=ck-j+sj-ck-t+st*ct-j+sj;/每次计算完U的元素后才能计算L的元素

11、，所以上下两个循环不能合并在一起if(csk=0)printf("a%d%d=0,Unable to show the solution of equationsn",k,k);exit(1);if(j<min_2(k+s,N-1)+1)/因为矩阵L没有N行，因此这里要加一个判断，符合条件才能进行下面的循环语句for(t=max_3(0,j-r+1,k-s);t<k;t+)/求L矩阵的行中各个元素的循环cj-k+s+1k=cj-k+s+1k-cj-t+s+1t*ct-k+sk;cj-k+s+1k=cj-k+s+1k/csk;void solve(double c

12、501,double b,double x)/第1个参数为分解完的LU矩阵，第2个参数b为已知的右端值，第3个参数x为求得的解向量存储位置int i,t;for(i=0;i<N;i+)/此循环求出向量yxi=bi;for(t=max_2(0,i-r);t<i;t+)xi=xi-ci-t+st*xt;for(i=N-1;i>=0;i-)for(t=i+1;t<min_2(i+s,N-1)+1;t+)xi=xi-ci-t+st*xt;xi=xi/csi;double inverse_power_method(double offset)/带原点偏移的反幂法，返回值为特征值i

13、nt i,j,k=0;double b=0,b0=0,e;for(i=0;i<5;i+)if(i!=2)for(j=0;j<501;j+)LUij=Aij;elsefor(j=0;j<501;j+)LUij=Aij-offset;LU_decomposition(LU);/LU分解，只分解一次即可init_u(u);dok+;b0=b;for(i=0;i<501;i+)yi=ui/module_value_u(u);solve(LU,y,u);b=0;for(i=0;i<501;i+)b=b+yi*ui;b=1/b;b=b+offset;while(fabs(e=

14、b-b0)>1e-12);printf("offset=%1.11e =%1.11e e=%1.11e k=%d n",offset,b,e,k);/分别为特征值，迭代误差，迭代次数return b;int max_2(int a,int b)if(a<b)a=b;return a;int max_3(int a,int b,int c)if(a<b)a=b;if(a<c)a=c;return a;int min_2(int a,int b)if(a>b)a=b;return a;3、 计算结果如下所示在幂法和反幂法中偏移量、特征值、误差及计算迭代次数也在计算过程中打印出来4、 迭代过程中的现象、原因及解决方法迭代的初始向量对计算结果可能产生极为重要的影响，如果选取初始向量不恰当，则可能得到错误的结果。以最大（或最小）特征值也即按模最大特征值为例：分别取一下两组初始向量：1.2.当初始向量为时可以得到按模最大的特征向量为：当初始向量为时可以得到按模最