C++ でローレンツ・アトラクタの計算をしてみました。
今回は、微分方程式の近似解法に Runge-Kutta(ルンゲ=クッタ)法を使用します。

前回は、微分方程式の近似解法に Euler(オイラー)法を使用しています。

0. 前提条件

  • Debian GNU/Linux 10.6 (64bit) での作業を想定。
  • GCC 9.2.0 (G++ 9.2.0) (C++17) でのコンパイルを想定。

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

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

2. Runge-Kutta(ルンゲ=クッタ)法とは

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

3. C++ ソースコードの作成

  • パラメータ p, r, b の値はコマンドライン引数で指定する。
  • 計算部分と実行部分とでファイルを分けている。
calc.hpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
#ifndef LORENZ_ATTRACTOR_EULER_CALC_HPP_
#define LORENZ_ATTRACTOR_EULER_CALC_HPP_

#include <vector>

namespace my_lib {

class Calc {
  double p;
  double r;
  double b;
  bool lorenz(const double[], double(&)[3]);  // 計算(各ステップ)

public:
  Calc(double p, double r, double b) : p(p), r(r), b(b) {}     // コンストラクタ
  bool lorenz_runge_kutta(std::vector<std::vector<double>>&);  // 計算
};

}  // namespace my_lib

#endif
calc.cpp
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
#include "calc.hpp"

#include <cmath>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <vector>

namespace my_lib {

// 定数
constexpr double kDt   = 1.0e-3;  // Differential interval
constexpr int    kStep = 100000;  // Time step count
constexpr double kX0   = 1.0;     // Initial value of x
constexpr double kY0   = 1.0;     // Initial value of y
constexpr double kZ0   = 1.0;     // Initial value of z

/**
 * @brief      計算(ローレンツ・アトラクタ(Runge-Kutta 法)
 *
 * @param[ref] データ配列(計算結果) rec: (vector<vector<double>>)
 * @return     真偽(true|false)(bool)
 */
bool Calc::lorenz_runge_kutta(std::vector<std::vector<double>>& res) {
  double xyz[] = {kX0, kY0, kZ0};
  double xyz_l_0[3];  // 計算用(LorenzAttractor)
  double xyz_l_1[3];  // 計算用(LorenzAttractor)
  double xyz_l_2[3];  // 計算用(LorenzAttractor)
  double xyz_l_3[3];  // 計算用(LorenzAttractor)
  double xyz_w[3];    // 計算用(作業用)
  unsigned int i;     // ループインデックス
  unsigned int j;     // ループインデックス

  try {
    for (i = 0; i < kStep; i++) {
      if (!lorenz(xyz, xyz_l_0)) return false;
      xyz_w[0] = xyz[0] + xyz_l_0[0] * kDt / 2.0;
      xyz_w[1] = xyz[1] + xyz_l_0[1] * kDt / 2.0;
      xyz_w[2] = xyz[2] + xyz_l_0[2] * kDt / 2.0;
      if (!lorenz(xyz_w, xyz_l_1)) return false;
      xyz_w[0] = xyz[0] + xyz_l_1[0] * kDt / 2.0;
      xyz_w[1] = xyz[1] + xyz_l_1[1] * kDt / 2.0;
      xyz_w[2] = xyz[2] + xyz_l_1[2] * kDt / 2.0;
      if (!lorenz(xyz_w, xyz_l_2)) return false;
      xyz_w[0] = xyz[0] + xyz_l_2[0] * kDt;
      xyz_w[1] = xyz[1] + xyz_l_2[1] * kDt;
      xyz_w[2] = xyz[2] + xyz_l_2[2] * kDt;
      if (!lorenz(xyz_w, xyz_l_3)) return false;
      for (j = 0; j < 3; ++j) {
        xyz[j] += (xyz_l_0[j] + 2 * xyz_l_1[j] + 2 * xyz_l_2[j] + xyz_l_3[j])
                * kDt / 6.0;
      }
      res.push_back({xyz[0], xyz[1], xyz[2]});
    }
  } catch (...) {
    return false;  // 計算失敗
  }

  return true;  // 計算成功
}

bool Calc::lorenz(const double xyz[], double(&xyz_l)[3]) {
  try {
    xyz_l[0] = -p * xyz[0] + p * xyz[1];
    xyz_l[1] = -xyz[0] * xyz[2] + r * xyz[0] - xyz[1];
    xyz_l[2] = xyz[0] * xyz[1] - b * xyz[2];
  } catch (...) {
    return false;
  }

  return true;
}

}  // namespace my_lib
lorenz_attractor_runge_kutta.cpp
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
/***********************************************************
  Lorenz attractor (Runge-Kutta method)

    DATE          AUTHOR          VERSION
    2020.10.14    mk-mode.com     1.00 新規作成

  Copyright(C) 2020 mk-mode.com All Rights Reserved.
***********************************************************/
#include "calc.hpp"

#include <cstdlib>   // for EXIT_XXXX
#include <iomanip>   // for setprecision
#include <iostream>
#include <string>
#include <vector>

int main(int argc, char* argv[]) {
  double p;
  double r;
  double b;
  std::vector<std::vector<double>> res;  // データ配列(計算結果)
  std::size_t i;                         // loop インデックス

  try {
    // コマンドライン引数のチェック
    if (argc < 4) {
      std::cerr << "[ERROR] Number of arguments is wrong!\n"
                << "[USAGE] ./lorenz_attractor_runge_kutta p r b"
                << std::endl;
      return EXIT_FAILURE;
    }

    // p, r, b の取得
    p = std::stod(argv[1]);
    r = std::stod(argv[2]);
    b = std::stod(argv[3]);

    // 計算用オプジェクトのインスタンス化
    my_lib::Calc calc(p, r, b);

    // 計算
    if (!calc.lorenz_runge_kutta(res)) {
      std::cout << "[ERROR] Failed to calculare!" << std::endl;
      return EXIT_FAILURE;
    }

    // 結果出力
    std::cout << std::fixed << std::setprecision(8);
    for (i = 0; i < res.size(); ++i) {
      std::cout << std::setw(14) << std::right << res[i][0]
                << std::setw(14) << std::right << res[i][1]
                << std::setw(14) << std::right << res[i][2]
                << std::endl;
    }
  } catch (...) {
      std::cerr << "EXCEPTION!" << std::endl;
      return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;
}

4. ソースコードのコンパイル

まず、以下のように Makefile を作成する。(行頭のインデントはタブ文字

Makefile
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
gcc_options = -std=c++17 -Wall -O2 --pedantic-errors

lorenz_attractor_euler: lorenz_attractor_runge_kutta.o calc.o
  g++ $(gcc_options) -o $@ $^

lorenz_attractor_euler.o : lorenz_attractor_runge_kutta.cpp
  g++ $(gcc_options) -c $<

calc.o : calc.cpp
  g++ $(gcc_options) -c $<

run : lorenz_attractor_runge_kutta
  ./lorenz_attractor_runge_kutta

clean :
  rm -f ./lorenz_attractor_runge_kutta
  rm -f ./*.o

.PHONY : run clean

そして、ビルド(コンパイル&リンク)。

1
$ make

5. 動作確認

コマンドライン引数に p, r, b の値を指定して実行する。
計算結果が出力される。

1
2
3
4
5
6
7
8
9
$ ./lorenz_attractor_runge_kutta 10 28 2.66666667

                      :

   14.12381298   18.68378846   29.33978919
   14.16899368   18.64485725   29.52533963
   14.21332423   18.60327057   29.71065962
   14.25678662   18.55901971   29.89569240
   14.29936294   18.51209735   30.08038071

計算結果をファイルに出力したければ、以下のようにする。

1
$ ./lorenz_attractor_runge_kutta 10 28 2.66666667 > data.txt

6. 結果確認

参考までに、出力された計算結果を GNUplot で描画してみた。

LORENZ_ATTRACTOR_RUNGE_KUTTA


以上。