用 Runge-Kutta 法解微分方程

日前從圖書館借了一本書 From Calculus to Chaos ，這本書主要是在介紹 動態系統 的一些特性，附了很多例子， 其中第4章有介紹如何用 Runge Kutta 方法來求得微分方程的數值解， 方法是出乎我意料之外的簡單，於是找了一個例子來練習，考慮下面的 Duffing equation

為了用 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 函式庫。