北航数值分析大作业一.docxVIP

数值分析—计算实习作业一学院：机械工程学院专业：材料加工工程姓名：暴一品学号： SY12071342012-10-29一、算法设计方案观察矩阵A，结构为带状，且与主对角线相邻的两个带的值b和c都是常数。从而可以用带原点平移的幂法或反幂法计算λ1和λ501。所以算法的设计方案如下：1.求按模最大的特征值，并记为max_eigenvalue，算法如下所示2.平移矩阵得到A’=A-max_eigenvalueI，再次用幂法，这次求出的A’的按模大的特征值pymax_eigenvalue就是与步骤1求出的特征值相差最大的特征值。即两者一个为最大的特征值，另一个为最小的特征值。3.根据max_eigenvalue和pymax_eigenvalue的正负性，直接确定λ1，和λ501。4.对原矩阵A用反幂法，求出其按模最小的特征值，记为s_eigenvalue，此即λs。在反幂法的求解过程中，每迭代一次都要求满足解线性方程组Auk=yk-1。本题中矩阵A上半带宽为2，下半带宽也为2 。故选择采用三角分解法求解方程组：先将原矩阵改写成5行501列的矩阵C（不存储A的0元素）A的带内元素aij=c中的元素ci-j+3。再对C矩阵做LU分解。对于k=1,2，…，n，执行 [j=k,k+1,…,min(k+2,n)][i=k+1,k+2,…,min(k+2,n);kn]求解Lx=b，Uuk=x（数组b先是存放原方程组右端向量yk-1，后来存放中间向量x） （i=2,3,…,n） (i=n-1,n-2,…,1)5.对k=1，2，……39执行：先根据题中给出的公式算出μk，再将矩阵平移A”=A-μk，对矩阵A”运用反幂法（线性方程组的解法同上），就可以求出与μk最接近的特征值λik，保存在数组py_eigenvalue中。6.根据公式 其中λmax和λmin分别是矩阵A的按模最大和按模最小的特征值。将max_eigenvalue和s_eigenvalue分别代入λmax和λmin，就可以求出cond(A) 2。7.对矩阵A做三角分解，A=LU。则有，其中u（i，i）表示上三角阵U的第i行i列的对角线上元素。二、源程序代码#include iostream#include cmath#include iomanipusing namespace std;const double JD=1e-12;const double B=0.16;const double C=-0.064;double a[501],c[5][501]={0};void sanjiaofj(){int ci,neici,i,j,k,t;double sum;for (i=0;i501;i++) c[2][i]=a[i];for (i=2;i501;i++) c[0][i]=C;for (i=1;i501;i++) c[1][i]=B;for (i=0;i500;i++) c[3][i]=B;for (i=0;i499;i++) c[4][i]=C;for (k=0;k501;k++){ci=k+2;if (ci500) ci=500;for (j=k;j=ci;j++){neici=0;if ((k-2)neici) neici=k-2;if ((j-2)neici) neici=j-2;sum=0;for (t=neici;t=k-1;t++) sum+=c[k-t+2][t]*c[t-j+2][j];c[k-j+2][j]-=sum;}if (k500){for (i=k+1;i=ci;i++){neici=0;if ((k-2)neici) neici=k-2;if ((i-2)neici) neici=i-2;sum=0;for (t=neici;t=k-1;t++) sum+=c[i-t+2][t]*c[t-k+2][k];c[i-k+2][k]=(c[i-k+2][k]-sum)/c[2][k];}}}}void mifa(double *bold,double pingyi){int k=0,i,r;double u[501]={0},y[501];double umax,umo,bi=1,b;for (i=0;i501;i++) u[i]=1;//对初始向量赋值if (pingyi!=0) for (i=0;i501;i++) a[i]-=pingyi;while ((k==0)||(k==1)||(biJD)){r=0;umax=u[r];for (i=1;i501;i++) if (fabs(u[i])fabs(umax)){umax=u[i];r=i;}umo=fabs(umax);for (i=0;i501;i++) y[i]=u[i]/umo;for