{
"cells": [
{
"cell_type": "code",
"execution_count": 53,
"metadata": {},
"outputs": [],
"source": [
"using Plots, ComplexPhasePortrait, ApproxFun, SingularIntegralEquations, DifferentialEquations, LinearAlgebra\n",
"gr();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# M3M6: Methods of Mathematical Physics 2018\n",
"\n",
"$$\n",
"\\def\\dashint{{\\int\\!\\!\\!\\!\\!\\!-\\,}}\n",
"\\def\\infdashint{\\dashint_{\\!\\!\\!-\\infty}^{\\,\\infty}}\n",
"\\def\\D{\\,{\\rm d}}\n",
"\\def\\E{{\\rm e}}\n",
"\\def\\dx{\\D x}\n",
"\\def\\dt{\\D t}\n",
"\\def\\dz{\\D z}\n",
"\\def\\C{{\\mathbb C}}\n",
"\\def\\R{{\\mathbb R}}\n",
"\\def\\CC{{\\cal C}}\n",
"\\def\\HH{{\\cal H}}\n",
"\\def\\I{{\\rm i}}\n",
"\\def\\qqqquad{\\qquad\\qquad}\n",
"\\def\\qqfor{\\qquad\\hbox{for}\\qquad}\n",
"\\def\\qqwhere{\\qquad\\hbox{where}\\qquad}\n",
"\\def\\Res_#1{\\underset{#1}{\\rm Res}}\\,\n",
"\\def\\sech{{\\rm sech}\\,}\n",
"\\def\\vc#1{{\\mathbf #1}}\n",
"$$\n",
"\n",
"Dr. Sheehan Olver\n",
"
\n",
"s.olver@imperial.ac.uk\n",
"\n",
"
\n",
"Website: https://github.com/dlfivefifty/M3M6LectureNotes\n",
"\n",
"# Solution Sheet 1\n",
"\n",
"\n",
"----------\n",
"\n",
"## Problem 1.1\n",
"\n",
"### 1.\n",
"Use fundamental theorem of algebra: a polynomial is a constant times a product of terms like $z - \\lambda_k$, where $\\lambda_k$ are the roots. In this case, the roots are $a$ times the quartic-root of $-1$, hence this gives us:\n",
"\\begin{align*}\n",
"z^4 + a^4 = (z-a \\E^{\\I \\pi /4})(z-a \\E^{3\\I \\pi /4})(z-a \\E^{5\\I \\pi /4})(z-a \\E^{7\\I \\pi /4}) \n",
"\\end{align*}\n",
"We are only interested in the root $a \\E^{\\I \\pi /4}$, thus we simplify the expression\n",
"Therefore \n",
"$$\n",
"{z^3 \\sin z \\over z^4 + a^4} = {z^3 \\sin z \\over (z-a \\E^{3\\I \\pi /4})(z-a \\E^{5\\I \\pi /4})(z-a \\E^{7\\I \\pi /4})} {1 \\over z-a \\E^{\\I \\pi /4}} = {a^3\\E^{3\\I \\pi /4} \\sin(a \\E^{\\I \\pi /4}) \\over a^3 (\\E^{\\I \\pi / 4} - \\E^{3\\I \\pi /4})(\\E^{\\I \\pi / 4} - \\E^{5\\I \\pi /4})(\\E^{\\I \\pi / 4} - \\E^{7\\I \\pi /4})} {1 \\over z-a \\E^{\\I \\pi /4}} + O(1)\n",
"$$\n",
"Therefore, \n",
"$$\n",
"\\Res_{z = a \\E^{\\I \\pi /4}} {z^3 \\sin z \\over z^4 + a^4} = { \\E^{3\\I \\pi /4} \\sin(a \\E^{\\I \\pi /4}) \\over \\E^{3\\I \\pi /4} (1 - \\E^{\\I \\pi \\over 2})(1 - \\E^{\\I \\pi })(1 - \\E^{3\\I \\pi/2 })}= { \\sin(a \\E^{\\I \\pi /4}) \\over (1 - \\I) (2) (1+ \\I)}= { \\sin(a \\E^{\\I \\pi /4}) \\over 4}\n",
"$$\n",
"\n",
"Let's check our work: we compare the numerically calculated residue to the formula we have derived:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.5378838853348212 + 0.07544036746694016im, 0.5378838853348215 + 0.0754403674669402im)"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a = 2.0\n",
"γ = Circle(a*exp(im*π/4), 0.1)\n",
"f = Fun(z -> z^3*sin(z)/(z^4+a^4), γ)\n",
"sum(f)/(2π*im), sin(a*exp(im*π/4))/4"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### 2. \n",
"\n",
"We have \n",
"$$\n",
"(z^2-1)^2 = (z-1)^2(z+1)^2\n",
"$$\n",
"Thus this is a slightly more challenging since it has a double pole. But we can expand using Geometric series:\n",
"$$\n",
" {z+1 \\over (z^2-1)^2} = {1 \\over (z-1)^2} {1 \\over 2 - (1-z)} = {1 \\over (z-1)^2} {1 \\over 2} (1 + (1-z)/2 +O(1-z)^2) = {1 \\over 2(z-1)^2} - {1 \\over 4 (z-1)} + O(1)\n",
"$$\n",
"Thus the residue is the negative-first Laurent coefficient, namely $-{1 \\over 4}$.\n",
"\n",
"We again check our work:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-0.2500000000000023 - 1.1916485920484316e-16im"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"γ = Circle(1, 0.1)\n",
"f = Fun(z -> (z+1)/(z^2-1)^2, γ)\n",
"sum(f)/(2π*im) # almost equals -1/4"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### 3. \n",
"\n",
"$$ {z^2 \\E^z \\over z^3 - a^3} = {z^2 \\E^z \\over (z-a) (z^2 + az + a^2)}\n",
"$$\n",
"We thus need only evaluate the extra term at $z=a$:\n",
"$$\n",
"\\Res_{z = a} {z^2 \\E^z \\over z^3 - a^3} = {\\E^a \\over 3}\n",
"$$\n",
"\n",
"Let's check:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(2.4630186996435506 + 4.738982805534393e-16im, 2.46301869964355)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a = 2.0\n",
"\n",
"γ = Circle(a, 0.1)\n",
"f = Fun(z -> z^2*exp(z)/(z^3-a^3), γ)\n",
"sum(f)/(2π*im), exp(a)/3"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 1.2\n",
"\n",
"### 1.\n",
"Change of variables $z = \\E^{\\I \\theta}$, $\\dz = \\I \\E^{\\I \\theta} \\D\\theta = \\I z \\D\\theta$, $\\cos \\theta = {z + z^{-1} \\over 2}$ gives\n",
"$$\n",
"\\int_0^{2 \\pi} {\\D\\theta \\over 5-4 \\cos \\theta} = - \\I \\oint {\\dz \\over 5z - 2z^2 - 2} = \\I \\oint {\\dz \\over (z-2)(2z-1)} = -\\pi \\Res_{z=1/2} {1 \\over (z-2)(z-1/2)} = {2 \\over 3} \\pi \n",
"$$\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"8.881784197001252e-16"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"θ = Fun(0 .. 2π)\n",
"sum(1/(5-4cos(θ))) -2π/3"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 2. \n",
"Use $\\cos 2 \\theta = { \\E^{2 \\I \\theta} + \\E^{-2 \\I \\theta} \\over 2} = {z^2 + z^{-2} \\over 2}$ to get\n",
"$$\n",
"\\int_0^{2 \\pi} {\\cos 2 \\theta \\D\\theta \\over 5+4 \\cos \\theta} = - {\\I \\over 4} \\oint { (z^4 + 1)\\dz \\over z^2 (z+1/2)(z+2)} = {\\pi \\over 2} \\left( \\Res_{z=-1/2} + \\Res_{z=0}\\right) {z^4+1 \\over z^2(z+2)(z+1/2)} = {\\pi \\over 6}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.1666666666666669"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"θ = Fun(0 .. 2π)\n",
"sum(cos(2θ)/(5-4cos(θ)))/π"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 3. \n",
"\n",
"Because the integrand is analytic and $O(z^{-2})$ in the upper half plane, we can use the residue theorem in the upper half plane using\n",
"$$\n",
" {1 \\over (z^2+1) (z^2 + 4)} = {1 \\over (z+\\I)(z-\\I) (z+2\\I)(z-2\\I)}\n",
"$$\n",
"\n",
"This has two poles in the upper half plane:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"phaseplot(-3..3, -3..4, z-> 1/((z^2+1)*(z^2+4)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\\begin{align*}\n",
"\t\\int_{-\\infty}^\\infty {1 \\over (x^2+1)(x^2+4)} \\dx &=\n",
" 2 \\pi \\I \\left( \\Res_{z = \\I} + \\Res_{z = 2\\I}\\right){1 \\over (z^2+1)(z^2+4)} \\\\\n",
" & = 2 \\pi \\I \\left( {1 \\over 2 \\I 3 \\I (-\\I)} + {1 \\over 3 \\I \\I 4\\I} \\right) = \\pi/ 6\n",
"\\end{align*}\n",
"\n",
"We can check the result numerically:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.523598775598299"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = Fun( Line())\n",
"sum(1/((x^2+1)*(x^2+4)))"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.5235987755982988"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"π/6"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 4. \n",
"\n",
"Again, decays like $O(z^{-2})$ in upper half plane so we can use residue calculus. This integrand has poles at $z = \\I$ and $z = 3 \\I$:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"phaseplot(-4..4, -4..4, z-> (z^2 - z + 2) / (z^4 + 10z^2 +9))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The residues are $(-1-\\I)/16$ and $(3-7\\I)/48$ giving the answer\n",
"$$\n",
"5\\pi \\over 12\n",
"$$\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 5.\n",
"$$\n",
"\t\\int_{-\\infty}^\\infty {1 \\over x + \\I } \\dx\n",
"$$\n",
"\n",
"Trick question: it's undefined because the integral doesn't decay fast enough. But what if I had asked for \n",
"$$\n",
"\t\\dashint_{-\\infty}^\\infty {1 \\over x + \\I } \\dx?\n",
"$$\n",
"\n",
"We can't use residue theorem since it doesn't decay fast enough, but we can use, with a contour $C_R = \\{ R \\E^{\\I theta} : 0 \\leq \\theta \\leq \\pi \\}$\n",
"$$\n",
"\\oint_{[-M, M] \\cup C_R} {1 \\over z + \\I} = 0\n",
"$$\n",
"Further, by direct substitution, we have \n",
"$$\n",
"\\int_{C_R} {1 \\over z + \\I}\\dz = \\I \\int_0^\\pi R {\\E^{\\I \\theta} \\over R \\E^{\\I \\theta} + \\I} \\D \\theta \n",
"$$\n",
"Letting $R \\rightarrow \\infty$, the integrand tends to one uniformly hence \n",
"$$\n",
" \\int_{C_R} {1 \\over z + \\I}\\dz \\rightarrow \\I \\int_0^\\pi \\D \\theta = \\I \\pi.\n",
"$$\n",
"Therefore, we have\n",
"$$\n",
"\t\\dashint_{-\\infty}^\\infty {1 \\over x + \\I } \\dx = \\lim_{M\\rightarrow \\infty} \\int_{-M}^M {1 \\over x + \\I} \\dx = - \\I \\pi .\n",
" $$\n",
" \n",
" Indeed:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-8.881784197001252e-16 - 3.13959265425659im"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = Fun(-1000 .. 1000)\n",
"sum(1/(x+im))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 6.\n",
"$$\n",
"\\int_{-\\infty}^\\infty {\\sin 2x \\over x^2 + x + 1} \\dx = \\Im \\int_{-\\infty}^\\infty {\\E^{ 2 \\I x} \\over x^2 + x + 1} \\dx \n",
"$$\n",
"This can be deformed in the upper half plane with a pole at ${-1 + \\I \\sqrt{3} \\over 2}$, using residue calculus gives us\n",
"$$\n",
" -{2 \\pi \\over \\sqrt 3} {\\sin 1 \\over \\E^\\sqrt{3}}\n",
"$$\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-0.5400548830723747"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = Fun(-100 .. 100)\n",
"sum(sin(2x)/(x^2+x+1))"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-0.5400553569742235"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"-2π/sqrt(3) * sin(1)/exp(sqrt(3))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"### 7.\n",
"$$\n",
"\t\\int_{-\\infty}^\\infty {\\cos x \\over x^2 + 4} \\dx = \\Re \t\\int_{-\\infty}^\\infty {\\E^{\\I x}\\over x^2 + 4} \\dx\n",
"$$\n",
"and residue calculus gives ${\\pi \\over 2 \\E^2}$\n"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.21254026836702827"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M = 200\n",
"x = Fun(-M .. M)\n",
"sum(cos(x)/(x^2+4)) # converges if we make M even bigger"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.21258416579381814"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"π/(2*exp(2))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"### 8.\n",
"$$\n",
"\t\\int_{-\\infty}^\\infty {x \\sin x \\over x^2 + 1} \\dx = { \\pi \\over \\E}\n",
" $$\n",
"using Residue calculus. You need to appeal to Jordan's lemma to argue that it can still be done even with only $O(x^{-1})$ decay."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.1560943440671996"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M = 2000\n",
"x = Fun(-M .. M)\n",
"sum(x*sin(x)/(x^2+1)) # Converges if we make M even bigger"
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.1557273497909217"
]
},
"execution_count": 49,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"π/ℯ"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 9. \n",
"$$\n",
"\\int_{-\\infty}^\\infty {\\cos a x - \\cos b x \\over x^2} \\dx \\qqwhere a,b > 0\n",
"$$\n",
"\n",
"We have for $f(x) = {\\E^{\\I a x} - \\E^{\\I b x} \\over x^2}$\n",
"$$\n",
"\\Re f(x) = {\\cos a x - \\cos b x \\over x^2}\n",
"$$\n",
"\n",
"Note that, since $\\cos x = 1 + x^2/2 + O(x^4)$, the integrand is fine near zero:\n",
"$$\n",
"{\\cos a x - \\cos b x \\over x^2} = {(a -b) \\over 2} + O(x^2)\n",
"$$\n",
"But $f(x)$ has a pole:\n",
"$$\n",
" {\\E^{\\I a x} - \\E^{\\I b x} \\over x^2} = {\\I (a-b) \\over x} + O(1)\n",
"$$\n",
"To rectify this, we need to be a bit more careful. First note that\n",
"$$\n",
"\\int_{\\infty}^\\infty {\\cos a x - \\cos b x \\over x^2} \\dx = \\lim_{\\epsilon \\rightarrow 0} \\left(\\int_{-\\infty}^{-\\epsilon} + \\int_\\epsilon^\\infty \\right){\\cos a x - \\cos b x \\over x^2} \\dx = \\Re \\dashint_{-\\infty}^\\infty f(x) \\dx\n",
"$$\n",
"Then we construct a contour avoiding zero as follows:"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M = 10\n",
"ε = 1.0\n",
"\n",
"plot(Segment(-M, -ε) ∪ Segment(ε, M);label=\"[-M,e] U [e,M]\")\n",
"plot!(Arc(0.,ε, (π,0.)); label=\"C_e\")\n",
"plot!(Arc(0., M, (0,π)); label = \"C_M\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that $\\oint_\\gamma f(z) \\dz = 0$, \n",
"$$\\int_{C_\\epsilon} {\\E^{\\I a z } - \\E^{\\I b z} \\over z^2} \\dz =\n",
" \\int_\\pi^0 { (b-a) \\E^{\\I \\theta} + O(\\epsilon) \\over \\E^{ \\I \\theta}} \\D\\theta \\rightarrow (a-b) \\pi\n",
"$$\n",
"Also, as the integrand is $O(z^{-2})$ the integral over $C_M$ vanishes as $M \\rightarrow \\infty$. We therefore get\n",
"$$\n",
" \\dashint f(x)\\dx = (b-a)\\pi\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4.703235784846408"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ε =0.001\n",
"M = 600.0\n",
"x = Fun(Segment(-M , -ε) ∪ Segment(ε, M))\n",
"a = 2.3; b = 3.8\n",
"sum((cos(a*x) - cos(b*x))/x^2) # Converges if we make M bigger"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4.71238898038469"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"π*(b-a)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 10.\n",
"Use binomial formula\n",
"\\begin{align*}\n",
"\\int_0^{2 \\pi} (\\cos \\theta)^n \\D \\theta &= {1 \\over 2^n \\I} \\oint (z+z^{-1})^n {\\dz \\over z} \\\\\n",
" &= {1 \\over 2^n \\I} \\sum_{k=0}^n {n! \\over k! (n-k)!} \\oint z^k z^{k-n} {\\dz \\over z} \\\\\n",
" &= {1 \\over 2^n \\I} \\sum_{k=0}^n {n! \\over k! (n-k)!} \\oint z^{2k-n-1}\\dz\n",
"\\end{align*}\n",
"\n",
"We only have a residue of $2k-n-1 = -1$, that is, if $2k = n$. If $n$ is odd, this can't happen (duh! the integral is symmetric with respect to $\\theta$). If it's even, then we have\n",
"$$\n",
"\\int_0^{2 \\pi} (\\cos \\theta)^n \\D \\theta = {\\pi \\over 2^{n-1} } {n! \\over 2 (n/2)!} \n",
" $$"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.356194490192353"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"θ = Fun(0 .. 2π)\n",
"n = 4;\n",
"sum(cos(θ)^n)"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.356194490192345"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"π*factorial(1.0n)/(2^(n-1)*2*factorial(n/2))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 2.1 \n",
"\n",
"By integrating around a rectangular contour with vertices at $\\pm R$ and $\\pi \\I \\pm R$ and letting $R \\rightarrow \\infty$, show that:\n",
"\t$$\\int_0^\\infty \\sech x \\dx = {\\pi \\over 2}$$\n",
"where $\\sech x = {2 \\over \\E^{-x} + \\E^x}$.\n",
"\n",
"\n",
"Recall $\\sech x = {2 \\over \\E^{-x} + \\E^x} $. This shows that $\\sech(-x) = \\sech x$\n",
"But we also have \n",
"$\n",
"\\sech(x + \\I \\pi) = {2 \\over \\E^{-x-\\I \\pi} + \\E^{x+ \\I \\pi}} = {2 \\over -\\E^{-x} - \\E^{x}} = -\\sech x\n",
"$\n",
"Thus we have\n",
"$$\n",
" 4\\int_0^\\infty \\sech x \\dx = \\left[ \\int_{-\\infty}^\\infty + \\int_{\\infty+\\I\\pi}^{-\\infty + \\I \\pi} \\right] \\sech z \\dz\n",
"$$\n",
"\n",
"We can approximate this using\n",
"$$\n",
"\\left[ \\int_{-R}^R + \\int_R^{R+\\I \\pi} + \\int_{R+\\I \\pi}^{-R+\\I\\pi} + \\int_{-R+\\I \\pi} \\right ] \\sech z \\dz = 2 \\pi \\I \\Res_{z = {\\I \\pi \\over 2}} \\sech z = 2 \\pi \n",
"$$\n",
"since, for $z_0 = {\\I \\pi \\over 2}$, we have\n",
"$$\n",
"\\sech z = {1 \\over \\cos \\I z} = {1 \\over - \\I \\sin \\I z_0 (z- z_0) + O(z-z_0)^2} \n",
"= -{\\I \\over (z- z_0)} + O(1)\n",
"$$\n",
"\n",
"Finally, we need to show that the limit as $R \\rightarrow \\infty$ tends to the right value. In this case, it follows since \n",
"$$\n",
" \\left|\\int_R^{R+\\I \\pi} \\sech z \\dz \\right| \\leq {2 \\pi \\E^{-R} \\over 1 - \\E^{-2R}} \\rightarrow 0\n",
"$$\n",
"(and by symmetry for $\\int_{-R+\\I \\pi}^{-R}$.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 2.2 \n",
"\n",
"Show that the Fourier transform of $\\sech x$ satisfies\n",
"\t$$\n",
"\t\\int_{-\\infty}^\\infty \\E^{\\I k x } \\sech x \\dx = \\pi \\sech {\\pi k \\over 2}\n",
"\t$$\n",
" \n",
"Define \n",
"$$f(z) = \\E^{\\I k z } \\sech z = {2 \\E^{(1+\\I k) z} \\over \\E^{2 z } + 1}$$\n",
"In this case, we have the symmetry\n",
"$$\n",
"f(x + \\I \\pi) = - \\E^{- k \\pi} \\E^{\\I k x} \\sech x = - \\E^{-k \\pi} f(x)\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.00032444937189257726 + 0.0003756543878221788im"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"k = 2.0\n",
"f = z -> exp(im*k*z)*sech(z)\n",
"\n",
"-exp(-k*π)f(2.0)"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0003244493718925772 + 0.00037565438782217884im"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f(2.0+im*π)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In other words, we have\n",
"$$\n",
"(1 + \\E^{-k \\pi}) \\int_{-\\infty}^\\infty f(x) \\dx =\\left(\\int_{-\\infty}^\\infty + \\int_{\\infty+\\I \\pi}^{-\\infty + \\I \\pi}\\right) f(z) \\dz \n",
"$$\n",
"By similar logic as above, we can show that the integral over the rectangular contour converges to this.\n",
"\n",
"Again, the only pole inside is at $z = {\\I \\pi \\over 2}$, where the residue is $-\\I \\E^{-\\pi k \\over 2}$. Thus we have\n",
"$$\n",
"\\int_{-\\infty}^\\infty f(x) \\dx = {2 \\pi \\E^{-\\pi k \\over 2} \\over 1 + \\E^{-k \\pi}} = \\pi \\sech{\\pi k \\over 2}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 3.1\n",
"\n",
"We have \n",
"$$\n",
"\\rho(A) \\subseteq B(1,3) \\cup B(2,3) \\cup B(4,1)\n",
"$$\n",
"where $B(z_0,r)$ is the ball of radius $r$ around $z_0$.\n",
"\n",
"Here's a depiction:"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"drawdisk!(z0, R) = plot!(θ-> real(z0) + R[1]*cos(θ), θ-> imag(z0) + R[1]*sin(θ), 0, 2π, fill=(0,:red), α = 0.2, legend=false)\n",
"\n",
"A = [1 2 -1; -2 2 1; 0 1 4]\n",
"\n",
"λ = eigvals(A)\n",
"p = plot()\n",
"drawdisk!(1,3)\n",
"drawdisk!(2,3)\n",
"drawdisk!(4,1)\n",
"scatter!(complex.(λ); label=\"eigenvalues\")\n",
"scatter!(complex.(diag(A)); label=\"diagonals\")\n",
"p"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 3.2 \n",
"\n",
"We get \n",
"$$\n",
"\\rho(A) \\subseteq B(1,2) \\cup B(2,3) \\cup B(4,2)\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"λ = eigvals(A)\n",
"p = plot()\n",
"drawdisk!(1,2)\n",
"drawdisk!(2,3)\n",
"drawdisk!(4,2)\n",
"scatter!(complex.(λ); label=\"eigenvalues\")\n",
"scatter!(complex.(diag(A)); label=\"diagonals\")\n",
"p"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 3.3\n",
"\n",
"Because the spectrum live in the intersection of the two estimates, the sharpest bound is\n",
"$$\n",
"\\rho(A) \\subseteq B(2,3) \n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"λ = eigvals(A)\n",
"p = plot()\n",
"drawdisk!(2,3)\n",
"scatter!(complex.(λ); label=\"eigenvalues\")\n",
"scatter!(complex.(diag(A)); label=\"diagonals\")\n",
"p"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Thus we can take $2 + 3 \\E^{\\I \\theta}$ as the contour."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 4.1\n",
"\n",
"Note that in the scalar case $u'' = a u$ we have the solution\n",
"$$\n",
"u(t) = u_0 \\cosh \\sqrt{a} t + v_0 {\\sinh \\sqrt{a} t \\over \\sqrt a}\n",
"$$\n",
"Write \n",
"$$A = Q \\begin{pmatrix} \\lambda_1 \\\\ & \\ddots \\\\\n",
" && \\lambda_n \\end{pmatrix} Q^\\top $$\n",
"where $\\lambda_k > 0$ \n",
"and then the solution has the form, where $\\gamma$ is a contour surrounding the eigenvalues and to the right of zero: \n",
"\\begin{align*}\n",
" \\vc u(t) = Q \\begin{pmatrix} \\cosh \\sqrt{\\lambda_1} t \\\\ & \\ddots \\\\\n",
" && \\cosh \\sqrt{\\lambda_n} t \\end{pmatrix} Q^\\top \\vc u_0 + \n",
" Q \\begin{pmatrix} {\\sinh \\sqrt{\\lambda_1} t \\over \\sqrt{\\lambda_1}}\\\\ & \\ddots \\\\\n",
" && {\\sinh \\sqrt{\\lambda_n} t \\over \\sqrt{\\lambda_n}} \\end{pmatrix} Q^\\top \\vc v_0 \\\\\n",
" = {1 \\over 2 \\pi \\I} Q \\begin{pmatrix} \\oint_\\gamma {\\cosh \\sqrt{z} t \\dz \\over z - \\lambda_1}\\\\ & \\ddots \\\\\n",
" && \\oint_\\gamma {\\cosh \\sqrt{z} t \\dz \\over z - \\lambda_n}\\end{pmatrix} Q^\\top \\vc u_0 + \n",
" {1 \\over 2 \\pi \\I} Q \\begin{pmatrix} \\oint_\\gamma {\\sinh \\sqrt{z} t \\over \\sqrt z} {\\dz \\over z - \\lambda_1} \\\\ & \\ddots \\\\\n",
" && \\oint_\\gamma {\\sinh \\sqrt{z} t \\over \\sqrt z} {\\dz \\over z - \\lambda_n} \\end{pmatrix} Q^\\top \\vc v_0 \\\\\n",
"= {1 \\over 2 \\pi \\I} \\oint_\\gamma \\cosh \\sqrt{z} t Q \\begin{pmatrix} (z - \\lambda_1)^{-1}\\\\ & \\ddots \\\\\n",
" && (z - \\lambda_n)^{-1}\\end{pmatrix} Q^\\top \\vc u_0 \\dz + {1 \\over 2 \\pi \\I} \\oint_\\gamma {\\sinh \\sqrt{z} t \\over \\sqrt z} Q \\begin{pmatrix} (z - \\lambda_1)^{-1}\\\\ & \\ddots \\\\\n",
" && (z - \\lambda_n)^{-1}\\end{pmatrix} Q^\\top \\vc v_0 \\dz \\\\\n",
"= {1 \\over 2 \\pi \\I} \\oint_\\gamma \\cosh \\sqrt{z} t (z I - A)^{-1} \\vc u_0 \\dz + {1 \\over 2 \\pi \\I} \\oint_\\gamma {\\sinh \\sqrt{z} t \\over \\sqrt z} (z I - A)^{-1} \\vc v_0 \\dz \n",
"\\end{align*}\n",
"\n",
"Here we verify the formulae numerically:"
]
},
{
"cell_type": "code",
"execution_count": 54,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3.323277889766459e-14"
]
},
"execution_count": 54,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n = 5\n",
"A = randn(n,n) \n",
"A = A+ A' + 10I\n",
"\n",
"λ, Q = eigen(A)\n",
"\n",
"norm(A - Q*Diagonal(λ)*Q')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Time-stepping solution:"
]
},
{
"cell_type": "code",
"execution_count": 57,
"metadata": {},
"outputs": [],
"source": [
"u₀ = randn(n)\n",
"v₀ = randn(n)\n",
"uv = solve(ODEProblem((uv,_,t) -> [uv[n+1:end]; A*uv[1:n]], [u₀; v₀], (0.,2.)); reltol=1E-10);"
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Array{Float64,1}:\n",
" 431.3487825666393\n",
" 842.6721648972994\n",
" -510.6730776173628\n",
" -1275.0672015889802\n",
" -852.8144701538115"
]
},
"execution_count": 58,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = 2.0\n",
"uv(t)[1:n]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Solution via diagonalization:"
]
},
{
"cell_type": "code",
"execution_count": 60,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Array{Float64,1}:\n",
" 431.34878253723247\n",
" 842.672164825094 \n",
" -510.67307756879427\n",
" -1275.0672015062562 \n",
" -852.8144700629356 "
]
},
"execution_count": 60,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Q*Diagonal(cosh.(sqrt.(λ) .* t))*Q'*u₀ + Q*Diagonal(sinh.(sqrt.(λ) .* t) ./ sqrt.(λ))*Q'*v₀"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Solution via elliptic integrals. We chose the ellipse to surround all the spectrum of our particular $A$ with eigenvalues:"
]
},
{
"cell_type": "code",
"execution_count": 63,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Array{Complex{Float64},1}:\n",
" 431.3489776293208 - 2.1947119676805203e-14im\n",
" 842.6728232388965 - 8.751061717472864e-15im \n",
" -510.67355346183 + 2.8015565557209588e-14im\n",
" -1275.067724375493 - 6.622136600059223e-14im \n",
" -852.8154291281783 + 2.365862961710633e-14im "
]
},
"execution_count": 63,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"periodic_rule(n) = 2π/n*(0:(n-1)), 2π/n*ones(n)\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\n",
"\n",
"function ellipse_f(f, 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]*f(z[j])*inv(z[j]*I - A)\n",
" end\n",
" ret/(2π*im)\n",
"end\n",
"\n",
"\n",
"n = 50\n",
"ellipse_f(z -> cosh(sqrt(z)*t), A, n, 10.0, 7.0, 2.0)*u₀ +\n",
" ellipse_f(z -> sinh(sqrt(z)*t)/sqrt(z), A, n, 10.0, 7.0, 2.0)*v₀"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## Problem 4.2\n",
"\n",
"I put the restriction in because of the $\\sqrt z$ term, which look like it is not analytic on $(-\\infty,0]$. However, this restriction was NOT necessary, since in fact $\\cosh \\sqrt z t$ and ${\\sinh \\sqrt z t \\over \\sqrt z}$ are entire:"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"phaseplot(-6..6, -3..3, z -> cosh(sqrt(z)))"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"phaseplot(-3..3, -3..3, z -> sinh(sqrt(z))/sqrt(z))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" This follows from Taylor series, though I prefer the following argument: we have\n",
"$$\n",
"\\cosh z = {\\E^{z} + \\E^{-z} \\over 2}\n",
"$$\n",
"hence $\\cosh{\\I t } = \\cos t = \\cosh{-\\I t}$. Therefore, on the possible branch cut $\\cosh \\sqrt z$ is continuous (hence analytic):\n",
"$$\n",
"\\cosh \\sqrt{x}_+ = \\cosh \\I \\sqrt{|x|} = \\cosh -\\I \\sqrt{|x|} = \\cosh\\sqrt{x}_-\n",
"$$\n",
"Similarly, \n",
"$$\n",
"\\sinh z = {\\E^z - \\E^{-z} \\over 2}\n",
"$$\n",
"implies $\\sinh{\\I t} = \\I \\sin t = - \\I \\sin(-t) = -\\sinh{-\\I t}$, which gives us continuity:\n",
"$$\n",
"{\\sinh \\sqrt{x}_+ \\over \\sqrt{x}_+} = {\\sinh \\I \\sqrt{|x|} \\over \\I \\sqrt{|x|}} = {\\sinh(-\\I \\sqrt{|x|)} \\over -\\I \\sqrt{|x|}} ={\\sinh \\sqrt{x}_- \\over \\sqrt{x}_-}\n",
"$$\n",
"furthermore, they are both bounded at zero, hence analytic there too.\n",
"\n",
"Here's a numerical example:"
]
},
{
"cell_type": "code",
"execution_count": 64,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5.93500670014108e-15"
]
},
"execution_count": 64,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n = 5\n",
"A = randn(n,n) \n",
"λ, V = eigen(A)\n",
"\n",
"norm(A - V*Diagonal(λ)*inv(V))"
]
},
{
"cell_type": "code",
"execution_count": 65,
"metadata": {},
"outputs": [],
"source": [
"u₀ = randn(n)\n",
"v₀ = randn(n)\n",
"uv = solve(ODEProblem((uv,_,t) -> [uv[n+1:end]; A*uv[1:n]], [u₀; v₀], (0.,2.)); reltol=1E-10);"
]
},
{
"cell_type": "code",
"execution_count": 66,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Array{Float64,1}:\n",
" 3.803344206258081 \n",
" 0.9252585031733016\n",
" -1.193805803087017 \n",
" 5.72798765024654 \n",
" 1.7188451287601578"
]
},
"execution_count": 66,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = 2.0\n",
"uv(t)[1:n]"
]
},
{
"cell_type": "code",
"execution_count": 67,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Array{Complex{Float64},1}:\n",
" 3.803344206066801 + 8.317464758586041e-16im \n",
" 0.9252585029740441 + 2.2092966791530562e-16im\n",
" -1.1938058030645906 + 9.308804664865466e-17im \n",
" 5.727987650140507 + 7.402447551422406e-16im \n",
" 1.7188451286631823 - 1.5268203762514324e-16im"
]
},
"execution_count": 67,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"V*Diagonal(cosh.(sqrt.(λ) .* t))*inv(V)*u₀ + \n",
" V*Diagonal(sinh.(sqrt.(λ) .* t) ./ sqrt.(λ))*inv(V)*v₀"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here's the solution using an elliptic integral:"
]
},
{
"cell_type": "code",
"execution_count": 68,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Array{Complex{Float64},1}:\n",
" 3.8033442060667997 - 1.0942835705045212e-16im\n",
" 0.9252585029740352 + 2.0727119042774893e-16im\n",
" -1.1938058030645888 - 2.8947259388740214e-16im\n",
" 5.7279876501405145 - 1.9746496661590837e-16im\n",
" 1.7188451286631805 + 8.925757831726994e-17im "
]
},
"execution_count": 68,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n = 100\n",
"ellipse_f(z -> cosh(sqrt(z)*t), A, n, 0.0, 3.0, 3.0)*u₀ +\n",
" ellipse_f(z -> sinh(sqrt(z)*t)/sqrt(z), A, n, 0.0, 3.0, 3.0)*v₀ "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 5.1\n",
"\n",
"\n",
"\n",
"We have \n",
"$$\n",
"{1 \\over n} \\sum_{j=0}^{n-1} g(\\theta_j) = \\sum_{k=0}^\\infty g_k {1 \\over n} \\sum_{j=0}^{n-1} \\E^{\\I k \\theta_j}\n",
"$$\n",
"Define the $n$-th root of unity as $\\omega = e^{2 \\pi \\I \\over n}$ (that is $\\omega^n = 1$), and simplify \n",
"$$\n",
"\\sum_{j=0}^{n-1} \\E^{\\I k \\theta_j} =\\sum_{j=0}^{n-1} e^{{2 \\pi j \\I k \\over n}} =\\sum_{j=0}^{n-1} \\omega^{kj}= \\sum_{j=0}^{n-1} (\\omega^k)^j\n",
"$$\n",
"If $k$ is a multiple of $n$, then $\\omega^k = 1$, and this sum is equal to $n$. If $k$ is not a multiple of $n$, use Geometric series:\n",
"$$\n",
"\\sum_{j=0}^{n-1} (\\omega^k)^j = {\\omega^{nk} - 1 \\over \\omega^k -1} = {1^k - 1 \\over \\omega^k-1} = 0\n",
"$$\n",
"## Problem 5.2 \n",
"\n",
"\n",
"From lecture 4, we have $|f_k| \\leq M_r r^{-|k|}$ for any $1 \\leq r < R$, where $M_r$ is the supremum of $f$ in an annulus $\\{ z : r^{-1} < |z| < r \\}$. Thus from the previous part we have (using geometric series)\n",
"$$\n",
"\\left| {1 \\over n} \\sum_{j=0}^{n-1} g(\\theta_j) - { 1\\over 2 \\pi}\t \\int_0^{2\\pi} g(\\theta) \\D \\theta \\right| \\leq \\sum_{K=1}^\\infty |f_{Kn}| + |f_{-Kn}| \\leq 2M \\sum_{K=1}^\\infty r^{-Kn} = 2M_r {r^{-n} \\over 1 - r^{-n}}.\n",
"$$\n",
"This is an upper bound that decays exponentially fast. \n",
"\n",
"## Problem 5.3\n",
"\n",
"Note that $f(z) = 2 z/(4z - z^2 - 1)$ satisfies $f(\\E^{\\I \\theta}) = g(\\theta)$. This has two poles at $2 \\pm \\sqrt 3$:"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.1102230246251565e-16 + 0.0im"
]
},
"execution_count": 43,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f = z -> 2z/(4z-z^2-1)\n",
"f(exp(0.1im)) - 1/(2-cos(0.1))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here's a phase plot showing the location of poles:"
]
},
{
"cell_type": "code",
"execution_count": 69,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 69,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"phaseplot(-5..5, -5..5, f)"
]
},
{
"cell_type": "code",
"execution_count": 72,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(3.732050807568877, 0.2679491924311228, 0.2679491924311227)"
]
},
"execution_count": 72,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"2+sqrt(3), 2- sqrt(3), 1/(2+sqrt(3))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that $2+\\sqrt 3 = 1/(2-\\sqrt 3)$ so in the previous result, we take $R = 2+\\sqrt 3$. For any $1 \\leq r < 2+\\sqrt 3$ we have \n",
"$$\n",
"M_r = {2 \\over 4 - r - r^{-1}}\n",
"$$\n",
"\n",
"Thus we get the upper bounds\n",
"$$\n",
" {4 \\over 4 - r - r^{-1}} {r^{-n} \\over 1 - r^{-n}}\n",
"$$\n",
"\n",
"Let's see how sharp it is:"
]
},
{
"cell_type": "code",
"execution_count": 73,
"metadata": {},
"outputs": [],
"source": [
"periodic_rule(n) = 2π/n*(0:(n-1)), 2π/n*ones(n)\n",
"\n",
"g = θ -> 1/(2 - cos(θ))\n",
"\n",
"Q = sum(Fun(g, 0 .. 2π))\n",
"\n",
"\n",
"err = Float64[ (\n",
" (θ, w) = periodic_rule(n);\n",
" sum(w.*g.(θ)) - Q\n",
" ) for n=1:30];"
]
},
{
"cell_type": "code",
"execution_count": 78,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"