{
"cells": [
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"ellipse_rule (generic function with 1 method)"
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Plots, ComplexPhasePortrait, Interact, ApproxFun, Statistics, SpecialFunctions, LinearAlgebra\n",
"gr();\n",
"\n",
"periodic_rule(n) = 2π/n*(0:(n-1)), 2π/n*ones(n)\n",
"\n",
"function circle_rule(n, r) \n",
" θ = periodic_rule(n)[1]\n",
" r*exp.(im*θ), 2π*im*r/n*exp.(im*θ)\n",
"end\n",
"\n",
"function ellipse_rule(n, a, b) \n",
" θ = periodic_rule(n)[1]\n",
" a*cos.(θ) + b*im*sin.(θ), 2π/n*(-a*sin.(θ) + im*b*cos.(θ))\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\def\\I{{\\rm i}}\n",
"\\def\\E{{\\rm e}}\n",
"\\def\\D{{\\rm d}}\n",
"$$\n",
"\n",
"# M3M6: Methods of Mathematical Physics\n",
"\n",
"Dr. Sheehan Olver\n",
"
\n",
"s.olver@imperial.ac.uk\n",
"\n",
"\n",
"# Lecture 4: Trapezium rule, Fourier series and Laurent series\n",
"\n",
"\n",
"This lecture we cover\n",
"\n",
"1. Periodic and complex Trapezium rule \n",
" - Application: Calculating matrix exponentials\n",
" - Gershorwin's theorem\n",
"1. Laurent series and Fourier series\n",
" - Application: Decay of Fourier coefficients"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Periodic and complex Trapezium rule\n",
"\n",
"Quadrature rules are pairs of _nodes_ $x_0,\\ldots,x_{n-1}$ and weights $w_0,\\ldots,w_{n-1}$ to approximate integrals\n",
"$$\n",
"\\int_a^b f(x) dx \\approx \\sum_{j=0}^{n-1} w_j f(x_j)\n",
"$$\n",
"In this lecture we construct quadrature rules on complex contours $\\gamma$ to approximate contour integrals.\n",
"\n",
"\n",
"The trapezium rule gives an easy approximation to integrals. On $[0,2\\pi)$ for periodic $f(\\theta)$, we have a simplified form:\n",
"\n",
"**Definition (Periodic trapezium rule)** The _periodic trapezium rule_ is the approximation\n",
"\n",
"$$\\int_0^{2 \\pi} f(\\theta) d \\theta \\approx {2\\pi \\over n} \\sum_{j=0}^{n-1} f(\\theta_k)$$\n",
"\n",
"for $\\theta_j = {2 \\pi j \\over n}$.\n",
"\n",
"The periodic trapezium rule is amazingly accurate for smooth, periodic functions:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f = θ -> 1/(2 + cos(θ))\n",
"\n",
"errs = [((x, w) = periodic_rule(n); abs(sum(w.*f.(x)) - sum(Fun(f, 0 .. 2π)))) for n = 1:45];\n",
"plot(errs.+eps(); yscale=:log10, title=\"exponential convergence of n-point trapezium rule\", legend=false, xlabel=\"n\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The accuracy in integration is remarkable as the trapezoidal interpolant does not accurately approximate $f$: we can see clear differences between the functions here:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n=20\n",
"(x, w) = periodic_rule(n)\n",
"plot(Fun(f, 0 .. 2π); label = \"integrand\")\n",
"plot!(x, f.(x); label = \"trapezium approximation\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Complex Trapezium rule\n",
"\n",
"We can use the map that defines a contour to construct an approximation to integrals over closed contour $\\gamma$:\n",
"\n",
"**Definition (Complex trapezium rule)** The _complex trapezium rule_ on a contour $\\gamma$ (mapped from $[0,2\\pi)$) is the approximation\n",
"\n",
"$$\\oint_\\gamma f(z) dz \\approx \\sum_{j=0}^{n-1} w_j f(z_j)$$\n",
"\n",
"for\n",
"\n",
"$$z_j = \\gamma(\\theta_j) \\qquad\\hbox{and}\\qquad w_j = {2\\pi \\over n} \\gamma'(\\theta_j)$$\n",
"\n",
"\n",
"*Example (Circle trapezium rule)* On a circle $\\{r e^{i \\theta} : 0 \\leq \\theta < 2 \\pi\\}$, we have \n",
"\n",
"$$\\oint_\\gamma f(z) dz \\approx \\sum_{j=0}^{n-1} w_j f(z_j)$$\n",
"\n",
"for $z_j = r e^{i \\theta_j}$ and $w_j = {2 \\pi i r \\over n} e^{i \\theta_j}$.\n",
"\n",
"Here we plot the quadrature points:"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ζ, w = circle_rule(20, 1.0)\n",
"scatter(ζ; title=\"quadrature points\", legend=false, ratio=1.0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Circle trapezium rule is surprisingly accurate for analytic functions:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-9.814371537686384e-14 - 1.3111040031432708e-14im"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ζ, w = circle_rule(20, 1.0)\n",
"\n",
"f = z -> cos(z)\n",
"\n",
"z = 0.1+0.2im\n",
"\n",
"sum(f.(ζ)./(ζ .- z).*w)/(2π*im) - f(z)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Application: Numerical differentiation\n",
"\n",
"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:"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f = z -> gamma(z)\n",
"fp = z -> gamma(z)polygamma(0,z) # exact derivative \n",
"x = 1.2\n",
"fp_fd = [(h=2.0^(-n); (f(x+h)-f(x))/h) for n = 1:50]\n",
"plot(abs.(fp_fd .- fp(x)); yscale=:log10, legend=false, title = \"error of finite-difference with h=2^(-n)\", xlabel=\"n\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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! "
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"trap_fp = [((ζ, w) = circle_rule(n, 0.5); \n",
" ζ .+= x; # circle around x\n",
" sum(f.(ζ)./(ζ .- x).^2 .*w)/(2π*im)) for n=1:50]\n",
"\n",
"plot(abs.(trap_fp .- fp.(x)); yscale=:log10, \n",
" title=\"Error of trapezium rule applied to Cauchy integral formula\", xlabel=\"n\", legend=false)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can take things further and use this to calculate the higher order derivatives, with some care taken for choosing the radius:"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-7.993605777301127e-15 + 3.675487826103639e-16im"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"k=100\n",
"r = 1.0k\n",
"g = Fun( ζ -> exp(ζ)/(ζ - 0.1)^(k+1), Circle(0.1,r))\n",
"factorial(1.0k)/(2π*im) * sum(g) - exp(0.1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[Bournemann 2011](https://www-m3.ma.tum.de/foswiki/pub/M3/Allgemeines/FolkmarBornemannPublications/FoCM_Stability_Cauchy_Integrals.pdf) investigates this further and optimizes the radius."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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:"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"θ = periodic_rule(100)[1]\n",
"plot(θ, real(f.(0.6exp.(im*θ) .+ x)./(0.5exp.(im*θ))))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*Example (Ellipse trapezium rule)* On an ellipse $\\{a \\cos \\theta + b i \\sin \\theta : 0 \\leq \\theta < 2 \\pi\\}$ we have \n",
"$$\\oint_\\gamma f(z) dz \\approx \\sum_{j=0}^{n-1} w_j f(z_j)$$\n",
"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)$.\n",
"\n",
"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\n",
"$$\n",
"f(z) = 1/(25z^2 + 1)\n",
"$$\n",
"which has poles at $\\pm \\I/5$. Using an ellipse allows us to design a contour that avoids these singularities:"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"scatter([1/5im,-1/5im]; label=\"singularities\")\n",
"θ = range(0; stop=2π, length=2000)\n",
"a = 2; b= 0.1\n",
"plot!(a * cos.(θ) + im*b * sin.(θ); label=\"ellipse\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Thus we can still use Cauchy's integral formula:"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = 0.1\n",
"f = z -> 1/(25z^2 + 1)\n",
"\n",
"f_ellipse = [((z, w) = ellipse_rule(n, a, b); sum(f.(z)./(z.-x).*w)/(2π*im)) for n=1:1000]\n",
"plot(abs.(f_ellipse .- f(x)); yscale=:log10, title=\"convergence of n-point ellipse approximation\", legend=false, xlabel=\"n\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Taylor series versus Cauchy integral formula\n",
"\n",
"\n",
"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$:"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f = z -> sqrt(z)\n",
"\n",
"function sqrtⁿ(n,z,z₀) \n",
" ret = sqrt(z₀)\n",
" c = 0.5/ret*(z-z₀)\n",
" for k=1:n\n",
" ret += c\n",
" c *= -(2k-1)/(2*(k+1)*z₀)*(z-z₀)\n",
" end\n",
" ret\n",
"end\n",
"\n",
"\n",
"z₀ = 0.3\n",
"n = 40\n",
"phaseplot(-2..2, -2..2, z -> sqrtⁿ.(n,z,z₀))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here we see that the approximation is vallid on the expected ellipse:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(ζ, w) = ellipse_rule(20, 4.0, 1.0);\n",
"ζ .= ζ .+ 4.5;\n",
"f_c = z -> sum(f.(ζ)./(ζ.-z).*w)/(2π*im)\n",
"phaseplot(-2..10, -2 .. 2, f_c)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-1.3987033753437572e-11 + 2.459960261444439e-16im"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(ζ, w) = ellipse_rule(100, 4.0, 1.0);\n",
"ζ .= ζ .+ 4.5;\n",
"f_c = z -> sum(f.(ζ)./(ζ.-z).*w)/(2π*im)\n",
"f_c(5.3) - sqrt(5.3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Application: calculating matrix exponentials\n",
"\n",
"Let $A \\in {\\mathbb C}^{d \\times d}$, ${\\mathbf u}_0 \\in {\\mathbb C}^d$ and consider the constant coefficient linear ODE\n",
"\n",
"$$\n",
" {\\mathbf u}'(t) = A {\\mathbf u}(t)\\qquad\\hbox{and}\\qquad {\\mathbf u}(0) = {\\mathbf u}_0(0)\n",
" $$\n",
" \n",
"The solution to this is given by the _matrix exponential_\n",
"$$\n",
" {\\mathbf u}(t) = \\exp(A t) {\\mathbf u}_0\n",
" $$\n",
" where the matrix exponential is defined by it's Taylor series:\n",
"$$\n",
" \\exp(A) = \\sum_{k=0}^\\infty {A^k \\over k!}\n",
" $$\n",
"This has stability problems, so a more convenient form is as follows:\n",
"\n",
"**Theorem (Cauchy integral representation for matrix exponential)**\n",
"Let $A \\in {\\mathbb C}^{n \\times n}$ be a diagonalizable matrix:\n",
"$$\n",
" A = V \\Lambda V^{-1} \\qquad\\hbox{for}\\qquad\\Lambda =\n",
" \\begin{pmatrix}\\lambda_1 \\cr &\\ddots \\cr && \\lambda_d\\end{pmatrix}\n",
"$$\n",
"Let $\\gamma$ be a contour that surrounds the spectrum of $A$. Then we have\n",
"$$\n",
" \\exp(A) = {1 \\over 2 \\pi i} \\oint_\\gamma e^z (z I - A)^{-1} dz\n",
" $$\n",
" \n",
"**Proof** Note by definition\n",
"$$\n",
"\\begin{align*}\n",
"\\exp(A) &= \\sum_{k=0}^\\infty {A^k \\over k!} = V \\sum_{k=0}^\\infty {\\Lambda^k \\over k!}V^{-1} = V \\begin{pmatrix} \\E^{\\lambda_1}\\\\ & \\ddots \\\\ && \\E^{\\lambda_d} \\end{pmatrix} V^{-1} \\\\\n",
"&= {1 \\over 2\\pi \\I} V \\begin{pmatrix} \\oint_\\gamma {\\E^\\zeta \\over \\zeta - \\lambda_1 } \\D\\zeta \\\\ & \\ddots \\\\ && \\oint_\\gamma {\\E^\\zeta \\over \\zeta - \\lambda_d } \\D\\zeta\\end{pmatrix} V^{-1} \n",
"\\end{align*}\n",
"$$\n",
"It was important here that $\\gamma$ surrounded all $\\zeta_j$.\n",
"\n",
"We now take out the integration from the matrix, the easiest way to see this is to apply the Trapezium rule approximation:\n",
"$$\n",
"\\begin{align*}\n",
"V\\begin{pmatrix} \\oint_\\gamma {\\E^\\zeta \\over \\zeta - \\lambda_1 } \\D\\zeta \\\\ & \\ddots \\\\ && \\oint_\\gamma {\\E^\\zeta \\over \\zeta - \\lambda_d } \\D\\zeta\\end{pmatrix} V^{-1} &= V\\lim_{n \\rightarrow \\infty} \\begin{pmatrix} \\sum_{j=1}^n {w_j \\E^\\zeta_j \\over \\zeta_j - \\lambda_1 } \\\\ & \\ddots \\\\ && \\sum_{j=1}^n {w_j \\E^\\zeta_j \\over \\zeta_j - \\lambda_d } \\end{pmatrix} V^{-1} \\\\\n",
"&= \\lim_{n \\rightarrow \\infty} V\\sum_{j=1}^n w_j \\E^\\zeta_j \\begin{pmatrix} {1 \\over \\zeta_j - \\lambda_1 } \\\\ & \\ddots \\\\ && {1 \\over \\zeta_j - \\lambda_d } \\end{pmatrix} V^{-1} \\\\\n",
"&= \\lim_{n \\rightarrow \\infty} \\sum_{j=1}^n w_j \\E^\\zeta_j V(I \\zeta_j - \\Lambda)^{-1}V^{-1} \\\\\n",
"&= \\oint_\\gamma \\E^\\zeta V(I \\zeta - \\Lambda)^{-1} V^{-1} \\D \\zeta \n",
"= \\oint_\\gamma \\E^\\zeta (I \\zeta - V\\Lambda V^{-1})^{-1} \\D \\zeta \\\\\n",
"& = \\oint_\\gamma \\E^\\zeta (I \\zeta - A)^{-1} \\D \\zeta \\\\\n",
"\\end{align*}\n",
"$$\n",
"\n",
"\n",
"■"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"_Demonstration_ we use this formula alongside the complex trapezium rule to calculate matrix exponentials. Begin by creating a random symmetric matrix (which only has real eigenvalues):"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Array{Float64,1}:\n",
" -1.1701131879499997\n",
" -0.8307847323168219\n",
" -0.2011088438314149\n",
" 1.7390356067088892\n",
" 2.1354582101565316"
]
},
"execution_count": 42,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = randn(5,5)\n",
"A = A + A'\n",
"λ = eigvals(A)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can now by hand create a circle that surrounds all the eigenvalues:"
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 49,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"z,w = circle_rule(100,maximum(abs.(λ))+1)\n",
"\n",
"plot(z)\n",
"scatter!(λ,zeros(5); label = \"eigenvalues of A\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here we wrap this up into a function `circle_exp` that calculates the matrix exponential:"
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4.240422855934442e-13"
]
},
"execution_count": 50,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function circle_exp(A, n, z₀, r)\n",
" z,w = circle_rule(n,r)\n",
" z .+= z₀\n",
"\n",
" ret = zero(A)\n",
" for j=1:n\n",
" ret += w[j]*exp(z[j])*inv(z[j]*I - A)\n",
" end\n",
"\n",
" ret/(2π*im)\n",
"end\n",
" \n",
"circle_exp(A, 100, 0, 8.0) -exp(A) |>norm"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this case, it is beneficial to use an ellipse:"
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.3594516845940306e-13"
]
},
"execution_count": 51,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function ellipse_exp(A, n, z₀, a, b)\n",
" z,w = ellipse_rule(n,a,b)\n",
" z .+= z₀\n",
"\n",
" ret = zero(A)\n",
" for j=1:n\n",
" ret += w[j]*exp(z[j])*inv(z[j]*I - A)\n",
" end\n",
" ret/(2π*im)\n",
"end\n",
"\n",
"\n",
"ellipse_exp(A, 50, 0, 8.0, 5.0) -exp(A) |>norm"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For matrices with large negative eigenvalues (For example, discretisations of the Laplacian), complex quadrature can lead to much better accuracy than Taylor series:"
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4.813209010780333e-8"
]
},
"execution_count": 52,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function taylor_exp(A,n)\n",
" ret = Matrix(I, size(A))\n",
" for k=1:n\n",
" ret += A^k/factorial(1.0k)\n",
" end\n",
" ret\n",
"end\n",
"\n",
"B = A - 20I\n",
"\n",
"taylor_exp(B, 200) -exp(B) |>norm"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can use an ellpise to surround the spectrum:"
]
},
{
"cell_type": "code",
"execution_count": 53,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 53,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"scatter(complex.(eigvals(B)))\n",
"plot!(ellipse_rule(50,8,5)[1] .- 20)"
]
},
{
"cell_type": "code",
"execution_count": 54,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"7.24987088647028e-22"
]
},
"execution_count": 54,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"norm(ellipse_exp(B, 50, -20.0, 8.0, 5.0) - exp(B))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Gershgorin circle theorem\n",
"\n",
"If we only know $A$, how do we know how big to make the contour? Gershgorin's circle theorem gives the answer:\n",
"\n",
"**Theorem (Gershgorin)** Let $A \\in {\\mathbb C}^{d \\times d}$ and define \n",
"$$R_k = \\sum_{j=1 \\atop j \\neq k}^d |a_{kj}| $$\n",
"Then \n",
"$$\n",
"\\rho(A) \\subset \\bigcup_{k=1}^d \\bar B(a_{kk}, R_k)\n",
"$$\n",
"where $\\bar B(z_0, r)$ is the closed disk of radius $r$ and $\\rho(A)$ is the set of eigenvalues.\n",
"\n",
"**Proof **\n",
"\n",
"\n",
"■\n",
"\n",
"_Demonstration_ Here we apply this to a particular matrix:"
]
},
{
"cell_type": "code",
"execution_count": 55,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3×3 Array{Int64,2}:\n",
" 1 2 3\n",
" 1 5 2\n",
" -4 1 6"
]
},
"execution_count": 55,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = [1 2 3; 1 5 2; -4 1 6]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following calculates the row sums:"
]
},
{
"cell_type": "code",
"execution_count": 56,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3×1 Array{Int64,2}:\n",
" 5\n",
" 3\n",
" 5"
]
},
"execution_count": 56,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"R = sum(abs.(A - Diagonal(diag(A))),dims=2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Gershgorin's theorem tells us that the spectrum lies in the union of the circles surrounding the diagonals:"
]
},
{
"cell_type": "code",
"execution_count": 57,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 57,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"drawcircle!(z0, R) = plot!(θ-> real(z0) + R[1]*cos(θ), θ-> imag(z0) + R[1]*sin(θ), 0, 2π, fill=(0,:red), α = 0.2, legend=false)\n",
"\n",
"λ = eigvals(A)\n",
"p = plot()\n",
"for k = 1:size(A,1)\n",
" drawcircle!(A[k,k], R[k])\n",
"end\n",
"scatter!(complex.(λ); label=\"eigenvalues\")\n",
"scatter!(complex.(diag(A)); label=\"diagonals\")\n",
"p"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can therefore use this to choose a contour big enough to surround all the circles. Here's a fairly simplistic construction for our case where everything is real:"
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 58,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"z₀ = (maximum(diag(A) .+ R) + minimum(diag(A) .- R)) /2 # average edges of circle\n",
"r = max(abs.(diag(A) .- R .- z₀)..., abs.(diag(A) .+ R .- z₀)...)\n",
"\n",
"plot!(Circle(z₀, r))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Fourier series and Laurent series \n",
"\n",
"**Definition (Fourier series)** On $[-\\pi, \\pi)$, _Fourier series_ is an expansion of the form\n",
"$$\n",
" g(\\theta) = \\sum_{k=-\\infty}^\\infty g_k e^{i k \\theta}\n",
"$$\n",
"where\n",
"$$\n",
"g_k = {1\\over 2 \\pi} \\int_{-\\pi}^\\pi g(\\theta) e^{-i k \\theta} d \\theta\n",
"$$\n",
"\n",
"**Definition (Laurent series)** In the complex plane, _Laurent series around $z_0$_ is an expansion of the form \n",
"$$ \n",
" f(z) = \\sum_{k=-\\infty}^\\infty f_k (z-z_0)^k \n",
"$$ \n",
"\n",
"**Lemma (Fourier series = Laurent series)** On a circle in the complex plane \n",
"$$\n",
" \\gamma_r = \\{z_0 + re^{i \\theta} : -\\pi \\leq \\theta < \\pi \\},$$ \n",
"Laurent series of $f(z)$ around $z_0$ converges for $z \\in C_r$ if the Fourier series of $g(\\theta) = f(z_0 + r e^{i \\theta})$ converges, and the coefficients are given by\n",
"\n",
"$$ \n",
"f_k = {g_k \\over r^k} = {1 \\over 2 \\pi i} \\oint_{\\gamma_r} {f(\\zeta) \\over (\\zeta - z_0)^{k+1}} d \\zeta.\n",
"$$\n",
"\n",
"**Proof** This follows immediately from the change of variables $z = r \\E^{\\I \\theta} + z_0$. ■"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Interestingly, analytic properties of $f$ can be used to show decaying properties in $g$:\n",
"\n",
"**Theorem (Decay in Fourier/Laurent coefficients)** Suppose $f(z)$ is analytic in a closed annulus $A_{r,R}$ around $0$:\n",
"$$A_r(z_0) = \\{z : r ≤ | z| ≤ R\\}$$\n",
"Then for all $k$\n",
"$$|f_k | \\leq M\\min\\left\\{{1 \\over R^k} , {1 \\over r^k}\\right\\}$$\n",
"where $M = \\sup_{z \\in A} |f(z)|$.\n",
"\n",
"**Proof**\n",
"This is a simple application of the ML lemma: \n",
"$$\n",
"|f_k| = {1 \\over 2 \\pi } \\left|\\oint_{\\gamma_1} {f(\\zeta) \\over \\zeta^{k+1}} \\D \\zeta\\right| = {1 \\over 2 \\pi}\\left|\\oint_{\\gamma_r} {f(\\zeta) \\over \\zeta^{k+1}} \\D \\zeta\\right| \\leq \\sup_{\\zeta \\in \\gamma_r} |f(\\zeta)| R^{-k} \\leq M R^{-k}\n",
"$$\n",
"which similar applies deforming to $\\gamma_r$.\n",
"■\n",
"\n",
"\n",
"\n",
"_Demonstration_ The Fourier coefficients of\n",
"$$g(\\theta) = {1 \\over 2 - \\cos \\theta}$$\n",
"satisfies for $k \\geq 0$\n",
"$$\n",
" |g_k| \\leq {2 \\over 4 - R -R^{-1}} R^{-k}\n",
"$$\n",
"for all $R \\leq 2 + \\sqrt{3}$:"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"