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

2020年6月29日 (月)

Python+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)を使う(その3) Aizawa Attractor 

今回はみかん袋(あみあみ)を丸めたようなAizawa attractor。

http://analog-ontology.blogspot.com/2015/06/the-aizawa-attractor-iv.html

dx/dt = (z-b)*x - d*y
dy/dt = d*x+(z-b)*y
dz/dt = c+a*z-(z^3)/3 - (x^2 +y^2)*(1+e*z)+f*z*x^3

のような形で、計算結果はこちら。

Aizawa

プログラムはこんな感じで。

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


def aizawa(t, x, a, b, c, d, e, f): #odeintのときとt,xの並びが逆
    x_dot = (x[2]-b)*x[0] - d*x[1]
    y_dot = d*x[0]+(x[2]-b)*x[1]
    z_dot = c+a*x[2]-(x[2]**3)/3 - (x[0]**2 +x[1]**2)*(1+e*x[2])+f*x[2]*x[0]**3
    return [x_dot, y_dot, z_dot]

t0=0
tmax = 1000
N=100000

x0=[0.1, 0.0, 0.0]
solver=ode(aizawa)
solver.set_integrator('dop853')
solver.set_initial_value(x0,t0) #なぜか関数と並びが逆
solver.set_f_params(0.95,0.7,0.6,3.5,0.25,0.1)

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], lw=0.5)
ax.set_xlabel("X Axis")
ax.set_ylabel("Y Axis")
ax.set_zlabel("Z Axis")
ax.set_title("Aizawa Attractor(DOP853)")

plt.show()

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

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

コメント

コメントを書く

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

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

最近の記事

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