Python+Scipyでルンゲクッタ8次のDOP853(Dormand&Prince)を使う(その1) まずは準備編
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がすごく誤差が小さい。
ではこれを使っていろいろ計算してみる(続く)。
« Pythonの3Dのscatter+cmap+SciPyの特殊関数で水素原子の波動関数の電子分布を描いてみる。 | トップページ | Python+Scipyでルンゲクッタ8次のDOP853(Dormand&Prince)を使う(その2) やはり最初はローレンツ方程式。 »
「学問・資格」カテゴリの記事
- 高周波(RF・マイクロ波・ミリ波・5G)関連ニュース2021年2月16日 IEEE Microwave Magazineの特集はオールデジタルのRFID、Microwave JournalはEバンド ミリ波通信に衛星や気球を使う話、アメリカの半導体企業がバイデンに投資を迫る、(2021.02.17)
- カオスを生じる電気回路、Chua’s circuitをLTspiceで回路シミュレーションしてみる。(2021.02.19)
- Labyrinth Chaos(迷宮カオス)を生むThomas-Rössler方程式のパラメータbを色々変えて、Python+Scipyでルンゲクッタ8次のDOP853(Dormand&Prince)を使って計算してGIFアニメ(2021.02.16)
- フィッツヒュー・南雲 (FitzHugh-Nagumo) 方程式をPython+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)で計算。(2021.02.23)
- 「水晶振動子の等価回路計算」をカシオの高精度計算サイトkeisan.casio.jpの自作式としてUP! インピーダンスの大きさと位相がグラフ化できる。(2021.02.12)
« Pythonの3Dのscatter+cmap+SciPyの特殊関数で水素原子の波動関数の電子分布を描いてみる。 | トップページ | Python+Scipyでルンゲクッタ8次のDOP853(Dormand&Prince)を使う(その2) やはり最初はローレンツ方程式。 »
コメント