最終更新: 2022/7/5

一次元の精度保証付き数値積分 (特異点なし)

柏木雅英

1. はじめに

1次元の有界閉区間 [a,b] 上の数値積分 \( \displaystyle I = \int_a^b f(x) dx \) を精度保証付きで計算する。被積分関数がベキ級数演算が可能な関数の組み合わせで書かれていて、積分区間内に特異点が無い場合に計算可能。今のところ、

ベキ級数演算を利用する方法
複合Newton-Cotes公式を利用する方法

の2つが提供されている。将来は、Gauss-Legendre公式等も提供される予定。

2. ベキ級数演算を利用した精度保証付き数値積分

ベキ級数演算を利用しており、基本的にアルゴリズムは、

ベキ級数演算についての6節

に従っている。

インクルードファイルは、defint.hppを使う。また下請けとして、psa.hpp, interval.hpp, またはそれらの下請けを使う。

参考ファイルは、test/test-defint.cc, example/defint-example.ccなど。
KVライブラリのwebデモの中の 1変数の数値積分で生成されるファイルも参考になるかも。

例えば、数値積分 \( \displaystyle \int_{-1}^1 \frac{1-x^2}{1+x^2} dx \) を計算する問題を考える。このとき、次のようなプログラムを作成すればよい (sample1d.cc)。

#include <iostream>
#include <kv/defint.hpp>

typedef kv::interval<double> itv;

struct Func {
    template <class T> T operator() (const T& x) {
        return (1 - x * x) / (1 + x * x);
    }
};

int main() {
    std::cout.precision(17);

    std::cout << kv::defint(Func(), itv(-1), itv(1), 12, 10) << "\n";
    std::cout << kv::defint_autostep(Func(), itv(-1), itv(1), 12) << "\n";
}

例のように、()演算子をオーバーロードしてT型を受け取ってT型を返す関数オブジェクトで解きたい問題 (被積分関数) を定義する。

精度保証付き定積分を行う関数は、defintとdefint_autostepの2種類がある。

defintは、積分区間を等間隔に分割して積分する関数で、引数は

  defint(被積分関数, 区間の左端, 区間の右端、内部で用いる級数の次数, 分割数)

の意味。戻り値が計算結果。この例だと12次のpsaを用いて[-1,1]を10分割している。
区間の両端点は区間であり、interval型である必要がある (端点が円周率のような場合を想定している)。

defint_autostepは、区間の分割を自動的に行う。引数は、

  defint_autostep(被積分関数, 区間の左端, 区間の右端、内部で用いる級数の次数, 許容誤差(省略可))

の意味。戻り値が計算結果。

このプログラムをコンパイルして実行すると、

[1.1415926535894358,1.1415926535902046]
[1.1415926535897868,1.1415926535898005]

と表示される。

なお、端点を表す区間型をinterval<double>でなく interval<dd>にするなどすれば、それに応じた高精度計算を行うことが出来る。すなわち、defintやdefint_autostepはintervalの内部で使う型を切り替えられるtemplate関数として作られている。

defint_autostepの許容誤差は、全体の許容誤差ではなく、自動分割するときの 1ステップ毎の許容誤差(の目安)。 defaultは区間の端点の型のmachine epsilin。ステップ幅の自動調節の詳細は、ベキ級数演算についての10節の通り。

3. Newton-Cotesの公式による数値積分

複合Newton-Cotes公式による数値積分を行う。アルゴリズムは、

精度保証付き二重積分についてのpp.27-28

に従っている。誤差項の計算には積分区間上での高階微分の評価が必要だが、それには

ベキ級数演算についての第4節

の方法と、

大域最適化

の方法を併用している。

インクルードファイルは、defint-newtoncotes.hppを使う。また下請けとして、 psa.hpp, interval.hpp, highderiv.hpp, optimize.hpp, またはそれらの下請けを使う。

参考ファイルは、test/test-defint-newtoncotes.cc 。

例えば、数値積分 \( \displaystyle \int_{-1}^1 \frac{1-x^2}{1+x^2} dx \) を計算する問題を考える。このとき、次のようなプログラムを作成すればよい (sample1d-nc.cc)。

#include <iostream>
#include <kv/defint-newtoncotes.hpp>

typedef kv::interval<double> itv;

struct Func {
    template <class T> T operator() (const T& x) {
        return (1 - x * x) / (1 + x * x);
    }
};

int main() {
    std::cout.precision(17);

    std::cout << kv::defint_newtoncotes(Func(), itv(-1), itv(1), 2, 20) << "\n";
    std::cout << kv::defint_newtoncotes(Func(), itv(-1), itv(1), 6, 0) << "\n";
}

defint_newtoncotesを呼び出す。引数は

  defint_newtoncotes(被積分関数, 区間の左端, 区間の右端, 分割区間の複合数, 分割数, 誤差項の計算時の分割数)

の意味。戻り値が計算結果。

分割区間の複合数 : 1から7までのint型で指定する。1なら複合台形公式、 2なら複合Simpson則。ここに指定した数の小区間を組にして計算を行う。公式の次数は、複合数1〜7に対してそれぞれ2,4,4,6,6,8,8次となる。このような次数の上がり方をするため、ここに指定する数は偶数がお勧め。
分割数 : ここに指定した数の等間隔区間に与えられた区間を分割して計算する。「分割区間の複合数」の倍数でなければならないが、そうでない場合は自動で増やして計算する (例えば、複合数3で分割数10ならば、分割数は12まで自動的に増やされる)。なお、分割数に0を指定することが出来て、その場合は丸め誤差の大きさと誤差項の大きさから最も精度が高くなるような分割数を自動推定して計算する。自動推定のアルゴリズムは、精度保証付き二重積分についてのpp.34-36の方法に従っている。
誤差項の計算時の分割数 : 誤差項をType-I PSAを用いて計算するときの区間分割数。省略可能で、defaultでは 1 (分割しない)。多く分割すれば誤差項が小さくなるが、それだけ計算時間がかかる。なお、計算時にゼロ除算が起こるなど (分数を含む被積分関数の場合分母の区間評価のover-estimationでよく起こる) 、計算が正常に行えない場合はここに指定した分割数よりも多く分割されることがある。

このプログラムをコンパイルして実行すると、

[1.1406859472725467,1.1422226139392174]
[1.1415926535897606,1.1415926535898265]

と表示される。

4. おわりに

前者は Type-II PSAの応用例の一つとしてとりあえず作ってみたもので、一般的な関数では後者の方法の方が高速で使い勝手がいいものと思われる。しかし、局所的に鋭いピークを持つような関数では、後者の方法では誤差評価が大きくなってしまい、前者の自動ステップ幅調整が機能して前者の方がよいとも考えられる。

そのうち組織的な比較をしてみたい。