北航数值分析大作业第二次QR分解.docVIP

目标：使用带双步位移的QR分解法求矩阵的全部特征值，并对其中的每一个实特征值求相应的特征向量。已知： i,j 1,2,……,10 算法： 输入需要求解的矩阵 首先输入需要求解的矩阵。即为下述程序中的void creaA 子程序。 对上述生成的矩阵进行拟上三角化 为了减少求特征值和特征向量过程中的计算量，在对矩阵进行QR分解前先进行拟上三角化。由于矩阵的QR分解不改变矩阵的结构，所以具有拟上三角形状的矩阵的QR分解可以大大极少计算量。在下述程序中用void hess 子程序来实现。 对拟上三角化后的矩阵进行QR分解。 为了直观的了解普通的QR分解过程及结果，下述程序中用void AQR（）子程序来对拟上三角化过后的A阵进行QR分解，并输出Q阵R阵RQ阵。 对拟上三角化后的矩阵进行带双步位移的QR分解。 子程序void QRmeth（）实现对拟上三角化后的A阵进行带双步位移的QR分解，得出特征值并输出，并用子程序void gauss（）对其中的实数特征值进行求解，得出对应的特征向量。 输出相关结果。 源程序： #include #include #include #define L 2500 #define n 11 #define e 0.000000000001 int i,j,s,p,k,ik; double A[n][n],q[n][n],r[n][n],rq[n][n],I[n][n]; double P[n],W[n],u[n],Q[n]; double dr,cr,hr,ar,tr; int nR,nC; double s1,t,x,tzR[n],tzC[2][n],sum,M[n][n],v[n]; void creaA ////生成矩阵A for i 1;i 11;i++ for j 1;j 11;j++ if j! i A[i][j] sin 0.5*i+0.2*j ; else A[i][j] 1.5*cos i+1.2*j ; void hess for s 1;s n-2;s++ for ar 0.0,i s+2;i n;i++ ar+ A[i][s]*A[i][s];//printf ar %1.12e\n,ar ; if ar 0 continue; else ar+ A[s+1][s]*A[s+1][s]; dr sqrt ar ;//printf dr %1.12e\n,dr ; if A[s+1][s] 0 cr -dr; else cr dr; hr cr*cr-cr*A[s+1][s]; for i 1;i s;i++ ////生成u向量 u[i] 0.0; u[s+1] A[s+1][s]-cr; for i s+2;i n;i++ u[i] A[i][s]; for j 1;j n;j++ //P for P[j] 0.0,i 1;i n;i++ P[j]+ A[i][j]*u[i]/hr; for tr 0.0,i 1;i n;i++ //tr tr+ P[i]*u[i]/hr; for i 1;i n;i++ //Q for Q[i] 0.0,j 1;j n;j++ Q[i]+ A[i][j]*u[j]/hr; for i 1;i n;i++ //W W[i] Q[i]-tr*u[i]; for i 1;i n;i++ //生成Ar+1 for j 1;j n;j++ A[i][j]- W[i]*u[j]+u[i]*P[j] ; for i 1;i n;i++ for j 1;j n;j++ if fabs A[i][j] e A[i][j] 0.0; printf 拟上三角化后A n-1 :\n ; for i 1;i n;i++ for j 1;j n;j++ printf %1.12e,,A[i][j] ; printf \n ; void AQR double u[n],w[n],F[n]; for s 1;s n;s++ for p 1;p n;p++ r[s][p] A[s][p]; for s 1;s n-1;s++ for dr 0.0,i s;i n;i++ dr+ r[i][s]*r[i][s]; dr sqrt dr ; if dr 0 continue; else if A[s][s] 0 cr -dr; else cr dr; hr cr*cr-cr*r[s][s]; for i 1;i s;i++ //u u[i] 0; u[s] r[s][s]-cr; for i s+1;i n;i++ u[i] r[i][s]; for i 1;i n;i++ //F for F[i] 0.0,j 1;j n;j++ F[