幂法,反幂法求解矩阵最大最小特征值及其对应的特征向量(DOC)

数值计算解矩阵的按模最大最小特征值及对应的特征向量

一.幂法

1. 幂法简介：

当矩阵A 满足一定条件时，在工程中可用幂法计算其主特征值(按模最大)及其特征向量。矩阵A 需要满足的条件为： (1) 的特征值为A i n λλλλ,0||...||||2

1≥≥≥>

(2) 存在n 个线性无关的特征向量，设为n x x x ,...,,21 1.1计算过程：

i n

i i i u x

x αα,1

)

0()0(∑==，有对任意向量不全为0，则有

1

11111221

12111

1

1

11

1

011)()(...u u a u a u λu λαu αA x A Ax x k n n k n k k n

i i

k i i n

i i i k )(k (k))(k αλλλλλα++++=+=+++≈?

?

????+++======∑∑ 可见，当||1

2

λλ越小时，收敛越快；且当k 充分大时，有1)11

11)111

11λαλαλ=??????==+++(k )(k k

(k k )(k x x u x u x ，对应的特征向量即是)(k x 1+。

2 算法实现

.

,, 3,,1 , ).5()

5(,,,,||).4();max(,).3()

(max(;0,1).2(,).1()()

()

(停机否则输出失败信息转置若转否则输出若计算最大迭代次数，误差限，初始向量输入矩阵βλβεβλβλε←+←<<-←←=

←←k k N k y x Ay x x abs x y k N x A k k k

3 matlab 程序代码

function [t,y]=lpowerA,x0,eps,N) % t 为所求特征值，y是对应特征向量

k=1;

z=0; % z 相当于λ

y=x0./max(abs(x0)); % 规范化初始向量

x=A*y; % 迭代格式

b=max(x); % b 相当于β

if abs(z-b)

t=max(x);

return;

end

while abs(z-b)>eps && k

k=k+1;

z=b;

y=x./max(abs(x));

x=A*y;

b=max(x);

end

[m,index]=max(abs(x)); % 这两步保证取出来的按模最大特征值

t=x(index); % 是原值，而非其绝对值。

end

4 举例验证

选取一个矩阵A，代入程序，得到结果，并与eig(A)的得到结果比较，再计算A*y-t*y，验证y是否是对应的特征向量。结果如下：

结果正确，表明算法和代码正确，然后利用此程序计算15阶Hilb矩阵，与eig(A)的得到结果比较，再计算A*y-t*y，验证y是否是对应的特征向量。设

置初始向量为x0=ones(15,1)，结果显示如下

可见，结果正确。得到了15阶Hilb 矩阵的按模最大特征值和对应的特征向量。

二.反幂法

1.反幂法简介及其理论

在工程计算中，可以利用反幂法计算矩阵按模最小特征值及其对应特征向量。其基本理论如下，与幂法基本相同：

x x A x A x x Ax λ

λλ1

)(11

==?=--则，

由，可知，A 和A -1的特征值互为倒数，

求A 按模最小特征值即求A -1的按模最大特征值,取倒数即为A 的按模最小特征

值所以算法基本相同，区别就是在计算

)

1()()1()(1-)1()()1()1(;,+++++===k k k k k k k k x LU A y Ax y A x Ay x x 分解，来计算做对具体计算时，变换为而是时，不是令

2.算法实现

.

, );

4( ,))

(max(, ,1,).7(),

7(,,,1

,||).6( ),max().5(),,().4().3(,

))

(max(,0 ,1).2(,,,,).1(000停机否则输出失败信息转置若否则转停机输出若解方程组作三角分解置最大迭代次数误差限初始向量输入矩阵x abs x

y k k N k y x z Ux y Lz y LUx LU

A x abs x

y k N x A ←←+←<<-←←====←←←λλλ

ελλμ

λμλε

3 matlab 程序代码

function [s,y]=invpower(A,x0,eps,n) % s 为按模最小特征值，y 是对应特征向量 k=1; r=0; % r 相当于0λ

y=x0./max(abs(x0)); % 规范化初始向量 [L,U]=lu(A); z=L\y; x=U\z; u=max(x);

s=1/u; % 按模最小为A -1按模最大的倒数.

if abs(u-r)

while abs(u-r)>eps && k

y=x./max(abs(x)); z=L\y; x=U\z;

u=max(x); end

[m,index]=max(abs(x)); % 这两步保证取出来的按模最大特征值

s=1/x(index); % 是原值，而非其绝对值。 end

4 举例验证

同幂法一样，选取一个矩阵A ，代入程序，得到结果，并与eig(A)的得到结果比较，再计算 A*y-t*y ，验证y 是否是对应的特征向量。

可见结果正确，然后利用此程序计算15阶Hilb 矩阵，eig(A)的得到结果比较，再计算 A*y-s*y ，验证y 是否是对应的特征向量。设置初始向量为

x0=ones(15,1)，结果显示如下

可见，结果真确。得到了15阶Hilb矩阵的按模最大特征值和对应的特征向量。

三. 计算条件数

矩阵A的条件数等于A的范数与A的逆的范数的乘积，即cond(A)=‖A‖·‖A^(-1)‖，对应矩阵的3种范数，可以定义3种条件数。函数cond(A,1)、cond(A)或cond(A inf)是判断矩阵病态与否的一种度量，条件数越大表明矩阵的病态程度越大.

的最大和最小特征值

矩阵

分别为

范数，即

这里我们选择矩阵的

A

A

A

cond

T

2

1

2

1,

,

)

(

2λ

λ

λ

λ

=

, 而如果A为对称矩阵，如Hilb矩阵，A

A T的最大最小特征值，分别为A的最大最小特征值的平方。所以cond(A) 为A的最大最小特征值得比值。对于本例中的15阶Hilb矩阵来说，利用上面计算结果得其条件数(选择第二种条件数）为：3.0934e+017；这与直接利用cond(A)得到的结果：2.5083e+017在同一

数量级，再次表明了上述算得得最大最小特征值的正确性，同时又表明Hilb 矩阵是病态矩阵。

四.Aitken 商加速法

1.简介与原理

{}{}.

,??:2)( 2)(,,

;0lim ,122

1122

1222

1211加速法这种方法称为逼近用有充分大时当线性收敛即且收敛与若Aitken a a a

a a a a a a a x x x x x x y a

a a

a a a a a k a c a a a

a a a k k k

k k k k k n

n n n n n n k k k k k k

k k k =+---≈?+---

=--≈--≠=--+++++++++++++∞→

同幂法和反幂法计算最大和最小特征值类似，如果计算最大特征值，则迭代格式为)()

1(k k Ay x

=+；计算最小特征值时，迭代格式为

)1()()(1)1(,+-+==k k k k Ax y y A x 即。

2.算法实现

计算按模最大特征值算法如下：

.

,,),

3()

max(

,,1 , , ,,).6(),6(,,,).5(,

2)().4(,)max( , ).3(,

)

max( ,0.1 ,0 ,0 ,1).2(,,,),().1(0120100

122

0102010停机输出失败信息否则转置若否则转停机输出若计算置计算置最大迭代次数误差限初始向量输入y x x

k k N k y x Ay x x x

y k N x a A ij ??+???<<-+---=?=======λλααααλελλααααααλαλααε

类似幂法和反幂法可以写出按模最小特征值算法，此处不再赘述。

3.matlab 程序代码

function [r,y]=aitken(A,x0,eps,n) % r按模最大特征值,y为对应特征向量

k=1;

α

a0=0; % a 相当于0

α

a1=1; % a1 相当于1

r0=1; % 相当于2中的0λ

y=x0./max(abs(x0)); % 规范化初始向量

x=A*y;

α

a2=max(abs(x)); % a2相当于2

r=a0-(a1-a0)^2/(a2-2*a1+a0); % 相当于λ

if (a2-2*a1+a0)==0 % 若上式中分母为0，则迭代失败，返回disp "初始向量迭代失败"

return;

end

if abs(r-r0)

end

while abs(r-r0)>eps && k

k=k+1;

a0=a1;

a1=a2;

r0=r;

y=x./max(abs(x));

x=A*y; % 迭代格式

a2=max(abs(x));

if (a2-2*a1+a0)==0 % 若分母为0，则迭代失败，返回

return;

end

r=a0-(a1-a0)^2/(a2-2*a1+a0);

[m,index]=max(abs(eig(A))); % 以下代码保证取出来的按模最大特征值aa=eig(A); % 是原值，而非其绝对值。

if aa(index)>0 ||aa(index)==0

r=r;

e lse

r=-r;

end

end

end

类似可得按模最小特征值和特征向量的代码如下：与上面类似，所不同的只是迭代格式不同.

function [r,y]=invaitken(A,x0,eps,n)

k=1;

a0=0;

a1=1;

r0=1;

y=x0./max(abs(x0));

[L,U]=lu(A); % 迭代格式的不同

z=L\y;

x=U\z;

a2=max(abs(x));

r=a0-(a1-a0)^2/(a2-2*a1+a0);

i f (a2-2*a1+a0)==0

disp "初始向量迭代失败"

return;

end

if abs(r-r0)

end

while abs(r-r0)>eps && k

k=k+1;

a=b;

b=c;

r0=r;

y=x./max(abs(x));

z=L\y;

x=U\z;

a2=max(abs(x));

if (a2-2*a1+a0)==0

return;

end

r=a0-(a1-a0)^2/(a2-2*a1+a0);

end

[m,index]=min(abs(eig(A))); % 以下代码保证取出来的按模最大特征值aa=eig(A); % 是原值，而非其绝对值。

if aa(index)>0 ||aa(index)==0

r=1/r;

else

r=-1/r;

end

end

4.计算Hilb矩阵特征值

此处不再举例，而是直接应用于15阶Hilb矩阵，初始向量选为ones(15,1)，结果如下，并将结果与幂法和反幂法得到结果比较

这与幂法得到的特征值和特征向量一致，表明算法和代码正确；同理，最小特征值结果如下：

这与反幂法得到的结果一致，表明结果正确。

五，对称矩阵的Rayleigh 商加速法 1.简介与原理

商的为关于则称为对称矩阵，设ayleigh )(,0R A x

x Ax

x x R x A T T =≠

原理如下：

)

max(,)()

()()()()max(0)()

()

(1)

()()()

1()()

()()

1()()

()

(i i k k k k k T k k T k k k k k k k i i i i x x y

y R Rayleigh y y x y y R Ay x

x x y X X AX X →

→????

?????====≠++λλλ其中商加速法。

这称为左乘上式有：用的特征向量，即为设

2. 算法实现

.

, );

3( ,))

(max(, ,1,).5(),5(,,,,||).4(,*).3(,

))(max(,0 ,1).2(,

,,,).1(000停机否则输出失败信息转置若否则转停机输出若置最大迭代次数误差限初始向量输入矩阵x abs x

y r r k k N k y r r r y

y x

y r y A x x abs x

y r k N x A T

T ←←+←<<-←←←←←εε

3. Matlab 程序代码

function [r,y]=rayleigh(A,x0,eps,n) % r 是特征值，y 是特征向量 k=1; r0=0;

y=x0./max(abs(x0));

x=A*y; % 迭代格式计算新的x r=dot(y,x)/dot(y,y); % Reyleigh 商 if abs(r-r0)

while abs(r-r0)>eps && k

y=x./max(abs(x)); x=A*y;

r=dot(y,x)/dot(y,y); end end

类似得计算按模最小特征值的Rayleigh 商加速法，如下： function [r,y]=invrayleigh(A,x0,eps,n) k=1;r0=0;

y=x0./max(abs(x0));

[L,U]=lu(A); % 迭代格式不同 z=L\y; x=U\z;

r=dot(y,x)/dot(y,y); if abs(r-r0)

while abs(r-r0)>eps && k

r0=r;

y=x./max(abs(x));

z=L\y;

x=U\z;

r=dot(y,x)/dot(y,y);

end

r=1/r;

end

4. 计算Hilb矩阵特征值

此处不再举例，而是直接应用于15阶Hilb矩阵，初始向量选为ones(15,1)，结果如下，并将结果与幂法和反幂法得到结果比较

这与幂法得到结果一致，表明算法和代码正确。

同理，最小特征值如下：

与反幂法得到结果一致，表明代码和算法正确。