はじめに
数値的に解いたとき、カオス的な解を示す微分方程式としてよく知られるものの1つにローレンツ方程式があります。式自体はシンプルですが、解の振る舞いがとても興味深いと思います。
インターネットにたくさん資料が存在しているので、自身でも実装してみました。やっぱり勉強で手を動かすことは重要だと思います。
微分方程式
ローレンツ方程式は変数$x,y,z$とパラメータ$\sigma,\rho,\beta$で構成され、以下で示されます。
\frac{dx}{dt} = \sigma (y-x)
\frac{dy}{dt} = x(\rho-z) - y \\
\frac{dz}{dt} = xy - \beta z
プログラム
数値解法はオイラー法を用います。
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def lorenz_eq(x, y, z, dt, sigma=10, beta=8/3, rho=28):
xn = x + (sigma*(y - x))*dt
yn = y + (x*(rho - z) - y)*dt
zn = z + (x*y - beta*z)*dt
return xn, yn, zn
def plot_3d(x, y, z):
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.plot(x, y, z)
plt.show()
dt = 0.001
n = 50000
x, y, z = np.empty(n), np.empty(n), np.empty(n)
x[0], y[0], z[0] = (0, 1, 1.05)
for i in range(n-1):
x[i+1], y[i+1], z[i+1] = lorenz_eq(x[i], y[i], z[i], dt)
plot_3d(x, y, z)
描画できました！
おわりに
ローレンツ方程式を数値的に解き、描画しました。
どういう結果か分かっているとはいえ、やっぱり実際に手を動かしてみると、簡単な式と実装で複雑な振る舞いが得られると面白いですよね。
他にもロジスティック写像など簡単な方程式でカオス的な挙動を示すものがあるので、それも実装したいと思います。
読んでいただきありがとうございました。
参考
https://disassemble-channel.com/python-lorentz-equation/
http://www.den.t.u-tokyo.ac.jp/ad_prog/ode/
https://qiita.com/shokishimada/items/acdb1e736d8f441d27d1