计算实习报告1

数值分析（B）大作业（一）姓名：学号：电话：1、算法设计：①求、和的值：：表示矩阵的按模最小特征值，为求得直接对待求矩阵A应用反幂法即可。、：若矩阵A的特征值满足关系且，要求、及时，可按如下方法求解：对矩阵A用幂法，求得按模最大的特征值。按平移量对矩阵A进行原点平移得矩阵，对矩阵B用反幂法求得B的按模最小特征值。则：，即为所求。②求和A的与数最接近的特征值（k=0，1，…39）：求矩阵A的特征值中与P最接近的特征值的大小，采用原点平移的方法：先求矩阵B=A-PI对应的按模最小特征值，则+P即为矩阵A与P最接近的特征值。在本次计算实习中则是先求平移矩阵，对该矩阵应用反幂法求得，则与最接近的A的特征值为：重复以上过程39次即可求得（k=0，1，…39）的值。③求A的（谱范数）条件数和行列式：在（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"//****************************************************************************////在存储带状矩阵时，下面的几个量在程序中反复用到，为方便编程故把它们定义成宏.////M：转换后的矩阵的行数，M=R+S+1。////N：转换后的矩阵的列数，与原矩阵列数相等。////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：用反幂法求矩阵的按模最小特征值时，矩阵的数据会被更改。////因此实际计算中，使用源矩阵A的拷贝arrayA。////u：用来存放迭代向量，初始化后的u里存储的是初始迭代向量。////u_k：u_k=λ(1)+k*(λ(501)-λ(1))/40，用来存储A的与数。////cond(A)：用来存储A的条件数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;/*调用函数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;/*调用反幂法，求变换矩阵的按模最小特征值*/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的与数最接近的特征值的大小////****************************************************************************//printf("/***********************************************************/\n");printf("和A的与数最接近的特征值λ(i_k)分别为：\n");/*调用求和A的与数最接近的特征值的大小子函数*/Solution_Yushu(MatrixA,u,namda1,namda501);//****************************************************************************////求矩阵A的条件数cond(A)和行列式的值det(A)////****************************************************************************//printf("/***********************************************************/\n");cond_A=fabs(max(namda1,namda501)/namdas);printf("矩阵A的条件数：cond(A)=%.12e\n",cond_A);printf("A的行列式的值为：det(A)=%.12e\n",det_A);getch();return0;}/*创建源矩阵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各元素，测试创建矩阵是否正确*//*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)*//*经过转转换后的带状矩阵A，在与向量的乘法中按元素位置，在左上角、中间和右下角三个地方分别表现出不同的下标变化规律。*/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时，是第一次进入反幂法求解子函数，要进行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]：用来传递经过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的与数最接近的特征值的大小*//*A的与数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);}}/*求两个数中的最大、最小值*//*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的与数最接近的特征值λ(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