文档介绍：数值分析
课
程
设
计
QR方法求矩阵全部特征值
问题复述
用算法求矩阵特征值:
(i) (ii)
要求:
根据算法原理编制求(i)与(ii)中矩阵全部特征值的程序并输出
计算结果(要求误差)
直接用现有的数学软件求(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矩阵结构上的特点,(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化为上三角阵。
数据结构描述
主要数据成员
说明
double A[n][n]
存放矩阵A
double Q[n][n],
存放QR分解式的正交矩阵Q
double R[n][n]
存放QR分解式的上三角阵R
double p[n][n]
Givens矩阵p
double I[n][n]
N阶单位阵
double V[n][n]
存放Q矩阵的转置
double T[n][n]
初等反射阵T
double eps
精度
double max
最大值
double det
存放行列式的值
int count
存放迭代次数
主要函数成员
说明
double Det(double