Python+Scipyでルンゲクッタ8次のDOP853(Dormand&Prince)を使う(その2) やはり最初はローレンツ方程式。
昨日で大体使い方が分かったので、まずは一番メジャーなローレンツ方程式か。
dop853で計算してみる。本当はrtolとatolも設定しないといけないが手抜きで。
おなじみの形が得られた。プログラムはこんな感じで。
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でデフォルトで計算してみた結果も。やっぱり結構違うね。
« Python+Scipyでルンゲクッタ8次のDOP853(Dormand&Prince)を使う(その1) まずは準備編 | トップページ | Qualcommのリファレンスデザインは富士通が作ってるが、全く売れる見込みもないということでスマホが出始めのころの国内メーカの変化に対応できなくて没落していった事例をちょっとだけ思い出す、、、 »
「学問・資格」カテゴリの記事
- 時代に逆行してCOBOL(GnuCOBOL)を学んでみる(6) 31桁の10進演算で複素数を計算するライブラリを作ってみる。外部ファイルのサブルーチンを呼び出すCALL文とUSINGの練習。(2023.06.07)
- 原著で持っているダンズモアさんのHandbook of Microwave Component Measurements With Advanced VNA Techniques(2nd ed.) が邦訳されてマイクロ波部品測定ハンドブック(第2版)として出版されてる!Pozarさんのマイクロ波工学に続き、最近名著の翻訳が増えてるのはいい傾向!(2023.06.03)
- 平均が0でなく相関がある正規分布の比(商)Z=X/Yの分布 Fieller-Hinkleyの式(平均が0ならコーシー分布)をカシオの高精度計算サイトkeisan.casio.jpにあげていたが式にミスがあったのでやり直し。Pythonでモンテカルロして全然違う値になったのでもっと新しい論文探した。(2023.06.02)
- 時代に逆行してCOBOL(GnuCOBOL)を学んでみる(5) 中二病のような名前のブルースカイカタストロフィ(Blue-sky catastrophe)を生じるGavrilov Shilnikov modelをドルマン・プリンス法(ode45)で計算する。テキストベース(アスキーアート)で!(2023.06.01)
- TensorFlow(Keras)でモデルをsaveで保存してload_modelで読み込むときに突然エラー(utf8で読めないとかディレクトリが存在しないとか)が出始めた。なんで?といろいろやると、単に日本語がフォルダ名に使われているときにエラーになるだけだった…(Windowsネイティブの場合)(2023.05.28)
« Python+Scipyでルンゲクッタ8次のDOP853(Dormand&Prince)を使う(その1) まずは準備編 | トップページ | Qualcommのリファレンスデザインは富士通が作ってるが、全く売れる見込みもないということでスマホが出始めのころの国内メーカの変化に対応できなくて没落していった事例をちょっとだけ思い出す、、、 »
コメント