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