« Pythonの3Dのscatter+cmap+SciPyの特殊関数で水素原子の波動関数の電子分布を描いてみる。 | トップページ | Python+Scipyでルンゲクッタ8次のDOP853(Dormand&Prince)を使う(その2) やはり最初はローレンツ方程式。 »

2020年6月26日 (金)

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オーダー。

Sin_dop853

では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くらい。

Sin_odeint

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

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

« Pythonの3Dのscatter+cmap+SciPyの特殊関数で水素原子の波動関数の電子分布を描いてみる。 | トップページ | Python+Scipyでルンゲクッタ8次のDOP853(Dormand&Prince)を使う(その2) やはり最初はローレンツ方程式。 »

学問・資格」カテゴリの記事

コメント

コメントを書く

(ウェブ上には掲載しません)

« Pythonの3Dのscatter+cmap+SciPyの特殊関数で水素原子の波動関数の電子分布を描いてみる。 | トップページ | Python+Scipyでルンゲクッタ8次のDOP853(Dormand&Prince)を使う(その2) やはり最初はローレンツ方程式。 »

最近の記事

2021年3月
  1 2 3 4 5 6
7 8 9 10 11 12 13
14 15 16 17 18 19 20
21 22 23 24 25 26 27
28 29 30 31      
フォト
無料ブログはココログ