2025-11-10 Transformed Quadrature¶

Singular integrals and Tanh-Sinh quadrature
Finite element integration and mapped elements
Adaptive integration
Integration in multiple dimensions

[16]:

using LinearAlgebra
using Plots
using Polynomials
default(lw=4, ms=5, legendfontsize=12, xtickfontsize=12, ytickfontsize=12)

# Vandermonde with Legendre polynomials
function vander_legendre(x, n=nothing)
    if isnothing(n)
        n = length(x) # Square by default
    end
    m = length(x)
    Q = ones(m, n)
    if n > 1
        Q[:, 2] = x
    end
    for k in 1:n-2
        Q[:, k+2] = ((2*k + 1) * x .* Q[:, k+1] - k * Q[:, k]) / (k + 1)
    end
    Q
end

# And a utility for points distributed via cos
CosRange(a, b, n) = (a + b)/2 .+ (b - a)/2 * cos.(LinRange(-pi, 0, n))

# Some test functions
F_expx(x) = exp(2x) / (1 + x^2)
f_expx(x) = 2*exp(2x) / (1 + x^2) - 2x*exp(2x)/(1 + x^2)^2

F_dtanh(x) = tanh(x)
f_dtanh(x) = cosh(x)^-2

integrands = [f_expx, f_dtanh]
antiderivatives = [F_expx, F_dtanh]
tests = zip(integrands, antiderivatives)

# Plotting utils for accuracy
function plot_accuracy(fint, tests, ns; ref=[1,2])
    a, b = -2, 2
    p = plot(xscale=:log10, yscale=:log10, xlabel="n", ylabel="error")
    for (f, F) in tests
        Is = [fint(f, a, b, n=n) for n in ns]
        Errors = abs.(Is .- (F(b) - F(a)))
        scatter!(ns, Errors, label=f)
    end
    for k in ref
        plot!(ns, ns.^(-1. * k), label="\$n^{-$k}\$")
    end
    p
end

[16]:

plot_accuracy (generic function with 1 method)

FastGaussQuadrature.jl ¶

There are packages to help us here.

[2]:

using Pkg
pkg"add FastGaussQuadrature"
using FastGaussQuadrature

n = 10
x, q = gausslegendre(n)
scatter(x, q, label="Gauss-Legendre", ylabel="weight", xlims=(-1, 1))
scatter!(gausslobatto(n)..., label="Gauss-Lobatto")

    Updating registry at `~/.julia/registries/General.toml`
   Resolving package versions...
  No Changes to `~/.julia/environments/v1.11/Project.toml`
  No Changes to `~/.julia/environments/v1.11/Manifest.toml`

[2]:

Other reading: Trefethen, Six Myths of Polynomial Interpolation and Quadrature

[21]:

@time gausslegendre(10_000_000);

  0.418783 seconds (16 allocations: 305.187 MiB, 64.74% gc time)

Singular Integrands¶

What should we do when the integrand contains a singularity?

[22]:

plot([sqrt log x->.5*x^(-.5)], xlim=(0, 2), ylim=(-2, 2), label=["\$sqrt(x)\$" "\$log(x)\$" "\$1/(2 sqrt(x))\$"])

[22]:

Error is higher than we would expect with other functions.

[23]:

function fint_gauss(f, a, b, n)
    x, w = gausslegendre(n)
    x = (a+b)/2 .+ (b-a)/2*x
    w *= (b - a)/2
    w' * f.(x)
end

[23]:

fint_gauss (generic function with 1 method)

[24]:

plot(3:4:100,
    n -> abs(fint_gauss(x -> .5*x^(-.5), 0, 1, n) - 1),
    marker=:auto, yscale=:log10, xscale=:log10, label="error")
plot!(n -> 1/n, label="\$1/n\$")

[24]:

Tanh-Sinh quadrature¶

When functions have singularities near the endpoints, it is usually more efficient to integrate via a change of variables. Suppose we have a strictly monotone differentiiable function \(\phi : \left( -\infty, \infty \right) \rightarrow \left( -1, 1 \right)\). Then with \(x = \phi \left( s \right)\), our integral transforms as

\[\int_{-1}^1 f \left( x \right) \, dx = \int_{-\infty}^\infty f \left( \phi \left( s \right) \right) \phi ' \left( s \right) \, ds\]

The tanh-sinh method uses a transformation such that \(\phi ' \left( s \right) \rightarrow 0\) faster than the singularity \(f \left( \phi \left( s \right) \right)\) grows, such that the integrand goes to \(0\) at finite \(s\).

[25]:

tanhsinh(s) = tanh(π/2*sinh(s))

# Here's the transformation
function dtanhsinh(s)
    ds = 1
    t = π/2 * sinh(s)
    dt = π/2 * cosh(s) * ds
    (1 - tanh(t)^2) * dt
end

# And lets plot the function and integrands
p = plot([tanhsinh], color=:black, label="tanhsinh(s)",
    xlims=(-3, 3),
    xlabel="s", title="tanh-sinh function and integrands")
for f in integrands
    plot!(s -> f(tanhsinh(s))*dtanhsinh(s), label="$f ∘ tanhsinh")
end
p

[25]:

Implementation¶

The function below implements tanh-sinh quadrature on the interval \(\left( -1, 1 \right)\). Given the number of points, we need to choose both the limits of integration (we can’t afford to integrate all the way to infinity) and the spacing. Here we make an arbitrary choice to integrate on the interval \(\left( -L, L \right)\) where \(L = \log \left( n \right)\). The grid spacing thus scales as \(h \approx 2 \log \left( n \right) / n\).

Modify the quadrature so it can be used to integrate on an arbitrary interval \(\left( a, b \right)\).

[26]:

function fint_tanhsinh(f, a, b; n=9)
    L = log(n)
    h = 2 * L / (n - 1)
    s = LinRange(-L, L, n)
    x = tanhsinh.(s)
    w = h * dtanhsinh.(s)
    ## Challenge: modify the weights w and points x to integrated on (a, b), not (-1, 1)
    w' * f.(x)
end

[26]:

fint_tanhsinh (generic function with 1 method)

[27]:

plot_accuracy(fint_tanhsinh, tests, 9:4:60, ref=[2,3,4])
plot!(xscale=:identity,)

[27]:

[10]:

# If you complete the challenge above
@assert fint_tanhsinh(log, 0, 1, n=20) ≈ -1
println("Tests pass")

DomainError with -0.9999999999999508:
log was called with a negative real argument but will only return a complex result if called with a complex argument. Try log(Complex(x)).

Stacktrace:
  [1] throw_complex_domainerror(f::Symbol, x::Float64)
    @ Base.Math ./math.jl:33
  [2] _log
    @ ./special/log.jl:295 [inlined]
  [3] log(x::Float64)
    @ Base.Math ./special/log.jl:261
  [4] _broadcast_getindex_evalf
    @ ./broadcast.jl:678 [inlined]
  [5] _broadcast_getindex
    @ ./broadcast.jl:651 [inlined]
  [6] getindex
    @ ./broadcast.jl:610 [inlined]
  [7] macro expansion
    @ ./broadcast.jl:973 [inlined]
  [8] macro expansion
    @ ./simdloop.jl:77 [inlined]
  [9] copyto!
    @ ./broadcast.jl:972 [inlined]
 [10] copyto!
    @ ./broadcast.jl:925 [inlined]
 [11] copy(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(log), Tuple{Vector{Float64}}})
    @ Base.Broadcast ./broadcast.jl:897
 [12] materialize
    @ ./broadcast.jl:872 [inlined]
 [13] fint_tanhsinh(f::Function, a::Int64, b::Int64; n::Int64)
    @ Main ./In[8]:8
 [14] top-level scope
    @ In[10]:2

Finite Element methods¶

Need to integrate the product of basis functions over every element
Choose a quadrature (e.g. Gauss) to minimize number of points needed for sufficient accuracy
Bad effects if insufficient points

Rough/singularities in middle of domain¶

If a function has rough or singular behavior in the interior of the domain, we will get low accuracy results

[11]:

using Pkg
pkg"add SpecialFunctions"
using SpecialFunctions

f_heaviside(x) = 1.0 * (x > 0)
F_heaviside(x) = max(x, 0)
f_mollify(x, sigma=.02) = 1/(sigma * sqrt(2*π)) * exp(-.5*(x/sigma)^2)
F_mollify(x, sigma=.02) = .5*erf(x/(sigma*sqrt(2)))
rough_tests = [(f_heaviside, F_heaviside), (f_mollify, F_mollify)]
plot([f_heaviside, f_mollify], xlim=(-1, 1), label=["heavyside" "mollify"])

   Resolving package versions...
  No Changes to `~/.julia/environments/v1.11/Project.toml`
  No Changes to `~/.julia/environments/v1.11/Manifest.toml`

[11]:

[12]:

plot_accuracy(fint_tanhsinh, rough_tests, 5:4:800, ref=[2,3,4]); plot!(legend=:bottomleft)

[12]:

Adaptive integration¶

Adaptive integration is a process where we adaptively subdivide the integration into smaller pieces until the numerical integral converges.

[28]:

using Pkg
pkg"add HCubature"
using HCubature

@show F_heaviside(1 - .3) - F_heaviside(-1 - .3)
hquadrature(x -> f_heaviside(x - 0.3), -1, 1, maxevals=1000)

   Resolving package versions...
  No Changes to `~/.julia/environments/v1.11/Project.toml`
  No Changes to `~/.julia/environments/v1.11/Manifest.toml`

F_heaviside(1 - 0.3) - F_heaviside(-1 - 0.3) = 0.7

[28]:

(0.700000008169331, 9.264907000048497e-9)

[14]:

hquadrature(f_mollify, -1, 1, maxevals=1000)

[14]:

(1.0000000000000002, 6.773756033519694e-10)

[15]:

F_mollify(1) - F_mollify(-1)

[15]:

1.0

Quadrature/cubature for multiple dimensions¶

Here is one such package: quadpy