Python - ローレンツ・アトラクタ（Runge-Kutta 法）！

[ プログラミング, 数学 ] [ Ruby ]

先日、 Ruby でローレンツ・アトラクタを計算＆描画しました。

そして、前回、 Python で微分方程式の近似解法に Euler（オイラー）法を使用して、ローレンツ・アトラクタを計算＆描画しました。

今回は、微分方程式の近似解法に Runge-Kutta（ルンゲ＝クッタ）法を使用して、計算＆描画してみました。（Python で）

0. 前提条件

  • Python 3.6.4 （エンコード：UTF-8）での作業を想定。

1. ローレンツ方程式／アトラクタとは

  • 「ローレンツ方程式」とは、気象学者「エドワード・Ｎ・ローレンツ(Edward N. Lorenz)」が作成した力学系方程式をより単純化した、次のような非線形微分方程式。
    パラメータ p, r, b をほんの少し変えるだけで、これらの方程式から得られる軌跡は大きく異なったものになる。
  • 「ローレンツ方程式」は、カオス理論を学習する際に序盤で登場する方程式で、カオス研究の先駆的なもの。
  • 「アトラクタ」とは、ある力学系がそこに向かって時間発展する集合のことで、カオス理論における研究課題の一つ。
  • 「ローレンツ・アトラクタ」とは、ストレンジ・アトラクタの一種。
  • 「ローレンツ・アトラクタ」は、言い換えれば、「ローレンツ方程式のカオスのストレンジ・アトラクタ」である。

2. Runge-Kutta（ルンゲ＝クッタ）法とは

  • Euler 法よりは計算に時間がかかるが、その分、精度も高い。
  • 実際の研究等では、 Euler 法ではなく Runge-Kutta 法が使用されることが多い。
  • ここでは、 Runge-Kutta 法の詳細については説明しない。

3. Python スクリプトの作成

  • 数値演算ライブラリ NumPy は使用しない。（使用するほどでもないので）
  • グラフ描画には “matplotlib” ライブラリを使用。（要インストール）
  • パラメータ p, r, b のデフォルト値は 10, 28, 8/3 としている。
    変更したければ、 lorenz メソッド呼び出し時に引数で指定する。
  • 今回使用するのは4次の Runge-Kutta 法。
  • 実際、ローレンツ・アトラクタを計算する際に計算の粗い Euler 法を使用することはあまりなく、 Runge-Kutta 法で計算することが多い。
lorenz_attractor_runge_kutta.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82

#! /usr/local/bin/python3.6
"""
Lorenz attractor (Runge-Kutta method)
"""
import sys
import traceback
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


class LorenzAttractorRungeKutta:
    DT            = 1e-3     # Differential interval
    STEP          = 100000   # Time step count
    X_0, Y_0, Z_0 = 1, 1, 1  # Initial values of x, y, z

    def __init__(self):
        self.res = [[], [], []]

    def exec(self):
        """ Loranz attractor (Runge-Kutta method) execution """
        try:
            xyz = [self.X_0, self.Y_0, self.Z_0]
            for _ in range(self.STEP):
                k_0 = self.__lorenz(xyz)
                k_1 = self.__lorenz([
                    x + k * self.DT / 2 for x, k in zip(xyz, k_0)
                ])
                k_2 = self.__lorenz([
                    x + k * self.DT / 2 for x, k in zip(xyz, k_1)
                ])
                k_3 = self.__lorenz([
                    x + k * self.DT for x, k in zip(xyz, k_2)
                ])
                for i in range(3):
                    xyz[i] += (k_0[i] + 2 * k_1[i] + 2 * k_2[i] + k_3[i]) \
                            * self.DT / 6.0
                    self.res[i].append(xyz[i])
            self.__plot()
        except Exception as e:
            raise

    def __lorenz(self, xyz, p=10, r=28, b=8/3.0):
        """ Lorenz equation

        :param  list xyz
        :param  float  p
        :param  float  r
        :param  float  b
        :return list xyz
        """
        try:
            return [
                -p * xyz[0] + p * xyz[1],
                -xyz[0] * xyz[2] + r * xyz[0] - xyz[1],
                xyz[0] * xyz[1] - b * xyz[2]
            ]
        except Exception as e:
            raise

    def __plot(self):
        """ Protting """
        try:
            fig = plt.figure()
            ax = Axes3D(fig)
            ax.set_xlabel("x")
            ax.set_ylabel("y")
            ax.set_zlabel("z")
            ax.set_title("Lorenz attractor (Runge-Kutta method)")
            ax.plot(self.res[0], self.res[1], self.res[2], color="red", lw=1)
            #plt.show()
            plt.savefig("lorenz_attractor_runge_kutta.png")
        except Exception as e:
            raise


if __name__ == '__main__':
    try:
        obj = LorenzAttractorRungeKutta()
        obj.exec()
    except Exception as e:
        traceback.print_exc()
        sys.exit(1)

4. Ruby スクリプトの実行

まず、実行権限を付与。

1

$ chmod +x lorenz_attractor_runge_kutta.py

そして、実行。

1

$ ./lorenz_attractor_runge_kutta.py

5. 結果確認

Python スクリプトと同じディレクトリ内に “lorenz_attractor_runge_kutta.png” という画像ファイルが作成されるので、確認してみる。

LORENZ_ATTRACTOR_RUNGE_KUTTA

以上。