Euler法の結果をgnuplotでスッと出したかった

#目的
数年前に講義で習った常微分方程式の数値解法.
Euler法は簡単に理解できるし, 書くのも簡単で, 精度が求められない場合にはそこそこ便利.
当時はC言語で配列に結果を格納して, csvで書き出して, EXCELで開いてグラフにして…みたいなことをしていた.
ちょっとプログラムを見返していて, このなんとも面倒な部分をスマートに書けないものかと, ちょっとだけ書き直してみた.

#準備するもの
gnuplotをインストールしてください.
今回はver. 5.2.7で試してあります.
あとはC++のコンパイラと標準ライブラリだけでできる…はず.

#コード

EulerCalculator.hpp

#pragma once
#include <vector>

class alignas(8) EulerCalculator final{
/*定数定義*/
private:
	constexpr static double_t	DEFAULT_X_MAX_		= 10.0;
	constexpr static size_t		DEFAULT_LOOP_NUM_	= 1000;
	constexpr static double_t	DEFAULT_Y_INIT_		= 1.0;
	
/*毎回コンストラクタの作り方に頭を抱える*/
public:
	EulerCalculator() noexcept{
		init(DEFAULT_Y_INIT_, DEFAULT_X_MAX_, DEFAULT_LOOP_NUM_);
	}
	EulerCalculator(const double_t y_init, const double_t x_maximum = DEFAULT_X_MAX_, const size_t loop_num = DEFAULT_LOOP_NUM_) noexcept(false){
		init(y_init, x_maximum, loop_num);
	}

	EulerCalculator(const EulerCalculator&)						= delete;
	EulerCalculator(const EulerCalculator&&)					= delete;
	const EulerCalculator& operator=(const EulerCalculator&)	= delete;
	const EulerCalculator& operator=(const EulerCalculator&&)	= delete;
	~EulerCalculator() noexcept									= default;

/*ここからがメンバ*/
private:
	double_t				step_	= 0.0;
	double_t				y_init_ = DEFAULT_Y_INIT_;

	/*この辺はグラフ描画で使用*/
	/*メンバ用意しなくてもよかったかな*/
	double_t				x_min_ 	= 0.0;
	double_t				x_max_ 	= DEFAULT_X_MAX_;
	double_t				y_min_ 	= 0.0;
	double_t				y_max_ 	= DEFAULT_Y_INIT_;

	/*計算結果格納用コンテナ*/
	std::vector<double_t>	x_;
	std::vector<double_t>	y_;

/*以下メソッド*/
public:
	const double_t getMinX() const noexcept{
		return(x_min_);
	}
	const double_t getMaxX() const noexcept{
		return(x_max_);
	}
	const double_t getMinY() const noexcept{
		return(y_min_);
	}
	const double_t getMaxY() const noexcept{
		return(y_max_);
	}

	std::void_t<>	calc() 				noexcept(false);
	bool			writeDat() const 	noexcept(false);

private:
	std::void_t<>	init(const double_t y_init_, const double_t x_maximum, const size_t loop_num) noexcept(false);
	/*関数の名前の付け方が下手*/
	const double_t	dy_dxFunc(const double_t x, const double_t y) const noexcept(false);
};

EulerCulculator.cpp

#include "EulerCalculator.hpp"
#include <algorithm>
#include <fstream>

using namespace std;

void_t<> EulerCalculator::calc() noexcept(false){
	for(size_t i(0), end = y_.size() - 1; i < end; ++i){
		x_[i + 1] = x_[i] + step_;
		/*ここがEuler法*/
		y_[i + 1] = y_[i] + step_ * dy_dxFunc(x_[i], y_[i]);
	}

	/*yの最大・最小を更新*/
	y_min_ = *min_element(y_.begin(), y_.end());
	y_max_ = *max_element(y_.begin(), y_.end());

	return;
}

bool EulerCalculator::writeDat() const noexcept(false){
	ofstream dat_file("Output.dat");
	if(!dat_file.is_open()){
		return(false);
	}

	for(size_t i(0), end(y_.size()); i < end; ++i){
		dat_file << x_[i] << " " << y_[i] << endl;
	}

	dat_file.close();

	return(true);
}

void_t<> EulerCalculator::init(const double_t y_init_, const double_t x_maximum, const size_t loop_num) noexcept(false){
	step_ = x_maximum / static_cast<double_t>(loop_num);
	x_max_ = x_maximum;

	x_.reserve(loop_num + 1);
	x_.resize(loop_num + 1, 0.0);

	y_.reserve(loop_num + 1);
	y_.resize(loop_num + 1, y_init_);

	return;
}

const double_t EulerCalculator::dy_dxFunc(const double_t x, const double_t y) const noexcept(false){
	return(-y);
}

Main.cpp

include "EulerCalculator.hpp"

using namespace std;

int main(int argc, char** argv){
	EulerCalculator calculator;
	calculator.calc();
	if(!calculator.writeDat()){
		return(EXIT_FAILURE);
	}

	/*グラフの描画範囲を取得*/
	double_t x_min = calculator.getMinX();
	double_t x_max = calculator.getMaxX();
	double_t y_min = calculator.getMinY();
	double_t y_max = calculator.getMaxY();

#ifdef _WIN32
	/*ファイルの置き場所が置き場所なのでPathを""で囲う必要があった*/
	const char* GNUPLOT_PATH = "\"C:/Program Files/gnuplot/bin/gnuplot.exe\"";

	/*プロセスオープン*/
	FILE* gnuplot = _popen(GNUPLOT_PATH, "w");
	if(gnuplot == nullptr){
		return(EXIT_FAILURE);
	}

	/*後はgnuploのコマンドをfprintfで入力していく*/
	fprintf(gnuplot, "set size square\n");
	fprintf(gnuplot, "set xrange [%lf:%lf]\n", x_min, x_max);
	fprintf(gnuplot, "set yrange [%lf:%lf]\n", y_min, y_max);
	fprintf(gnuplot, "unset key\n");
	fprintf(gnuplot, "set title 'Euler Method'\n");

	fprintf(gnuplot, "plot 'Output.dat'w p pt 7 ps 0.5\n");
	/*結果を吐き出させる*/
	fflush(gnuplot);

	/*すぐ消えてしまうので待ってもらう*/
	system("pause");

	/*FIXME: たまに文字列がeとxitで切れてしまっている?*/
	fprintf(gnuplot, "exit\n");
	fflush(gnuplot);
	_pclose(gnuplot);
#endif

	return(EXIT_SUCCESS);
}

#結果
上のプログラム通りならきっと結果は下みたいな感じ.

見返せば, クラス作る必要あったかな…と思ったり.
コード自体に難しいことはない…はず.
簡単なことしかやってないからね.

#一部を変えてみる
コードを変えるならこんな風に.
以下, 変更部分のみ抜粋.

EulerCalculator.cpp

const double_t EulerCalculator::dy_dxFunc(const double_t x, const double_t y) const noexcept(false){
	return((x * x + y * y) / (2.0 * x * y + 1.0));
}

Main.cpp

int main(int argc, char** argv){
	EulerCalculator calculator(5.0, 15.0, 5000);
	calculator.calc();
	if(!calculator.writeDat()){
		return(EXIT_FAILURE);
	}

	/*以下略*/

すると結果は以下のようになる.

#雑記
なんとなくお分かりかもしれませんが, $y$の初期値は$x=0$のときのものです.
$x$の初期値だったり範囲だったりもきちんと入力できた方がいいな, と一通り書き終わって気が付きました.
$\dfrac{\mathrm{d}y}{\mathrm{d}x}=\ln{x}$
のようなものは, このままじゃ解けないでエラー吐くんだろうなと.

どうしても_popen使用で付随するFILE*やらfprintfやらをモダンに書きたいのであるが…ううむ.

今度はJavaScriptかPython辺りで書こうか.