A numerical method of computing oscillatory integral related to hyperfunction theory

In this paper, we propose a numerical method of computing an integral whose integrand is a slowly decaying oscillatory function. In the proposed method, we consider a complex analytic function in the upper-half complex plane, which is defined by an integral of the Fourier-Laplace transform type, and we obtain the desired integral by the analytic continuation of this analytic function onto the real axis using a continued fraction. We also remark that the proposed method is related to hyperfunction theory, a theory of generalized functions based on complex function theory. Numerical examples show the effectiveness of the proposed method.

READ FULL TEXT VIEW PDF

page 1

page 2

page 3

page 4

09/25/2019

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

In this paper, we propose a numerical method of computing a Hadamard fin...
09/19/2019

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

In this paper, we propose a numerical method for computing Hadamard fini...
03/23/2022

Efficient extraction of resonant states in systems with defects

We introduce a new numerical method to compute resonances induced by loc...
01/07/2021

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

Material response of passive media to external influences is described b...
11/18/2019

Unsteady aerodynamics of porous aerofoils

We extend unsteady thin aerofoil theory to model aerofoils with generali...
01/31/2020

Revisiting integral functionals of geometric Brownian motion

In this paper we revisit the integral functional of geometric Brownian m...
10/12/2020

Computation of the Complex Error Function using Modified Trapezoidal Rules

In this paper we propose a method for computing the Faddeeva function w(...

1 Introduction

In this paper, we consider the computation of an integral

(1)

where is a slowly decaying oscillatory function, such as

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

(2)

which is a complex analytic function in the upper-half complex plane under some conditions on . We obtain the function in the upper-half plane, and we compute the analytic continuation of on the real axis . Then, we obtain the desired integral by

(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 . From (3), the desired integral is a boundary value on of the function analytic in the upper half plane , and, in this sense, the integral 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

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

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 is a slowly decaying oscillatory function and propose a numerical method for computing it. As mentioned in the previous section,

  1. we obtain a function , which is a complex analytic function in the upper half plane under some conditions on ,

  2. and we obtain the desired integral by the analytic continuation of onto the real axis.

The details of each step are as follows.

  1. (Computation of in the upper half plane ) We obtain the analytic function as a Taylor series

    where is a complex constant given by the user such that . The Taylor coefficient is given by the integral

    (4)

    This integral involves an exponentially decaying factor , where we remark that , and it is easy to compute it by conventional numerical integration formulas.

  2. (Analytic continuation of ) We obtain the analytic continuation of onto the real axis 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

    into a continued fraction

    (5)

    using the so-called quotient difference (QD) algorithm (see Theorem 12.4c of [3]). Namely, we generate the sequences and by the recurrence relations

    and, then, we obtain

    Once we obtain the coefficients of the continued fraction (5), we can compute the continued fraction by

    where and are the polynomials given by the recurrence relations

    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

(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 in a subdomain of the upper half plane and a complex analytic function in a subdomain of the lower half plane

where the analytic functions are called defining functions of the hyperfunction . We also use the notation

where

For example, the Dirac delta function is defined as a hyperfunction by

and the Heaviside step function is defined as a hyperfunction by

where the complex logarithmic function 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 of the Dirac delta function and the Heaviside step function . From these figures, we find easily that the difference between the boundary values on the real axis of the defining functions gives the hyperfunctions.

x y z

x y z

Figure 1: The graphs of the defining functions of the Dirac delta function and the Heaviside step function .

More precisely, hyperfunctions are defined as follows. Let be an open interval on the real axis and be a complex neighborhood of , that is, a complex domain including as a closed subset. We denote the set of complex analytic functions in by , 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 by , which can be regarded as a linear subspace of . Then, we can consider the quotient space . We call an element of , that is, an equivalence class a hyperfunction on and denote it by

where we call a defining function of the hyperfunction . We also use the notation

where

In other words, a hyperfunction on is a complex analytic function in , where another complex analytic function in is identified with if is analytic in .

We define the value of a hyperfunction at a point by

if the limit on the right hand side exists. We do not define the value 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 is a complex analytic function in the upper half plane , we can regard the desired oscillatory integral as the value at of the hyperfunction whose defining functions are and .

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

by the proposed method. We took the center of the Taylor series as , and we computed the Taylor coefficients 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 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 in (1) of the proposed method applied to the integrals (1-8).

Here, we point out an interesting fact. Since we took , the integral in (4) which gives the Taylor coefficients 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

(7)

where and

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 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 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.