C++ source code to calculate a Lorenz attractor by Runge-Kutta's method.

calc.cpp

#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

calc.hpp

#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

lorenz_attractor_runge_kutta.cpp

/***********************************************************
  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;
}