{
"cells": [
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"using Plots, ApproxFun, SingularIntegralEquations, DifferentialEquations, ComplexPhasePortrait\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",
"\\def\\pr(#1){\\left({#1}\\right)}\n",
"\\def\\br[#1]{\\left[{#1}\\right]}\n",
"\\def\\set#1{\\left\\{{#1}\\right\\}}\n",
"\\def\\ip<#1>{\\left\\langle{#1}\\right\\rangle}\n",
"\\def\\iip<#1>{\\left\\langle\\!\\langle{#1}\\right\\rangle\\!\\rangle}\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 2\n",
"\n",
"\n",
"## Problem 1.1\n",
"\n",
"\n",
"### 1.\n",
"\n",
"Take as an initial guess\n",
"$$\n",
" \\phi_1(z) = {\\sqrt{z-1} \\sqrt{z+1} \\over 2\\I(1 + z^2)}\n",
"$$\n",
"This satisfies for $-1 < x < 1$\n",
"$$\n",
"\\phi_1^+(x) -\\phi_1^-(x) = { \\sqrt{1-x^2} \\over 2(1 + x^2)} - {- \\sqrt{1-x^2} \\over 2(1 + x^2)} = {\\sqrt{1 - x^2} \\over 1+x^2}\n",
"$$\n",
"\n",
"Further, as $z \\rightarrow \\infty$, \n",
"$$\n",
"\\phi_1(z) \\sim {z \\over \\I(1+ z^2)} \\rightarrow 0\n",
"$$\n",
"\n",
"The catch is that it has poles at $\\pm \\I$:\n",
"\\begin{align*}\n",
" \\phi_1(z) = -{\\sqrt{\\I -1} \\sqrt{\\I+1} \\over 4} { 1 \\over z - \\I} + O(1) \\\\\n",
" \\phi_1(z) = {\\sqrt{-\\I -1} \\sqrt{-\\I+1} \\over 4} { 1 \\over z + \\I} + O(1) \n",
"\\end{align*}\n",
"Thus it follows that\n",
"$$\n",
"\\phi(z) = \\phi_1(z) + {\\sqrt{\\I -1} \\sqrt{\\I+1} \\over 4} { 1 \\over z - \\I} - {\\sqrt{-\\I -1} \\sqrt{-\\I+1} \\over 4} { 1 \\over z + \\I}\n",
"$$\n",
"is \n",
" 1. Analytic at $\\pm \\I$ and off $[-1,1]$ (Analyticity)\n",
" 2. $\\phi(\\infty) = 0$ (Decay)\n",
" 3. Has weaker than pole singularities (Regularity)\n",
" 4. Satisfies\n",
"$$\n",
"\\phi_+(x) - \\phi_-(x) = {\\sqrt{1-x^2} \\over 1+x^2}\n",
"$$\n",
" \n",
"By Plemelj II, this must be the Cauchy transform.\n",
"\n",
"_Demonstration_ We will see experimentally that it correct. First we do a phase plot to make sure we satisfy (Analyticity):"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"φ = z -> sqrt(z-1)sqrt(z+1)/(2im*(1+z^2)) + \n",
" sqrt(im-1)sqrt(im+1)/4*1/(z-im) - \n",
" sqrt(-im-1)sqrt(-im+1)/4*1/(z+im)\n",
"\n",
"phaseplot(-3..3, -3..3, φ)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can also see from the phase plot (Regularity): we have weaker than pole singularities, otherwise we would have at least a full,counter clockwise colour wheel. We can check decay as well:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0005177682933717976 + 0.000517765612545622im"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"φ(200.0+200.0im)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we compare it numerically it to `cauchy(f, z)` which is implemented in SingularIntegralEquations.jl:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.05303535516221752 + 0.05036581190871381im"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"φ(2.0+2.0im)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0530353551622175 + 0.050365811908713795im"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = Fun()\n",
"cauchy(sqrt(1-x^2)/(1+x^2), 2.0+2.0im)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 2.\n",
"\n",
"Recall that\n",
"$$\n",
"\\psi(z) = {\\log(z-1) - \\log(z+1) \\over 2 \\pi \\I}\n",
"$$\n",
"satisfies\n",
"$$\n",
"\\psi_+(x) - \\psi_-(x) = 1\n",
"$$\n",
"Therefore, consider\n",
"$$\n",
"\\phi_1(z) = {\\psi(z) \\over 2 + z}\n",
"$$\n",
"This has the right jump, but has an extra pole at $z = -2$: for $x < -1$ we have\n",
"$$\n",
"\\phi_1(x) = {\\log_+(x-1) - \\log_+(x+1) \\over 2 \\pi \\I} {1 \\over 2 + x} = \n",
" {\\log(1-x) - \\log(-1-x) \\over 2 \\pi \\I} {1 \\over 2 + x}\n",
"$$\n",
"hence we arrive at the solution\n",
"$$\n",
"\\phi_1(z) - {\\log 3 \\over 2 \\pi \\I (2+z)}\n",
"$$\n",
"We can verify that $\\phi_1(\\infty) = 0$."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"φ = z -> (log(z-1)-log(z+1)) / ((2π*im)*(2+z)) - log(3)/(2π*im*(2+z))\n",
"\n",
"phaseplot(-3..3, -3..3, φ)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0416136101650096 + 0.04191488519537722im"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"φ(2.0+2.0im)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.041613610165009585 + 0.04191488519537721im"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cauchy(1/(2+x), 2.0+2.0im)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 4. \n",
"\n",
"We first calculate the Cauchy transform of $f(x) = x/\\sqrt{1-x^2}$:\n",
"$$\n",
"\\phi(z) = {\\I z \\over 2 \\sqrt{z -1} \\sqrt{z+1}} - {\\I \\over 2}\n",
"$$\n",
"This vanishes at $\\infty$ and has the correct jump. We then have\n",
"$$\n",
"-\\I{\\cal H}f(x) = \\phi^+(x) + \\phi^-(x) = -\\I\n",
"$$\n",
"This implies that \n",
"$$\n",
"\\dashint_{-1}^1 {t \\over (t-x) \\sqrt{1-t^2}} \\dt = \\pi \\HH f(x) = \\pi\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3.141592653589793"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f = x/sqrt(1-x^2)\n",
"π*hilbert(f, 0.1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 1.2\n",
"\n",
"### 1.2.1\n",
"\n",
"From Problem 1.1 part 4, we have a solution: \n",
"$$\n",
"\\phi(z) = -{ z \\over 2 \\sqrt{z -1} \\sqrt{z+1}} + {1 \\over 2}\n",
"$$\n",
"All other solutions are then of the form:\n",
"$$\n",
"\\phi(z) + {C \\over \\sqrt{z-1} \\sqrt{z+1}}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"#11 (generic function with 1 method)"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"C = randn()\n",
"φ = z -> -z/(2*sqrt(z-1)*sqrt(z+1))+1/2 + C/(sqrt(z-1)*sqrt(z+1))"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.0 + 0.0im"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"φ(0.1+0.0im)+φ(0.1-0.0im)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.1249130872220166e-8"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"φ(1E8)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 1.2.2\n",
"\n",
"$\\psi(z) = \\phi(z) - 1$ satifies \n",
"$$\n",
"\\psi_+(x) + \\psi_-(x) = -2, \\psi(\\infty) = 0\n",
"$$\n",
"hence we know that $\\psi(z) = { z \\over \\sqrt{z -1} \\sqrt{z+1}} - 1 + {C \\over \\sqrt{z-1} \\sqrt{z+1}}$, giving \n",
"$$\n",
"\\psi(z) = { z \\over \\sqrt{z -1} \\sqrt{z+1}} + {C \\over \\sqrt{z-1} \\sqrt{z+1}} + 1\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"#13 (generic function with 1 method)"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"C = randn()\n",
"φ = z -> z/(sqrt(z-1)*sqrt(z+1)) + C/(sqrt(z-1)*sqrt(z+1)) "
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0 + 0.0im"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"φ(0.1+0.0im)+φ(0.1-0.0im)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.9999999983329383"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"φ(1E9)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 1.2.3 \n",
"\n",
"For $f(x) = \\sqrt{1-\\diamond^2}$, we use the formula\n",
"$$\n",
"\\phi(z) = {\\I \\over \\sqrt{z - 1} \\sqrt{z+1}} \\CC[\\sqrt{1-\\diamond^2} f](z) + {C \\over \\sqrt{z-1}\\sqrt{z+1}} = {\\I \\over \\sqrt{z - 1} \\sqrt{z+1}} \\CC[1-\\diamond^2](z) + {C \\over \\sqrt{z-1}\\sqrt{z+1}}\n",
"$$\n",
"We already know $\\CC1(z)$, and we can deduce $\\CC[\\diamond^2]$ as follows: try\n",
"$$\n",
"\\phi_1(z) = z^2 \\CC1(z) = z^2 {\\log(z-1) - \\log(z+1) \\over 2 \\pi \\I}\n",
"$$\n",
"this has the right jump, but blows up at $\\infty$ like:\n",
"$$\n",
"x^2 (\\log(x-1) - \\log(x+1)) = x^2 (\\log(1-1/x) - \\log(1+1/x) )\n",
"= -2 x + O(x^{-1})\n",
"$$\n",
"using\n",
"$$\n",
"\\log z = (z-1) - {1\\over 2}(z-1)^2 + O(z-1)^3\n",
"$$\n",
"Thus we have\n",
"$$\n",
"\\CC[\\diamond^2](z) = {z^2 (\\log(z-1) - \\log(z+1)) + 2 z \\over 2 \\pi \\I}\n",
"$$\n",
"and\n",
"$$\n",
"\\phi(z) = {\\I \\over \\sqrt{z - 1} \\sqrt{z+1}} {(1-z^2)(\\log(z-1) - \\log(z+1)) - 2 z \\over 2 \\pi \\I} + {C \\over \\sqrt{z-1}\\sqrt{z+1}}\n",
"$$\n",
"\n",
"_Demonstration_ Here we see that the Cauchy transform of $x^2$ has the correct formula:"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.028322293739596095 + 0.024377589786690298im"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"z = 2.0+2.0im\n",
"cauchy(x^2, z)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.028322293739596032 + 0.02437758978669024im"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(z^2*(log(z-1)-log(z+1))+2z)/(2π*im)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We now see that $\\phi$ has the right jumps:"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"#15 (generic function with 1 method)"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"C = randn()\n",
"φ = z -> im/(sqrt(z-1)*sqrt(z+1)) * ((1-z^2)*(log(z-1)-log(z+1))-2z)/(2π*im) + C/(sqrt(z-1)sqrt(z+1))"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.1102230246251565e-16 + 0.0im"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"φ(0.1+0.0im) + φ(0.1-0.0im) - sqrt(1-0.1^2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, it vanishes at infinity:"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-1.0949779687205364e-6 - 0.0im"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"φ(1E5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 1.2.4\n",
"\n",
"Let $f(x) = {1 \\over 1+x^2}$. From Problem 1.1 part 1 we know \n",
"$$\n",
"\\CC[\\sqrt{1-\\diamond^2}] f(z) = {\\sqrt{z-1} \\sqrt{z+1} \\over 2\\I(1 + z^2)} + {\\sqrt{\\I -1} \\sqrt{\\I+1} \\over 4} { 1 \\over z - \\I} - {\\sqrt{-\\I -1} \\sqrt{-\\I+1} \\over 4} { 1 \\over z + \\I}\n",
"$$\n",
"hence from the solution formula we have\n",
"$$\n",
"\\phi(z) = {1 \\over 2(1 + z^2)} + {\\sqrt{\\I -1} \\sqrt{\\I+1} \\I \\over 4\\sqrt{z-1} \\sqrt{z+1}} { 1 \\over z - \\I} - {\\sqrt{-\\I -1} \\sqrt{-\\I+1} \\I \\over 4\\sqrt{z-1} \\sqrt{z+1}} { 1 \\over z + \\I} + {C \\over \\sqrt{z-1} \\sqrt{z+1}}\n",
"$$\n",
"\n",
"But we want something stronger: that $\\phi(z) = O(z^{-2})$. To accomplish this, we need to choose $C$. Fortunately, I made the problem easy as every term apart from the last one is already $O(z^{-2})$, so choose $C = 0$:"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-0.2071067812011923 + 0.0im"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"φ = z -> 1/(2*(1+z^2)) + \n",
" sqrt(im-1)sqrt(im+1)*im/(4sqrt(z-1)sqrt(z+1))*1/(z-im) - \n",
" sqrt(-im-1)sqrt(-im+1)*im/(4sqrt(z-1)sqrt(z+1))*1/(z+im)\n",
"\n",
"φ(1E5)*1E5^2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see also that it has the right jump:"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.9900990099009901 + 0.0im"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"φ(0.1+0.0im) + φ(0.1-0.0im)"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.9900990099009901"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"1/(1+0.1^2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 1.3\n",
"\n",
"1. From the Hilbert formula, we know that the general solution of $\\HH u = f$ is \n",
"$$\n",
" u(x) = {-1 \\over \\sqrt{1 - x^2}}\\HH \\left[{ f(\\diamond) \\sqrt{1-\\diamond^2} }\\right](x) - {C \\over \\sqrt{1-x^2}}\n",
"$$ \n",
"Plugging in $f(x) = x/\\sqrt{1-x^2}$ means we need to calculate\n",
"$$\n",
"\\HH \\left[{\\diamond}\\right](x)\n",
"$$\n",
"We do so by first finding the Cauchy transform. Consider\n",
"$$\n",
"\\phi_1(z) = z \\CC 1(z) = z {\\log(z-1) - \\log(z+1) \\over 2 \\pi \\I}\n",
"$$\n",
"This has the right jump:\n",
"$$\n",
"\\phi_1^+(x) - \\phi_1^-(x) = x\n",
"$$\n",
"but doesn't decay at $\\infty$:\n",
"\\begin{align*}\n",
"x {\\log(x-1) - log(x+1) \\over 2 \\pi \\I} = x {\\log x + \\log(1-1/x) - \\log x -\\log(1+1/x) \\over 2 \\pi \\I} \\\\\n",
"= -x { 2 \\over x 2 \\pi \\I} = -{1 \\over \\I \\pi}\n",
"\\end{align*}\n",
"But this means that\n",
"$$\n",
"\\phi(z) = \\phi_1(z) + {1 \\over \\I \\pi} = z {\\log(z-1) - \\log(z+1) \\over 2 \\pi \\I} + {1 \\over \\I \\pi}\n",
"$$\n",
"Decays and has the right jump, hence is $\\CC[\\diamond](z)$. "
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.013174970881571598 - 0.00098617598822645im, 0.013174970881571569 - 0.0009861759882264232im)"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = Fun()\n",
"z = 2.0+2.0im\n",
"cauchy(t, z), z*(log(z-1)-log(z+1))/(2π*im) + 1/(im*π)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Therefore, we have\n",
"$$\n",
"\\HH[\\diamond](x) = \\I (\\CC^+ + \\CC^-) \\diamond(x) = x {\\log(1-x) - \\log(1+x) \\over \\pi} +{2 \\over \\pi}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.6302322257442835, 0.6302322257442834)"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = 0.1\n",
"hilbert(t,x), x*(log(1-x)-log(1+x))/(π) + 2/π"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Therefore, we get\n",
"$$\n",
"u(x) = - {x(\\log(1-x) - \\log(1+x))+2 \\over \\pi\\sqrt{1-x^2}} - {C \\over \\sqrt{1-x^2}}\n",
"$$\n",
"This can be verified in Mathematica via\n",
"```mathematica\n",
"NIntegrate[-((\n",
" x (Log[1 - x] - Log[1 + x]) + \n",
" 2)/(π Sqrt[1 - x^2] (x - 0.1))), {x, -1, 0.1, 1}, \n",
" PrincipalValue -> True]/π\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"2. Following the procedure of multiplying $\\CC[\\sqrt{1-\\diamond^2}](z)$ by $1/(2+z)$ and subtracting off the pole at $z=-2$, we first find:\n",
"$$\n",
"\\CC\\left[{\\sqrt{1-\\diamond^2} \\over 2+\\diamond}\\right](z) = {\\sqrt{z-1}\\sqrt{z+1} - z \\over 2\\I(2+z)} -{\\sqrt{-2-1}_+ \\sqrt{-2+1}_+ +2 \\over 2\\I(2+z)} = {\\sqrt{z-1}\\sqrt{z+1} - z \\over 2\\I(2+z)} +{ \\sqrt{3} -2 \\over 2\\I(2+z)} \n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.03230545315801244 + 0.032449695183223826im, 0.032305453158012455 + 0.032449695183223784im)"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = Fun()\n",
"z = 2.0+2.0im\n",
"cauchy(sqrt(1-t^2)/(2+t), z), (sqrt(z-1)sqrt(z+1)-z)/(2im*(2+z)) + (sqrt(3)-2)/(2im*(z+2))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Therefore, calculating $\\I(\\CC^+ + \\CC^-)$ we find that\n",
"$$\n",
"\\HH\\left[{\\sqrt{1-\\diamond^2} \\over 2+\\diamond}\\right](z) = -{x\\over 2+x} +{ \\sqrt{3} -2 \\over 2+x} \n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(-0.17521390115767732, -0.17521390115767752)"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = 0.1\n",
"hilbert(sqrt(1-t^2)/(2+t), x), -x/(2+x)+(sqrt(3)-2)/(2+x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Thus the general solution is\n",
"$$\n",
"u(x) = {1 \\over \\sqrt{1-x^2}} \\left( {x\\over 2+x} -{ \\sqrt{3} -2 \\over 2+x} + C \\right)\n",
"$$\n",
"We need to choose $C$ so this is bounded at the right-endpoint: In other words, \n",
"$$\n",
"u(x) = {1 \\over \\sqrt{1-x^2}} \\left( {x\\over 2+x} -{ \\sqrt{3} -2 \\over 2+x} +{ \\sqrt{3} -3 \\over 3} \\right)\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"u = (t/(2+t)-(sqrt(3)-2)/(2+t)+(sqrt(3)-3)/3)/sqrt(1-t^2)\n",
"plot(u; ylims=(-5,5))"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.4761904761904754, 0.47619047619047616)"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = 0.1\n",
"hilbert(u,x) , 1/(2+x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 2.1\n",
"\n",
"Doing the change of variables $\\zeta = b s$ we have\n",
"$$\n",
"\\log(ab) = \\int_1^{ab} {\\D \\zeta \\over \\zeta} = \\int_{1/b}^a {\\D s \\over s}\n",
"$$\n",
"if $\\gamma$ does not surround the origin, we have\n",
"$$\n",
"0 = \\oint_\\gamma {\\D s \\over s} = \\left[\\int_1^{1/b} + \\int_{1/b}^a + \\int_a^1\\right] {\\D s \\over s}\n",
"$$\n",
"which implies \n",
"$$\n",
"\\log(ab) = \\left[-\\int_a^1 -\\int_1^{1/b} \\right] {\\D s \\over s} = \\log a - \\log {1 \\over b} = \\log a + \\log b\n",
"$$\n",
"Here's a picture:"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"log(a * b) = 1.6094379124341003 - 0.9272952180016122im\n",
"log(a) + log(b) = 1.6094379124341003 - 0.9272952180016123im\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a = 1.0+2.0im\n",
"b = -1.0-2.0im\n",
"\n",
"@show log(a*b)\n",
"@show log(a) + log(b)\n",
"\n",
"\n",
"scatter([real(1/b)], [imag(1/b)]; label=\"1/b\")\n",
"scatter!([real(a)], [imag(a)]; label=\"a\")\n",
"scatter!([0.0], [0.0]; label=\"0\")\n",
"plot!(Segment(1, 1/b) ∪ Segment(1/b, a) ∪ Segment(a, 1); label=\"contour\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If it surrounds the origin counbter-clockwise, that is, it has positive orientation, we have $2\\pi \\I = \\oint_\\gamma {\\D s \\over s}$, which shoes that \n",
"$$\n",
"\\log(ab) = 2 \\pi \\I - \\left[\\int_a^1 +\\int_1^{1/b} \\right] {\\D s \\over s} = \\log a + \\log b + 2\\pi \\I\n",
"$$\n",
"and a similar result when counter clockwise."
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"log(a * b) = 1.8444397270569681 - 2.819842099193151im\n",
"(log(a) + log(b)) - (2π) * im = 1.844439727056968 - 2.819842099193151im\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a = 1.0+2.0im\n",
"b = -2.0+2.0im\n",
"\n",
"@show log(a*b)\n",
"@show log(a) + log(b) - 2π*im\n",
"scatter([real(1/b)], [imag(1/b)]; label=\"1/b\")\n",
"scatter!([real(a)], [imag(a)]; label=\"a\")\n",
"scatter!([0.0], [0.0]; label=\"0\")\n",
"plot!(Segment(1, 1/b) ∪ Segment(1/b, a) ∪ Segment(a, 1); label=\"contour\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If the contour passes through the origin, there are three possibility: \n",
"1. $[a,1]$ contains zero, hence $a < 0$ \n",
"2. $[1,1/b]$ contains zero, hence $b < 0$\n",
"3. $[1/b, a]$ contains zero, which can only be true if $a b < 0$ by considering the equation of the line segment.\n",
"\n",
"1\\. In the case where $a < 0$ and $b < 0$ (and hence $a b > 0$), perturbing $a$ above and $b$ below or vice versa avoids $\\gamma$ winding around zero, so we have\n",
"$$\n",
"\\log(a b) = \\log_+ a + \\log_- b = \\log_- a + \\log_+ b = \\log_+ a + \\log_+ b - 2 \\pi \\I = \\log_- a + \\log_- b + 2 \\pi \\I\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"log(a * b) = 1.791759469228055\n",
"log(a + 0.0im) + log(b - 0.0im) = 1.791759469228055 + 0.0im\n",
"log(a - 0.0im) + log(b + 0.0im) = 1.791759469228055 + 0.0im\n",
"log(a - 0.0im) + log(b - 0.0im) + (2π) * im = 1.791759469228055 + 0.0im\n",
"(log(a + 0.0im) + log(b + 0.0im)) - (2π) * im = 1.791759469228055 + 0.0im\n"
]
}
],
"source": [
"a = -2.0\n",
"b = -3.0\n",
"\n",
"@show log(a*b)\n",
"@show log(a+0.0im) + log(b-0.0im)\n",
"@show log(a-0.0im) + log(b+0.0im)\n",
"@show log(a-0.0im) + log(b-0.0im) + 2π*im\n",
"@show log(a+0.0im) + log(b+0.0im) - 2π*im;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the case where $a < 0$ and $b > 0$, then $ a b < 0$, but we can perturb $a$ above/below to get\n",
"$$\n",
"\\log_\\pm(a b) = \\log_\\pm a + \\log b\n",
"$$\n",
"(and by symmetry, the equivalent holds for $b < 0$ and $a > 0$.)"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"log(a * b + 0.0im) = 1.791759469228055 + 3.141592653589793im\n",
"log(a + 0.0im) + log(b) = 1.791759469228055 + 3.141592653589793im\n",
"log(a * b - 0.0im) = 1.791759469228055 - 3.141592653589793im\n",
"log(a - 0.0im) + log(b) = 1.791759469228055 - 3.141592653589793im\n"
]
}
],
"source": [
"a = -2.0\n",
"b = 3.0\n",
"\n",
"@show log(a*b +0.0im)\n",
"@show log(a+0.0im) + log(b);\n",
"\n",
"@show log(a*b -0.0im)\n",
"@show log(a-0.0im) + log(b);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the case where $a < 0$, if $\\Im b > 0$ we can perturb $a$ below so that $\\gamma$ does not contain zero, giving us\n",
"$$\n",
"\\log(ab) = \\log_- a + \\log b\n",
"$$\n",
"similarly, if $\\Im b < 0$ we can perturb $a$ above."
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"log(a * b) = 1.8444397270569681 - 2.819842099193151im\n",
"log(a - 0.0im) + log(b) = 1.8444397270569683 - 2.819842099193151im\n",
"log(a * b) = 1.8444397270569681 - 2.819842099193151im\n",
"log(a + 0.0im) + log(b) = 1.8444397270569683 + 3.4633432079864352im\n"
]
}
],
"source": [
"a = -2.0\n",
"b = 3.0 + im\n",
"\n",
"@show log(a*b)\n",
"@show log(a-0.0im) + log(b);\n",
"\n",
"b = 3.0 + im;\n",
"@show log(a*b)\n",
"@show log(a+0.0im) + log(b);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"2\\. In this case, swap the role of $a$ and $b$ and use the answers for $a < 0$. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"3\\. Finally, we have the case $a b < 0$ and neither $a$ nor $b$ is real. Note that\n",
"$$\n",
"ab = (a_x + \\I a_y) (b_x + \\I b_y) = a_x b_x - a_y b_y + \\I(a_x b_y + a_yb_x)\n",
"$$\n",
"It follows if $b_x > 0$ we have\n",
"$$\n",
"(ab)_+ = a_+ b\n",
"$$\n",
"and if $b_x < 0$ we have\n",
"$$\n",
"(ab)_+ = a_- b\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"log(a * b + eps() * im) = 0.6931471805599453 + 1.5707963267948966im\n",
"log((a + eps() * im) * b) = 0.6931471805599453 + 1.5707963267948968im\n"
]
},
{
"data": {
"text/plain": [
"0.6931471805599453 + 1.5707963267948968im"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a = 1.0 + 1.0im\n",
"b = 1.0 + 1.0im\n",
"@show log(a*b + eps()im) \n",
"@show log((a+eps()im)*b)"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"log(a * b + eps() * im) = 0.6931471805599453 + 3.141592653589793im\n",
"log((a - eps() * im) * b) = 0.6931471805599452 + 3.141592653589793im\n"
]
},
{
"data": {
"text/plain": [
"0.6931471805599452 + 3.141592653589793im"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a = 1.0 + 1.0im\n",
"b = -1.0 + 1.0im\n",
"@show log(a*b + eps()im) \n",
"@show log((a-eps()im)*b)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can use this perturbation to reduce to the previous cases. For example, if $a = 1 + \\I$ and $b = -1 + \\I$, pertubing $ab$ above causes $a$ to be perturbed above, which causes the contour to surround the origin clockwise, hence we have\n",
"$$\n",
"\\log_+(ab) = \\log(a)_+b = \\log a b - 2 \\pi \\I\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"log(a * b - eps() * im) = 0.6931471805599453 - 3.141592653589793im\n",
"(log(a) + log(b)) - (2π) * im = 0.6931471805599453 - 3.141592653589793im\n"
]
}
],
"source": [
"a = 1.0 + 1.0im\n",
"b = -1.0 + 1.0im\n",
"@show log(a*b - eps()im) \n",
"@show log(a)+log(b)-2π*im;"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"## Problem 2.2\n",
"\n",
"Use the contour $\\gamma(t) = 1 + t(1-z)$ to reduce it to a normal integral:\n",
"$$\\overline{\\log z } = \\overline{\\int_1^z {1 \\over \\zeta} \\D \\zeta} = \\overline{\\int_0^1 {(z-1) \\over 1+(z-1) t} \\dt} = \\int_0^1 {(\\bar z-1) \\over 1+(\\bar z-1) t} \\dt = \\int_1^{\\bar z} {\\D \\zeta \\over \\zeta} = \\log \\bar z.$$\n",
"We then have, since the contour from $1$ to $1/(\\bar z)$ to $z$ never surrounds the origin since both $\\Im z$ and $\\Im 1/(\\bar z)$ have the same sign, we have \n",
"$$\n",
"2 \\Re \\log z = \\log z + \\overline{\\log z} = \\log z + \\log \\bar z = \\log z \\bar z = \\log |z|^2 = 2 \\log |z|\n",
"$$\n",
"On the other hand, we have, where the contour of integration is chosen to be to the right of zero and then we do the change of variables $\\zeta = |z| \\E^{\\I \\theta}$\n",
"$$\n",
"2 \\Im \\log z = \\log z - \\log \\bar z = \\int_{\\bar z}^z {\\D \\zeta \\over \\zeta} = \\I \\int_{-\\arg z}^{\\arg z} \\D \\theta = 2 \\I \\arg z\n",
"$$\n",
"\n",
"\n",
"## Problem 2.3\n",
"\n",
"We first show that it is analytic on $(-\\infty,0)$. To do this, we need to show that the limit from above equals the limit from below: for $x < 0$ we have\n",
"$$\\log_1^+ x -\\log_1^- x = \\log_+x - \\log_- x +2 \\pi \\I = 0$$\n",
"Then for $x > 0$ and using $\\log_1^\\pm(x) = \\lim_{\\epsilon\\rightarrow 0} \\log(x\\pm \\I \\epsilon)$ we find\n",
"$$\n",
"\\log_1^+(x) - \\log_1^-x = \\log x- \\log x- 2\\pi \\I = -2 \\pi \\I\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"_Demonstration_ Here we see that the following is the analytic continuation:"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"#19 (generic function with 1 method)"
]
},
"execution_count": 43,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"log1 = z -> begin\n",
" if imag(z) > 0 \n",
" log(z)\n",
" elseif imag(z) == 0 && real(z) < 0\n",
" log(z + 0.0im)\n",
" elseif imag(z) < 0\n",
" log(z) + 2π*im\n",
" else\n",
" error(\"log1 not defined on real axis\")\n",
" end\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 48,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"phaseplot(-3..3, -3..3, log1)"
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.6931471805599453 + 3.141592653589793im"
]
},
"execution_count": 49,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"log1(-2.0+0.0im)"
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.6931471805599453 + 3.141592653589793im"
]
},
"execution_count": 50,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"log1(-2.0-0.0im)"
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.6931471805599453 + 3.141592653589793im"
]
},
"execution_count": 51,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"log1(-2.0)"
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.6931471805599453 + 1.1102230246251565e-16im"
]
},
"execution_count": 52,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"log1(2.0+eps()im)"
]
},
{
"cell_type": "code",
"execution_count": 53,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.6931471805599453 + 6.283185307179586im"
]
},
"execution_count": 53,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"log1(2.0-eps()im)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 3.1\n",
"\n",
"1. Because it's absolutely integrable, we can exchange derivatives and integrals to determine\n",
"$$\n",
"{\\D \\CC f \\over \\D z} = {1 \\over 2 \\pi \\I} \\oint {f(\\zeta)\\over (\\zeta-z)^2} \\D \\zeta\n",
"$$\n",
"2. There are two different possible approaches:\n",
" - the subtract and add back in technique: since $f$ is analytic for $z$ near $\\zeta$, we can write\n",
"$$\n",
"\\CC f(z) = {1 \\over 2 \\pi \\I} \\oint {f(\\zeta)-f(z) \\over (\\zeta-z)^2} \\D \\zeta + f(z) \\CC 1(z)\n",
"$$\n",
"Therefore\n",
"$$\n",
"\\CC^+ f(\\zeta) - \\CC^- f(\\zeta) = f(\\zeta)( \\CC^+ 1(\\zeta) - \\CC^- 1(\\zeta))\n",
"$$\n",
"But we know (using Cauchy's integral formula / Residue calculus)\n",
"$$\n",
"\\CC 1(z) = \\begin{cases} 1 & |z| < 1 \\\\\n",
" 0 & |z| > 0\n",
" \\end{cases}\n",
"$$\n",
"hence $(\\CC^+-\\CC^-)1(\\zeta) = 1$\n",
" - Since $f$ is analytic, we have for any radius $R > 1$ but inside the annulus\n",
"$$\n",
"\\CC^+ f(\\zeta) = {1 \\over 2 \\pi \\I} \\oint_{|z| = R} {f(\\zeta) \\over \\zeta -z} \\D\\zeta\n",
"$$\n",
"Similarly, for $\\CC^-f(\\zeta)$ with any radius $r < 1$ but inside the annulus. Therfore,\n",
"$$\n",
"\\CC^+ f(\\zeta) - \\CC^- f(\\zeta) = {1 \\over 2 \\pi \\I} \\left[\\oint_{|z| = R} - \\oint_{|z| = r} \\right] {f(\\zeta) \\over \\zeta -z} \\D\\zeta\n",
"$$\n",
"Deforming the contour and using Cauchy-integral formula gives the result. \n",
"3. This follow since ${1 \\over \\zeta - z} \\rightarrow 0$ uniformly.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 3.2\n",
"\n",
"Suppose we have another solution $\\phi$ and consider $\\psi(z) = \\phi(z) - \\CC f(z)$. Then on the circle we have\n",
"$$\n",
"\\psi_+(\\zeta) - \\psi_-(\\zeta) = \\phi_+(\\zeta) - \\CC_+f(\\zeta) - \\phi_-(\\zeta) + \\CC_+f(\\zeta) = f(\\zeta)-f(\\zeta) = 0\n",
"$$\n",
"Thus $\\psi$ is entire, and since it decays at infinity, it must be zero by Liouville's theorem.\n",
"\n",
"## Problem 3.3\n",
"\n",
"When $k \\geq 0$, we have from 3.1 and 3.2\n",
"$$\n",
"\\CC[\\diamond^k](z) =\\begin{cases}\n",
" z^k & |z| < 0 \\\\\n",
" 0 & |z| > 0\n",
" \\end{cases}\n",
"$$\n",
"when $k < 0$ since\n",
"$$\\CC[\\diamond^k]^+(\\zeta) - \\CC[\\diamond^k]^-(\\zeta) = \\zeta^k - 0 = \\zeta^k$$. we similarly have \n",
"$$\n",
"\\CC[\\diamond^k](z) =\\begin{cases}\n",
"0 & |z| < 0 \\\\\n",
" -z^k & |z| > 0\n",
" \\end{cases}\n",
"$$\n",
"\n",
"Therefore,\n",
"$$\n",
"\\Im \\CC^-[\\diamond^k](\\zeta) =\\begin{cases}\n",
"0 & k \\geq 0 \\\\\n",
" -{\\zeta^k - \\zeta^{-k} \\over 2 \\I} & k < 0\n",
" \\end{cases}\n",
"$$\n",
"and\n",
"$$\n",
"\\Re \\CC^-[\\diamond^k](\\zeta) =\\begin{cases}\n",
"0 & k \\geq 0 \\\\\n",
" -{\\zeta^k + \\zeta^{-k} \\over 2} & k < 0\n",
" \\end{cases}\n",
"$$\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 3.4\n",
"\n",
"Express the solution outside the circle as\n",
"$$\n",
"v(x,y) = \\Im ( \\E^{-\\I \\theta} z + \\CC f(z))\n",
"$$\n",
"for a to-be-determined $f$. On the circle, this reduces to\n",
"$$\n",
"\\Im \\CC^- f(\\zeta) = -\\cos \\theta {\\zeta - \\zeta^{-1} \\over 2 \\I} + \\sin \\theta {\\zeta + \\zeta^{-1} \\over 2}\n",
"$$\n",
"Unlike the real case, we can include imaginary coefficients, thus the solution is\n",
"$$\n",
"f(\\zeta) = (\\cos \\theta + \\I \\sin \\theta) \\zeta^{-1}\n",
"$$\n",
"and thus the full solution is\n",
"$$\n",
"v(x,y) = \\Im ( \\E^{-\\I \\theta} z - \\E^{\\I \\theta} z^{-1}))\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 54,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"#21 (generic function with 1 method)"
]
},
"execution_count": 54,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"θ = 0.1\n",
"v = (x,y) -> x^2 + y^2 < 1 ? 0 : imag(exp(-im*θ) * (x+im*y) + exp(im*θ) * (x+im*y)^(-1))"
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 58,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"xx = yy = -5:0.01:5\n",
"\n",
"contour(xx, yy, v.(xx', yy); nlevels=100)\n",
"plot!(Circle(); color=:black, legend=false)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 4.1\n",
"\n",
"$z^\\alpha$ has the limits $z_\\pm^\\alpha = \\E^{\\pm \\I \\pi \\alpha} |z|^\\alpha$, thus choose $\\alpha = -{\\theta \\over 2\\pi}$ where if we take $0 < \\theta < 2\\pi$ we have $0 < \\alpha < 1$ (the case $\\theta = 0$ and $\\theta = \\pi$ are covered by the Cauchy transform, that is ). Then consider\n",
"$$\n",
"\\kappa(z) = (z-1)^{-\\alpha} (z+1)^{\\alpha-1}\n",
"$$\n",
"which has weaker than pole singularities and satisfies $\\kappa(z) \\sim z^{-1}$.\n",
"For $-1 < x < 1$ it has the right jump\n",
"$$\n",
"\\kappa_+(x) = (x-1)_+^\\alpha (x+1)^{1 - \\alpha} = \\E^{\\I \\pi \\alpha} (1-x)^\\alpha (x+1)^{1 - \\alpha} = \\E^{2\\I \\pi \\alpha} (x-1)_-^\\alpha (x+1)^{1 - \\alpha} = \\E^{2 \\I \\pi \\alpha} \\kappa_-(x)= \\E^{\\I \\theta} \\kappa_-(x)\n",
"$$\n",
"and for $x < -1$ it has the jump\n",
"$$\n",
"\\kappa_+(x) = (x-1)_+^\\alpha (x+1)_+^{1 - \\alpha} = \\E^{\\I \\pi \\alpha}\\E^{\\I \\pi (1-\\alpha)} (1-x)^\\alpha (-1-x)^{1 - \\alpha} = \\kappa_-(x)\n",
"$$\n",
"hence $\\kappa$ is analytic. \n",
"\n",
"We need to show this times a constant spans the entire space. Suppose we have another solution $\\tilde \\kappa$ and consider $r(z) = {\\tilde \\kappa(z) \\over \\kappa(z)}$. Note by construction that $\\kappa$ has no zeros. Then\n",
"$$\n",
"r_+(x) = {\\tilde\\kappa_+(x) \\over \\kappa_+(x)} = {\\tilde\\kappa_-(x) \\over \\kappa_-(x)} = r_-(x)\n",
"$$\n",
"hence $r$ is analytic on $(-1,1)$. It has weaker than pole singularities because $\\kappa(z)^{-1}$ is actually bounded at $\\pm 1$. Therefore $r$ is bounded and entire, and thus must be a constant $r(z) \\equiv r$, and thence $\\tilde \\kappa(z) = r \\kappa(z)$."
]
},
{
"cell_type": "code",
"execution_count": 59,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0 + 1.1102230246251565e-16im"
]
},
"execution_count": 59,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"θ =2.3\n",
"α = -θ/(2π)\n",
"κ = z -> (z-1)^(-α)*(z+1)^(α-1)\n",
"κ(0.1+0.0im) - exp(im*θ)*κ(0.1-0.0im)"
]
},
{
"cell_type": "code",
"execution_count": 60,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.00982876598532333"
]
},
"execution_count": 60,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"κ(100.0)"
]
},
{
"cell_type": "code",
"execution_count": 62,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 62,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"phaseplot(-3..3, -3..3, κ)"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"## Problem 4.2\n",
"\n",
"We want to mimic the solution of $\\phi_+(x) + \\phi_-(x)$. So take\n",
"$$\n",
"\\phi(z) = \\kappa(z) \\CC\\br[{f \\over \\kappa_+}](z) =\\E^{-\\I \\theta/2} (z-1)^{-\\alpha}(z+1)^{\\alpha-1} \\CC[f (1-x)^{\\alpha}(1+x)^{1-\\alpha}](z)\n",
"$$\n",
"This has the jump\n",
"\\begin{align*}\n",
"\\phi_+(x) - \\E^{\\I \\theta}\\phi_-(x) = \\kappa_+(z) \\CC_+\\br[{f \\over \\kappa_+}](x) - \\E^{\\I \\theta}\\kappa_-(z)\\CC_-\\br[{f \\over \\kappa_+}](x) = \\kappa_+(x) \\pr({\\CC_+\\br[{f \\over \\kappa_+}](x) - \\CC_-\\br[{f \\over \\kappa_+}](x) }) = f(x)\n",
"\\end{align*}\n",
"\n",
"Thus the general solution is $\\phi(z) + C \\kappa(z)$."
]
},
{
"cell_type": "code",
"execution_count": 63,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"8.881784197001252e-16 - 2.220446049250313e-16im"
]
},
"execution_count": 63,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"θ =2.3\n",
"α = -θ/(2π)\n",
"κ = z -> (z-1)^(-α)*(z+1)^(α-1)\n",
"\n",
"\n",
"x = Fun()\n",
"κ₊ = exp(im*θ/2)*(1-x)^(-α)*(x+1)^(α-1)\n",
"\n",
"f = Fun(exp)\n",
"\n",
"z = 2+im\n",
"φ = z -> κ(z)*cauchy(f/κ₊, z)\n",
"\n",
"φ(0.1+0.0im)-exp(im*θ)*φ(0.1-0.0im) - f(0.1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 4.3\n",
"\n",
"Note for $x < 0$\n",
"$$\n",
"x_+^{\\I \\beta} = \\E^{\\I \\beta \\log_+ x} = \\E^{\\I \\beta \\log_- x - 2\\pi y} = \\E^{-2 \\pi \\beta} x_-^{\\I \\beta}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 64,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-3.3881317890172014e-21 + 1.0842021724855044e-19im"
]
},
"execution_count": 64,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"β = 2.3;\n",
"x = -2.0\n",
"(x+0.0im)^(im*β) - exp(-2π*β)*(x-0.0im)^(im*β)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We actually have bounded (oscillatory) growth near zero since \n",
"$$\n",
"|\\E^{\\I \\beta \\log z}| = |\\E^{\\I \\beta \\log |z|} \\E^{-\\beta \\arg z}| = \\E^{-\\beta \\arg z}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 67,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 67,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"xx = yy = -1:0.01:1\n",
"contourf(xx, yy, abs.((xx' .+ im.*yy).^(im*β)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Thus if we write $c = r \\E^{\\I \\theta}$ for $0 < \\theta < 2 \\pi$ and define $\\alpha = -{\\theta \\over 2 \\pi} + \\I{\\log r \\over 2 \\pi}$ we can write the solution to 4.1 as\n",
"$$\n",
"\\kappa(z) = (z-1)^{-\\alpha} (z+1)^{\\alpha-1}\n",
"$$\n",
"\n",
"The same arguments as before then proceed. and the solution to 4.3 is \n",
"$$\n",
"\\phi(z) = \\kappa(z) \\CC\\br[{f \\over \\kappa_+}](z)+ C \\kappa(z)\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 68,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0 + 1.1102230246251565e-16im"
]
},
"execution_count": 68,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"θ =2.3\n",
"α = -θ/(2π)\n",
"κ = z -> (z-1)^(-α)*(z+1)^(α-1)\n",
"κ(0.1+0.0im) - exp(im*θ)*κ(0.1-0.0im)"
]
},
{
"cell_type": "code",
"execution_count": 69,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.220446049250313e-16 + 2.220446049250313e-16im"
]
},
"execution_count": 69,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"r = 2.4\n",
"θ = 2.1\n",
"c = r*exp(im*θ)\n",
"\n",
"α = -θ/(2π) + im*log(r)/(2π)\n",
"\n",
"κ = z -> (z-1)^(-α)*(z+1)^(α-1)\n",
"\n",
"\n",
"κ(0.1+0.0im)-c*κ(0.1-0.0im)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 5.1\n",
"\n",
"1. It is a product of functions analytic off $(-\\infty,1]$ hence is analytic off $(-\\infty,1]$, and we just have to check that it has no jump on $(-\\infty,-1)$ and $(-a,a)$. This follows via, for $x < -1$:\n",
"$$\n",
"\\kappa_+(x) = {1 \\over \\I^4 \\sqrt{1-x}\\sqrt{-1-x}\\sqrt{a-x}\\sqrt{-a-x}} = {1 \\over (-\\I)^4 \\sqrt{1-x}\\sqrt{-1-x}\\sqrt{a-x}\\sqrt{-a-x}} = \\kappa_-(x)\n",
"$$\n",
"and for $-a < x < a$ we have\n",
"$$\n",
"\\kappa_+(x) = {1 \\over \\I^2 \\sqrt{1-x}\\sqrt{x+1}\\sqrt{a-x}\\sqrt{x+a}} = {1 \\over (-\\I)^2 \\sqrt{1-x}\\sqrt{1+x}\\sqrt{a-x}\\sqrt{-a-x}} = \\kappa_-(x)\n",
"$$\n",
"2. This follows via the usual arguments: for $a < x < 1$ we have:\n",
"$$\n",
"\\kappa_+(x) = {1 \\over \\I \\sqrt{1-x}\\sqrt{x+1}\\sqrt{x-a}\\sqrt{x+a}} = -\\kappa_-(x)\n",
"$$\n",
"and for $-1 < x < -a$ we have\n",
"$$\n",
"\\kappa_+(x) = {1 \\over \\I^3 \\sqrt{1-x}\\sqrt{x+1}\\sqrt{x-a}\\sqrt{x+a}} = -\\kappa_-(x)\n",
"$$\n",
"3. This has at most square singularitie4s which are weaker than poles\n",
"4. $\\kappa(z) = {1 \\over z^2 \\sqrt{1-1/z}\\sqrt{1-a/z} \\sqrt{1+a/z} \\sqrt{1+1/z}} \\sim {1 \\over z^2} \\rightarrow 0$. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 5.2\n",
"\n",
"Ah, this is a trick question! Note that $z \\kappa(z) \\sim z^{-1} = O(z)$ and satisfies all the other properties. Thus consider any other solution $\\tilde \\kappa(z)$ and write\n",
"$$\n",
"r(z) = {\\tilde \\kappa(z) \\over \\kappa(z) }\n",
"$$\n",
"This has trivial jumps and hence is entire: for example, on $(a,1)$ we have\n",
"$$\n",
"r_+(x) = {\\tilde \\kappa_+(x) \\over \\kappa_+(x) } = {-\\tilde \\kappa_-(x) \\over -\\kappa_-(x) } = r_-(x)\n",
"$$\n",
"But since $\\kappa \\sim O(z^{-2})$ we only know that $\\kappa$ has at most $O(z)$ growth, hence it can be any first degree polynomial. Therefore, the space of all solutions is in fact two-dimensional: $\\psi(z) = (A + Bz) \\kappa(z)$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 5.3\n",
"\n",
"Here we mimick the usual solution techniques and propose:\n",
"$$\n",
"\\phi(z) = \\kappa(z) \\CC\\br[{f \\over \\kappa_+}](z) + (A+Bz) \\kappa(z)\n",
"$$\n",
"A quick check confirms it has the right jumps:\n",
"$$\n",
"\\phi_+(x) = \\kappa_+(x) \\CC_+\\br[{f \\over \\kappa_+}](x) + (A+Bx) \\kappa_+(x) =\\kappa_+(x) \\left({f(x) \\over \\kappa_+(x)} + \\CC_-\\br[{f \\over \\kappa_+}](x)\\right) - (A+Bx) \\kappa_-(x) = -\\phi_-(x) + f(x)\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"## Problem 6.1\n",
"\n",
"Let's first do a plot and histogram. Here we see the right scaling is $N^{1/4}$:"
]
},
{
"cell_type": "code",
"execution_count": 71,
"metadata": {},
"outputs": [],
"source": [
"V = x -> x^4\n",
"Vp = x -> 4x^3\n",
"\n",
"N = 100\n",
"λ⁰ = randn(N) # initial location\n",
"prob = ODEProblem((λ,t,_) -> Float64[sum(1 ./(λ[k] .- λ[[1:k-1;k+1:end]])) - Vp(λ[k]) for k=1:N], λ⁰, (0.0, 10.0))\n",
"λ = solve(prob; reltol=1E-6);"
]
},
{
"cell_type": "code",
"execution_count": 73,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"┌ Info: Saved animation to \n",
"│ fn = /Users/sheehanolver/Documents/Teaching/Imperial/M3M6 Methods of Mathematical Physics/2018/tmp.gif\n",
"└ @ Plots /Users/sheehanolver/.julia/packages/Plots/rmogG/src/animation.jl:90\n"
]
},
{
"data": {
"text/html": [
"![](\"tmp.gif?0.021308146046194576)
\" />"
],
"text/plain": [
"Plots.AnimatedGif(\"/Users/sheehanolver/Documents/Teaching/Imperial/M3M6 Methods of Mathematical Physics/2018/tmp.gif\")"
]
},
"execution_count": 73,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@gif for t=0.0:0.05:7.0\n",
" scatter(λ(t)/N^(1/4) ,zeros(N); label=\"charges\", xlims=(-2,2), title=\"t = $t\")\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The limiting distribution has the following form:"
]
},
{
"cell_type": "code",
"execution_count": 74,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 74,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"histogram(λ(10.0)/N^(1/4); nbins=30, normalize=true, label=\"histogram of charges\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"We want to solve\n",
"$$\n",
"{\\D \\lambda_k \\over \\D t} = \\sum_{j=1 \\atop j \\neq k}^N {1 \\over \\lambda_k -\\lambda_j} - 4\\lambda_k^3\n",
"$$\n",
"Rescale via $\\mu_k = {\\lambda_k \\over N^{1/4}}$ gives\n",
"$$\n",
"0 = N^{1/4} {\\D \\mu_k \\over \\D t} = N^{-1/4} \\sum_{j=1 \\atop j \\neq k}^N {1 \\over \\mu_k -\\mu_j} - 4 N^{3/4} \\mu_k^3\n",
"$$\n",
"or in other words\n",
"$$\n",
"0 = N^{-1/2} {\\D \\mu_k \\over \\D t} = {1 \\over N} \\sum_{j=1 \\atop j \\neq k}^N {1 \\over \\mu_k -\\mu_j} - 4 \\mu_k^3\n",
"$$\n",
"We can now formally let $N\\rightarrow \\infty$ to get our equation\n",
"$$\n",
"\\dashint_{-b}^b {w(t) \\over x-t} \\dt = 4 x^3\n",
"$$\n",
"where I've used symmetry to assume that the interval is symmetric. \n",
"We want to find $w$ and $b$ so that this equations holds true and $w$ is a bounded probability density:\n",
"1. $w(x) >0$ for $-b < x < b$\n",
"2. $\\int w(x) \\dx = 1$\n",
"3. $w$ is bounded\n",
"\n",
"\n",
"Our equation is equivalent to\n",
"$$\n",
"\\HH_{[-b,b]} w(x) = -{4 x^3 \\over \\pi}\n",
"$$\n",
"recall the inversion formula\n",
"$$\n",
" u(x) = {-1 \\over \\sqrt{b^2 - x^2}}\\HH \\left[{ f(\\diamond) \\sqrt{b^2-\\diamond^2} }\\right](x) - {C \\over \\sqrt{b^2-x^2}}\n",
"$$ \n",
"In our case $f(x) = -{4 x^3 \\over \\pi}$ and we use\n",
"$$\n",
"\\sqrt{z-b}\\sqrt{z+b} = z\\sqrt{1-b/z}\\sqrt{1+b/z} = z - {b^2 \\over 2 z} -{b^4 \\over 8z^3} + O(z^{-4})\n",
"$$\n",
"to determine\n",
"$$\n",
"2 \\I \\CC \\left[{ \\diamond^3 \\sqrt{b^2-\\diamond^2} }\\right](x) = z^3(\\sqrt{z-b}\\sqrt{z+b} -z + {b^2 \\over 2 z} + {b^4 \\over 8z^3}) = z^3\\sqrt{z-b}\\sqrt{z+b} -z^4 + {b^2z^2 \\over 2} + {b^4 \\over 8}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 75,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"71.6687198577313 + 40.090510492433154im"
]
},
"execution_count": 75,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"b = 5\n",
"x = Fun(-b .. b)\n",
"(2im)cauchy(x^3*sqrt(b^2-x^2),z)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Therefore,\n",
"\\begin{align*}\n",
"u(x) = {4\\I \\over \\pi \\sqrt{b^2 - x^2}}(\\CC^+ + \\CC^-) \\left[{ \\diamond^3 \\sqrt{b^2-\\diamond^2} }\\right](x) - {C \\over \\sqrt{b^2-x^2}} \\\\\n",
"{4 \\over \\pi \\sqrt{b^2 - x^2}}(-x^4+{b^2 x^2 \\over 2} + {b^4 \\over 8})- {C \\over \\sqrt{b^2-x^2}} \n",
"\\end{align*}\n",
"We choose $C$ so this is bounded, in particular, we get get the solution\n",
"$$\n",
"u(x) = {4 \\over \\pi \\sqrt{b^2 - x^2}}(-x^4+{b^2 x^2 \\over 2} + {b^4 \\over 2})\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 76,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-0.0012732395447351292"
]
},
"execution_count": 76,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"u = 4/(π*sqrt(b^2-x^2))*(-x^4 + b^2*x^2/2 + b^4/2)\n",
"hilbert(u, 0.1)"
]
},
{
"cell_type": "code",
"execution_count": 77,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-0.001273239544735163"
]
},
"execution_count": 77,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"-4*0.1^3/π"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"At least it looks right, we just need to get the right b:"
]
},
{
"cell_type": "code",
"execution_count": 78,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 78,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(u; label = \"b = $b\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We want to choose $b$ now so that this integrates to $1$. For example, this choice of $b$ is horrible:"
]
},
{
"cell_type": "code",
"execution_count": 79,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"937.5"
]
},
"execution_count": 79,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sum(u)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There's a nice trick: If $-2\\pi \\I \\CC u(z) \\sim {1 \\over z}$ then $\\int_{-b}^b u(x) \\dx = 1$. We know since\n",
"$$\n",
"{1 \\over \\sqrt{1+z}} = 1 - {z \\over 2} + {3z^2 \\over 8} - {5z^3 \\over 16} + O(z^4)\n",
"$$\n",
"\\begin{align*}\n",
" {-z^4 + {b^2z^2 \\over 2} + {b^4 \\over 2} \\over \\sqrt{z-b}\\sqrt{z+b}} &= {-z^3 + {b^2z \\over 2} + {b^4 \\over 2 z} \\over \\sqrt{1-b/z}\\sqrt{1+b/z}} \\\\\n",
" &= (-z^3 + {b^2z \\over 2} )(1 + {b \\over 2 z} + {3 b^2 \\over 8 z^2} + {5 b^3 \\over 16z^3} )(1 - {b \\over 2 z} + {3 b^2 \\over 8 z^2} - {5 b^3 \\over 16z^3}) + O(z^{-1}) \\\\\n",
" &=-z^3 + \\pr({b^2 \\over 2} +{b^2 \\over 4} + {b^2 \\over 2} - {3 b^2 \\over 8} -{3 b^2 \\over 8} ) z +O(z^{-1})\\\\\n",
" &= -z^3 + O(z^{-1})\n",
" \\end{align*}\n",
"Thus we know\n",
"$$\n",
"\\CC u(z) = {2\\I \\over \\pi} z^3 +{2 \\I \\over \\pi} {-z^4 + {b^2z^2 \\over 2} + {b^4 \\over 2} \\over \\sqrt{z-b}\\sqrt{z+b}}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 80,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"40.55106915762885 + 9.075717966750503im"
]
},
"execution_count": 80,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cauchy(u, z)"
]
},
{
"cell_type": "code",
"execution_count": 81,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.4908359390683472 - 0.0im"
]
},
"execution_count": 81,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"z =100.0im;\n",
"2im/π*z^3 + 2im/π*(-z^4 + b^2*z^2/2 +b^4/2)/(sqrt(z-b)sqrt(z+b))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"taking this one term further we find \n",
"$$\n",
" {-z^4 + {b^2z^2 \\over 2} + {b^4 \\over 2} \\over \\sqrt{z-b}\\sqrt{z+b}} = -z^3 +{3 b^4 \\over 8 z} + O(z^{-2})\n",
" $$\n",
"Hence we want to choose $b$ so that\n",
"$$\n",
"-2\\I\\pi \\CC u(z) = {3 b^4 \\over 2 z} \\sim {1 \\over z}\n",
"$$\n",
"in other words, $b = ({2 \\over 3})^{1/4}$"
]
},
{
"cell_type": "code",
"execution_count": 82,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.9999999999999997"
]
},
"execution_count": 82,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"b = (2/3)^(1/4)\n",
"x = Fun(-b .. b)\n",
"u = 4/(π*sqrt(b^2-x^2))*(-x^4 + b^2*x^2/2 + b^4/2)\n",
"sum(u)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And it worked!"
]
},
{
"cell_type": "code",
"execution_count": 83,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 83,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"histogram(λ(10.0)/N^(1/4); nbins=25, normalize=true, label=\"histogram of charges\")\n",
"plot!(u; label=\"u\",xlims=(-2,2))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.0.2",
"language": "julia",
"name": "julia-1.0"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.0.2"
}
},
"nbformat": 4,
"nbformat_minor": 2
}