Library

QuadratureRule{T,N}

A quadrature rule on the reference interval $[0,1]$, that is, an approximation of the form

\[\int_0^1 f(x) \, dx \approx \sum_{i=1}^{N} b_i \, f(c_i) ,\]

with nodes $c_i$ and weights $b_i$. All rules in this package are normalised to $[0,1]$, so that the weights sum to one. An integral over a general interval $[a,b]$ is obtained by the affine change of variables

\[\int_a^b f(x) \, dx = (b-a) \int_0^1 f \big( a + (b-a) \, \xi \big) \, d\xi .\]

Type parameters

• T: element type of the nodes and weights, e.g. Float64 or BigFloat.
• N: number of nodes, encoded in the type so that it is available at compile time.

Fields

• order::Int: the order $p$ of the rule; a rule of order $p$ integrates polynomials of degree $\le p-1$ exactly (see order).
• nodes::Vector{T}: the nodes $c_i \in [0,1]$.
• weights::Vector{T}: the weights $b_i$.

Examples

A rule is usually obtained from one of the constructors listed in the Quadrature Rules section rather than built by hand:

julia> quad = TrapezoidalQuadrature()
QuadratureRule{Float64, 2}(2, [0.0, 1.0], [0.5, 0.5])

julia> quad(x -> x^2)
0.5

See also nodes, weights, order, nnodes.

nnodes(quad::QuadratureRule)

Return the number of nodes of quad.

This is the type parameter N of QuadratureRule and therefore known at compile time.

julia> nnodes(GaussLegendreQuadrature(3))
3

nodes(quad::QuadratureRule)

Return the vector of nodes $c_i \in [0,1]$ of quad.

julia> nodes(TrapezoidalQuadrature())
2-element Vector{Float64}:
 0.0
 1.0

order(quad::QuadratureRule)

Return the order $p$ of quad.

A quadrature rule has order $p$ if it integrates all polynomials of degree $\le p-1$ exactly, that is, if the moment conditions

\[\sum_{i} b_i \, c_i^{k} = \frac{1}{k+1} , \qquad k = 0, 1, \dots, p-1 ,\]

are satisfied. The order is sharp for every rule in this package: the rule fails for some polynomial of degree $p$. It follows that two rules with the same nodes and weights also have the same order, so rules that coincide compare equal across families, for instance the three-node Clenshaw-Curtis and Lobatto-Legendre rules, both of which are Simpson's rule.

julia> order(GaussLegendreQuadrature(3))     # exact up to degree 5
6

julia> order(LobattoLegendreQuadrature(3))   # exact up to degree 3
4

julia> order(ClenshawCurtisQuadrature(3))    # the same rule, hence the same order
4

weights(quad::QuadratureRule)

Return the vector of weights $b_i$ of quad.

As all rules are normalised to the interval $[0,1]$, the weights sum to one.

julia> weights(TrapezoidalQuadrature())
2-element Vector{Float64}:
 0.5
 0.5

eachindex(quad::QuadratureRule)

Return an iterator over the indices of the nodes and weights of quad, so that all pairs $(b_i, c_i)$ can be visited in a loop.

julia> quad = TrapezoidalQuadrature();

julia> for i in eachindex(quad)
           println(nodes(quad)[i], " ", weights(quad)[i])
       end
0.0 0.5
1.0 0.5

eltype(quad::QuadratureRule)

Return the element type T of the nodes and weights of quad.

julia> eltype(GaussLegendreQuadrature(BigFloat, 3))
BigFloat

==(quad1::QuadratureRule, quad2::QuadratureRule)

Compare two quadrature rules by value, i.e., check that their orders, nodes and weights agree. Rules of different element type may compare equal; use isequal to distinguish them.

isequal(quad1::QuadratureRule, quad2::QuadratureRule)

Check that two quadrature rules are equal by value and of the same type, that is, that they agree in their element type T and their number of nodes N in addition to their orders, nodes and weights.

isapprox(quad1::QuadratureRule, quad2::QuadratureRule; kwargs...)

Check that two quadrature rules have the same order and approximately equal nodes and weights. All keyword arguments are forwarded to isapprox for the nodes and weights.

This is the appropriate comparison for rules computed by different algorithms or in different working precisions, e.g. with and without fast=true.

RiemannQuadratureLeft(T=Float64)

The left Riemann sum on $[0,1]$, i.e., the one-node rule

\[\int_0^1 f(x) \, dx \approx f(0) .\]

It is the interpolatory rule for the single node $c_1 = 0$ and therefore exact for constants only, giving order 1.

julia> RiemannQuadratureLeft()
QuadratureRule{Float64, 1}(1, [0.0], [1.0])

RiemannQuadratureRight(T=Float64)

The right Riemann sum on $[0,1]$, i.e., the one-node rule

\[\int_0^1 f(x) \, dx \approx f(1) .\]

It is the interpolatory rule for the single node $c_1 = 1$ and therefore exact for constants only, giving order 1.

julia> RiemannQuadratureRight()
QuadratureRule{Float64, 1}(1, [1.0], [1.0])

MidpointQuadrature(T=Float64)

The midpoint rule on $[0,1]$, i.e., the one-node rule

\[\int_0^1 f(x) \, dx \approx f \big( \tfrac{1}{2} \big) .\]

Although it uses a single node, the symmetry of the node about the centre of the interval makes it exact for linear functions as well, so its order is 2 rather than 1. It is the one-node Gauss-Legendre rule, cf. GaussLegendreQuadrature.

julia> MidpointQuadrature()
QuadratureRule{Float64, 1}(2, [0.5], [1.0])

TrapezoidalQuadrature(T=Float64)

The trapezoidal rule on $[0,1]$, i.e., the two-node rule

\[\int_0^1 f(x) \, dx \approx \tfrac{1}{2} \, f(0) + \tfrac{1}{2} \, f(1) .\]

It integrates the linear interpolant through the endpoints and is thus exact for linear functions, giving order 2. It is the two-node Lobatto-Legendre rule, cf. LobattoLegendreQuadrature.

julia> TrapezoidalQuadrature()
QuadratureRule{Float64, 2}(2, [0.0, 1.0], [0.5, 0.5])

GaussLegendreQuadrature(s; IT=BigFloat, fast=false)
GaussLegendreQuadrature(T, s; IT=_default_arithmetic(T), fast=false)

The Gauss-Legendre quadrature rule with s nodes on the interval $[0,1]$.

The nodes are the roots of the Legendre polynomial $P_s$ mapped to $[0,1]$. Being the free choice of both nodes and weights, they achieve the maximal possible degree of exactness $2s-1$, so the rule has order $2s$. All nodes lie in the interior of the interval and all weights are positive.

The weights are computed from the closed form

\[w_i = \frac{1}{P_s'(x_i)^2} \int_{-1}^{+1} \left( \frac{P_s(x)}{x - x_i} \right)^2 dx ,\]

where the integrand is the square of the (unnormalised) Lagrange basis polynomial associated with $x_i$, evaluated by exact polynomial division and integration in the arithmetic IT. The nodes and weights are finally shifted and scaled to $[0,1]$, and are also available on their own as gauss_legendre_nodes and gauss_legendre_weights.

Arguments

• T: element type of the resulting rule, Float64 if omitted.
• s: number of nodes.
• IT: arithmetic in which nodes and weights are computed. For a numeric T this defaults to BigFloat, so that the result is accurate to full precision in T; for any other T, in particular a symbolic one, it defaults to T itself, so that the rule is computed exactly, cf. QuadratureRules._default_arithmetic.
• fast: if true, take the nodes and weights directly from FastGaussQuadrature.gausslegendre, which is much faster but computes in double precision, so the result is only accurate to about Float64 precision.

julia> GaussLegendreQuadrature(2)
QuadratureRule{Float64, 2}(4, [0.2113248654051871, 0.7886751345948129], [0.5, 0.5])

julia> GaussLegendreQuadrature(3)(x -> x^5)   # exact up to degree 5
0.16666666666666669

See also LobattoLegendreQuadrature, gauss_legendre_nodes.

LobattoLegendreQuadrature(s; IT=BigFloat, fast=false)
LobattoLegendreQuadrature(T, s; IT=_default_arithmetic(T), fast=false)

The Lobatto-Legendre quadrature rule with s nodes on the interval $[0,1]$.

In contrast to GaussLegendreQuadrature, both endpoints of the interval are included among the nodes. This constraint costs two degrees of exactness: only the $s-2$ interior nodes are free, so the maximal degree of exactness is $2s-3$ and the rule has order $2s-2$. In exchange, the rule can be applied when the values at the endpoints are needed anyway, which is why it is common in finite element methods, collocation schemes and variational integrators.

The interior nodes are the roots of $P_{s-1}'$ and the weights are given in closed form by

\[w_i = \frac{2}{s \, (s-1) \, \big[ P_{s-1}(x_i) \big]^2} ,\]

which holds for the endpoints as well. Nodes and weights are computed in the arithmetic IT and then shifted and scaled to $[0,1]$, and are also available on their own as lobatto_legendre_nodes and lobatto_legendre_weights.

Throws an ErrorException for s == 1. See GaussLegendreQuadrature for the meaning of T, IT and fast.

julia> LobattoLegendreQuadrature(2) == TrapezoidalQuadrature()
true

julia> LobattoLegendreQuadrature(3)
QuadratureRule{Float64, 3}(4, [0.0, 0.5, 1.0], [0.16666666666666666, 0.6666666666666666, 0.16666666666666666])

The three-node rule is Simpson's rule, and it is exact up to degree $2 \cdot 3 - 3 = 3$.

RadauLegendreQuadrature(s, endpoint; IT=BigFloat, fast=false)
RadauLegendreQuadrature(T, s, ::Val{endpoint}; IT=_default_arithmetic(T), fast=false)

The Radau-Legendre quadrature rule with s nodes on the interval $[0,1]$.

It sits between GaussLegendreQuadrature, which prescribes no node, and LobattoLegendreQuadrature, which prescribes both endpoints: exactly one endpoint is included among the nodes, selected by endpoint as :left (the node 0) or :right (the node 1). This single constraint costs one degree of exactness, so the maximal degree of exactness is $2s-2$ and the rule has order $2s-1$. All weights are positive.

The one-sided constraint is what makes the Radau nodes the natural choice for stiffly accurate collocation: the Radau IA and Radau IIA Runge-Kutta methods are built on the :left and :right nodes respectively.

The free nodes are the roots of $P_{s-1} + P_s$ other than the prescribed endpoint, and the weights are given in closed form by

\[w_i = \frac{1 \mp x_i}{s^2 \, \big[ P_{s-1}(x_i) \big]^2} ,\]

with the upper sign for :left and the lower one for :right, valid at the prescribed endpoint as well. Nodes and weights are computed in the arithmetic IT and then shifted and scaled to $[0,1]$, and are also available on their own as radau_legendre_nodes and radau_legendre_weights.

Unlike the Lobatto rules, the Radau rules are defined for s == 1, where the single node is the prescribed endpoint and the rule reduces to a Riemann sum.

See GaussLegendreQuadrature for the meaning of T, IT and fast, and radau_legendre_nodes for endpoint.

julia> RadauLegendreQuadrature(1, :left) == RiemannQuadratureLeft()
true

julia> RadauLegendreQuadrature(1, :right) == RiemannQuadratureRight()
true

julia> RadauLegendreQuadrature(2, :right)
QuadratureRule{Float64, 2}(3, [0.3333333333333333, 1.0], [0.75, 0.25])

julia> RadauLegendreQuadrature(3, :right)(x -> x^4)   # exact up to degree 2*3-2
0.19999999999999996

See also GaussLegendreQuadrature and LobattoLegendreQuadrature.

ClenshawCurtisQuadrature(s; IT=BigFloat)
ClenshawCurtisQuadrature(T, s; IT=_default_arithmetic(T))

The Clenshaw-Curtis quadrature rule with s nodes on the interval $[0,1]$.

The nodes are the Chebyshev points of the second kind, which include both endpoints of the interval. The weights are those of the interpolatory rule through them: the integrand is expanded in a Chebyshev series, whose coefficients follow from the values at the nodes by a discrete cosine transform, and the series is integrated term by term. Since $\int_{-1}^{+1} T_{2j}(x) \, dx = -2/(4j^2-1)$ and the odd terms integrate to zero, this yields

\[w_k = \frac{c_k}{n} \left( 1 - \sum_{j=1}^{\lfloor n/2 \rfloor} \frac{b_j}{4 j^2 - 1} \cos ( j \vartheta_k ) \right) , \qquad \vartheta_k = \frac{2 \pi k}{n} , \qquad n = s-1 ,\]

with $c_k = 1$ at the endpoints and $2$ otherwise, and $b_j = 1$ for the last term of an even sum and $2$ otherwise; a further factor $1/2$ maps the rule to $[0,1]$. This is the explicit form given by Reid; see the Clenshaw-Curtis section of the manual for the derivation and for the literature. Nodes and weights are also available on their own as clenshaw_curtis_nodes and clenshaw_curtis_weights.

Being an interpolatory rule on s nodes, it is exact for polynomials of degree $\le s-1$. For odd s it gains one further degree, because the monomial of degree s that it would otherwise fail on is odd about the midpoint of the interval and its error therefore cancels. The order is s for even s and s+1 for odd s, and it is sharp: the rule is not exact one degree beyond. All weights are positive, which is what guarantees convergence for every continuous integrand.

Although its order is only about half that of GaussLegendreQuadrature with the same number of nodes, Clenshaw-Curtis converges at essentially the same rate for integrands that are not analytic in a sizable neighbourhood of the interval, and its nodes are given in closed form, requiring no root finding.

Throws an ErrorException for s == 1.

Arguments

• T: element type of the resulting rule, Float64 if omitted.
• s: number of nodes.
• IT: arithmetic in which nodes and weights are computed, BigFloat for a numeric T and T itself otherwise, cf. QuadratureRules._default_arithmetic.

julia> ClenshawCurtisQuadrature(3)     # Simpson's rule
QuadratureRule{Float64, 3}(4, [0.0, 0.5, 1.0], [0.16666666666666666, 0.6666666666666666, 0.16666666666666666])

julia> ClenshawCurtisQuadrature(3) == LobattoLegendreQuadrature(3)
true

julia> ClenshawCurtisQuadrature(5)(x -> x^4)
0.19999999999999996

See also LobattoChebyshevQuadrature, which is the same rule, and GaussChebyshevQuadrature, which is the analogous rule on the Chebyshev points of the first kind.

GaussChebyshevQuadrature(s; IT=BigFloat)
GaussChebyshevQuadrature(T, s; IT=_default_arithmetic(T))

The Gauss-Chebyshev quadrature rule with s nodes on the interval $[0,1]$, also known as Fejér's first rule.

The nodes are the Chebyshev points of the first kind and the weights are those of the interpolatory rule through them, obtained by integrating the Chebyshev interpolant of f term by term:

\[w_i = \frac{2}{s} \left( 1 - 2 \sum_{j=1}^{\lfloor s/2 \rfloor} \frac{\cos ( 2 j \theta_i )}{4 j^2 - 1} \right) , \qquad \theta_i = \frac{(2i-1) \, \pi}{2s} ,\]

scaled by $1/2$ for the interval $[0,1]$. Being an interpolatory rule on s nodes, it is exact for polynomials of degree $\le s-1$. For odd s it gains one further degree, because the monomial of degree s that it would otherwise fail on is odd about the midpoint of the interval and its error therefore cancels. The order is s for even s and s+1 for odd s, and it is sharp. All weights are positive.

Arguments

• T: element type of the resulting rule, Float64 if omitted.
• s: number of nodes.
• IT: arithmetic in which nodes and weights are computed, BigFloat for a numeric T and T itself otherwise, cf. QuadratureRules._default_arithmetic. As for ClenshawCurtisQuadrature, the weight sum costs $O(s^2)$ operations, so a lower working precision is faster but accumulates round-off in the intermediate terms; see the note there.

julia> quad = GaussChebyshevQuadrature(3);

julia> order(quad)          # odd s, so exact up to degree 3
4

julia> quad(x -> x^2)
0.33333333333333337

julia> GaussChebyshevQuadrature(1) == MidpointQuadrature()
true

See also ClenshawCurtisQuadrature and ChebyshevQuadrature.

LobattoChebyshevQuadrature(s; IT=BigFloat)
LobattoChebyshevQuadrature(T, s; IT=_default_arithmetic(T))

The Lobatto-Chebyshev quadrature rule with s nodes on the interval $[0,1]$.

Its nodes are the Chebyshev points of the second kind, which include both endpoints of the interval. These are precisely the Clenshaw-Curtis nodes, and since an interpolatory rule is uniquely determined by its nodes, this rule is the Clenshaw-Curtis rule: LobattoChebyshevQuadrature(s) == ClenshawCurtisQuadrature(s). The implementation therefore delegates to ClenshawCurtisQuadrature, where the weights and the resulting order, s for even s and s+1 for odd s, are documented.

Throws an ErrorException for s == 1.

julia> LobattoChebyshevQuadrature(3) == ClenshawCurtisQuadrature(3)
true

julia> LobattoChebyshevQuadrature(3)
QuadratureRule{Float64, 3}(4, [0.0, 0.5, 1.0], [0.16666666666666666, 0.6666666666666666, 0.16666666666666666])

ChebyshevQuadrature(s, kind; kwargs...)
ChebyshevQuadrature(T, s, ::Val{kind}; kwargs...)

Umbrella constructor for the Chebyshev quadrature rules with s nodes on $[0,1]$, dispatching on the kind of the Chebyshev points:

• kind = 1 gives GaussChebyshevQuadrature, i.e., Fejér's first rule,
• kind = 2 gives LobattoChebyshevQuadrature, i.e., the Clenshaw-Curtis rule.

Keyword arguments are forwarded to the selected rule; both accept IT, the arithmetic in which nodes and weights are computed. An unsupported keyword raises a MethodError rather than being silently ignored.

julia> ChebyshevQuadrature(4, 1) == GaussChebyshevQuadrature(4)
true

julia> ChebyshevQuadrature(4, 2) == LobattoChebyshevQuadrature(4)
true

julia> ChebyshevQuadrature(8, 1; IT=Float64) ≈ ChebyshevQuadrature(8, 1)
true

TanhSinhQuadrature(n; IT=BigFloat)
TanhSinhQuadrature(T, n; IT=BigFloat)

The tanh-sinh quadrature rule of level n on the interval $[0,1]$.

Also known as the double-exponential formula of Takahasi and Mori [4]. It is not an interpolatory rule on a prescribed set of nodes, but the trapezoidal rule applied after the change of variables

\[x = \tanh \left( \frac{\pi}{2} \sinh t \right) , \qquad t \in \mathbb{R} ,\]

which maps $\mathbb{R}$ onto $(-1,+1)$. Since $dx/dt$ decays like $\exp ( - \tfrac{\pi}{2} e^{|t|} )$, the transformed integrand vanishes double-exponentially at both ends whatever the integrand does at the endpoints, and the infinite trapezoidal sum can be truncated after a handful of terms. With $h = 2^{-n}$ the nodes and weights on $[0,1]$ are

\[c_k = \frac{1}{1 + e^{-\pi \sinh (k h)}} , \qquad b_k = \frac{\pi h}{4} \, \frac{\cosh (k h)}{\cosh^2 \big( \tfrac{\pi}{2} \sinh (k h) \big)} , \qquad k \in \mathbb{Z} ,\]

symmetric about the centre node $c_0 = 1/2$. See the Tanh-Sinh section of the manual for the derivation and for the literature.

The rule is exact for no polynomial, not even for a constant: the weights sum to one only up to the truncation error. Its order is therefore reported as 0, and it is the one rule in this package whose accuracy is not described by a polynomial degree of exactness. What it does instead is converge like $\mathcal{O} ( e^{-cN/\log N} )$ in the number of nodes $N$, which in practice means that each level roughly doubles the number of correct digits until the precision of T is reached. For Float64 that point is reached at level 3; a BigFloat rule at the default precision needs level 5.

Arguments

• T: element type of the resulting rule, Float64 if omitted. Unlike the other families, this rule accepts floating point types only, since the truncation below is defined in terms of rounding in T.
• n: the level, i.e. the number of halvings of the step size, so that $h = 2^{-n}$. This is not a node count — the number of nodes follows from the truncation below and grows like $2^n$, e.g. 13, 25, 51, 101, 203 for levels 1 to 5 at T=Float64.
• IT: arithmetic in which nodes and weights are computed, BigFloat by default. It must be a floating point type as well, as every candidate node and weight is rounded to T from it to test the truncation criterion.

The grid is truncated at the first k whose weight rounds to zero in T or whose node rounds to an endpoint — of either $[0,1]$ or $[-1,+1]$, so that neither representation degenerates. A pair that is indistinguishable from its predecessor in T has its weight folded into that predecessor instead of being added, which leaves the quadrature sum unchanged and keeps the nodes strictly increasing. Consequently no node ever coincides with an endpoint, and integrands that are singular there may be passed in directly.

Throws an ErrorException for n < 1 and an ArgumentError for a non-floating point T.

julia> quad = TanhSinhQuadrature(3);

julia> nnodes(quad)
51

julia> order(quad)
0

julia> quad(x -> log(x)) ≈ -1                     # ∫₀¹ log x dx, singular at x = 0
true

julia> quad(x -> exp(x)) ≈ exp(1) - 1
true

julia> isfinite(quad(x -> 1 / sqrt(x * (1 - x))))  # no node sits on an endpoint
true

See also tanh_sinh_nodes and tanh_sinh_weights for the nodes and weights alone, and GaussLegendreQuadrature, which is the better choice for an integrand that is smooth up to and including the endpoints.

QuadratureInterval

Supertype of the intervals on which nodes and weights can be requested, namely UnitInterval and SymmetricInterval.

Every node and weight function takes an interval keyword argument of this type, defaulting to UnitInterval(). It selects the interval the result lives on and nothing else: the two conventions describe the same quadrature rule, related by the affine map $c_i = (x_i + 1)/2$ and $b_i = w_i / 2$.

UnitInterval()

The interval $[0,1]$, on which the weights sum to $1$.

This is the reference interval of QuadratureRule and the default for every node and weight function. See QuadratureInterval and SymmetricInterval.

julia> gauss_legendre_nodes(2; interval = UnitInterval())
2-element Vector{Float64}:
 0.2113248654051871
 0.7886751345948129

SymmetricInterval()

The interval $[-1,+1]$, on which the weights sum to $2$.

This is the interval on which the classical theory is formulated, being where the Legendre and Chebyshev polynomials are defined, and the one the closed forms of most families are computed on. It is symmetric about the origin; the nodes need not be, and for a Radau rule deliberately are not. See QuadratureInterval and UnitInterval.

julia> gauss_legendre_nodes(2; interval = SymmetricInterval())
2-element Vector{Float64}:
 -0.5773502691896257
  0.5773502691896257

gauss_legendre_nodes(s; kwargs...)
gauss_legendre_nodes(T, s; IT=_default_arithmetic(T), interval=UnitInterval())

The s Gauss-Legendre nodes, i.e., the roots of the Legendre polynomial $P_s$ mapped to interval. They lie in the interior of the interval.

On the default UnitInterval they are equivalently the roots of the shifted Legendre polynomial

\[P_s (2x - 1) ,\]

which is the form in which they are usually stated in the Runge-Kutta literature, $2x-1$ being the inverse of the map to $[0,1]$.

The roots are computed on $[-1,+1]$ in the arithmetic IT, mapped to interval in that same arithmetic and converted to T only at the very end. For a numeric IT they are obtained by refining the double precision approximations of FastGaussQuadrature.gausslegendre with Newton's method, so that IT=BigFloat, the default whenever T is a numeric type, yields nodes that are accurate to full BigFloat precision independently of T. For any other IT, in particular a symbolic one, they are instead obtained exactly as the eigenvalues of the companion matrix of $P_s$; this is the default whenever T is not a numeric type, so that gauss_legendre_nodes(T, s) with a symbolic T returns exact expressions.

Arguments

• T: element type of the returned vector, Float64 if omitted.
• s: number of nodes.
• IT: arithmetic in which the roots are computed, BigFloat for a numeric T and T itself otherwise, cf. QuadratureRules._default_arithmetic.
• interval: UnitInterval for $[0,1]$, the default and the interval of GaussLegendreQuadrature, or SymmetricInterval for $[-1,+1]$, related by $c_i = (x_i + 1)/2$.

julia> gauss_legendre_nodes(2)
2-element Vector{Float64}:
 0.2113248654051871
 0.7886751345948129

julia> gauss_legendre_nodes(2; interval = SymmetricInterval())
2-element Vector{Float64}:
 -0.5773502691896257
  0.5773502691896257

gauss_legendre_weights(s; kwargs...)
gauss_legendre_weights(T, s; IT=_default_arithmetic(T), interval=UnitInterval())

The s Gauss-Legendre weights belonging to the nodes returned by gauss_legendre_nodes for the same interval. All of them are positive.

They are computed from the closed form

\[w_i = \frac{1}{P_s'(x_i)^2} \int_{-1}^{+1} \left( \frac{P_s(x)}{x - x_i} \right)^2 dx ,\]

where the integrand is the square of the (unnormalised) Lagrange basis polynomial associated with $x_i$, evaluated by exact polynomial division and integration in the arithmetic IT. Being formulated on $[-1,+1]$, where the weights sum to $2$, it is the SymmetricInterval weights that are primary here; those on $[0,1]$ are obtained as $b_i = w_i / 2$ and sum to $1$.

In terms of the UnitInterval nodes $c_i$ themselves the same weights read

\[b_i = \bigg( \frac{dP}{dx} (c_i) \bigg)^{-2} \int \limits_0^1 \bigg( \frac{P(x)}{x - c_i} \bigg)^2 dx , \qquad P(x) = P_s (2x-1) ,\]

the two forms agreeing because the substitution multiplies the integral by $2$ — four from the squared denominator, one half from $dx$ — and $P'^2$ by $4$, which leaves the halving. See gauss_legendre_nodes for the arguments.

julia> gauss_legendre_weights(2)
2-element Vector{Float64}:
 0.5
 0.5

julia> gauss_legendre_weights(2; interval = SymmetricInterval())
2-element Vector{Float64}:
 1.0
 1.0

lobatto_legendre_nodes(s; kwargs...)
lobatto_legendre_nodes(T, s; IT=_default_arithmetic(T), interval=UnitInterval())

The s Lobatto-Legendre nodes, i.e., the two endpoints of the interval together with the $s-2$ roots of $P_{s-1}'$ mapped to interval.

Rather than differentiating the Legendre polynomial, the interior nodes are obtained as the roots of the $(s-2)$-nd derivative of $(1-x^2)^{s-1}$, which has the same roots by Rodrigues' formula. That derivative has degree $s$, and its $s$ roots are therefore the whole node set, endpoints included. The roots are computed in the arithmetic IT, either Newton-refined from the double precision approximations of FastGaussQuadrature.gausslobatto or, for a non-numeric IT such as a symbolic one, exactly from the companion matrix.

On the default UnitInterval the corresponding polynomial is

\[\frac{d^{\,s-2}}{dx^{\,s-2}} \big( (x - x^2)^{s-1} \big) ,\]

since $x - x^2 = (1 - \xi^2)/4$ under $x = (\xi+1)/2$; this is the form in which the Lobatto nodes are usually stated in the Runge-Kutta literature.

The endpoints are set to exactly $\mp 1$ before the mapping, so the first and last node come out as exactly -1 and +1 on the symmetric interval and exactly 0 and 1 on the unit one.

Throws an ErrorException for s == 1, as a Lobatto rule needs at least the two endpoints. See gauss_legendre_nodes for the arguments.

julia> lobatto_legendre_nodes(3)
3-element Vector{Float64}:
 0.0
 0.5
 1.0

julia> lobatto_legendre_nodes(3; interval = SymmetricInterval())
3-element Vector{Float64}:
 -1.0
  0.0
  1.0

lobatto_legendre_weights(s; kwargs...)
lobatto_legendre_weights(T, s; IT=_default_arithmetic(T), interval=UnitInterval())

The s Lobatto-Legendre weights belonging to the nodes returned by lobatto_legendre_nodes for the same interval. All of them are positive.

They are given in closed form by

\[w_i = \frac{2}{s \, (s-1) \, \big[ P_{s-1}(x_i) \big]^2} ,\]

which holds for the interior nodes and for the two endpoints alike. Being formulated on $[-1,+1]$, where the weights sum to $2$, it is the SymmetricInterval weights that are primary here; those on $[0,1]$ are obtained as $b_i = w_i / 2$ and sum to $1$.

In terms of the UnitInterval nodes $c_j$ the halving cancels the numerator, so the same weights read

\[b_j = \frac{1}{s \, (s-1) \, \big[ P_{s-1}(2 c_j - 1) \big]^2} ,\]

which is the form in which they are usually tabulated for the Lobatto Runge-Kutta methods.

Throws an ErrorException for s == 1. See gauss_legendre_nodes for the arguments.

julia> lobatto_legendre_weights(3)
3-element Vector{Float64}:
 0.16666666666666666
 0.6666666666666666
 0.16666666666666666

julia> lobatto_legendre_weights(3; interval = SymmetricInterval())
3-element Vector{Float64}:
 0.3333333333333333
 1.3333333333333333
 0.3333333333333333

radau_legendre_nodes(s, endpoint; kwargs...)
radau_legendre_nodes(T, s, ::Val{endpoint}; IT=_default_arithmetic(T),
                     interval=UnitInterval())

The s Radau-Legendre nodes, i.e., one prescribed endpoint of the interval together with the $s-1$ free nodes that maximise the degree of exactness, mapped to interval.

Which endpoint is prescribed is selected by endpoint:

• :left includes the left end of the interval, the classical Gauss-Radau convention and the one underlying the Radau IA Runge-Kutta methods,
• :right includes the right end, the convention underlying the Radau IIA methods.

There is deliberately no default for endpoint: the two variants are not interchangeable, and silently choosing one would be easy to overlook.

The left nodes are the $s$ roots of $P_{s-1} + P_s$ on $[-1,+1]$, one of which is exactly $-1$. They are computed in the arithmetic IT, either Newton-refined from the double precision approximations of FastGaussQuadrature.gaussradau or, for a non-numeric IT such as a symbolic one, exactly from the companion matrix, and the prescribed endpoint is then set to exactly $-1$. The right nodes are obtained by reflection, $x \mapsto -x$, which makes the two variants exact mirror images of each other. As the prescribed endpoint is exact, so is the corresponding node after the mapping.

Note that only the interval of SymmetricInterval is symmetric: the Radau nodes are asymmetric by construction on either interval, which is the whole point of the family.

Equivalently, and as the Radau nodes are usually stated in the Runge-Kutta literature, on the default UnitInterval they are the roots of

\[\frac{d^{\,s-1}}{dx^{\,s-1}} \big( x^s (x - 1)^{s-1} \big) \qquad \text{for } \texttt{:left} , \qquad \frac{d^{\,s-1}}{dx^{\,s-1}} \big( x^{s-1} (x - 1)^s \big) \qquad \text{for } \texttt{:right} ,\]

each of degree $s$. These agree with $P_{s-1} + P_s$ through the Rodrigues formula for the Jacobi polynomial $P^{(0,1)}_{s-1}$,

\[P^{(0,1)}_{s-1} (x) \; \propto \; \frac{1}{1+x} \, \frac{d^{\,s-1}}{dx^{\,s-1}} \big( (1-x)^{s-1} (1+x)^{s} \big) ,\]

whose numerator is what the differentiated product above becomes under the map to $[-1,+1]$. That numerator is therefore proportional to $P_{s-1} + P_s$ itself, of degree $s$, and the division by $1+x$ is precisely what strips the prescribed endpoint off to leave the degree $s-1$ Jacobi polynomial whose roots are the free nodes. The two expressions are mirror images of one another under $x \mapsto 1-x$, as the node sets are.

Arguments

• T: element type of the returned vector, Float64 if omitted.
• s: number of nodes.
• endpoint: :left or :right, the endpoint included among the nodes.
• IT: arithmetic in which the roots are computed, BigFloat for a numeric T and T itself otherwise, cf. QuadratureRules._default_arithmetic.
• interval: UnitInterval for $[0,1]$, the default and the interval of RadauLegendreQuadrature, or SymmetricInterval for $[-1,+1]$.

julia> radau_legendre_nodes(2, :left)
2-element Vector{Float64}:
 0.0
 0.6666666666666666

julia> radau_legendre_nodes(2, :right)
2-element Vector{Float64}:
 0.3333333333333333
 1.0

julia> radau_legendre_nodes(2, :left; interval = SymmetricInterval())
2-element Vector{Float64}:
 -1.0
  0.3333333333333333

radau_legendre_weights(s, endpoint; kwargs...)
radau_legendre_weights(T, s, ::Val{endpoint}; IT=_default_arithmetic(T),
                       interval=UnitInterval())

The s Radau-Legendre weights belonging to the nodes returned by radau_legendre_nodes for the same endpoint and interval. All of them are positive.

They are given in closed form by

\[w_i = \frac{1 \mp x_i}{s^2 \, \big[ P_{s-1}(x_i) \big]^2} ,\]

with the upper sign for endpoint = :left and the lower one for endpoint = :right. Like the corresponding Lobatto formula this holds for the free nodes and for the prescribed endpoint alike, where $P_{s-1}(\mp 1)^2 = 1$ reduces it to $2/s^2$. Being formulated on $[-1,+1]$, where the weights sum to $2$, it is the SymmetricInterval weights that are primary here; those on $[0,1]$ are obtained as $b_i = w_i / 2$ and sum to $1$.

In terms of the UnitInterval nodes $c_i$ the same weights read

\[b_i = \frac{1 \mp (2 c_i - 1)}{2 \, s^2 \, \big[ P_{s-1}(2 c_i - 1) \big]^2} ,\]

with the signs as above, so that the prescribed endpoint again carries $1/s^2$.

See radau_legendre_nodes for the arguments.

julia> radau_legendre_weights(2, :left)
2-element Vector{Float64}:
 0.25
 0.75

julia> radau_legendre_weights(2, :left; interval = SymmetricInterval())
2-element Vector{Float64}:
 0.5
 1.5

chebyshev_nodes(s, kind; kwargs...)
chebyshev_nodes(T, s, ::Val{kind}; IT=_default_arithmetic(T), interval=UnitInterval())

The s Chebyshev nodes of the first (kind = 1) or second (kind = 2) kind, in ascending order, mapped to interval.

On $[-1,+1]$ the nodes of the first kind are the roots of the Chebyshev polynomial $T_s$,

\[x_i = \cos \left( \frac{(2i-1) \, \pi}{2s} \right) , \qquad i = 1, \dots, s ,\]

which lie strictly inside the interval. They are evaluated in the equivalent sin form used in the implementation, which is more accurate near the ends of the interval. The nodes of the second kind are the extrema of $T_{s-1}$,

\[x_i = \cos \left( \frac{(i-1) \, \pi}{s-1} \right) , \qquad i = 1, \dots, s ,\]

and include both endpoints.

The closed forms are evaluated on $[-1,+1]$ in the arithmetic IT, mapped to interval in that same arithmetic and converted to T only at the very end. For a numeric T the default is IT=BigFloat, so that the closed forms are evaluated to full precision whatever T is and the returned values are correctly rounded; see ClenshawCurtisQuadrature for the trade-off involved in choosing a lower working precision. For any other T the default is IT=T, so that with a symbolic T the closed forms are evaluated exactly, $\pi$ included; cf. QuadratureRules._default_arithmetic.

Nodes of the second kind require s ≥ 2 and throw an ErrorException otherwise.

Arguments

• T: element type of the returned vector, Float64 if omitted.
• s: number of nodes.
• kind: 1 or 2, cf. above.
• IT: arithmetic in which the closed forms are evaluated, BigFloat for a numeric T and T itself otherwise, cf. QuadratureRules._default_arithmetic.
• interval: UnitInterval for $[0,1]$, the default and the interval of ChebyshevQuadrature, or SymmetricInterval for $[-1,+1]$.

julia> chebyshev_nodes(3, 1)
3-element Vector{Float64}:
 0.06698729810778067
 0.5
 0.9330127018922193

julia> chebyshev_nodes(3, 1; interval = SymmetricInterval())
3-element Vector{Float64}:
 -0.8660254037844386
  0.0
  0.8660254037844386

chebyshev_weights(s, kind; kwargs...)
chebyshev_weights(T, s, ::Val{kind}; IT=_default_arithmetic(T), interval=UnitInterval())

The s interpolatory weights belonging to the Chebyshev nodes of the first (kind = 1) or second (kind = 2) kind, i.e. to chebyshev_nodes for the same kind and interval. All of them are positive.

For kind = 1 these are the weights of Fejér's first rule,

\[b_i = \frac{1}{s} \left( 1 - 2 \sum_{j=1}^{\lfloor s/2 \rfloor} \frac{\cos (2 j \vartheta_i)}{4 j^2 - 1} \right) , \qquad \vartheta_i = \frac{(2i-1) \pi}{2s} ,\]

already normalised to $[0,1]$, where they sum to $1$. For kind = 2 the nodes coincide with the Clenshaw-Curtis nodes, so the weights are those of clenshaw_curtis_weights, to which this function delegates. Both formulas carry the normalisation to $[0,1]$, so unlike the Legendre families it is the UnitInterval weights that are primary here; those on $[-1,+1]$ are obtained as $w_i = 2 b_i$ and sum to $2$.

See chebyshev_nodes for the arguments. kind = 2 requires s ≥ 2.

julia> chebyshev_weights(3, 1)
3-element Vector{Float64}:
 0.2222222222222222
 0.5555555555555556
 0.2222222222222222

julia> chebyshev_weights(3, 1; interval = SymmetricInterval())
3-element Vector{Float64}:
 0.4444444444444444
 1.1111111111111112
 0.4444444444444444

gauss_chebyshev_nodes(s; kwargs...)
gauss_chebyshev_nodes(T, s; IT=_default_arithmetic(T), interval=UnitInterval())

The s Gauss-Chebyshev nodes, i.e., the Chebyshev nodes of the first kind, cf. chebyshev_nodes.

On the unit interval these are the nodes of GaussChebyshevQuadrature.

gauss_chebyshev_weights(s; kwargs...)
gauss_chebyshev_weights(T, s; IT=_default_arithmetic(T), interval=UnitInterval())

The s Gauss-Chebyshev weights, i.e., the weights of Fejér's first rule, cf. chebyshev_weights.

On the unit interval these are the weights of GaussChebyshevQuadrature.

lobatto_chebyshev_nodes(s; kwargs...)
lobatto_chebyshev_nodes(T, s; IT=_default_arithmetic(T), interval=UnitInterval())

The s Lobatto-Chebyshev nodes, i.e., the Chebyshev nodes of the second kind, cf. chebyshev_nodes. They include both endpoints of the interval.

These nodes coincide with clenshaw_curtis_nodes, and correspondingly LobattoChebyshevQuadrature coincides with ClenshawCurtisQuadrature.

Requires s ≥ 2.

lobatto_chebyshev_weights(s; kwargs...)
lobatto_chebyshev_weights(T, s; IT=_default_arithmetic(T), interval=UnitInterval())

The s Lobatto-Chebyshev weights, cf. chebyshev_weights.

As the nodes coincide with the Clenshaw-Curtis nodes, so do the weights: these are identical to clenshaw_curtis_weights, and on the unit interval they are the weights of LobattoChebyshevQuadrature.

Requires s ≥ 2.

clenshaw_curtis_nodes(s; kwargs...)
clenshaw_curtis_nodes(T, s; IT=_default_arithmetic(T), interval=UnitInterval())

The s Clenshaw-Curtis nodes, i.e., the Chebyshev nodes of the second kind, cf. chebyshev_nodes. They include both endpoints of the interval.

Identical to lobatto_chebyshev_nodes. On the unit interval these are exactly the nodes of ClenshawCurtisQuadrature at the same T and IT. Requires s ≥ 2.

julia> clenshaw_curtis_nodes(5)
5-element Vector{Float64}:
 0.0
 0.14644660940672624
 0.5
 0.8535533905932737
 1.0

clenshaw_curtis_weights(s; kwargs...)
clenshaw_curtis_weights(T, s; IT=_default_arithmetic(T), interval=UnitInterval())

The s Clenshaw-Curtis weights, belonging to the nodes returned by clenshaw_curtis_nodes for the same interval. All of them are positive.

They are the weights of the interpolatory rule through the Chebyshev points of the second kind, given explicitly on $[0,1]$, where they sum to $1$, by

\[b_k = \frac{c_k}{2n} \left( 1 - \sum_{j=1}^{\lfloor n/2 \rfloor} \frac{b_j}{4 j^2 - 1} \cos ( j \vartheta_k ) \right) , \qquad \vartheta_k = \frac{2 \pi k}{n} , \qquad n = s-1 ,\]

with $c_k = 1$ at the endpoints and $2$ otherwise, and $b_j = 1$ for the last term of an even sum and $2$ otherwise. The formula carries the normalisation to $[0,1]$, so it is the UnitInterval weights that are primary here; those on $[-1,+1]$ are obtained as $w_i = 2 b_i$ and sum to $2$. See ClenshawCurtisQuadrature for the derivation, the literature and the arguments.

Throws an ErrorException for s == 1.

julia> clenshaw_curtis_weights(3)
3-element Vector{Float64}:
 0.16666666666666666
 0.6666666666666666
 0.16666666666666666

julia> clenshaw_curtis_weights(3; interval = SymmetricInterval())
3-element Vector{Float64}:
 0.3333333333333333
 1.3333333333333333
 0.3333333333333333

tanh_sinh_nodes(n; kwargs...)
tanh_sinh_nodes(T, n; IT=BigFloat, interval=UnitInterval())

The tanh-sinh nodes of level n, in ascending order, mapped to interval.

On $[-1,+1]$ they are the images

\[x_k = \tanh \left( \frac{\pi}{2} \sinh (k h) \right) , \qquad h = 2^{-n} ,\]

of the equidistant grid $t = k h$ under the tanh-sinh transformation. They lie strictly inside the interval and cluster double-exponentially fast towards its ends. On $[0,1]$ the transformation is the logistic sigmoid,

\[c_k = \frac{1 + x_k}{2} = \frac{1}{1 + e^{-\pi \sinh (k h)}} , \qquad h = 2^{-n} ,\]

which is the form the nodes are actually computed in, since it yields the small distance of the outer nodes from the endpoints without cancellation. Unlike every other family it is therefore the UnitInterval nodes that are primary here.

The number of nodes is not a free parameter: the grid is truncated where the nodes can no longer be distinguished from $\pm 1$ in the target precision, so it depends on n, on T and, through the weights, on IT. See TanhSinhQuadrature for the arguments and for the truncation criterion.

For the same reason T must be a floating point type here, whereas the other families also accept exact and symbolic ones: it is the rounding in T that decides where the infinite sum is cut off, so there is nothing to compute for a type that does not round. Anything else throws an ArgumentError.

On the unit interval these are exactly the nodes of TanhSinhQuadrature at the same n, T and IT.

julia> tanh_sinh_nodes(1)
13-element Vector{Float64}:
 2.1470805279391204e-14
 5.562167559007666e-9
 1.1261403769203567e-5
 0.0012425717713878065
 0.024316017963626528
 0.1628642538757821
 0.5
 0.8371357461242179
 0.9756839820363735
 0.9987574282286122
 0.9999887385962308
 0.9999999944378325
 0.9999999999999786

julia> x = tanh_sinh_nodes(1; interval = SymmetricInterval());

julia> x == -reverse(x)
true

tanh_sinh_weights(n; kwargs...)
tanh_sinh_weights(T, n; IT=BigFloat, interval=UnitInterval())

The tanh-sinh weights of level n, belonging to the nodes returned by tanh_sinh_nodes for the same interval, i.e. the trapezoidal weights $h$ times the Jacobian of the tanh-sinh transformation,

\[b_k = \frac{\pi h}{4} \, \frac{\cosh (k h)}{\cosh^2 \big( \tfrac{\pi}{2} \sinh (k h) \big)} , \qquad h = 2^{-n} .\]

They are symmetric about the centre weight bit for bit, and all of them are positive. The formula carries the normalisation to $[0,1]$, so it is the UnitInterval weights that are primary here; those on $[-1,+1]$ are obtained as $w_i = 2 b_i$.

Unlike every other family in this package these weights attain their sum — $1$ on the unit interval, $2$ on the symmetric one — only up to the truncation error, since the rule is not exact even for the constant. See TanhSinhQuadrature, whose weights these are, for the arguments and for the truncation criterion.

julia> b = tanh_sinh_weights(1);

julia> b == reverse(b)
true

julia> isapprox(sum(b), 1; atol = 1E-5)
true

julia> isapprox(sum(tanh_sinh_weights(1; interval = SymmetricInterval())), 2; atol = 1E-5)
true

The arithmetic in which nodes and weights are computed by default for the element type T.

Types from the numeric tower — floating point, rational, integer and complex ones — are computed in BigFloat and rounded to T at the very end, so that the result is accurate to full precision in T. Every other type is computed in itself, on the assumption that it does its own arithmetic exactly, so that a symbolic T yields nodes and weights in closed form.

Should that assumption not hold for some element type, IT=BigFloat recovers the behaviour of the numeric tower for it.

The Legendre polynomial $P_j(x)$ of degree $j$ on the interval $[-1,+1]$, evaluated by the three-term recurrence

\[j \, P_j (x) = (2j-1) \, x \, P_{j-1} (x) - (j-1) \, P_{j-2} (x) , \qquad P_0 = 1 , \quad P_1 = x .\]

The recurrence is used rather than the equivalent Rodrigues formula

\[P_j (x) = \frac{1}{j! \, 2^j} \, \frac{d^j}{dx^j} \big( x^2 - 1 \big)^j ,\]

because it costs $O(j)$ operations — it is carried upwards in a loop — instead of building and differentiating a polynomial of degree $2j$. The Rodrigues form is what identifies the $(s-2)$-nd derivative of $(1-x^2)^{s-1}$ used for the Lobatto nodes as an antiderivative of $P_{s-1}$; that its roots are the whole Lobatto node set is the separate, Jacobi-Rodrigues argument given in the manual, cf. lobatto_legendre_nodes.

Works for any x supporting arithmetic, including a Polynomial — see QuadratureRules._legendre_polynomial — and symbolic types.

Legendre polynomial P_s(x) of degree s on the interval [-1..+1].

Refine the approximate roots x₀ of the polynomial p with Newton's method.

The iteration is carried out in the arithmetic of p's coefficients and runs until the correction stops decreasing, i.e., until the roots are accurate to the working precision. This is used to compute quadrature nodes in arbitrary precision from double precision initial guesses. All roots are assumed to be real and simple.

The roots of the polynomial p, all assumed to be real and simple.

For coefficients from the numeric tower the double precision approximations x₀() are refined with Newton's method in the arithmetic of p. For any other coefficient type, in particular a symbolic one, the roots are computed exactly as the eigenvalues of the companion matrix; this needs eigvals for that type, which a computer algebra system supplies and which resolves the roots into radicals. In both cases they come out in no particular order, so the caller sorts them.

x₀ is a thunk rather than a vector so that the initial guess, which the exact branch has no use for, is only computed where it is actually needed.

Nodes and weights of the tanh-sinh rule of level n on the interval $[0,1]$, computed and returned in the arithmetic IT and truncated to what the target type T can resolve. This is the common back end of TanhSinhQuadrature, tanh_sinh_nodes and tanh_sinh_weights, so that all three truncate identically.

Substituting $x = \tanh ( \tfrac{\pi}{2} \sinh t )$ and applying the trapezoidal rule with step $h = 2^{-n}$ at $t = k h$ gives, after the move to $[0,1]$,

\[\delta_k = \frac{1 - \tanh u_k}{2} = \frac{1}{1 + e^{2 u_k}} , \qquad b_k = \frac{\pi h}{4} \, \frac{\cosh (k h)}{\cosh^2 u_k} , \qquad u_k = \frac{\pi}{2} \sinh (k h) ,\]

where $\delta_k$ is the distance of the outer pair of nodes from the endpoints, $c_{\pm k} = \tfrac{1}{2} \pm (\tfrac{1}{2} - \delta_k)$. It is written through the logistic form on the right to avoid the cancellation that $1 - \tanh u_k$ would suffer for large $u_k$ — and it is those small numbers that the rule is all about.

The pairs are collected for $k > 0$ and the rule is assembled symmetrically about the centre node $1/2$, so that the returned weights are symmetric bit for bit and the nodes ascending by construction.

Level n uses $h = 2^{-n}$, so the nodes of level n are contained in those of level n+1. The single loop over $k = 1, 2, \dots$ at the finest step is equivalent to refining level by level: the truncation criterion is a threshold in $t$, shared by all levels, so the union of the level grids is exactly the multiples of $2^{-n}$ below that threshold.

T must be a floating point type: unlike the other rules, tanh-sinh has no exact variant, as it is the rounding in T that decides where the infinite sum is cut off. So must IT, which is not checked, since every candidate is rounded from IT to T along the way.

Shift and scale nodes from the interval [-1,+1] to the interval [0,1].

Shift and scale nodes from the interval [0,1] to the interval [-1,+1].

Scale weights from the interval [-1,+1] to the interval [0,1].

Scale weights from the interval [0,1] to the interval [-1,+1].

Scale nodes and weights from the interval [-1,+1] to the interval [0,1].

Scale nodes and weights from the interval [0,1] to the interval [-1,+1].