0.499988 -0.5 0.707115
此结果说明迭代6次，求得误差为1.16685?10-9的矩阵A的所有特征值?1=0.585786，?2=3.41421，?3=2 及对应的特征向量
{0.500012，0.707107，0.499988}，
{-0.5，0.707107，-0.5}，
{-0.707098，0.0000120684，0.707115}
本题矩阵A的3个特征值为{2-sqrt[2],2+sqrt[2], 2}。线性代数告诉我们，不同特征值对应的特征向量是正交的，用如上求出的特征向量验证，可以得到
-17
{0.500012，0.707107，0.499988}.{-0.5，0.707107，-0.5}= -5.55112?10，
{0.500012，0.707107，0.499988}.{-0.707098，0.0000120684，0.707115}=-5.55112?10
{-0.5，0.707107，-0.5}.{-0.707098，0.0000120684，0.707115}=-1.66533?10
可见求出的特征向量有很好的正交性，这也是Jacobi方法的优点
-16
-17
QR 方 法
求镜面反射矩阵算法
1．输入向量x，维数n 2. 计算a ?sgn(x1)||x||2 3. x1 ?x1 +a 4. b ? a.x1
5. 计算H?E-b-1xTx
注：向量x=( x1 ,x2 ,… ,xn)T,E为n阶单位矩阵。该算法算出的镜面反射矩阵H左乘给定的非零向量x=( x1，x2,???? , xn)T后，获得只有第一个分量不为零的向量 Hx=-?e 1
T
其中e1=( 1，0,???? , 0), ? =?||x|| 2。
求Householder矩阵（镜面反射矩阵）程序
x=Input[“输入非零向量x=”]
n=Length[x];
a=Sign[x[[1]]]*Sqrt[Sum[x[[i]]^2,{i,1,n}]]; x[[1]]=x[[1]]+a; b=a*x[[1]];
h=Table[0,{n},{n}];
Do[h[[i,j]]=x[[i]]*x[[j]],{i,1,n},{j,1,n}]; h=IdentityMatrix[n]-h/b; MatrixForm[%]
说明：本程序用于求把向量x变为?||x||2e1的镜面反射矩阵H。程序执行后，先通过键盘输入向量x，程序即可给出镜面反射矩阵H。
程序中变量说明 h：存放镜面反射矩阵H x:存放向量x
求Householder矩阵（镜面反射矩阵）例题
设向量x={3,5,1,1}，试构造一个镜面反射矩阵H，使Hx=-||x||2e1, e1 ={1,0,0,0}. 解：执行求镜面反射矩阵H程序后在输入窗口中按输入向量:
{3,5,1,1}，然后用鼠标点击窗口的“OK”按扭，得如下输出结果。
T
T
-1/2 -5/6 -1/6 -1/6
-5/6 29/54 -5/54 5/54
-1/6 -5/54 53/54 -1/54
-1/6 -5/54 -1/54 53/54
这个结果就是所求的镜面反射矩阵H.
QR分解算法
1. Q?E, R?A
2. For k=1,2,..,n-1
取出矩阵R的第k列的ak+1 ,k,ak+2 k,ak+3 k,…,an k元素构造(n-k+1)?(n-k+1)阶镜面反射矩阵H(k)
将n阶单位矩阵的右下脚子阵用H(k) 取代构造n?n阶正交矩阵Qk
做Q? QkR , Q? QQk
注： E为n阶单位矩阵。该算法算出矩阵A的QR分解A=QR，其中Q存放正交矩阵, R存放上三角矩阵。
用Householder变换求矩阵QR分解程序
Clear[x,h,a];
h[x_]:=Module[{a,b,v,t,n},n=Length[x];
a=Sign[x[[1]]]*Sqrt[Sum[x[[i]]^2,{i,1,n}]]; v=Table[0,{n},{n}]; t=Table[0,{n}]; t=x;
t[[1]]=x[[1]]+a;
b=a*t[[1]];
Do[v[[i,j]]=t[[i]]*t[[j]],{i,1,n},{j,1,n}]; IdentityMatrix[n]-v/b];
a=Input[“输入矩阵A=”] n=Input[“输入矩阵A的阶n=”] hh=IdentityMatrix[n]; Do[
b=Table[a[[i,k]],{i,k,n}]; h1=h[b];
h2=IdentityMatrix[n];
Do[h2[[i+k-1,j+k-1]]=h1[[i,j]],{i,1,n-k+1},{j,1,n-k+1}]; a=h2.a; hh=hh.h2; Print[\ Print[\ {k,1,n-1}]
说明：本程序用于求矩阵A的QR分解。程序执行后，先通过键盘输入矩阵A，程序即可给出每次的镜面反射矩阵H和对应的正交矩阵。
程序中变量说明
h[x]：求出向量x的镜面反射矩阵H
a: 存放矩阵A和QR分解过程中的矩阵R
hh: 存放矩阵A阵QR分解过程中的正交矩阵Q
QR分解例题与实验
例5．设矩阵
??3?A??4???12??14130?313?19??135??13?11??13??求A的QR分解。
解：执行QR分解程序后在输入窗口中按输入向量
{{3,-14/13,-19/13},{4,0,5/13},{12,-3/13,-11/13}}，3
然后用鼠标点击窗口的“OK”按扭，得如下输出结果。
R={{-13, 6/13, 1}, {0, 5/13, 1}, {0, 12/13, 1}}
Q={{-3/13, -4/13,-12/13},{-4/13, 12/13, -3/13}, {-12/13, -3/13, 4/13}} R={{-13, 6/13, 1}, {0, -1, -17/13}, {0, 0, -7/13}} Q={{-3/13,164/169,-12/169},{-4/13,-24/169,-159/169},
{-12/13,-33/169,56/169}}
故得到A的QR分解为
??3?A??4???12??14130??31319??3???13?13?54????13??1311??12???13???1316416924?16933?16912????13169??159???0169??56????0169???613?10?1??17??13?7???13?
QR方法算法
1.输入矩阵A、最大迭代次数nmax 2.k?0
3.将A做QR分解得A= Q0 R0 4.交换乘积次序得A1 =R0 Q0
5.判断A1 是否本质上一致收敛于上三角形矩阵或分快三角形矩阵，如果收敛，则停止，否则做A?A1，k?k+1,转3
QR算法程序
h[x_]:=Module[{a,b,v,t,n},n=Length[x];
a=N[Sign[x[[1]]]*Sqrt[Sum[x[[i]]^2,{i,1,n}]],10]; v=Table[0,{n},{n}]; t=Table[0,{n}]; t=x;
t[[1]]=x[[1]]+a; b=a*t[[1]];
Do[v[[i,j]]=t[[i]]*t[[j]],{i,1,n},{j,1,n}]; IdentityMatrix[n]-v/b]; a=Input[“输入矩阵A=”]
n=Input[“输入矩阵A的阶n=”] nmax=Input[“输入迭代最大次数nmax=”] Do[ hh=IdentityMatrix[n];
Do[ b=Table[a[[i,k]],{i,k,n}]; h1=h[b];
h2=IdentityMatrix[n];
Do[h2[[i+k-1,j+k-1]]=h1[[i,j]],{i,1,n-k+1},{j,1,n-k+1}]; a=h2.a; hh=hh.h2, {k,1,n-1}];
Print[\迭代次数m=\ a=a.hh;
Print[\ {mm,1,nmax}]
说明：本程序用于求一般矩阵A的全部特征值。程序执行后，先通过键盘输入矩阵A、矩阵的阶数、迭代最大次数，程序即可给出每次迭代的矩阵序列，通过观察最后给出的矩阵
可以得到矩阵A的全部特征值。
程序中变量说明：
h[x]：求出向量x的镜面反射矩阵H函数
a:存放矩阵A和QR分解过程中的R和QR方法的迭代序列Ak hh: 存放矩阵A阵QR分解过程中的正交矩阵Q nmax: 存放QR方法的迭代允许的最大次数。 注：最大次数nmax可以调控。
QR方法例题与实验
例6．设矩阵
?5?1A???0??0?2020?5?321?1??2??3???2?用 QR方法求A的全部特征值。
解：执行QR方法程序后在输入窗口中按输入向量
{{5,-2,-5,-1}, {1,0,-3,2}, {0,2,2,-3}, {0,0,1,-2}}，4，20 然后用鼠标点击窗口的“OK”按扭，得如下输出结果。
迭代次数m=1 a=
4.61538 -5.95115 1.42869 1.25978 0.399704 1.94017 -2.52456 1.66248 1.18016?10-16 2.46957 -0.853801 3.19062 5.11023?10-17 6.6913?10-18 0.303869 -0.701754 迭代次数m=2 a=
4.11828 -3.46451 5.39263 -0.151083 0.299921 0.687205 -3.30064 -4.2376 3.64519?10-17 0.930157 1.11718 -0.711995 8.69753?10-18 4.45096?10-17 0.22015 -0.922664 迭代次数m=3
a=
3.87093 0.996828 6.20813 0.862851 0.0959 -0.352723 -2.44633 3.34172 2.14451?10 2.39185 2.56005 -2.23815 1.3935?10-18 3.94583?10-17 0.0700648 -1.07825 …………………………………………………………………. 迭代次数m=18
a=
4.00003 -6.2063 0.580683 -0.845242 0.0000160962 2.41472 -2.23612 -2.00627
-17
1.06607?10 2.68388 -0.414755 -3.49264 3.81442?10 6.06345?10 3.57803?10 -1. 迭代次数m=19 a=
4.00001 -3.71942 5.00213 0.845252 0.0000145278 1.07369 -3.89046 3.9383 6.8727?10-21 1.02954 0.926304 0.844604 1.77561?10 3.26111?10 2.58357?10 -1.
迭代次数m=20
a=
3.99999 0.777338 6.18475 -0.845265 5.4027?10-6 -0.426111 -2.59369 -3.4272 2.76352?10-21 2.32634 2.42612 2.1161
2.12415?10 2.46193?10 7.68638?10 -1.
观察输出结果可以明显看到矩阵序列虽然开始收敛规律不明显，但随着迭代次数的增加，其本质上收敛特点变得越来越明显，如在矩阵的a21，a31 ，a41 ，a42 ，a43对应的数列具有明显的收敛于零的特点。当迭代20次后，这些数相应变为5.4027?10-6，2.76352?10，2.12415?10 ，2.46193?10 ，7.68638?10 ，由此可以知道矩阵A有2个实特征值3.99999 和-1和两个复特征值，它们由二阶主子阵
-0.426111 -2.59369 2.32634 2.42612 的特征值得出，键入命令
Eigenvalues[{{-0.426111,-2.59369},{2.32634,2.42612}}]
可以得出两个复特征值{1. + 2. I, 1. - 2. I}。本题的4个特征值为{4., 1. + 2. I, 1. - 2. I, -1.}，可见收敛效果相当满意。
-21
-28
-22
-8
-28
-22
-8
-27
-22
-7
-27
-22
-7
-20
(k)(k)(k)(k)(k)