PythonでRungeKutta

Python微分方程式を解きます。アルゴリズムはRungeKutta法を使います。
忘れた頃に定期的に必要になってくるのでここに関数化して書いておきます。
下にコードがあるのでコピペして使ってください。

numpyをimportしないと使えません。

以下のsolve_RK()を使うと \frac{d y}{d t} = f(t,y) のような一階微分方程式をyについて解くことができます。
この微分方程式に従ってRungeKutta法でtに関して発展させて行きます。

solve_RK() 説明

solve_RK(t0,tmax,N,func,y0)

引数

  • t0 : tの初期値
  • tmax : tの最大値
  • N : 計算ステップ数
  • func : f(t,y) pythonは関数がオブジェクトとして引数に使えます。
  • y0 : yの初期値

戻り値

t発展の経過がlist型で返ってきます。 ndarrayではありません。

コード

コピペして使ってください

def solve_RK(t0,tmax,N,func,y0):
    dt = (tmax-t0)/(N-1.)
    tlist = np.arange(t0,tmax,dt)
    y=[y0]
    i = 0
    for t in tlist:
        k1=dt*func(t,y[i])
        k2=dt*func(t+dt/2.0,y[i]+k1/2.0)
        k3=dt*func(t+dt/2.0,y[i]+k2/2.0)
        k4=dt*func(t+dt,y[i]+k3)
        y.append(y[i] + (k1+2.0*k2+2.0*k3+k4)/6.0)
        i += 1
    return y

使用例

 f(t,y) = \frac{dy(t)}{dt} = y

tの範囲と初期条件
t:0\rightarrow10, y(0) = 1

厳密解  y(t) = e^{t}

コード例

RunkeKuttaを関数かしたおかげで1行でRungeKuttaができてます。これからはこれをコピペして物理の時間発展などの微分方程式を解いていこう。

参考:

JupyterをGistで共有する。ブログに貼れるぞ。 - ましろのログ