カオスの軌跡のプロット

[やること]
微分方程式で表現されるカオスをPythonで解いて、軌跡をプロットします。

[カオスについて]
カオスは式で明確に定義されているのに、初期値の選び方で将来の状態が予測できない現象のことです。
参考文献[1][2]がわかりやすいです

[結果1]
ローレンツ方程式

f:id:wada0421514:20200329140013p:plain

p,r,b=10,28,8/3

初期値 x,y,z=1,1,1
f:id:wada0421514:20200329135713p:plain
初期値 x,y,z=1.0001,1,1
f:id:wada0421514:20200329135717p:plain
初期値 x,y,z=1.001,1,1
f:id:wada0421514:20200329135721p:plain
初期値 x,y,z=1.01,1,1
f:id:wada0421514:20200329135727p:plain

横軸に時間、縦軸にxを取ったプロットです。
初期値が少し変わるだけで軌跡が大きく変わることがわかります。

f:id:wada0421514:20200329140511p:plain


[結果2]
レスラー方程式
f:id:wada0421514:20200329140830g:plain

a,b,c=0.2,1.0,5.7

初期値 x,y,z=1.0,0.0,0.0
f:id:wada0421514:20200329140958p:plain
初期値 x,y,z=8.01,0.0,0.0
f:id:wada0421514:20200329141003p:plain
初期値 x,y,z=8.02,0.0,0.0
f:id:wada0421514:20200329141007p:plain

[プログラム]
微分方程式を解くライブラリodeintについては、参考文献[3]がわかりやすいです

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

#レスラー方程式
def rossler(var,t,a,b,c):
    x_next=-var[1]-var[2]
    y_next=var[0]+a*var[1]
    z_next=b*var[0]-(c-var[0])*var[2]
    return np.array([x_next,y_next,z_next])

#ローレンツ方程式
def lorentz(var,t,p,r,b):
    x_next=-p*var[0]+p*var[1]
    y_next=-var[0]*var[2]+r*var[0]-var[1]
    z_next=var[0]*var[1]-b*var[2]
    return np.array([x_next,y_next,z_next])

if __name__ == "__main__":

    #レスラー方程式のパラメータと初期値
    #abc=(0.2,1.0,5.7)
    #xyz_first=[1.0,0.0,0.0]

    #ローレンツ方程式のパラメータと初期値
    prb=(10,28,8/3)
    xyz_first=[1.01,1,1]

    #ステップ数
    t=np.linspace(0,1000,100000)

    #微分方程式を解く
    #var=odeint(chaos,xyz_first,t,args=abc)
    var=odeint(lorentz,xyz_first,t,args=prb)

    #グラフの枠を作っていく
    fig = plt.figure()
    ax = Axes3D(fig)

    #軸にラベルを付ける
    ax.set_xlabel("X")
    ax.set_ylabel("Y")
    ax.set_zlabel("Z")

    #プロットして表示
    ax.plot(var[:, 0], var[:, 1], var[:, 2], marker="o", linestyle='None')
    plt.show()


[参考文献]
[1]カオスアトラクタとは
http://www.gifu-nct.ac.jp/elec/deguchi/sotsuron/hisaki/node9.html#eqresura
[2]カオス
http://www.aoni.waseda.jp/yuuka/Sim/Chaos.html
[3]Python運動方程式を解く(odeint)
https://qiita.com/binaryneutronstar/items/ad5efa27fd626826846f