北航数值分析计算实习题目二矩阵QR分解.doc
《北航数值分析计算实习题目二矩阵QR分解.doc》由会员分享，可在线阅读，更多相关《北航数值分析计算实习题目二矩阵QR分解.doc（13页珍藏版）》请在冰豆网上搜索。
数值分析
实习二
院（系）名称
航空科学与工程学院
专业名称
动力工程及工程热物理
学号
SY0905303
学生姓名
解立垚
1.题目
试用带双步位移QR的分解法求矩阵A=[aij]10*10的全部特征值，并对其中的每一个实特征值求相应的特征向量。
已知。
说明：
1、求矩阵特征值时，要求迭代的精度水平为。
2、打印以下内容：
算法的设计方案；
全部源程序（要求注明主程序和每个子程序的功能）；
矩阵A经过拟上三角话之后所得的矩阵；
对矩阵进行QR分解方法结束后所得的矩阵；
矩阵A的全部特征值，和A的相应于实特征值的特征向量；
其中如果是实数，则令
3、采用e型输出数据，并且至少显示12位有效数字。
2.算法设计方案
本题采用带双步位移的QR分解方法。
为了使程序简洁，自定义类Xmatrix，其中封装了所需要的函数方法。
在Xmatrix类中封装了运算符重载的函数，即定义了矩阵的加、减、乘、除、数乘运算及转置运算（T（））。
同时为了避免传递数组带来的额外内存开销，使用引用（&）代替值传递，以节省内存空间，避免溢出.
（1）此程序的主要部分为Xmatrix中的doubleQR（）方法，具体如下：
Step1：
使用矩阵拟上三角化的算法将A化为拟上三角阵A（n-1）（此处调用Xmatrix中的preQR（）方法）
Step2：
令，其中k为迭代次数。
Step3：
如果，则得到A的一个特征值，令，gotoStep4;否则gotoStep5.
Step4：
如果，则得到A的一个特征值，gotoStep11;如果，则gotoStep11;如果，则gotoStep3;
Step5（Step6）：
如果，则得到A的两个特征值（为右下角两阶子阵对应的特征方程的两个根。
），gotoStep11，否则gotoStep7：
Step7：
如果，则得到A的两个特征值（定义同上），令，gotoStep4；否则gotoStep8.
Step8：
如果，则输出“error”否则gotoStep9.
Step9（Step10）：
对A（m×m阶）进行QR分解，此处调用Xmatrix中的方法QR（）.令k=k+1，gotoStep3.
Step11：
若到此步，则认为程序已经计算出全部特征值或已超过最大迭代次数k，退出此函数，相当于return0.
（2）在求特征向量对应的特征值时，采用列主元素的高斯消元法，来解方程组（A-λI）x=0，以得出某一组特征向量。
在消元法的过程中，首先进行行交换，把aik（i=k,k+1,…,n）中绝对值最大的元素交换到第k行的主对角线上，然后再使用第二种初等行变换进行消元。
由于A-λI为非满秩矩阵，故这样变换得到的矩阵ann=0，令xn=1,即可解出某一个特征向量。
具体算法在Xmatrix:
:
Get_Eigenvector中，如下：
消元过程:
对于k=1,2,…,n-1执行
Step1:
选行号ik,使
Step2:
交换以及所含的数值
Step3:
对于计算
回代过程:
（3）拟上三角化的过程在Xmatrix:
:
QR（）中，具体算法如下：
记，并记的第r列至第n列的元素为
对执行：
Step1:
若全为零，则令，转（5）；否则转
（2）。
Step2:
计算
（若，则取）
Step3:
令
Step4:
计算
（4）QR分解的过程在Xmatrix:
:
QR（）中，具体算法如下：
记
对于r=1，2，…，m-1执行
Step1:
若全为零，则令，转（5）；否则转
（2）。
Step2:
计算
（若）
Step3:
令
Step4:
计算
此算法执行完后，就得到。
3.全部源程序
#include
#include
#include
usingnamespacestd;
classXMatrix //定义类Xmatrix
{
private:
public:
intdim_a,dim_b; //类成员dim_a和dim_b为矩阵的行数和列数
doublea[11][11]; //类成员a数组用来存储矩阵元素
XMatrix（inti=1,intj=1）:
dim_a（i）,dim_b（j） //类的构造函数，传递行数和列数
{
}
voidShowIt（）const //此方法用来显示矩阵元素
{
for（inti=1;i<=dim_a;i++）
for（intj=1;j<=dim_b;j++）
{
printf（"%1.12f",a[i][j]）;
if（j==dim_b）printf（"\n"）;
}
printf（"\n"）;
}
voidShowIt_complex（）const //此方法用来显示复数特征值
{
for（inti=1;i<=dim_a;i++）
{
printf（"%1.12f",a[i][1]）;
if（a[i][2]>0）printf（"+%1.12fi",a[i][2]）;
if（a[i][2]<0）printf（"%1.12fi",a[i][2]）;
printf（"\n"）;
}
printf（"\n"）;
}
XMatrixoperator*（XMatrix&a1）
//重载乘法运算,定义矩阵与矩阵的乘法.若不能相乘,则显示multiplyerror
{
inti,j,n;
XMatrixtemp（dim_a,a1.dim_b）;
if（dim_b!
=a1.dim_a）cout<";
for（i=1;i<=temp.dim_a;i++）
for（j=1;j<=temp.dim_b;j++）
{
temp.a[i][j]=0;
for（n=1;n<=dim_b;n++）temp.a[i][j]=temp.a[i][j]+a[i][n]*a1.a[n][j];
}
returntemp;
}
XMatrixoperator+（XMatrix&b）
//重载加法运算,定义矩阵与矩阵的加法运算,若不能相加,则显示adderror.
{
inti,j;
XMatrixtemp（dim_a,dim_b）;
if（dim_a!
=b.dim_a||dim_b!
=b.dim_b）cout<";
for（i=1;i<=temp.dim_a;i++）
for（j=1;j<=temp.dim_b;j++） temp.a[i][j]=a[i][j]+b.a[i][j];
returntemp;
}
XMatrixoperator*（doublec）
//重载乘法运算,定义矩阵与实数相乘的运算
{
inti,j;
XMatrixtemp（dim_a,dim_b）;
for（i=1;i<=dim_a;i++）
for（j=1;j<=dim_b;j++）
temp.a[i][j]=c*a[i][j];
returntemp;
}
XMatrixoperator-（XMatrix&e）
//重载减法运算,定义矩阵与矩阵相减的运算
{
return（*this）+（e*（-1））;
}
XMatrixT（） //此方法用来转置矩阵
{
inti,j;
XMatrixtemp（dim_b,dim_a）;
for（i=1;i<=temp.dim_a;i++）
for（j=1;j<=temp.dim_b;j++）temp.a[i][j]=a[j][i];
returntemp;
}
boolcheck_data（inti,intr） //此方法用来判断air是否全都等于0
{
intj;
doublesum=0;
for（j=i;j<=dim_a;j++）
if（abs（a[j][r]）>1e-12）returnfalse;
returntrue;
}
doublegetdata（） //此方法用来获得一阶方阵的唯一元素的值
{
if（dim_a==1&&dim_b==1）returna[1][1];else{cout<<"geterror!
";return0;};
}
intpreQR（） //此方法用来拟上三角化矩阵
{
intr,j;
XMatrixu（dim_a,1）,p（dim_a,1）,q（dim_a,1）,w（dim_a,1）,t（1,1）;
doubled,c,h;
for（r=1;r<=dim_a-2;r++）
{
if（check_data（r+2,r））continue;
else
{
d=0;
for（j=r+1;j<=dim_a;j++）d=sqrt（d*d+a[j][r]*a[j][r]）;
c=（a[r+1][r]>0）?
（-d）:
d;h=c*c-c*a[r+1][r];
for（j=1;j<=dim_a;j++）u.a[j][1]=（j<=r）?
0:
a[j][r];
u.a[r+1][1]=a[r+1][r]-c;
p=this->T（）*u*（1/h）;
q=（*this）*u*（1/h）;
t=p.T（）*u*（1/h）;
w=q-u*t.getdata（）;
（*this）=（*this）-（w*u.T（））-（u*p.T（））;
}
}
return0;
}
intQR（intm,XMatrix&C） //此方法对A进行QR分解
{
intr,j;
doubled,c,h;
XMatrixu（m,1）,v（m,1）,p（m,1）,q（m,1）,w（m,1）,t（1,1）;
for（r=1;r<=m-1;r++）
{
if（check_data（r+1,r））continue;
else
{
d=0;
for（j=r;j<=m;j++）d=sqrt（d*d+a[j][r]*a[j][r]）;
c=（a[r][r]>0）?
（-d）:
d;h=c*c-c*a[r][r];
for（j=1;j<=m;j++）u.a[j][1]=（j0:
a[j][r];
u.a[r][1]=a[r][r]-c;
v=this->T（）*u*（1/h）;
（*this）=（*this）-u*v.T（）;
p=C.T（）*u*（1/h）;
q=C*u*（
升级会员 