Reference

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

Numerically integrate f(x) over an arbitrary interval [a, b] and return the integral value I and an estimated error E. 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. Optionally, one can divide the integral interval [a, b, c...], which returns ∫ f(x)dx in [a, b] + ∫f(x)dx in [b, c[1]] + ⋯. It is worth noting that each endpoint allows discontinuity or singularity.

The integrand f can also return any value other than a scalar, as far as +, -, multiplication by real values, and LinearAlgebra.norm, are implemented. For example, Vector or Array of numbers are acceptable although, unfortunately, it may not be very performant.

Examples

julia> using DoubleExponentialFormulas

julia> using LinearAlgebra: norm

julia> f(x) = 2/(1 + x^2);

julia> I, E = quadde(f, -1, 1);

julia> I ≈ π
true

julia> E ≤ sqrt(eps(Float64))*norm(I)
true

julia> I, E = quadde(x -> [1/(1 + x^2), 2/(1 + x^2)], 0, Inf);

julia> I ≈ [π/2, π]
true

julia> E ≤ sqrt(eps(Float64))*norm(I)
true

julia> I, E = quadde(x -> 1/sqrt(abs(x)), -1, 0, 1);  # singular point at x = 0

julia> I ≈ 4
true

julia> E ≤ sqrt(eps(Float64))*norm(I)
true

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

Numerically integrate f(x) over an arbitrary interval [a, b] and return the integral value I and an estimated error E. The quaddeo function is specialized for the decaying oscillatory integrands,

f(x) = f₁(x)sin(ωx + θ),

where f₁(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ω). 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. Optionally, one can divide the integral interval [a, b, c...], which returns ∫ f(x)dx in [a, b] + ∫f(x)dx in [b, c[1]] + ⋯. It is worth noting that each endpoint allows discontinuity or singularity.

The integrand f can also return any value other than a scalar, as far as +, -, multiplication by real values, and LinearAlgebra.norm, are implemented. For example, Vector or Array of numbers are acceptable although, unfortunately, it may not be very performant.

Examples

julia> using DoubleExponentialFormulas

julia> using LinearAlgebra: norm

julia> f(x) = sin(x)/x;

julia> I, E = quaddeo(f, 1.0, 0.0, 0.0, Inf);

julia> I ≈ π/2
true

julia> E ≤ sqrt(eps(Float64))*norm(I)
true

QuadDE(T::Type{<:AbstractFloat}; maxlevel::Integer=10, h0::Real=one(T)/8)

A callable object to integrate an function over an arbitrary interval, which uses the double exponential formulas also known as the tanh-sinh quadrature and its variants. It utilizes the change of variables to transform the integrand into a form well-suited to the trapezoidal rule.

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, see below.

The type T represents the accuracy of interval. The integrand should accept values x<:T as its parameter.

I, E = (q::QuadDE{T})(f::Function, a::Real, b::Real, c::Real...;
                      atol::Real=zero(T),
                      rtol::Real=atol>0 ? zero(T) : sqrt(eps(T)))
                      where {T<:AbstractFloat}

Numerically integrate f(x) over an arbitrary interval [a, b] and return the integral value I and an estimated error E. 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. Optionally, one can divide the integral interval [a, b, c...], which returns ∫f(x)dx in [a, b] + ∫f(x)dx in [b, c[1]] + ⋯. It is worth noting that each endpoint allows discontinuity or singularity.

The integrand f can also return any value other than a scalar, as far as +, -, multiplication by real values, and LinearAlgebra.norm, are implemented. For example, Vector or Array of numbers are acceptable although, unfortunately, it may not be very performant.

Examples

julia> using DoubleExponentialFormulas

julia> using LinearAlgebra: norm

julia> qde = QuadDE(Float64);

julia> f(x) = 2/(1 + x^2);

julia> I, E = qde(f, -1, 1);

julia> I ≈ π
true

julia> E ≤ sqrt(eps(Float64))*norm(I)
true

julia> g(x) = [1/(1 + x^2), 2/(1 + x^2)];

julia> I, E = qde(g, 0, Inf);

julia> I ≈ [π/2, π]
true

julia> E ≤ sqrt(eps(Float64))*norm(I)
true

julia> h(x) = 1/sqrt(abs(x));  # singular point at x = 0

julia> I, E = qde(h, -1, 0, 1);

julia> I ≈ 4
true

julia> E ≤ sqrt(eps(Float64))*norm(I)
true

QuadES(T::Type{<:AbstractFloat}; maxlevel::Integer=10, h0::Real=one(T)/8)

A callable object to integrate a function over the range [0, ∞) using the exp-sinh quadrature. It utilizes the change of variables to transform the integrand into a form well-suited to the trapezoidal rule.

QuadES 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, see below.

The type T represents the accuracy of interval. The integrand should accept values x<:T as its parameter.

I, E = (q::QuadES)(f::Function;
                   atol::Real=zero(T),
                   rtol::Real=atol>0 ? zero(T) : sqrt(eps(T)))
                   where {T<:AbstractFloat}

Numerically integrate f(x) over the interval [0, ∞) and return the integral value I and an estimated error E. 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.

The integrand f can also return any value other than a scalar, as far as +, -, multiplication by real values, and LinearAlgebra.norm, are implemented. For example, Vector or Array of numbers are acceptable although, unfortunately, it may not be very performant.

Examples

julia> using DoubleExponentialFormulas

julia> using LinearAlgebra: norm

julia> qes = QuadES(Float64);

julia> f(x) = 2/(1 + x^2);

julia> I, E = qes(f);

julia> I ≈ π
true

julia> E ≤ sqrt(eps(Float64))*norm(I)
true

julia> I, E = qes(x -> [1/(1 + x^2), 2/(1 + x^2)]);

julia> I ≈ [π/2, π]
true

julia> E ≤ sqrt(eps(Float64))*norm(I)
true

QuadSS(T::Type{<:AbstractFloat}; maxlevel::Integer=10, h0::Real=one(T)/8)

A callable object to integrate a function over the range (-∞, ∞) using the sinh-sinh quadrature. It utilizes the change of variables to transform the integrand into a form well-suited to the trapezoidal rule.

QuadSS 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, see below.

The type T represents the accuracy of interval. The integrand should accept values x<:T as its parameter.

I, E = (q::QuadSS)(f::Function;
                   atol::Real=zero(T),
                   rtol::Real=atol>0 ? zero(T) : sqrt(eps(T)))
                   where {T<:AbstractFloat}

Numerically integrate f(x) over the interval (-∞, ∞) and return the integral value I and an estimated error E. 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.

The integrand f can also return any value other than a scalar, as far as +, -, multiplication by real values, and LinearAlgebra.norm, are implemented. For example, Vector or Array of numbers are acceptable although, unfortunately, it may not be very performant.

Examples

julia> using DoubleExponentialFormulas

julia> using LinearAlgebra: norm

julia> qss = QuadSS(Float64);

julia> f(x) = 2/(1 + x^2);

julia> I, E = qss(f);

julia> I ≈ 2π
true

julia> E ≤ sqrt(eps(Float64))*norm(I)
true

julia> I, E = qss(x -> [1/(1 + x^2), 2/(1 + x^2)]);

julia> I ≈ [π, 2π]
true

julia> E ≤ sqrt(eps(Float64))*norm(I)
true

QuadTS(T::Type{<:AbstractFloat}; maxlevel::Integer=10, h0::Real=one(T)/8)

A callable object to integrate a function over the range [-1, 1] using the tanh-sinh quadrature. It utilizes the change of variables to transform the integrand into a form well-suited to the trapezoidal rule.

QuadTS 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, see below.

The type T represents the accuracy of interval. The integrand should accept values x<:T as its parameter.

I, E = (q::QuadTS)(f::Function;
                   atol::Real=zero(T),
                   rtol::Real=atol>0 ? zero(T) : sqrt(eps(T)))
                   where {T<:AbstractFloat}

Numerically integrate f(x) over the interval [-1, 1] and return the integral value I and an estimated error E. 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.

The integrand f can also return any value other than a scalar, as far as +, -, multiplication by real values, and LinearAlgebra.norm, are implemented. For example, Vector or Array of numbers are acceptable although, unfortunately, it may not be very performant.

Examples

julia> using DoubleExponentialFormulas

julia> using LinearAlgebra: norm

julia> qts = QuadTS(Float64);

julia> f(x) = 2/(1 + x^2);

julia> I, E = qts(f);

julia> I ≈ π
true

julia> E ≤ sqrt(eps(Float64))*norm(I)
true

julia> I, E = qts(x -> [1/(1 + x^2), 2/(1 + x^2)]);

julia> I ≈ [π/2, π]
true

julia> E ≤ sqrt(eps(Float64))*norm(I)
true