北航研究生数值分析编程大作业1
北航数值分析大作业 第一题 幂法与反幂法
数 值 分 析（B ） 大 作 业（一）姓名： 学号： 电话：1、算法设计：①求1λ、501λ和s λ的值：s λ：s λ表示矩阵的按模最小特征值，为求得s λ直接对待求矩阵A 应用反幂法即可。
1λ、501λ：若矩阵A 的特征值满足关系 1n λλ<<且1n λλ≠，要求1λ、及501λ时，可按如下方法求解： a ． 对矩阵A 用幂法，求得按模最大的特征值1m λ。
b ． 按平移量1m λ对矩阵A 进行原点平移得矩阵1m BA I λ=+，对矩阵B 用反幂法求得B 的按模最小特征值2m λ。
c ． 321m m m λλλ=-则：113min(,)m m λλλ=，13max(,)n m m λλλ=即为所求。
②求和A 的与数5011140k k λλμλ-=+最接近的特征值ik λ（k=0，1，…39）：求矩阵A 的特征值中与P 最接近的特征值的大小，采用原点平移的方法：先求矩阵 B=A-PI 对应的按模最小特征值k β，则k β+P 即为矩阵A 与P 最接近的特征值。
在本次计算实习中则是先求平移矩阵k B A I μ=-，对该矩阵应用反幂法求得s λ，则与k μ最接近的A 的特征值为：s P λ+重复以上过程39次即可求得ik λ（k=0，1，…39）的值。
③求A 的（谱范数）条件数2cond()A 和行列式det A ：在（1）中用反幂法求矩阵A 的按模最小特征值时，要用到Doolittle 分解方法，在Doolittle 分解完成后得到的两个矩阵分别为L 和U ，则A 的行列式可由U 阵求出，即：det(A)=det(U)。
求得det(A)不为0，因此A 为非奇异的实对称矩阵，则： max 2()scond A λλ=，max λ和s λ分别为模最大特征值与模最小特征值。
2、程序源代码：#include "Stdio.h"#include "Conio.h"#include "math.h"//****************************************************************************// // 在存储带状矩阵时，下面的几个量在程序中反复用到，为方便编程故把它们定义成宏.// // M ：转换后的矩阵的行数，M=R+S+1。
北航研究生数值分析作业第一题
北航研究⽣数值分析作业第⼀题北航研究⽣数值分析作业第⼀题：⼀、算法设计⽅案1.要求计算矩阵的最⼤最⼩特征值，通过幂法求得模最⼤的特征值，进⾏⼀定判断即得所求结果；2.求解与给定数值接近的特征值，可以该数做漂移量，新数组特征值倒数的绝对值满⾜反幂法的要求，故通过反幂法即可求得；3.反幂法计算时需要⽅程求解中间过渡向量，需设计Doolite分解求解；4.|A|=|B||C|，故要求解矩阵的秩，只需将Doolite分解后的U矩阵的对⾓线相乘即为矩阵的Det。
算法编译环境：vlsual c++6.0需要编译函数：幂法，反幂法，Doolite分解及⽅程的求解⼆、源程序如下：#include#include#include#includeint Max(int value1,int value2);int Min(int value1,int value2);void Transform(double A[5][501]);double mifa(double A[5][501]);void daizhuangdoolite(double A[5][501],double x[501],double b[501]); double fanmifa(double A[5][501]); double Det(double A[5][501]);/***定义2个判断⼤⼩的函数，便于以后调⽤***/int Max(int value1,int value2){return((value1>value2)?value1:value2);}int Min(int value1,int value2){return ((value1}/*****************************************//***将矩阵值转存在⼀个数组⾥，节省空间***/void Transform(double A[5][501],double b,double c){int i=0,j=0;A[i][j]=0,A[i][j+1]=0;for(j=2;j<=500;j++)A[i][j]=c;i++;j=0;A[i][j]=0;for(j=1;j<=500;j++)A[i][j]=b;i++;for(j=0;j<=500;j++)A[i][j]=(1.64-0.024*(j+1))*sin(0.2*(j+1))-0.64*exp(0.1/(j+1)); i++;for(j=0;j<=499;j++)A[i][j]=b;A[i][j]=0;i++;for(j=0;j<=498;j++)A[i][j]=c;A[i][j]=0,A[i][j+1]=0;}/***转存结束***///⽤于求解模最⼤的特征值，幂法double mifa(double A[5][501]){int s=2,r=2,m=0,i,j;double b2,b1=0,sum,u[501],y[501];for (i=0;i<=500;i++){u[i] = 1.0;}do{sum=0;if(m!=0)b1=b2;m++;for(i=0;i<=500;i++)sum+=u[i]*u[i];for(i=0;i<=500;i++)y[i]=u[i]/sqrt(sum);for(i=0;i<=500;i++){u[i]=0;for(j=Max(i-r,0);j<=Min(i+s,500);j++)u[i]=u[i]+A[i-j+s][j]*y[j];}b2=0;for(i=0;i<=500;i++)b2=b2+y[i]*u[i];}while(fabs(b2-b1)/fabs(b2)>=exp(-12));return b2;}//带状DOOLITE分解，并且求解出⽅程组的解void daizhuangdoolite(double A[5][501],double x[501],double b[501]) { int i,j,k,t,s=2,r=2;double B[5][501],c[501];for(i=0;i<=4;i++){for(j=0;j<=500;j++)B[i][j]=A[i][j];}for(i=0;i<=500;i++)c[i]=b[i];for(k=0;k<=500;k++){for(j=k;j<=Min(k+s,500);j++){for(t=Max(0,Max(k-r,j-s));t<=k-1;t++)B[k-j+s][j]=B[k-j+s][j]-B[k-t+s][t]*B[t-j+s][j]; }for(i=k+1;i<=Min(k+r,500);i++){for(t=Max(0,Max(i-r,k-s));t<=k-1;t++)B[i-k+s][k]=B[i-k+s][k]-B[i-t+s][t]*B[t-k+s][k]; B[i-k+s][k]=B[i-k+s][k]/B[s][k];}}for(i=1;i<=500;i++)for(t=Max(0,i-r);t<=i-1;t++)c[i]=c[i]-B[i-t+s][t]*c[t];x[500]=c[500]/B[s][500];for(i=499;i>=0;i--){x[i]=c[i];for(t=i+1;t<=Min(i+s,500);t++)x[i]=x[i]-B[i-t+s][t]*x[t];x[i]=x[i]/B[s][i];}}//⽤于求解模最⼤的特征值，反幂法double fanmifa(double A[5][501]){int s=2,r=2,m=0,i;double b2,b1=0,sum=0,u[501],y[501];for (i=0;i<=500;i++){u[i] = 1.0;}do{if(m!=0)b1=b2;m++;sum=0;for(i=0;i<=500;i++)sum+=u[i]*u[i];for(i=0;i<=500;i++)y[i]=u[i]/sqrt(sum);daizhuangdoolite(A,u,y);b2=0;for(i=0;i<=500;i++)b2+=y[i]*u[i];}while(fabs(b2-b1)>=fabs(b1)*exp(-12));return 1/b2;}//⾏列式的LU分解，U的主线乘积即位矩阵的DET double Det(double A[5][501]) {int i,j,k,t,s=2,r=2;for(k=0;k<=500;k++){for(j=k;j<=Min(k+s,500);j++){for(t=Max(0,Max(k-r,j-s));t<=k-1;t++)A[k-j+s][j]=A[k-j+s][j]-A[k-t+s][t]*A[t-j+s][j];}for(i=k+1;i<=Min(k+r,500);i++){for(t=Max(0,Max(i-r,k-s));t<=k-1;t++)A[i-k+s][k]=A[i-k+s][k]-A[i-t+s][t]*A[t-k+s][k];A[i-k+s][k]=A[i-k+s][k]/A[s][k];}}double det=1;for(i=0;i<=500;i++)det*=A[s][i];return det;}void main(){double b=0.16,c=-0.064,p,q;int i,j;double A[5][501];Transform(A,b,c); //进⾏A的赋值cout.precision(12); //定义输出精度double lamda1,lamda501,lamdas;double k=mifa(A);if(k>0) //判断求得最⼤以及最⼩的特征值.如果K>0，则它为最⼤特征值值，//并以它为偏移量再⽤⼀次幂法求得新矩阵最⼤特征值，即为最⼤ //与最⼩的特征值的差{lamda501=k;for(i=0;i<=500;i++)A[2][i]=A[2][i]-k;lamda1=mifa(A)+lamda501;for(i=0;i<=500;i++)A[2][i]=A[2][i]+k;}else //如果K<=0，则它为最⼩特征值值，并以它为偏移量再⽤⼀次幂法//求得新矩阵最⼤特征值，即为最⼤与最⼩的特征值的差{lamda1=k;for(i=0;i<=500;i++)A[2][i]=A[2][i]-k;lamda501=mifa(A)+lamda1;for(i=0;i<=500;i++)A[2][i]=A[2][i]+k;}lamdas=fanmifa(A);FILE *fp=fopen("result.txt","w");fprintf(fp,"λ1=%.12e\n",lamda1);fprintf(fp,"λ501=%.12e\n",lamda501);fprintf(fp,"λs=%.12e\n\n",lamdas);fprintf(fp,"\t要求接近的值\t\t\t实际求得的特征值\n");for(i=1;i<=39;i++) //反幂法求得与给定值接近的特征值{p=lamda1+(i+1)*(lamda501-lamda1)/40;for(j=0;j<=500;j++)A[2][j]=A[2][j]-p;q=fanmifa(A)+p;for(j=0;j<=500;j++)A[2][j]=A[2][j]+p;fprintf(fp,"µ%d: %.12e λi%d: %.12e\n",i,p,i,q);}double cond=fabs(mifa(A)/fanmifa(A));double det=Det(A);fprintf(fp,"\ncond(A)=%.12e\n",cond);fprintf(fp,"\ndetA=%.12e\n",det);}三、程序运⾏结果λ1=-1.069936345952e+001λ501=9.722283648681e+000λs=-5.557989086521e-003要求接近的值实际求得的特征值µ1: -9.678281104107e+000 λi1: -9.585702058251e+000µ2: -9.167739926402e+000 λi2: -9.172672423948e+000µ3: -8.657198748697e+000 λi3: -8.652284007885e+000µ4: -8.146657570993e+000 λi4: -8.0934********e+000µ5: -7.636116393288e+000 λi5: -7.659405420574e+000µ6: -7.125575215583e+000 λi6: -7.119684646576e+000µ7: -6.615034037878e+000 λi7: -6.611764337314e+000µ8: -6.104492860173e+000 λi8: -6.0661********e+000µ9: -5.593951682468e+000 λi9: -5.585101045269e+000µ10: -5.0834********e+000 λi10: -5.114083539196e+000µ11: -4.572869327058e+000 λi11: -4.578872177367e+000µ12: -4.062328149353e+000 λi12: -4.096473385708e+000µ13: -3.551786971648e+000 λi13: -3.554211216942e+000µ14: -3.0412********e+000 λi14: -3.0410********e+000µ15: -2.530704616238e+000 λi15: -2.533970334136e+000µ16: -2.020*********e+000 λi16: -2.003230401311e+000µ17: -1.509622260828e+000 λi17: -1.503557606947e+000µ18: -9.990810831232e-001 λi18: -9.935585987809e-001µ19: -4.885399054182e-001 λi19: -4.870426734583e-001µ20: 2.200127228676e-002 λi20: 2.231736249587e-002µ21: 5.325424499917e-001 λi21: 5.324174742068e-001µ22: 1.043083627697e+000 λi22: 1.052898964020e+000µ23: 1.553624805402e+000 λi23: 1.589445977158e+000µ24: 2.064165983107e+000 λi24: 2.060330427561e+000µ25: 2.574707160812e+000 λi25: 2.558075576223e+000µ26: 3.0852********e+000 λi26: 3.080240508465e+000µ27: 3.595789516221e+000 λi27: 3.613620874136e+000µ28: 4.106330693926e+000 λi28: 4.0913********e+000µ29: 4.616871871631e+000 λi29: 4.603035354280e+000µ30: 5.127413049336e+000 λi30: 5.132924284378e+000µ31: 5.637954227041e+000 λi31: 5.594906275501e+000µ32: 6.148495404746e+000 λi32: 6.080933498348e+000µ33: 6.659036582451e+000 λi33: 6.680354121496e+000µ34: 7.169577760156e+000 λi34: 7.293878467852e+000µ35: 7.680118937861e+000 λi35: 7.717111851857e+000µ36: 8.190660115566e+000 λi36: 8.225220016407e+000µ37: 8.701201293271e+000 λi37: 8.648665837870e+000µ38: 9.211742470976e+000 λi38: 9.254200347303e+000µ39: 9.722283648681e+000 λi39: 9.724634099672e+000cond(A)=1.925042185755e+003detA=2.772786141752e+118四、分析如果初始向量选择不当,将导致迭代中X1的系数等于零.但是，由于舍⼊误差的影响,经若⼲步迭代后,.按照基向量展开时,x1的系数可能不等于零。
北航数值分析实验报告
北航数值分析实验报告篇一：北航数值分析报告第一大题《数值分析》计算实习报告第一大题学号：DY1305姓名：指导老师：一、题目要求已知501*501阶的带状矩阵A，其特征值满足?1?2...?501。
试求：1、?1，?501和?s的值；2、A的与数?k??1?k?501??140最接近的特征值?ik（k=1，2，...，39）；3、A的（谱范数）条件数c nd(A)2和行列式de tA。
二、算法设计方案题目所给的矩阵阶数过大，必须经过去零压缩后进行存储和运算，本算法中压缩后的矩阵A1如下所示。
?0?0?A1??a1??b??c0b a2bcc bb c............c bb ccb a500b0a 3...a499c?b??a501??0?0??由矩阵A的特征值满足的条件可知?1与?501之间必有一个最大，则采用幂法求出的一个特征值必为其中的一个：当所求得的特征值为正数，则为?501；否则为?1。
在求得?1与?501其中的一个后，采用带位移的幂法则可求出它们中的另一个，且位移量即为先求出的特征值的值。
用反幂法求得的特征值必为?s。
由条件数的性质可得，c nd(A)2为模最大的特征值与模最小的特征值之比的模，因此，求出?1，?501和?s的值后，则可以求得c nd(A)2。
北航数值分析计算实习1
《数值分析》计算实习题目110091013 劳云杰一、算法设计方案根据提示的算法，首先使用幂法求出按模最大的特征值λt1，再根据已求出的λt1用带原点平移的幂法求出另一个特征值λt2，比较两个λ的大小，根据已知条件，可以得出λ1和λ501.至于λs，由于是按模最小的特征值，使用反幂法求之，由于反幂法需要解线性方程组，故对矩阵进行Doolittle分解。
再通过带原点平移的反幂法求跟矩阵的与数最接近的特征值。
对非奇异的矩阵A，根据条件数定义，取λt1/λs的绝对值，两个特征值在之前步骤中均以求得。
由于对矩阵进行了Doolittle分解，所以矩阵的行列式det A可由分解得出的上三角阵U 的对角线上元素相乘求得。
为了使A的所有零元素都不存储，使用书本25页的压缩存储法对A进行存储，在计算时通过函数在数组C中检索A中元素即可。
由于A是501*501矩阵，C应取为5*501矩阵。
由于数据不大，为了方便起见，在程序中取502*502矩阵或者502向量，C也取为6*502矩阵。
程序编写参考《数值分析》颜庆津著和[C数值算法].(美国)W ILLIAM.H．P RESS.扫描版。
二、全部源程序#include <stdio.h>#include <math.h>#define XS 1.0e-12//精度水平void fz_a();//对矩阵A赋值double js(int,int);//在压缩矩阵中检索A的元素double mf(double);//幂法double fmf(double);//反幂法int lu(double);//Doolittle分解int jfc(double[],double[]);//解方程int max(int,int);int min(int,int);double (*u)[502]=new double[502][502];//上三角阵double (*l)[502]=new double[502][502];//单位下三角阵double a[6][502];//压缩存储矩阵int max(int x,int y)//比大小函数×2{ return (x>y?x:y);}int min(int x,int y)//精度关系，比较下标用{ return (x<y?x:y);}int main(){printf("请耐心等待,先看看中间过程吧～\n");int i,k;double ldt1,ldt2,ld1,ld501,lds,mu[40],det;double ld[40];fz_a();//对A赋值ldt1=mf(0);//幂法求模最大的特征值ldt2=mf(ldt1);//以第一次求得的特征值进行平移ld1=ldt1<ldt2?ldt1:ldt2;//大的就是λ501ld501=ldt1<ldt2?ldt2:ldt1;lu(0);lds=fmf(0);//反幂法求λsdet=1;//初始化行列式for(i=1;i<=501;i++)det=det*u[i][i];//用U的对角元素求行列式for(k=1;k<=39;k++){mu[k]=ld1+k*(ld501-ld1)/40;//与数lu(mu[k]);ld[k]=fmf(mu[k]);}printf("\n 列出结果\n");printf("λ1=%1.12e λ501=%1.12e\n",ld1,ld501);printf("λs=%1.12e \n",lds);printf("cond(A)=%1.12e \n",fabs(ldt1/lds));printf("detA=%1.12e \n",det);for(k=1;k<=39;k++)//列出跟与数最接近特征值{printf("λi%d=%1.12e\t",k,ld[k]);if(k%2==0)printf("\n");}//界面友好性delete []u;delete []l;getchar();return 0;}void fz_a()//对A赋值{int i;for(i=3;i<=501;i++)a[1][i]=a[5][502-i]=-0.064;//原A矩阵的cfor(i=2;i<=501;i++)a[2][i]=a[4][502-i]=0.16;//原A矩阵的bfor(i=1;i<=501;i++)a[3][i]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i);//原对角线元素}double js(int i,int j)//对压缩矩阵检索A的元素{if(abs(i-j)<=2)return a[i-j+3][j];else return 0;}double mf(double offset)//幂法{int i,x1;double u[502],y[502];double beta=0,prebeta=-1000,yita=0;//用幂法的第一种迭代方法for(i=1;i<=501;i++) //用到了2-范数u[i]=1,y[i]=0;for(int k=1;k<=10000;k++)//对迭代次数进行限制{yita=0;for(i=1;i<=501;i++)yita=sqrt(yita*yita+u[i]*u[i]);for(i=1;i<=501;i++)y[i]=u[i]/yita;for(x1=1;x1<=501;x1++){u[x1]=0;for(int x2=1;x2<=501;x2++)u[x1]=u[x1]+((x1==x2)?(js(x1,x2)-offset):js(x1,x2))*y[x2];}prebeta=beta;beta=0;for(i=1;i<=501;i++)beta=beta+y[i]*u[i];if(fabs((prebeta-beta)/beta)<=XS){printf("offset=%f lb=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};}//满足误差条件后，迭代终止，并输出平移量，误差和迭代次数return(beta+offset);//加上平移量，方便比较}double fmf(double offset)//反幂法{ int i;double u[502],y[502];double beta=0,prebeta=0,yita=0;for(i=1;i<=501;i++)u[i]=1,y[i]=0; //相关量初始化for(int k=1;k<=10000;k++)//限制迭代次数{yita=0;for(i=1;i<=501;i++)yita=sqrt(yita*yita+u[i]*u[i]);for(i=1;i<=501;i++)y[i]=u[i]/yita;jfc(u,y);prebeta=beta;beta=0;for(i=1;i<=501;i++)beta=beta+y[i]*u[i];beta=1/beta;if(fabs((prebeta-beta)/beta)<=XS){printf("offset=%f lb=%f err=%ek=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};}//满足误差条件后，迭代终止，并输出平移量，误差和迭代次数return(beta+offset);}int lu(double offset)//Doolittle分解{int i,j,k,t;double sum;//中间量for(k=1;k<=501;k++)for(j=1;j<=501;j++){u[k][j]=l[k][j]=0;if(k==j)l[k][j]=1;}//对LU矩阵初始化for(k=1;k<=501;k++)//对式（2.12）的程序实现{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+l[k][t]*u[t][j];//j=k,k+1,……,nu[k][j]=((k==j)?(js(k,j)-offset):js(k,j))-sum;}if(k==501)continue;for(i=k+1;i<=min(k+2,501);i++)//i=k+1,……，n{sum=0;for(t=max(1,max(i-2,k-2));t<=(k-1);t++)sum=sum+l[i][t]*u[t][k];l[i][k]=(((i==k)?(js(i,k)-offset):js(i,k))-sum)/u[k][k];}}return 0;}int jfc(double x[],double b[])//解方程{int i,t;double y[502];double sum;y[1]=b[1];for(i=2;i<=501;i++){sum=0;for(t=max(1,i-2);t<=i-1;t++)sum=sum+l[i][t]*y[t];y[i]=b[i]-sum;}x[501]=y[501]/u[501][501];for(i=500;i>=1;i--){sum=0;for(t=i+1;t<=min(i+2,501);t++)sum=sum+u[i][t]*x[t];x[i]=(y[i]-sum)/u[i][i];}return 0;}三、结果λ1=-1.070011361502e+001λ501=9.724634098777e+000λs=-5.557910794230e-003cond(A)=1.925204273902e+003detA=2.772786141752e+118λi1=-1.018293403315e+001 λi2=-9.585707425068e+000 λi3=-9.172672423928e+000λi4=-8.652284007898e+000 λi5=-8.0934********e+000 λi6=-7.659405407692e+000λi7=-7.119684648691e+000 λi8=-6.611764339397e+000 λi9=-6.0661********e+000λi10=-5.585101052628e+000 λi11=-5.114083529812e+000 λi12=-4.578872176865e+000λi13=-4.096470926260e+000 λi14=-3.554211215751e+000 λi15=-3.0410********e+000 λi16=-2.533970311130e+000 λi17=-2.003230769563e+000 λi18=-1.503557611227e+000 λi19=-9.935586060075e -001 λi20=-4.870426738850e -001 λi21=2.231736249575e -002 λi22=5.324174742069e -001 λi23=1.052898962693e+000 λi24=1.589445881881e+000 λi25=2.060330460274e+000 λi26=2.558075597073e+000 λi27=3.080240509307e+000 λi28=3.613620867692e+000 λi29=4.0913********e+000 λi30=4.603035378279e+000 λi31=5.132924283898e+000 λi32=5.594906348083e+000 λi33=6.080933857027e+000 λi34=6.680354092112e+000 λi35=7.293877448127e+000 λi36=7.717111714236e+000 λi37=8.225220014050e+000 λi38=8.648666065193e+000 λi39=9.254200344575e+000四、讨论迭代初始向量的选取对计算结果的影响1.在反幂法中取迭代向量u[1]=1,u[i]=0，i=2,……，501，最后得出的结果中λs=2.668886923785e -002，cond(A)也随之改变成4.009204556274e+0022.在幂法中取迭代向量u[1]=1,u[i]=2,i=2,……,501,最后得出的结果不变。
北航数值分析大作业一
《数值分析B》大作业一SY1103120 朱舜杰一．算法设计方案：1.矩阵A的存储与检索将带状线性矩阵A[501][501]转存为一个矩阵MatrixC[5][501] .由于C语言中数组角标都是从0开始的，所以在数组MatrixC[5][501]中检索A的带内元素a ij的方法是：A的带内元素a ij=C中的元素c i-j+2,j2.求解λ1，λ501，λs①首先分别使用幂法和反幂法迭代求出矩阵按摸最大和最小的特征值λmax和λmin。
λmin即为λs；如果λmax>0,则λ501=λmax；如果λmax<0,则λ1=λmax。
②使用带原点平移的幂法（mifa（）函数），令平移量p=λmax，求出对应的按摸最大的特征值λ，max，如果λmax>0,则λ1=λ，max+p；如果λmax<0,则λ501=λ，max+p。
3.求解A的与数μk=λ1+k（λ501-λ1）/40的最接近的特征值λik （k=1,2，…，39）。
使用带原点平移的反幂法，令平移量p=μk，即可求出与μk最接近的特征值λik。
4.求解A的（谱范数）条件数cond（A）2和行列式d etA。
①cond（A）2=|λ1/λn|，其中λ1和λn分别是矩阵A的模最大和最小特征值。
②矩阵A的行列式可先对矩阵A进行LU分解后，detA等于U所有对角线上元素的乘积。
二．源程序#include<stdio.h>#include<iostream.h>#include<stdlib.h>#include<math.h>#include<float.h>#include<iomanip.h>#include<time.h>#define E 1.0e-12 /*定义全局变量相对误差限*/int max2(int a,int b) /*求两个整型数最大值的子程序*/{if(a>b)return a;elsereturn b;}int min2(int a,int b) /*求两个整型数最小值的子程序*/{if(a>b)return b;elsereturn a;}int max3(int a,int b,int c) /*求三整型数最大值的子程序*/{ int t;if(a>b)t=a;else t=b;if(t<c) t=c;return(t);}void assignment(double array[5][501]) /*将矩阵A转存为数组C[5][501]*/{int i,j,k;//所有元素归零for(i=0;i<=4;){for(j=0;j<=500;){array[i][j]=0;j++;}i++;}//第0,4行赋值for(j=2;j<=500;){k=500-j;array[0][j]=-0.064;array[4][k]=-0.064;j++;}//第1,3行赋值for(j=1;j<=500;){k=500-j;array[1][j]=0.16;array[3][k]=0.16;j++;}//第2行赋值for(j=0;j<=500;){ k=j;j++;array[2][k]=(1.64-0.024*j)*sin((double)(0.2*j))-0.64*exp((double)(0.1/j));}}double mifa(double u[501],double array[5][501],double p) /*带原点平移的幂法*/ {int i,j; /* u[501]为初始迭代向量*/double a,b,c=0; /* array[5][501]为矩阵A的转存矩阵*/double y[501]; /*p为平移量*/for(;;){a=0;b=0;/*选用第一种迭代格式*///求ηk-1for(i=0;i<=500;i++){a=a+u[i]*u[i];}a=sqrt(a);//求y k-1for(i=0;i<=500;i++){y[i]=u[i]/a;}//求u kfor(i=0;i<=500;i++){u[i]=0;for(j=max2(i-2,0);j<=min2(i+2,500);j++){u[i]+=array[i-j+2][j]*y[j];}u[i]=u[i]-p*y[i]; /*引入平移量*/}//求βkfor(i=0;i<=500;i++){b+=y[i]*u[i];}if(fabs((b-c)/b)<=E) /*达到精度水平，迭代终止*/break;c=b;}return (b+p); /*直接返回A的特征值*/}void chuzhi(double a[]) /*用随机数为初始迭代向量赋值*/ {int i;srand((int)time(0));for(i=0;i<=500;i++){a[i]=(10.0*rand()/RAND_MAX); /*生成0~10的随机数*/}}void chuzhi2(double a[],int j) /*令初始迭代向量为e i*/{int i;for(i=0;i<=500;i++){a[i]=0;}a[j]=1;}void LU(double array[5][501]) /*对矩阵A进行Doolittle分解*/{ /*矩阵A转存在C[5][501]中*/int j,k,t; /*分解结果L，U分别存在C[5][501]的上半部与下半部*/ for(k=0;k<=500;k++){for(j=k;j<=min2((k+2),500);j++){for(t=max3(0,k-2,j-2);t<=(k-1);t++){array[k-j+2][j]-=array[k-t+2][t]*array[t-j+2][j];}}if(k<500)for(j=k+1;j<=min2((k+2),500);j++){for(t=max3(0,k-2,j-2);t<=(k-1);t++){array[j-k+2][k]-=array[j-t+2][t]*array[t-k+2][k];}array[j-k+2][k]=array[j-k+2][k]/array[2][k];}}}double fmifa(double u[501],double array[5][501],double p){ /*带原点平移的反幂法*/ int i,j;double a,b,c=0;double y[501];//引入平移量for(i=0;i<=500;i++){array[2][i]-=p;}//先将矩阵Doolittle分解LU(array);for(;;){a=0;b=0;//求ηk-1for(i=0;i<=500;i++){a=a+u[i]*u[i];}a=sqrt(a);//求y k-1for(i=0;i<=500;i++){y[i]=u[i]/a;}//回带过程，求解u kfor(i=0;i<=500;i++){u[i]=y[i];}for(i=1;i<=500;i++){for(j=max2(0,(i-2));j<=(i-1);j++){u[i]-=array[i-j+2][j]*u[j];}}u[500]=u[500]/array[2][500];for(i=499;i>=0;i--){for(j=i+1;j<=min2((i+2),500);j++){u[i]-=array[i-j+2][j]*u[j];}u[i]=u[i]/array[2][i];}//求βkfor(i=0;i<=500;i++){b+=y[i]*u[i];}if(fabs((b-c)/b)<=E) /*达到精度要求，迭代终止*/break;c=b;}return (p+(1/b)); /*直接返回距离原点P最接近的A的特征值*/ }//主函数main(){ int i;double d1,d501,ds,d,a;double u[501];double MatrixC[5][501];printf(" 《数值分析》计算实习题目第一题\n");printf(" SY1103120 朱舜杰\n");//将矩阵A转存为MatrixCassignment(MatrixC);//用带原点平移的幂法求解λ1，λ501chuzhi(u);d=mifa(u,MatrixC,0);chuzhi(u);a=mifa(u,MatrixC,d);if(d<0){d1=d;d501=a;}else{d501=d;d1=a;}printf("λ1=%.12e\n",d1);printf("λ501=%.12e\n",d501);//用反幂法求λschuzhi(u);ds=fmifa(u,MatrixC,0);printf("λs=%.12e\n",ds);//用带原点平移的反幂法求λikfor(i=1;i<=39;i++){a=d1+(i*(d501-d1))/40;assignment(MatrixC);chuzhi(u);d=fmifa(u,MatrixC,a);printf("与μ%02d=%+.12e最接近的特征值λi%02d=%+.12e\n",i,a,i,d);}//求A的条件数d=fabs((d1/ds));printf("A的（谱范数）条件数cond<A>2=%.12e\n",d);//求detAassignment(MatrixC);LU(MatrixC);a=1;for(i=0;i<=500;i++){a*=MatrixC[2][i];}printf("行列式detA=%.12e\n",a);//测试不同迭代初始向量对λ1计算结果的影响。
北航数值分析大作业一
北京航空航天大学数值分析大作业一学院名称自动化专业方向控制工程学号ZY*******学生姓名许阳教师孙玉泉日期2021 年11月26 日设有501501⨯的实对称矩阵A ，⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=5011A a b c b c c b c b a其中，064.0,16.0),501,,2,1(64.0)2.0sin()024.064.1(1.0-==⋅⋅⋅=--=c b i e i i a ii 。
矩阵A 的特征值为)501,,2,1(⋅⋅⋅=i i λ，并且有||min ||,501150121i i s λλλλλ≤≤=≤⋅⋅⋅≤≤1λ，501λ和s λ的值。
A 的与数4015011λλλμ-+=kk 最接近的特征值)39,,2,1(⋅⋅⋅=k k i λ。
A 的(谱范数)条件数2)A (cond 和行列式detA 。
一 方案设计1 求1λ，501λ和s λ的值。
s λ为按模最小特征值，||min ||5011i i s λλ≤≤=。
可使用反幂法求得。
1λ，501λ分别为最大特征值及最小特征值。
可使用幂法求出按模最大特征值，如结果为正，即为501λ，结果为负，那么为1λ。
使用位移的方式求得另一特征值即可。
2 求A 的与数4015011λλλμ-+=kk 最接近的特征值)39,...,2,1(=k k i λ。
题目可看成求以k μ为偏移量后，按模最小的特征值。
即以k μ为偏移量做位移，使用反幂法求出按模最小特征值后，加上k μ，即为所求。
3 求A 的(谱范数)条件数2)(A cond 和行列式detA 。
矩阵A 为非奇异对称矩阵，可知，||)(min max2λλ=A cond(1-1)其中m ax λ为按模最大特征值，min λ为按模最小特征值。
detA 可由LU 分解得到。
因LU 均为三角阵，那么其主对角线乘积即为A 的行列式。
二 算法实现1 幂法使用如下迭代格式：⎪⎪⎩⎪⎪⎨⎧⋅===⋅⋅⋅=------||max |)|sgn(max ||max /),,(111111)0()0(10k k k k k k k k Tn u u Ay u u u y u u u β任取非零向量 (2-1)终止迭代的控制理论使用εβββ≤--||/||1k k k ， 实际使用εβββ≤--||/||||||1k k k(2-2)由于不保存A 矩阵中的零元素，只保存主对角元素a[501]及b,c 值。
北航数值分析计算实习第一题编程
i − t + s +1,t t − k + s +1, k t = max(1,i − r ,k − s )
∑c
c
) / cs +1, k
[i = k + 1, k + 2,⋯ , min( k + r , n); k < n]
（2） 求解 Ly = b,Ux = y (数组 b 先是存放原方程右端向量，后来存放中间向量 y)
0 b a2
b c
c b a3 b c
⋯ ⋯ ⋯ ⋯ ⋯
c b a499 b c
c b a500 b 0
c ⎤ b ⎥ ⎥ a501 ⎥ ⎥ 0 ⎥ 0 ⎥ ⎦
在数组 C 中检索矩阵 A 的带内元素 aij 的方法是： A 的带内元素 aij =C 中的元素 ci − j + s +1, j
2
数值分析计算实习题目一
i −1
bi := bi −
பைடு நூலகம்
i − t + s +1,t t t = max(1,i − r )
∑c
b
(i = 2,3,⋯ , n)
xn := bn / cs +1, n
min( i + s )
xi := (bi −
t = i +1
∑c
i −t + s +1,t t
x ) / cs +1,i
(i = n − 1, n − 2,⋯ ,1)
3、Doolittle 分解求解 n 元带状线性方程组（doolittle（）函数）
按照上述对带状矩阵 A 的存储方法和元素 aij 的检索方法，并且把三角分解的结果 ukj 和 lik 分 别存放在 akj 和 aik 原先的存储单元内，那么用 Doolittle 分解法求解 n 元带状线性方程组的算法 可重新表述如下（其中“:=”表示赋值） ： （1） 作分解 A = LU 。 对于 k=1，2， ……，n 执行
北航数值分析-实习作业1(C语言详细注释)
《数值分析》计算实习作业《一》北航第一题 设有501501⨯的矩阵⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=501500499321a bc b a b cc b a b ccb a bc c b a b c b a A其中.064.0,16.0);501,2,1(64.0)2.0sin()024.064.1(1.0-===--=c b i e i i a i i 矩阵的特征值)501,,2,1( =i i λ满足||min ||,501150121i i s λλλλλ≤≤=<<<试求1. 5011，λλ和s λ的值2. 的与数4015011λλκλμ-+=k 最接近的特征值)39,,2,1( =K κλi3. 的（谱范数）条件数2)A (cond 和行列式A det 要求1. 算法的设计方案（A 的所有零元素都不能存储）2. 全部源程序（详细注释）。
变量为double ，精度-1210=ε，输出为e 型12位有效数字3. 特征值s 5011，，λλλ和)39,,2,1( =K κλi 以及A cond det ,)A (2的值4. 讨论迭代初始向量的选取对计算结果的影响，并说明原因解答：1. 算法设计对于s λ满足||min ||5011i i s λλ≤≤=，所以s λ是按模最小的特征值，直接运用反幂法可求得。
对于5011，λλ，一个是最大的特征值，一个是最小的特征值，不能确定两者的绝对值是否相等，因此必须首先假设||||5011λλ≠，然后运用幂法，看能否求得一个特征值，如果可以求得一个，证明A 是收敛的，求得的结果是正确的，然后对A 进行带原点平移的幂法，偏移量是前面求得的特征值，可以求得另一个特征值，最后比较这两个特征值，较大的特征值是501λ，较小的特征值就是1λ。
如果在假设的前提下，无法运用幂法求得按模最大的特征值，即此时A 不收敛，则需要将A 进行带原点平移的幂法，平移量可以选取1，再重复上述步骤即可求得两个特征值。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值分析大作业一、算法设计方案1、矩阵初始化矩阵[]501501⨯=ij a A 的下半带宽r=2,上半带宽s=2，设置矩阵[][]5011++s r C ，在矩阵C 中检索矩阵A 中的带内元素ij a 的方法是：j s j i ij c a ,1++-=。
这样所需要的存储单元数大大减少，从而极大提高了运算效率。
2、利用幂法求出5011λλ，幂法迭代格式：0111111nk k k k kk T k k k u R y u u Ay y u ηηβ------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩非零向量 当1210/-≤-k k βββ时，迭代终止。
首先对于矩阵A 利用幂法迭代求出一个λ，然后求出矩阵B ，其中I A B λ-=(I 为单位矩阵），对矩阵B 进行幂法迭代，求出λ'，之后令λλλ+'=''，比较的大小与λλ''，大者为501λ，小者为1λ。
3、利用反幂法求出ik s λλ,反幂法迭代格式：0111111nk k k k kk T k k k u R y u Au y y u ηηβ------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩非零向量 当1210/-≤-k k βββ时，迭代终止，1s k λβ=。
每迭代一次都要求解一次线性方程组1-=k k y Au ，求解过程为：（1）作分解LU A =对于n k ,...,2,1=执行[][]s k n r k k k i c c c c c n s k k k j c cc c k s ks k t k s k r i t t s t i k s k i k s k i js j t k s j r k t t s t k j s j k j s j k <+++=-=++=-=+++----=++-++-++-++----=++-++-++-∑∑);,min(,...,2,1/)(:),min(,...,1,:,1,11),,1max(,1,1,1,11),,1max(,1,1,1（2）求解y Ux b Ly ==,（数组b 先是存放原方程组右端向量，后来存放中间向量y))1,...,2,1(/)(:/:),...,3,2(:,1),min(1.1.11),1max(,1--=-===-=+++-++-+--=++-∑∑n n i c x c b x c b x n i b c b b i s t n s i i t t s t i i i ns n n ti r i t t s t i i i使用反幂法，直接可以求得矩阵按模最小的特征值s λ。
求与数)39,...,2,1(4015011=-+=k k k λλλμ最接近的特征值ik λ，对矩阵IA k μ-实行反幂法，即可求出对应的k k ik k k μλλβλ+==,/1。
4、求出A 的条件数和行列式根据max 2()scond A λλ=，其中分子分母分别对应按模最大和最小的特征值。
det()A 的计算：由于A LU =,其中L 为下三角矩阵，且对角线元素为1，故det()1L =，所以有A LU U ==，又U 为上三角矩阵，故det()U 为对其对角线上各元素的乘积，最后可得det()det()A U =。
二、程序源代码（1）定义所需要的函数：#include <stdio.h>#include <conio.h>#include <math.h>#define N 501#define R 2#define S 2int min(int a,int b); // 求最小值int max(int a,int b,int c); // 求最大值double Fan_two(double x[N]);//计算二范数void FenjieLU(double (*C)[N]);//解线性方程组的LU分解过程void Solve(double (*C)[N], double *b,double *x);//解线性方程组的求解过程double PowerMethod(double C[][N],double u[N],double y[N],double bta,double D);//幂法double InversePowerMethod(double C[][N],double u[N],double y[N],double bta,double D);//反幂法};（2）程序的主函数，Main.cpp代码如下：void main(){double C[R+S+1][N];double u[N];double y[N];double miu[39];double C1[R+S+1][N];double bta = 1.0;double Namda1,Namda501,NamdaS;double Namda[39];double CondA2;double detA = 1.0;double D = 1.0e-12;int i, j, k;FILE * fp;fp = fopen("Namda.txt","w");//对数组进行初始化//int i, j;for (i = 0; i < N; i++){u[i] = 1;}for (i = 0;i< R + S + 1;i++){for (j = 0;j< N;j++){if (i==0||i==4){C[i][j]=-0.064;}else if (i==1||i==3){C[i][j]=0.16;}else if (i==2){C[i][j]=(1.64-0.024*(j+1))*sin(0.2*(j+1))-0.64*exp(0.1/(j+1));}}}//幂法求Namda1//Namda1 = PowerMethod(C, u, y, bta, D);printf("\n================================================\n"); printf("Namda1 = %12.11e", Namda1);printf("\n================================================\n"); //幂法求Namda501//bta = 1.0;for (i = 0; i < R + S + 1; i++){for (j = 0; j < N; j++){if (i == 2)C1[i][j] = C[i][j] -Namda1;elseC1[i][j] = C[i][j];}}Namda501 = algorism.PowerMethod(C1, u, y, bta, D) +Namda1;printf("\n================================================\n"); printf("Namda501 = %12.11e", Namda501);printf("\n================================================\n"); //反幂法求NamdaS//bta = 1.0;NamdaS = InversePowerMethod(C, u, y, bta, D);printf("\n================================================\n");printf("NamdaS = %12.11e", NamdaS);printf("\n================================================\n");//反幂法求Namda[k]//printf("\n================================================\n");for (k = 0; k < 39; k++){miu[k] = Namda1 + (k + 1) * (Namda501 - Namda1) / 40.0;bta = 1.0;for (i = 0; i < R + S + 1; i++){for (j = 0; j < N; j++){if (i == 2)C1[i][j] = C[i][j] - miu[k];elseC1[i][j] = C[i][j];}}Namda[k] = InversePowerMethod(C1, u, y, bta, D) + miu[k];fprintf(fp,"与%12.11e最接近的特征值为:%12.11e\n",miu[k],Namda[k]);}printf("求与miu[k]最接近的Namda[k]的计算结果已经输出到文件Namda.txt 中");printf("\n================================================\n");//求A的谱范数//printf("\n================================================\n");printf("A的谱范数为:%12.11e", sqrt(Namda501));printf("\n================================================\n");//求A的条件数//CondA2 = fabs( Namda1 / NamdaS);printf("\n================================================\n");printf("A的谱范数的条件数Cond(A)2为:%12.11e",CondA2);printf("\n================================================\n");//求det(A)2的值//for (j = 0; j < N; j++)detA *= C[2][j];printf("\n================================================\n");printf("行列式A的值为:%12.11e",detA);printf("\n================================================\n");fclose(fp);_getch();return;}（3）成员函数的实现int min(int a,int b){return a < b ? a : b;}int max(int a,int b,int c){int temp;temp = a > b ? a : b;return temp > c ? temp : c;}double Fan_two(double x[N]){double sum = 0.0;int i;for (i = 0; i < N; i++){sum += pow(x[i],2);}return sqrt(sum);}void FenjieLU(double (*C)[N]){double sum = 0;int i, j, k,t;for (k = 0; k < N; k++){j = k;i = k + 1;while (1){if (j == min(k + S + 1, N))break;for (t = max(0, k - R, j - S); t <= k - 1; t++){sum += C[k-t+S][t] * C[t-j+S][j];}C[k-j+S][j] = C[k-j+S][j] - sum;sum = 0.0;j++;if (k == N-1)break;if (i == min(k + R + 1, N))break;for (t = max(0, i - R,k - S); t <= k - 1; t++){sum += C[i-t+S][t] * C[t-k+S][k];}C[i-k+S][k] = (C[i-k+S][k] - sum) / C[S][k];sum = 0;i++;}}}void Solve(double (*C)[N], double *b,double *x){double sum = 0;int i, t;sum = 0;for (i = 1; i < N; i++){for (t = max(0, i - R); t <= i - 1; t++){sum += C[i-t+S][t] * b[t];}b[i] = b[i] - sum;sum = 0;}x[N-1] = b[N-1] / C[S][N-1];for (i = N - 2; i >= 0; i--){for (t = i+1; t <= min(i + S, N - 1); t++){sum += C[i-t+S][t] * x[t];}x[i] = (b[i] - sum) / C[S][i];sum = 0;}}double PowerMethod(double C[][N],double u[N],double y[N],double bta,double D) {double ita;double sum = 0;double temp = 0.0;int i,j,k = 0;while (fabs(bta - temp) / fabs(bta) > D){temp = bta;ita = Fan_two(u);for (i = 0; i < N; i++){y[i] = u[i] / ita;}for (i = 0; i < N; i++){for (j = max(0,i - R); j < min(i + S + 1,N); j++){sum += C[i - j + S][j] * y[j];}u[i] = sum;sum = 0;}for (i = 0; i < N; i++){sum += y[i] * u[i];}bta = sum;sum = 0;k++;}return bta;}double InversePowerMethod(double C[][N],double u[N],double y[N],double bta,double D){double TC[R+S+1][N];double ty[N];double ita;double sum = 0;double temp = 0.0;int i,j,k = 0;FenjieLU(C);while (abs(1/bta - 1/temp) / abs(1/bta) > D){temp = bta;ita = Fan_two(u);for (i = 0; i < N; i++){y[i] = u[i] / ita;}//用到临时存储数组TC[][]和ty[][]是因为函数Solve执行过程中会改变A[][]和y[][]for (i = 0; i < R + S + 1; i++){for (j = 0; j < N; j++)TC[i][j] = C[i][j];}for (i = 0; i < N; i++)ty[i] = y[i];Solve(C, y, u);for (i = 0; i < R+S+1; i++){for (j = 0; j < N; j++)C[i][j] = TC[i][j];}for (i = 0; i < N; i++)y[i] = ty[i];for (i = 0; i < N; i++){sum += y[i] * u[i];}bta = sum;sum = 0;k++;}bta = 1.0 / bta;return bta;}三、程序运行结果下图为主程序运行结果的结果输出在Namda.txt文件中，结果如下：其中ik四、分析迭代初始向量对计算结果的影响选择不同的初始向量[]u N可能会得到不同的特征值。