最終更新: 2019/11/21

数値積分の精度保証

柏木 雅英

1. はじめに

数値積分の精度保証を行うプログラムがいくつか用意されている。 など。

2. 1次元の数値積分

ベキ級数演算を利用しており、基本的にアルゴリズムは、 に従っている。 積分区間が有限区間で積分区間内に特異点が無ければ原理的には 大体計算できるはず。

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

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

例えば、関数 (1-x2) / (1+x2) を-1から1まで 数値積分する問題を考える。 このとき、次のようなプログラムを作成すればよい (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.1415926535897686,1.1415926535898205]
と表示される。

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

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

3. 長方形領域及び三角形領域での二重積分

インクルードファイルは、doubleintegral.hppを使う。 また下請けとして、psa.hpp, interval.hpp, またはそれらの下請けを使う。 参考ファイルはtest-doubleintegral.cc。
KVライブラリのwebデモの中の 二重積分(長方形領域)二重積分(三角形領域)、 で生成されるファイルも参考になるかも。

使っているアルゴリズムは、 精度保証付き二重積分の方法まとめ (Version: 2019/11/20) の第2節に説明がある。

使い方は、

  doubleintegral(被積分関数, x1, x2, y1, y2, 次数, 分割数)
のような感じ。戻り値が計算結果。

x1, x2, y1, y2はinterval型で、次の図のような積分範囲になる。

次数は用いる級数の次数、分割数をnとすると、x, yそれぞれをn等分して 計算する。

なお、

  doubleintegral(被積分関数, x1, x2, y1, y2, 次数, 分割数, true)
のように末尾にフラグtrueを設定すると、

のような三角形上の積分になる。

  doubleintegral_triangle(被積分関数, x1, y1, x2, y2, x3, y3, 次数, 分割数)
で、点(x1, y1)、点(x2, y2)、点(x3, y3)で構成される三角形上の 積分を行う。これは内部で変数変換してdoubleintegral(..., true)を呼び出している。

4. 端点に特異性のある1次元の数値積分

端点に特異性のある関数の自動積分も試験的に作っている。 アルゴリズムは、 のようなものを使用している。

現在のところ、

が出来るようになっている。test-defint-singular.ccは使い方の例。 まだ試作品で関数名等がこなれてなく、使い方は今後変更される可能性がある。

lobachevsky関数やgamma関数の計算の内部で、特異関数の積分機能を使っている。

4.1 ファイル構成

下請けとして、defint.hppまたはそれらの下請けを使う。

4.2 除去可能特異点を端点に持つ関数の数値積分

sin(x)/x を0から1まで積分する例を示す。

この例は、原点で0/0になっているので普通に解かせるとゼロ除算が起きてしまう。 実際には分母分子をそれぞれTaylor展開すると原点では分母分子ともに 定数項が0になっていて、「約分」すると正常に計算できる。 そこで、これらを積分するには、

の2種類を用意する必要がある。

[defint_singular]

まず一つ目の例を示す。 (sample-singular.cc)。
#include <iostream>
#include <kv/defint-singular.hpp>

typedef kv::interval<double> itv;

struct Sinc {
    template <class T> T operator()(const T& x) {
        return sin(x) / x;
    }
};

struct Sinc_s {
    template <class T> T operator()(const T& x) {
        return div_tn(sin(x), 1);
        // return div_reduce(sin(x), x, 1);
    }
};

int main() {
    std::cout.precision(17);
    itv r1, r2;

    r1 = kv::defint_singular(Sinc_s(), (itv)0., (itv)0.125, 12);
    r2 = kv::defint_autostep(Sinc(), (itv)0.125, (itv)1., 12);

    std::cout << r1 << "\n";
    std::cout << r2 << "\n";
    std::cout << r1 + r2 << "\n";
}
これを実行すると、
[0.12489154390467224,0.12489154390467227]
[0.82119152646164261,0.82119152646727767]
[0.94608307036631478,0.94608307037194995]
となる。 Sincは元の被積分関数、Sinc_sはそれを特異点(0)で「約分」した関数である。 Sinc_sの中で使われているdiv_tn(p, n)は、pをベキ級数型、nを自然数として、 pをtnで割る関数。 div_reduce(p1, p2, n)は、p1, p2をベキ級数型、nを自然数として、 p1, p2をともにtnで割ってから除算を行う関数。 よって、 のどちらでもよい。

defint_singular関数は、

    defint_singular(f, start, end, order)
      f: 被積分関数
      start: 積分区間の左端
      end: 積分区間の右端
      order: 内部で使うベキ級数の次数
のような使い方。 積分区間の左端(start)でTaylor展開を作成し、積分値を計算する。 fにはdiv_tnまたはdiv_reduceで強制的に約分を行う関数を渡すが、 このとき約分を行える (上の例では分母のsin(x)の展開の定数項が0になっている) ことを保証するのはユーザの役割である。

この例では0から0.125までをdefint_singularで約分した関数を0中心のTaylor展開で 積分し、0.125から1までは通常の方法で元の被積分関数を自動積分している。

[defint_singular_autostep]

先の例だと、特異点周りで特別扱いする領域とそれ以外の領域を 手動で分割していたが、最適な分割位置は必ずしも自明でない。 そこで、分割位置を自動で決めてくれる関数defint_singular_autostepも用意した。 使用例を示す (sample-singular2.cc)。
#include <iostream>
#include <kv/defint-singular.hpp>

typedef kv::interval<double> itv;

struct Sinc {
    template <class T> T operator()(const T& x) {
        return sin(x) / x;
    }
};

struct Sinc_s {
    template <class T> T operator()(const T& x) {
        return div_tn(sin(x), 1);
    }
};

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

    std::cout<< kv::defint_singular_autostep(Sinc(), Sinc_s(), (itv)0., (itv)1., 12) << "\n";
}
これを実行すると、
[0.94608307036684913,0.94608307036909978]
となる。

defint_singular_autostep関数は、

    defint_singular_autostep(f1, f2, start, end, order, epsilon)
      f1: 元の被積分関数
      f2: 特異点除去を行った被積分関数
      start: 積分区間の左端
      end: 積分区間の右端
      order: 内部で使うベキ級数の次数
      epsilon: 1ステップ毎の誤差上限
のような使い方。epsilonは省略可(machine epsilonになる)。 最初はstartの点でf2を使って展開し、誤差がepsilon程度になるように step幅を調整して積分する。それ以降は、f1に対してdefint_autostepを 呼び出してるだけ。

本5節の端点特異性を持つ関数に対する積分関数全てについて、非autostep版と autostep版を両方作成した。通常autostep版を使えば良いだろう。 以降は、autostep版のみを例に説明する。

以下は、x2/(1-cos(x)) を0から1まで積分する例 (sample-singular3.cc)。

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

typedef kv::interval<double> itv;

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

struct RemSing_s {
    template <class T> T operator()(const T& x) {
        return div_reduce(x * x, 1 - cos(x), 2);
    }
};

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

    std::cout<< kv::defint_singular_autostep(RemSing(), RemSing_s(), (itv)0., (itv)1., 12) << "\n";
}
実行すると、
[2.057270784424849,2.0572707844408656]
と表示される。

4.3 ベキ型特異点を端点に持つ関数の数値積分

[defint_power]

defint_powerは、積分区間を[a,b]として、 (x-a)pg(x) (pは非整数) の形の積分を行う。

sqrt(x)cos(x)を0から1まで積分する例を示す (sample-singular4.cc)。

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

typedef kv::interval<double> itv;

struct Func1 {
    template <class T> T operator() (const T& x) {
        return sqrt(x) * cos(x);
    }
};

struct Func1_s {
    template <class T> T operator() (const T& x) {
        return cos(x);
    }
};
int main() {
    std::cout.precision(17);

    std::cout << kv::defint_power_autostep(Func1(), Func1_s(), (itv)0., (itv)1., 12, (itv)0.5) << "\n";
}
実行すると、
[0.53120268308451157,0.5312026830845168]
となる。

このように、

の2つを渡す。

defint_power関数は、

    defint_power(g, start, end, order, p)
      g: (x-start)pg(x) のg
      start: 積分区間の左端
      end: 積分区間の右端
      order: 内部で使うベキ級数の次数
      p: ベキ数
defint_power_autostep関数は、
    defint_power_autostep(h, g, start, end, order, p, epsilon)
      h: 元の被積分関数 (x-start)pg(x)
      g: (x-start)pg(x) のg
      start: 積分区間の左端
      end: 積分区間の右端
      order: 内部で使うベキ級数の次数
      p: ベキ数
      epsilon: 1ステップ毎の誤差上限
という使い方。

[defint_power2]

defint_power2は、積分区間を[a,b]として、 (f(x))p (pは非整数) の形の積分を行う。 ただし、f(a) = 0 とする。

sqrt(sin(x))を0から1まで積分する例を示す (sample-singular5.cc)。

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

typedef kv::interval<double> itv;

struct SqrtSin {
    template <class T> T operator()(const T& x) {
        return sqrt(sin(x));
    }
};

struct SqrtSin_s {
    template <class T> T operator()(const T& x) {
        return sin(x);
    }
};

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

    std::cout << kv::defint_power2_autostep(SqrtSin(), SqrtSin_s(), (itv)0., (itv)1., 12, (itv)0.5) << "\n";
}
実行すると、
[0.64297763465831214,0.64297763465831726]
となる。

この例は、f(a)=0かつf'(a)≠0 、すなわちaという解の多重度は1だったが、 多重度が1より大きい場合、すなわち n > 1 があって f(i)(a)=0 (i < n)、f(n)(a)≠0の場合、 そのn (multiplicity) を引数に指定する。

そのような例として、 sqrt(1-cos(x))を0から1まで積分する例を示す (sample-singular6.cc)。 1-cos(x)のTaylor展開はx^2/2 - x^4/24 + ... なので、多重度は2である。

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

typedef kv::interval<double> itv;

struct Func2 {
    template <class T> T operator()(const T& x) {
        return sqrt(1 - cos(x));
    }
};

struct Func2_s {
    template <class T> T operator()(const T& x) {
        return 1 - cos(x);
    }
};

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

    std::cout << kv::defint_power2_autostep(Func2(), Func2_s(), (itv)0., (itv)1., 12, (itv)0.5, 2) << "\n";
}
[0.34624880247392419,0.34624880250131529]

defint_power2関数は、

    defint_power2(f, start, end, order, p, multiplicity)
      f: f(x)p のf
      start: 積分区間の左端
      end: 積分区間の右端
      order: 内部で使うベキ級数の次数
      p: ベキ数
      multiplicity: fn(a)≠0となる最小のn (default: 1)
defint_power2_autostep関数は、
    defint_power2_autostep(h, f, start, end, order, p, multiplicity, epsilon)
      h: 元の被積分関数 f(x)p
      f: f(x)p のf
      start: 積分区間の左端
      end: 積分区間の右端
      order: 内部で使うベキ級数の次数
      p: ベキ数
      multiplicity: fn(a)≠0となる最小のn (default: 1)
      epsilon: 1ステップ毎の誤差上限
という使い方。

[defint_power3]

defint_power3は、積分区間を[a,b]として、 (f(x))p g(x) (pは非整数) の形の積分を行う。 ただし、f(a) = 0 とする。 要は、f(x)=xならdefint_powerと同じで、g(x)=1ならdetint_power2と同じ。

sqrt(sin(x)) * cos(x) を0から1まで積分する例を示す (sample-singular7.cc)。

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

typedef kv::interval<double> itv;

struct Func3_f {
    template <class T> T operator()(const T& x) {
        return sin(x);
    }
};

struct Func3_g {
    template <class T> T operator()(const T& x) {
        return cos(x);
    }
};

struct Func3 {
    template <class T> T operator()(const T& x) {
        return sqrt(sin(x)) * cos(x);
    }
};

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

    std::cout << kv::defint_power3_autostep(Func3(), Func3_f(), Func3_g(), (itv)0., (itv)1, 12, (itv)0.5) << "\n";
}
実行すると、
[0.51459724773239412,0.51459724773239902]
となる。

defint_power3関数は、

    defint_power3(f, g, start, end, order, p, multiplicity)
      f: f(x)pg(x)のf
      g: f(x)pg(x)のg
      start: 積分区間の左端
      end: 積分区間の右端
      order: 内部で使うベキ級数の次数
      p: ベキ数
      multiplicity: fn(a)≠0となる最小のn (default: 1)
defint_power3_autostep関数は、
    defint_power3_autostep(h, f, g, start, end, order, p, multiplicity, epsilon)
      h: 元の被積分関数 f(x)pg(x)
      f: f(x)pg(x)のf
      g: f(x)pg(x)のg
      start: 積分区間の左端
      end: 積分区間の右端
      order: 内部で使うベキ級数の次数
      p: ベキ数
      multiplicity: fn(a)≠0となる最小のn (default: 1)
      epsilon: 1ステップ毎の誤差上限
という使い方。

4.4 log特異点を端点に持つ関数の数値積分

[defint_log]

defint_logは、積分区間を[a,b]として、 log(x-a)g(x) の形の積分を行う。

log(x)/(x+1)を0から1まで積分する例を示す (sample-singular8.cc)。

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

typedef kv::interval<double> itv;

struct LogF {
    template <class T> T operator()(const T& x) {
        return log(x) / (x + 1);
    }
};

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

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

    std::cout << kv::defint_log_autostep(LogF(), LogF_s(), (itv)0., (itv)1., 12) << "\n";
}
実行すると
[-0.82246703342412398,-0.82246703342410998]
となる。

defint_log関数は、

    defint_log(g, start, end, order)
      g: log(x-start)g(x)のg
      start: 積分区間の左端
      end: 積分区間の右端
      order: 内部で使うベキ級数の次数
defint_log_autostep関数は、
    defint_log_autostep(h, g, start, end, order, epsilon)
      h: 元の被積分関数 log(x-start)g(x)
      g: log(x-start)g(x)のg
      start: 積分区間の左端
      end: 積分区間の右端
      order: 内部で使うベキ級数の次数
      epsilon: 1ステップ毎の誤差上限
という使い方。

[defint_log2]

defint_log2は、積分区間を[a,b]として、 log(f(x)) の形の積分を行う。 ただし、f(a) = 0 とする。

log(sin(x))を0から1まで積分する例を示す (sample-singular9.cc)。

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

typedef kv::interval<double> itv;

struct LogSin {
    template <class T> T operator()(const T& x) {
        return log(sin(x));
    }
};

struct LogSin_s {
    template <class T> T operator()(const T& x) {
        return sin(x);
    }
};

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

    std::cout << kv::defint_log2_autostep(LogSin(), LogSin_s(), (itv)0., (itv)1., 12) << "\n";
}
実行すると、
[-1.0567202059915891,-1.0567202059915823]
となる。

この例は、f(a)=0かつf'(a)≠0だったが、もっと高い次数まで0が続く場合も ある。n回微分で初めてf(n)(a)≠0になる場合は、 そのn (multiplicity) を指定することも出来る。

そのような例として、 log(1-cos(x))を0から1まで積分する例を示す (sample-singular10.cc)。

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

typedef kv::interval<double> itv;

struct LogCos {
    template <class T> T operator()(const T& x) {
        return log(1-cos(x));
    }
};

struct LogCos_s {
    template <class T> T operator()(const T& x) {
        return 1-cos(x);
    }
};

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

    std::cout << kv::defint_log2_autostep(LogCos(), LogCos_s(), (itv)0., (itv)1., 12, 2) << "\n";
}
実行すると、
[-2.7210654452814907,-2.7210654452814773]
となる。

defint_log2関数は、

    defint_log2(f, start, end, order, multiplicity)
      f: log(f(x))のf
      start: 積分区間の左端
      end: 積分区間の右端
      order: 内部で使うベキ級数の次数
      multiplicity: fn(a)≠0となる最小のn (default: 1)
defint_log2_autostep関数は、
    defint_log2_autostep(h, f, start, end, order, multiplicity, epsilon)
      h: 元の被積分関数 log(f(x))
      f: log(f(x))g(x)のf
      start: 積分区間の左端
      end: 積分区間の右端
      order: 内部で使うベキ級数の次数
      multiplicity: fn(a)≠0となる最小のn (default: 1)
      epsilon: 1ステップ毎の誤差上限
という使い方。

[defint_log3]

defint_log3は、積分区間を[a,b]として、 log(f(x))g(x) の形の積分を行う。 ただし、f(a) = 0 とする。 要は、f(x)=xならdefint_logと同じで、g(x)=1ならdetint_log2と同じ。

log(sin(x))cos(x)を0から1まで積分する例を示す (sample-singular11.cc)。

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

typedef kv::interval<double> itv;

struct LogSinCos {
    template <class T> T operator()(const T& x) {
        return log(sin(x)) * cos(x);
    }
};

struct LogSinCos_f {
    template <class T> T operator()(const T& x) {
        return sin(x);
    }
};

struct LogSinCos_g {
    template <class T> T operator()(const T& x) {
        return cos(x);
    }
};

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

    std::cout << kv::defint_log3_autostep(LogSinCos(), LogSinCos_f(), LogSinCos_g(), (itv)0., (itv)1., 12) << "\n";
}
実行すると、
[-0.98671202916248546,-0.98671202916247868]
となる。

defint_log3関数は、

    defint_log3(f, g, start, end, order, multiplicity)
      f: log(f(x))g(x)のf
      g: log(f(x))g(x)のg
      start: 積分区間の左端
      end: 積分区間の右端
      order: 内部で使うベキ級数の次数
      multiplicity: fn(a)≠0となる最小のn (default: 1)
defint_log3_autostep関数は、
    defint_log3_autostep(h, f, g, start, end, order, multiplicity, epsilon)
      h: 元の被積分関数 log(f(x))g(x)
      f: log(f(x))g(x)のf
      g: log(f(x))g(x)のg
      start: 積分区間の左端
      end: 積分区間の右端
      order: 内部で使うベキ級数の次数
      multiplicity: fn(a)≠0となる最小のn (default: 1)
      epsilon: 1ステップ毎の誤差上限
という使い方。

4.5 右端に特異点を持つ場合

今までの例題は全て積分区間の左端に特異点を持つ例だったが、 右端に特異点を持つ場合のために、関数名の末尾に"_r"を付けたversionを用意した。 という感じ。

sqrt(1-x2)を0から1まで積分する例を示す (sample-singular12.cc)。

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

typedef kv::interval<double> itv;

struct Quadrant {
    template <class T> T operator() (const T& x) {
        return sqrt(1 - x * x);
    }
};

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

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

    std::cout << kv::defint_power2_autostep_r(Quadrant(), Quadrant_s(), (itv)0., (itv)1., 12, (itv)0.5) << "\n";
}
実行すると、
[0.78539816339744561,0.78539816339745206]
となる。

5. 辺に特異性を持つ関数の二重積分

辺にベキ型の特異性のある関数の二重積分。 アルゴリズムは、 の第3節に解説がある。

インクルードファイルは、doubleint-singular.hpp。 下請けとしてdefint-singular.hppと、その下請けも使う。

参考ファイルは、test/test-doubleint-singular.cc。

基本的に、pが非整数で、 \[ \int_a^b \int_c^d f(x,y)^p g(x,y) dx dy \] のような形の積分を行うことができる。

最初に、 \[ \int_0^{0.125} \int_0^{0.125} \sqrt{x \cos{y}} \cos{(xy)} dx dy \] を計算する例を示す。

これは、平方根の中が辺x=0で0になっているので、辺に特異性がある。 これを計算するには、

#include <kv/doubleint-singular.hpp>

struct Func {
    template <class T> T operator() (const T& x, const T& y) {
        return sqrt(x * cos(y)) * (cos(x * y));
    }
};

struct Func_f {
    template <class T> T operator() (const T& x, const T& y) {
        return x * cos(y);
    }
};

struct Func_g {
    template <class T> T operator() (const T& x, const T& y) {
        return cos(x * y);
    }
};

typedef kv::interval<double> itv;

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

    std::cout << kv::doubleintegral_power3(Func_f(), Func_g(), itv(0.), itv(0.125), itv(0.), itv(0.125), 12, itv(0.5), 1, 0) << "\n";
}
のようなプログラムを作ればよい (sample-doublesingular1.cc)。

これを実行すると、

[0.0036779864914043189,0.003677986491404342]
と表示される。 関数doubleintegral_power3は、

doubleintegral_power3(f, g, start1, end1, start2, end2, order, power, multiplicity1, multiplicity2)

のような引数を取る。戻り値は精度保証された積分値。 それぞれの引数の意味は以下の通り。

この例だと、x cos(y)はx=0で0になっており、それは単根なので、 multiplicity1=1、yについては特異でないので、multiplicity2=0と指定している。

次は、y=0について特異な例。 \[ \int_0^{0.125} \int_0^{0.125} \sqrt{(\cos{x}) y} \cos{(xy)} dx dy \] を計算する例を示す。

プログラムは以下の通り。 (sample-doublesingular2.cc)。

#include <kv/doubleint-singular.hpp>

struct Func {
    template <class T> T operator() (const T& x, const T& y) {
        return sqrt(cos(x) * y) * (cos(x * y));
    }
};

struct Func_f {
    template <class T> T operator() (const T& x, const T& y) {
        return cos(x) * y;
    }
};

struct Func_g {
    template <class T> T operator() (const T& x, const T& y) {
        return cos(x * y);
    }
};

typedef kv::interval<double> itv;

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

    std::cout << kv::doubleintegral_power3(Func_f(), Func_g(), itv(0.), itv(0.125), itv(0.), itv(0.125), 12, itv(0.5), 0, 1) << "\n";
}
実行すると、
[0.0036779864914043176,0.0036779864914043438]
と表示される。multiplicityの指定が0,1となっている以外はほとんど前と同じ。

次は、x=0、y=0両方の辺について特異な例。 \[ \int_0^{0.125} \int_0^{0.125} \sqrt{x y} \cos{(xy)} dx dy \] を計算する例を示す。

プログラムは以下の通り。 (sample-doublesingular3.cc)。

#include <kv/doubleint-singular.hpp>

struct Func {
    template <class T> T operator() (const T& x, const T& y) {
        return sqrt(x * y) * (cos(x * y));
    }
};

struct Func_f {
    template <class T> T operator() (const T& x, const T& y) {
        return x * y;
    }
};

struct Func_g {
    template <class T> T operator() (const T& x, const T& y) {
        return cos(x * y);
    }
};

typedef kv::interval<double> itv;

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

    std::cout << kv::doubleintegral_power3(Func_f(), Func_g(), itv(0.), itv(0.125), itv(0.), itv(0.125), 12, itv(0.5), 1, 1) << "\n";
}
これを実行すると、次のように表示される。
[0.00086803609297475211,0.00086803609297475852]
x,y両方について特異なので、multiplicityの指定が1,1となっている。

次は、多重度が2の例。 \[ \int_0^{0.125} \int_0^{0.125} ((1-\cos{x}) \cos{y})^{1/3} \cos{xy} dx dy \] を計算する例を示す。

プログラムは以下の通り。 (sample-doublesingular4.cc)。

#include <kv/doubleint-singular.hpp>

struct Func {
    template <class T> T operator() (const T& x, const T& y) {
        return pow((1-cos(x)) * cos(y), 1/T(3)) * (cos(x * y));
    }
};

struct Func_f {
    template <class T> T operator() (const T& x, const T& y) {
        return (1-cos(x)) * cos(y);
    }
};

struct Func_g {
    template <class T> T operator() (const T& x, const T& y) {
        return cos(x * y);
    }
};

typedef kv::interval<double> itv;

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

    std::cout << kv::doubleintegral_power3(Func_f(), Func_g(), itv(0.), itv(0.125), itv(0.), itv(0.125), 12, 1./itv(3.), 2, 0) << "\n";
}
これを実行すると、次のように表示される。
[0.0018582185546727903,0.0018582185546728117]
1-cos(x)はx=0で多重度2の解になっているので、multiplicityは2,0と 指定している。

最後に、逆側の辺で特異になる例。 \[ \int_0^{0.125} \int_0^{0.125} \sqrt{\sin{0.125-x}\cos{y}} \cos{xy} dx dy \] を計算する例を示す。

プログラムは以下の通り。 (sample-doublesingular5.cc)。

#include <kv/doubleintegral.hpp>
#include <kv/doubleint-singular.hpp>

struct Func {
    template <class T> T operator() (const T& x, const T& y) {
        return sqrt(sin(0.125-x) * cos(y)) * (cos(x * y));
    }
};

struct Func_f {
    template <class T> T operator() (const T& x, const T& y) {
        return sin(0.125-x) * cos(y);
    }
};

struct Func_g {
    template <class T> T operator() (const T& x, const T& y) {
        return cos(x * y);
    }
};

typedef kv::interval<double> itv;

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

    std::cout << kv::doubleintegral_power3_r1(Func_f(), Func_g(), itv(0.), itv(0.125), itv(0.), itv(0.125), 12, itv(0.5), 1, 0) << "\n";
}
これを実行すると、次のように表示される。
[0.0036759640850438091,0.0036759640850438352]
この例はxについて特異だが反対側の辺x=0.125で特異になっている。 その場合は、doubleintegral_power3_r1を呼び出す。 r1は、start1=0でなくend1=0で特異と見なして計算する。 内部では関数を反転してdoubleintegral_power3を呼び出しているだけ。 同様に、doubleintegral_power3_r2は、第2変数について反転し、 end2=0で特異の場合の計算を行う。 また、doubleintegral_power3_r12は、第1,2変数両方について 反転し、end1=0, end2=0の両方で特異の場合の計算を行う。

6. 一点で特異性を持つ関数の二重積分

一点にベキ型の特異性のある関数の二重積分。 アルゴリズムは、 の第4節に解説がある。

インクルードファイルは、doubleint-singular.hpp。 下請けとしてdefint-singular.hppと、その下請けも使う。

参考ファイルは、test/test-doubleint-singular.cc。

pdfファイル中に説明があるように、一点でベキ型の特異性がある関数については、 まず三角形の1つの頂点で特異になるように領域を分割する必要がある。 ここでは、1つの頂点で特異になるような簡単な例を示す。

点(0,0)、点(0.1,0)、点(0.1,0.1)を頂点とする三角形上の積分で、 原点(0,0)が特異となっている、 \[ \int_0^{0.1} \int_0^x \sqrt{x + y} \, dy dx \] を計算することを考える。

プログラムは以下の通り。 (sample-doublesingular6.cc)。

#include <kv/doubleint-singular.hpp>

struct Func {
    template <class T> T operator() (const T& x, const T& y) {
        return sqrt(x + y);
    }
};

struct Func_f {
    template <class T> T operator() (const T& x, const T& y) {
        return x + y;
    }
};

struct Func_g {
    template <class T> T operator() (const T& x, const T& y) {
        return T(1);
    }
};

typedef kv::interval<double> itv;

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

    std::cout << kv::doubleintegral_singular_point(Func_f(), Func_g(), itv(0.), itv(0.), itv("0.1"), itv(0.), itv("0.1"), itv("0.1"), 12, itv(0.5), 1) << "\n";;
}
これを実行すると、次のように表示される。
[0.0015418651332575329,0.0015418651333636226]
このように、doubleintegral_singular_pointを呼び出せばよい。

関数doubleintegral_singular_pointは、

doubleintegral_singular_point(f, g, x0, y0, x1, y1, x2, y2, order, power, multiplicity, div)

のような引数を取る。戻り値は精度保証された積分値。 それぞれの引数の意味は以下の通り。

7. その他

Type-II PSAの応用例の一つとしてとりあえず作ってみたもので、 数値積分の精度保証のベストな方法とは言えないかもしれない。 従来型の数値積分公式を用い、誤差公式の定数を評価する部分のみに 高階自動微分 (PSA) を使うのが王道かも。 しかし、それだと誤差定数をグローバルに評価することになるので、 変化の激しい関数を長い区間に渡って積分する場合は、本手法も 有力と言えそう。

2次元の方は全くの試作品で、速度、精度ともに不十分。 そもそも2次元の積分の精度保証を行っている例そのものが少ないのではないかと思い、 一応作ってみた。内部では、psa型の係数にpsa型を持たせたものを用いている。

端点特異性への対応は、分類が多すぎて面倒なのでもう少し自動化出来ないかと 考えているがなかなか難しい。log(f(x))g(x)ph(x)の形は 現状では出来ないので作りたいと考えているが、なんというか切りがないな。

2次元の特異積分は、応用例はたくさんありそうだが、現状は遅すぎて なかなか使ってもらえない気がする。