数值分析大作业QR分解.docx
《数值分析大作业QR分解.docx》由会员分享，可在线阅读，更多相关《数值分析大作业QR分解.docx（8页珍藏版）》请在冰点文库上搜索。
题目:
设计求取实数矩阵A的所有特征值及其特征向量的数值算法，并以矩阵
为例进行具体的求解计算。
一、算法分析：
一般而言，求取实数矩阵所有特征值的方法有雅克比法和QR分解法，两者都是变换法。
其中雅克比法只能求解对称矩阵的全部特征值和特征向量，而QR则可用于更一般的矩阵特征值的求解，结合反幂法可进而求出各个特征向量。
二、算法设计：
1、原始实矩阵拟上三角化
为了减少求特征值和特征向量过程中的计算量，对生成的矩阵A进行拟上三角化，得到拟上三角化矩阵A’
记A
（1）=A，并记A（r）的第r列到第n列的元素。
对于r=1,2,…,n-2执行
（1）若全为零，则令A（r+1）=A（r），转（5）；否则转
（2）。
（2）计算
令
（3）令，
（4）计算
（5）继续
算法执行完后，就得到与原矩阵A相似的拟上三角矩阵A（n-1）。
2、拟上三角矩阵QR分解的求原矩阵的全部特征值
记Ak是对拟上三角矩阵进行QR算法，产生的矩阵序列，A0是起始拟上三角矩阵，
对于k=1,2,…,n-1执行：
（1）分解
选取旋转矩阵，使，从而使第一列次对角元；再选取旋转矩阵，使，从而使第一列次对角元……如此进行下去，最多经过n-1步，必然变为上三角矩阵，即
必为正交矩阵，且满足
（2）计算
（3）上述两过程反复进行直到变为近似舒尔矩阵，对角线元素即为A的近似特征值。
3、带位移的QR方法求实矩阵全部特征值
一般情况下，QR分解的收敛速度比较慢，因此可通过仿乘幂法将原矩阵先进行一定长度的位移，可显著加速QR算法所得矩阵序列的收敛。
（A是原始矩阵的拟上三角矩阵）
（1）分解，其中即位移量，一般取A的某一特征值的近似值；实际计算通常取为的右下角元素，或取为右下角矩阵特征之中接近者。
（2）计算。
（3）重复过程
（1）
（2）直到变为近似舒尔矩阵，对角线元素即为A的近似特征值。
4、反幂法求实矩阵的特征向量
记通过QR分解得到A的某一特征值的近似值为p，反幂法步骤如下：
（1）三角分解。
（2）对k=1，2，…，做
①当时，令
当时，解
②解
③求绝对值最大的元素
④令
⑤当或小于规定的误差界时，停止计算，即为所求的特征向量，即为对应特征值的更为精确的取值。
三、程序设计：
1、对生成的矩阵A进行拟上三角化
利用hohd函数对A进行householder变换，得到A的拟上三角矩阵。
2、对拟上三角化后的矩阵进行QR分解和带位移的QR分解，求矩阵的全部特征值
在绝对误差界为的条件下，利用qrfenjie函数对拟上三角矩阵进行QR分解，利用dp_qrfenjie函数进行带位移的QR分解，比较两者的收敛速度。
3、反幂法求更精确的特征值和特征向量
利用vect函数和
（2）得到的特征值的近似值计算更为精确的特征值和对应的特征向量，是绝对误差界为。
4、输出相关结果。
四、源程序：
1、hohd函数
function[A]=hohd（a）
n=length（a）;
fori=1:
n-2
b=a（i+1:
n,i）;
ifb
（1）>=0
c=-sqrt（sum（b.^2））;
else
c=sqrt（sum（b.^2））;
end
rho=sqrt（2*c*（c-b
（1）））;
u1=（b-c*eye（n-i,1））/rho;
H1=eye（n-i）-2*u1*u1';
H=eye（n-i）;
H（i+1:
n,i+1:
n）=H1;
a=H*a*（H'）;
end
A=a
2、qrfenjie函数
functionA=qrfenjie（a）
n=length（a）;
A=a;
fori=1:
100
theta（n-1）=0;
c（n-1）=0;
s（n-1）=0;
Q=1;
forj=1:
n-1
theta（j）=atan（A（j+1,j）/A（j,j））;
c（j）=cos（theta（j））;
s（j）=sin（theta（j））;
P=eye（n）;
P（j,j）=c（j）;
P（j+1,j）=-s（j）;
P（j,j+1）=s（j）;
P（j+1,j+1）=c（j）;
invP=eye（n）;
invP（j,j）=c（j）;
invP（j+1,j）=s（j）;
invP（j,j+1）=-s（j）;
invP（j+1,j+1）=c（j）;
Q=Q*invP;
R=P*A;
A=R;
end
A=R*Q;
tor=max（abs（diag（A）-diag（a）））;
iftor>5*10^-5
a=A;
else
break;
end
end
i,Ak=A,lanmda=diag（A）'
3、dp_qrfenjie函数
functionA=dp_qrfenjie（a）
n=length（a）;
A=a;
fori=1:
100
A=a-a（n,n）*eye（n）;
theta（n-1）=0;
c（n-1）=0;
s（n-1）=0;
Q=1;
forj=1:
n-1
theta（j）=atan（A（j+1,j）/A（j,j））;
c（j）=cos（theta（j））;
s（j）=sin（theta（j））;
P=eye（n）;
P（j,j）=c（j）;
P（j+1,j）=-s（j）;
P（j,j+1）=s（j）;
P（j+1,j+1）=c（j）;
invP=eye（n）;
invP（j,j）=c（j）;
invP（j+1,j）=s（j）;
invP（j,j+1）=-s（j）;
invP（j+1,j+1）=c（j）;
Q=Q*invP;
R=P*A;
A=R;
end
A=R*Q+a（n,n）*eye（n）;
tor=max（abs（diag（A）-diag（a）））;
iftor>5*10^-5
a=A;
else
break;
end
end
i,Ak=A,lanmda=diag（A）'
4、vect函数
functionz=vect（a0,p）
n=length（a0）;
Z=zeros（n,n+2）;
forx=1:
n
a=a0-p（x）*eye（n）;
r=zeros（n,n）;
r（1,1:
n）=a（1,1:
n）;
l=eye（n）;
l（2:
n,1）=a（2:
n,1）/a（1,1）;
fori=2:
n
forj=i:
n
M=0;
fork=1:
i-1
M=l（i,k）*r（k,j）+M;
end
r（i,j）=a（i,j）-M;
end
ififorj=i+1:
n
N=0;
fork=1:
i-1
N=l（j,k）*r（k,i）+N;
end
l（j,i）=（a（j,i）-N）/r（i,i）;
end
else
break;
end
end
R=r;
L=l;
u=ones（n,1）;
mo=0;
invr=inv（R）;
invl=inv（L）;
fori=1:
100
y=invr*u;
p1=max（y）;
p2=max（-y）;
ifp1>p2
m=p1;
else
m=-p2;
end
z=y/m;
tor=abs（m-mo）;
iftor<5*10^-7
break;
else
mo=m;
end
u=invl*z;
end
Z（x,3:
n+2）=z';
Z（x,1）=i;
Z（x,2）=p（x）+1/m;
end
Z
五、运行结果：
1、执行命令：
>>clear,clc,a=[2001;0-1-24;0-213;1431];hohd（a）;
得到原始实数矩阵的拟上三角矩阵A。
A=
2.0000-1.00000.00000.0000
-1.00001.00005.0000-0.0000
0.00005.0000-2.2000-0.4000
0.0000-0.0000-0.40002.2000
2、执行命令：
>>clear,clc,A=[2-100;-1150;05-2.2-0.4;00-0.42.2];qrfenjie（A）;
得到迭代次数i，舒尔矩阵Ak，以及A的特征值矩阵lamda。
i=
38
Ak=
-5.90680.0315-0.00000.0000
0.03154.8972-0.0000-0.0000
0.0000-0.00002.21380.0007
0.00000.00000.00071.7958
lanmda=
-5.90684.89722.21381.7958
3、执行命令：
>>clear,clc,A=[2-100;-1150;05-2.2-0.4;00-0.42.2];dp_qrfenjie（A）;
得到用带位移的QR分解法所用的迭代次数i，舒尔矩阵Ak，以及A的特征值矩阵lamda。
i=
8
Ak=
-5.9068-0.0056-0.0000-0.0000
-0.00564.89730.00000.0000
-0.00000.00001.79580.0000
0.00000.00000.00002.2138
lanmda=
-5.90684.89731.79582.2138
4、执行命令：
>>clear,clc,a=[2001;0-1-24;0-213;1431],p=[-5.90684.89722.21381.7958];vect（a,p）;
得到Z矩阵，Z矩阵每一行的第一个数是迭代次数，第二个数是特征值，后面的数为特征向量。
Z=
Columns1through4
5.000000000000000-5.9068479421191640.1124195215529501.000000000000000
5.0000000000000004.8973023267404870.3451486545848380.505132030845888
5.0000000000000002.213757601733807-0.282974649968150-0.697610774659887
5.0000000000000001.7957880136448701.000000000000000-0.324049010823671
Columns5through6
0.675655845769651-0.888884062644966
0.5105418495906991.000000000000000
1.000000000000000-0.060487982528655
0.044562197436887-0.204211986355130
六、结果分析：
1、QR算法是求取一般矩阵所有特征值的一种非常有效的数值方法。
2、由2、3两步的运行结果看出带位移的QR算法与一般的QR算法相比具有更好的收敛效果，极大地减少了运算量。
3、通过反幂法只需比较小的计算量就可以得到精度更高的特征值，以及对应的特征向量。
升级会员 