Python+Scipyでルンゲクッタ8次のDOP853(Dormand&Prince)を使う（その１） まずは準備編

PythonでWikipediaに乗っているカオスの写像リストを描いて遊んでいたが、そろそろネタ切れ、、、

そこで写像から常微分方程式に移る。とはいってもodeintを使うと面白くないので、このサイトでいつもやっている（ExcelのVBAに移植した）

DOP853（Dormand&Prince, ルンゲクッタ8次）を使ってみよう。

(参考リンク)

 Dormand-Prince(ルンゲ・クッタ8次）の有名なルーチンDOP853をExcel VBAに移植

Dormand-Prince(ルンゲ・クッタ8次）DOP853をExcel VBAに移植(説明編)

では準備として単振動する

dx/dt=y

dy/dt=-x

を計算してみよう。これを100周期（2π*100)繰り返して厳密解のcos(t),sin(t)と比較してみる。

Pythonのプログラムはこんな感じ。odeintを使っている方は、そもそも関数のtとxの並びが反転しているとか、1ステップずつ計算しないといけないとか、初期値の与え方がまた並びが反転しているとか謎が多い、、、

import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt

def f1(t, x): #odeintのときとt,xの並びが逆
    x_dot = x[1]
    y_dot = -x[0]
    return [x_dot, y_dot]

t0=0
tmax = 2.*np.pi*100
N=10000

x0=[0., 1.0]
solver=ode(f1)
solver.set_integrator('dop853')
solver.set_initial_value(x0,t0) #なぜか関数と並びが逆

t=np.linspace(0, tmax, N)
sol= np.empty((N, 2))
sol[0] = x0

k=1
while solver.successful() and solver.t < tmax:
    solver.integrate(t[k])
    sol[k] = solver.y
    k+= 1

# Plot
fig = plt.figure(figsize=(12,12))

plt.xlabel("Time")
plt.ylabel("Numerical - Exact")
plt.title("DOP853")
plt.grid()
plt.plot(t,sol[:,0]-np.sin(t))
plt.plot(t,sol[:,1]-np.cos(t))

plt.show()

結果はこちら。100周期繰り返しても誤差が1e-14オーダー。

ではodeintを使うと？

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def f2(x, t):
    x_dot = x[1]
    y_dot = -x[0]
    return [x_dot, y_dot]

t0=0
tmax = 2.*np.pi*100
N=10000

t=np.linspace(0, tmax, N)

x0=[0., 1.0]
x=odeint(f2, x0, t)
# Plot
fig = plt.figure(figsize=(12,12))
plt.xlabel("Time")
plt.ylabel("Numerical - Exact")
plt.title("odeint")
plt.grid()
plt.plot(t,sol[:,0]-np.sin(t))
plt.plot(t,sol[:,1]-np.cos(t))

plt.plot(t,x[:,0]-np.sin(t))
plt.plot(t,x[:,1]-np.cos(t))
plt.show()

簡単なのだが誤差は0.00001くらい。

やっぱりこの手のStiffじゃない問題ならDOP853がすごく誤差が小さい。

ではこれを使っていろいろ計算してみる（続く）。