北航数值分析大作业2..docVIP

一、算法设计方案计算矩阵A对i,j=1到10执行i,j=1, …,n。设置标志量flag=0，flag=0表示A[i][r]全为零，flag=1表示A[i][r]不全为零。对r=1到n-2执行 (1) if |air|Epsilon,flag=1,转(2)； (2) 计算 dr=sqrt(A[i][r]* A[i][r]) cr=--sgn(A[r+1][r]) hr = cr*(cr-A[r+1][r]); (3) 对i=r+1到n执行ur[i] = A[i][r]; ur[r+1] = ur[r+1]-cr； (4) 计算 pr=A’*ur/hr qr=A*ur/hr tr=pr’*ur/hr wr=qr-tr*ur A=A-wr*ur-ur*pr’ 3.对A(n-1)做双步位移QR分解 (1) 记A1=A(n-1)=A[i][j]，令k=0,m=n; (2) if | A[m][m-1]=Epsilon|,则得到A的一个特征值Lamda[m][0]=A[m][m]，置m=m-1,转(3);否则转(4); (3) if m=1,则得到A的一个特征值Lamda[m][0]=A[1][1]，转(10);if m=0，则直接转(10)；if m1,则转(2); (4) 求二阶子阵的两个特征值，即计算二次方程的两个根， x*x- s1[0]*x+ s2[0]=0 （*） 其中s1[0] = A[m-1][m-1]+A[m][m]; s2[0]= A[m-1][m-1]*A[m][m]-A[m-1][m]*A[m][m-1] （该过程可用求解二次方程子程序实现，注意讨论根的情况）； (5) if m=2,则得到A的两个特征值，转(10)；否则转(7)； (6) if |A[m-1][m-2]) |=Epsilon,则得到A的两个特征根，置m=m-2，转(3);否则转(7)； (7) if k=L,printf “Fail to get all the eigenvalues of matrix A!”；否则转(8)； (8) 对i,j=1,…,m，计算 s= A[m-1][m-1]+A[m][m]; t= A[m-1][m-1]*A[m][m]-A[m-1][m]*A[m][m-1]; Mk[i][j] = Mk[i][j]-s*A[i][j]+t*I[i][j] Mk=Qk*Rk（对Mk作QR分解，可调用QR分解子程序） Ak+1=Qk’*Ak*Qk (9) 置k=k+1，转(2)； (10) printf “Get all the eigenvalues of matrix A! Congratulations!”，停止计算。 Mk的QR分解算法： 记B[][N]=Mk[][N],C[][N]=A[][N]，置flag=0(B[i][r]=0)。对r=1到m-1执行 (1) if |B[i][r]|=Epsilon,flag=1，转(2); (2) 计算 dr=sqrt(A[i][r]* A[i][r]) cr=--sgn(A[r+1][r]) hr = cr*(cr-A[r+1][r]); (3) 对i=r+1到m执行ur[i] = A[i][r]; ur[r+1] = ur[r+1]-cr； (4)计算 vr=Br’*ur/hr Br+1=Br-ur*vr’ pr=Crur/hr qr=Crr/hr tr=pr’*ur/hr wr=qr-tr*ur Cr+1=Cr-wr*ur-ur*pr’ (5) continue. 求二元一次方程算法： 置根Lamda[][3]={0}。对方程（*），令t1=s1[0];t2=s2[0]; 计算 delt = t1*t1-4*t2; sqdelt = sqrt(fabs(delt)); (2) if delt=0.0，则s1[0] = (t1+sqdelt)/2.0，s2[0] = (t1-sqdelt)/2.0;否则转(3); (3)置Lamda[m][2]=1; Lamda[m-1][2]=1; s1[0]=t1/2.0; s1[1]=sqdelt/2.0; s2[0]=t1/2.0; s2[1]=-sqdelt/2.0; 4.用列主元Gauss消去法求解实特征根对应的特征向量 (1) if Lamda[m][2]=0，转(2); (2) 调用GaussOptimal子程序； (3) printf “When the eigenvalue is %4.12e, the corresponding eigenvector is\n,Lamda[i][0]” GaussOptimal算法 记a1=aij，令b[N]