為了用 Runge Kutta 方法,要先把上式改寫為
Duffing equation 是二階非線性時變微分方程,二階是因為 \ddot{x},
非線性是因為 x^3,而時變是因為 cos \Omega t。我是用
Python 來寫,用
Matplotlib 來做繪圖的工作,
注意在安裝 Matplotlib 之前需先安裝
Numpy。
所使用的參數是 k=0.05, A=7.5, \Omega=1,
上圖的曲線是 x 對 t 做圖,兩條曲線的初始條件 (x, y) 分別為 (3.003, 3) 跟 (3, 3),
即兩個初始條件的差別只在第一個元素,而且只差了 0.003,可是從上圖可知大概在 20 秒之後,
兩條曲線開始有明顯的差異,這就是混沌系統 (
chaotic system
) 的特性之一,
系統對於初始條件的變化非常敏感,初始條件的微小變動會造成系統動態的巨大改變。
程式碼如下:
(我是用 GeSHi 做 syntax highlight,程式碼在
chaos.py )
from pylab import *
from math import cos
k = 0.05
Omega = 1
A = 7.5
def f(t, state):
[x, d_x] = state
y = d_x
d_y = -k * d_x - x*x*x + A*cos(Omega*t)
return array([y, d_y])
def Runge_Kutta(a, b, h, initial):
z = []
state = initial
hh = h/2
for t in arange(a, b+h, h):
z.append(state[0])
k1 = f(t, state)
k2 = f(t+hh, state+hh*k1)
k3 = f(t+hh, state+hh*k2)
k4 = f(t+h, state+h*k3)
state += h/6 * (k1 + 2*k2 + 2*k3 + k4)
return z
a = 0
b = 50
h = 0.01
z = Runge_Kutta(a, b, h, [3, 3])
y = Runge_Kutta(a, b, h, [3.003, 3])
t = [i for i in arange(a, b+h, h)]
plot(t, z, t, y)
xlabel('t')
ylabel('x')
show()
程式碼很簡單,我只介紹兩個 Numpy 提供的函式,分別是 arange() 跟 array(),
Python 內建 range() 的參數型態只能是整數,Numpy 提供的 arange() 則允許參數型態是浮點。
學過線性代數都知道,向量空間中最基本的兩個運算分別是純量乘法跟向量加法,比如
2 [1, 2, 3] = [2, 4, 6]
[1, 2, 3] + [1, 0, 1] = [2, 2, 4]
可是 Python 內建的 list 型態,並沒提供這兩種運算,但是如果是用 array() 所宣告的型態,
則可直接使用這兩種運算。
有了 Numpy 提供的 arange() 跟 array (),程式碼是不是看起來就跟
維基百科
所列的演算法相差無幾呢?這就是 Python 程式語言的迷人之處,
趕快來學
Python 吧!
註:我是在
Yukuan's Blog
發現 Matplotlib 這個好用的 Python 函式庫。