求实对称三对角矩阵的特征值和特征向量

求实对称三对角矩阵的特征值和特征向量（一）
摘要
在特征值计算问题上，QR方法具有里程碑意义。

QR 方法是一种变换方法，是计算一般矩阵（中小型矩阵）全部特征值问题的最有效方法之一。

QR方法具有收敛快，算法稳定等特点.由于特征值和特征向量能从本质上揭露矩阵的某些重要性质，因而得到它们的精确解十分重要，但其计算一直是很繁琐的数学问题。

特别是当矩阵的阶数较高时，计算量非常大，且不易求其精确解。

关键词：特征值；特征向量；QR分解
Solve Real Symmetry Three Diagonal Matrix Eigenvalue And
Eigenvector
ABSTRACT
Values in the feature, the QR method has milepost sense. QR method is a transformation method, is the calculation of the general matrix ( small and medium-sized matrix ) one of the most effective methods of eigenvalue problems. The QR method has fast convergence, algorithm stability. Because the eigenvalues and eigenvectors can reveal some important properties of matrix from the nature, and thus obtain their exact solutions is very important, but the calculation is very complicated mathematical problems. Especially when the high rank of matrix, the calculation is very large, and is not easy to find the exact solution.
Key words:eigenvalue; eigenvector; QR decomposition
目录
1 绪论 (1)
1.1 问题重述 (1)
1.2研究方法 (1)
2 QR方法 (3)
2.1 QR分解的概念 (3)
2.2 Givens方法 (3)
2.3豪斯霍尔德方法（镜像变换） (5)
2.2.1 Householder 矩阵和Householder变换 (5)
2.2.2QR算法 (6)
3 QR算法C实现过程 (8)
3.1主要参数 (8)
3.2组成模块 (8)
3.3程序改错 (8)
4 测试运行 (11)
参考文献……………………………………………………………………………….…….. 附录…………………………………………………………………………….……………..
1 绪论
1.1 问题重述
（1）用你所熟悉的计算机语言编制利用QR 方法求实对称三对角矩阵全部特征值和特征向量的通用子程序。

（2）利用你所编制的子程序求如下矩阵（从70到80阶）
⎥⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎣⎡=4114114114 A （1.1） 的全部特征值和特征向量。

1.2研究方法
在特征值计算问题上，QR 方法具有里程碑意义。

在1955年的时候，人们还觉得特征值的计算是十分困扰的问题，到1965年它的计算——基于QR 方法的程序已经完全成熟。

直到今天QR 方法仍然是特征值计算的有效方法之一。

设A 是n n ⨯阶矩阵，如果数λ和n 维非零向量x 满足
x Ax λ= （1.2）
则λ称为矩阵A 的一个特征值，x 称为矩阵A 的属于λ的特征向量。

由于特征值和特征向量能从本质上揭露矩阵的某些重要性质，因而得到它们的精确解十分重要，但其计算一直是很繁琐的数学问题。

特别是当矩阵的阶数较高时，计算量非常大，且不易求其精确解。

故在工程技术上，计算矩阵的特征值和特征向量主要使用数值解法，得到其在某一精度水平上的近似解。

常用的算法有：幂法、反幂法、Jacobi 方法和QR 方法。

通过这次课程设计实现用QR 方法求解实对称三对角矩阵的全部特征值和特征向量。

QR 方法是一种变换方法，是计算一般矩阵（中小型矩阵）全部特征值问题的最有效方法之一。

QR 方法具有收敛快，算法稳定等特点.
对矩阵A 进行拟上三角化得到)1(-n A 后，使用带双步位移的QR 方法的迭代公式为：
k k T k k k k k k k k k n Q A Q A QR M R Q M tI
sA A M A A ==+-==+-12)
1(1)
( 分解作对 （1.3）
2 QR 算法应用
2.1 QR 分解的概念
定理1.1[1]
如果实（复）非奇异矩阵A 能化成正交矩阵Q 与实非奇异上三角矩阵R 的乘积， 即
=QR A （2.1） 则称式（2.1）是A 的QR 分解。

引理1.1 任何实的非奇异n 阶矩阵A 可以分解成正交矩阵Q 和上三角矩阵R 的乘积，且除去相差一个对角线元素之绝对值全等于1的对角矩阵因子D 外，分解式（2.1）是唯一的。

定理1.2[1]
设A 为m n ⨯复矩阵（m n ≥），且n 个列向量线性无关，则A 具有分解
R U A = （2.2）
其中U 是m n ⨯复矩阵，且满足H U U ＝Ｉ ，R 是n 阶复非奇异上三角矩阵，且除去相差一个对角元素的模全为1的对角矩阵因子外，分解式（2.2）唯一的。

2.2 Givens 方法
我们知初等旋转变换的性质，即用
ij R 在乘矩阵A 时，仅影响A 的第i 行和第j 行，且选适当的ij R ，就可以消去A 的一个非零元素。

一般地说，作一次旋转可以消去一个
非零元素。

如果在作下一次旋转时不会影响前面已化为零的元素,即不会重新又变成非零,那么,借助于初等旋转矩阵将A 约化成上三角矩阵就有希望。

在平面解析几何中，将向量x 沿顺时针旋转角度θ后变为向量y 时的旋转变换为
cos sin sin cos y x Tx θθθ
θ⎛⎫== ⎪-⎝⎭ （2.3）
由于旋转变换不改变向量的模，即
2222Tx x =， （2.4）
此即,,Tx Tx x x <>=<>，所以T 是正交变换，从而T 是正交矩阵，且
det 1T =-。

（2.5）
一般地说，在n 维欧氏空间n
R 中引入旋转变换如下。

定义1.2[2] 设实数c 与s 满足221c s += ，称 111
11ij c s T s c ⎛⎫ ⎪ ⎪ ⎪ ⎪
⎪ ⎪= ⎪- ⎪ ⎪ ⎪ ⎪ ⎪⎝⎭()i j < （2.4）
为Givens 矩阵（或初等旋转矩阵（Elementary Rotation Matrix ）），有时也记为(,)ij ij T T c s =，由Givens 矩阵所确定的线性变换称为Givens 变换（或初等旋转变换（Elementary Rotation Transformation ））.
当221c s +=时，必有角度θ，使得cos ,sin c s θθ==
性质1.2.1[2] Givens 矩阵是正交矩阵且有
1[(,)][(,)](,)
T ij ij ij T c s T c s T c s -==- （2.5） det[(,)]1ij T c s = （2.6）
性质1.2.2[2] 设1212(,,...,),(,,...,)T T
n ij n x y T x εεεηηη=== 则有
(,)
i i j j j k k c s s c k i j ηεεηεεηε=+⎧⎪=-+⎨⎪=≠⎩ （2.7）
式（2.7）表明，当22
0i j εε+≠ 时，选取
第i 行 第j 行
2
22
2i
i j i j c s εεεεε==++ （2.8）
就可使
22
0,0
i i j j ηεεη=+>=。

（2.9）
定理1.3[2] 设12(,,...,)0T n x εεε=≠，则存在有限多个Givens 矩阵的乘积,记为T ,
使得
12Tx x e =，其中1(1,0,...,0)T e 。

（2.10） 推论1[2].设非零向量n x R ∈及单位列向量n z R ∈ ，则存在有限多个Givens 矩阵的乘积，记为T ，使得 2Tx x z =。

（2.11）
2.3豪斯霍尔德方法 （镜像变换）
2.3.1 Householder 矩阵 和Householder 变换
定义3.1[2] 设n u R ∈，且21u = 称
2T I uu H =- （2.12）
为Householder 矩阵（或初等反射矩阵（Elementary Reflection Matrix ））,由Householder 矩阵所确定的初等变换称为Householder 变换（或初等反射变换
（（Elementary Reflection Transformation ））.
Householder 矩阵具有如下性质：
（1）T H =H （对称矩阵）；
（2）T I H H = （正交矩阵）；
（3）2I H = （对合矩阵）；
（4）1-H =H （自逆矩阵）；
（5）det 1H =-
使用式（2.12）可直接验证性质(1)到性质(4)
设3u R ∈ 为一个单位列向量，s 是过原点且与u 垂直的平面，312,R νννν∀∈=+,其中
12,s s νν∈⊥ ，则11112T uu ννννH =-=,(因为u s ⊥ ，而1s ν∈ ，故1u ν⊥ ，从而1
0T u ν=) 2222222222()2()2T T T uu uu ku ku u u ννννννννH =-=-=-=-=-
（因为u s ⊥，故u s ⊥∈ ，而2s ν⊥∈，故2ν与u 平行，从而2ku ν=）。

这样ν经变换后的
镜像νH 是ν关于s 对称向量1212νννννH =H +H =- 。

(6)n
x R ∀∈ ，记y x =H ，则有 2()2()T T y x x uu x x u x u =H =-=- （2.13） T T T T y y x x x x =H H = （2.14）
式（2.13）第二个等号成立的理由是矩阵乘法满足结合律，而T u x 为实数。

式（2.14）表明
222y x x =H = 。

定理1.4[2] 设n x R ∈为非零列向量 ，n z R ∈ 为单位列向量，则存在Householder 矩
阵H ，使得2x x z
H = 。

2.3.2 QR 算法
假设A 为非奇异矩阵，则由定理1.1知A 有QR 分解，先说
QR A =
令
1A =A ，对1A 作QR 分解，得
交换1Q 与1R 的次序 ，得
211Q R A =
再将2A 作QR 分解 ，得
222Q R A =
又交换2Q 与2R 的次序 ，得
111Q R A =
322R Q A =
如此反复作下去，得迭代序列{}k A 如下
11Q R (1,2,...)
Q k k k k k k k R +A =A ⎧⎪A ==⎨⎪A =⎩ （2.15）
称式（2.15）为矩阵A 的QR 算法（Algorithm ） 由式（2.15）易知 ，1Q k k k R -=A ，故
11Q Q k k k k -+A =A ，这表明矩阵序列{}k A 中任何两个相邻矩阵都是相似的，从而每个k A 都与矩阵A 相似。

3 QR解法C实现过程
3.1 主要参数
求实对称三对角对称矩阵的全部特征值及特征向量可利用变型QR方法计算，下面是主要参数：
n-矩阵的阶数；
b-长度为n的数组，返回时存放三对角阵的主对角线元素；
c-长度为n的数组，返回时前n-1个元素存放次对角线元素；
q-长度为n*n的数组，若存放单位矩阵，则返回实对称三对角矩阵的特征向量组；
若存放Householder变换矩阵，则返回实对称矩阵A的特征向量组；
a-长度为n*n的数组，存放n阶实对称矩阵；
3.2 组成模块
程序由两个函数构成，分别是主函数和功能函数[3]：
int main(){}
主要负责各变量的初始化；
int sstq(int n,double b[],double c[],double q[],double eps,int l){}
主要负责通过多个循环实现求出实对称三对角矩阵的全部特征值和特征向量；
3.3 程序改错
错误一如图3.1所示：
图3.1 建立数组出错
改正： static double b[N],c[N],q[N]这一语句出错，将其改为static double b[N],c[N],q[N][N]后不报错。

错误二如图3.2所示：
图3.2 数组A初始化出错改正：
for(i=0;i<N;i++)
for(j=0;j<N;j++)
{
if(j==i)
{
a[i][j]=4;
q[i][j]=1;
}
if(j==i-1||j==i+1)
a[i][j]=1;
else a[i][j]=0;
}
改为：
for(i=0;i<N;i++)
for(j=0;j<N;j++)
{
if(j==i)
{
a[i][j]=4;
q[i][j]=1;
}
if(j==i-1||j==i+1)
a[i][j]=1;
}后不报错，错在赋值出错。

4 测试运行
首先，我用三阶实对称三对角矩阵测试运行
4.000000 1.0000000.0000001.000000 4.000000 1.0000000.000000
1.000000
4.000000⎡⎤
⎢⎥⎢⎥⎢⎥⎣⎦
运行结果：
图4.1 三阶实对称三对角矩阵运行结果
显示结果正确，然后按照老师的要求用70阶实对称三对角矩阵测试运行：
图4.2 70阶矩阵的测试运行结果截图(局部)
由于dos界面的显示能力有限，不能显示70阶矩阵的全部结果，故通过cmd运行exe 文件并将其结果导出至一个文本文件中便可分析其对错。

结果如下图所示：
图4.3 导出结果至文本（1）
图4.4 导出结果至文本（2）
参考文献
[1] 杨大地,王开荣.数值分析[M].北京:科学出版社,2006.
[2] 颜庆津.数值分析[M].北京:北京航空航天大学出版社,2006.
[3] [美] Ivor Horton著,李颂华、康会光译.Visual C++ 2005入门经典[M].北京:清华出版社,2007.
附录源程序清单
#include "stdio.h"
#include "math.h"
#define N 70
int sstq(int n,double b[],double c[],double q[],double eps,int l) {
int i,j,k,m,it,u,v;
double d,f,h,g,p,r,e,s;
c[n-1]=0.0;
d=0.0;
f=0.0;
for (j=0; j<=n-1; j++)
{
it=0;
h=eps*(fabs(b[j])+fabs(c[j]));
if (h>d)
d=h;
m=j;
while ((m<=n-1)&&(fabs(c[m])>d))
m=m+1;
if (m!=j)
{
do
{
if (it==l)
printf("fail\n");
return(-1);
}
it=it+1;
g=b[j];
p=(b[j+1]-g)/(2.0*c[j]);
r=sqrt(p*p+1.0);
if (p>=0.0)
b[j]=c[j]/(p+r);
else
b[j]=c[j]/(p-r);
h=g-b[j];
for (i=j+1; i<=n-1; i++)
b[i]=b[i]-h;
f=f+h;
p=b[m];
e=1.0;
s=0.0;
for (i=m-1; i>=j; i--)
{
g=e*c[i];
h=e*p;
if (fabs(p)>=fabs(c[i]))
{
e=c[i]/p;
r=sqrt(e*e+1.0);
c[i+1]=s*p*r;
s=e/r;
e=1.0/r;
else
{
e=p/c[i];
r=sqrt(e*e+1.0);
c[i+1]=s*c[i]*r;
s=1.0/r;
e=e/r;
}
p=e*b[i]-s*g;
b[i+1]=h+s*(e*g+s*b[i]);
for (k=0; k<=n-1; k++)
{
u=k*n+i+1;
v=u-1;
h=q[u];
q[u]=s*q[v]+e*h;
q[v]=e*q[v]-s*h;
}
}
c[j]=s*p;
b[j]=e*p;
}
while (fabs(c[j])>d);
}
b[j]=b[j]+f;
}
for (i=0; i<=n-1; i++)
{
k=i; p=b[i];
{
j=i+1;
while ((j<=n-1)&&(b[j]<=p))
{
k=j;
p=b[j];
j=j+1;
}
}
if (k!=i)
{
b[k]=b[i];
b[i]=p;
for (j=0; j<=n-1; j++)
{
u=j*n+i;
v=j*n+k;
p=q[u];
q[u]=q[v];
q[v]=p;
}
}
}
return(l);
}
main()
{
int i,j,k,l=60;
double eps = 0.000001;
static double b[N],c[N],q[N][N];
static double a[N][N];
for(i=0;i<N;i++)
for(j=0;j<N;j++)
{
if(j==i)
{
a[i][j]=4;
q[i][j]=1;
}
if(j==i-1||j==i+1)
a[i][j]=1;
}
for(i=0;i<N;i++)
{
b[i] = 4;
c[i] = 1;
}
k=sstq(N,b,c,q,eps,l);
k = 1;
printf("您输入的矩阵实对称三对角矩阵A:\n");
for(i=0;i<N;i++)
{
for(j=0;j<N;j++)
printf("%6lf\t",a[i][j]);
printf("\n");
}
printf("\n");
if(k>0)
{
求实对称三对角矩阵的特征值和特征向量（一）printf("全部特征值:\n");
for(i=0;i<N;i++)
printf("%6lf\t",b[i]);
printf("\n\n");
printf("对应每个特征值的特征向量(按列看):\n");
for(i=0;i<N;i++)
{
for(j=0;j<N;j++)
printf("%6lf\t",q[i][j]);
printf("\n");
}
printf("\n");
}
}。