{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true,
"scrolled": false
},
"outputs": [],
"source": [
"using ApproxFun, Plots, ComplexPhasePortrait, ApproxFun, SingularIntegralEquations,\n",
" SpecialFunctions, DualNumbers, SO\n",
"\n",
"using SingularIntegralEquations.HypergeometricFunctions\n",
"gr();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# M3M6: Methods of Mathematical Physics\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\\ds{\\D s}\n",
"\\def\\C{{\\mathbb C}}\n",
"\\def\\R{{\\mathbb R}}\n",
"\\def\\H{{\\mathbb H}}\n",
"\\def\\CC{{\\cal C}}\n",
"\\def\\HH{{\\cal H}}\n",
"\\def\\FF{{\\cal F}}\n",
"\\def\\I{{\\rm i}}\n",
"\\def\\Ei{{\\rm Ei}\\,}\n",
"\\def\\qqqquad{\\qquad\\qquad}\n",
"\\def\\qqand{\\qquad\\hbox{and}\\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\\acos{\\,{\\rm acos}\\,}\n",
"\\def\\erfc{\\,{\\rm erfc}\\,}\n",
"\\def\\vc#1{{\\mathbf #1}}\n",
"\\def\\ip<#1,#2>{\\left\\langle#1,#2\\right\\rangle}\n",
"\\def\\br[#1]{\\left[#1\\right]}\n",
"\\def\\norm#1{\\left\\|#1\\right\\|}\n",
"\\def\\half{{1 \\over 2}}\n",
"\\def\\fL{f_{\\rm L}}\n",
"\\def\\fR{f_{\\rm R}}\n",
"\\def\\HF(#1,#2,#3,#4){{}_2F_1\\!\\left({#1,#2 \\atop #3}; #4 \\right)}\n",
"\\def\\questionequals{= \\!\\!\\!\\!\\!\\!{\\scriptstyle ? \\atop }\\,\\,\\,}\n",
"$$\n",
"\n",
"Dr Sheehan Olver\n",
"
\n",
"s.olver@imperial.ac.uk\n",
"\n",
"Office Hours: 3-4pm Mondays, 11-12am Thursdays, Huxley 6M40\n",
"
\n",
"Website: https://github.com/dlfivefifty/M3M6LectureNotes\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"# Lecture 27: Singular points of ordinary differential equations\n",
"\n",
"\n",
"2. Singular points of ODEs\n",
" - regular singular point\n",
" - irregular singular point\n",
" - singular points at ∞"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Singular points\n",
"\n",
"We will consider the first order ODE\n",
"$$\n",
"u' + a(z) u = 0\n",
"$$\n",
"and the second order ODE\n",
"$$\n",
"u'' + a(z) u' + b(z) u = 0\n",
"$$\n",
"where $a$ and $b$ are rational. \n",
"\n",
"The simplest point is an ordinary point, where we know that $u$ is analytic:\n",
"\n",
"**Definition (ordinary point)** $z_0$ is a ordinary point of the ODE if $a(z)$ and $b(z)$ are analytic at $z_0$.\n",
"\n",
"\n",
"### Regular singular point\n",
"\n",
"Suppose $a$ has a simple pole, and w.l.o.g. take $z_0=0$, that is: \n",
"$$\n",
"a(z) = {a_{-1} \\over z} + a_0 + a_1 z + a_2 z^2 + \\cdots \n",
"$$\n",
"where $a_0 + a_1 z + a_2 z^2 + \\cdots $ is analytic.\n",
"Take as an ansatz to the first order ODE\n",
"$$u(z) = z^\\alpha ( 1 + z v(z)) = z^\\alpha + z^{\\alpha+1} v$$\n",
"where $\\alpha = a_{-1}$. Note that $c u(z)$ is a solution of $u'(z) = a(z) u(z)$: the choice of $1$ above makes the solution unique. \n",
"\n",
"\n",
"Note that\n",
"$$\n",
"u'(z) = \\alpha z^{\\alpha-1} + z^\\alpha ((\\alpha+1) v + z v')\n",
"$$\n",
"The ODE $u'(z) = a(z) u(z)$ becomes an inhomogenous equation:\n",
"\\begin{align*}\n",
"z^\\alpha ((\\alpha+1) v + z v' - z a(z) v) &= z^{\\alpha} (a(z) - {\\alpha \\over z} ) & \\Leftrightarrow\\\\\n",
"z v'(z) + (\\alpha+1 - z a(z)) v(z) &= a_0 + a_1 z + a_2 z^2 + \\cdots\n",
"\\end{align*}\n",
"In operator form, expressing \n",
"$$\n",
"v(z) = v_0 + v_1z + v_2 z^2 = (1,z,z^2,\\ldots)\\begin{pmatrix}v_0\\\\v_1\\\\\\vdots \\end{pmatrix}\n",
"$$\n",
"we have\n",
"\\begin{align*}\n",
"z v' &=(1,z,z^2,\\ldots) \\begin{pmatrix} 0 \\cr & 1 \\cr && 2 \\cr &&& 3 \\cr &&&&\\ddots \\end{pmatrix} \\begin{pmatrix}v_0\\\\v_1\\\\\\vdots \\end{pmatrix} \\\\\n",
"(\\alpha+1 - z a(z)) v(z) &= (1,z,z^2,\\ldots) \\begin{pmatrix} 1 \\cr -a_0 & 1 \\cr -a_1 & -a_0& 1 \\cr -a_2 & -a_1&-a_0& 1 \\cr \\vdots&\\ddots&\\ddots&\\ddots&\\ddots \\end{pmatrix} \\begin{pmatrix}v_0\\\\v_1\\\\\\vdots \\end{pmatrix}\n",
"\\end{align*}\n",
"Thus the equation for our unknown coefficients becomes:\n",
"$$\n",
"\\begin{pmatrix} 1 \\cr -a_0 & 2 \\cr -a_1 & -a_0& 3 \\cr -a_2 & -a_1&-a_0& 4 \\cr \\vdots&\\ddots&\\ddots&\\ddots&\\ddots \\end{pmatrix} \\begin{pmatrix}v_0\\\\v_1\\\\\\vdots \\end{pmatrix} = \\begin{pmatrix}a_0\\\\a_1\\\\\\vdots \\end{pmatrix}\n",
"$$\n",
"This equation induces decay in $|v_k|$ just like last lecture, hence we know its analytic with the same radius of convergence as $a$.\n",
"\n",
"Here is an example: as we can see, `v(z)` is analytic in a disk of radius `r` where `r` is dictated by the function `a`. "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"r = 0.9 # radius we are solving in, can be up to but not including 1\n",
"S = Taylor(Circle(r))\n",
"α = 0.3\n",
"a = z -> α/z + 1/(1-z)\n",
"ã = Fun( z -> a(z) - α/z, S)\n",
"za = Fun( z -> z*a(z), S)\n",
"z = Fun(S)\n",
"D = Derivative() : S → S\n",
"L = z*D + (α+1-za)\n",
"v = L \\ ã\n",
"\n",
"# v is analytic in disk of radius r\n",
"phaseplot(v, (-1,1), (-1,1))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here is a plot of `u`, which shows that it has a branch cut:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# u has a branch \n",
"u = z -> z^α + z^(α+1)*v(z)\n",
"phaseplot(u, (-1,1), (-1,1))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can confirm it solves the ODE:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4.884981308350689e-15 - 7.521914234669674e-17im"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"up = dualpart(u(dual(0.1,1.0))) # uses DualNumbers.jl to calculate derivative\n",
"\n",
"up - a(0.1)u(0.1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Definition (first order regular singular point)** $z_0$ is a regular singular point of the first order ODE if $(z-z_0)a(z)$ is analytic at $z_0$.\n",
"\n",
"For second order equations, we will arrive at the following analogue:\n",
"\n",
"**Definition (second order regular singular point)** $z_0$ is a regular singular point of the second order ODE if $(z-z_0) a(z)$ and $(z-z_0)^2b(z)$ are analytic at $z_0$.\n",
"\n",
"Assume \n",
"$$\n",
"a(z) = {a_{-1} \\over z} + a_0 + a_1 z + \\cdots\n",
"$$\n",
"and\n",
"$$\n",
"b(z) = {b_{-2} \\over z^2} + {b_{-1} \\over z} + b_0 + b_1 z + \\cdots\n",
"$$\n",
"To determine the right ansatz, consider\n",
"$$\n",
"u'' + {a_{-1} \\over z} u' + {b_{-2} \\over z^2} u'' = 0\n",
"$$\n",
"which has solutions $z^\\alpha$ solving the indicial equation:\n",
"\n",
"**Definition (indicial equation)**\n",
"$$\n",
"\\alpha(\\alpha-1) + a_{-1} \\alpha + b_{-2} = 0\n",
"$$\n",
"\n",
"Now take as an ansatz:\n",
"$$\n",
" u(z) = z^\\alpha (1 + z v(z)) = z^\\alpha + z^{\\alpha+1} v(z)\n",
"$$\n",
"Differentiating we have\n",
"\\begin{align*}\n",
"u'(z) &= \\alpha z^{\\alpha-1} + z^\\alpha((\\alpha+1) v + z v') \\\\\n",
"u''(z) &= \\alpha(\\alpha-1) z^{\\alpha-2} + z^{\\alpha-1}( \\alpha(\\alpha+1) v + \\alpha z v' + (\\alpha+1) z v' + z v' + z^2 v'' ) \\\\\n",
" &= \\alpha(\\alpha-1) z^{\\alpha-2} + z^{\\alpha-1}( \\alpha(\\alpha+1) v + 2(\\alpha+1) z v' + z^2 v'' ) \n",
"\\end{align*}\n",
"Thus $v$ solves \n",
"$$\n",
"z^2 v''+ z (2(\\alpha+1) + z a(z)) v' + ( \\alpha(\\alpha+1) + z a(z) (\\alpha+1) + z^2 b(z)) v = - {\\alpha(\\alpha-1) \\over z} - a(z) \\alpha -z b(z) = \\underbrace{- \\alpha a_0 - b_{-1} - (\\alpha a_1 + b_0) z - (\\alpha a_2 + b_1) z^2 + \\cdots}_{f(z) = f_0 + f_1 z + \\cdots}\n",
"$$\n",
"\n",
"The question is when can forward substitution proceed. Let's write out the operators, first we have:\n",
"$$\n",
"z^2 v'' \\equiv \\begin{pmatrix} 0 \\cr & 0 \\cr &&2 \\cr &&& 6 \\cr &&&&12 \\cr &&&&&\\ddots \\end{pmatrix}\n",
"$$\n",
"Then we have\n",
"$$\n",
"2 (\\alpha+1) z v' \\equiv \\begin{pmatrix} 0\\cr & 2 (\\alpha+1) \\cr &&&4 (\\alpha+1) \\cr &&&& 6(\\alpha+1) \\cr &&&&&8 (\\alpha+1) \\cr &&&&&&\\ddots \\end{pmatrix}\n",
"$$\n",
"and\n",
"$$\n",
"z^2 a(z) v' = (z a(z))(z v') \\equiv \\begin{pmatrix} a_{-1} \\\\ a_0 & a_{-1} \\\\ a_1 & a_0 & a_{-1} \\\\ a_2 & a_1 & a_0 & a_{-1} \\\\ \\vdots &\\ddots&\\ddots&\\ddots&\\ddots \\end{pmatrix} \\begin{pmatrix} 0 \\cr & 1 \\cr &&2 \\cr &&& 3 \\cr &&&&4 \\cr &&&&&\\ddots \\end{pmatrix}\n",
"$$\n",
"Finally since\n",
"\\begin{align*}\n",
"\\alpha(\\alpha+1) + z (\\alpha+1) a(z) + z^2 b(z) &= \\alpha(\\alpha-1) + \\alpha a_{-1} + b_{-2} + 2 \\alpha + a_{-1} + ((\\alpha+1) a_0 + b_{-1}) z + ((\\alpha+2) a_1 + b_0) z^2 + \\cdots \\\\\n",
"& = 2 \\alpha + a_{-1} + ((\\alpha+1) a_0 + b_{-1}) z + ((\\alpha+2) a_1 + b_0) z^2 + \\cdots\n",
"\\end{align*}\n",
"we have\n",
"$$\n",
"(\\alpha(\\alpha+1) + z (\\alpha+1) a(z) + z^2 b(z) ) v \\equiv \n",
"\\alpha \\begin{pmatrix} \n",
"2 \\alpha + a_{-1} \\\\ \n",
"(\\alpha+1) a_0 + b_{-1} & 2 \\alpha + a_{-1} \\\\ \n",
"(\\alpha+1) a_1 + b_0 & (\\alpha+1) a_0 + b_{-1} & 2 \\alpha + a_{-1} \\\\ \n",
"(\\alpha+1) a_2 + b_1 & (\\alpha+1) a_1 + b_0 & (\\alpha+1) a_0 + b_{-1} & 2 \\alpha + a_{-1} \\\\ \n",
"\\vdots &\\ddots&\\ddots&\\ddots&\\ddots \\end{pmatrix}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We need to only worry about the diagonal to determine if forward substitution proceeds, in this case we have\n",
"$$\n",
"\\begin{pmatrix}\n",
"2 \\alpha + a_{-1} \\\\\n",
"\\times & 4 \\alpha + 2 + 2a_{-1} \\\\\n",
"\\times & \\times &6 \\alpha + 3a_{-1} + 6 \\\\\n",
"\\times & \\times & \\times & 8 \\alpha + 4a_{-1} + 12 \\\\\n",
"\\times & \\times & \\times & \\times & 10 \\alpha + 5a_{-1} + 20\\\\\n",
"\\vdots & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots\n",
"\\end{pmatrix}\n",
"$$\n",
"Thus the first condition for this to give us a solution is that $2 (k+1)\\alpha + (k+1)a_{-1} + k(k+1) \\neq 0$ for all $k=1,2,\\ldots$.\n",
"This can be thought of another way: note that plugging $\\alpha+k+1$ in to the indicial equation gives us:\n",
"$$\n",
"(\\alpha+k+1) (\\alpha+k) + a_{-1} (\\alpha+k+1) + b_{-2} = (\\alpha) (\\alpha-1) + a_{-1} \\alpha + b_{-2}+ 2 (k+1)\\alpha + (k+1)a_{-1} + k(k+1) = 2 (k+1)\\alpha + (k+1)a_{-1} + k(k+1) \n",
"$$\n",
"Thus forward substitution fails if and only $\\alpha+k+1$ is the second root to the indicial equation for some $k=0,1,2,\\ldots$.\n",
"\n",
"**Remark** If this condition is satisfied, then $v$ is analytic near the singular point, which follows by bounding the forward elimination as before, but we'll skip this for now.\n",
"\n",
"**Remark** If $\\alpha+k+1$ is an integer, then using the ansatz $z^{\\alpha+k+1}(1 + z v(z))$ proceeds without problem. To get the second solution, we need to introduce logarithms, but we'll skip this for now."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We now check the results of solving\n",
"$$\n",
"z^2 v''+ z (2(\\alpha+1) + z a(z)) v' + ( \\alpha(\\alpha+1) + z a(z) (\\alpha+1) + z^2 b(z)) v = - {\\alpha(\\alpha-1) \\over z} - a(z) \\alpha -z b(z)\n",
"$$\n",
"numerically.\n",
"\n",
"We first choose a solutino to the an indicial equation:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"α = 1.3830951894845303\n"
]
},
{
"data": {
"text/plain": [
"2.7755575615628914e-16"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"S = Taylor()\n",
"z = Fun(S)\n",
"D = Derivative() : S\n",
"\n",
"a₋₁ = -0.6\n",
"b₋₂ = 0.3\n",
"\n",
"a = z -> a₋₁/z + 1/(2-z)\n",
"ã = Fun( z -> a(z) - a₋₁/z, S)\n",
"za = Fun( z -> z*a(z), S)\n",
"\n",
"b = z -> b₋₂/z^2 + 0.4/z + 1/(3-z)\n",
"zb̃ = Fun(z -> z*b(z) - b₋₂/z, S)\n",
"z²b = Fun( z -> z^2*b(z), S)\n",
"\n",
"α = ((1-a₋₁) + sqrt((a₋₁-1)^2-4b₋₂))/2\n",
"@show α\n",
"α*(α-1) + α*a₋₁ + b₋₂"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can solve for `v`, which is analytic in a neighbourhood of 0:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"L = z^2*D^2 + z*(2(α+1) + za)*D + (α*(α+1) + za*(α+1) + z²b)\n",
"v = L \\ (-α*ã - zb̃)\n",
"\n",
"phaseplot(v, (-3,3), (-3,3))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Thus one solution is:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"u = z -> z^α + z^(α+1)*v(z)\n",
"\n",
"phaseplot(u, (-3,3), (-3,3))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's double check it satisfies the original ODE:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3.3654295491028256e-12"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# dual numbers don't work to second degree so let's check numerically: \n",
"uf = Fun(u, 0.1 .. 0.2)\n",
"norm(uf'' + Fun(a, space(uf))*uf' + Fun(b, space(uf))*uf)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The second, linearly independent solution can be found taking the other $\\alpha$:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"α̃ = 0.21690481051546984\n"
]
},
{
"data": {
"text/plain": [
"1.6653345369377348e-16"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"α̃ = ((1-a₋₁) - sqrt((a₋₁-1)^2-4b₋₂))/2\n",
"@show α̃\n",
"α̃*(α̃-1) + α̃*a₋₁ + b₋₂"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"L = z^2*D^2 + z*(2(α̃+1) + za)*D + (α̃*(α̃+1) + za*(α̃+1) + z²b)\n",
"ṽ = L \\ (-α̃*ã - zb̃)\n",
"\n",
"phaseplot(ṽ, (-3,3), (-3,3))"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ũ = z -> z^α + z^(α+1)*ṽ(z)\n",
"\n",
"phaseplot(ũ, (-3,3), (-3,3))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Any general solution to the ODE is a linear combination of these two."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"### Irregular singular points\n",
"\n",
"Every isolated singularity that is not regular is an _irregular singular point_: even something as simple as ${\\D u / \\dz} = {u \\over z^2}$ has an irregular singular point. Irregular singular points correspond to essential singularities: in this case the solution is\n",
"$$\n",
"u(z) = \\E^{-1/z}\n",
"$$\n",
"\n",
"### Singular points at infinity\n",
"\n",
"To determine the type of singularity at $\\infty$ we do our usual change of variables $u(z) = v(1/z)$. For example, if\n",
"$$\n",
"{\\D u \\over \\dz} = a(z) u\n",
"$$\n",
"Then \n",
"$$\n",
"{\\D v \\over \\dz} = -z^{-2}a(z^{-1}) v\n",
"$$\n",
"**Definition (first order ordinary point at ∞)** $\\infty$ is an ordinary point of a first order ODE if $z^2a(z)$ is analytic at $\\infty$. \n",
"\n",
"**Definition (first order regular singular point at ∞)** $\\infty$ is a regular singular point of a first order ODE if $za(z) $ is analytic at $\\infty$. \n",
"\n",
"For second order, we see that $v$ satisfies\n",
"\\begin{align*}\n",
"v'(z) &= -z^{-2} u'(z^{-1}) \\\\\n",
"v''(z) &=2 z^{-3} u'(z^{-1}) + z^{-4} u''(z^{-1}) =-2 z^{-1} v'(z)+ z^{-4} u''(z^{-1})\n",
"\\end{align*}\n",
"which gives us the ODE\n",
"$$\n",
"u'' + a(z) u' + b(z) u \\Rightarrow\n",
"z^4v'' + (2 z^3 - a z^2) v' + b v\n",
"$$\n",
"\n",
"\n",
"**Definition (second order regular singular point at ∞)** $\\infty$ is a regular singular point of a first order ODE if $za(z)$ and $z^2 b(z)$ are analytic at $\\infty$. \n",
"\n",
"\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.6.0",
"language": "julia",
"name": "julia-0.6"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "0.6.1"
}
},
"nbformat": 4,
"nbformat_minor": 2
}