計算實習報告1

數(shù)值分析（B）大作業(yè)（一）姓名：學號：電話：1、算法設計：①求、和的值：：表示矩陣的按模最小特征值，為求得直接對待求矩陣A應用反冪法即可。、：若矩陣A的特征值滿足關系且，要求、及時，可按如下方法求解：對矩陣A用冪法，求得按模最大的特征值。按平移量對矩陣A進行原點平移得矩陣，對矩陣B用反冪法求得B的按模最小特征值。則：，即為所求。②求和A的與數(shù)最接近的特征值（k=0，1，…39）：求矩陣A的特征值中與P最接近的特征值的大小，采用原點平移的方法：先求矩陣B=A-PI對應的按模最小特征值，則+P即為矩陣A與P最接近的特征值。在本次計算實習中則是先求平移矩陣，對該矩陣應用反冪法求得，則與最接近的A的特征值為：重復以上過程39次即可求得（k=0，1，…39）的值。③求A的（譜范數(shù)）條件數(shù)和行列式：在（1）中用反冪法求矩陣A的按模最小特征值時，要用到Doolittle分解方法，在Doolittle分解完成后得到的兩個矩陣分別為L和U，則A的行列式可由U陣求出，即：det(A)=det(U)。求得det(A)不為0，因此A為非奇異的實對稱矩陣，則：，和分別為模最大特征值與模最小特征值。2、程序源代碼：#include"Stdio.h"#include"Conio.h"#include"math.h"http://****************************************************************************////在存儲帶狀矩陣時，下面的幾個量在程序中反復用到，為方便編程故把它們定義成宏.////M：轉換后的矩陣的行數(shù)，M=R+S+1。////N：轉換后的矩陣的列數(shù)，與原矩陣列數(shù)相等。////S：上半帶寬。////R：下半帶寬。////Epsilon：用來指定求解精度。////****************************************************************************//#defineM5#defineN501#defineR2#defineS2#defineEpsilon0.000000000001voidCreat_MatrixA(doublearray[M+1][N+1]);voidLoad_MatrixA(doublearrayA[M+1][N+1],doublearrayB[M+1][N+1]);voidLoad_vectoru(doubleu[N+1]);doubleMifa(doubleu[N+1],doublearray[M+1][N+1]);doubleFanmifa(doubleu[N+1],doublearrayA[N+1][N+1]);voidSolution_Yushu(doublearray[N+1][N+1],doubleu[N+1],doublenamda1,doublenamda501);voidDoolittle_Dai(doublearray[M+1][N+1]);voidBack_Doolittle_Dai(doublearray[M+1][N+1],doubley[N+1],doublex[N+1]);doublemax(doublea,doubleb);doublemin(doublea,doubleb);main(void){//****************************************************************************////定義主程序中用到的變量////****************************************************************************////MatrixA：用來存儲源矩陣A。////arrayA：用反冪法求矩陣的按模最小特征值時，矩陣的數(shù)據(jù)會被更改。////因此實際計算中，使用源矩陣A的拷貝arrayA。////u：用來存放迭代向量，初始化后的u里存儲的是初始迭代向量。////u_k：u_k=λ(1)+k*(λ(501)-λ(1))/40，用來存儲A的與數(shù)。////cond(A)：用來存儲A的條件數(shù)cond(A)=λ(max)/λ(s)。////det_A：用來存儲A的行列式的值。////****************************************************************************//doubleMatrixA[M+1][N+1],arrayA[M+1][N+1],u[N+1];doublenamdas=0,namda1=0,namda2=0,namda501=0;doublecond_A,det_A=1;inti;/*調(diào)用函數(shù)Creat_MatrixA生成帶狀矩陣A*/Creat_MatrixA(MatrixA);//****************************************************************************////求矩陣A的最小、最大和按模最小的特征值////****************************************************************************//printf("/***********************************************************/\n");printf("矩陣A的兩端特征值及按模最小特征值分別為：\n");/*求λ(s)：加載矩陣A，對A用反冪法*/Load_MatrixA(MatrixA,arrayA);Load_vectoru(u);namdas=Fanmifa(u,arrayA);/*利用反冪法進行Doolittle分解得到的U陣求矩陣A的行列式的值：det(A)=det(U)*/for(i=1;i<=N;i++) det_A*=arrayA[S+1][i];/*重新加載矩陣A，對A用冪法求得A的按模最大特征值λ(m1)*/Load_MatrixA(MatrixA,arrayA);Load_vectoru(u);namda1=Mifa(u,arrayA);/*計算B=A+λ(m1)I，并對B用反冪法求得λ(m2)*/for(i=1;i<=N;i++) arrayA[S+1][i]+=namda1;/*調(diào)用反冪法，求變換矩陣的按模最小特征值*/namda2=Fanmifa(u,arrayA);namda501=namda2-namda1;printf("λ(1)=%.12e\nλ(501)=%.12e\nλ(s)=%.12e\n",min(namda1,namda501),max(namda1,namda501),namdas);//****************************************************************************////求和A的與數(shù)最接近的特征值的大小////****************************************************************************//printf("/***********************************************************/\n");printf("和A的與數(shù)最接近的特征值λ(i_k)分別為：\n");/*調(diào)用求和A的與數(shù)最接近的特征值的大小子函數(shù)*/Solution_Yushu(MatrixA,u,namda1,namda501);//****************************************************************************////求矩陣A的條件數(shù)cond(A)和行列式的值det(A)////****************************************************************************//printf("/***********************************************************/\n");cond_A=fabs(max(namda1,namda501)/namdas);printf("矩陣A的條件數(shù)：cond(A)=%.12e\n",cond_A);printf("A的行列式的值為：det(A)=%.12e\n",det_A);getch();return0;}/*創(chuàng)建源矩陣A,按帶狀特征存儲*/voidCreat_MatrixA(doublearray[M+1][N+1]){inti,j;for(i=0;i<=M;i++)for(j=0;j<=N;j++)array[i][j]=0;for(j=3;j<=N;j++)array[1][j]=-0.064;for(j=2;j<=N;j++)array[2][j]=0.16;for(j=1;j<=N;j++)array[3][j]=(1.64-0.024*j)*sin(0.2*j)-0.64*exp(0.1/j);for(j=1;j<=N-1;j++)array[4][j]=0.16;for(j=1;j<=N-2;j++)array[5][j]=-0.064;/*輸出矩陣A各元素，測試創(chuàng)建矩陣是否正確*//*for(i=1;i<=M;i++)for(j=1;j<=N;j++)printf("%lf",array[i][j]);*/}/*生成矩陣A的拷貝，避免計算過程中對源矩陣A的修改*/voidLoad_MatrixA(doublearrayA[M+1][N+1],doublearrayB[M+1][N+1]){inti,j;for(i=0;i<=M;i++)for(j=0;j<=N;j++)arrayB[i][j]=arrayA[i][j];/*測試輸出矩陣A的拷貝，看加載是否正確*//*for(i=1;i<=M;i++)for(j=1;j<=N;j++)printf("%f",arrayB[i][j]);*/}/*給出迭代初始向量u0*/voidLoad_vectoru(doubleu[N+1]){inti;for(i=0;i<=200;i++)u[i]=1;for(i=201;i<=N;i++)u[i]=0;/*測試輸出迭代初始向量u*//*for(i=1;i<=N;i++)printf("%f",u[i]);*/}/*冪法*/doubleMifa(doubleu[N+1],doublearray[M+1][N+1]){/*定義變量：ita表示η，beta表示β(也即是按模最大的λ值)*/doubleita=0,beta=0,y[N+1],temp=0,namda0=-1;inti,j,flag;flag=0;while(fabs(beta-namda0)>Epsilon||flag<900){ namda0=beta; ita=beta=0;for(i=1;i<=N;i++)ita+=u[i]*u[i];ita=sqrt(ita);/*計算向量y：y_k-1=u_k-1/η_k-1*/for(i=1;i<=N;i++)y[i]=u[i]/ita;/*計算向量u：u_k=A(y_k-1)*//*經(jīng)過轉轉換后的帶狀矩陣A，在與向量的乘法中按元素位置，在左上角、中間和右下角三個地方分別表現(xiàn)出不同的下標變化規(guī)律。*/for(i=1;i<=N;i++) {if(i<=2) {for(j=1;j<=i+2;j++)temp+=array[i-j+S+1][j]*y[j];u[i]=temp;temp=0; }if(i>=3&&i<=499) {for(j=i-2;j<=i+2;j++)temp+=array[i-j+S+1][j]*y[j];u[i]=temp;temp=0; }if(i>=500) {for(j=i-2;j<=N;j++)temp+=array[i-j+S+1][j]*y[j];u[i]=temp;temp=0; }}for(i=1;i<=N;i++)beta+=y[i]*u[i]; flag++;}/*測試輸出用冪法求得的λ*//*printf("λ=%.12e\n",beta);*/returnbeta;}/*反冪法*/doubleFanmifa(doubleu[N+1],doublearray[N+1][N+1]){doubleita=0,beta=0,y[N+1],namda0=-1;/*在這里ita表示η，beta表示β*/inti,flag;flag=0;while(fabs(beta-namda0)>Epsilon||flag<900){namda0=beta; ita=beta=0;for(i=1;i<=N;i++)ita+=u[i]*u[i];ita=sqrt(ita);for(i=1;i<=N;i++)y[i]=u[i]/ita;/*y_k-1=u_k-1/η_k-1*//*用flag作為是否進行Doolittle分解的標志。當flag=0時，是第一次進入反冪法求解子函數(shù)，要進行Doolittle分解。Doolittle分解完成后隨即進行回代求解：Au_k=y_k-1*/if(flag==0)Doolittle_Dai(array);Back_Doolittle_Dai(array,y,u);for(i=1;i<=N;i++)beta+=y[i]*u[i]; flag++;}/*測試用反冪法輸出求得的λ*//*printf("λ=%.12e\n",1/beta);*/return1/beta;}/*帶狀矩陣的Doolittle分解*/voidDoolittle_Dai(doublearray[M+1][N+1]){intk,i,j,t;doubletemp;for(k=1;k<=N;k++){for(j=k;j<=min(k+S,N);j++) {temp=0;for(t=max(max(1,k-R),j-S);t<=k-1;t++)temp+=array[k-t+S+1][t]*array[t-j+S+1][j];array[k-j+S+1][j]=array[k-j+S+1][j]-temp; }for(i=k+1;i<=min(k+R,N);i++) {temp=0;for(t=max(max(1,i-R),k-S);t<=k-1;t++)temp+=array[i-t+S+1][t]*array[t-k+S+1][k];array[i-k+S+1][k]=(array[i-k+S+1][k]-temp)/array[S+1][k]; }}}/*帶狀矩陣Doolittle分解的回代過程*//*array[M+1][N+1]：用來傳遞經(jīng)過Dolittle分解后的矩陣A。b[N+1]：用作向量y的拷貝，避免在回代過程中造成對y的修改。*/voidBack_Doolittle_Dai(doublearray[M+1][N+1],doubley[N+1],doublex[N+1]){inti,t;doubletemp,b[N+1];for(i=0;i<=N;i++)b[i]=y[i];for(i=2;i<=N;i++){temp=0;for(t=max(1,i-R);t<=i-1;t++)temp+=array[i-t+S+1][t]*b[t];b[i]=b[i]-temp;}x[N]=b[N]/array[S+1][N];for(i=N-1;i>=1;i--){temp=0;for(t=i+1;t<=min(i+S,N);t++)temp+=array[i-t+S+1][t]*x[t];x[i]=(b[i]-temp)/array[S+1][i];}}/*求和A的與數(shù)最接近的特征值的大小*//*A的與數(shù)uk=λ(1)+k*(λ(501)-λ(1))/40。求矩陣A與uk(k=1,2,……39)最接近的特征值，先求A的原點平移矩陣A-ukI,再求A-ukI按模最小特征值，最后將該特征值加上反向平移量uk。*/voidSolution_Yushu(doubleMatrixA[M+1][N+1],doubleu[N+1],doublenamda1,doublenamda501){inti,k;doubleu_k,namda,array[M+1][N+1];for(k=1;k<=39;k++){/*平移量不同，反冪法的對象A-ukI就不一樣，因此要重新加載A陣*/Load_MatrixA(MatrixA,array);Load_vectoru(u);u_k=namda1+(namda501-namda1)*k/40;for(i=1;i<=N;i++)array[S+1][i]-=u_k;namda=Fanmifa(u,array);namda=u_k+namda;printf("λ(i_%d)=%.12e\n",k,namda);}}/*求兩個數(shù)中的最大、最小值*//*max(a,b)：返回a、b中的較大值。min(a,b)：返回a、b中的較小值。*/doublemin(doublea,doubleb){if(a<=b)returna;elsereturnb;}doublemax(doublea,doubleb){if(a>=b)returna;elsereturnb;}3、程序運行結果：/***********************************************************/矩陣A的兩端特征值及按模最小特征值分別為：λ(1)=-1.070011361514e+001λ(501)=9.724634099672e+000λ(s)=-5.557910794214e-003/***********************************************************/和A的與數(shù)最接近的特征值λ(i_k)分別為：λ(i_1)=-1.018293403315e+001λ(i_2)=-9.581110544848e+000λ(i_3)=-9.172672423928e+000λ(i_4)=-8.652284007898e+000λ(i_5)=-8.093483808675e+000λ(i_6)=-7.659405407692e+000λ(i_7)=-7.119684648691e+000λ(i_8)=-6.611764339397e+000λ(i_9)=-6.066103226595e+000λ(i_10)=-5.585101052628e+000λ(i_11)=-5.114083529812e+000λ(i_12)=-4.578872176865e+000λ(i_13)=-4.096470926260e+000λ(i_14)=-3.554211215751e+000λ(i_15)=-3.041090018133e+000λ(i_16)=-2.533970311130e+000λ(i_17)=-2.003230769564e+000λ(i_18)=-1.5035576