北航数值分析大作业第二题Word格式.docx
《北航数值分析大作业第二题Word格式.docx》由会员分享，可在线阅读，更多相关《北航数值分析大作业第二题Word格式.docx（26页珍藏版）》请在冰豆网上搜索。
2.记，令。
3.如果，则得到的一个特征值，置（降阶），转4；
否则转5。
4.如果，则得到的一个特征值，转11；
如果，则转3。
5.求2阶子阵
的两个特征值和，即计算二次方程
的两个根和。
6.如果，则得到的两个特征值和，转11；
否则转7。
7.如果，则得到的两个特征值和，置（降阶），转4；
否则转8
8.如果，则计算终止，未得到的全部特征值；
否则转9。
9.记，计算
10.置，转3。
11.的全部特征值已计算完毕，停止计算。
其中，的分解与的计算用下列算法实现：
记。
1.若全为零，则令,转5；
此算法执行完后，就得到。
（4）计算Q,R，一边求R*Q矩阵。
记
1.若全为零，则令转5；
否则转2.
2.计算
3.令。
4.计算
5.继续
当此算法执行完毕后，就得到正交矩阵和上三角矩阵
6.然后计算出矩阵
（5）用列主元素Gauss消去法计算矩阵对应于实特征值的特征向量，具体算法如下：
1.消元过程
（1）选行号，使。
（2）交换与所含的数值。
（3）对于计算
2.回代过程
最终得到的向量的即为对应于实特征值的特征向量。
二、源程序
#include<
iostream.h>
math.h>
iomanip.h>
#definen10
#defineE1.0e-12
voidHouseholder（doublea[n][n]）;
//拟上三角化函数
doublesgn（doublea）;
//符号函数
voidQRfenjie（doublea[n][n],doubleQ[n][n],doubleR[n][n]）;
voidQR（doublea[n][n],doubleL[n][2]）;
//带双步位移的QR分解
voidMxM（doubleM[n][n],doubleA[n][n],doubleB[n][n],intm）;
//矩阵相乘
voidsolve（doublea[n][n],doubles1[2],doubles2[2],intm）;
//解方程函数
voidGauss（doublea[n][n],doublex[n]）;
//定义列主元高斯消去法函数
voidmain（）
{
inti,j,k;
doublea[n][n],qr[n][n],q[n][n],r[n][n],L[n][2],x[n];
for（i=0;
i<
n;
i++）//录入要进行变换的矩a，并对q初始赋值。
{
for（j=0;
j<
j++）
{
if（i==j）
a[i][j]=1.5*cos（i+1+1.2*（j+1））;
else
a[i][j]=sin（0.5*（i+1）+0.2*（j+1））;
}
}
//cout<
<
endl;
Householder（a）;
//调用拟上三角化函数
cout<
endl<
"
对矩阵A拟上三角化前七列的结果:
i++）
//k=0;
//为了显示方便，每行显示三个元素，使用k来实现
//cout<
矩阵A的第"
i+1<
行元素为：
7;
if（fabs（a[i][j]）<
E）
a[i][j]=0;
cout<
setprecision（12）<
setiosflags（ios:
:
scientific）<
a[i][j]<
"
;
//k++;
//cout<
//if（k%3==0）//判断每行是否到达三个元素
cout<
对矩阵A拟上三角化后三列的结果:
for（j=7;
QRfenjie（a,q,r）;
QR（a,L）;
对矩阵A进行QR分解后所得矩阵前七列的结果:
对矩阵A进行QR分解后所得矩阵三列的结果:
/*cout<
对矩阵A进行QR分解后所得矩阵:
k=0;
k++;
if（k%3==0）
cout<
}*/
j++）
for（k=0,qr[i][j]=0;
k<
k++）
qr[i][j]=qr[i][j]+r[i][k]*q[k][j];
R*Q相乘的前七列:
R*Q的第"
if（fabs（qr[i][j]）<
qr[i][j]<
//k++;
//if（k%3==0）
//cout<
R*Q相乘的后三列:
矩阵A的全部特征值为"
if（L[i][1]==0）
λ"
="
L[i][0]<
else
+"
L[i][1]<
i"
实特征值对应的特征向量分别是：
for（k=0;
k++）
if（L[k][1]==0）//判断特征值是否为实特征值
for（i=0;
i++）//重新录入矩阵A
{
for（j=0;
{
if（i==j）
a[i][j]=1.5*cos（i+1+1.2*（j+1））-L[k][0];
else
a[i][j]=sin（0.5*（i+1）+0.2*（j+1））;
}
}
Gauss（a,x）;
k+1<
L[k][0]<
对应的特征向量是:
for（j=0;
x[j]<
elsecontinue;
}
}
voidHouseholder（doublea[n][n]）//拟上三角化函数
intm=0;
doubled,c,h,t;
doubleu[n],p[n],q[n],w[n];
n-2;
for（j=i+2;
if（a[j][i]<
=E）
m=m+1;
if（m==（n-2-i））
continue;
for（j=i+1,d=0;
d=d+a[j][i]*a[j][i];
d=sqrt（d）;
c=-1*sgn（a[i+1][i]）*d;
h=c*c-c*a[i+1][i];
for（j=i+2;
u[j]=a[j][i];
升级会员 