A numerical method of computing oscillatory integral related to hyperfunction theory

1 Introduction

In this paper, we consider the computation of an integral

I = \int_{0}^{\infty} f (x) d x,

(1)

where $f (x)$ is a slowly decaying oscillatory function, such as

\int_{0}^{\infty} \frac{x cos x}{1 + x^{2}} d x, \int_{0}^{\infty} J_{ν} (x) d x .

It is difficult to compute this integral by conventional methods such as the DE formula [8]. In this paper, we propose a numerical method for computing integrals of this type. In the proposed method, we consider an integral of the Fourier-Laplace transform type

F [f] (ζ) = \int_{0}^{\infty} f (x) e^{i ζ x} d x,

(2)

which is a complex analytic function in the upper-half complex plane $Im ζ > 0$ under some conditions on $f (x)$ . We obtain the function $F [f] (ζ)$ in the upper-half plane, and we compute the analytic continuation of $F [f] (ζ)$ on the real axis $R$ . Then, we obtain the desired integral by

I = \int_{0}^{\infty} f (x) d x = lim ϵ \to 0 + F [f] (i ϵ) .

(3)

Our method is related to hyperfunction theory, a theory of generalized functions proposed by Sato [7]. Roughly speaking, a hyperfunction is given by the boundary values of functions which are analytic in subdomains of $C ∖ R$ . From (3), the desired integral $I$ is a boundary value on $R$ of the function $F [f] (ζ)$ analytic in the upper half plane $Im ζ > 0$ , and, in this sense, the integral $I$ can be regarded as a value of a hyperfunction.

Previous studies related to this paper are as follows. Toda and Ono proposed a method of computing (1) by the limit

I = lim n \to \infty \int_{0}^{\infty} f (x) exp (- 2^{- n} x) d x,

which is obtained by the DE formula and the Richardson extrapolation [9]. Sugihara improved Toda and Ono’s method and proposed a method of computing (1) by the limit

I = lim n \to \infty \int_{0}^{\infty} f (x) exp (- 2^{- n} x^{2}) d x,

which is obtained by the DE formula and the Richardson extrapolation. Ooura and Mori proposed a DE-type formula designed especially for Fourier transform type integrals using a unique DE transform

[5] and improved it for integrands with singularities in the complex plane [6]. The author proposed a method of computing Fourier transform in a way similar to that of this paper [4].

The remainder of this paper is structured as follows. In Section 2, we propose a numerical method for computing the oscillatory integral (1). In Section 3, we give a brief review of hyperfunction theory and mention a relation between the proposed method and hyperfunction theory. In Section 5, we give the summary and concluding remarks of this paper and refer to problems for future study.

2 Numerical method

We consider an integral (1), where $f (x)$ is a slowly decaying oscillatory function and propose a numerical method for computing it. As mentioned in the previous section,

we obtain a function $F [f] (ζ)$ , which is a complex analytic function in the upper half plane $Im ζ > 0$ under some conditions on $f (x)$ ,
and we obtain the desired integral $I$ by the analytic continuation of $F [f] (ζ)$ onto the real axis.

The details of each step are as follows.

(Computation of $F [f] (ζ)$ in the upper half plane $Im ζ > 0$ ) We obtain the analytic function $F [f] (ζ)$ as a Taylor series

$F [f] (ζ) = \infty \sum n = 0 c_{n} (ζ - ζ_{0})^{n},$

where $ζ_{0}$ is a complex constant given by the user such that $Im ζ_{0} > 0$ . The Taylor coefficient $c_{n}$ is given by the integral

$c_{n} = \frac{1}{n!} {\frac{d^{n}}{d ζ^{n}} F [f] (ζ) ∣ ∣ ∣}_{ζ = ζ_{0}} = \frac{1}{n!} \int_{0}^{\infty} (i x)^{n} f (x) e^{i ζ_{0} x} d x,$ (4)

This integral involves an exponentially decaying factor $e^{i ζ_{0} x}$ , where we remark that $Im ζ_{0} > 0$ , and it is easy to compute it by conventional numerical integration formulas.
(Analytic continuation of $F [f] (ζ)$ ) We obtain the analytic continuation of $F [f] (ζ)$ onto the real axis $R$ using its continued fraction expansion since the continued fraction expansion of an analytic function has a larger convergence region compared with that of its Taylor series. We can convert a Taylor series

$F [f] (ζ) = \infty \sum n = 0 c_{n} (ζ - ζ_{0})^{n}$

into a continued fraction

$F [f] (ζ) = \frac{a_{0}}{1 + \frac{a_{1} (ζ - ζ_{0})}{1 + \frac{a_{2} (ζ - ζ_{0})}{1 + ⋱}}}$ (5)

using the so-called quotient difference (QD) algorithm (see Theorem 12.4c of [3]). Namely, we generate the sequences $e_{k}^{(n)}$ $(n = 0, 1, 2, \dots; k = 0, 1, 2, \dots)$ and $q_{k}^{(n)}$ $(n = 0, 1, 2, \dots; k = 1, 2, \dots)$ by the recurrence relations

$e_{0}^{(n)} = 0, q_{1}^{(n)} = \frac{c_{n + 1}}{c_{n}} (n = 0, 1, 2, \dots),$

$\begin{matrix} e_{k}^{(n)} = q_{k}^{(n + 1)} - q_{k}^{(n)} + e_{k - 1}^{(n + 1)}, q_{k + 1}^{(n)} = \frac{e_{k}^{(n + 1)}}{e_{k}^{(n)}} q_{k}^{(n + 1)} (n = 0, 1, 2, \dots; k = 1, 2, \dots), \end{matrix}$

and, then, we obtain

$F [f] (ζ) = \frac{c_{0}}{1 - \frac{q_{1}^{(0)} (ζ - ζ_{0})}{1 - \frac{e_{1}^{(0)} (ζ - ζ_{0})}{1 - \frac{q_{2}^{(0)} (ζ - ζ_{0})}{1 - \frac{e_{2}^{(0)} (ζ - ζ_{0})}{1 - ⋱}}}}} .$

Once we obtain the coefficients $a_{n}$ of the continued fraction (5), we can compute the continued fraction by

$F [f] (ζ) = lim k \to \infty \frac{P_{k} (ζ)}{Q_{k} (ζ)},$

where $P_{k} (ζ)$ and $Q_{k} (ζ)$ are the polynomials given by the recurrence relations

$P_{- 1} = 0, Q_{- 1} = 1, P_{0} = c_{0}, Q_{0} = 1,$

${\begin{matrix} P_{k} (ζ) = a_{k} (ζ - ζ_{0}) P_{k - 2} (ζ) + P_{k - 1} (ζ) Q_{k} (ζ) = a_{k} (ζ - ζ_{0}) Q_{k - 2} (ζ) + Q_{k - 1} (ζ) \end{matrix} (k = 1, 2, \dots) .$

In the procedure of the QD algorithm, we use multiple precision arithmetic since the algorithm is numerically unstable.

Finally, we obtain the desired oscillatory integral by

I = \int_{0}^{\infty} f (x) d x = F [f] (0) = lim ϵ \to 0 + F [f] (i ϵ) .

(6)

3 Relation with hyperfunction theory

In this section, we give a brief review of hyperfunction theory and its relation with the proposed method. For the detail of hyperfunction theory, see the textbook [2].

Intuitively, a hyperfunction is defined by the difference between the boundary values of a complex analytic functions $F_{+} (z)$ in a subdomain of the upper half plane $Im z > 0$ and a complex analytic function $F_{-} (z)$ in a subdomain of the lower half plane $Im z < 0$

f (x) = F_{+} (x + i 0) - F_{-} (x - i 0),

where the analytic functions $F_{\pm} (z)$ are called defining functions of the hyperfunction $f (x)$ . We also use the notation

f (x) = [F (z)] = F (x + i 0) - F (x - i 0),

where

F (z) = {\begin{matrix} F_{+} (z) & Im z > 0 F_{-} (z) & Im z < 0. \end{matrix}

For example, the Dirac delta function $δ (x)$ is defined as a hyperfunction by

δ (x) = - \frac{1}{2 π i} (\frac{1}{x + i 0} - \frac{1}{x - i 0}),

and the Heaviside step function $Y (x)$ is defined as a hyperfunction by

Y (x) = - \frac{1}{2 π i} {log (- (x + i 0)) - log (- (x - i 0))}

where the complex logarithmic function $log z$ is the principal value, that is, the branch such that it takes a real value on the positive real axis. Figure 1 shows the graphs of the defining functions $F (z)$ of the Dirac delta function $δ (x)$ and the Heaviside step function $Y (x)$ . From these figures, we find easily that the difference between the boundary values on the real axis of the defining functions $F (z)$ gives the hyperfunctions.

Figure 1: The graphs of the defining functions $F (z)$ of the Dirac delta function $δ (x)$ and the Heaviside step function $Y (x)$ .

More precisely, hyperfunctions are defined as follows. Let $I$ be an open interval on the real axis $R$ and $D$ be a complex neighborhood of $I$ , that is, a complex domain including $I$ as a closed subset. We denote the set of complex analytic functions in $D ∖ I$ by $O (D ∖ I)$ , which can be regarded as a complex linear space by the usual summation and multiplication by complex numbers, and the set of complex analytic functions in $D$ by $O (D)$ , which can be regarded as a linear subspace of $O (D ∖ I)$ . Then, we can consider the quotient space $B (I) \equiv O (D ∖ I) / O (D)$ . We call an element of $B (I)$ , that is, an equivalence class $[F (z)]$ $(F (z) \in O (D ∖ I))$ a hyperfunction on $I$ and denote it by

f (x) = [F (z)] = F (x + i 0) - F (x - i 0),

where we call $F (z)$ a defining function of the hyperfunction $f (x)$ . We also use the notation

f (x) = [F_{+} (z), F_{-} (z)] = F_{+} (x + i 0) - F_{-} (x - i 0),

where

F_{\pm} (z) = F |_{D_{\pm}} (z) with D_{\pm} = {z \in C | \pm Im z > 0} .

In other words, a hyperfunction $f (x)$ on $I$ is a complex analytic function $F (z)$ in $D ∖ I$ , where another complex analytic function $˜ F (z)$ in $D ∖ I$ is identified with $F (z)$ if $˜ F (z) - F (z)$ is analytic in $D$ .

We define the value of a hyperfunction $f (x) = F_{+} (x + i 0) - F_{-} (x - i 0)$ at a point $x_{0} \in I$ by

f (x_{0}) = lim ϵ \to 0 + {F_{+} (x_{0} + i ϵ) - F_{-} (x_{0} - i ϵ)}

if the limit on the right hand side exists. We do not define the value $f (x_{0})$ if the limit does not exist.

We return to the discussion on the method of computing the oscillatory integral (1). From (6) and the fact that $F [f] (ζ)$ is a complex analytic function in the upper half plane $Im ζ > 0$ , we can regard the desired oscillatory integral $I$ as the value at $ξ = 0$ of the hyperfunction $f (x) = [F_{+} (ζ), F_{-} (ζ)]$ whose defining functions are $F_{+} (ζ) = F [f] (ζ)$ and $F_{-} (ζ) = 0$ .

4 Numerical examples

In this section, we discuss show some numerical examples that illustrate the effectiveness of the proposed method. All the computations were performed using programs coded in C++ with 100 decimal digit precision, where we use the multiple precision arithmetic library exflib [1].

We computed the oscillatory integrals

	$(1)$	$\int_{0}^{\infty} \frac{cos (x / 2) - cos x}{x} d x = log 2 = 0.6 9315 \dots,$
	$(2)$	$\int_{0}^{\infty} log x cos x d x = - \frac{π}{2} = - 1.57080 \dots,$
	$(3)$	$\int_{0}^{\infty} J_{0} (x) d x = 1,$
	$(4)$	$\int_{0}^{\infty} \frac{x J_{0} (x)}{x^{2} + 1} d x = K_{0} (1) = 0. 42102 \dots,$
	$(5)$	$\int_{0}^{\infty} \frac{J_{0} (x)}{(x^{2} + 1)^{1 / 2}} d x = K_{0} (\frac{1}{2}) I_{0} (\frac{1}{2}) = 0.98310 \dots,$
	$(6)$	$\int_{0}^{\infty} log x J_{0} (x) d x = - γ - log 2 = - 1.2704 \dots,$
	$(7)$	$\int_{0}^{\infty} \frac{x J_{1} (\sqrt{x^{2} + 1})}{\sqrt{x^{2} + 1}} d x = J_{0} (1) = 0.76520 \dots,$
	$(8)$	$\int_{0}^{\infty} \frac{Y_{0} (x)}{x^{2} + 1} d x = - K_{0} (1) = - 0 .42102 \dots .$

by the proposed method. We took the center $ζ_{0}$ of the Taylor series $F [f] (ζ) = \sum_{n = 0}^{\infty} c_{n} (ζ - ζ_{0})^{n}$ as $ζ_{0} = i$ , and we computed the Taylor coefficients $c_{0}, c_{1}, c_{2}, \dots, c_{100}$ given by the integral (4), which were computed using the DE formula. Table 1 show the relative error and the number of functional evaluations for $f (x)$ in (1) of the proposed method. Table 1 show that we can compute oscillatory integrals with high accuracy by the proposed method.

integral	(1)	(2)	(3)	(4)
relative error	5.4E-26	6.2E-35	3.8E-34	1.4E-36
number of functional evaluations	917	964	957	927
integral	(5)	(6)	(7)	(8)
relative error	1.3E-35	3.8E-36	1.1E-33	2.1E-37
number of functional evaluations	954	958	927	947

Table 1: The relative errors and the number of functional evaluations for

f (x)

in (1) of the proposed method applied to the integrals (1-8).

Here, we point out an interesting fact. Since we took $ζ_{0} = i$ , the integral in (4) which gives the Taylor coefficients $c_{n}$ is given by

which does not involve any oscillatory function. It means that we can obtain oscillatory integrals without computing oscillatory integrals. In other words, we can obtain Fourier transform type integrals by computing Laplace transform type integrals.

As a comparison with the results by a previous method, we computed some of the integrals by the Euler transform for computing alternating series. Namely, we express the desired oscillatory integral as

\int_{0}^{\infty} f (x) d x = \infty \sum k = 0 (- 1)^{k} ∣ ∣ ∣ \int_{a_{k}}^{a_{k + 1}} f (x) d x ∣ ∣ ∣,

(7)

where $0 = a_{0} < a_{1} < a_{2} < \dots$ and

\int_{a_{k}}^{a_{k + 1}} f (x) d x \cdot \int_{a_{k + 1}}^{a_{k + 2}} f (x) d x < 0 (k = 0, 1, 2, \dots) .

We computed each integral on the right hand side of (7) by the 100-point Gauss-Legendre formula, and we obtained the desired integral by applying the Euler transform to the alternating series (7). Table 2 show the relative errors of the method using the Euler transform applied to the integrals (2,3,4,5). The number of functional evaluations for $f (x)$ in (1) are 5000 in all the computations of the integrals. Comparing the results in Table 1 and 2, the proposed method is far superior to the method using the Euler transform.

integral	(2)	(3)	(4)	(5)
relative error	6.3E-25	1.9E-25	3.1E-25	1.4E-25

Table 2: The relative errors of the method using the Euler transform applied to the integrals (2,3,4,5). The number of functional evaluations for

f (x)

in (1) are 5000 in all the computations of the integrals.

5 Summary

In this paper, we proposed a numerical method of computing an integral involving a slowly decaying oscillatory function. In the proposed method, we introduce a complex analytic function in the upper half plane and obtain the desired integral by the analytic continuation of this analytic function onto the real axis. The analytic function is obtained as a Taylor series by computing its coefficients, which are given by integrals involving exponentially decaying functions and easily evaluated by the conventional numerical integration formula. The analytic continuation is performed by converting the Taylor series into a continued fraction using the quotient difference (QD) algorithm. Then, the desired integral is obtained as the value of the analytic function at the origin. We also remarked that the proposed method is closely related to hyperfunction theory in the sense that the desired oscillatory integral can be regarded as the boundary value of the hyperfunction whose defining function is the above complex analytic function.

As presented in this paper, hyperfunctions, that is, the boundary values of analytic functions appear in many studies of science and engineering. In the evaluation of hyperfunctions, the techniques of analytic continuation are necessary. Therefore, the author believes that the study and the computation of hyperfunctions together with analytic continuation will be important and useful in applied mathematics.

Problems for future study related to this paper are, for example, theoretical error estimate of the proposed method and the reduction of the computational cost of analytic continuation, which needs multiple precision arithmetic in the QD algorithm applied to the conversion of a Taylor series into a continued fraction.

References

[1] H. Fujiwara Exflib information. Note: http://www-an.acs.i.kyoto-u.ac.jp/~fujiwara/exflib/ Cited by: §4.
[2] U. Graf (2010) Introduction to hyperfunctions and their integral transforms — an applied and computational approach. Birkha̋user, Basel. Cited by: §3.
[3] P. Henrici (1977) Applied and computational complex analysis. Vol. 2, John Wiley & Sons, New York. Cited by: item 2.
[4] H. Ogata (2018) A numerical method of fourier transform based on hyperfunction theory. Note: arXiv:1808.03460v1 [math.NA] Cited by: §1.
[5] T. Ooura and M. Mori (1991) The double exponential formula for oscilattory functions over the half infinite interval. J. Comput. Appl. Math. 38, pp. 353–360. Cited by: §1.
[6] T. Ooura and M. Mori (1999) A robust double exponential formula for fourier type integrals. J. Comput. Appl. Math. 112, pp. 229–241. Cited by: §1.
[7] M. Sato (1959) THeory of hyperfunctions. J. Fac. Sci. Univ. Tokyo, Sect. 1A Math 8, pp. 139–193. Cited by: §1.
[8] H. Takahasi and M. Mori (1978) Double exponential formulas for numerical integration. Publ. Res. Inst. Math. Sci., Kyoto Univ. 339, pp. 721–741. Cited by: §1.
[9] H. Toda and H. Ono (1978) Some remarks for efficient usage of the double exponential formulas (in japanese). Kokyuroku, Res. Inst. Math. Sci. Kyoto Univ. 339, pp. 74–109. Cited by: §1.

A numerical method of computing oscillatory integral related to hyperfunction theory

Numerical method for computing Hadamard finite-part integrals with a non-integral power singularity at an endpoint

A numerical method for Hadamard finite-part integrals with an integral power singularity at an endpoint

Efficient extraction of resonant states in systems with defects

Reconstructing Stieltjes functions from their approximate values: a search for a needle in a hay stack

Unsteady aerodynamics of porous aerofoils

Revisiting integral functionals of geometric Brownian motion

Computation of the Complex Error Function using Modified Trapezoidal Rules

1 Introduction

2 Numerical method

3 Relation with hyperfunction theory

4 Numerical examples

5 Summary

References

	$e_{0}^{(n)} = 0, q_{1}^{(n)} = \frac{c_{n + 1}}{c_{n}} (n = 0, 1, 2, \dots),$
	$\begin{matrix} e_{k}^{(n)} = q_{k}^{(n + 1)} - q_{k}^{(n)} + e_{k - 1}^{(n + 1)}, q_{k + 1}^{(n)} = \frac{e_{k}^{(n + 1)}}{e_{k}^{(n)}} q_{k}^{(n + 1)} (n = 0, 1, 2, \dots; k = 1, 2, \dots), \end{matrix}$

	$P_{- 1} = 0, Q_{- 1} = 1, P_{0} = c_{0}, Q_{0} = 1,$
	${\begin{matrix} P_{k} (ζ) = a_{k} (ζ - ζ_{0}) P_{k - 2} (ζ) + P_{k - 1} (ζ) Q_{k} (ζ) = a_{k} (ζ - ζ_{0}) Q_{k - 2} (ζ) + Q_{k - 1} (ζ) \end{matrix} (k = 1, 2, \dots) .$

x $Re z$ y $Im z$ z $Re F (z)$	x $Re z$ y $Im z$ z $Re F (z)$
$δ (x)$	$Y (x)$

A numerical method of computing oscillatory integral related to hyperfunction theory

Related Research

Numerical method for computing Hadamard finite-part integrals with a non-integral power singularity at an endpoint

A numerical method for Hadamard finite-part integrals with an integral power singularity at an endpoint

Efficient extraction of resonant states in systems with defects

Reconstructing Stieltjes functions from their approximate values: a search for a needle in a hay stack

Unsteady aerodynamics of porous aerofoils

Revisiting integral functionals of geometric Brownian motion

Computation of the Complex Error Function using Modified Trapezoidal Rules

1 Introduction

2 Numerical method

3 Relation with hyperfunction theory

4 Numerical examples

5 Summary

References