数値積分を実行する

数値積分を実行する。

被積分関数までは解析的に書けているけれど、
積分が解析的に実行できず、さらにその値が他の計算に必要な場合。
しかも高精度に。

参考：ニュートン・コーツの公式(Wikipedia)

# include <functional>
# include <iostream>
# include <cmath>

enum IntegralMethod{
    Trapezoidal = 0,
    Simpson = 1,
    Boole = 2,
};

template<IntegralMethod>
class Integrater{
public:
    /** 一次元定積分 @f$ \int_a^b f(x) dx @f$ を積分する */
    static double calc(int N /** 空間刻み */,
                       double a /** 下限 */, double b /** 上限 */,
                       std::function<double(double)> f /** 積分する関数 */);
    /** 一次元定積分 @f$ \int_a^b f(x) dx @f$ を積分する(精度保証付き)。
     *  ただし打ち消して0になる場合は使えない。 */
    static double accurate(double a /** 下限 */, double b /** 上限 */,
                           std::function<double(double)> f /** integrand */,
                           double threshold=1e-9 /** 収束の閾値 */,
                           std::ostream* report = nullptr /** 収束の過程を表示 */)
    {
        int N = 128;
        double result_p = calc(N, a, b, f);
        double result_n;
        while(true){
            N *= 2;
            result_n = calc(N, a, b, f);
            if(report){
                *report << N << " " << result_n << std::endl;
            }
            if(std::abs(result_n - result_p) / std::abs(result_n) < threshold){
                break;
            }
            result_p = result_n;
        }
        return result_n;
    }
};

template<>
double Integrater<Trapezoidal>::calc(int N, double a, double b, std::function<double(double)> f) {
    long double sum = 0.0;
    double dx = (b - a) / N;
    for(int i=0;i<N;++i){
        double x_l = a + dx*i;
        double x_r = a + dx*(i+1);
        sum += f(x_l) + f(x_r);
    }
    return sum * dx / 2;
}

template<>
double Integrater<Simpson>::calc(int N, double a, double b, std::function<double(double)> f) {
    long double sum = 0.0;
    double dx = (b - a) / N;
    for(int i=0;i<N;++i){
        double x_l = a + dx*i;
        double x_r = a + dx*(i+1);
        double x_m = (x_l + x_r) / 2;
        sum += f(x_l) + 4*f(x_m) + f(x_r);
    }
    return sum * dx / 6;
}

template<>
double Integrater<Boole>::calc(int N, double a, double b, std::function<double(double)> f) {
    long double sum = 0.0;
    double dx = (b - a) / N;
    for(int i=0;i<N;++i){
        double x_l = a + dx*i;
        double x_r = a + dx*(i+1);
        double x_m = (x_l + x_r) / 2;
        double x_lm = (x_l + x_m) / 2;
        double x_rm = (x_r + x_m) / 2;
        sum += 7*f(x_l) + 32*f(x_lm) + 12*f(x_m) + 32*f(x_rm) + 7*f(x_r);
    }
    return sum * dx / (7+32+12+32+7);
}

台形公式(Trapezoidal)で行くつもりだったけど、
高精度を出そうとすると、台形公式は非常に効率がわるい。
そこで高次のアルゴリズムとしてSimpson, Booleを実装した。

最初はaccurateの引数にstd::functionで計算手法(Trapezoidal,...)を貰おうかと思ったけど、
この方がすっきりする。calcをvirtualにして継承するよりは、こちらの方が好み。

地味なこだわりは

std::ostream* report = nullptr

でstreamをもらう事。これなら必ずしも出力する必要がない。

if(report){
...
}

の評価は通常非常に軽いはず。

もっといい実装キボンヌ