QR方法求矩阵全部特征值Word文件下载.doc
《QR方法求矩阵全部特征值Word文件下载.doc》由会员分享，可在线阅读，更多相关《QR方法求矩阵全部特征值Word文件下载.doc（16页珍藏版）》请在冰豆网上搜索。
一、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<
stdio.h>
math.h>
#definen3
#defineeps1E-5
doubleDet（doubleL[n][n]）{//用高斯列主元方法求行列式
doubledet=1,t;
for（intk=0;
k<
n-1;
k++）{
doublemax=L[k][k];
intIk=k;
for（intj=k+1;
j<
n;
j++）
if（max<
L[j][k]）{
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<
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;
for（intj=0;
printf（"
%.6f"
H[i][j]）;
intIsZero（doublea[],intj）{//判断数组是否全为0
j;
i++）
if（a[i]!
=0）
return0;
return1;
intsgn（doubley）{//符号函数
if（y>
0）return1;
elseif（y==0）return0;
elsereturn-1;
voidHessenberg（doubleA[n][n]）{//将矩阵化为上Hessenberg矩阵
原矩阵为:
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;
n-2;
for（i=0;
n-k-1;
for（j=0;
if（i==j）I[i][j]=1;
elseI[i][j]=0;
//定义单位阵I[][]
doublemax=fabs（A[k+1][k]）;
//求最大值
fabs（A[i+k+1][k]））
max=fabs（A[i+k+1][k]）;
c[i]=A[i+k+1][k]/max;
//标准化数组
if（IsZero（c,n-k-1））//判断数组是否全为0
continue;
//数组为0，则这一步不需要约化
for（i=0,t=0.0;
t+=c[i]*c[i];
v[k]=sgn（A[k+1][k]）*sqrt（t）;
//
u[0]=c[0]+v[k];
for（j=1;
u[j]=c[j];
//求矩阵u[][]
w=v[k]*（c[0]+v[k]）;
R[i][j]=I[i][j]-u[i]*u[j]/w;
if（i==j）T[i][j]=1;
elseT[i][j]=0;
//定义矩阵T[][]为单位阵
T[i+k+1][j+k+1]=R[i][j];
//初等反射阵T[][]
i++）//矩阵T[][]左乘矩阵A[][]
for（m=0,s=0.0;
m<
m++）
s+=T[i][m]*A[m][j];
B[i][j]=s;
i++）//矩阵T[][]右乘矩阵A[][]
for
升级会员 