« Python+Scipyでルンゲクッタ8次のDOP853(Dormand&Prince)を使う(その1) まずは準備編 | トップページ | Qualcommのリファレンスデザインは富士通が作ってるが、全く売れる見込みもないということでスマホが出始めのころの国内メーカの変化に対応できなくて没落していった事例をちょっとだけ思い出す、、、 »

2020年6月27日 (土)

Python+Scipyでルンゲクッタ8次のDOP853(Dormand&Prince)を使う(その2) やはり最初はローレンツ方程式。

昨日で大体使い方が分かったので、まずは一番メジャーなローレンツ方程式か。

dop853で計算してみる。本当はrtolとatolも設定しないといけないが手抜きで。

Lorenz_dop853_20200623205801

おなじみの形が得られた。プログラムはこんな感じで。

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


def lorenz(t, x, s, r, b): #odeintのときとt,xの並びが逆
    x_dot = s*(x[1] - x[0])
    y_dot = r*x[0] - x[1] - x[0]*x[2]
    z_dot = x[0]*x[1] - b*x[2]
    return [x_dot, y_dot, z_dot]

t0=0
tmax = 100
N=100000

x0=[0.1, 0.1, 0.1]
solver=ode(lorenz)
solver.set_integrator('dop853')
solver.set_initial_value(x0,t0) #なぜか関数と並びが逆
solver.set_f_params(10,28,8/3)

t=np.linspace(0, tmax, N)
sol= np.empty((N, 3))
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))
ax = fig.gca(projection='3d')

ax.plot(sol[:,0], sol[:,1], sol[:,2])
ax.set_xlabel("X Axis")
ax.set_ylabel("Y Axis")
ax.set_zlabel("Z Axis")
ax.set_title("Lorenz Attractor(DOP853)")

plt.show()

---

odeintでデフォルトで計算してみた結果も。やっぱり結構違うね。

Lorenz_odeint

 

 

« Python+Scipyでルンゲクッタ8次のDOP853(Dormand&Prince)を使う(その1) まずは準備編 | トップページ | Qualcommのリファレンスデザインは富士通が作ってるが、全く売れる見込みもないということでスマホが出始めのころの国内メーカの変化に対応できなくて没落していった事例をちょっとだけ思い出す、、、 »

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

コメント

コメントを書く

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

« Python+Scipyでルンゲクッタ8次のDOP853(Dormand&Prince)を使う(その1) まずは準備編 | トップページ | Qualcommのリファレンスデザインは富士通が作ってるが、全く売れる見込みもないということでスマホが出始めのころの国内メーカの変化に対応できなくて没落していった事例をちょっとだけ思い出す、、、 »

最近の記事

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      
フォト
無料ブログはココログ