数値計算の基本をC++で書いていきます.

boostやgnuplotを使用しているコードのコンパイルには, boostやgnuplotのインストールがそれぞれ必要です.

変更した方が良い箇所があれば, この記事かTwitterの@Yonono_01まで連絡ください.

目次

常微分方程式 (Ordinaly Differential Equation)

常微分方程式 (Ordinaly Differential Equation)

常微分方程式を数値的に解くとは, ある初期値から始めて, ある刻み幅で増加する時の値を順々に求めていくことをいう.
一階常微分方程式を数値的に解く方法として, Euler法やRunge-Kutta法などがある.

Euler法 (Euler Method)

Euler法は一階常微分方程式

\frac{dy}{dx} = f(x,y)

と初期条件

y(x_0) = y_0

について, $x_{i+1}=x_i+dh$に対応する$y_{i+1}$の値を

y_{i+1}=y_i+f(x,y)dh

と近似する.
Euler法は直線での近似であるので, 精度はあまり良くない.

例えば, 初期条件$v(t_0)=v_0, x(t_0)=x_0$で自由落下する物体の運動を, 刻み幅$dt$で増加する時刻$t_1, t_2, \cdots$に対して順に求める.
この物体の運動方程式は

\frac{d^2x}{dt^2} = g

のようになる.
この二階常微分方程式を次の連立一階常微分方程式

\frac{dv}{dt} = g \\
\frac{dx}{dt} = v

としてEuler法を用いる.

#include <iomanip>  // for std::ios::fixed, std::setiosflags, std::setprecision
#include <iostream> // for std::cout
#include <string>   // for std::string

namespace{
    static auto constexpr MAXBUFSIZE = 32;
    static auto constexpr g = 9.80665;

    template <typename T, bool check_positive_num = true>
    T input_parameter(std::string const & str){
        T param;
        while (true){
		    std::cout << str;
            std::cin >> param;

		    std::cin.clear();
		    std::cin.ignore(MAXBUFSIZE, '\n');

            // 初期条件として<=0の場合に再入力が要らない時は<type, false>
            if constexpr (check_positive_num){
                if (!std::cin.fail() && param > static_cast<T>(0)){ 
			        break;
		        }
            } 
            else{
                if (!std::cin.fail()){
                    break;
                }
            }
	    }

        return param;
    }

    // motion equation dv/dt
    double dvdt(double v);

    // motion equation dx/dt
    double dxdt(double v);

    void Euler(double v0, double x0, double dt);
}

int main(){
    auto const v0 = input_parameter<double, false> ("Enter the initial velocity v0[m/s]\n");
    auto const x0 = input_parameter<double> ("Enter the initial position x0[m]\n");
    auto const dt = input_parameter<double> ("Enter the partition width dt[s]\n");

    Euler(v0, x0, dt);
}

namespace{
    void Euler(double v0, double x0, double dt){
        // t, v, xの初期化
        double t = 0;
        auto v = v0;
        auto x = x0;

        std::cout << std::setprecision(7) << std::setiosflags(std::ios::fixed);

        // x<=0になるまで計算
        while (x > 0){
            std::cout << "when the time is " << t << "[s], the velocity is " << v << "[m/s] and the position x is " << x << "[m]" << std::endl;
            
            // Euler method
            t += dt;
            v += dvdt(v) * dt;
            x -= dxdt(v) * dt;
        }
    }

    double dvdt(double v){
        return g;
    }

    double dxdt(double v){
        return v;
    }
}

前進Euler法の誤差

前進Euler法による数値計算が$1$次精度であることをグラフを描いて確認する.
空気抵抗のある質点の運動方程式は適当な変数変換により無次元化すると

\frac{dv}{dt} = 1 - v

となる.
この微分方程式を前進Euler法で計算する.

#include <array>    // for std::array
#include <cmath>    // for std::exp
#include <fstream>  // for std::ofstream
#include <iostream> // for std::cerr, std::endl
#include <utility>  // for std::make_pair, std::move, std::pair
#include <vector>   // for std::vector
#include <boost/format.hpp> // for boost::format

namespace{
    using mypair = std::pair<std::vector<double>, std::vector<double> >;

    static auto constexpr V0 = 0.0;
    static auto constexpr TEND = 1.0;
    static std::array<double, 3> const dtarray = {0.1, 0.01, 0.001};
    static auto const size = static_cast<int>(dtarray.size());

    double dvdt(double v);

    std::vector<double> error_assesment(std::vector<double> const &vendvec);

    std::vector<double> forward_Euler(std::vector<double> const &stependvec);

    std::vector<double> make_stependvec();

    bool output_result(std::vector<double> const &dterrvec);
}

int main(){
    auto const stependvec = make_stependvec();

    auto const vendvec = forward_Euler(stependvec);

    auto const dterrvec = error_assesment(vendvec);

    if (!output_result(dterrvec)){
        std::cerr << "output file not open" << std::endl;
        return -1;
    }
}

namespace{
    double dvdt(double v){
        return 1 - v;
    }

    std::vector<double> error_assesment(std::vector<double> const &vendvec){
        std::vector<double> vexactvec(size, 1.0 - exp(-TEND));
        std::vector<double> dterrvec;

        for (auto i = 0; i < size; i++){
           dterrvec.push_back((vendvec[i] - vexactvec[i]) / vexactvec[i]);
        }

        return dterrvec;
    }

    std::vector<double> forward_Euler(std::vector<double> const &stependvec){
        std::vector<double> vendvec(size, 0.0);

        for (auto i = 0; i < size; i++){
            for (auto j = 0; j < stependvec[i]; j++){
                vendvec[i] += dvdt(vendvec[i]) * dtarray[i];
            }
        }

        return vendvec;
    }

    std::vector<double> make_stependvec(){
        std::vector<double> stependvec;

        for (auto i = 0; i < size; i++){
            stependvec.push_back(TEND / dtarray[i]);
        }

        return stependvec;
    }

    bool output_result(std::vector<double> const &dterrvec){
        std::ofstream ofs;
        ofs.open("data_forward_Euler.txt");

        if (!ofs){
            return false;
        }

        for (auto i = 0; i < size; i++){
            ofs << boost::format("%.7f %.7f\n") % dtarray[i] % dterrvec[i];
        }

        return true;
    }
}

出力ファイルdata_forward_Euler.txtからデータを読み取り, gnuplotを用いてグラフを描き, forward_Euler.pngに保存する.

set title 'forward Euler method'

set xlabel 'dt'
set ylabel 'error'
set key left top

set logscale x
set logscale y

set xr[1.0e-4:1.0]
set yr[1.0e-4:1.0]

set grid

plot 'data_forward_Euler.txt' using 1:2 title 'forward Euler'
replot x title "\Delta t"

set terminal pngcairo
set output "forward_Euler.png"
replot

forward_Euler.pngは次のようになる.

Adams Boshforth法 (Adams Bashforth Method)

空気抵抗のある質点の自由落下の運動をAdams Bashforth法で計算し, 理論解との誤差を見る.
微分方程式

\frac{dv}{dt} = 1 - v

の理論解は

v(t) = 1 - \exp(-t)

である.

#include <fstream>  // for std::ofstream
#include <iostream> // for std::cerr, std::cout
#include <utility>  // for std::make_pair, std::move, std::pair
#include <vector>   // for std::vector
#include <boost/format.hpp> // for boost::format

namespace{
    using mypair = std::pair<std::vector<double>, std::vector<double> >;

    static auto constexpr V0 = 0.0;
    static auto constexpr DT = 0.1;
    static auto constexpr TEND = 1.0;
    static auto constexpr STEPEND = TEND / DT;

    // Adams Bashforth method
    mypair Adams_Bashforth();

    // open output file
    bool output_result(mypair const &tvpair);

    // motion equation
    double dvdt(double v);
}

int main(){
    auto const tvpair = Adams_Bashforth();

    if (!output_result(tvpair)){
        std::cerr << "output file not open" << std::endl;
        return -1;
    }
}

namespace{
    double dvdt(double v){
        return 1 - v;
    }

    mypair Adams_Bashforth(){
        std::vector<double> tvec(STEPEND + 2, 0.0), vvec(STEPEND + 2, 0.0);

        // vの初期条件
        vvec[0] = V0;

        // 2nd order Runge Kutta for vvec[1]
        auto const k1 = DT * dvdt(vvec[0]);
        auto const k2 = DT * dvdt(vvec[0] + k1 * 0.5);
        vvec[1] = vvec[0] + k2;

        for (auto i = 1; i < STEPEND; i++){
            tvec[i] = static_cast<double>(i) * DT; 
            
            // Adams Bashforth method
            vvec[i + 1] = vvec[i] + (DT / 2) * (-dvdt(vvec[i - 1]) + 3 * dvdt(vvec[i]));      
        }

        return std::make_pair(std::move(tvec), std::move(vvec));
    }

    bool output_result(mypair const &tvpair){
        std::ofstream ofs;
        ofs.open("data_Adams_Bashforth.txt");

        if (!ofs){
            return false;
        }

        auto const [tvec, vvec] = tvpair;
        for (int i = 0; i <= STEPEND; i++){
            ofs << boost::format("%.7f %.7f\n") % tvec[i] % vvec[i];
        }

        return true;
    }
}

出力ファイルdata_Adams_Bashforth.txtからデータを読み取り, gnuplotを用いてグラフを描き, Adams_Borshforth.pngに保存する.

set title 'Adams Bashforth method'

set xlabel 'time'
set ylabel 'velocity'
set key left top

set xr[0.0:1.1]
set yr[0.0:0.7]

set xtics 0.1
set ytics 0.1

set grid

plot 'data_Adams_Bashforth.txt' using 1:2 title 'Adams Bashforth'
replot 1 - exp(-x) title 'theoritical solution 1 - exp(-t)'

set terminal pngcairo
set output 'Adams_Bashforth.png'
replot

Adams_Borshforth.pngは次のようになる.

Runge Kutta法 (Runge-Kutta Method)

Euler法は$y(t+dt)$をTaylor展開した時の$1$次の項まで考えたものと捉えることができる.
高階の微分を一回の微分関数$f(x,y)$のみで表す式をRunge Kutta法という.
しかし, 特に断らない限り, $4$次のRunge Kutta法の次の形を単にRunge Kutta法という.

Runge Kutta法は一階常微分方程式

\frac{dy}{dx} = f(x,y)

と初期条件

y(x_0) = y_0

について, $x_{i+1}=x_i+dh$に対応する$y_{i+1}$の値を

y_{i+1} = y(x_i)+\frac{dh}{6}(k_1+2k_2+2k_3+k_4) \\
k_1 = f(x,y) \\
k_2 = f(x+\frac{dh}{2}, y+dh\frac{k_1}{2}) \\
k_3 = f(x+\frac{dh}{2}, y+dh\frac{k_2}{2}) \\
k_4 = f(x+dh, y+dhk_3)

と近似する.

例えば, 初期条件$v(t_0)=v_0$で粘性抵抗を受けながら自由落下する物体の運動を, 刻み幅$dt$で増加する時刻$t_1, t_2, \cdots$に対して順に求める.
この物体の運動方程式は

\frac{dv}{dt} = g - cv

のようになる (ただし, $c$は粘性抵抗係数).

#include <iomanip>  // for std::ios::fixed, std::setiosflags, std::setprecision
#include <iostream> // for std::cout
#include <string>   // for std::string

namespace{
    static auto constexpr MAXBUFSIZE = 32;
    static auto constexpr c = 0.1;
    static auto constexpr g = 9.80665;

    template <typename T, bool check_positive_num = true>
    T input_parameter(std::string const & str){
        T param;
        while (true){
		    std::cout << str;
            std::cin >> param;

		    std::cin.clear();
		    std::cin.ignore(MAXBUFSIZE, '\n');

            // 初期条件として<=0の場合に再入力が要らない時は<type, false>
            if constexpr (check_positive_num){
                if (!std::cin.fail() && param > static_cast<T>(0)){ 
			        break;
		        }
            } 
            else{
                if (!std::cin.fail()){
                    break;
                }
            }
	    }

        return param;
    }

    // motion equation
    double f(double v);

    void Runge_Kutta(int n, double tmax, double v0);
}

int main(){
    // 分割数, 最大時間, 初期条件の入力
    auto const n = input_parameter<int> ("Enter the partition number n\n");
    auto const tmax = input_parameter<double> ("Enter the maximum time t_max[s]\n");
    auto const v0 = input_parameter<double, false> ("Enter the initial velocity v_0[m/s]\n");
    
    Runge_Kutta(n, tmax, v0);
}

namespace{
    double f(double v){
        return g - c * v;
    }

    void Runge_Kutta(int n, double tmax, double v0){
        // 分割幅
        auto const dt = tmax / static_cast<double>(n);
        // vの初期化
        auto v = v0;

        std::cout << std::setprecision(7) << std::setiosflags(std::ios::fixed);
        for (auto i = 0; i <= n; i++){
            auto const t = static_cast<double>(i) * dt;

            std::cout << "when the time is " << t << "[s], velocity is " << v << "[m/s]" << std::endl;

            // Runge Kutta method
            auto const k1 = f(v);
            auto const k2 = f(v + k1 * dt * 0.5);
            auto const k3 = f(v + k2 * dt * 0.5);
            auto const k4 = f(v + k3 * dt);
            v += dt / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4);            
        }
    }
}

物体の運動が高次元であれば, 配列として考えればよい.

その他の解法

Euler法は前進差分で求めたので, 前進Euler法(forward Euler method)と呼ばれることもある. 後退差分で求める方法を後退Euler法(backward Euler method)と呼ぶ.

$2$次のRunge-Kutta法として, 係数を$1$としたものをHeun法(Heun method)や改良Euler法(improved Euler method)という.
係数を$1/2$としたものは中点法(midpoint method)や修正Euler法(modified Euler method)と呼ばれる.

boost::odeint

常微分方程式の数値解析ライブラリであるboost::odeintを用いる.
実際にC++で常微分方程式を解く時はboost::odeintを用いる.

例えば, $2$次元のLotka Volterra方程式

\frac{dx}{dt} = -x (\alpha - \beta y) \\
\frac{dy}{dt} = -y (\gamma - \delta x)

をあるパラメタ$\alpha, \beta, \gamma, \delta$に対して解いてみる.

#include <array>    // for std::array
#include <fstream>  // for std::ofstream
#include <iostream> // for std::cerr, std::endl
#include <boost/format.hpp>         // for boost::format
#include <boost/numeric/odeint.hpp> // for boost::numeric::odeint

namespace{
    using state = std::array<double, 2>;

    static auto constexpr DT = 0.05;
    static auto constexpr T0 = 0.0;
    static auto constexpr TMAX = 5.0;

    class Lotka_Volterra_system final{
        double const alpha_;
        double const beta_;
        double const delta_;
        double const gamma_;

        public:
            Lotka_Volterra_system(double alpha, double beta, double delta, double gamma)
                : alpha_(alpha), beta_(beta), delta_(delta), gamma_(gamma)
            {}
            Lotka_Volterra_system() = delete;
            ~Lotka_Volterra_system() = default;

            void operator()(state const &x, state &dx, double t) const; //関数内でalpha等を変更しないためのconst
    }; 
}

int main(){   
    Lotka_Volterra_system System(2.0, 3.0, 5.0, 4.0);
    state State = {1.0, 0.5};

    //6th-Adams Bashforth Moulton method
    boost::numeric::odeint::adams_bashforth_moulton<6, state> Stepper;

    // open output file
    std::ofstream ofs("data_Lotka_Volterra.txt");

    if (!ofs){
        std::cerr << "output file not open" << std::endl;
        return -1;
    }

    //time = T0からtime = TMAXまで時間発展をDT刻みで計算
    boost::numeric::odeint::integrate_const(
        Stepper,
        System,
        State,
        T0,
        TMAX,
        DT,
        [&ofs](state const &x, double t) {ofs << boost::format("%.2f %.7f %.7f\n") % t % x[0] % x[1];}
    );
}

namespace{
    void Lotka_Volterra_system::operator()(state const &x, state &dx, double t) const{
        dx[0] = x[0] * (alpha_ - beta_ * x[1]);
        dx[1] = - x[1] * (gamma_ - delta_ * x[0]);
    }
}

出力ファイルdata_Lotka_Volterra.txtからデータを読み取り, gnuplotを用いてグラフを描き, Lotka_Volterra.pngに保存する.

set title 'boost::odeint Lotka Volterra'

set xlabel 'time'
set ylabel 'x and y'
set key left top

set xr[0:5.1]
set yr[0.3:1.2]

set grid

plot 'data_Lotka_Volterra.txt' using 1:2 w lp title 'x'
replot "data_Lotka_Volterra.txt" using 1:3 w lp title 'y'

set terminal pngcairo
set output 'Lotka_Volterra.png'
replot

Lotka_Volterra.pngは次のようになる.

補遺 boostついて

boostはC++のライブラリであり, 様々な機能がある.
この記事では常微分方程式の数値計算に特化したboost::odeintについて説明する.
ここで紹介していない便利なコマンドもたくさんあるので, 参考文献を適宜参考にしてほしい.