北航数值分析大作业一

《北航数值分析大作业一》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

1、数值分析—计算实习作业一学院：机械工程学院专业：材料加工工程姓名：暴一品学号：SY12071342012-10-29一、算法设计方案观察矩阵A，结构为带状，且与主对角线相邻的两个带的值b和c都是常数。从而可以用带原点平移的幂法或反幂法计算λ1和λ501。所以算法的设计方案如下：1.求按模最大的特征值，并记为max_eigenvalue，算法如下所示2.平移矩阵得到A’=A-max_eigenvalueI，再次用幂法，这次求出的A’的按模大的特征值pymax_eigenvalue就是与步骤1求出的特征值相差最大的特征值。即两者一个为最大的特征值，另一个为最小的特征值。3.根据max_ei

2、genvalue和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);k

3、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的第

4、i行i列的对角线上元素。二、源程序代码#include#include#includeusingnamespacestd;constdoubleJD=1e-12;constdoubleB=0.16;constdoubleC=-0.064;doublea[501],c[5][501]={0};voidsanjiaofj(){intci,neici,i,j,k,t;doublesum;for(i=0;i<501;i++)c[2][i]=a[i];for(i=2;i<501;i++)c[0][i]=C;for(i=1;i<501;i++)c

5、[1][i]=B;for(i=0;i<500;i++)c[3][i]=B;for(i=0;i<499;i++)c[4][i]=C;for(k=0;k<501;k++){ci=k+2;if(ci>500)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(k<500){for(i=k+1;i<=ci;i++){neic

6、i=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];}}}}voidmifa(double*bold,doublepingyi){intk=0,i,r;doubleu[501]={0},y[501];doubleumax,umo,bi=1,b;for(i=0;i<501;i++)u[i]=1;//对初始向量赋值if(pingyi!=0)for

7、(i=0;i<501;i++)a[i]-=pingyi;while((k==0)

8、

9、(k==1)

10、

11、(bi>JD)){r=0;umax=u[r];for(i=1;i<501;i++)if(fabs(u[i])>fabs(umax)){umax=u[i];r=i;}umo=fabs(umax);for(i=0;i<501;i++)y[i]=u[i]/umo;for(i=0;i<501;i++){u[i]=0;if(i==0)u[i]=a[i]*y