Skip to content

Instantly share code, notes, and snippets.

@komasaru

komasaru/calc.cpp

Created Oct 15, 2020
Embed
What would you like to do?
C++ source code to calculate a Lorenz attractor by Euler's method.
#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 計算(ローレンツ・アトラクタ(Euler 法)
*
* @param[ref] データ配列(計算結果) rec: (vector<vector<double>>)
* @return 真偽(true|false)(bool)
*/
bool Calc::lorenz_euler(std::vector<std::vector<double>>& res) {
double xyz[] = {kX0, kY0, kZ0};
double xyz_l[3]; // 計算用(LorenzAttractor)
unsigned int i; // ループインデックス
unsigned int j; // ループインデックス
try {
for (i = 0; i < kStep; ++i) {
if (!lorenz(xyz, xyz_l)) return false;
for (j = 0; j < 3; ++j) xyz[j] += kDt * xyz_l[j];
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
#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_euler(std::vector<std::vector<double>>&); // 計算
};
} // namespace my_lib
#endif
/***********************************************************
Lorenz attractor (Euler method)
DATE AUTHOR VERSION
2020.09.18 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_euler 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_euler(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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment