常微分方程式の初期問題の数値計算をするC++コード

ねこさんは昔から常微分方程式の数値計算がそれはそれは大好きで、学生時代もそれにのめり込んで気が付いたらプログラムばっかり書いてる人生になってしまったのですが、「常微分方程式を計算するプログラムを綺麗に書きたい！」って言ってるうちに泥沼にはまり込んでしまったともいえる訳で orz

んで最近になってboost::fusionっていうステキライブラリの存在を知りまして、ろくにtemplate metaprogrammingなんて解らないくせに、「これで常微分方程式のスキームを書いたらきっとステキだー」と思ってしまって、試してみたら直ぐに蹴躓いて、恥ずかしながら某スレに質問して戴いた答えを元にどうにか書いてみました。

#include <iostream> #include <boost/fusion/tuple.hpp> #include <boost/fusion/algorithm.hpp> namespace ode { namespace euler { template <typename FLT=double> struct scheme {     FLT dt;     template<typename T>     T operator()(const T& x0, const T& x1) const {         return x0+x1*dt;     }     template<typename T>     struct result;     template<typename T>     struct result<scheme<FLT>(T,T)> {         typedef T type;     };     scheme(const FLT& dt=1.0) : dt(dt) {} }; } template <typename FLT, template <typename> class Sc, class Op, class X> void solve(Sc<FLT>& prop, const Op& op, X& x) {     X x0 = x;     op(x);     x = boost::fusion::transform(x0,x,prop); } } typedef boost::fusion::tuple<double,double> state; std::ostream& operator<<(std::ostream& os, const state& x_) {     double t,x;     boost::fusion::tie(t,x) = x_;     return os << t << " " << x; } struct Logistic {     double a,b;     Logistic(const double& a, const double& b) : a(a), b(b) {}     template <class X>     void operator()(X& x_) const {         double t,x;         boost::fusion::tie(t,x) = x_;         t = 1.0;         x = a*x - b*x*x;         x_ = X(t,x);     } }; int main() {     state x(0.0,0.01);     ode::euler::scheme<double> euler(0.01);     Logistic logi(1.0,1.0);     std::cout << std::scientific;     for (int i=0; i<1000; ++i) {         ode::solve(euler,logi,x);         std::cout << x << std::endl;     }     return 0; }

ベタに書くよりはオーバーヘッドがあるだろうし、あんまり速くないかも(計測しろよ)しれないけど、簡単にコードが書けるようになりそうなのが良いかな。とりあえず簡単にEuler法で書いてるけど、多分インターフェースそのままにRunge-Kuttaでも何でも拡張できると思うです…やってみてないから確証はないですけど orz。時間ができたらそのうち書いてみるです。