In [37]:
using Plots, ComplexPhasePortrait, Interact, ApproxFun, Statistics, SpecialFunctions, LinearAlgebra
gr();

periodic_rule(n) = 2π/n*(0:(n-1)), 2π/n*ones(n)

function circle_rule(n, r) 
    θ = periodic_rule(n)[1]
    r*exp.(im*θ), 2π*im*r/n*exp.(im*θ)
end

function ellipse_rule(n, a, b) 
    θ = periodic_rule(n)[1]
    a*cos.(θ) + b*im*sin.(θ), 2π/n*(-a*sin.(θ) + im*b*cos.(θ))
end


Out[37]:
ellipse_rule (generic function with 1 method)
$$ \def\I{{\rm i}} \def\E{{\rm e}} \def\D{{\rm d}} $$

M3M6: Methods of Mathematical Physics

Dr. Sheehan Olver
s.olver@imperial.ac.uk

Lecture 4: Trapezium rule, Fourier series and Laurent series

This lecture we cover

  1. Periodic and complex Trapezium rule
    • Application: Calculating matrix exponentials
    • Gershorwin's theorem
  2. Laurent series and Fourier series
    • Application: Decay of Fourier coefficients

Periodic and complex Trapezium rule

Quadrature rules are pairs of nodes $x_0,\ldots,x_{n-1}$ and weights $w_0,\ldots,w_{n-1}$ to approximate integrals $$ \int_a^b f(x) dx \approx \sum_{j=0}^{n-1} w_j f(x_j) $$ In this lecture we construct quadrature rules on complex contours $\gamma$ to approximate contour integrals.

The trapezium rule gives an easy approximation to integrals. On $[0,2\pi)$ for periodic $f(\theta)$, we have a simplified form:

Definition (Periodic trapezium rule) The periodic trapezium rule is the approximation

$$\int_0^{2 \pi} f(\theta) d \theta \approx {2\pi \over n} \sum_{j=0}^{n-1} f(\theta_k)$$

for $\theta_j = {2 \pi j \over n}$.

The periodic trapezium rule is amazingly accurate for smooth, periodic functions:


In [2]:
f = θ ->  1/(2 + cos(θ))

errs = [((x, w) = periodic_rule(n); abs(sum(w.*f.(x)) - sum(Fun(f, 0 .. 2π)))) for n = 1:45];
plot(errs.+eps(); yscale=:log10, title="exponential convergence of n-point trapezium rule", legend=false, xlabel="n")


Out[2]:
0 10 20 30 40 10 - 15 10 - 10 10 - 5 10 0 exponential convergence of n-point trapezium rule n

The accuracy in integration is remarkable as the trapezoidal interpolant does not accurately approximate $f$: we can see clear differences between the functions here:


In [3]:
n=20
(x, w) = periodic_rule(n)
plot(Fun(f, 0 .. 2π); label = "integrand")
plot!(x, f.(x); label = "trapezium approximation")


Out[3]:
0 1 2 3 4 5 6 0.4 0.5 0.6 0.7 0.8 0.9 1.0 integrand trapezium approximation

Complex Trapezium rule

We can use the map that defines a contour to construct an approximation to integrals over closed contour $\gamma$:

Definition (Complex trapezium rule) The complex trapezium rule on a contour $\gamma$ (mapped from $[0,2\pi)$) is the approximation

$$\oint_\gamma f(z) dz \approx \sum_{j=0}^{n-1} w_j f(z_j)$$

for

$$z_j = \gamma(\theta_j) \qquad\hbox{and}\qquad w_j = {2\pi \over n} \gamma'(\theta_j)$$

Example (Circle trapezium rule) On a circle $\{r e^{i \theta} : 0 \leq \theta < 2 \pi\}$, we have

$$\oint_\gamma f(z) dz \approx \sum_{j=0}^{n-1} w_j f(z_j)$$

for $z_j = r e^{i \theta_j}$ and $w_j = {2 \pi i r \over n} e^{i \theta_j}$.

Here we plot the quadrature points:


In [27]:
ζ, w = circle_rule(20, 1.0)
scatter(ζ; title="quadrature points", legend=false, ratio=1.0)


Out[27]:
-1.0 -0.5 0.0 0.5 1.0 -1.0 -0.5 0.0 0.5 1.0 quadrature points Re(x) Im(x)

The Circle trapezium rule is surprisingly accurate for analytic functions:


In [5]:
ζ, w = circle_rule(20, 1.0)

f = z -> cos(z)

z = 0.1+0.2im

sum(f.(ζ)./(ζ .- z).*w)/(2π*im) - f(z)


Out[5]:
-9.814371537686384e-14 - 1.3111040031432708e-14im

Application: Numerical differentiation

Calculating high-order derivatives using limits is numerically unstable. Here is a demonstration using finite-difference: making $h$ small does not increase the accuracy after a certain point:


In [30]:
f = z -> gamma(z)
fp = z -> gamma(z)polygamma(0,z) # exact derivative 
x = 1.2
fp_fd = [(h=2.0^(-n);  (f(x+h)-f(x))/h) for n = 1:50]
plot(abs.(fp_fd .- fp(x)); yscale=:log10, legend=false, title = "error of finite-difference with h=2^(-n)", xlabel="n")


Out[30]:
0 10 20 30 40 50 10 - 8 10 - 6 10 - 4 10 - 2 error of finite-difference with h=2^(-n) n

But the previous formula tells us that we can reduce a derivative to a contour integral. The example above shows that it's still numerically unstable, but we can deform the integration contour, to make it stable!


In [31]:
trap_fp = [((ζ, w) = circle_rule(n, 0.5); 
        ζ .+= x; # circle around x
       sum(f.(ζ)./(ζ .- x).^2 .*w)/(2π*im)) for n=1:50]

plot(abs.(trap_fp .- fp.(x)); yscale=:log10, 
       title="Error of trapezium rule applied to Cauchy integral formula", xlabel="n", legend=false)


Out[31]:
0 10 20 30 40 50 10 - 15 10 - 10 10 - 5 10 0 Error of trapezium rule applied to Cauchy integral formula n

We can take things further and use this to calculate the higher order derivatives, with some care taken for choosing the radius:


In [32]:
k=100
r = 1.0k
g = Fun( ζ -> exp(ζ)/(ζ - 0.1)^(k+1), Circle(0.1,r))
factorial(1.0k)/(2π*im) * sum(g) - exp(0.1)


Out[32]:
-7.993605777301127e-15 + 3.675487826103639e-16im

Bournemann 2011 investigates this further and optimizes the radius.

The exponential convergence of the complex trapezium rule is a consequence of $f(\gamma(t))$ being 2π-periodic, as we will see later in a few lectures:


In [33]:
θ = periodic_rule(100)[1]
plot(θ, real(f.(0.6exp.(im*θ) .+ x)./(0.5exp.(im*θ))))


Out[33]:
0 1 2 3 4 5 6 -3 -2 -1 0 1 2 y1

Example (Ellipse trapezium rule) On an ellipse $\{a \cos \theta + b i \sin \theta : 0 \leq \theta < 2 \pi\}$ we have $$\oint_\gamma f(z) dz \approx \sum_{j=0}^{n-1} w_j f(z_j)$$ for $z_j = a \cos \theta_j + b i \sin \theta_j$ and $w_j = {2 \pi \over n} (-a \sin \theta_j + i b \cos \theta_j)$.

We can use the ellipse trapezium rule in place of the circle trapzium rule and still achieve accurate results. This gives us flexibility in avoiding singularities. Consider $$ f(z) = 1/(25z^2 + 1) $$ which has poles at $\pm \I/5$. Using an ellipse allows us to design a contour that avoids these singularities:


In [34]:
scatter([1/5im,-1/5im]; label="singularities")
θ = range(0; stop=2π, length=2000)
a = 2; b= 0.1
plot!(a * cos.(θ) + im*b * sin.(θ); label="ellipse")


Out[34]:
-2 -1 0 1 2 -0.2 -0.1 0.0 0.1 0.2 Re(x) Im(x) singularities ellipse

Thus we can still use Cauchy's integral formula:


In [35]:
x = 0.1
f = z -> 1/(25z^2 + 1)

f_ellipse = [((z, w) = ellipse_rule(n, a, b); sum(f.(z)./(z.-x).*w)/(2π*im)) for n=1:1000]
plot(abs.(f_ellipse .- f(x)); yscale=:log10, title="convergence of n-point ellipse approximation", legend=false, xlabel="n")


Out[35]:
0 250 500 750 1000 10 - 15 10 - 10 10 - 5 10 0 convergence of n-point ellipse approximation n

Taylor series versus Cauchy integral formula

The Taylor series gives a polynomial approximation to $f$. The Cauchy's integral formula discretisation gives a rational approximation, which is more adaptiple and does not require knowning the derivatives of $f$:


In [36]:
f = z -> sqrt(z)

function sqrtⁿ(n,z,z₀) 
    ret = sqrt(z₀)
    c = 0.5/ret*(z-z₀)
    for k=1:n
        ret += c
        c *= -(2k-1)/(2*(k+1)*z₀)*(z-z₀)
    end
    ret
end


z₀ = 0.3
n = 40
phaseplot(-2..2, -2..2, z -> sqrtⁿ.(n,z,z₀))


Out[36]:
-2 -1 0 1 2 -2 -1 0 1 2

Here we see that the approximation is vallid on the expected ellipse:


In [13]:
(ζ, w) = ellipse_rule(20, 4.0, 1.0);
ζ  .= ζ .+ 4.5;
f_c = z ->  sum(f.(ζ)./(ζ.-z).*w)/(2π*im)
phaseplot(-2..10, -2 .. 2, f_c)


Out[13]:
-2 0 2 4 6 8 10 -2 -1 0 1 2