使用Matlab实现:幂法、反幂法(原点位移)

使用Matlab实现:幂法、反幂法(原点位移)

幂法

例题

使用幂法,计算下面矩阵的主特征值及对应的特征向量。

A = [ 2 − 1 0 − 1 2 − 1 0 − 1 2 ] A= \left[\begin{array}{ccc} 2 & -1 & 0\\ -1 & 2 & -1\\ 0 & -1 & 2 \end{array}\right] A=210121012

实现

幂法的数学迭代公式为:

取 v ( 0 ) ≠ 0 , α ≠ 0 , 令 u ( 0 ) = v ( 0 ) 取 v^{(0)} \neq 0,\alpha \neq 0, 令 u^{(0)} = v^{(0)} v(0)̸=0α̸=0u(0)=v(0)

v ( 1 ) = A v ( 0 ) = A u ( 0 ) v^{(1)} = A v^{(0)} = A u^{(0)} v(1)=Av(0)=Au(0)

u ( 1 ) = v ( 1 ) m a x ( v ( 1 ) ) = A v ( 0 ) m a x ( A v ( 0 ) ) u^{(1)} = \frac{v^{(1)}}{max(v^{(1)})} = \frac{A v^{(0)}}{max(A v^{(0)})} u(1)=max(v(1))v(1)=max(Av(0))Av(0)

v ( 2 ) = A u ( 1 ) = A v ( 1 ) m a x ( A v ( 0 ) ) = A 2 v ( 0 ) m a x ( A 2 v ( 0 ) ) v^{(2)} = A u^{(1)} = \frac{Av^{(1)}}{max(A v^{(0)})} = \frac{A^2 v^{(0)}}{max(A^2 v^{(0)})} v(2)=Au(1)=max(Av(0))Av(1)=max(A2v(0))A2v(0)

u ( 2 ) = v ( 2 ) m a x ( v ( 2 ) ) = A 2 v ( 0 ) m a x ( A 2 v ( 0 ) ) u^{(2)} = \frac{v^{(2)}}{max(v^{(2)})} = \frac{A^2 v^{(0)}}{max(A^2 v^{(0)})} u(2)=max(v(2))v(2)=max(A2v(0))A2v(0)

依次类推。

使用Matlab实现:

format long g;
v0 = [1;1;1];
u0 = [1;1;1];
A = [2,-1,0;-1,2,-1;0,-1,2];
v = A * u0;
u = v / norm(v, inf);
i = 0;
while norm(u - u0, inf) >= 1e-5
    u0 = u;
    v = A * u0;
    u = v / norm(v, inf);
    i ++;
end;
norm(v, inf)
i
u

求解得:

i = 8 , u ( 9 ) = ( 0.70711 , − 1 , 0.70711 ) T , λ = m a x ( v ( 9 ) ) = 3.41422 i = 8,u^{(9)} = (0.70711, -1, 0.70711)^T,\lambda = max(v^{(9)}) = 3.41422 i=8u(9)=(0.70711,1,0.70711)Tλ=max(v(9))=3.41422

反幂法

例题

已知下列矩阵有特征值 λ \lambda λ 的近似值 p = 4.3 p = 4.3 p=4.3 ,用原点位移的反幂法,求对应的特征向量 u u u ,并改善 λ \lambda λ

A = [ 3 0 − 10 − 1 3 4 0 1 − 2 ] A= \left[\begin{array}{ccc} 3 & 0 & -10\\ -1 & 3 & 4\\ 0 & 1 & -2 \end{array}\right] A=3100311042

实现

反幂法,是基于幂法,推导出来的一个求最小特征值的方法。由于特征值的性质,可以用来求某个近似值的精确特征值。

其数学迭代方法如下:

首先,对矩阵 A − p I A - p I ApI进行三角分解,便于以后求方程组的解。

A − p I = L U A - p I = L U ApI=LU

v ( k ) v^{(k)} v(k) 时,相当于求两个三角方程组。

{ L y ( k ) = u ( k − 1 ) U v ( k ) = y ( k ) \begin{cases} L y^{(k)} = u^{(k - 1)}\\ U v^{(k)} = y^{(k)}\\ \end{cases} {Ly(k)=u(k1)Uv(k)=y(k)

又有:

u ( k ) = v ( k ) m a x ( v ( k ) ) u^{(k)} = \frac{v^{(k)}}{max(v^{(k)})} u(k)=max(v(k))v(k)

使用Matlab实现:

A = [3,0,-10;-1,3,4;0,1,-2];
I = eye(3,3);
p = 4.3;
u0 = [1;1;1];
v = inv(A - p * I) * u0;
u = v / norm(v, inf);
i = 0;
while norm(u - u0, inf) > 1e-5
    u0 = u;
    v = inv(A - p * I) * u0;
    u = v / norm(v, inf);
    i ++;
end;
i
u
x = p + 1 / norm(v, inf)

求解得:

i = 6 , u ( 7 ) = ( − 0.96606 , 1 , 0.15210 ) T , λ = p + 1 m a x ( v ( k ) ) = 4.57447 i = 6,u^{(7)} = (-0.96606, 1, 0.15210)^T,\lambda = p + \frac{1}{max(v^{(k)})} = 4.57447 i=6u(7)=(0.96606,1,0.15210)Tλ=p+max(v(k))1=4.57447

  • 27
    点赞
  • 181
    收藏
    觉得还不错? 一键收藏
  • 7
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 7
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值