彈簧振子數值解

從上圖的結構，我們可以得到運動學方程

這是一組（還有兩個端點的方程）二階微分方程，我們可以得到方程的特解為

將其帶入運動學方程可得

整理可以得到如下形式

下面，我們就用數值方法來計算方程組的解

from

numpy

import

zeros

，

empty

from

pylab

import

plot

，

show

# Constants

N

=

20

C

=

1。0

m

=

1。0

k

=

6。0

omega

=

2。0

alpha

=

2

*

k

-

m

*

omega

*

omega

# Set up the initial values of the arrays

A

=

zeros

（［

N

，

N

］，

float

）

for

i

in

range

（

N

-

1

）：

A

［

i

，

i

］

=

alpha

A

［

i

，

i

+

1

］

=

-

k

A

［

i

+

1

，

i

］

=

-

k

A

［

0

，

0

］

=

alpha

-

k

A

［

N

-

1

，

N

-

1

］

=

alpha

-

k

v

=

zeros

（

N

，

float

）

v

［

0

］

=

C

# Perform the Gaussian elimination

for

i

in

range

（

N

-

1

）：

# Divide row i by its diagonal element

A

［

i

，

i

+

1

］

/=

A

［

i

，

i

］

v

［

i

］

/=

A

［

i

，

i

］

# Now subtract it from the next row down

A

［

i

+

1

，

i

+

1

］

-=

A

［

i

+

1

，

i

］

*

A

［

i

，

i

+

1

］

v

［

i

+

1

］

-=

A

［

i

+

1

，

i

］

*

v

［

i

］

# Divide the last element of v by the last diagonal element

v

［

N

-

1

］

/=

A

［

N

-

1

，

N

-

1

］

# Backsubstitution

x

=

empty

（

N

，

float

）

x

［

N

-

1

］

=

v

［

N

-

1

］

for

i

in

range

（

N

-

2

，

-

1

，

-

1

）：

x

［

i

］

=

v

［

i

］

-

A

［

i

，

i

+

1

］

*

x

［

i

+

1

］

# Make a plot using both dots and lines

plot

（

x

）

plot

（

x

，

“ro”

）

show

（）

#最近在嘗試數值求解晶體的聲子色散譜和態密度，沒有找到可用程式碼，於是就從彈簧振子開始探索。今天終於碼完正確計算fcc結構聲子譜和dos的程式碼，大快人心。先把這個彈簧振子的小程式搬上來。