北航数值分析计算实习题目二矩阵QR分解.doc
- 文档编号:118378
- 上传时间:2022-10-03
- 格式:DOC
- 页数:13
- 大小:562.50KB
北航数值分析计算实习题目二矩阵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]=(j 0: 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*(
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航 数值 分析 计算 实习 题目 矩阵 QR 分解
![提示](https://static.bdocx.com/images/bang_tan.gif)