北航数值分析第二次大作业

《北航数值分析第二次大作业》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

1、数值分析第二次大作业一、算法的设计方案首先构造矩阵A，利用Householder矩阵对矩阵A作相似变换，把A化为拟上三角矩阵A(n-1)，其算法如书上第61页所示。使用QR分解法对矩阵A(n-1)进行QR分解，其算法如书上第59页所示，进而可以求出所得矩阵的Q、R、RQ阵。然后对A(n-1)进行带双步位移的QR分解求解矩阵的全部特征值。这里采用这么几步进行：第一步：判断是否，若不是，则进入第四步。若是，则为特征值，m=m-1，若m=1，进入第二步，若m=2进入第三步，否则进入第四步。第二步：m=1，则为特

2、征值，转向结束步。第三步：m=2，则可以求出A的两个特征值和，转向结束步。第四步：判断是否，若不是，进入第五步。若是，则得到A的两个特征值和，令m=m-2，若m=1，进入第二步，若m=2进入第三步，否则进入第一步。第五步：判断是否达到循环上限，若达到上限，则结束，否则进入第六步。第六步：对A进行双步移位QR分解，这里的算法如书上第64页所示，分解之后转向计数步。计数步：对循环次数进行计数，并转向第一步。结束步：显示所求得的特征值。最后对实特征值利用列主元高斯消元法求解其对应的特征向量，其算法见书第17页。

3、二、对程序的几点说明1、对A求拟上三角化矩阵A(n-1)，这里采用一个单独的函数NSSjiao()来求出。2、对矩阵A(n-1)进行QR分解采用一个单独的函数QRfenjie()求出。1、对矩阵A(n-1)进行带双步位移的QR分解采用一个单独的函数SBQR()求出。2、对实特征值利用列主元高斯消元法求解其对应的特征向量在这里采用一个单独的函数vector()来实现。3、所有的结果都用C语言文件指令使之显示在一个txt文件中，方便结果的复制与粘贴。A(n-1)输出到NSSjiao.txt中，A(n-1)QR

4、分解所得矩阵的Q、R、RQ阵输出到q.txt中，A的特征值输出到tezhenzhi.txt中，A实特征值对应的特征向量输出到vector.txt中。三、源程序如下：#include"stdio.h"#include"math.h"#defineN10+1#definen10#defineL1000doublea[N][N];doubleRe[N]={0};doubleIm[N]={0};intsgn(doublea);//符号函数voidNSSjiao();//拟上三角化函数voidQRfenjie();

5、//QR分解函数voidSBQR();//双步QR分解函数voidvector();//求特征向量的函数voidmain(){inti,j;for(i=1;i<=n;i++){for(j=1;j<=n;j++){if(i!=j)a[i][j]=sin(0.5*i+0.2*j);elsea[i][j]=1.5*cos(i+1.2*j);}}//赋值，矩阵初始化NSSjiao();QRfenjie();SBQR();vector();}intsgn(doublea){if(a>=0)return(1);els

6、ereturn(-1);}voidNSSjiao(){intr,i,j;intcnt=0;doubledr=0;doublecr=0;doublehr=0;doubletr=0;doubleur[N]={0};doublepr[N]={0};doubleqr[N]={0};doublewr[N]={0};for(r=1;r<=n-2;r++){for(i=r+2;i<=n;i++){if(a[i][r]==0)cnt++;}if(cnt!=n-r-1){for(i=r+1;i<=n;i++)dr+=a[i

7、][r]*a[i][r];dr=sqrt(dr);cr=-sgn(a[r+1][r])*dr;hr=cr*cr-cr*a[r+1][r];for(i=1;i<=r;i++)ur[i]=0;for(i=r+2;i<=n;i++)ur[i]=a[i][r];ur[r+1]=a[r+1][r]-cr;for(i=1;i<=n;i++)for(j=1;j<=n;j++)pr[i]+=a[j][i]*ur[j]/hr;for(i=1;i<=n;i++)for(j=1;j<=n;j++)qr[i]+=a[i][j]*

8、ur[j]/hr;for(i=1;i<=n;i++)tr+=pr[i]*ur[i]/hr;for(i=1;i<=n;i++)wr[i]=qr[i]-tr*ur[i];for(i=1;i<=n;i++)for(j=1;j<=n;j++)a[i][j]-=(wr[i]*ur[j]+ur[i]*pr[j]);dr=0;cr=0;hr=0;tr=0;for(i=1;i<=n;i++){pr[i]=0;qr[i]=0;}}}FILE*f