北航數(shù)值分析大作業(yè)一

1、數(shù)值分析大作業(yè)一1、 算法設(shè)計方案由于算法要求所有零元素都不存儲，因此采用帶狀矩陣來存儲矩陣,主要子程序為LU分解子程序、求解方程組子程序、帶偏移量冪法求特征值子程序、帶偏移量反冪法求特征值子程序具體算法設(shè)計思路如下：利用冪法求按模最大的特征值使用帶原點平移的冪法求的按模最大的特征值，則為另外一特征值比較與的大小，較大值即為，較小者即為使用反冪法求的按模最小的特征值，利用LU分解求取每一次迭代后的 使用帶偏移量的反冪法，偏移量分別為，即可求得分別與之相距最近的特征值為非奇異的是對稱矩陣，因此，其中和分別為按模最大的特征值和按模最小的特征值本題中LU分解采用doolittle方法，分解后所得矩陣

2、即為對A進(jìn)行初等行列式變換后所得矩陣，切變換過程中行列式值不變，因此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;/*子函數(shù)聲明*/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);/帶原點偏移的反冪法，返回值為特征值，子函數(shù)中打印出偏移量，求得的特征值，誤差，迭代次數(shù)void LU_d

4、ecomposition(double c501);/參數(shù)c為矩陣LU5501首地址，程序進(jìn)行完后，保存分解后的L和U縮減后的矩陣void solve(double c501,double b,double x);/第1個參數(shù)為LU分解完的LU矩陣，第2個參數(shù)b為已知的右端值，第3個參數(shù)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)/初始化迭代初始向量,自變量為初始向量的數(shù)組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表示迭代次數(shù)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);/分別為特征值，迭代誤差，迭代次數(shù)return (b);void LU_decomposition(double c501)/參數(shù)c為矩陣LU5501首地址，程序進(jìn)行完后，保存分解后的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矩陣的行中各個元素的循環(huán)ck-j+sj=ck-j+sj-ck-t+st*ct-j+sj;/每次計算完U的元素后才能計算L的元素

11、，所以上下兩個循環(huán)不能合并在一起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行，因此這里要加一個判斷，符合條件才能進(jìn)行下面的循環(huán)語句for(t=max_3(0,j-r+1,k-s);t<k;t+)/求L矩陣的行中各個元素的循環(huán)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個參數(shù)為分解完的LU矩陣，第2個參數(shù)b為已知的右端值，第3個參數(shù)x為求得的解向量存儲位置int i,t;for(i=0;i<N;i+)/此循環(huán)求出向量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);/分別為特征值，迭代誤差，迭代次數(shù)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、 計算結(jié)果如下所示在冪法和反冪法中偏移量、特征值、誤差及計算迭代次數(shù)也在計算過程中打印出來4、 迭代過程中的現(xiàn)象、原因及解決方法迭代的初始向量對計算結(jié)果可能產(chǎn)生極為重要的影響，如果選取初始向量不恰當(dāng)，則可能得到錯誤的結(jié)果。以最大（或最?。┨卣髦狄布窗茨Ｗ畲筇卣髦禐槔悍謩e取一下兩組初始向量：1.2.當(dāng)初始向量為時可以得到按模最大的特征向量為：當(dāng)初始向量為時可以得到按模最