インクルードファイルは、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(f, g, start1, end1, start2, end2, order, power, multiplicity1, multiplicity2)
のような引数を取る。戻り値は精度保証された積分値。 それぞれの引数の意味は以下の通り。
次は、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の両方で特異の場合の計算を行う。
インクルードファイルは、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)
のような引数を取る。戻り値は精度保証された積分値。 それぞれの引数の意味は以下の通り。