DoubleExponentialFormulas
Introduction
This package provides functions for onedimensional numerical integration using the double exponential formula [1,2] also known as the tanhsinh 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.
Halfinfinite 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, [35]
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

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

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

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/03770427(91)90181I