数值分析大作业二.doc

《数值分析大作业二.doc》由会员分享，可在线阅读，更多相关《数值分析大作业二.doc（17页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

1、数值分析大作业二一、算法的设计总体思路：采用带双步位移的QR分解法求解矩阵A的特征值，给定精度水平为E，以及最大的迭代次数L=1000；采用高斯分解法求解矩阵A特征值对应的特征向。1、对矩阵A进行拟上三角化。记=A，并记的第r至第n列的元素为若(i=r+1,r+2,n;j=r,r+1,n)。对于r=1,2,.,n-2执行(1) 若(i=r+2,r+3,n)全为0，则令，转(5),否则转(2)。(2)计算=-sgn()(若=0，取=)(3)令。(4)计算(5)继续。当此算法执行完后，就得到与原矩阵A相似的拟上三角矩阵。2、对矩阵A进行QR分解，即求得矩阵Q、R及RQ。记=A，并记=，令=I.对于

2、r=1,2，n-1，执行(1)若(i=r+1,r+2,n)全为0，则令, 。转（5），否则转（2）。(2)计算=-sgn()(若=0，取=)(3)令。(4)计算(5)继续。当此算法执行完后，就得到正交矩阵和上三角矩阵R=且有A=QR。3、对矩阵A进行带双步位移的QR分解法，求得其全部特征值，这也就是矩阵A的特征值。具体算法如下：(1)使用矩阵的拟上三角化算法把矩阵A化为拟上三角矩阵；给定精度水平E1，则转(3)。(5)求二阶子阵的两个特征值和，即计算二次方程的两个根和。(6)如果m=2，则得到A的两个特征值和，转(11)，否则转(7)。(7)如果，则得到A的两个特征值和，置m：=m-2，转(4

3、)，否则则转(8)。(8)如果k=L，则计算终止，未得到A的全部特征值，否则转(9)。(9)记（），计算（I是m阶单位矩阵）（对做QR分解）(10)置k=k+1，转(3)。(11)A的全部特征值已经计算完毕，停止计算。其中的QR分解和的计算用下列算法实现:记。对r=1,2，m-1执行(1)若全为零，则令，。(2)计算，(若=0，则取=)(3)令。(4)计算(5)继续。此算法执行完后，就得到。4、最后再通过列主元高斯消去法，得到矩阵A的特征向量。先对原矩阵A作原点平移，平移量为它的实数特征值，再对平移之后的矩阵A作QR分解。解方程组AX=0求特征向量，即RX=0。因为R为上三角阵，令特征向量的最

4、后一个元素为1，然后再进行迭代求解。这样，就求出了所有的值。二、全部源程序#include#include#define n 10#define E 1.0e-12 /*精确度为*/#define L 1000 /*定义迭代次数为1000次*/int twostep(double A111+1,double Vector3,FILE *fp);void upper_triangular_matrix(double A1111);void QR (double C11,double B11,int m);void solve_equation(double s1,double s2,double

5、 Vector3,int m);void Gauss(double a111,double x111,double Vector3,int p);void EigenVector(double A11,double x11,double Vector3,FILE *fp);/*矩阵A的拟上三角化*/void upper_triangular_matrix(double A1111)int r,k,j,m;double dr,cr,tr,hr;double judge;double sum, sum1,sum2,sum3;double ur11,pr11,qr11,wr11;for(r=1;r=

6、n-2;r+)judge=0; for (k=r+2;k=n;k+)judge=judge+Akr*Akr;if(judge=0)continue;elsesum=0;for(j=r+1;j=n;j+)sum+=Ajr*Ajr;dr=sqrt(double (sum);if(Ar+1r=0)cr=dr;elsecr=-dr;hr=cr*cr-cr*Ar+1r;for(m=1;m=n;m+)if(m=r)urm=0;else if(m=(r+1)urm=Amr-cr;elseurm=Amr;for(m=1;m=n;m+) sum1=0; sum2=0;for(k=1;k=n;k+) sum1+=

7、(double)(Akm*urk); prm=sum1/hr;for(k=1;k=n;k+)sum2+=(double)(Amk*urk); qrm=sum2/hr; sum3=0;for(m=1;m=n;m+) sum3+=(double)(prm*urm);tr=(double)(sum3/hr);for(m=1;m=n;m+)wrm=qrm-tr*urm; for(m=1;m=n;m+) for(k=1;k=n;k+) Amk=Amk-wrm*urk-urm*prk; /*矩阵QR分解*/void QR(double C11,double B11,int m) int i,j,r; in

8、t flag=0; double cr,dr,hr,tr,temp; double pr11,qr11,ur11,vr11,wr11; for(r=1;r=m-1;r+) cr=0.0;dr=0.0;hr=0.0;tr=0.0;temp=0.0; for(i=1;i=m;i+) pri=0.0;qri=0.0;uri=0.0;vri=0.0;wri=0.0; for(i=r+1;i=E) flag=1; break; if(flag=1) for(i=r;iE) cr = -dr; else cr = dr; hr = cr*(cr-Brr); for(i=r;i=m;i+) uri = Bi

9、r; urr = urr-cr; for(i=1;i=m;i+) for(j=1;j=m;j+) vri = vri+Bji*urj; pri = pri+Cji*urj; qri = qri+Cij*urj; vri = vri/hr; pri = pri/hr; qri = qri/hr; for(i=1;i=m;i+) tr = tr+pri*uri; tr = tr/hr; for(i=1;i=m;i+) wri = qri-tr*uri; for(i=1;i=m;i+) for(j=1;j=0.0) s10 = (t1+sqdelt)/2.0; s20 = (t1-sqdelt)/2

10、.0; /*不存在实根,设Vector2为复数标志量*/ else Vectorm2=1;Vectorm-12=1; s10=t1/2.0;s11=sqdelt/2.0; s20=t1/2.0;s21=-sqdelt/2.0; /*带双步位移的QR分解*/int twostep(double A11,double Vector3,FILE *fp) int i,j,k,m,l; int flag=1; double s12=0.0,s22=0.0,s,t; double Mkn+1n+1=0.0,In+1n+1=0.0; k=0;m=n; for(i=1;i=n;i+) Iii=1.0; wh

11、ile(1) STEP1: if(fabs(Amm-1)=E) Vectorm0=Amm; m=m-1;STEP2: if(m=1) Vectorm0=Amm; return(flag); else if(m=0) return(flag); else goto STEP1; s10 = Am-1m-1+Amm;s11=0; s20 = Am-1m-1*Amm-Am-1m*Amm-1;s21=0; solve_equation(s1,s2,Vector,m); if(m=2) for(i=0;i2;i+) Vectormi = s1i; Vectorm-1i = s2i; return(fla

12、g); if(fabs(Am-1m-2)=E) for(i=0;i2;i+) Vectormi = s1i; Vectorm-1i = s2i; m=m-2; goto STEP2; if(k=L) printf(计算失败！不能得到所有特征值！n); flag=0; return(flag); s = Am-1m-1+Amm; t = Am-1m-1*Amm-Am-1m*Amm-1; for(i=1;i=m;i+) for(j=1;j=m;j+) Mkij=0; for(i=1;i=m;i+) for(j=1;j=m;j+) for(l=1;l=m;l+) Mkij = Mkij+Ail*Al

13、j; for(i=1;i=m;i+) for(j=1;j=m;j+) Mkij = Mkij-s*Aij+t*Iij; QR(A,Mk,m); k=k+1; /*列主元Gauss消去法*/void Gauss(double a111,double x11,double Vector3,int p) int i,j,k,maxindex; double max,mik,temp=0; double A1n+1n+1=0,bn+1=0; for(i=1;i=n;i+) for(j=1;j=n;j+) A1ij = a1ij;for(i=1;i=n;i+)A1ii = A1ii-Vectorp0;/

14、*消元过程*/for(k=1;k=n-1;k+)max=A1kk;maxindex=k; for(i=k+1;ifabs(max) maxindex=i; max=A1ik; if(maxindexk)for(j=k;j=n;j+) temp = A1kj;A1kj = A1maxindexj;A1maxindexj = temp;for(i=k+1;i=n;i+)mik = A1ik/A1kk;for(j=k+1;j=1;k-) temp=0.0; for(j=k+1;j=n;j+) temp = temp+A1kj*xpj; xpk = -temp/A1kk; /*特征向量归一化*/ te

15、mp=0.0; for(i=1;i=n;i+) temp += pow(xpi,2); temp = sqrt(temp); for(i=1;i=n;i+) xpi=xpi/temp;/*求实根特征值对应特征向量*/void EigenVector(double An+1,double xn+1,double Vector3,FILE *fp) int i,j; for(i=1;i=n;i+) /*求实特征值的特征向量*/ if(Vectori2=0) Gauss(A,x,Vector,i); printf(n当特征值为%.12e时,对应的特征向量为n,Vectori0); fprintf(f

16、p,n当特征值为%.12e时,对应的特征向量为n,Vectori0); for(j=1;j=n;j+) printf(%.12en,xij); fprintf(fp,%.12en,xij); printf(n); fprintf(fp,n); /*主函数程序*/void main() int i=0,j=0,flag=0; double an+1n+1=0.0,P1111=0.0,Vector113=0.0,x1111=0.0,a11111=0.0; FILE *fp=fopen(RESULT.txt,wt);/由于需要显示的内容过多而需要产生一个txt文件来全部显示/*以下为对矩阵A的初始化

17、赋值*/ for(i=1;i=n;i+) for(j=1;j=n;j+) if(i=j) aij = 1.5*cos(i+1.2*j); Pij = 1.0; else aij = sin(0.5*i+0.2*j); /*调用函数upper_triangular_matrix(a),对A拟上三角化*/ upper_triangular_matrix(a); /*输出拟上三角矩阵A(n-1)*/ fprintf(fp,n拟上三角化矩阵A(n-1)n); printf(n*数值分析大作业二*n); printf(n*拟上三角化矩阵A(n-1)*n);for(i=1;i=n;i+)for (int

18、j=1;j=n;j+)a1ij=aij;fprintf(fp,A%2d%2d=%.12en,i,j,aij);if(j+1)%2=0) printf(n);printf( A%2d%2d=%.12e ,i,j,aij);fprintf(fp,n);printf(n);/*执行带双步位移QR分解，返回非零值时计算成功*/ flag=twostep(a,Vector,fp); /*输出对A(n-1)进行带双步位移的QR分解后的矩阵*/ printf(n矩阵A(n-1)进行QR分解后所得矩阵n); fprintf(fp,n矩阵A(n-1)进行QR分解后所得矩阵n);for(i=1;i=n;i+)fo

19、r (int j=1;j=n;j+)fprintf(fp,A%2d%2d=%20.12en,i,j,aij);printf(A%2d%2d=%20.12en,i,j,aij);fprintf(fp,n);printf(n);/*输出矩阵A的特征值*/ if(flag) fprintf(fp,n矩阵A的所有特征值nRi Ii 复数标志量n); printf(n矩阵A的所有特征值nRi Ii 复数标志量n); for(i=1;i=n;i+) printf(%.12e %.12e %6gn,Vectori0,Vectori1,Vectori2);if(j+1)%2=0) printf(n); fprintf(fp,%.12e %.12e %.12e ,Vectori0,Vectori1,Vectori2); else printf(n特征值求解失败！n); /*调用函数求解实特征值对应特征向量*/ EigenVector(a1,x,Vector,fp); fclose(fp); 三、矩阵A经过拟上三角化后所得的矩阵四、对矩阵进行QR分解方法结束后所得的矩阵Q,R以及乘积矩阵RQ.五、矩阵A的全部特征值六、A的相应于实特征值的特征向量