【Python】ローレンツ方程式を実装して可視化する

Posted: , Category: Python , 物理学

ローレンツ方程式 (Lorentz equation) は、非線形微分方程式(nonlinear ordinary differential equation, non-linear ODE)で表される方程式で、特定のパラメータと初期値を与えることで、カオス的な振る舞いをすることがよく知られています。

今回は、このローレンツ方程式をPythonで実装し、上図のように図示する方法を解説します。

常微分方程式をPythonで実装して可視化するといった流れについて理解できると思います。

本記事の内容
  • ローレンツ方程式について解説
  • Pythonでローレンツアトラクタを実装
  • ローレンツアトラクタをアニメーションで表示する

ローレンツ方程式の定義

まず最初にローレンツ方程式の定義を示します。

ローレンツ方程式は、次の3つの非線形微分方程式(non-linear ODE)で表現されます。

ローレンツ方程式の定義
\begin{equation}
\begin{split}
\frac{dx}{dt} &=  s ( y  -x)\\
\frac{dy}{dt} &=  rx - y  -x z\\
\frac{dz}{dt} &= xy - bz
\end{split}
\end{equation}

ここで、3つの変数$s, r, b$はパラメータである。

ローレンツアトラクタ

(1)と3つのパラメータ$s, r, b$で定義されるローレンツ方程式において、所定の値が与えられたときに、ローレンツ・アトラクタとなります。ローレンツアトラクタはカオス的な振る舞いをすることが知られています。

カオス的な振る舞いとは、「決定論的非周期的流れ」であり、パラメータと初期位置が定まると、その後の動きが決定するにもかかわらず、その挙動が周期的にならず、非周期的な挙動をします。

なかなかイメージも湧かないと思うので、ローレンツ方程式を実装し、ローレンツアトラクタの挙動を確認してみましょう。

ローレンツ方程式を実装する

(1)式で示したローレンツ方程式を、Pythonを用いて実装してみます。

ここで実装している内容は、他の微分方程式を実装したり可視化したりする際に、よく利用できるので、自分で手を動かしてみることをお勧めします。

今回は、ローレンツ方程式を実装して、可視化するところまでやりたいので、Pythonの行列計算ライブラリであるnumpyと3D可視化用のモジュールであるmatplotlibとそのモジュールであるAxes3Dを読み込みます。

import numpy as np
import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

続いて、(1)のローレンツ方程式の関数lorenzと、lorenz関数から時間発展した x, y, z の値を得るスクリプトはこのようになります。

# ローレンツ方程式
def lorenz(x, y, z, s=10, r=28, b=8/3):
    dot_x = s * (y - x)
    dot_y = r * x - y - x * z
    dot_z = x * y - b * z
    return dot_x, dot_y, dot_z

# 時間感覚とステップ数
dt = 0.01
n_steps = 5000

# データの軌跡を保存する配列
xs = np.empty(n_steps + 1, dtype="float64")
ys = np.empty(n_steps + 1, dtype="float64")
zs = np.empty(n_steps + 1, dtype="float64")

xs[0], ys[0], zs[0] = (0., 1., 1.05)

for i in range (n_steps):
    dot_x, dot_y, dot_z = lorenz(xs[i ], ys[i ], zs[i] )
    xs[i + 1] = xs[i] +  ( dot_x * dt )
    ys[i + 1] = ys[i] +  ( dot_y * dt )
    zs[i + 1] = zs[i] +  ( dot_z * dt )

非線形の常微分方程式の解を数値計算で得る場合、時刻$i$の値とその時の傾きから、次の時刻$i+1$の値をえる実装になっています。

値を得ることができたらプロットを行います。

from matplotlib import colors
cm = plt.cm.get_cmap('RdYlBu')
cs = colors.Normalize(0, n_steps)

fig = plt.figure(figsize=(7, 7))
ax = Axes3D(fig)

ax.xaxis.pane.fill = False
ax.yaxis.pane.fill = False
ax.zaxis.pane.fill = False

for i in range(n_steps):
    ax.plot(xs[i:i+2], ys[i:i+2], zs[i:i+2], color=plt.cm.viridis(cs(i)), lw=0.5)
    
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")

plt.show()

3次元プロットし、さらに、時間発展と共に色が変わるようにしています。これを図示すると、次のようなグラフを得ることができます。

これがいわゆる、ローレンツアトラクタですね。綺麗に得ることができています。

ローレンツアトラクタのアニメーションを実装する

最後に、先ほど描写したローレンツアトラクタのアニメーションを実装してみましょう。

コードの全文は以下のようになります。解像度をあげているので、描写に少し時間がかかりますが、実行すると次のようなアニメーションを得ることができます。

from matplotlib.animation import FuncAnimation
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d, Axes3D
from matplotlib import colors

def lorenz(x, y, z, s=10, r=28, b=8/3):
    dot_x = s * (y - x)
    dot_y = r * x - y - x * z
    dot_z = x * y - b * z
    return dot_x, dot_y, dot_z


# 時間感覚とステップ数
dt = 0.01
n_steps = 5000

# データの軌跡を保存する配列
xs = np.empty(n_steps + 1, dtype="float64")
ys = np.empty(n_steps + 1, dtype="float64")
zs = np.empty(n_steps + 1, dtype="float64")

xs[0], ys[0], zs[0] = (0., 1., 1.05)

for i in range (n_steps):
    dot_x, dot_y, dot_z = lorenz(xs[i ], ys[i ], zs[i] )
    xs[i + 1] = xs[i ] +  ( dot_x * dt )
    ys[i + 1] = ys[i ] +  ( dot_y * dt )
    zs[i + 1] = zs[i ] +  ( dot_z * dt )
    
cm = plt.cm.get_cmap('RdYlBu')
cs = colors.Normalize(0, n_steps)

fig = plt.figure(figsize=(5, 5))
ax = Axes3D(fig)

ax.xaxis.pane.fill = False
ax.yaxis.pane.fill = False
ax.zaxis.pane.fill = False

ax.set_xlim([-20, 20])
ax.set_ylim([-30, 30])
ax.set_zlim([0, 60])

def animate(i):
    artist = ax.plot(xs[i:i+2], ys[i:i+2], zs[i:i+2], color=plt.cm.viridis(cs(i)), lw=0.5)
    return artist

ani = FuncAnimation(fig, animate, n_steps, interval=10, blit=True)

ani.save('lorenz.gif', writer="imagemagick",dpi=100)

【広告】
統計学的にあなたの悩みを解決します。
仕事やプライベートでお悩みの方は、ベテラン占い師 蓮若菜にご相談ください。

機械学習と情報技術