1、數(shù)值分析計算實習(xí)題目第一題：1. 算法設(shè)計方案(1) 1， 501 和 s的值。1)首先通過冪法求出按模最大的特征值t1，然后根據(jù) t1 進行原點平移求出另一特征值t2，比較兩值大小， 數(shù)值小的為所求最小特征值 1，數(shù)值大的為是所求最大特征值501。2)使用反冪法求 s，其中需要解線性方程組。因為A 為帶狀線性方程組，此處采用 LU分解法解帶狀方程組。( 2)與 k = 1+k 501 1 最接近的特征值 ik 。k 1 40 ik通過帶有原點平移的反冪法求出與數(shù) k 最接近的特征值 ik 。(3) cond(A) 2和 detA 。1) cond(A)2 = 1 ，其中 1和 n 分別是按模

2、最大和最小特征值。2)利用步驟( 1)中分解矩陣 A 得出的 LU矩陣， L為單位下三角陣 ,U 為上三角陣，其 中 U 矩陣的主對角線元素之積即為 detA 。由于 A 的元素零元素較多，為節(jié)省儲存量，將 A 的元素存為 6501的數(shù)組中，程序中 采用 get_an_element() 函數(shù)來從小數(shù)組中取出 A 中的元素。2. 全部源程序#include #include void init_a();/ 初始化 Adouble get_an_element(int,int);/ 取 A 中的元素函數(shù)double powermethod(double);/ 原點平移的冪法double inve

3、rsepowermethod(double);/ 原點平移的反冪法int presolve(double);/ 三角 LU 分解int solve(double ,double );/ 解方程組int max(int,int);int min(int,int);double (*u)502=new double502502;/ 上三角 U 數(shù)組double (*l)502=new double502502;/ 單位下三角 L 數(shù)組double a6502;/ 矩陣 Aint main()int i,k;double lambdat1,lambdat2,lambda1,lambda501,lam

4、bdas,mu40,det;double lambda40; init_a();/ 初始化 A lambdat1=powermethod(0); lambdat2=powermethod(lambdat1); lambda1=lambdat1lambdat2?lambdat1:lambdat2; presolve(0);lambdas=inversepowermethod(0);det=1; for(i=1;i=501;i+) det=det*uii;for (k=1;k=39;k+) muk=lambda1+k*(lambda501-lambda1)/40; presolve(muk); l

5、ambdak=inversepowermethod(muk);printf( 所有特征值如下 n);printf( =%1.11e =%1.11en,lambda1,lambda501); printf( s=%1.n11,leambdas);printf(cond(A)=%1.11en,fabs(lambdat1/lambdas); printf(detA=%1.11e n,det);for (k=1;k=39;k+)printf( i%d=%1.11e ,k,lambdak);if(k % 3=0) prinntf();delete u;delete l;/ 釋放堆內(nèi)存return 0;v

6、oid init_a()/ 初始化 Aint i;for (i=3;i=501;i+) a1i=a5502-i=-0.064;for (i=2;i=501;i+) a2i=a4502-i=0.16;for (i=1;i=501;i+) a3i=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i); double get_an_element(int i,int j)/ 從 A 中節(jié)省存儲量的提取元素方法 if (fabs(i-j)=2) return ai-j+3j;else return 0;double powermethod(double offset)/

7、冪法int i,x1;double u502,y502;double beta=0,prebeta=-1000,yita=0;for (i=1;i=501;i+)ui=1,yi=0;/ 設(shè)置初始向量 ufor (int k=1;k=10000;k+)yita=0;for (i=1;i=501;i+) yita=sqrt(yita*yita+ui*ui);for (i=1;i=501;i+) yi=ui/yita;for (x1=1;x1=501;x1+)ux1=0;for (int x2=1;x2=501;x2+)ux1=ux1+(x1=x2)?(get_an_element(x1,x2)-o

8、ffset):get_an_element(x1,x2)*yx2; prebeta=beta;beta=0;for (i=1;i=501;i+) beta=beta+ yi*ui;if (fabs(prebeta-beta)/beta)=1e-12) printf(offset=%f lambda=%f k=%dn,offset,(beta+offset),fabs(prebeta-beta)/beta),k);break;/ 輸出中間過程 ,包括偏移量 ,誤差 ,迭代次數(shù)return (beta+offset);double inversepowermethod(double offset)

9、/ 反冪法int i;double u502,y502;double beta=0,prebeta=0,yita=0;for (i=1;i=501;i+)ui=1,yi=0; / 設(shè)置初始向量 ufor (int k=1;k=10000;k+)yita=0;for (i=1;i=501;i+) yita=sqrt(yita*yita+ui*ui);for (i=1;i=501;i+) yi=ui/yita;solve(u,y); prebeta=beta;beta=0;for (i=1;i=501;i+) beta=beta+ yi*ui;beta=1/beta;if (fabs(prebet

10、a-beta)/beta)=1e-12) printf(offset=%f lambda=%f k=%dn,offset,(beta+offset),fabs(prebeta-beta)/beta),k);break;/ 輸出中間過程 ,包括偏移量 ,誤差 ,迭代次數(shù)return (beta+offset);err=%eerr=%eint presolve(double offset)/ 三角 LU 分解int i,k,j,t;double sum;for (k=1;k=501;k+)for (j=1;j=501;j+) ukj=lkj=0; if (k=j) lkj=1; / 初始化 LU

11、矩陣for (k=1;k=501;k+)for (j=k;j=min(k+2,501);j+)sum=0;for (t=max(1,max(k-2,j-2) ; t=(k-1) ; t+) sum=sum+lkt*utj;ukj=(k=j)?(get_an_element(k,j)-offset):get_an_element(k,j)-sum; if (k=501) continue;for (i=k+1;i=min(k+2,501);i+)sum=0;for (t=max(1,max(i-2,k-2);t=(k-1);t+) sum=sum+lit*utk;lik=(i=k)?(get_a

12、n_element(i,k)-offset):get_an_element(i,k)-sum)/ukk; return 0;int solve(double x,double b)/ 解方程組int i,t;double y502;double sum;y1=b1;for (i=2;i=501;i+)sum=0;for (t=max(1,i-2);t=1;i-)sum=0;for (t=i+1;ty?x:y);int min(int x,int y)return (xy?x:y);3. 計算結(jié)果 結(jié)果如下圖所示部分中間結(jié)果：給出了偏移量 (offset) ，誤差 (err) ，迭代次數(shù) (k)

13、4. 討論迭代初始向量的選取對計算結(jié)果的影響,并說明原因使用 ui=1( i=1,2,.,501)作為初始向量進行迭代， 可得出以上結(jié)果。 經(jīng)過 Mathematica 計算驗證結(jié)果正確。現(xiàn)修改初始向量 u1=1 ，ui=0，(i=2,3,.,501)。得出結(jié)果此結(jié)果與正確結(jié)果相差較多。令初始向量 um=1 ，un=0 ，(m=1,2,.,250 n=251,252,501) ，得出結(jié)果：此結(jié)果也與正確結(jié)果相差較多。但與上次結(jié)果相比，更加靠近準(zhǔn)確值一些。再增加初始向量 u 中等于 1 的元素個數(shù)，可以發(fā)現(xiàn)其結(jié)果更加靠近準(zhǔn)確值。經(jīng)驗證， 只有當(dāng)不為 0 元素的個數(shù)達到比較高的一個值時， 才能得到精確的收斂結(jié)果， 與元素的絕對 值大小無關(guān)。分析算法，設(shè) 1為按模最大特征值，觀察式子u0 1x1