Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う(その11) 中二病用語、ブルースカイ・カタストロフィを生じるGavrilov Shilnikov modelを計算してGIFアニメに。
千葉逸人さんが命名者が中二病だと言ったので有名な?Blue-Sky catastropheを生じるGavrilov Shilnikov modelをPythonとScipyのDOP853で計算してみる。
スカラーペディアはご本人が書いているので詳しい。
http://www.scholarpedia.org/article/Blue-sky_catastrophe
こんな式です。
計算結果はこちら:クリックでアニメが始まります。
ソースはこんな感じで。
%matplotlib notebook
import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation
from matplotlib.animation import PillowWriter
INTERVAL = 0.001 # [s] interval time
FRAME_INTERVAL = 1000 * INTERVAL # [msec] interval between frames
FPS = 1000/FRAME_INTERVAL # frames per second
def bluesky(t, x,myu,eps): #odeintのときとt,xの並びが逆
x_dot = x[0]*(2.0+myu-10.0*(x[0]**2+x[1]**2))+x[2]**2+x[1]**2+2.0*x[1]
y_dot = -x[2]**3-(1+x[1])*(x[2]**2+x[1]**2+2.0*x[1])-4.0*x[0]+myu*x[1]
z_dot = (1+x[1])*(x[2]**2)+x[0]**2-eps
return [x_dot, y_dot, z_dot]
t0=0
tmax = 1000
N=20000
#x0=[0.012277918,-2.356078578,0.018241293]
x0=[0,0,0]
solver=ode(bluesky)
solver.set_integrator('dop853',atol=1.E-12,rtol=1.E-12)
solver.set_initial_value(x0,t0) #なぜか関数と並びが逆
solver.set_f_params(0.456,0.0357)
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
def update(i):
plt.cla() #現在描写されているグラフを消去
ax.plot(sol[:i,0], sol[:i,1], sol[:i,2],c='Blue',lw=1,alpha=0.8)
ax.scatter(sol[i,0], sol[i,1], sol[i,2],c='Blue')
ax.set_xlim(-1.5,1.5)
ax.set_ylim(-1.5,1.5)
ax.set_zlim(-1.5,1.5)
ax.set_xlabel("X Axis")
ax.set_ylabel("Y Axis")
ax.set_zlabel("Z Axis")
ax.set_title("Gavrilov Shilnikov model : Blue-sky catastrophe(DOP853)")
plt.grid()
# Plot
fig = plt.figure(figsize=(12,12))
ax = fig.gca(projection='3d')
ax.view_init(elev=10, azim=60)
ani = FuncAnimation(fig, update, frames=range(0,len(t),100), interval=FRAME_INTERVAL, blit=True)
plt.show()
ani.save("bluesky.gif", writer = 'pillow',fps=FPS)
« Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う(その10) 7体問題の昴問題(the Pleadis) | トップページ | AX(伊坂幸太郎さん)を読んだ。恐妻家の殺し屋の連作短編集だが、途中で驚く、、、そして最後にあの伏線が効いてるとは、、、 »
「学問・資格」カテゴリの記事
- 高周波・RFニュース 2024年12月11日 5G AmericasがセルラーネットワークでのAIのホワイトペーパー、GSAが2024年の5Gレビュー、NordicがIoTに適したThingy:91 X発表、Wi-Fi 6GHzが世界でどうなっているか、Huaweiが5G GaNパワーアンプで有利、Samsung Galaxy M34 5G分解(2024.12.11)
- 高周波・RFニュース2024年12月9日 iFixitがDJI Neo分解、TechInsightsがApple Pencil Pro分解、QualcommのNeurIPS 2024でのAI技術発表、IntelのIEDM 2024での発表、 Nokiaの7GHz帯の6G、Analog DevicesのPhased Array Antennaのホワイトペーパー、ZDTが史上二番目の売上高(2024.12.09)
- 高周波・RFニュース 2024年12月6日 NGMNが無線パフォーマンス評価フレームワーク発行、5GAAがC-V2Xのロードマップ発行、Marvellの3nm 1.6Tbps PAM4インターコネクト、Nokiaの2.4Tbps光伝送、Silicon Labsの低消費電力モジュール、Xiaomi 14T Pro分解動画(2024.12.06)
- 高周波回路シミュレータQucsStudioがuSimmicsに名称変更し、バージョンも4.8.3から5.8にアップデートされた。Qucsと区別するためだそうだ。また、Pythonの高周波用ライブラリscikit-rfもv1.5.0にバージョンアップされていた(2024.12.04)
- 日経サイエンス2025年1月号の特集 和算再発見の佐藤賢一さんの記事「算聖 関孝和の実像」に出てきた矢高に対する円弧の2乗の近似式をカシオの高精度計算サイトkeisan.casio.jpの自作式として作った。ものすごい精度であることがよくわかる。(2024.12.03)
« Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う(その10) 7体問題の昴問題(the Pleadis) | トップページ | AX(伊坂幸太郎さん)を読んだ。恐妻家の殺し屋の連作短編集だが、途中で驚く、、、そして最後にあの伏線が効いてるとは、、、 »
コメント