# Reference

## Index

`DoubleExponentialFormulas.quadde`

`DoubleExponentialFormulas.quaddeo`

`DoubleExponentialFormulas.QuadDE`

`DoubleExponentialFormulas.QuadES`

`DoubleExponentialFormulas.QuadSS`

`DoubleExponentialFormulas.QuadTS`

`DoubleExponentialFormulas.quadde`

— Method.```
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
```

`DoubleExponentialFormulas.quaddeo`

— Method.```
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
```

`DoubleExponentialFormulas.QuadDE`

— Type.`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 formula*s 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
```

`DoubleExponentialFormulas.QuadES`

— Type.`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.

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

`DoubleExponentialFormulas.QuadSS`

— Type.`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.

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

`DoubleExponentialFormulas.QuadTS`

— Type.`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.

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

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