振動積分のためのDE公式の変数変換の計算

はじめに

本記事では振動積分 \[\begin{aligned} \int_0^\infty \sin(x) f(x) \,\mathrm{d}{x}, \qquad \int_0^\infty \cos(x) f(x) \,\mathrm{d}{x} \end{aligned}\] のためのDE公式(Ooura and Mori 1999)で用いられる変数変換とその正弦・余弦・微分を計算するときの注意事項を示す．

具体的には, \(\int_0^\infty \sin(x)f(x)\,\mathrm{d}{x}\)の計算では \[\begin{aligned} x_{\sin,h}(t) = \frac{\pi}{h} \phi(t), \qquad \phi(t) :=\frac{t}{1 - \mathrm{e}^{u(t)}}, \qquad u(t) :=-2t - \alpha(1-\mathrm{e}^{-t}) - \beta (\mathrm{e}^t-1) \end{aligned}\] という変換を用いる．なお，\(x_{\sin,h}(t)\)の微分は \[\begin{aligned} x_{\sin,h}'(t) = \frac{\pi}{h} \phi'(t) = \frac{\pi}{h} \left[ \frac{1}{1 - \mathrm{e}^{u(t)}} - \frac{t u'(t) \mathrm{e}^{u(t)}}{(1-\mathrm{e}^{u(t)})^2} \right] = \frac{\pi}{h} \frac{1 -\mathrm{e}^{u(t)} + tu'(t)\mathrm{e}^{u(t)}}{(1-\mathrm{e}^{u(t)})^2} \end{aligned}\] である． また，\(\int_0^\infty \cos(x)f(x)\,\mathrm{d}{x}\)の計算では \[\begin{aligned} x_{\cos,h}(t) = \frac{\pi}{h}\phi\left(t - \frac{h}{2}\right) \end{aligned}\] という変換を用いる．なお，\(x_{\cos,h}(t)\)の微分は\(x_{\cos,h}'(t) = \pi\phi'(t-h/2)/h\)である． 本記事では\(t = kh\; (k\in\mathbb{Z})\)に対して計算する際の注意事項を示す．

\(\phi(t)\)の計算

\(t\approx 0\)での計算

\(\phi(t)\)の\(t=0\)におけるTaylor展開は \[\begin{aligned} & \phi(t) = c_0 + c_1t + \frac{c_{2\mathrm{n}}}{c_{2\mathrm{d}}}t^2 + \mathcal{O}(t^3),\\ & \quad c_0 = \frac{1}{2 +\alpha + \beta},\\ & \quad c_1 = \frac{\alpha^2 + 2 \alpha\beta + 5 \alpha +\beta^2 + 3\beta + 4}{2(2 +\alpha + \beta)^2},\\ & \quad c_{2\mathrm{n}} = \alpha^{4} + 4 \alpha^{3} \beta + 8 \alpha^{3} + 6 \alpha^{2} \beta^{2} + 24 \alpha^{2} \beta + 25 \alpha^{2} + 4 \alpha \beta^{3} + 24 \alpha \beta^{2}\\ & \qquad + 38 \alpha \beta + 28 \alpha + \beta^{4} + 8 \beta^{3} + 25 \beta^{2} + 28 \beta + 16\\ & \quad c_{2\mathrm{d}} = 12(\alpha+\beta+2)^3 \end{aligned}\] であるので，\(|t| \ll 1\)であるときは \[\begin{aligned} \phi(t) \approx c_0 + c_1 t + c_2 t^2 \end{aligned}\] と計算するのが良いだろう．

\(\phi'(t)\)の計算

\(t\approx 0\)での計算

\(\phi(t)\)の\(t=0\)におけるTaylor展開は \[\begin{aligned} & \frac{h}{\pi} x_h'(t) = c_0 + \frac{c_{1\mathrm{n}}}{c_{1\mathrm{d}}}t + \frac{c_{2\mathrm{n}}}{c_{2\mathrm{d}}}t^2 + \mathcal{O}(t^3),\\ % & \quad c_0 = \frac{\alpha^{2} + 2 \alpha \beta + 5 \alpha + \beta^{2} + 3 \beta + 4}{2 (\alpha^{2} + 2 \alpha \beta + 4 \alpha + \beta^{2} + 4 \beta + 4)},\\ % & \quad c_{1\mathrm{n}} = \alpha^{4} + 4 \alpha^{3} \beta + 8 \alpha^{3} + 6 \alpha^{2} \beta^{2} + 24 \alpha^{2} \beta + 25 \alpha^{2} + 4 \alpha \beta^{3} + 24 \alpha \beta^{2} + 38 \alpha \beta + 28 \alpha\\ & \qquad + \beta^{4} + 8 \beta^{3} + 25 \beta^{2} + 28 \beta + 16,\\ % & \quad c_{1\mathrm{d}} = 6 \alpha^{3} + 18 \alpha^{2} \beta + 36 \alpha^{2} + 18 \alpha \beta^{2} + 72 \alpha \beta + 72 \alpha + 6 \beta^{3} + 36 \beta^{2} + 72 \beta + 48,\\ % & \quad c_{2\mathrm{n}} = - \alpha^{5} - 3 \alpha^{4} \beta - 8 \alpha^{4} - 2 \alpha^{3} \beta^{2} - 16 \alpha^{3} \beta - 24 \alpha^{3} + 2 \alpha^{2} \beta^{3} - 36 \alpha^{2} \beta - 36 \alpha^{2}\\ & \qquad + 3 \alpha \beta^{4} + 16 \alpha \beta^{3} + 36 \alpha \beta^{2} - 12 \alpha + \beta^{5} + 8 \beta^{4} + 24 \beta^{3} + 36 \beta^{2} + 12 \beta,\\ % & \quad c_{2\mathrm{d}} = 8 \alpha^{4} + 32 \alpha^{3} \beta + 64 \alpha^{3} + 48 \alpha^{2} \beta^{2} + 192 \alpha^{2} \beta + 192 \alpha^{2} + 32 \alpha \beta^{3} + 192 \alpha \beta^{2}\\ & \qquad + 384 \alpha \beta + 256 \alpha + 8 \beta^{4} + 64 \beta^{3} + 192 \beta^{2} + 256 \beta + 128 \end{aligned}\] であるので，\(|t|\ll 1\)であるときは \[\begin{aligned} \phi'(t) = c_0 + \frac{c_{1\mathrm{n}}}{c_{1\mathrm{d}}}t + \frac{c_{2\mathrm{n}}}{c_{2\mathrm{d}}}t^2 \end{aligned}\] と計算するのが良いだろう．

\(t<0\)での計算

\(t<0\)のとき\(u(t)>0\)となる．このため倍精度では\(t < -8\)において\(\mathrm{e}^{u(t)}=\texttt{Inf}\)が起こり得る． そこで，分母分子を\(\mathrm{e}^{u(t)}\)で割れば，そのような状況下での計算結果を強制的0に持っていくことができる． \[\begin{aligned} \phi'(t) = \frac{1 -\mathrm{e}^{u(t)} + tu'(t)\mathrm{e}^{u(t)}}{(1-\mathrm{e}^{u(t)})^2} = \frac{\mathrm{e}^{-u(t)}- 1 + tu'(t)}{\mathrm{e}^{-u(t)}(1-\mathrm{e}^{u(t)})^2} = \frac{\mathrm{e}^{-u(t)}- 1 + tu'(t)}{(2\sinh(u(t)/2))^2} \end{aligned}\]

\(\sin(x_{\sin,h}(kh))\;(k\in\mathbb{Z})\)の計算

\(kh \approx 0\)のとき

\(x_{\sin,h}(kh) \approx 0\)なので，\(x_{\sin,h}(kh)\)を\(\phi(t)\)のTaylor展開で計算して，素直に\(\sin(x_{\sin,h}(kh))\)を計算するのが良いであろう．

\(k > 0\)のとき

\(k>0\)のとき，\(u(kh) < 0\)なので \[\begin{aligned} \sin(x_{\sin,h}(kh)) = \sin\left(\frac{k \pi}{1 - \mathrm{e}^{u(kh)}}\right) = \sin\left(k\pi \left(1 + \frac{\mathrm{e}^{u(kh)}}{1 - \mathrm{e}^{u(kh)}}\right)\right) = (-1)^k \sin\left(\frac{k\pi \mathrm{e}^{u(kh)}}{1 - \mathrm{e}^{u(kh)}}\right) \end{aligned}\] と計算するのが良いだろう．

\(k < 0\)のとき

\(k<0\)のとき，\(u(kh) > 0\)なので \[\begin{aligned} \sin(x_{\sin,h}(kh)) = \sin\left(\frac{k \pi}{1 - \mathrm{e}^{u(kh)}}\right) \end{aligned}\] と計算するのが良いだろう．

\(\cos(x_{\cos,h}(kh))\;(k\in\mathbb{Z})\)の計算

\((k-1/2)h \approx 0\)のとき

\(x_{\cos,h}(kh)\)を\(\phi((k-1/2)h)\)のTaylor展開で計算して，素直に\(\cosh(\pi x_{\cos,h}(kh)/h)\)を計算するのが良いであろう．

\(k > 0\)のとき

\(k>0\)のとき，\(u(kh) < 0\)なので \[\begin{aligned} \cos(x_{\cos,h}(kh)) = \cos\left(\frac{(k-1/2) \pi}{1 - \mathrm{e}^{u(kh)}}\right) = \cos\left(\left(k+\frac{1}{2}\right)\pi \left(1 + \frac{\mathrm{e}^{u(kh)}}{1 - \mathrm{e}^{u(kh)}}\right)\right) = (-1)^k \sin\left(\frac{(k-1/2)\pi \mathrm{e}^{u(kh)}}{1 - \mathrm{e}^{u(kh)}}\right) \end{aligned}\] と計算するのが良いだろう．

\(k < 0\)のとき

\(k<0\)のとき，\(u(kh) > 0\)なので \[\begin{aligned} \cos(x_{\cos,h}(kh)) = \cos\left(\frac{(k-1/2)\pi}{1 - \mathrm{e}^{u(kh)}}\right) \end{aligned}\] と計算するのが良いだろう．

Ooura, Takuya, and Masatake Mori. 1999. “A Robust Double Exponential Formula for Fourier-Type Integrals.” Journal of Computational and Applied Mathematics 112 (1): 229–41. https://doi.org/https://doi.org/10.1016/S0377-0427(99)00223-X.