C++で数値計算 (数値積分)

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

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

目次

• 数値積分 (Numerical Integral)
  • 台形積分 (Trapezoidal Integral)
  • Simpson則による積分 (Simpson Rule)
  • その他の解法

数値積分 (Numerical Integral)

関数$y=f(x)$と$x$軸とある区間で囲まれた面積を, 微小区間$dh$に分割して近似値として求める.
数値積分の方法としては台形積分やSimpson則などがある.

台形積分 (Trapezoidal Integral)

積分区間を微小な幅に分割し, 面積を台形として近似する.
区分求積法における長方形で近似する方法よりも精度が高い.

$f(x)=x^2$の定積分を台形積分の方法を用いて求める.

#include <cassert>  // for assert
#include <iomanip>  // for std::ios::fixed, std::setiosflags, std::setprecision
#include <iostream> // for std::cout, std::endl
#include <utility>  // for std::make_pair, std::pair 

namespace{
    static auto constexpr MAXBUFSIZE = 32;

    std::pair<double, double> integ_interval();

    int partition_num();

    double trapezoidal_integral(int n, double x0, double x1);

    double y(double x);
}

int main(){
    auto const [low_end, up_end] = integ_interval();
    auto const n = partition_num();
    auto const res = trapezoidal_integral(n, low_end, up_end);

    std::cout << std::setprecision(7) << std::setiosflags(std::ios::fixed);
    std::cout << "When the partition number is " << n << ", the numerical integral in ["
    << low_end << "," << up_end << "] of y = x^2 is " << res << std::endl; 
}

namespace{
    std::pair<double, double> integ_interval(){
        double low_end, up_end;

        while (true){
            // 積分の下端と上端の入力
		    std::cout << "Enter endpoints of for the integration interval" << std::endl;
		    std::cin >> low_end >> up_end;
            
		    if (!std::cin.fail() && low_end < up_end){
			    break;
		    }

		    std::cin.clear();
		    std::cin.ignore(MAXBUFSIZE, '\n');
	    }
    
        return std::make_pair(low_end, up_end);
    }

    int partition_num(){
        int n;

        while (true){
            // 分割数の入力
		    std::cout << "Enter the partition number" << std::endl;
            std::cin >> n;

		    if (!std::cin.fail()) {
			    break;
		    }
            
		    std::cin.clear();
		    std::cin.ignore(MAXBUFSIZE, '\n');
	    }

        return n;
    }

    double trapezoidal_integral(int n, double x0, double x1){
        // x0 >= x1の時, 強制終了
        assert(x0 < x1);

        // 分割幅
        auto const dh = (x1 - x0) / static_cast<double>(n); 
        // 面積の初期化
        double sum = 0.0; 

        // trapezoidal integral
        for (auto i = 1; i < n; i++) {
            auto const x = x0 + static_cast<double>(i) * dh;
            sum += y(x);
        }

        return dh * ((y(x0) + y(x1)) / 2.0 + sum);
    }

    double y(double x){
        return x * x;
    }
}

Simpson則による積分 (Simpson Rule)

$f(x)=\sin(x)$の定積分をSimpson則を用いて求める.

#include <cassert>  // for assert
#include <cmath>    // for std::sin
#include <iomanip>  // for std::ios::fixed, std::setiosflags, std::setprecision
#include <iostream> // for std::cout, std::endl
#include <utility>  // for std::make_pair, std::pair 

namespace{
    static auto constexpr MAXBUFSIZE = 32;

    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;
    }

    std::pair<double, double> integ_interval();

    double Simpson(int n, double x0, double x1);

    double f(double x);
}

int main(){
    auto const [low_end, up_end] = integ_interval();
    auto const n = input_parameter<int> ("Enter the partition number\n");
    auto const res = Simpson(n, low_end, up_end);

    std::cout << std::setprecision(7) << std::setiosflags(std::ios::fixed);
    std::cout << "When the partition number is " << n << ", the numerical integral in ["
    << low_end << "," << up_end << "] of y = sin(x) is " << res << std::endl; 
}

namespace{
    std::pair<double, double> integ_interval(){
        double low_end, up_end;

        while (true){
            // 積分の下端と上端の入力
		    std::cout << "Enter endpoints of for the integration interval" << std::endl;
		    std::cin >> low_end >> up_end;
            
		    if (!std::cin.fail() && low_end < up_end){
			    break;
		    }

		    std::cin.clear();
		    std::cin.ignore(MAXBUFSIZE, '\n');
	    }
    
        return std::make_pair(low_end, up_end);
    }

    double Simpson(int n, double x0, double x1){
        // x0 >= x1の時, 強制終了
        assert(x0 < x1);

        // 分割幅
        auto const dh = (x1 - x0) / static_cast<double>(n); 
        // 面積の初期化
        double sum = 0.0; 

        // Simpson's rule
        for (auto i = 1; i <= 2 * n; i++) {
            auto const x = x0 + 0.5 * static_cast<double>(i) * dh;

            if (i == 0 || i == 2 * n){
                sum += f(x);
            }
            else if (i % 2 == 0){
                sum += 2 * f(x);
            }
            else{
                sum += 4 * f(x);
            }
        }

        return (dh / 6) * sum;
    }

    double f(double x){
        return std::sin(x);
    }
}

その他の解法

$1$変数の数値積分としてはNewton Cotesの公式が主に挙げられる.
Newton-Cotesの公式とは, 積分区間を等間隔に分割して被積分関数の値を考えて数値積分を行う手法の総称である.
例えば, 区分求積法の定義にある中点法や台形積分, Simpson則がある.
Newton Cotesの公式以外では, Romberg積分やGauss求積法などがある.