北航数值分析计算实习题目二矩阵QR分解Word文档下载推荐.doc
《北航数值分析计算实习题目二矩阵QR分解Word文档下载推荐.doc》由会员分享,可在线阅读,更多相关《北航数值分析计算实习题目二矩阵QR分解Word文档下载推荐.doc(13页珍藏版)》请在冰豆网上搜索。
本题采用带双步位移的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:
计算
(若,则取)
令
Step4:
(4)QR分解的过程在Xmatrix:
记
对于r=1,2,…,m-1执行
(若)
此算法执行完后,就得到。
3.全部源程序
#include<
stdio.h>
math.h>
#include<
iostream>
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("
voidShowIt_complex()const //此方法用来显示复数特征值
{
%1.12f"
a[i][1]);
if(a[i][2]>
0)printf("
+%1.12fi"
a[i][2]);
if(a[i][2]<
%1.12fi"
}
XMatrixoperator*(XMatrix&
a1)
//重载乘法运算,定义矩阵与矩阵的乘法.若不能相乘,则显示multiplyerror
inti,j,n;
XMatrixtemp(dim_a,a1.dim_b);
if(dim_b!
=a1.dim_a)cout<
<
endl<
"
multiplyerror!
;
for(i=1;
=temp.dim_a;
for(j=1;
=temp.dim_b;
temp.a[i][j]=0;
for(n=1;
n<
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<
adderror!
j++) temp.a[i][j]=a[i][j]+b.a[i][j];
XMatrixoperator*(doublec)
//重载乘法运算,定义矩阵与实数相乘的运算
temp.a[i][j]=c*a[i][j];
}
XMatrixoperator-(XMatrix&
e)
//重载减法运算,定义矩阵与矩阵相减的运算
return(*this)+(e*(-1));
XMatrixT() //此方法用来转置矩阵
XMatrixtemp(dim_b,dim_a);
j++)temp.a[i][j]=a[j][i];
boolcheck_data(inti,intr) //此方法用来判断air是否全都等于0
intj;
doublesum=0;
for(j=i;
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++)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++)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;
=m-1;
if(check_data(r+1,r))continue;
else
{
for(j=r;
=m;
c=(a[r][r]>
h=c*c-c*a[r][r];
r)?
u.a[r][1]=a[r][r]-c;
v=this->
(*this)=(*this)-u*v.T();
p=C.T()*u*(1/h);
q=C*u*(