QR方法求矩阵全部特征值.doc
- 文档编号:134283
- 上传时间:2022-10-04
- 格式:DOC
- 页数:16
- 大小:423KB
QR方法求矩阵全部特征值.doc
《QR方法求矩阵全部特征值.doc》由会员分享,可在线阅读,更多相关《QR方法求矩阵全部特征值.doc(16页珍藏版)》请在冰豆网上搜索。
数值分析
课
程
设
计
QR方法求矩阵全部特征值
问题复述
用算法求矩阵特征值:
(i)(ii)
要求:
(1)根据算法原理编制求(i)与(ii)中矩阵全部特征值的程序并输出
计算结果(要求误差)
(2)直接用现有的数学软件求(i),(ii)的全部特征值,并与
(1)的结果比较。
问题分析
QR方法是目前公认为计算矩阵全部特征值的最有效的方法。
它适用于任一种实矩。
QR方法的原理是利用矩阵的正交分解产生一个与矩阵A相似的矩阵迭代序列,这个序列将收敛于一个上三角阵或拟上三角阵,从而求得原矩阵A的全部特征值。
QR是一个迭代算法,每一步迭代都要进行QR分解,再作逆序的矩阵乘法。
要是对矩阵A直接用QR方法,计算量太大,效率不高。
因此,一个实际的QR方法计算过程包括两个步骤,首先要对矩阵A用初等相似变换约化为上Hessenberg矩阵,再用QR方法求上Hessenberg矩阵的全部特征值。
示意如下:
对B矩阵的约化只需将每列次对角线上的元素约化为0。
因此常常用平面旋转阵(Givens变换阵)来进行约化。
一、QR方法原理及收敛性
设已实现了QR分解,记
其中是正交阵,是上三角阵。
因为,用对作正交相似变换有
可改写为
显然只是的QR分解因子阵的逆序相乘,而且与原矩阵有相同的特征值。
对矩阵再重复以上过程并继续下去,可以得到一个与原矩阵A有相同特征值的矩阵序列。
其过程如下:
记
对作正交分解
作矩阵,~
对作正交分解
作矩阵,~~,~
重复以上过程可得一般的形式为
对作正交分解
构成矩阵序列(k=1,2……)
~A
从矩阵A开始得到一个矩阵序列
这个矩阵序列中每一个矩阵都与原矩阵相似,即都有与A相同的特征值。
这个矩阵序列在实质上收敛于依次以为对角元的上三角阵。
具体可以表示为
其中
的极限不一定存在。
二、用正交相似变换约化矩阵为上Hessenberg阵
用Householder变换可以将一个向量指定的某个分量以下的各分量变为0。
我们只要求消掉A的次对角线以下的元素,即将A约化为上Hessenberg阵。
为了使变换前后矩阵的特征值不变,需要用Householder矩阵对A作相似变换,即用正交阵同时左乘和右乘A时,原来已变为0的元素不再改变。
若设是Householder矩阵,用它对A的第一列元素的变换示意如下:
依次对A的各列进行类似的变换,一共要进行次变换,最终可以得到一个与原矩阵A有相同特征值的上Hessenberg阵。
以上约化A为上Hessenberg阵的过程可以用一系列Householder矩阵来实现。
其中,对于每一个有
经过步约化就可得到一个上Hessenberg阵(A的第列不需要约化)
三、Hessenberg阵的QR算法
设矩阵,其特征值都是实数。
若已用Householder变换约化为上Hessenberg阵
对已得到的上Hessenberg阵可用QR变换,经过迭代过程约化为上三角形矩阵以求出A的特征值。
只要A的特征值是实数,将B约化为上三角形矩阵总是可以做到的。
由于B矩阵结构上的特点,对B矩阵的约化只需将每列次对角线上的元素约化为0.可用平面旋转阵(Givens变换阵)来进行约化。
n阶方阵为平面旋转阵
。
还是一个非对称的正交阵,有,也是一个平面旋转阵。
有以下几种作用:
1、左乘向量只改变x的第i个和第j个分量。
现构造对x作变换使的结果将x的第j个分量约化为0。
令,有
调整角可使。
若记,按下式选取
于是
有。
2、对非零的n维向量x连续左乘,可将x的第i+1到第n个分量都约化为零;即
其中
3、用对矩阵A作变换得到的结论是
左乘A只改变A的第i,j行;
右乘A只改变A的第i,j列;
用对A作正交相似变换改变了A的i行和j行以及i列和j列。
用一系列连续左乘矩阵A,可以将矩阵A化为上三角阵。
数据结构描述
主要数据成员
说明
doubleA[n][n]
存放矩阵A
doubleQ[n][n],
存放QR分解式的正交矩阵Q
doubleR[n][n]
存放QR分解式的上三角阵R
doublep[n][n]
Givens矩阵p
doubleI[n][n]
N阶单位阵
doubleV[n][n]
存放Q矩阵的转置
doubleT[n][n]
初等反射阵T
doubleeps
精度
doublemax
最大值
doubledet
存放行列式的值
intcount
存放迭代次数
主要函数成员
说明
doubleDet(doubleL[n][n])
用高斯列主元方法求行列式
intNon_singularMatrix(doubleL[n][n])
判断是否是非奇异矩阵
voidDisp(doubleH[n][n])
输出矩阵
intIsZero(doublea[],intj)
判断数组是否全为0
intsgn(doubley)
符号函数
voidHessenberg(doubleA[n][n])
将矩阵化为上Hessenberg矩阵
intIsHessenberg(doubleE[n][n])
判断是否是上Hessenberg矩阵
voidQRAlgorithm(doubleA[n][n])
QR算法求特征值
voidSeekEigenvalue(doubleA[n][n])
判断是否满足QR算法条件,满足则进行QR方法求特征值
算法的描述(流程图)
源程序
C程序为:
#include
#include
#definen3
#defineeps1E-5
doubleDet(doubleL[n][n]){//用高斯列主元方法求行列式
doubledet=1,t;
for(intk=0;k doublemax=L[k][k];intIk=k; for(intj=k+1;j if(max max=L[j][k]; Ik=j; } if(max==0)return0; if(Ik! =k){ for(intj=k;j t=L[Ik][j]; L[Ik][j]=L[k][j]; L[k][j]=t; } det*=-1; } for(inti=k+1;i t=L[i][k]/L[k][k]; L[i][k]=t; for(intj=k+1;j L[i][j]=L[i][j]-t*L[k][j]; } det*=L[k][k]; } if(L[n-1][n-1]==0)return0; else{ printf("矩阵A[][]的行列式为: %.5f\n",det*L[n-1][n-1]); returndet*L[n-1][n-1]; } } intNon_singularMatrix(doubleL[n][n]){//判断是否是非奇异矩阵 printf("判断矩阵A[][]的行列式是否为0? \n"); if(Det(L)! =0)return1; return0; } voidDisp(doubleH[n][n]){//输出矩阵 for(inti=0;i for(intj=0;j printf("%.6f",H[i][j]); printf("\n"); } } intIsZero(doublea[],intj){//判断数组是否全为0 for(inti=0;i if(a[i]! =0) return0; return1; } intsgn(doubley){//符号函数 if(y>0)return1; elseif(y==0)return0; elsereturn-1; } voidHessenberg(doubleA[n][n]){//将矩阵化为上Hessenberg矩阵 printf("原矩阵为: \n");Disp(A);//输出原矩阵 doubleT[n][n],B[n][n],C[n][n],c[n-1],v[n-1],u[n-1],R[n-1][n-1],I[n-1][n-1],t,w,s; inti,j,k,m; for(k=0;k for(i=0;i for(j=0;j if(i==j)I[i][j]=1; elseI[i][j]=0;//定义单位阵I[][] doublemax=fabs(A[k+1][k]);//求最大值 for(i=0;i if(max max=fabs(A[i+k+1][k]); } for(i=0;i c[i]=A[i+k+1][k]/max;//标准化数组 if(IsZero(c,n-k-1))//判断数组是否全为0 continue;//数组为0,则这一步不需要约化 for(i=0,t=0.0;i t+=c[i]*c[i]; v[k]=sgn(A[k+1][k])*sqrt(t);// u[0]=c[0]+v[k]; for(j=1;j u[j]=c[j];//求矩阵u[][] w=v[k]*(c[0]+v[k]); for(i=0;i for(j=0;j R[i][j]=I[i][j]-u[i]*u[j]/w; for(i=0;i for(j=0;j if(i==j)T[i][j]=1; elseT[i][j]=0;//定义矩阵T[][]为单位阵 } for(i=0;i for(j=0;j T[i+k+1][j+k+1]=R[i][j];//初等反射阵T[][] for(i=0;i for(j=0;j for(m=0,s=0.0;m s+=T[i][m]*A[m][j]; B[i][j]=s; } for(i=0;i for
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- QR 方法 矩阵 全部 特征值