The first 11 roots of the Airy function in Float64 precision https://mathworld.wolfram.com/AiryFunctionZeros.html

Examples

julia> zeros = airy.(FastGaussQuadrature.AIRY_ROOTS);

julia> all(zeros .< 1e-14)
true

First twenty roots of Bessel function $J_0$ in Float64. https://mathworld.wolfram.com/BesselFunctionZeros.html

Examples

julia> zeros = besselj0.(FastGaussQuadrature.BESSELJ0_ROOTS);

julia> all(zeros .< 1e-14)
true

Values of Bessel function $J_1$ on first ten roots of Bessel function J_0.

Examples

julia> roots = approx_besselroots(0,10);

julia> (besselj1.(roots)).^2 ≈ FastGaussQuadrature.BESSELJ1_ON_BESSELJ0_ROOTS
true

10-point Chebyshev type 2 integration matrix computed using Chebfun.

Used for numerical integration in boundary asymptotics for Gauss-Jacobi.

10-point Chebyshev type 2 differentiation matrix computed using Chebfun.

Used for numerical differentiation in boundary asymptotics for Gauss-Jacobi.

Coefficients of Chebyshev series approximations for the zeros of the Bessel functions

\[j_{\nu, s} \approx \sum_{k}^{n} C_{k,s} T_k(\frac{\nu-2}{3})\]

where $j_{\nu, s}$ is a $s$-th zero of Bessel function $J_\nu$, $\{T_k\}$ are Chebyshev polynomials and $C_{k, s}$ is the coefficients.

approx_besselroots(ν::Real, n::Integer) -> Vector{Float64}

Return the first $n$ roots of Bessel function. Note that this function is only 12-digits of precision.

\[J_{\nu}(x) = \sum_{m=0}^{\infty}\frac{(-1)^j}{\Gamma(\nu+j+1)j!} \left(\frac{x}{2}\right)^{2j+\nu}\]

Examples

julia> ν = 0.3;

julia> roots = approx_besselroots(ν, 10);

julia> zeros = (x -> besselj(ν, x)).(roots);

julia> all(zeros .< 1e-12)
true

Evaluate the orthonormal associated Laguerre polynomial with positive leading coefficient, as well as its derivative, in the point x using the recurrence relation.

Evaluate the interior asymptotic formula at x = cos(t). Assumption:

• length(t) == n ÷ 2

Evaluate the boundary asymptotic formula at x = cos(t). Assumption:

• length(t) == n ÷ 2

gausschebyshevt(n::Integer) -> x, w  # nodes, weights

Return nodes x and weights w of Gauss-Chebyshev quadrature of the 1st kind.

\[\int_{-1}^{1} \frac{f(x)}{\sqrt{1-x^2}} dx \approx \sum_{i=1}^{n} w_i f(x_i)\]

Examples

julia> x, w = gausschebyshevt(3);

julia> f(x) = x^4;

julia> I = dot(w, f.(x));

julia> I ≈ 3π/8
true

gausschebyshevu(n::Integer) -> x, w  # nodes, weights

Return nodes x and weights w of Gauss-Chebyshev quadrature of the 2nd kind.

\[\int_{-1}^{1} f(x)\sqrt{1-x^2} dx \approx \sum_{i=1}^{n} w_i f(x_i)\]

Examples

julia> x, w = gausschebyshevu(3);

julia> f(x) = x^4;

julia> I = dot(w, f.(x));

julia> I ≈ π/16
true

gausschebyshevv(n::Integer) -> x, w  # nodes, weights

Return nodes x and weights w of Gauss-Chebyshev quadrature of the 3rd kind.

\[\int_{-1}^{1} f(x)\sqrt{\frac{1+x}{1-x}} dx \approx \sum_{i=1}^{n} w_i f(x_i)\]

Examples

julia> x, w = gausschebyshevv(3);

julia> f(x) = x^4;

julia> I = dot(w, f.(x));

julia> I ≈ 3π/8
true

gausschebyshevw(n::Integer) -> x, w  # nodes, weights

Return nodes x and weights w of Gauss-Chebyshev quadrature of the 4th kind.

\[\int_{-1}^{1} f(x)\sqrt{\frac{1-x}{1+x}} dx \approx \sum_{i=1}^{n} w_i f(x_i)\]

Examples

julia> x, w = gausschebyshevw(3);

julia> f(x) = x^4;

julia> I = dot(w, f.(x));

julia> I ≈ 3π/8
true

gausshermite(n::Integer; normalize = false) -> x, w  # nodes, weights

Return nodes x and weights w of Gauss-Hermite quadrature.

\[\int_{-\infty}^{+\infty} f(x) \exp(-x^2) dx \approx \sum_{i=1}^{n} w_i f(x_i)\]

The option normalize=true instead computes

\[\int_{-\infty}^{+\infty} f(x) \phi(x) dx \approx \sum_{i=1}^{n} w_i f(x_i),\]

where $\phi$ is the standard normal density function.

Examples

julia> x, w = gausshermite(3);

julia> f(x) = x^4;

julia> I = dot(w, f.(x));

julia> I ≈ 3(√π)/4
true

gaussjacobi(n::Integer, α::Real, β::Real) -> x, w  # nodes, weights

Return nodes x and weights w of Gauss-Jacobi quadrature for exponents α and β.

\[\int_{-1}^{1} f(x) (1-x)^\alpha(1+x)^\beta dx \approx \sum_{i=1}^{n} w_i f(x_i)\]

Examples

julia> x, w = gaussjacobi(3, 1/3, -1/3);

julia> f(x) = x^4;

julia> I = dot(w, f.(x));

julia> I ≈ 268π/729(√3)
true

gausslaguerre(n::Integer, α::Real) -> x, w  # nodes, weights

Return nodes x and weights w of generalized Gauss-Laguerre quadrature.

\[\int_{0}^{+\infty} f(x) x^\alpha e^{-x} dx \approx \sum_{i=1}^{n} w_i f(x_i)\]

Examples

julia> x, w = gausslaguerre(3, 1.0);

julia> f(x) = x^4;

julia> I = dot(w, f.(x));

julia> I ≈ 120
true

Optionally, a reduced quadrature rule can be computed. In that case, only those points and weights are computed for which the weight does not underflow in the floating point precision type. Supply the optional argument reduced = true.

Though the code is generic, heuristical choices on the choice of the algorithm are based on achieving machine precision accuracy only for Float64 type. In case the default choice of algorithm does not achieve the desired accuracy, the user can manually invoke the following routines:

• gausslaguerre_GW: computation based on Golub-Welsch
• gausslaguerre_rec: computation based on Newton iterations applied to evaluation using the recurrence relation
• gausslaguerre_asy: the asymptotic expansions

gausslaguerre(n::Integer) -> x, w  # nodes, weights

Return nodes x and weights w of Gauss-Laguerre quadrature.

\[\int_{0}^{+\infty} f(x) e^{-x} dx \approx \sum_{i=1}^{n} w_i f(x_i)\]

Examples

julia> x, w = gausslaguerre(3);

julia> f(x) = x^4;

julia> I = dot(w, f.(x));

julia> I ≈ 24
true

Calculate Gauss-Laguerre nodes and weights from the eigenvalue decomposition of the Jacobi matrix.

Compute the Gauss-Laguerre rule using explicit asymptotic expansions for the nodes and weights. Optional parameters are:

• reduced: compute a reduced quadrature rule, discarding all points and weights as soon as the weights underflow
• T: the order of the expansion. Set T=-1 to determine the order adaptively depending on the size of the terms in the expansion
• recompute: if a crude measure of the error is larger than a tolerance, the point and weight are recomputed using the (slower) recursion+newton approach, yielding more reliable accurate results.

Compute Gauss-Laguerre rule based on the recurrence relation, using Newton iterations on an initial guess.

gausslegendre(n::Integer) -> x, w  # nodes, weights

Return nodes x and weights w of Gauss-Legendre quadrature.

\[\int_{-1}^{1} f(x) dx \approx \sum_{i=1}^{n} w_i f(x_i)\]

Examples

julia> x, w = gausslegendre(3);

julia> f(x) = x^4;

julia> I = dot(w, f.(x));

julia> I ≈ 2/5
true

gausslobatto(n::Integer) -> x, w  # nodes, weights

Return nodes x and weights w of Gauss-Lobatto quadrature.

\[\int_{-1}^{1} f(x) dx \approx \sum_{i=1}^{n} w_i f(x_i)\]

Examples

julia> x, w = gausslobatto(4);

julia> f(x) = x^4;

julia> I = dot(w, f.(x));

julia> I ≈ 2/5
true

Note that the both ends of nodes are fixed at -1 and 1.

julia> x, w = gausslobatto(4);

julia> x[1], x[end]
(-1.0, 1.0)

gaussradau(n::Integer) -> x, w  # nodes, weights

Return nodes x and weights w of Gauss-Radau quadrature.

\[\int_{-1}^{1} f(x) dx \approx \sum_{i=1}^{n} w_i f(x_i)\]

Examples

julia> x, w = gaussradau(3);

julia> f(x) = x^4;

julia> I = dot(w, f.(x));

julia> I ≈ 2/5
true

Note that the first node is fixed at -1.

julia> x, w = gaussradau(3);

julia> x[1]
-1.0