## DoubleExponentialFormulas.jl

One-dimensional numerical integration using the double exponential formula
Author machakann
Popularity
6 Stars
Updated Last
1 Year Ago
Started In
December 2019

# DoubleExponentialFormulas

## Introduction

This package provides functions for one-dimensional numerical integration using the double exponential formula [1,2] also known as the tanh-sinh quadrature and its variants.

## Instllation

Press `]` on a Julia REPL to enter the Pkg mode and run the following command.

``````add https://github.com/machakann/DoubleExponentialFormulas.jl.git
``````

## Usage

### Handy interface with Float64 precision

``````    I, E = quadde(f::Function, a::Real, b::Real, c::Real...;
atol::Real=zero(Float64),
rtol::Real=atol>0 ? zero(Float64) : sqrt(eps(Float64)))
``````

The `quadde` function provides a handy way to integrate a function `f(x)` over an arbitrary interval.

```using DoubleExponentialFormulas
using LinearAlgebra: norm

f(x) = 1/(1 + x^2)
I, E = quadde(f, -1, 1)

I ≈ π/2    # true
E ≤ sqrt(eps(Float64))*norm(I)  # true```

The above example computes `∫ 1/(1+x^2) dx in [-1, 1]`. The `I` is the obtained integral value and the `E` is an estimated numerical error. The `E` is not exactly equal to the difference from the true value. However, one can expect that the integral value `I` is converged if `E <= max(atol, rtol*norm(I))` is true. Otherwise, the obtained `I` would be unreliable; the number of repetitions exceeds the `maxlevel` before converged.

Half-infinite intervals and the infinite interval are also valid, as far as the integral is convergent.

```# Computes ∫ 1/(1+x^2) dx in [0, ∞)
I, E = quadde(x -> 1/(1 + x^2), 0, Inf)
I ≈ π/2    # true

# Computes ∫ 1/(1+x^2) dx in (-∞, ∞)
I, E = quadde(x -> 1/(1 + x^2), -Inf, Inf)
I ≈ π      # true```

Optionally, one can divide the integral interval [a, b, c, ...], which returns `∫f(x)dx in [a, b] + ∫f(x)dx in [b, c] + ⋯`. It is worth noting that discontinuity or singularity is allowed at the endpoints.

```# Computes ∫ 1/sqrt(|x|) dx in (-∞, ∞)
# The integrand has a singular point at x = 0
I, E = quadde(x -> 1/sqrt(abs(x)), -1, 0, 1)
I ≈ 4    # true```

### Optimized numerical integrators

User can get an optimized integrators, for example, for better accuracy; `QuadDE` will provides the functionality.

``````qde = QuadDE(BigFloat; h0=one(BigFloat)/8, maxlevel=10)
qde(x -> 2/(1 + x^2), -1,  1)
``````

User can specify the required precision as a type (`T<:AbstractFloat`), the starting step size `h0` and the maximum number of repetition `maxlevel`. The `h0` and `maxlevel` shown above are the default values, so it can be omitted. `QuadDE` instance is an callable object which has the same interface of `quadde`, actually `quadde` is an alias to `QuadDE(Float64)(...)` with a precalculated instance.

`QuadDE` tries to calculate integral values `maxlevel` times at a maximum; the step size of a trapezoid is started from `h0` and is halved in each following repetition for finer accuracy. The repetition is terminated when the difference from the previous estimation gets smaller than a certain threshold. The threshold is determined by the runtime parameters, `atol` or `rtol`.

Using smaller `h0` may help if the integrand `f(x)` includes fine structure, such as spikes, in the integral interval. However, it seems that the subdivision of the interval would be more effective in many cases. Try subdivision first, and then think of an optimized integrator.

See documentation for more details.

### Numerical integrator for decaying oscillatory integrands

``````    quaddeo(f::Function, ω::Real, θ::Real, a::Real, b::Real;
h0::Real=one(ω)/5, maxlevel::Integer=12,
atol::Real=zero(ω),
rtol::Real=atol>0 ? zero(atol) : sqrt(eps(typeof(atol))))
``````

The `quaddeo` function is specialized for the decaying oscillatory integrands, [3-5]

``````f(x) = g(x)sin(ωx + θ),
``````

where `g(x)` is a decaying algebraic function. `ω` and `θ` are the frequency and the phase of the oscillatory part of the integrand. If the oscillatory part is `sin(ωx)`, then `θ = 0.0`; if it is `cos(ωx)` instead, then `θ = π/(2ω)`.

```using DoubleExponentialFormulas
using LinearAlgebra: norm

f(x) = sin(x)/x;
I, E = quaddeo(f, 1.0, 0.0, 0.0, Inf);
I ≈ π/2  # true```

## References

1. Takahasi, H.; Mori, M. Double Exponential Formulas for Numerical Integration. Publ. Res. Inst. Math. Sci. 1973, 9 (3), 721–741. 10.2977/prims/1195192451.

2. Mori, M.; Sugihara, M. The Double-Exponential Transformation in Numerical Analysis. J. Comput. Appl. Math. 2001, 127 (1–2), 287–296. 10.1016/S0377-0427(00)00501-X.

3. Ooura, T.; Mori, M. The double exponential formula for oscillatory functions over the half infinite interval. J. Comput. Appl. Math. 1991, 38, 353–360. 10.1016/0377-0427(91)90181-I

4. http://www.kurims.kyoto-u.ac.jp/~ooura/intde-j.html

5. http://www.kurims.kyoto-u.ac.jp/~ooura/intdefaq-j.html