数值分析第二次大作业

《数值分析》计算实习报告

第二题

院系：机械工程及自动化学院

学号：

姓名：

2017年11月

一、题目要求

试求矩阵A =[a ij ]10×10的全部特征值，并对其中的每一个实特征值求相应的特征向量，已知

a ij ={

sin (0.5i +0.2j ) i ≠j

1.52cos (i +1.2j ) i =j

(i,j =1,2, (10)

说明：

1.用带双步位移的QR 方法求矩阵特征值，要求迭代的精度水平为ε=10−12。

2.打印以下内容： （1）全部源程序；

（2）矩阵A 经过拟上三角化后所得的矩阵A (n−1)； （3）对矩阵A (n−1)实行QR 方法迭代结束后所得的矩阵；

（4）矩阵A 的全部特征值λi =(R i ,I i ) (i =1,2,⋯,10)，其中R i =Re(λi )，I i = Im(λi ) 。若λi 是实数，则令I i =0； （5）A 的相应于实特征值的特征向量。

3.采用e 型数输出实型数，并且至少显示12位有效数字。

二、算法设计思路和方案

1. 将矩阵A 拟上三角化得到矩阵A (n−1)

为了减少计算量，一般先利用Householder 矩阵对矩阵A 作相似变换，把A 化为拟上三角矩阵A (n−1)，然后用QR 方法计算A (n−1)的全部特征值，而A (n−1)的特征值就是A 的特征值。具体算法如下：

记(1)A A =，()r A 的第r 列至第n 列的元素为(r)(1,2,

,;,1,,)ij a i n j r r n ==+。

对于1,2,,2r n =-执行

（1）若()

(2,3,,)r ir

a i r r n =++全为零，则令(1)()r r A A +=，转(5)；否则转（2）。

（2）计算

r d =

()()

1,sgn r r r r r c a d +=- (若()1,0r r r a +=，则取r r c d =)

2()

1,r r r r r r h c c a +=-

（3）令()

()

()()1,2,0,,0,,,,T r r r n r r r

r r r

nr

u a

c a

a

R ++=-∈。

（4）计算

()()(1)()///r T r r r r r r r

T

r r r r

r r r r

r r T T r r r r

p A u h q A u h t p u h q t u A A u u p ωω+====-=--

（5）继续。

当此算法执行完后，就得到与原矩阵A 相似的拟上三角矩阵A (n−1)。 2.用带双步位移的QR 方法求矩阵A 全部特征值

为了加速收敛，采用带双步位移的QR 方法求解A 的全部特征值，具体算法如下：

（1）由步骤1得到拟上三角矩阵(1)n A -；给定精度水平0ε> 和最大迭代次数L 。

（2）记(1)(1)1=n ij n n

A A a -⨯⎡⎤=⎣⎦

，令k = 1,m = n 。

（3）如果()

,1k m m a ε-≤，则得到A 的一个特征值()k mm a ，置m:=m-1（降阶）

，转（4）；否则转（5）。

（4）如果2m ≤，则得到矩阵A 的一个特征值()11k a ，或两个特征值12,s s ，转

（9）；否则转（3）。

（5）如果()1,2k m m a ε--≤，则得到A 的两个特征值12,s s ，置m : = m – 2（降阶），转（4）；否则转（6）。

（6）如果k = L ，则终止计算，未得到A 的全部特征值；否则转（7）。

（7）记()

(1,)k k ij m m A a i j m ⨯⎡⎤=≤≤⎣⎦，计算

()()

1,1()()()()1,1,11,21)

() (k k m m mm

k k k k m m mm m m m m

k k k k k k k T k k k k

s a a t a a a a M A sA tI I m M Q R M QR A Q A Q ------+=+=-=-+==是阶单位矩阵对作分解 （8）置k : = k + 1，转（3）。

（9）A 的全部特征值已计算完毕，停止计算。 其中，k M 的QR 分解与1k A +的计算用下列算法实现：

记(1)()

11[],[],r k ij m m r ij m m k B M b B b C A ⨯⨯====。

对于1,2,,1r m =-执行

（1）若()()1,2,,r ir b i r r m =++全为零，则令11,r r r r B B C C ++==，转（5）

；否则转（2）。

（2）计算

()()()()

1,2()sgn 0,r r r r rr r r r r r r r r r rr

d c b d a c d h c c b +=

=-===-若则取 （3）令()()()()

1,0,,0,,,

,T

r r r m r rr r r r mr u b c b b R +=-∈。

（4）计算

11////T r r r r T r r r r T r r r r r r r r

T

r r r r

r r r r

T T r r r r r r

v B u h B B u v p C u h q C u h t p u h q t u C C u u p ωω++==-====-=--

（5）继续。

此算法执行完后，得到1k m A C +=。 3.用列主元素Gauss 消去法求特征向量

记λ为矩阵A 的实特征值，x 为对应的特征向量，则A x =λx ，即(A - λ)x = 0的解即为矩阵A 的特征向量。因此对于特征向量的求解可采用列主元素Gauss 消去法。

0A I λ-=，经检查由于矩阵A 无重特征值，所以每个特征值对应的特征向量的线性空间为1，矩阵A I λ-的秩为N-1，经过列主元素Gauss 消去法消元后，矩阵有且仅有最后一行全为0。此时设定1n x =-，再依次求出

()1,2,,1k x k n n =--，就可得到对应λ的其中一个特征向量12(,,

,)T n x x x x =，

对应λ的全部特征向量可表示为k x 。最后为了得到标准形式，将向量x 正交化。

具体算法如下。 1. 消元过程 对于1,2,

,1k n =-执行

（1）选行号k i ，使()()

max k k k i k ik k i n

a a ≤≤=。

（2）交换()k kj a 与()(,1,

,)k k i j a j k k n =+所含的数值。

（3）对于1,2,,i k k n =++计算