[Python] ローレンツアトラクタの曲線をmatplotlibで可視化する方法

ローレンツアトラクタは、カオス理論で有名な3次元の非線形力学系です。

Pythonでローレンツアトラクタを可視化するには、SciPyを使って微分方程式を解き、matplotlibでプロットします。

まず、SciPyのodeint関数を使ってローレンツ方程式を数値的に解きます。

次に、matplotlibのAxes3Dを使って3次元プロットを作成します。

ローレンツ方程式は次の形です：

$$\frac{dx}{dt}=\sigma (y–x),\frac{dy}{dt}=x(\rho –z)–y,\frac{dz}{dt}=xy–\beta z$$

ローレンツアトラクタとは

ローレンツアトラクタは、気象学者エドワード・ローレンツによって提唱されたカオス理論の一例です。

これは、非線形の常微分方程式から導かれる三次元の動的システムであり、初期条件に対する敏感な依存性を示します。

ローレンツアトラクタは、特に気象予測におけるカオス的な挙動を理解するための重要なモデルとされています。

アトラクタの形状は、複雑で美しい曲線を描き、カオス理論の視覚的な象徴として広く知られています。

このアトラクタは、システムの状態が時間とともにどのように変化するかを示すもので、カオス的な挙動を持つ多くの自然現象に関連しています。

必要なライブラリのインストール

ローレンツアトラクタを可視化するためには、いくつかのPythonライブラリが必要です。

以下に、各ライブラリのインストール方法を示します。

SciPyのインストール方法

SciPyは、科学技術計算のためのライブラリで、数値解法に使用します。

以下のコマンドを実行してインストールします。

pip install scipy

matplotlibのインストール方法

matplotlibは、データの可視化を行うためのライブラリです。

以下のコマンドを実行してインストールします。

pip install matplotlib

NumPyのインストール方法

NumPyは、数値計算を効率的に行うためのライブラリで、配列操作に使用します。

以下のコマンドを実行してインストールします。

pip install numpy

これらのライブラリをインストールすることで、ローレンツアトラクタの計算と可視化が可能になります。

ローレンツ方程式の実装

ローレンツアトラクタを計算するためには、ローレンツ方程式を定義し、初期条件やパラメータを設定する必要があります。

以下にその手順を示します。

ローレンツ方程式の定義

ローレンツ方程式は、以下の3つの常微分方程式から構成されています。

$$\begin{array}{rl}\frac{dx}{dt}& =\sigma (y–x)\\ \frac{dy}{dt}& =x(\rho –z)–y\\ \frac{dz}{dt}& =xy–\beta z\end{array}$$

ここで、$\sigma$、$\rho$、$\beta$はシステムのパラメータです。

初期条件とパラメータの設定

次に、ローレンツ方程式を解くための初期条件とパラメータを設定します。

以下のように設定します。

# パラメータの設定
sigma = 10.0
rho = 28.0
beta = 8/3
# 初期条件
initial_conditions = [1.0, 1.0, 1.0]  # [x0, y0, z0]

SciPyのodeintを使った数値解法

SciPyのodeint関数を使用して、ローレンツ方程式を数値的に解きます。

以下のコードでは、ローレンツ方程式を定義し、時間の範囲を設定して解を求めます。

import numpy as np
from scipy.integrate import odeint
# ローレンツ方程式の定義
def lorenz_system(state, t):
    x, y, z = state
    dxdt = sigma * (y - x)
    dydt = x * (rho - z) - y
    dzdt = x * y - beta * z
    return [dxdt, dydt, dzdt]
# 時間の範囲
t = np.linspace(0, 50, 10000)  # 0から50までの時間を10000分割
# 数値解法
solution = odeint(lorenz_system, initial_conditions, t)

このコードを実行すると、ローレンツ方程式の解がsolutionに格納されます。

これにより、ローレンツアトラクタの軌跡を計算する準備が整いました。

matplotlibを使ったローレンツアトラクタの可視化

ローレンツアトラクタを可視化するためには、matplotlibを使用して3Dプロットを作成します。

以下にその手順を示します。

matplotlibの3Dプロットの基本

matplotlibでは、3Dプロットを作成するためにmpl_toolkits.mplot3dモジュールを使用します。

まず、必要なライブラリをインポートし、3Dプロットの準備を行います。

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

Axes3Dを使った3次元プロットの作成

次に、Axes3Dを使用して3Dプロットを作成します。

以下のコードでは、ローレンツアトラクタの軌跡を描くための3Dプロットを設定します。

# 3Dプロットの作成
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# 軌跡のデータを取得
x = solution[:, 0]
y = solution[:, 1]
z = solution[:, 2]
# 3Dプロットにデータをプロット
ax.plot(x, y, z, lw=0.5)

ローレンツアトラクタの軌跡を描画する方法

上記のコードを実行すると、ローレンツアトラクタの軌跡が3D空間に描画されます。

プロットを表示するためには、plt.show()を呼び出します。

# プロットの表示
plt.title("Lorenz Attractor")
plt.xlabel("X-axis")
plt.ylabel("Y-axis")
ax.set_zlabel("Z-axis")
plt.show()

軌跡の色やスタイルのカスタマイズ

プロットの色やスタイルをカスタマイズすることも可能です。

例えば、軌跡の色を変更したり、線のスタイルを変更することができます。

以下のコードでは、色を青に設定し、線のスタイルを点線に変更しています。

# 軌跡の色とスタイルのカスタマイズ
ax.plot(x, y, z, color='blue', linestyle='--', lw=0.5)

これにより、ローレンツアトラクタの可視化がより魅力的になります。

プロットのカスタマイズを行うことで、データの視覚的な表現を向上させることができます。

ローレンツアトラクタのアニメーション化

ローレンツアトラクタをアニメーション化することで、時間の経過に伴う軌跡の変化を視覚的に表現できます。

以下に、matplotlibのFuncAnimationを使用したアニメーションの作成方法を示します。

matplotlibのFuncAnimationを使ったアニメーションの作成

FuncAnimationを使用して、ローレンツアトラクタのアニメーションを作成します。

以下のコードでは、アニメーションのフレームを更新する関数を定義し、アニメーションを生成します。

from matplotlib.animation import FuncAnimation
# アニメーションの更新関数
def update(frame):
    ax.cla()  # 現在のプロットをクリア
    ax.set_title("Lorenz Attractor")
    ax.set_xlabel("X-axis")
    ax.set_ylabel("Y-axis")
    ax.set_zlabel("Z-axis")
    ax.plot(x[:frame], y[:frame], z[:frame], color='blue', lw=0.5)
# アニメーションの作成
ani = FuncAnimation(fig, update, frames=len(t), interval=1)

アニメーションの保存方法

作成したアニメーションをファイルとして保存することも可能です。

以下のコードでは、アニメーションをMP4形式で保存します。

# アニメーションの保存
ani.save('lorenz_attractor.mp4', writer='ffmpeg', fps=30)

このコードを実行すると、lorenz_attractor.mp4というファイル名でアニメーションが保存されます。

アニメーションのカスタマイズ

アニメーションの見た目をカスタマイズすることもできます。

例えば、軌跡の色や線のスタイルを変更したり、アニメーションの速度を調整することが可能です。

以下のコードでは、軌跡の色を赤に変更し、アニメーションの速度を遅くしています。

# アニメーションの更新関数のカスタマイズ
def update(frame):
    ax.cla()  # 現在のプロットをクリア
    ax.set_title("Lorenz Attractor")
    ax.set_xlabel("X-axis")
    ax.set_ylabel("Y-axis")
    ax.set_zlabel("Z-axis")
    ax.plot(x[:frame], y[:frame], z[:frame], color='red', lw=0.5)
# アニメーションの作成
ani = FuncAnimation(fig, update, frames=len(t), interval=50)  # intervalを50msに設定

これにより、アニメーションの色や速度を調整し、視覚的に魅力的な表現を実現できます。

アニメーション化することで、ローレンツアトラクタの動的な特性をより効果的に伝えることができます。

応用例：ローレンツアトラクタのパラメータ変更

ローレンツアトラクタは、パラメータや初期条件の変更によってその挙動が大きく変わります。

以下に、これらの変更がどのように影響するかを示します。

パラメータ$\sigma$、$\rho$、$\beta$の変更による影響

ローレンツ方程式のパラメータである$\sigma$、$\rho$、$\beta$は、アトラクタの形状や動的特性に大きな影響を与えます。

例えば：

• $\sigma$(プルーム数)は、流体の対流の強さを示します。

大きくすると、アトラクタの複雑さが増し、よりカオス的な挙動を示します。

• $\rho$(レイリー数)は、流体の温度差を示し、これを変更すると、アトラクタの形状が変わります。

特に、$\rho$が特定の値(約24.74)を超えると、アトラクタが複雑な構造を持つようになります。

• $\beta$は、流体の粘性を示し、これを変更することで、アトラクタの収束速度や形状が変わります。

これらのパラメータを変更することで、異なるローレンツアトラクタを生成し、その挙動を観察することができます。

初期条件の変更による軌跡の変化

初期条件もローレンツアトラクタの挙動に大きな影響を与えます。

例えば、初期条件を以下のように変更することができます。

# 新しい初期条件
initial_conditions = [0.0, 1.0, 1.05]  # [x0, y0, z0]

このように初期条件を変更すると、アトラクタの軌跡が異なるパターンを描くことになります。

特に、カオス的なシステムでは、初期条件のわずかな違いが最終的な軌跡に大きな影響を与えることがあります。

他のカオス系との比較

ローレンツアトラクタは、他のカオス系と比較することで、その特性をより深く理解することができます。

例えば、ロスラーアトラクタやヘノン写像などの他のカオス的なシステムと比較することができます。

• ロスラーアトラクタは、以下の方程式で定義され、ローレンツアトラクタとは異なる形状を持ちます。

$$\begin{array}{rl}\frac{dx}{dt}& =-y–z\\ \frac{dy}{dt}& =x+ay\\ \frac{dz}{dt}& =b+z(x–c)\end{array}$$

ここで、$a$、$b$、$c$はパラメータです。

• ヘノン写像は、離散的なカオス系であり、以下のように定義されます。

$$\begin{array}{rl}{x}_{n+1}& =1–a{x}_n^2+y_n\\ y_{n+1}& =b{x}_n\end{array}$$

ここで、$a$と$b$はパラメータです。

これらのカオス系とローレンツアトラクタを比較することで、カオス理論の理解が深まり、異なるシステムの挙動を観察することができます。

応用例：他のカオスアトラクタの可視化

ローレンツアトラクタ以外にも、さまざまなカオスアトラクタが存在します。

ここでは、ロスラーアトラクタ、ヘノン写像、チュア回路の可視化方法を紹介します。

ロスラーアトラクタの可視化

ロスラーアトラクタは、以下の方程式で定義されます。

$$\begin{array}{rl}\frac{dx}{dt}& =-y–z\\ \frac{dy}{dt}& =x+ay\\ \frac{dz}{dt}& =b+z(x–c)\end{array}$$

ここで、$a$、$b$、$c$はパラメータです。

以下のコードを使用して、ロスラーアトラクタを可視化します。

# ロスラーアトラクタのパラメータ
a = 0.2
b = 0.2
c = 5.7
# ロスラー方程式の定義
def rossler_system(state, t):
    x, y, z = state
    dxdt = -y - z
    dydt = x + a * y
    dzdt = b + z * (x - c)
    return [dxdt, dydt, dzdt]
# 初期条件と時間の設定
initial_conditions = [1.0, 1.0, 1.0]
t = np.linspace(0, 100, 10000)
# 数値解法
rossler_solution = odeint(rossler_system, initial_conditions, t)
# ロスラーアトラクタの可視化
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(rossler_solution[:, 0], rossler_solution[:, 1], rossler_solution[:, 2], color='green', lw=0.5)
plt.title("Rossler Attractor")
plt.show()

ヘノン写像の可視化

ヘノン写像は、離散的なカオス系であり、以下のように定義されます。

$$\begin{array}{rl}{x}_{n+1}& =1–a{x}_n^2+y_n\\ y_{n+1}& =b{x}_n\end{array}$$

ここで、$a$と$b$はパラメータです。

以下のコードを使用して、ヘノン写像を可視化します。

# ヘノン写像のパラメータ
a = 1.4
b = 0.3
# 初期条件
x, y = 0, 0
iterations = 10000
x_values = []
y_values = []
# ヘノン写像の計算
for _ in range(iterations):
    x_new = 1 - a * x**2 + y
    y_new = b * x
    x_values.append(x_new)
    y_values.append(y_new)
    x, y = x_new, y_new
# ヘノン写像の可視化
plt.scatter(x_values, y_values, s=0.1, color='purple')
plt.title("Henon Map")
plt.xlabel("X-axis")
plt.ylabel("Y-axis")
plt.show()

チュア回路の可視化

チュア回路は、非線形回路の一例で、以下の方程式で表されます。

$$\begin{array}{rl}\frac{dx}{dt}& =\alpha (y–x)\\ \frac{dy}{dt}& =x–y+z\\ \frac{dz}{dt}& =-\beta y\end{array}$$

ここで、$\alpha$と$\beta$はパラメータです。

以下のコードを使用して、チュア回路を可視化します。

# チュア回路のパラメータ
alpha = 0.2
beta = 0.2
# チュア回路の方程式の定義
def chua_system(state, t):
    x, y, z = state
    dxdt = alpha * (y - x)
    dydt = x - y + z
    dzdt = -beta * y
    return [dxdt, dydt, dzdt]
# 初期条件と時間の設定
initial_conditions = [0.0, 0.0, 0.0]
t = np.linspace(0, 100, 10000)
# 数値解法
chua_solution = odeint(chua_system, initial_conditions, t)
# チュア回路の可視化
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(chua_solution[:, 0], chua_solution[:, 1], chua_solution[:, 2], color='orange', lw=0.5)
plt.title("Chua's Circuit")
plt.show()

これらのコードを実行することで、ロスラーアトラクタ、ヘノン写像、チュア回路の可視化が可能になります。

各カオスアトラクタの特性を観察することで、カオス理論の理解が深まります。

まとめ

この記事では、ローレンツアトラクタの可視化方法や、他のカオスアトラクタの実装について詳しく解説しました。

また、パラメータや初期条件の変更がアトラクタの挙動に与える影響についても触れました。

これらの知識を活用して、さまざまなカオスシステムの挙動を観察し、さらなる探求を行ってみてください。

興味を持った方は、実際にコードを実行してみたり、他のカオス系の可視化に挑戦してみることをお勧めします。