北航数值分析计算实习题目二矩阵QR分解.doc
- 文档编号:143134
- 上传时间:2023-04-28
- 格式: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*(1/h); t=p.T()*u*(1/h); w=q-u*(t.getdata()); C=C-w*u.T()-u*p.T(); } } return0; } intdoubleQR(XMatrix&lambda) //此方法为带双步位移的QR分解,其中调用QR() { intm,i,k,n1,n2; doubles1,t1,delta; i=1;k=0; m=dim_a; step3: if(abs(a[m][m-1])<=1e-10) { lambda.a[i][1]=a[m][m];lambda.a[i][2]=0;i++; m=m-1;gotostep4; } elsegotostep5; step4: if(m==1){lambda.a[i][1]=a[1][1];lambda.a[i][2]=0;i++;return0;} if(m==0)return0; if(m>1)gotostep3; step5: if(m==2) { delta=(a[m-1][m-1]+a[m][m])*(a[m-1][m-1]+a[m][m])-4*(a[m-1][m-1]*a[m][m]-a[m][m-1]*a[m-1][m]); if(delta<0) { lambda.a[i][1]=(a[m-1][m-1]+a[m][m])/2;lambda.a[i][2]=sqrt(-delta)/2; i++; lambda.a[i][1]=(a[m-1][m-1]+a[m][m])/2;lambda.a[i][2]=-(sqrt(-delta))/2;i++; } else { lambda.a[i][1]=(a[m-1][m-1]+a[m][m]+sqrt(delta))/2; lambda.a[i][2]=0;i++; lambda.a[i][1]=(a[m-1][m-1]+a[m][m]-sqrt(delta))/2; lambda.a[i][2]=0;i++; } return0; } elsegotostep7; step7: if(abs(a[m-1][m-2])<=1e-10) { delta=(a[m-1][m-1]+a[m][m])*(a[m-1][m-1]+a[m][m])-4*(a[m-1][m-1]*a[m][m]-a[m][m-1]*a[m-1][m]); if(delta<0) { lambda.a[i][1]=(a[m-1][m-1]+a[m][m])/2;lambda.a[i][2]=sqrt(-delta)/2;i++; lambda.a[i][1]=(a[m-1][m-1]+a[m][m])/2;lambda.a[i][2]=-(sqrt(-delta))/2;i++; } else { lambda.a[i][1]=(a[m-1][m-1]+a[m][m]+sqrt(delta))/2; lambda.a[i][2]=0;i++; lambda.a[i][1]=(a[m-1][m-1]+a[m][m]-sqrt(delta))/2; lambda.a[i][2]=0;i++; } m=m-2;gotostep4; } elsegotostep8; step8: if(k>=10000){cout<<"error"< step9: XMatrixM(m,m),I(m,m),C(dim_a,dim_b); C=(*this);C.dim_a=m;C.dim_b=m; for(n1=1;n1<=m;n1++) for(n2=1;n2<=m;n2++)I.a[n1][n2]=(n1==n2)? 1: 0;//initializeI s1=a[m-1][m-1]+a[m][m]; t1=a[m-1][m-1]*a[m][m]-a[m][m-1]*a[m-1][m]; M=C*C-C*s1+I*t1; M.QR(m,C); for(n1=1;n1<=C.dim_a;n1++) for(n2=1;n2<=C.dim_b;n2++) a[n1][n2]=C.a[n1][n2]; k++; gotostep3; } intGet_max_line_according_column(XMatrix&A,intk) //此方法找出按列主元素 { intj,temp=k; for(j=k+1;j<=A.dim_a;j++) if(abs(A.a[j][k])>abs(A.a[j-1][k]))temp=j; returntemp; } intExchange_Line(XMatrix&A,inta,intb,intk) //此方法用来交换矩阵的两行元素 { doubletemp; inti; for(i=k;i<=dim_b;i++) {temp=A.a[a][i]; A.a[a][i]=A.a[b][i]; A.a[b][i]=temp;} return0; } intGet_Eigenvector(doublelambda,XMatrix&ev) //此方法用来求矩阵的特征向量 { XMatrixA(dim_a,dim_b),I(dim_a,dim_b); intn1,n2,k,ik,i,j; doublem,sum; for(n1=1;n1<=dim_a;n1++) for(n2=1;n2<=dim_b;n2++)I.a[n1][n2]=(n1==n2)? 1: 0;//initializeI A=(*this)-I*lambda; for(n1=1;n1<=dim_a;n1++)A.a[n1][0]=0; for(k=1;k<=dim_a-1;k++) { ik=Get_max_line_according_column(A,k); if(k! =ik)Exchange_Line(A,k,ik,k); for(i=k+1;i<=dim_a;i++) { m=A.a[i][k]/A.a[k][k]; for(j=k;j<=dim_b;j++) A.a[i][j]=A.a[i][j]-m*A.a[k][j]; A.a[i][0]=A.a[i][0]-m*A.a[k][0]; } } ev.a[dim_a][1]=1; for(k=dim_a-1;k>=1;k--) { sum=0; for(j=k+1;j<=dim_b;j++) sum=sum+A.a[k][j]*ev.a[j][1]; ev.a[k][1]=(A.a[k][0]-sum)/A.a[k][k]; } return0; } }; voidmain() //主程序 { inti,j; XMatrixA(10,10),B(10,10),lambda(10,2),ev(10,1); //A(10x10)用来存储矩阵进行QR分解,B(10x10)用求特征向量,lambda(10x2),用来存放特征 //值,ev(10,1)用来存放特征向量. for(i=1;i<=10;i++) for(j=1;j<=10;j++) A.a[i][j]=(i==j)? (1.5*cos(i+1.2*j)): (sin(0.5*i+0.2*j)); B=A; A.preQR(); //拟上三角化 cout<<"对A进行拟上三角化结果"< A.ShowIt(); //输出拟上三角化结果 A.doubleQR(lambda); //带双步位移的QR分解 cout<<"对A进行QR分解结果"< A.ShowIt(); //输出QR分解结果 cout<<"A的全部特征值"< lambda.ShowIt_complex(); //输出A的全部特征值 for(i=1;i<=lambda.dim_a;i++) { if(lambda.a[i][2]==0) { B.Get_Eigenvector(lambda.a[i][1],ev); //求特征向量 cout<<"λ="< ev.ShowIt(); //输出特征向量 } } } 4.矩阵A经过拟上三角化后所得的矩阵A(n-1) 只保留3位有效数字的结果: 5.对A(n-1)进行QR分解方法结束后所得的矩阵 只保留3位有效数字的结果 6.矩阵A的全部特征值 7.A的相应于实特征值的特征向量 13
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航 数值 分析 计算 实习 题目 矩阵 QR 分解
![提示](https://static.bingdoc.com/images/bang_tan.gif)