最終更新: 2021/1/27

端点に特異性のある1次元の精度保証付き数値積分

柏木雅英

1. はじめに

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

端点特異性を持つ関数の精度保証付き数値積分 (jsiam2015での予稿 - 誤字修正あり)
端点特異性を持つ関数の精度保証付き数値積分 (jsiam2015でのスライド - 誤字修正あり)

のようなものを使用している。

現在のところ、

f(x)/g(x) (f(0)=0, g(0)=0) のような、見かけ上の特異点を含む関数
f(x)^p g(x) (f(0)=0, pは非整数) の形の関数
log(f(x))g(x) (f(0)=0) の形の関数

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

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

2. ファイル構成

defint-singular.hpp

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

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

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をtⁿで割る関数。 div_reduce(p1, p2, n)は、p1, p2をベキ級数型、nを自然数として、 p1, p2をともにtⁿで割ってから除算を行う関数。よって、

div_tn(sin(x), 1)
div_reduce(sin(x), x, 1)

のどちらでもよい。

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版のみを例に説明する。

以下は、x²/(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. ベキ型特異点を端点に持つ関数の数値積分

[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]

となる。

このように、

被積分関数(x-a)^pg(x)
被積分関数から(x-a)^pを除いたg(x)のみの関数

の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⁽ⁱ⁾(a)=0 (i < n)、f⁽ⁿ⁾(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: fⁿ(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: fⁿ(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: fⁿ(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: fⁿ(a)≠0となる最小のn (default: 1)
      epsilon: 1ステップ毎の誤差上限

という使い方。

5. 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⁽ⁿ⁾(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: fⁿ(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: fⁿ(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: fⁿ(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: fⁿ(a)≠0となる最小のn (default: 1)
      epsilon: 1ステップ毎の誤差上限

という使い方。

6. 右端に特異点を持つ場合

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

defint_power3に対してdefint_power3_r、
defint_power3_autostepに対してdefint_power3_autostep_r

という感じ。

sqrt(1-x²)を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]

となる。

7. おわりに

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