计算物理-chapter数值方法

矩阵对角化与本征值问题线性方程组：

AX

BAX

B本征值问题：KAKX

K

K

XKmj-1mjmj+1本征值问题：AX

X线性方程线组性：方A程X组：BKKKKmj-1Kmjmj+1......jjjjFNd

2

xdt

2dt

2dt

2d

2

x2d

x

m

d

2

xNdt

2

K

(x

x

)

K

(x

x

)F

m

K

(x1

x2

)

K

(x3

x2

)F2

m

21F

m1

Kx1

K

(x2

x1)

K

(xN

1

xN

)

KxN

.j

1j

1......jjjjd

2

x

Ndt

2d

2

xdt

2dt

22d

xd

2

xFN

m

K

(

x

x

)

K

(

x

x

)F

m

K

(

x1

x2)

K

(

x3

x2)2F

m

2dt

2F1

m 1

Kx1

K

(

x2

x1

)

K

(

xN

1

xN

)

Kx

N

.j

1j

1itx

j

Aje

...

0

0

1

N

N

A

Aj

...

Aj

2

...

0

...0

00

A

2

1

0

0

...1

2

1

0

...

01

2

1

......

A

A1

m2

...K2

A2

A1

n

n

nnxxanjaain

xi

aij

xi

a

x1

n1

an

2

ai1a2n

x

x2

2

a1n

x1

a1

ja2

j

a11

a1221a

a22

AX

X

0

0

0

n

xi

0

x

2

x1

an

2anjannan1ainai1a21a22a1n2na1

ja2

ja12

xij

a

a

a11

det(

A

I

)

0求矩阵特征值与特征向量的方法1.

乘幂法

(最大特征值)反幂法(最小特征值)Jacobi方法(对称矩阵)QR方法(更一般的算法)乘幂法：求矩阵的按模最大的特征值与相应的特征向量。基本思想：通过迭代产生向量序列，由此计算特征值和特征向量的近似值。乘幂法0

1A1

1特征值为:1=1.618032=

0.61803

11X

(0)做如下迭代：X

(k

1)

AX

(k

)

A

X

(2)

0

11

1

1

11

2

0

11

2

1

1

2

3

0

1

2

3

1

13

5

0

13

51

15

8X

(3)

A

X

(3)X

(4)1.66667

1.61.51.666671.61.6251.6251.61538

A

X

(1)X

(2)21.5

1318

1

0

15

8

(4)(5)X

A

X

(5)(6)X

A

X112X

(1)

A

X

(0)x(k

)

/

x(k1)1x(k

)

/

x(k

1)2

2

1

1

21

34

0 113

21(6)(7)X

A

X1.61538

1.61905

1

134

55

01

21

34

A

X

(7)

X

(8)1.61905

1.61765

1

A

X

(8)X

(9)1.61765

1.61818

1

189

144

0

155

89

(9)(10)X

A

X1.61808

1.61798

1

1144

233

0

189

144

(10)(11)X

A

X1.61798

1.61806

1 1

233

377

0 1144

233(11)(12)X

A

X1.61806

1.61803

1

1233

3770

1144

233(11)(12)X

A

X1.61806

1.61803特征值为:1=1.618032=

0.618031

10

1A任取初时向量X(0)R

xk

xk

2

1

k

n

iix

xk

x

n

AX

(k

)

Ak

1

xk

1

xk

1

2xk

1

1X

(k

1)

i

i

xlim

xkk

1k

为A按模最大特征值X(k+1)为对应的特征向量(k

0,1,2,...)maxk

11ink

1Y

(k

1)

X

(k

1)

/

mm

乘幂法的规范运算：任取初时向量X(0)R，通常取X(0)＝（1,1,…,1）T迭代过程为：X

(k

1)

AY

(k

)mk

1i|

x(k

1)

|例：用乘幂法求矩阵2

2

1 0

2

110A

0的按模最大的特征值和相应的特征向量。取：

X

(0)

(0,0,1)T

,

103.解：X

(0)

Step1.(1)

Tm1Y

(1)X

(0,

0.5,

1),(0,0,1)TY

(1)

AX(0)

(0,1,

2)T

,

m1

2,X

(2)

(0.2,

0.8,1)T

,m2Y

(2)Step2.

Y

(2)

AX

(1)

(0.5,

2,

2.5)T

,

m2

2.5,Y

(8)

AX

(7)

(2.7650948,

2.9981848,

2.9990924)m8

2.9990924X

(8)

(0.9219772,

0.9996973,

1)Y

(9)

AX

(8)

(2.8436517,2.9993946,2.9996973)m9

2.9996973由m

m

2.9996973

2.9990924

0.0006049

1039

8故1

2.9996973相应的特征向量为：u1

(2.8436517,2.9993946,2.9996973)精确解：1

3,2

2,1

1,1对应的本征向量:(1,1,1)T迭代解：1

2.9996973u1

(2.8436517

,2.9993946 ,2.9996973

)

T反幂法反幂法：计算矩阵按模最小的特征值及特征向量的方法基本思想：逆矩阵的特征值是原矩阵特征值的倒数，特征向量相同AX

X;A1

X

1

X反幂法算法：任取初时向量X(0)R，通常取X(0)＝（1,1,…,1）T迭代过程为：(k

0,1,2,...)ki|

x(k

)

|kmY

(k

)

X

(k

)

/

m1inmaxAX

(k

1)

Y

(k

)1/

mkX

(k

1)

A1Y

(k

)例：求如下矩阵最小特征值及对应的特征向量

4

1

9

3A

KY(k)X(k+1)0110.19047660.23809510.8

10.1809520.2761920.65517210.1740560.30377730.57297310.1701420.31943440.53263610.1682210.32711750.51425310.1673450.33061860.50615810.166960.3321670.50264910.1667930.33282980.50113710.1667210.333117A-1的按模最大特征值为0.333117A-1特征向量为(0.501137,1)TA的按模最小特征值为1/0.333117=3.0019A特征向量为(0.501137,

1)TJacobi方法（实对称矩阵的全部特征根与特征向量）定理：P为n阶可逆阵，则A与P－1AP相似，相似阵有相同的特征值；若A对称，则存在正交矩阵

TQ=I)，使得n

TQ AQ

12计算如下矩阵的特征值和相应的特征向量

01A

cos

sinsin

cosB

cos

BT

AB

sinsin1

01

cos

sin

010

2

sin

coscos

sincos

sincos

2

sin2222当=/4时：2

/

2 2

/

2 2

/

2 2

/

2

B1

00

1B AB

T1

010A

A对应于1=-1的特征向量为：A对应于2=1的特征向量为：2

/

2

v1

2

/

22

/

22

/

2v2

当=/4时：B

2

/

2 2

/

2 2

/

22

/

2

0

11

0B AB

A的特征值为1=-1,

2=-1T构造一系列特殊形式的正交阵Q1,...,Qn对A作正交变换使得对角元素

逐次增加，非对角元变小。当非对角元已经小得无足轻重时，可以近似认为对角元就是A的所有特征值。Jacobi法基本思路：Givens旋转变换：

1Q(

p,

q,)

1

sincossincosp列q列p行q行记：A

(a

)

,

B

QT

(

p,

q,)AQ(

p,

q,)

(b

)ijij则：

2qp

pqbi

b

a

cos

2

app

abpqcos2

abpp

a

pp

sin2

apqpq

sin

2sin

2pp

ppiin

asin

222cos,i

p,

qbbip

bpi

api

cos

aqi

sin,i

p,

q变换的目的是为了减少非对角元的分量，因此：2bpq

bqp

a

acos

2

app

sin

2

0pq2apq

app,

t

tans

a记s

0

,

t

2则

t

2ts

1

0

的按模较小根所以：s

0 ,

1sincos

d1

t

2

c1

t

2t12pppqb

a

cos

2

app

a

sin

2

0tan2

2

cot

2tan1

0

appa2apqcot

222222sin

2ppppppppqiqiippi

pisin

2a

a

aa

sin

a

pq

sin

2a

sin

a

cos

b

a

cos

a

sin,i

p,

qcos,i

p,

qb

bppqipiqiiqqiip

pi

pibpqb

bqp

a

apqtcaos

2pp

pqbpq

bqp

0

tapqbpp

appb

b

da

ca

,

i

p,

q,

i

p,

qbii

api

sin

ab

b

ca

dasin

22pp

pqpp

pppqpp

ppqi

piqi

iqapp

ab

acos

2b

asin

2

a

cos2

a

sin

2cos2

a

sin

2

a

sin

2pqb

asin

a

cos,i

p,

q

b

abbip

bpi

api

cos

aqi

sin,i

p,

qJacobi迭代算法：取p,q使i

j

max

aijapq，则A(k

1)

QT

(

p,

q,)A(k

)Q(

p,

q,)定理：若A对称，则A(k

1)

diag{

,,

}1

n解

记

A(0)=A,

取p=1,q=2,apq

=a12

=2,于是有(0)

(0)例：用Jacobi

方法计算对称矩阵的全部特征值6

21

4

2

2A

251

0.25122a(0)a(0)

a(0)s

11

22t

sgn(s)

/(|

s

|

1

s2

)

0.780776

0.7882061

t

21cos

0.6154121

t

2sints

0

,

1s

0

,

t

2

2ts

1

0t

所以再取p=2,q=3,apq(1)=a23(1)=2.020190,类似地可得10.615412

00.788206

00

1

0

0.788206

0

0.61541200cos

sincos01

pqsR

R

()

in1

10.9616

2.43844800

6.561552

2.020190

0.961

2.020190A(1)

RT

A(0)

R

0.631026

0.724794

2.43844800.631026

0.724794

8.3203860

4.241166A

(

2)

0.595192

2.18318500.595192

0

8.320386

0.2096140.209614

4.496424A

(3)0

0.0200488.377576

0.2086530.208653

4.496424

0

0.020048

2.125995A

(

4)4.485239

0.020019

0.001073

0.020019

2.125995

0.0010738.38876100A(5)

0

0.001072

2.1258250

0.0010728.388761

0.0000090.000009

4.485401A

(6)0

0

8.388761

0.0000090.000009

4.485401

00

2.125825A

(7)从而A的特征值可取为12.125825,

28.388761,

34.485401特征向量为旋转矩阵R=R1R2…的相应的各列矢量-

2.0-1.0

0.01.0

2.02.03.04.05.04.03.02.01.00.01.02.03.04.05.04.03.02.01.00.01.02.03.04.05.04.03.02.0-1.00.01.02.03.04.05.04.03.0-

2.0-1.00.01.02.03.04.05.04.0-

3.0-

4.0

-

3.05.0-1.0

-

2.0

4.0

5.0

4.0

3.0

2.0

1.0

0.0

-1.03.0

4.0

5.0

4.0

3.0

2.0

1.0

0.0-

2.0-1.00.0

1.0

2.0

3.0

4.03.05.04.33.02.01.00.0-1.0-