Fourier解析は科学技術研究で極めて重要です．

ところが，そのFourier解析におけるFourier変換

$$\mathscr{F}[f](\xi)=\int_{-\infty}^{\infty}f(x)\mathrm{e}^{-2\pi\mathrm{i}\xi x}\mathrm{d}x$$

は従来の数値積分公式で計算するのが難しいです．そこで，Fourier変換に対し佐藤超函数論を応用した数値計算法を考案しました[2]．

計算のカギは，Fourier変換を次のように表すことです．

$$\mathscr{F}[f](\xi)=\lim_{\epsilon\rightarrow 0+}\int_{-\infty}^{\infty}f(x)\mathrm{e}^{-2\pi\mathrm{i}\xi x-2\pi\epsilon|x|}\mathrm{d}x.$$

これは次のように書き換えられます．

$$\mathscr{F}[f](\xi)=\lim_{\epsilon\rightarrow 0+}\left\{\mathfrak{F}_+[f](x+\mathrm{i}\epsilon)-\mathfrak{F}_-[f](x-\mathrm{i}\epsilon)\right\},$$

$$\mathfrak{F}_{\pm}[f](\zeta) = \pm\int_0^{\infty}f(\mp x)\mathrm{e}^{\pm 2\pi\mathrm{i}\zeta x}\mathrm{d}x,\quad (\:\pm\mathrm{Im}\zeta>0\:).$$

ここに現れる関数\(\mathfrak{F}_{\pm}[f](\zeta)\)は上／下半平面\(\pm\mathrm{Im}\zeta>0\)で解析関数であり，その実軸上における境界値の差としてFourier変換が表されます．実はこれが，佐藤超函数論におけるFourier変換の定義です．我々の方法は，このFourier変換の定義に従ってFourier変換を数値計算しているのです．すなわち，解析関数\(\mathfrak{F}_{\pm}[f](\zeta)\)を数値的に求め，それらの境界値の差$$\mathfrak{F}_+[f](\xi+\mathrm{i}0)-\mathfrak{F}_-[f](\xi-\mathrm{i}0)$$を数値計算しているのです．

具体的な数値計算法は次のとおりです．

1. 解析関数\(\mathfrak{F}_{\pm}[f](\zeta)\)をTaylor級数の形$$\mathfrak{F}_{\pm}[f](\zeta)=\sum_{n=0}^{\infty}c_n^{(\pm)}(\zeta-\zeta_0^{(\pm)})^n$$で求めます．ここで，\(\zeta_0^{(\pm)} \ ( \:\pm\mathrm{Im}\zeta_0^{(\pm)}>0\:)\)は上／下半平面に適当に取った点です．\(c_n^{(\pm)}\)は\(\mathfrak{F}_{\pm}[f](\zeta)\)の\(\zeta_0^{(\pm)}\)における微分係数\((\div n!)\)であり，次で与えられます．$$c_n^{(\pm)}=\pm\frac{1}{n!}\int_0^{\infty}(\pm 2\pi\mathrm{i}x)^n f(\mp x)\mathrm{e}^{\pm 2\pi\mathrm{i}\zeta_0^{(\pm)}x}\mathrm{d}x \quad ( \: n = 0, 1, 2, \ldots \: ).$$ ここで，右辺の被積分関数には指数関数的減衰する因子\(\mathrm{e}^{-2\pi|\mathrm{Im}\zeta_0^{(\pm)}|x}\)が含まれているので，従来使われている数値積分公式で容易に計算することができます．
2. 1.で求めた解析関数\(\mathfrak{F}_{\pm}[f](\zeta)\)の実軸上での値\(\mathfrak{F}_{\pm}[f](\xi\pm\mathrm{i}0)\)が求まれば，Fourier変換が得られますが，Taylor級数のままだと収束円の範囲でしか\(\mathfrak{F}_{\pm}[f](\zeta)\)が計算できないので，これらを解析接続する必要があります．そのために，これらを連分数の形に変換します．$$\mathfrak{F}_{\pm}[f](\zeta)=\frac{a_0^{(\pm)}}{\displaystyle 1+\frac{\displaystyle a_1^{(\pm)}(\zeta-\zeta_0^{(\pm)})}{\displaystyle 1+\frac{a_2^{(\pm)}(\zeta-\zeta_0^{(\pm)})}{\displaystyle 1+\ddots}}}.$$ ここで，連分数の係数\(a_n^{(\pm)}\)は「商差法」というアルゴリズムを用いて，Taylor係数\(c_n^{(\pm)}\)から計算できます[1]．関数項連分数にすると一般にTaylor級数に比べて収束域が広くなるので，これにより\(\mathfrak{F}_{\pm}[f](\zeta)\)の解析接続が数値的に行なえます．そして，連分数は収束が速いので効率的に数値計算ができます．

（数値例）次の3つの関数\(f(x)\)に対し我々の方法でFourier変換を計算しました．

$$\mathrm{(1)}\quad \mathscr{F}[(1+x^2)^{-1}](\xi)=\pi\mathrm{e}^{-2\pi|\xi|},$$

$$\mathrm{(2)}\quad \mathscr{F}\mathscr[\tanh(\pi x)](\xi)=-\mathrm{i}\mathrm{cosech}(\pi\xi),$$

$$\mathrm{(3)}\quad \mathscr{F}\mathscr[\log|x|](\xi)=-\gamma-\frac{1}{2|\xi|}.$$

これら3つのFourier変換\(\mathscr{F}[f](\xi)\)の\(\xi=1\)における値を数値計算し，相対誤差および\(f(x)\)の計算回数を示したのが下の表です（[2]より引用）．ここでは，既存の杉原による方法[4]，大浦・森による方法[3]と比較しています．既存の2方法と比べて，我々の超函数論に基づく方法では効率的・高い精度でFourier変換が計算できていることがわかります．

\(f(x)\)		(1)	(2)	(3)
超函数論の方法[2]	相対誤差	\(2.5\times 10^{-39}\)	\(4.5\times 10^{-76}\)	\(8.9\times 10^{-39}\)
	\(f(x)\)の計算回数	\(1242\)	\(1202\)	\(1252\)
杉原の方法[4]	相対誤差	\(1.4\times 10^{-19}\)	\(9.1\times 10^{-20}\)	\(4.1\times 10^{-20}\)
	\(f(x)\)の計算回数	\(27072\)	\(53114\)	\(35280\)
大浦・森の方法[3]	相対誤差	\(6.7\times 10^{-63}\)	\(4.5\times 10^{-57}\)	\(1.2\times 10^{-69}\)
	\(f(x)\)の計算回数	\(1690\)	\(1618\)	\(1698\)

 

佐藤超函数論に基づく数値計算法のココロは，数値積分の場合と同様，実数計算の世界では計算が難しいものを，複素数計算の世界の計算が容易なところで計算し，その結果を実数計算の世界に反映させるところにあります．Fourier変換の場合，実数の世界で\(\mathscr{F}[f](\xi)\)を直接計算しようとするのは困難ですが，複素数の世界で解析関数\(\mathfrak{F}_{\pm}[f](\zeta)\)を計算するのは容易である．だから，複素数計算で\(\mathfrak{F}_{\pm}[f](\zeta)\)を計算し，それを実数の世界に解析接続して，$$\mathscr{F}[f](\xi)=\mathfrak{F}_+[f](\xi+\mathrm{i}0)-\mathfrak{F}_-[f](\xi-\mathrm{i}0)$$により所望のFourier変換を得るのです．

複素解析関数の特質は，Cauchyの積分定理などからわかるように，遠く離れた2点での値が密接に結びついているということです．そういう特質があるので，複素数の世界の計算しやすいところで計算して実数の世界に戻すという手法が可能になっています．その意味で，佐藤超函数論に基づく数値計算法は，2点の値が密接に結びついているという複素解析関数の特質を最大限に生かした方法と言えるでしょう．

参考文献

[1] P. Henrici: Applied and Computational Complex Analysis, Vol. 2, John Wiley & Sons, New York (1977).

[2] H. Ogata: Numerical calculation of Fourier transforms based on hyperfunction theory, J. Comput. Appl. Math., 378 (2020) 112924.

[3] T. Ooura and M. Mori: A robust double exponential formula for oscillatory functions over the half infinite interval, J. Comput. Appl. Math., 38 (1991) 353–360.

[4] M. Sugihara: Methods of numerical integration of oscillatory functions by the DE-formula with the Richardson extrapolation, J. Comput. Appl. Math., 17 (1987) 47–68.