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]:
Dr. Sheehan Olver
s.olver@imperial.ac.uk
This lecture we cover
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]:
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]:
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]:
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]:
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]:
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]:
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]:
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]:
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]:
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]:
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]:
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]: