数值分析第二次作业

1、数值分析第二次上机作业算法设计1. 矩阵A的拟上三角化对矩阵A 进行拟上三角化按照教材3.3.2给出的具体算法即可求出。2. 对A(n-1)执行QR方法第一步，并打印第一步做完后的矩阵RQ根据教材3.3.1给出的QR分解的具体算法，通过加入一个控制具体执行QR分解进行的步数的参数即可实现对A(n-1)执行QR方法第一步，将得到的R矩阵和Q矩阵相乘可得到矩阵RQ，将A(n-1)和RQ矩阵相减可得到两个矩阵的差。3. 用双步位移的QR 方法计算矩阵A的全部特征值根据教材3.3.3给出的用带双步位移的QR方法的求实矩阵A的全部特征值的具体算法，就可以求出所有的特征值。4. A的相应于实特征值的特征向

2、量求出矩阵A的所有特征值后，由公式AX=X可得(I-A)X=0，通过求解该方程组即可得到实特征值的特征向量。由于方程组的b向量为0，通过列主元素的Gauss消去法就可以快速的求解出向量的解。因为方程组的系数矩阵为奇异矩阵，需要给X向量的最后一个分量赋一个非零的值才能得到非零解。源程序(在CodeBlock C/C+集成开发环境下编译通过)#include #include using namespace std;#define MAXLEN 10#define THTA 1E-12#define ITIME 10000/矩阵的复数特征值typedef struct comp double re

3、al; double imag;Comp;/初始化矩阵void createMatrix(double AMAXLEN+1) for(int i=1;i=MAXLEN;i+) for(int j=1;j=MAXLEN;j+) if(i!=j) Aij=sin(0.5*i+0.2*j); else Aij=1.52*cos(i+1.2*j); /打印矩阵void printMatrix(double AMAXLEN+1) for(int i=1;i=MAXLEN;i+) for(int j=1;j=MAXLEN;j+) coutAij ; cout0) return 1; if(value=0)

4、 return -1;/矩阵的拟上三角化double (*triangulation(double AMAXLEN+1)MAXLEN+1 double (*tempA)MAXLEN+1=new doubleMAXLEN+1MAXLEN+1; for(int i=1;i=MAXLEN;i+) for(int j=1;j=MAXLEN;j+) tempAij=Aij; for(int r=1;r=MAXLEN-2;r+) bool canContinue=true; for(int i=r+2;i=MAXLEN;i+) if(tempAir!=0) canContinue=false; break

5、; if(canContinue) continue; double dr=0,cr=0,hr=0; for(int i=r+1;i=MAXLEN;i+) dr+=tempAir*tempAir; dr=sqrt(dr); cr=-sgn(tempAr+1r)*dr; hr=cr*cr-cr*tempAr+1r; double urMAXLEN+1,prMAXLEN+1,qrMAXLEN+1,wrMAXLEN+1,tr; for(int i=1;ir+1;i+) uri=0; for(int i=r+1;i=MAXLEN;i+) uri=tempAir; urr+1-=cr; tr=0; fo

6、r(int i=1;i=MAXLEN;i+) pri=0; qri=0; for(int j=1;j=MAXLEN;j+) pri+=tempAji*urj; qri+=tempAij*urj; pri=pri/hr; qri=qri/hr; tr+=pri*uri; tr=tr/hr; for(int i=1;i=MAXLEN;i+) wri=qri-tr*uri; for(int i=1;i=MAXLEN;i+) for(int j=1;j=MAXLEN;j+) tempAij=tempAij-wri*urj-uri*prj; for(int i=1;i=MAXLEN;i+) for(in

7、t j=1;j=MAXLEN;j+) if(fabs(tempAij)THTA) tempAij=0; return tempA;/矩阵相乘double (*matrixMut(double AMAXLEN+1,double BMAXLEN+1)MAXLEN+1 double (*result)MAXLEN+1=new doubleMAXLEN+1MAXLEN+1; for(int i=1;i=MAXLEN;i+) for(int j=1;j=MAXLEN;j+) resultij=0; for(int k=1;k=MAXLEN;k+) resultij+=Aik*Bkj; return re

8、sult;/矩阵转置double (*matrixTrans(double AMAXLEN+1)MAXLEN+1 double (*result)MAXLEN+1=new doubleMAXLEN+1MAXLEN+1; for(int i=1;i=MAXLEN;i+) for(int j=1;j=MAXLEN;j+) resultji=Aij; return result;/矩阵复制double (*matrixCopy(double AMAXLEN+1)MAXLEN+1 double (*result)MAXLEN+1=new doubleMAXLEN+1MAXLEN+1; for(int

9、i=1;i=MAXLEN;i+) for(int j=1;j=MAXLEN;j+) resultij=Aij; return result;/QR方法void QRMethod(double AMAXLEN+1,double (*&Q)MAXLEN+1,double (*&R)MAXLEN+1,int step) Q=new doubleMAXLEN+1MAXLEN+1; R=new doubleMAXLEN+1MAXLEN+1; for(int i=1;i=MAXLEN;i+) for(int j=1;j=MAXLEN;j+) if(i=j) Qij=1; else Qij=0; Rij=A

10、ij; for(int r=1;r=step;r+) bool canContinue=true; for(int i=r+1;i=MAXLEN;i+) if(Rir!=0) canContinue=false; break; if(canContinue) continue; double dr=0,cr=0,hr=0; for(int i=r;i=MAXLEN;i+) dr+=Rir*Rir; dr=sqrt(dr); cr=-sgn(Rrr)*dr; hr=cr*cr-cr*Rrr; double urMAXLEN+1,wrMAXLEN+1,prMAXLEN+1,tr; for(int

11、i=1;ir;i+) uri=0; for(int i=r+1;i=MAXLEN;i+) uri=Rir; urr=Rrr-cr; for(int i=1;i=MAXLEN;i+) wri=0; for(int j=1;j=MAXLEN;j+) wri+=Qij*urj; for(int i=1;i=MAXLEN;i+) for(int j=1;j=MAXLEN;j+) Qij=Qij-wri*urj/hr; for(int i=1;i=MAXLEN;i+) pri=0; for(int j=1;j=MAXLEN;j+) pri+=Rji*urj; pri=pri/hr; for(int i=

12、1;i=MAXLEN;i+) for(int j=1;j=MAXLEN;j+) Rij=Rij-uri*prj; /双步位移的QR方法bool twoStepQR(double AMAXLEN+1,int L,Comp flags) int k=1,m=MAXLEN;/第二步 int num=1; while(true) if(fabs(Amm-1)=THTA)/第三步 flagsnum.real=Amm; flagsnum+.imag=0; m=m-1; if(m=1)/第四步 flagsnum.real=Amm; flagsnum+.imag=0; return true; if(m1)

13、continue; else double a=1;/第五步 double b=-(Am-1m-1+Amm); double c=Am-1m-1*Amm-Am-1m*Amm-1; double delta=b*b-4*a*c; Comp s1,s2; if(delta=0) s1.real=(-b-sqrt(delta)/(2*a); s1.imag=0; s2.real=(-b+sqrt(delta)/(2*a); s2.imag=0; else s1.real=(-b)/(2*a); s1.imag=-sqrt(-delta)/(2*a); s2.real=(-b)/(2*a); s2.i

14、mag=sqrt(-delta)/(2*a); if(m=2)/第六步 flagsnum+=s1; flagsnum+=s2; return true; else if(fabs(Am-1m-2)1) continue; else if(k=L)/第八步 return false; else double s=Am-1m-1+Amm;/第九步 double t=Am-1m-1*Amm-Am-1m*Amm-1; double (*MK)MAXLEN+1=new doubleMAXLEN+1MAXLEN+1; double (*AK)MAXLEN+1=new doubleMAXLEN+1MAXLE

15、N+1; for(int i=1;i=MAXLEN;i+) for(int j=1;jm)|(jm) AKij=0; else AKij=Aij; double (*AK2)MAXLEN+1=matrixMut(AK,AK); for(int i=1;i=MAXLEN;i+) for(int j=1;jm)|(jm) MKij=0; continue; MKij=AK2ij-s*AKij; if(i=j) MKij+=t; double (*Q)MAXLEN+1,(*R)MAXLEN+1,(*QT)MAXLEN+1,(*temp)MAXLEN+1; QRMethod(MK,Q,R,MAXLEN

16、-1); QT=matrixTrans(Q); temp=matrixMut(QT,AK); delete A; A=matrixMut(temp,Q); delete AK2; delete MK; delete Q; delete R; delete QT; delete temp; delete AK; k+; /解线性方程组void solveEqua(double(*mat)MAXLEN+1,double x) for(int k=1;kMAXLEN;k+) int i=k; for(int j=k;jmatik) i=j; for(int j=k;j=MAXLEN;j+) doub

17、le temp=matkj; matkj=matij; matij=temp; for(int i=k+1;i=MAXLEN;i+) double mik=matik/matkk; for(int j=k+1;j=1;k-) xk=0; for(int j=k+1;j=MAXLEN;j+) xk-=matkj*xj; xk=xk/matkk; int main() double AMAXLEN+1MAXLEN+1,(*ssjResult)MAXLEN+1; createMatrix(A); cout.precision(11); cout.setf(ios:scientific|ios:upp

18、ercase); ssjResult=triangulation(A); cout矩阵A的拟上三角化endl; printMatrix(ssjResult); double (*Q)MAXLEN+1=NULL,(*R)MAXLEN+1=NULL; QRMethod(ssjResult,Q,R,1); coutendl; cout拟上三角矩阵A的QR分解第一步endl; coutR*Qendl; double (*QR)MAXLEN+1; QR=matrixMut(R,Q); printMatrix(QR); Comp flagsMAXLEN+1; coutendl; coutA(n-1)和R*

19、Q的差值endl; for(int i=1;i=MAXLEN;i+) for(int j=1;j=MAXLEN;j+) coutssjResultij-QRij ; coutendl; coutendl; if(twoStepQR(matrixCopy(ssjResult),ITIME,flags) for(int i=1;i=MAXLEN;i+) couti=flagsi.real0) cout+flagsi.imag; else if(flagsi.imag0) coutflagsi.imag; coutendl; coutendl; for(int i=1;i=MAXLEN;i+) if

20、(flagsi.imag=0) double tempMAXLEN+1MAXLEN+1; double xMAXLEN+1; for(int j=1;j=MAXLEN;j+) for(int k=0;k=MAXLEN;k+) tempjk=-Ajk; if(j=k) tempjk+=flagsi.real; solveEqua(temp,x); cout实特征值i=flagsi.real的特征向量endl; for(int j=1;j=MAXLEN;j+) coutxj ; coutendlendl; else coutFail!endl; return 0;程序输出1. 矩阵A的拟上三角化2

21、. 对A(n-1)执行QR方法第一步，并打印第一步做完后的矩阵RQ3. 用双步位移的QR 方法计算矩阵A的全部特征值4. A的相应于实特征值的特征向量计算结果1. 矩阵A的拟上三角化矩阵A得到的拟上三角化矩阵为：-8.E-001-9.E-002-1.E+000 -7.E-001 1.E-001-1.E+000 -8.E-002 9.E-001-6.E-0011.E-001-2.E+000 2.E+000 1.E+0003.E-001 2.E-001 2.E+000 1.E-001-1.E+000 6.E-001 -4.E-0010.E+000 1.E+000-1.E+000 -1.E+000-

22、6.E-001-1.E+000 2.E-001 6.E-001-1.E-001 3.E-0010.E+000 0.E+000-1.E+000 -1.E+000 1.E+000-1.E+000 1.E-001 4.E-001-1.E-001 1.E-0010.E+000 0.E+000 0.E+000 1.E+000 8.E-001 4.E-001 -4.E-002-4.E-001 2.E-001 -1.E-0010.E+000 0.E+000 0.E+000 0.E+000-7.E-001-1.E+000 -2.E-001-2.E-001 6.E-001 2.E-0010.E+000 0.E+

23、000 0.E+000 0.E+000 0.E+000-7.E-001 -7.E-003 9.E-001-1.E-001 2.E-0020.E+000 0.E+000 0.E+000 0.E+000 0.E+000 0.E+000 7.E-001-4.E-001 5.E-001 -1.E-0010.E+000 0.E+000 0.E+000 0.E+000 0.E+000 0.E+000 0.E+000 7.E-001 1.E-001 -3.E-0010.E+000 0.E+000 0.E+000 0.E+000 0.E+000 0.E+000 0.E+000 0.E+000-4.E-001

24、4.E-0012. 对A(n-1)执行QR方法第一步，并打印第一步做完后的矩阵RQ执行QR方法的第一步的到的R*Q矩阵为:1.E+000-3.E+000-1.E+000-3.E-002-2.E-001-1.E+00-1.E-001 8.E-001-4.E-0013.E-001-8.E-001 3.E-001 1.E+0008.E-001-8.E-002 2.E+0001.E-001-1.E+000 8.E-001-3.E-001-1.E+000 6.E-001-1.E+000-1.E+000-6.E-001-1.E+002.E-001 6.E-001-1.E-0013.E-0010.E+000

25、 0.E+000-1.E+000-1.E+000 1.E+000-1.E+0001.E-001 4.E-001-1.E-0011.E-0010.E+000 0.E+000 0.E+0001.E+000 8.E-001 4.E-001-4.E-002-4.E-001 2.E-001-1.E-0010.E+000 0.E+000 0.E+0000.E+000-7.E-001-1.E+000-2.E-001-2.E-001 6.E-0012.E-0010.E+000 0.E+000 0.E+0000.E+000 0.E+000-7.E-001-7.E-003 9.E-001-1.E-0012.E-0

26、020.E+000 0.E+000 0.E+0000.E+000 0.E+000 0.E+0007.E-001-4.E-001 5.E-001-1.E-0010.E+000 0.E+000 0.E+0000.E+000 0.E+000 0.E+0000.E+000 7.E-001 1.E-001-3.E-0010.E+000 0.E+000 0.E+0000.E+000 0.E+000 0.E+0000.E+000 0.E+000-4.E-0014.E-0013. 用双步位移的QR 方法计算矩阵A的全部特征值1=4.E-0022=6.E-0013=9.E-0014=-9.E-001 -1.E-

27、001i5=-9.E-001 +1.E-001i6=-1.E+0007=1.E+0008=-2.E+000 -8.E-001i9=-2.E+000 +8.E-001i10=3.E+0004. A的相应于实特征值的特征向量实特征值1=4.E-002的特征向量-4.E+000-4.E+000 8.E+000-7.E-001-8.E+000-2.E+0001.E+001-7.E+000-6.E+0001.E+000实特征值2=6.E-001的特征向量4.E+000 2.E+000 1.E+001-1.E+000-2.E+001 7.E+000-9.E+000 1.E+001 1.E+0011.E+000实特征值3=9.E-001的特征向量2.E+000 1.E+000-6.E-001-1.E+000-1.E+001 7.E+000-5.