北航数值分析大作业第二题.docx
《北航数值分析大作业第二题.docx》由会员分享，可在线阅读，更多相关《北航数值分析大作业第二题.docx（26页珍藏版）》请在冰豆网上搜索。
北航数值分析大作业第二题
数值分析第二次大作业
史立峰
SY1505327
一、方案
（1）利用循环结构将（i,j=1,2,……,10）进行赋值，得到需要变换的矩阵A；
（2）然后，对矩阵A利用Householder矩阵进行相似变换，把A化为上三角矩阵A（n-1）。
对A拟上三角化，得到拟上三角矩阵A（n-1），具体算法如下：
记A
（1）=A，并记A（r）的第r列至第n列的元素为。
对于执行
1.若全为零，则令A（r+1）=A（r）,转5；否则转2。
2.计算
3.令。
4.计算
5.继续。
（3）使用带双步位移的QR方法计算矩阵A（n-1）的全部特征值，也是A的全部特征值，具体算法如下：
1.给定精度水平和迭代最大次数。
2.记，令。
3.如果，则得到的一个特征值，置（降阶），转4；否则转5。
4.如果，则得到的一个特征值，转11；如果，则转3。
5.求2阶子阵
的两个特征值和，即计算二次方程
的两个根和。
6.如果，则得到的两个特征值和，转11；否则转7。
7.如果，则得到的两个特征值和，置（降阶），转4；否则转8
8.如果，则计算终止，未得到的全部特征值；否则转9。
9.记，计算
10.置，转3。
11.的全部特征值已计算完毕，停止计算。
其中，的分解与的计算用下列算法实现：
记。
对于执行
1.若全为零，则令,转5；否则转2。
2.计算
3.令。
4.计算
5.继续。
此算法执行完后，就得到。
（4）计算Q,R，一边求R*Q矩阵。
记
对于执行
1.若全为零，则令转5；否则转2.
2.计算
3.令。
4.计算
5.继续
当此算法执行完毕后，就得到正交矩阵和上三角矩阵
6.然后计算出矩阵
（5）用列主元素Gauss消去法计算矩阵对应于实特征值的特征向量，具体算法如下：
记
1.消元过程
对于执行
（1）选行号，使。
（2）交换与所含的数值。
（3）对于计算
2.回代过程
最终得到的向量的即为对应于实特征值的特征向量。
二、源程序
#include
#include
#include
#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{
for（j=0;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<Householder（a）;//调用拟上三角化函数
cout<"<for（i=0;i{
//k=0;//为了显示方便，每行显示三个元素，使用k来实现
//cout<<"矩阵A的第"<"<for（j=0;j<7;j++）
{
if（fabs（a[i][j]）a[i][j]=0;
cout<:
scientific）<//k++;
//cout<//if（k%3==0）//判断每行是否到达三个元素
}
cout<}
cout<<"对矩阵A拟上三角化后三列的结果:
"<for（i=0;i{
//k=0;//为了显示方便，每行显示三个元素，使用k来实现
//cout<<"矩阵A的第"<"<for（j=7;j{
if（fabs（a[i][j]）a[i][j]=0;
cout<:
scientific）<//k++;
//if（k%3==0）//判断每行是否到达三个元素
//cout<}
cout<}
cout<QRfenjie（a,q,r）;
QR（a,L）;
cout<"<for（i=0;i{
for（j=0;j<7;j++）
{
if（fabs（a[i][j]）a[i][j]=0;
cout<:
scientific）<}
cout<}
cout<<"对矩阵A进行QR分解后所得矩阵三列的结果:
"<for（i=0;i{
for（j=7;j{
if（fabs（a[i][j]）a[i][j]=0;
cout<:
scientific）<}
cout<}
/*cout<<"对矩阵A进行QR分解后所得矩阵:
"<for（i=0;i{
k=0;
cout<<"矩阵A的第"<"<for（j=0;j{
if（fabs（a[i][j]）a[i][j]=0;
cout<:
scientific）<k++;
if（k%3==0）
cout<}
cout<}*/
for（i=0;i{
for（j=0;jfor（k=0,qr[i][j]=0;kqr[i][j]=qr[i][j]+r[i][k]*q[k][j];
}
cout<"<for（i=0;i{
//k=0;
//cout<<"R*Q的第"<"<for（j=0;j<7;j++）
{
if（fabs（qr[i][j]）a[i][j]=0;
cout<:
scientific）<//k++;
//if（k%3==0）
//cout<}
cout<}
cout<"<for（i=0;i{
for（j=7;j{
if（fabs（qr[i][j]）a[i][j]=0;
cout<:
scientific）<}
cout<}
cout<for（i=0;i{
if（L[i][1]==0）
cout<<"λ"<else
cout<<"λ"<}
cout<"<for（k=0;k{
if（L[k][1]==0）//判断特征值是否为实特征值
{
for（i=0;i{
for（j=0;j{
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）;
cout<<"λ"<"<for（j=0;j{
cout<}
}
elsecontinue;
}
}
voidHouseholder（doublea[n][n]）//拟上三角化函数
{
inti,j,k;
intm=0;
doubled,c,h,t;
doubleu[n],p[n],q[n],w[n];
for（i=0;i{
for（j=i+2;j{
if（a[j][i]<=E）
m=m+1;
}
if（m==（n-2-i））
continue;
for（j=i+1,d=0;j{
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;j{
u[j]=a[j][i];
升级会员 