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

1、數(shù)值分析大作業(yè)一1、 算法設(shè)計(jì)方案由于算法要求所有零元素都不存儲(chǔ)，因此采用帶狀矩陣來(lái)存儲(chǔ)矩陣,主要子程序?yàn)長(zhǎng)U分解子程序、求解方程組子程序、帶偏移量?jī)绶ㄇ筇卣髦底映绦颉屏糠磧绶ㄇ筇卣髦底映绦蚓唧w算法設(shè)計(jì)思路如下：利用冪法求按模最大的特征值使用帶原點(diǎn)平移的冪法求的按模最大的特征值，則為另外一特征值比較與的大小，較大值即為，較小者即為使用反冪法求的按模最小的特征值，利用LU分解求取每一次迭代后的 使用帶偏移量的反冪法，偏移量分別為，即可求得分別與之相距最近的特征值為非奇異的是對(duì)稱矩陣，因此，其中和分別為按模最大的特征值和按模最小的特征值本題中LU分解采用doolittle方法，分解后所得矩陣

2、即為對(duì)A進(jìn)行初等行列式變換后所得矩陣，切變換過(guò)程中行列式值不變，因此2、 源程序#include <math.h>#include <stdio.h>#include <stdlib.h>#define N 501/矩陣為501*501的矩陣#define s 2/上半帶寬#define r 2/下半帶寬3個(gè)宏定義為方便LU分解及求解方程組過(guò)程/*全局變量定義*/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);/帶原點(diǎn)偏移的冪法，返回值為特征值double inverse_power_method(double offset);/帶原點(diǎn)偏移的反冪法，返回值為特征值，子函數(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個(gè)參數(shù)為L(zhǎng)U分解完的LU矩陣，第2個(gè)參數(shù)b為已知的右端值，第3個(gè)參數(shù)x為求得的解向量存儲(chǔ)位置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);/求相反方向另一個(gè)端點(diǎn)的特征值if(lambda_m>lambda40)/若大小反向，則交換兩個(gè)元素中的值，得到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)/帶原點(diǎn)偏移的冪法，返回值為特征值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矩陣的行中各個(gè)元素的循環(huán)ck-j+sj=ck-j+sj-ck-t+st*ct-j+sj;/每次計(jì)算完U的元素后才能計(jì)算L的元素

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