# 2022-04-15 Transformed Quadrature¶

## Last time¶

• Polynomial interpolation for integration

## Today¶

• Singular integrals and Tanh-Sinh quadrature

• Finite element integration and mapped elements

• Integration in multiple dimensions

using LinearAlgebra
using Plots
default(linewidth=4, legendfontsize=12)

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

CosRange(a, b, n) = (a + b)/2 .+ (b - a)/2 * cos.(LinRange(-pi, 0, n))

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

F_rsqrt(x) = 2 * sqrt(x + 1)
f_rsqrt(x) = 1 / sqrt(x + 1)

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

function plot_accuracy(fint, tests, ns; ref=[1,2])
a, b = -1, 1
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  plot_accuracy (generic function with 1 method)  # FastGaussQuadrature.jl¶ using FastGaussQuadrature n = 5 x, q = gausslegendre(n) scatter(x, q, label="Gauss-Legendre", ylabel="weight", xlims=(-1, 1)) scatter!(gausslobatto(n)..., label="Gauss-Lobatto") Trefethen, Six Myths of Polynomial Interpolation and Quadrature @time gausslegendre(1000000);   0.029118 seconds (10 allocations: 22.888 MiB, 31.66% gc time)  # Singular integrands¶ plot([sqrt log x->.5*x^(-.5)], xlim=(0, 2), ylim=(-2, 2)) 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  fint_gauss (generic function with 1 method)  plot(2:20, n -> abs(fint_gauss(x -> .5*x^(-.5), 0, 1, n) - 1), marker=:auto, yscale=:log10) plot!(n -> 1/n) # Tanh-Sinh quadrature: make everything smooth¶ 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 differentiable function $$\phi: (-\infty, \infty) \to (-1, 1)$$. Then with $$x = \phi(s)$$, our integral transforms as $\int_{-1}^1 f(x) \mathrm dx = \int_{-\infty}^\infty f(\phi(s)) \phi'(s) \mathrm d s .$ The tanh-sinh method uses a transformation such that $$\phi'(s) \to 0$$ faster than the singularity $$f(\phi(s))$$ grows, such that the integrand goes to 0 at finite $$s$$. tanhsinh(s) = tanh(pi/2*sinh(s)) function dtanhsinh(s) ds = 1 t = pi/2 * sinh(s) dt = pi/2 * cosh(s) * ds (1 - tanh(t)^2) * dt end 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 # Implementation¶

The function below implements tanh-sinh quadrature on the interval $$(-1,1)$$. 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 $$(-L, L)$$ where $$L = \log n$$. The grid spacing thus scales as $$h \approx 2 \log n / n$$.

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

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

fint_tanhsinh (generic function with 1 method)

plot_accuracy(fint_tanhsinh, tests, 9:4:60, ref=[2,3,4]) # If you complete the challenge above
@assert fint_tanhsinh(log, 0, 1, n=20) ≈ -1
println("Tests pass")

Tests pass