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

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

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

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

0. 前提条件

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

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

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

2. Euler（オイラー）法とは

  • 微分方程式の近似解法の中で計算が比較的簡単なものだが、その分、計算も粗い。
  • 実際の研究等で使用されることはほとんどない。
  • 近似解法の概念を理解するための一助にはなる。
  • ここでは、 Euler 法の詳細については説明しない。

3. Python スクリプト作成

  • 数値演算ライブラリ NumPy は使用しない。（使用するほどでもないので）
  • グラフ描画には “matplotlib” ライブラリを使用。（要インストール）
  • パラメータ p, r, b のデフォルト値は 10, 28, 8/3 としている。
    変更したければ、 lorenz メソッド呼び出し時に引数で指定する。
  • 実際、ローレンツ・アトラクタを計算する際に計算の粗い Euler 法を使用することはあまりない（ルンゲ＝クッタ法で計算することが多い）ので、あくまでも、参考程度に。
lorenz_attractor_euler.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

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


class LorenzAttractorEuler:
    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 (Euler method) execution """
        try:
            xyz = [self.X_0, self.Y_0, self.Z_0]
            for _ in range(self.STEP):
                l = self.__lorenz(xyz)
                for i in range(3):
                    xyz[i] += self.DT * l[i]
                    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 (Euler method)")
            ax.plot(self.res[0], self.res[1], self.res[2], lw=1)
            #plt.show()
            plt.savefig("lorenz_attractor_euler.png")
        except Exception as e:
            raise


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

4. Python スクリプトの実行

まず、実行権限を付与。

1

$ chmod +x lorenz_attractor_euler.py

そして、実行。

1

$ ./lorenz_attractor_euler.py

5. 結果確認

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

LORENZ_ATTRACTOR_EULER

以上。