数值分析大作业QR分解.docx
- 文档编号:78711
- 上传时间:2022-10-02
- 格式:DOCX
- 页数:8
- 大小:103.74KB
数值分析大作业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
ifi forj=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,
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 作业 QR 分解