Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う(その4) ローレンツ96モデル 36変数の方程式
ローレンツ方程式(ローレンツモデル)と言えば3変数のものが有名ですが多変数のローレンツ96というモデルもある。
https://en.wikipedia.org/wiki/Lorenz_96_model
F=8でカオスになるとか。
dx_i/dt=(x_i+1 - x_i-2)*x_i-1 -x_i +F
では計算。Wikipediaではodeintだがdop853で。やっぱりちょっと波形が違うな、、、
プログラムはこちら:
import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
N=36
F=8
def lorenz96(t, x): #odeintのときとt,xの並びが逆
"""Lorenz 96 model."""
# Compute state derivatives
d = np.zeros(N)
# First the 3 edge cases: i=1,2,N
d[0] = (x[1] - x[N-2]) * x[N-1] - x[0]
d[1] = (x[2] - x[N-1]) * x[0] - x[1]
d[N-1] = (x[0] - x[N-3]) * x[N-2] - x[N-1]
# Then the general case
for i in range(2, N-1):
d[i] = (x[i+1] - x[i-2]) * x[i-1] - x[i]
# Add the forcing term
d = d + F
# Return the state derivatives
return d
t0=0.
tmax=30.
dt=0.01
x0 = F * np.ones(N) # Initial state (equilibrium)
x0[19] += 0.01 # Add small perturbation to 20th variable
t = np.arange(t0, tmax, dt)
solver=ode(lorenz96)
solver.set_integrator('dop853')
solver.set_initial_value(x0,t0) #なぜか関数と並びが逆
sol= np.zeros((len(t),N))
sol[0] = x0
k=1
while solver.successful() and solver.t < tmax-dt:
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("Lorenz96 Attractor(DOP853)")
plt.show()
« Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う(その3) Aizawa Attractor | トップページ | Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う(その5) ピタゴラスの三体問題を計算する。rtolとatolを設定しないと無茶苦茶になる。 »
「学問・資格」カテゴリの記事
- 「数学がゲームを動かす! ゲームデザインから人工知能まで」を読んだ。面白い!パックマンのアルゴリズムやドラクエの計算式、ドンキーコングはベルレ法でジャンプ、カルマンフィルタ、遺伝的アルゴリズム、セガの線形代数本を書かれた方は理論物理出身など話題が豊富。(2025.05.13)
- 高周波・RFニュース 2025年5月12日 IEEE Microwave Magazineでミリ波ガラス基板・超伝導・液晶などの記事、Pythonの高周波ライブラリscikit-rfがv1.7.0に、Megamagneticsの希土類を使わないミリ波サーキュレータ、CoilcraftのRFインダクタホワイトペーパー(2025.05.12)
- 高周波・RFニュース 2025年5月9日 QorvoがスマートファクトリーにUWB利用の解説、KYOCERA AVXがUHF,VHF向け高方向性カプラ発表、Wireless Broadband Allianceが企業向けWi-Fi 7トライアル結果報告、NXPが第三世代イメージングレーダプロセッサ発表(2025.05.09)
- 高周波・RFニュース 2025年5月8日 6Gにはテラヘルツどころか7-14GHzも向いてないという意見、フィンランド企業が6G推進コンソーシアムRF ECO3発表、Infineonが7P7T内蔵7ch-LNA発表、SpirentがOctoboxにWi-Fi 6/7自動テスト追加(2025.05.08)
« Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う(その3) Aizawa Attractor | トップページ | Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う(その5) ピタゴラスの三体問題を計算する。rtolとatolを設定しないと無茶苦茶になる。 »
コメント