Double-Exponential transformation:
A quick review of a Japanese tradition^††thanks: This is a preliminary version of an article to be included in ICIAM 2023, Tokyo Intelligencer.

The DE formula was motivated by the fact that the trapezoidal rule is highly effective for integrals over the infinite interval $(-\infty,+\infty)$. For an integral

for example, we employ a change of variable $x=\phi(t)$ using some function $\phi(t)$ satisfying $\phi(-\infty)=-1$ and $\phi(+\infty)=1$, and apply the trapezoidal rule to the transformed integral

to obtain an infinite sum of discretization

A finite-term approximation to this infinite sum results in an integration formula

Such combination of the trapezoidal rule with a change of variables was conceived by several authors [2, 24, 25, 35] around 1970.

The error $I-I_{h}^{(N)}$ of the formula (3) consists of two parts, the error $E_{\rm D}\equiv I-I_{h}$ incurred by discretization (2) and the error $E_{\rm T}\equiv I_{h}-I_{h}^{(N)}$ caused by truncation of an infinite sum $I_{h}$ to a finite sum $I_{h}^{(N)}$.

The major findings of Takahasi and Mori consisted of two ingredients. The first was that the double exponential decay of the transformed integrand $f(\phi(t))\phi^{\prime}(t)$ achieves the optimal balance (or trade-off) between the discretization error $E_{\rm D}$ and the truncation error $E_{\rm T}$. The second finding was that a concrete choice of

is suitable for this purpose thanks to the double exponential decay shown in (1). With this particular function $\phi(t)$ the formula (3) reads

which is sometimes called “tanh-sinh quadrature.” The error of this formula is estimated roughly as

with some $C>0$. The DE formula has an additional feature that it is robust against end-point singularities of integrands.

The idea of the DE formula can be applied to integrals over other types of intervals of integration. For example,

Such formulas are also referred to as the double exponential formula. The DE formula is available in Mathematica (NIntegrate), Python library SymPy, Python library mpmath, C++ library Boost, Haskell package integration, etc.

Optimality of the DE transformation (4) was discussed already by Takahasi and Mori [37]. Numerical examples also support its optimality. Figure 1 (taken from [5]) shows the comparison of the DE transformation (4) against other transformations

for

that has integrable singularities at both ends of the interval of integration. The DE formula converges much faster than others. It is known that the tanh-rule (using $\phi(t)=\tanh t$) has the (rough) convergence rate $\exp(-C\sqrt{N})$, in contrast to $\exp(-CN/\log N)$ in (5) of the DE formula.

The optimality argument of [37], based on complex function theory, was convincing enough for the majority of scientists and engineers, but not perfectly satisfactory for theoreticians. Rigorous mathematical argument for optimality of the DE formula was addressed by Masaaki Sugihara [28, 29, 30] in the 1980–1990s in a manner comparable to Stenger’s framework [26] for optimality of the tanh rule. It is shown in [30] (also [42]) that the DE formula is optimal with respect to a certain class (Hardy space) of integrand functions.

In principle, for each class of integrand functions we may be able to find an optimal quadrature formula, and the optimal formula naturally depends on our choice of the admissible class of integrands. Thus the optimality of a quadrature formula is only relative. However, it was shown by Sugihara that no nontrivial class of integrand functions exists that admits a quadrature formula with smaller errors than the DE formula. We can interpret this fact as the absolute optimality of the DE formula.

For Fourier-type integrals like

the DE formula like (6) is not very successful. To cope with Fourier-type integrals, a novel technique, in the spirit of DE transformation, was proposed by Ooura and Mori [22, 23]. In [22] they proposed to use

($K>0$), which maps $(-\infty,+\infty)$ to $(0,+\infty)$ in such a way that (i) $\phi^{\prime}(t)\to 0$ double exponentially as $t\to-\infty$ and (ii) $\phi(t)\to t$ double exponentially as $t\to+\infty$. The proposed formula changes the variable by $x=M\phi(t)$ to obtain

to which the trapezoidal rule with equal mesh $h$ is applied, where $M$ and $h$ are chosen to satisfy $Mh=\pi$. The transformed integrand decays double-exponentially toward $t\to-\infty$ because of the factor $\phi^{\prime}(t)$ and also toward $t\to+\infty$ because $M\phi(t)$ for $t=kh$ (sample point of the trapezoidal rule) tends double-exponentially to $Mt=Mkh=k\pi$, at which sine function vanishes. Another (improved) transformation function

is given in [23], where $\beta=1/4$ and $\alpha=\beta/\sqrt{1+M\log(1+M)/(4\pi)}$.

In 1969, prior to the DE formula, a remarkable quadrature formula was proposed by Masao Iri, Sigeiti Moriguti, and Yoshimitsu Takasawa [2]. The formula is known today as the “IMT rule,” which name was introduced in [35] and used in [1].

For an integral

over $[0,1]$, the IMT rule applies the trapezoidal rule to the integral

resulting from the transformation by

where

is a normalizing constant to render $\phi(1)=1$.

The transformed integrand $g(t)=f(\phi(t))\phi^{\prime}(t)$ has the property that all the derivatives $g^{(j)}(t)$ $(j=1,2,\ldots)$ vanish at $t=0,1$. By the Euler–Maclaurin formula, this indicates that the IMT rule should be highly accurate. Indeed, it was shown in [2] via a complex analytic method that the error of the IMT rule can be estimated roughly as $\exp(-C\sqrt{N})$, which is much better than $N^{-4}$ of the Simpson rule, say, but not as good as $\exp(-CN/\log N)$ of the DE formula. Variants of the IMT rule have been proposed for possible improvement [4, 10, 21, 29], but it turned out that an IMT-type rule, transforming $\int_{0}^{1}{\rm d}x$ to $\int_{0}^{1}{\rm d}t$ rather than to $\int_{-\infty}^{+\infty}{\rm d}t$, cannot outperform the DE formula.

The Sinc approximation of a function $f(x)$ over $(-\infty,\infty)$ is given by

where $S(k,h)(x)$ is the so-called Sinc function defined by

and the step size $h$ is chosen appropriately, depending on $N$. The technique of variable transformation $x=\phi(t)$ is also effective in this context. By applying the formula (8) to $f(\phi(t))$ we obtain

or equivalently,

To approximate $f(x)$ over $[0,1]$, for example, we choose

etc. The methods using (9) and (10) are often called the SE- and DE-Sinc approximations, respectively. The error of the SE-Sinc approximation is roughly $\exp(-C\sqrt{N})$ and that of the DE-Sinc approximation is $\exp(-CN/\log N)$.

These approximation schemes are compared in Fig. 2 (taken from [34]) for function

over $[0,1]$. In Fig. 2, “Ordinary-Sinc” means the SE-Sinc approximation using (9), and the polynomial interpolation with the Chebyshev nodes is included for comparison.

Detailed theoretical analyses on the DE-Sinc method can be found in Sugihara [33] as well as Tanaka et al. [41] and Okayama et al. [16, 20]. An optimization technique is used to improve the DE-Sinc method in Tanaka and Sugihara [39].

Once a function approximation scheme is at hand, we can apply it to a variety of numerical problems. Indeed this is also the case with the DE-Sinc approximation as follows.

• •

  Indefinite integration by Muhammad and Mori [8], Tanaka et al. [40], and Okayama and Tanaka [19].

• •

  Initial value problem of differential equations by Nurmuhammad et al. [11] and Okayama [15].

• •

  Boundary value problem of differential equations by Sugihara [32], followed by Nurmuhammad et al. [12, 13] and Mori et al. [6].

• •

  Volterra integral equation by Muhammad et al. [9] and Okayama et al. [18].

• •

  Fredholm integral equation by Kobayashi et al. [3], Muhammad et al. [9], and Okayama et al. [17].