数值分析大作业二.doc

《数值分析》大作业二 一、算法的设计 总体思路：采用带双步位移的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. 对于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∈化为拟上三角矩阵；给定精度水平E，以及最大的迭代次数L=1000。 (2)记==，令k=1，m=n。 (3)如果，则得到A的一个特征值，置m：=m-1，转(4)，否则转(5)。 (4)如果m=1，则得到A 的一个特征值，转(11)；如果m=0，则直接转(11)；如果m1，则转(3)。 (5)求二阶子阵 的两个特征值和，即计算二次方程的两个根和。 (6)如果m=2，则得到A的两个特征值和，转(11)，否则转(7)。 (7)如果，则得到A的两个特征值和，置m：=m-2，转(4)，否则则转(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为上三角阵，令特征向量的最后一个元素为1，然后再进行迭代求解。这样，就求出了所有的值。 二、全部源程序 #includestdio.h #includemath.h #define n 10 #define E 1.0e-12 //***精确度为***// #define L 1000 //*******定义迭代次数为1000次*******// int twostep(double A[11][1+1],double Vector[][3],FILE *fp); void upper_triangular_matrix(double A[11][11]); void QR (double C[][11],double B[][11],int m); void solve_equation(double s1[],double s2[],double Vector[][3],int m); void Gauss(double a1[][11],double x[][111],double Vector[][3],int p); void EigenVector(double A[][11],double x[][11],double Vector[][3],FILE *fp); //*********************矩阵A的拟上三角化*********************// void upper_triangular_matrix(double A[11][11]) { int r,k,j,m; double dr,cr,tr,hr; double judge; double sum, sum1,sum2,sum3; double ur[11],pr[11],qr[11],wr[11]; for(r=1;r=n-2;r++) { judge=0; for (k=r+2;k=n;k++) judge=judge+A[k][r]*A[k][r]; if(judge==0) continue; else { sum=0; f