◆ _deoint_r()

半無限区間の積分 (フーリエ型積分) (自動積分) (二重指数関数(DE)公式) (リバースコミュニケーション版)

目的: 本ルーチンは二重指数関数(DE)公式を使用した自動積分によりフーリエ型積分を求める.

Ic = ∫ f(x) cosωx dx [a, +∞] または Is = ∫ f(x) sinωx dx [a, +∞]

ここで, f(x)は変数irevに従ってユーザーが必要な関数値を計算して与える.

DE公式では次の変換関数を使用することにより無限区間[-∞, +∞]の積分に変換して台形公式により積分値を計算する.
Ic: x = Mφ(t - π/2M)/ω

Is: x = Mφ(t)/ω,

φ(t)としては
φ1(t) = t/(1 - exp(-2t -α(1 - exp(-t)) - β(exp(t) - 1))),

β = 1/4, α = β/sqrt(1 + Mlog(1 + M)/4π)

または

φ2(t) = t/(1 - exp(-ksinh(t))), k = 2π

を使用する. 刻み幅hと定数Mは, Mh = π が成り立つように選ばれる.
詳細は下記文献を参照のこと.

引数

[in]	a	積分の下限値.
[in]	omega	cosまたはsinの中のωの値.
[in]	iflag	積分の種類を選択する. = 0: ∫ f(x) cosωx dx [0, +∞] = 1: ∫ f(x) sinωx dx [0, +∞] また, 次のようにして変換関数を選択することができる. iflag += 0: φ1(t) (推奨) iflag += 2: φ2(t)
[in]	eps	要求相対誤差. (0 < eps <= 1)
[out]	result	求められた積分値.
[out]	neval	f(x)の評価回数.
[out]	err	vの推定相対誤差. 多くの場合大きめに算出される. また, err > eps であっても正常終了として戻ることがある.
[out]	info	= 0: 正常終了 = -3: 入力パラメータ iflag の誤り (iflag = 0〜3 以外) = -4: 入力パラメータ eps の誤り (eps <= 0 または eps > 1) = 1: 収束しなかった
[out]	xx	irev = 1, 2: 再呼び出し時に与える関数値を求めるべき点を返す.
[in]	yy	irev = 1, 2: 再呼び出し時に関数値f(xx)を与えること.
[in,out]	irev	リバースコミュニケーションの制御変数. [in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない. [out] 0 以外の場合, 下記処理を行いirevを変更せずに再び本ルーチンを呼び出すこと. = 0: 処理終了. = 1, 2: xxにおける関数値を求めyyに設定する. yy以外の変数を変更してはならない.

参考文献

Ooura and Mori, "The double exponential formula for oscillatory functions over the half infinite interval", J. Comput. Appl. Math., 38 (1991), 353-360.
Ooura and Mori, "A robust double exponential formula for Fourier-type integrals", J. Comput. Appl. Math., 112 (1999), 229-241.