カオスの描写とPython

はじめに

微分方程式や漸化式を用いることで、過去や現在の状況から近未来の状況を予測することができる。しかし、かなり先の未来は、それらの数式によっては予想不可能な程不規則に変化する場合がある。それをカオスといい、初期条件が少し変化しただけで結果が大幅に変化してしまうことがある。今回は、質点の運動方程式（微分方程式）であるローレンツ方程式と、ロジスティック写像という漸化式を用いることで、複雑な模様を描写する。

ローレンツ方程式

導入

ローレンツ方程式は、気象における大気の循環や対流の運動を簡略したモデルの微分方程式であり、以下のような単純な式で表すことができる。

\begin{equation}
\left\{ \,
    \begin{aligned}
    & \frac{dx}{dt} = -px+py\\
    & \frac{dy}{dt} = -xz+rx \\
    & \frac{dz}{dt} = xy-bz
    \end{aligned}
\right.
\end{equation}

ただし、質点の座標を$(x,y,z)$とし、$p,r,b$を定数とする。

ここで、この単純な微分方程式は、初期条件や微分方程式の係数が少し変化するだけで、質点の運動の挙動が大幅に変化してしまう。また、運動の様子ですらかなり複雑である。

プログラム

以下のようなプログラムを書く。

python Lorenz.py

import numpy as np
import matplotlib.pyplot as plt
import math
import japanize_matplotlib

#定数設定
p=10
r=28
b=8/3
#初期条件 (x,v,a)=(位置,速度,加速度)
x=np.array([1,0,0])
#斜め上に発射
v=np.array([0,2,1])

#3次元のグラフ（散布図）を描くための設定
fig = plt.figure()
ax = fig.add_subplot(projection='3d')

#初期時間
t=0
#微小時間の時間間隔(小さくするほど精度は上がるが計算時間が比例して増大する)
delta_t=0.01
while(t<10):

    #位置の微小変化
    x=x+v*delta_t
    #微分方程式
    v=np.array([-p*x[0]+p*x[1],-x[0]*x[2]+r*x[0]-x[1],x[0]*x[1]-b*x[2]])

    #時間の更新
    t=t+delta_t
    #print(x)
    #粒子の位置を記録
    ax.scatter(x[0], x[1], x[2], color='blue')
#グラフを保存する
plt.savefig("Lolez.png")
#グラフの表示
plt.show()

これを実行すると以下のような画像が出力される。

このように、大気の循環のような複雑な軌道を示す。

ロジスティック写像（分岐図）

導入

以下のような数列を考える。

x_{n+1}=ax_n(1-x_n)

プログラム

ここで、$2.0<a<4.0$の範囲で定数$a$を動かしていった場合の$x$を描写するプログラムを書くと以下のようになる。

python logi.py

import numpy as np
import matplotlib.pyplot as plt
import math
import japanize_matplotlib

n=30000

x=0.8
a =np.linspace(2.0, 4.0, n)

for i in range(n):
  x=a*x*(1-x)
  plt.plot(a[i],x[i],"c.", markersize=1.7)

plt.savefig("logi.png")
plt.show()

これを実行すると以下のようになる。

このように、生物学の進化系統図のようなグラフを得ることができる。
ただし、このプログラムを見れば分かるが、漸化式の計算に再帰を使用せずforループを用いている。これは、再帰計算だと計算量が膨大になってしまうためである。

まとめ

今回はPythonを用いて、微分方程式や漸化式を解き、運動などを描写するプログラムを書いた。これによって、カオスな模様を描写することができた。このような模様は、初期条件や問題設定を少し変化させることによって多種多様に変化させることができると考えられる。