{
"cells": [
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [],
"source": [
"using Plots, ComplexPhasePortrait, ApproxFun, SingularIntegralEquations, OscillatoryIntegrals, QuadGK\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\\rR{{\\rm R}}\n",
"\\def\\rL{{\\rm L}}\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\\Ei{{\\rm Ei}\\,}\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 4\n",
"\n",
"\n",
"\n",
"## Problem 1.1\n",
"\n",
"We want to solve\n",
"$$\n",
"\\int_{-1}^1 \\log |x-t| u(t) \\dt = {1 \\over x^2 + 1}\n",
"$$\n",
"differentiating and multiplying though by $1/\\pi$ we have\n",
"$$\n",
"{1 \\over \\pi} \\dashint_{-1}^1 {u(t) \\over t-x} u(t) = {2 x \\over \\pi (1+x^2)^2}\n",
"$$\n",
"\n",
"This is an inverse Hilbert problem so we know\n",
"$$\n",
"u(x) = -{1 \\over \\sqrt{1-x^2} } H[\\sqrt{1-t^2} f(t)] + {C \\over \\sqrt{1-x^2}}\n",
"$$\n",
"for \n",
"$$\n",
"f(x) = {2 x \\over \\pi (1+x^2)^2}.\n",
"$$\n",
"\n",
"We determine the Cauchy transform of $\\sqrt{1-x^2} f(x)$ using the usual methods: start with the ansatz\n",
"$$\n",
"\\phi(z) = {\\sqrt{z-1} \\sqrt{z+1} \\over 2 \\I} {2 z \\over \\pi (1+z^2)^2}\n",
"$$\n",
"This decays at infinity so we just need to remove the poleas at $\\pm \\I$. Here we determine:\n",
"$$\n",
"\\phi(z) = {\\sqrt{\\I -1 } \\sqrt{\\I + 1} \\over \\pi} \\left(-{1 \\over 4 (z-\\I)^2} + {\\I \\over 8 (z-\\I)} + O(1) \\right)\n",
"$$\n",
"and\n",
"$$\n",
"\\phi(z) = {\\sqrt{-\\I -1 } \\sqrt{1-\\I } \\over \\pi} \\left({1 \\over 4 (z+\\I)^2} + {\\I \\over 8 (z+\\I)} + O(1) \\right)\n",
"$$\n",
"Telling us that\n",
"$$\n",
"C[\\sqrt{1-t^2} f(t)](z) = {\\sqrt{z-1} \\sqrt{z+1} \\over 2 \\I} {2 z \\over \\pi (1+z^2)^2} - {\\sqrt{\\I -1 } \\sqrt{\\I + 1} \\over \\pi} \\left(-{1 \\over 4 (z-\\I)^2} + {\\I \\over 8 (z-\\I)} \\right) - {\\sqrt{-\\I -1 } \\sqrt{1-\\I } \\over \\pi} \\left({1 \\over 4 (z+\\I)^2} + {\\I \\over 8 (z+\\I)} \\right)\n",
"$$\n",
"Thus we have\n",
"$$\n",
"H[\\sqrt{1-t^2} f](x) = \\I (C^+ + C^-)[\\sqrt{1-t^2} f](x) = - {2\\I \\sqrt{\\I -1 } \\sqrt{\\I + 1} \\over \\pi} \\left(-{1 \\over 4 (x-\\I)^2} + {2\\I \\over 8 (x-\\I)} \\right) - {2\\I \\sqrt{-\\I -1 } \\sqrt{1-\\I } \\over \\pi} \\left({1 \\over 4 (x+\\I)^2} + {\\I \\over 8 (x+\\I)} \\right)\n",
"$$\n",
"\n",
"In other words, for\n",
"$$\n",
"\\tilde u(x) = {2\\I \\sqrt{\\I -1 } \\sqrt{\\I + 1} \\over \\pi \\sqrt{1-x^2}} \\left(-{1 \\over 4 (x-\\I)^2} + {2\\I \\over 8 (x-\\I)} \\right) - {2\\I \\sqrt{-\\I -1 } \\sqrt{1-\\I } \\over \\pi \\sqrt{1-x^2}} \\left({1 \\over 4 (x+\\I)^2} + {\\I \\over 8 (x+\\I)} \\right)\n",
"$$\n",
"we have\n",
"$$\n",
"u(x) = \\tilde u(x) + {C \\over \\sqrt{1-x^2}}.\n",
"$$\n",
"\n",
"\n",
"To find $C$ we impose the condition that\n",
"$$\n",
"\\int_{-1}^1 \\log|t| u(t) \\dt = 1\n",
"$$\n",
"We thus need to determine $C$. Recall from Lecture 19 that\n",
"$$\n",
"{1 \\over \\pi} \\int_{-1}^1 {\\log(z-x) \\over \\sqrt{1-x^2}} \\dx = 2 \\log(\\sqrt{z-1}+\\sqrt{z+1}) - 2 \\log 2\n",
"$$\n",
"Thus for $x \\in [-1,1]$ we have\n",
"$$\n",
"{1 \\over \\pi} \\int_{-1}^1 {\\log|x-t| \\over \\sqrt{1-t^2}} \\dt = 2 \\Re(\\log(\\I \\sqrt{1-x}+\\sqrt{x+1}) - 2 \\log 2)\n",
"$$\n",
"or in particular \n",
"$$\n",
"{1 \\over \\pi} \\int_{-1}^1 {\\log|x| \\over \\sqrt{1-x^2}} \\dx = \\log(1/2)\n",
"$$\n",
"Thus we want to solve\n",
"$$\n",
"\\int_{-1}^1 \\log|t| \\tilde u(t) \\dt + C {\\log(1/2) \\over \\pi} = 1\n",
"$$\n",
"i.e.\n",
"$$\n",
"C = (1 - \\int_{-1}^1 \\log|t| \\tilde u(t) \\dt) {1 \\over \\pi \\log(1/2)}\n",
"$$\n",
"\n",
"\n",
"\n",
"_Check derivation_\n",
"\n",
"Let's check the derivation. First we can calculate $u$ numerically:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-2.220446049250313e-16"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"L = SingularIntegral(0) : JacobiWeight(-0.5,-0.5,Chebyshev())\n",
"\n",
"\n",
"x = Fun()\n",
"g = (1/(x^2+1))\n",
"u = (π*L) \\ g\n",
"\n",
"π*logkernel(u, 0.1) - g(0.1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Differentiating it satisfies $H u = g'/(-\\pi)$:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"true"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"hilbert(u,0.1) ≈ g'(0.1)/(-π) ≈ (2*0.1)/(π*(1+0.1^2)^2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Cauchy transform also works:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"true"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f = g'/(-π)\n",
"\n",
"z = 1+im\n",
"ψ = z -> sqrt(z-1)sqrt(z+1)/(2im) * 2z/(π*(1+z^2)^2) - \n",
" sqrt(im-1)sqrt(im+1)/π * (-1/(4*(z-im)^2) + im/(8(z-im))) - \n",
" sqrt(-im-1)sqrt(1-im)/π * (1/(4*(z+im)^2) + im/(8(z+im)))\n",
"\n",
"cauchy(sqrt(1-x^2)*f, z) ≈ ψ(z)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Therefore the Hilbert transform is given by:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"true"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"H = x-> - 2im*sqrt(im-1)sqrt(im+1)/π * (-1/(4*(x-im)^2) + im/(8(x-im))) - \n",
" 2im*sqrt(-im-1)sqrt(1-im)/π * (1/(4*(x+im)^2) + im/(8(x+im)))\n",
"\n",
"hilbert(sqrt(1-x^2)*f)(0.1) ≈ H(0.1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We thus have an expression for $u$, we are just missing the constant:"
]
},
{
"cell_type": "code",
"execution_count": 81,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"true"
]
},
"execution_count": 81,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ũ = x -> -H(x)/sqrt(1-x^2)\n",
"\n",
"C = u(0.0) + H(0.0)\n",
"u(0.1) ≈ ũ(0.1) + C/sqrt(1-0.1^2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can find $C$ interms of the relevant integral, which we call $D$:"
]
},
{
"cell_type": "code",
"execution_count": 92,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(-0.3247204711377926 + 0.0im, -0.32472047136213555 + 0.0im)"
]
},
"execution_count": 92,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"D = quadgk(x -> x == 0 ? 0.0 : ũ(x)*log(abs(x)),-1,1)[1]\n",
"C, (1-D)/(π*log(1/2))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Problem 1.2\n",
"\n",
"Since $\\int_{-1}^1 x \\dx = 0$, we have no logarithmic growth at infinity. Thus by the formula in Lecture 19 we find for\n",
"$$\n",
"U(x) = \\int_x^1 t \\dt = {1 - x^2 \\over 2}\n",
"$$\n",
"that\n",
"$$\n",
"\\int_{-1}^1 \\log|z-x| x \\dx = \\Re{-2 \\I \\pi C U(z)}\n",
"$$\n",
"\n",
"So we just have to work out the Cauchy transform. Try as an ansatz\n",
"$$\n",
"(1 - z^2 \\over 2) {\\log(z-1) - \\log(z+1) \\over 2 \\pi \\I}\n",
"$$\n",
"We only need to remove the growth at infinity. We do so via:\n",
"$$\n",
"\\log(z-1) - \\log(z+1) = -{2 \\over z} + O(z^{-3})\n",
"$$\n",
"telling us\n",
"$$\n",
"C U(z) = {1 - z^2 \\over 2} {\\log(z-1) - \\log(z+1) \\over 2 \\pi \\I} + {z \\over 2 \\pi \\I}\n",
"$$\n",
"and\n",
"$$\n",
"\\int_{-1}^1 \\log|z-x| x \\dx = \\Re \\left({1 - z^2 \\over 2} (\\log(z-1) - \\log(z+1)) - z \\right)\n",
"$$\n",
"\n",
"We can confirm the formula:"
]
},
{
"cell_type": "code",
"execution_count": 93,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"true"
]
},
"execution_count": 93,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"z = 2+im; π*logkernel(x,z) ≈ real((1-z^2)/2 * (log(z-1)-log(z+1)) - z)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 2\n",
"\n",
"### Problem 2.1\n",
"\n",
"From Lecture 19 we know that we want to solve\n",
"\n",
"1. $v_{xx} + v_{yy} =0$ for $z \\notin [-1,1] \\cup \\I$\n",
"2. $v(z) \\sim \\log|z|$ as $z \\rightarrow \\infty$\n",
"3. $v(z) \\sim \\log|z-\\I|$ as $z \\rightarrow \\I$\n",
"3. $v(x,0) = D$ for some unknown constant $D$ on $[-1,1]$. \n",
"\n",
"### Problem 2.2\n",
"\n",
"\n",
"We write\n",
"$$\n",
"v(x,y) = \\int_{-1}^1 u(t) \\log|z-t| \\dt + \\log|z-\\I|\n",
"$$\n",
"The behaviour at infinity requires that $\\int_{-1}^1 u(t) \\dt = 0$. We further have the singular integral equation\n",
"$$\n",
" \\int_{-1}^1 u(t) \\log|x-t| \\dt = C - \\log|x-\\I|\n",
"$$\n",
"which follows from $D = v(x,0)$\n",
"\n",
"We can actually solve this numerically:"
]
},
{
"cell_type": "code",
"execution_count": 150,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.18822640645958963"
]
},
"execution_count": 150,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ũ = SingularIntegral(0) \\ (-log(abs(x-im))/π)\n",
"u = ũ - sum(ũ)/(π*sqrt(1-x^2)) # ensure integrates to zero\n",
"v = z -> π*logkernel(u, z) + log(abs(z-im))\n",
"D = v(0)"
]
},
{
"cell_type": "code",
"execution_count": 151,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 151,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"xx = yy = -4:0.011:4\n",
"V = v.(xx' .+ im*yy)\n",
"\n",
"contour(xx, yy, V; nlevels=50)\n",
"plot!(domain(x); color=:black, legend=false)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"#### Problem 2.3\n",
"\n",
"As usual for logarithmic singular integral equations we want to solve\n",
"$$\n",
"\\dashint_{-1}^1 {u(t) \\over t-x} \\dt = f'(x) \n",
"$$\n",
"we need to be a bit careful differentiating $f$:\n",
"$$\n",
"{\\D \\over \\dx} \\log|x-\\I| = {\\D \\over \\dx} \\log\\sqrt{x^2 + 1} = { x \\over x^2 + 1}\n",
"$$\n",
"Now we do our usual game and solve for $u$ using the inverse Hilbert transform formula. That is, first calculate\n",
"$$\n",
"C[\\sqrt{1-x^2} {x \\over x^2+1}](z) = {z \\sqrt{z-1} \\sqrt{z+1} \\over 2\\I (z^2+1)} - {1 \\over 2 \\I} + {\\I \\sqrt{\\I - 1} \\sqrt{\\I+1} \\over 4 (z-\\I)} + {\\I \\sqrt{-\\I - 1} \\sqrt{-\\I+1} \\over 4 (z+\\I)}\n",
"$$\n",
"Therefore, we know that\n",
"$$\n",
"u(x) = {\\sqrt{\\I - 1} \\sqrt{\\I+1} \\over 2 \\pi (x-\\I)\\sqrt{1-x^2}} + {\\sqrt{-\\I - 1} \\sqrt{-\\I+1} \\over 2 \\pi (x+\\I) \\sqrt{1-x^2}} + {C \\over \\sqrt{1-x^2}}\n",
"$$\n",
"We can show (e.g. by finding its Cauchy transform and lookin at the asymptotic behaviour) that\n",
"$$\n",
"\\int_{-1}^1\\left[ {\\sqrt{\\I - 1} \\sqrt{\\I+1} \\over 2 \\pi (x-\\I)\\sqrt{1-x^2}} + {\\sqrt{-\\I - 1} \\sqrt{-\\I+1} \\over 2 \\pi (x+\\I) \\sqrt{1-x^2}} \\right] = -1\n",
"$$\n",
"Combined with the fact that\n",
"$$\n",
"\\int_{-1}^1 {1 \\over \\sqrt{1-x^2}} = \\pi\n",
"$$\n",
"we find that $C = 1/\\pi$.\n",
"\n",
"_Verification_ Let's check our work. First we see that the Hilbert transform of $u$ does indeed satisfy the specified equation (using the numerical `u` as calculated above):"
]
},
{
"cell_type": "code",
"execution_count": 154,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"true"
]
},
"execution_count": 154,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fp = x -> x/(x^2+1)\n",
"\n",
"π*hilbert(u, 0.1) ≈ fp(0.1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And the Cauchy transform of $\\sqrt{1-x^2} f'(x)$ satisfies the derived formula:"
]
},
{
"cell_type": "code",
"execution_count": 156,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"true"
]
},
"execution_count": 156,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"C = z -> z * sqrt(z-1)sqrt(z+1)/(2im*(z^2+1)) - 1/(2im) + im*sqrt(im-1)sqrt(im+1)/(4(z-im)) +\n",
" im*sqrt(-im-1)sqrt(-im+1)/(4(z+im))\n",
"cauchy(sqrt(1-x^2)fp(x), z) ≈ C(z)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Therefore we can invert to the Hilbert transform for $f'$:"
]
},
{
"cell_type": "code",
"execution_count": 166,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"true"
]
},
"execution_count": 166,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"w = -hilbert(sqrt(1-x^2)fp(x))/sqrt(1-x^2)\n",
"w̃ = sqrt(im-1)sqrt(im+1)/(2(x-im) * sqrt(1-x^2)) +\n",
" sqrt(-im-1)sqrt(-im+1)/(2(x+im) * sqrt(1-x^2))\n",
"hilbert(w,0.1) ≈ hilbert(w̃,0.1) ≈ fp(0.1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And we have recovered $u$ up to $C/\\sqrt{1-x^2}$:"
]
},
{
"cell_type": "code",
"execution_count": 180,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"true"
]
},
"execution_count": 180,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"C = 1/π\n",
"\n",
"u(0.1) ≈ w̃(0.1)/π + C/sqrt(1-0.1^2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 3\n",
"\n",
"### Problem 3.1 \n",
"\n",
"To be analytic at all we need decay at either $\\pm \\infty$, this has neither so is not defined."
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"### Problem 3.2\n",
"\n",
"It has exponential decay in the right-half plane, therefore\n",
"$$\n",
"\\E^{\\gamma x} f(x) = {\\E^{\\gamma x } \\over 1 + \\E^x}\n",
"$$\n",
"has exponential decay at both $\\pm \\infty$, provided $0 < \\gamma < 1$. Therefore, we can take the strip $0 < \\Im s < 1$. (Note in each case the contour for the inverse Fourier transform can be any contour in the domain of analyticity.)\n",
"\n",
"We can verify this by exact computation using Residue calculus: for $0 < \\Im s < 1$, we can integrate over a rectangle to get:\n",
"$$\n",
"\\left(\\int_{-R}^R + \\int_R^{2\\I \\pi + R} + \\int_{2 \\I \\pi + R}^{2\\I \\pi - R} + \\int_{2 \\I \\pi - R}^{-R} \\right) {\\E^{-\\I s x} \\over 1 + \\E^x} \\dx = 2 \\pi \\I \\Res_{z = \\I \\pi } {\\E^{-\\I s z} \\over 1 + \\E^z} = \n",
"- 2 \\pi \\I \\E^{\\pi s}\n",
"$$\n",
"Note that \n",
"$$\n",
"{\\E^{-\\I s (R + \\I t)} \\over 1 + \\E^{R + \\I t}} = \n",
"{\\E^{-\\I R \\Re s + R \\Im s + t} \\over 1 + \\E^{R + \\I t}} \\rightarrow 0\n",
"$$\n",
"and\n",
"$$\n",
"{\\E^{-\\I s (-R + \\I t)} \\over 1 + \\E^{R + \\I t}} = \n",
"{\\E^{\\I R \\Re s - R \\Im s + t} \\over 1 + \\E^{R + \\I t}} \\rightarrow 0\n",
"$$\n",
"uniformly in $t$ as $R \\rightarrow \\infty$, hence we deduce that\n",
"$$\n",
"\\left(\\int_{-\\infty}^\\infty + \\int_{2 \\I \\pi + \\infty}^{2\\I \\pi - \\infty}\\right) {\\E^{-\\I s x} \\over 1 + \\E^x} \\dx = \n",
"- 2 \\pi \\I \\E^{\\pi s}\n",
"$$\n",
"Now note that\n",
"$$\n",
"\\int_{2 \\I \\pi + \\infty}^{2\\I \\pi - \\infty} {\\E^{-\\I s t} \\over 1 + \\E^t} \\dt = \\int_{\\infty}^{-\\infty} {\\E^{-\\I s (x+2 \\I \\pi)} \\over 1 + \\E^x} \\dx = -\\E^{2 \\pi s} \\int_{-\\infty}^\\infty {\\E^{-\\I s x} \\over 1 + \\E^x} \\dx \n",
"$$\n",
"Therefore, we have\n",
"$$\n",
"\\int_{-\\infty}^\\infty {\\E^{-\\I s x} \\over 1 + \\E^x} \\dx = - 2 \\I \\pi {\\E^{\\pi s} \\over 1 -\\E^{2 \\pi s}} = \\I \\pi {\\rm csch}\\, \\pi x\n",
"$$\n",
"which has poles at $0$ and $\\I$:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"phaseplot(-3..3, -10..10, z -> 1/(1+exp(z))) #integrand"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"phaseplot(-3..3, -3..3, z -> im*π*csch(π*z)) # transform"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Problem 3.3 \n",
"\n",
"Here $\\E^{\\gamma x } f(x) = \\E^{(\\gamma+2) x}$ has decay at $+\\infty$ proved $\\gamma < -2$, hence we have the strip $\\Im s < -2$.\n",
"\n",
"Indeed, its Fourier transform is \n",
"$$\n",
"-{\\I \\over 2 \\I +s}\n",
"$$\n",
"by integration by parts.\n",
"\n",
"### Problem 3.4 \n",
"\n",
"Here it's $\\Im s > 0$: unlike 1.1, we now have decay at $x \\rightarrow \\infty$ since $f_{\\rm L}(x)$ is identically zero.\n",
"\n",
"It's Fourier transform is determinable by integration-by-parts:\n",
"$$\n",
"\\hat f(s) = \\int_{-\\infty}^0 x \\E^{-\\I s x} \\dx = {1 \\over \\I s} \\int_{-\\infty}^0\\E^{-\\I s x} \\dx = {1 \\over s^2}\n",
"$$\n",
"\n",
"### Problem 3.5\n",
"\n",
"The Fourier transforms are given above.\n",
"\n",
"### Problem 3.6\n",
"\n",
"$$\\int_{-\\infty}^\\infty \\delta(x) \\E^{\\I s x} \\dx = 1$$\n",
"\n",
"It's actually an entire function, but non-decaying. This is hinting at the relationship between smoothness of a function and decay of its Fourier transform, and vice-versa: since $\\delta(x)$ \"decays\" to all orders, we expect its Fourier transform to be entire, but since its n ot smooth at all, we expect no decay, so on a formal level we can predict the analyticity properties."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 4\n",
"\n",
"### Problem 4.1\n",
"\n",
"Note that\n",
"$$\n",
"K(z) = {3\\over 2} \\E^{-|x|} \\Rightarrow \\hat K(s) = {3 \\over 1+s^2}\n",
"$$\n",
"Provided $-1 < \\Im s < 1$, and \n",
"$$\n",
"\\widehat f_{\\rm R}(s) = -{\\I \\over s} - {\\alpha \\over s^2}\n",
"$$\n",
"for $\\Im s < 0$. Define\n",
"$$\n",
"h(s) = - \\widehat f_{\\rm R}(s) = {\\I \\over s} + {\\alpha \\over s^2}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5.773159728050814e-15 + 7.549516567451064e-15im"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"α = 0.3\n",
"x = Fun(0..100)\n",
"f = 1 + α*x\n",
"h = s -> (im/s + α/s^2)\n",
"\n",
"γ = -0.5 # we take the Fourier transform on R + im*γ\n",
"s = -0.5 + im*γ\n",
"\n",
"-sum(f*exp(-im*s*x)) - h(s)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Transforming the equation, we have\n",
"$$\n",
"\\Phi_+(s) - (1 + \\hat K(s)) \\Phi_-(s) = {\\I \\over s} + {\\alpha \\over s^2}\n",
"$$\n",
"where\n",
"$$\n",
"1 + \\hat K(s) = {4 + s^2 \\over 1 + s^2} = {(s-2 \\I) (s+2 \\I) \\over (s+\\I)(s-\\I)}\n",
"$$\n",
"This is very close to the the example we did in lectures, so we already know the homogenous solution:\n",
"$$\n",
"\\kappa(z) = \\begin{cases} {z + 2\\I \\over z + \\I} & \\Im z > \\gamma \\\\\n",
" {z - \\I \\over z - 2\\I} & \\Im z < \\gamma\n",
" \\end{cases}\n",
" $$\n",
"which is valid for $-1 < \\Im s < 0$.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\n",
"g = s -> (4+s^2)/(1+s^2)\n",
"\n",
"κ = z -> imag(z) > γ ? (z+im*2)/(z+im) :\n",
" (z-im)/(z-im*2)\n",
"\n",
"\n",
"phaseplot(-3..3, -3..3, κ)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-1.3322676295501878e-15 - 2.7755575615628914e-16im"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s = 0.1 + γ*im \n",
"κ₊ = κ(s + eps()*im)\n",
"κ₋ = κ(s - eps()*im)\n",
"\n",
"κ₊ - κ₋*g(s)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"We thus get the RH problem\n",
"$$\n",
"Y_+(s) - Y_-(s) = h(s)/\\kappa_+(s) = ({\\I \\over s} + {\\alpha \\over s^2}) {s + \\I \\over s + 2 \\I}\n",
"$$\n",
"We see this has poles at $0$ and $-2 \\I$, so using partial fraction expansion we get\n",
"$$\n",
"({\\I \\over s} + {\\alpha \\over s^2}) {s + \\I \\over s + \\sqrt{3} \\I} = \n",
"{\\alpha \\over 2 s^2}- {\\I (\\alpha-2) \\over 4 s} +{\\I (2+\\alpha) \\over 4 (s+2 \\I)}\n",
"$$\n",
"Therefore, splitting the poles between those above and below $\\gamma$, we have\n",
"$$\n",
"Y(z) = \\begin{cases} \n",
" {\\I (2+\\alpha) \\over 4 (z+2 \\I)} & \\Im z > \\gamma \\\\\n",
"-{\\alpha \\over 2 z^2}+{\\I (\\alpha-2) \\over 4 z} & \\Im z < \\gamma\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 91,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5.551115123125783e-16 - 4.440892098500626e-16im"
]
},
"execution_count": 91,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\n",
"s = 0.1 + γ*im \n",
"Y = z -> imag(z) > γ ? im*(2+α)/(4*(z+2im)) :\n",
" - α/(2z^2) + im*(α-2)/(4z)\n",
"\n",
"\n",
"Y₊ = Y(s + eps()*im)\n",
"Y₋ = Y(s - eps()*im)\n",
"\n",
"Y₊ - Y₋ - h(s)/κ₊"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We therefore have \n",
"\n",
"$$\n",
"\\Phi(z) = \\kappa(z) Y(z) = \\begin{cases} \n",
" {\\I (2+\\alpha) \\over 4 (z+ \\I)} & \\Im z > \\gamma \\\\\n",
"(-{\\alpha \\over 2 z^2}+{\\I (\\alpha-2) \\over 4 z} ) {z - \\I \\over z - 2\\I} & \\Im z < \\gamma\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 92,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"8.881784197001252e-16 - 1.1102230246251565e-15im"
]
},
"execution_count": 92,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Φ = z -> imag(z) > γ ? im*(2+α)/(4*(z+im)) :\n",
" (-α/(2z^2) + im*(α-2)/(4z))*(z-im)/(z-2im)\n",
"Φ₊ = Φ(s+eps()im) \n",
"Φ₋ = Φ(s-eps()im) \n",
"\n",
"Φ₊ - Φ₋*g(s) - h(s)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we recover the solution by inverting \n",
"$\\Phi_-$, using Residue calculus in the upper half plane: for $x > 0$ we have\n",
"\\begin{align*}\n",
"u(x) &= {1 \\over 2 \\pi} \\int_{-\\infty+ \\I \\gamma}^{\\infty + \\I \\gamma} (-{\\alpha \\over 2 z^2}+{\\I (\\alpha-2) \\over 4 z} ) {z - \\I \\over z - 2\\I} \\E^{\\I z x} \\dz \\\\\n",
"&= \\I (\\Res_{z = 0} + \\Res_{z = 2\\I}) (-{\\alpha \\over 2 z^2}+{\\I (\\alpha-2) \\over 4 z} ) {z - \\I \\over z - 2\\I} \\E^{\\I z x} = {1+x \\alpha \\over 4} - {\\alpha+1 \\over 4} \\E^{-2 x}\n",
"\\end{align*}\n",
"\n",
"Did it work? yes:"
]
},
{
"cell_type": "code",
"execution_count": 93,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(1.0300000000000011, 1.03)"
]
},
"execution_count": 93,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = Fun(0 .. 50)\n",
"\n",
"u = (1+t*α)/4 - (α-1)/4*exp(-2t)\n",
"\n",
"x = 0.1\n",
"\n",
"u(x) + 3/2*sum(exp(-abs(t-x))*u) , f(x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Problem 4.2\n",
"\n",
"Setting up the problem as above, we arrive at a degenerate RH problem:\n",
"$$\n",
"\\Phi_+(s) - g(s) \\Phi_-(s) = h(s)\n",
"$$\n",
"where \n",
"$$g(s) = \\widehat K(s) = {2\\alpha \\over \\alpha^2 +s^2}= {2 \\alpha \\over (s-\\I \\alpha)(s+\\I \\alpha)} $$\n",
"and\n",
"$$ \n",
"h(s) = {\\I \\over s} + {\\alpha \\over s^2} = \\I {s -\\I \\alpha \\over s^2}\n",
"$$\n",
"\n",
"Suppose we allow $\\kappa_-(s) \\sim s$ to have growth, then we can write\n",
"$$\n",
"\\kappa(z) = \\begin{cases} {1 \\over z+\\I \\alpha }& \\Im z > \\gamma \\\\\n",
" {z-\\I \\alpha \\over 2 \\alpha} & \\Im z < \\gamma\n",
"\\end{cases}\n",
"$$\n",
"so that \n",
"$$\n",
"\\kappa_+(s) = \\kappa_-(s) g(s)\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\n",
"α = 0.3\n",
"\n",
"\n",
"g = s -> (2α)/(α^2+s^2)\n",
"h = s -> (im/s + α/s^2)\n",
"\n",
"κ = z -> imag(z) > γ ? 1/(z + im*α) :\n",
" (z-im*α)/(2α)\n",
"\n",
"\n",
"phaseplot(-3..3, -3..3, κ)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.6645352591003757e-15 + 1.7763568394002505e-15im"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s = 0.1 + γ*im \n",
"κ₊ = κ(s + eps()*im)\n",
"κ₋ = κ(s - eps()*im)\n",
"\n",
"κ₊ - κ₋*g(s)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then we have \n",
"$$h(s)/\\kappa_+(s) = \\I {s^2 + \\alpha^2 \\over s^2} = \\I + \\I {\\alpha^2 \\over s^2}\n",
"$$\n",
"and then we can write\n",
"$$\n",
"Y(z) = \\begin{cases}\n",
"\\I & \\Im z > \\gamma \\\\\n",
"-{\\I \\alpha^2 \\over z^2} & \\Im z < \\gamma\n",
"\\end{cases}\n",
"$$\n"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.3877787807814457e-16 + 7.771561172376096e-16im"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\n",
"s = 0.1 + γ*im \n",
"Y = z -> imag(z) > γ ? im :\n",
" -im*α^2/s^2\n",
"\n",
"\n",
"Y₊ = Y(s + eps()*im)\n",
"Y₋ = Y(s - eps()*im)\n",
"\n",
"Y₊ - Y₋ - h(s)/κ₊"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Putting things together, we get\n",
"$$\n",
"\\Phi(z) = \\kappa(z) Y(z) = \\begin{cases} {\\I \\over z+\\I \\alpha }& \\Im z > \\gamma \\\\\n",
" -\\I {\\alpha^2 \\over z^2} {z-\\I \\alpha \\over 2 \\alpha} & \\Im z < \\gamma\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-3.9968028886505635e-15 + 4.218847493575595e-15im"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Φ = z -> imag(z) > γ ? im/(z + im*α) :\n",
" -im*α^2/z^2* (z-im*α)/(2α)\n",
"Φ₊ = Φ(s+eps()im) \n",
"Φ₋ = Φ(s-eps()im) \n",
"\n",
"Φ₊ - Φ₋*g(s) - h(s)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We now invert the Fourier transform of $\\Phi_-(s)$ using Jordan's lemma:\n",
"$$\n",
"u(x) = {1 \\over 2 \\pi} \\int_{-\\infty + \\I \\gamma}^{\\infty + \\I \\gamma} \\Phi_-(s) \\E^{\\I s x} \\D s = {\\alpha \\over 2}\\Res_{z = 0} {z- \\I \\alpha \\over z^2} \\E^{\\I z x} = {\\alpha \\over 2} (1+x \\alpha)\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 122,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.199040866595169e-14"
]
},
"execution_count": 122,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = Fun(0 .. 200)\n",
"\n",
"u = α*(1+t*α)/2\n",
"\n",
"x = 0.1\n",
"\n",
"sum(exp(-α*abs(t-x))*u) - (1 + α*x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 4.3\n",
"\n",
"1. From the same logic as 2.2, we know we need to solve\n",
"$$\n",
"\\Phi_+(s) - g(s) \\Phi_-(s) = h(s)\n",
"$$\n",
"where\n",
"$$\n",
"g(s) = 1 - {2 \\lambda \\over s^2 + 1} = {s^2 +1-2 \\lambda \\over s^2 + 1} = {(s- \\I\\gamma)(s+\\I \\gamma) \\over (s+\\I)(s-\\I)}\n",
"$$\n",
"and \n",
"$$\n",
"h(s) = {1 \\over s^2}\n",
"$$\n",
"where $-1 < \\Im s < 0$, let's say $\\Im s = \\delta$ because I annoyingly used $\\gamma$ in the statement of the problem. \n",
"Writing $s = t + \\I \\delta$, we see that\n",
"$$\n",
"g(s) = {t^2 +2 \\I \\delta t -\\delta^2 +\\gamma^2 \\over s^2 + 1}\n",
"$$\n",
"By ensuring its real part is positive, this has trivial winding number provided $\\gamma^2 = 1 - 2\\lambda > 0$, which is true for $0 < \\lambda < {1 \\over 2}$, and restricting the contour $s$ lives on to be $- {\\gamma} < \\delta < 0$. Factorizing the kernel we get\n",
"$$\n",
"\\kappa(z) = \\begin{cases}\n",
" {z+\\I \\gamma \\over z + \\I} & \\Im z > \\delta\\\\\n",
" {z-\\I \\over z-\\I \\gamma} & \\Im z < \\delta \n",
" \\end{cases}\n",
"$$\n",
"Thus we want to solve\n",
"$$\n",
"Y_+(s) - Y_-(s) = h(s) \\kappa_+(s)^{-1} = {s + \\I \\over s + \\I \\gamma} { 1 \\over s^2} = {1 \\over \\gamma s^2} - {\\I ( \\gamma - 1) \\over \\gamma^2 s} + {\\I \\over \\gamma^2 }{ \\gamma - 1 \\over s+ \\I \\gamma}\n",
"$$\n",
"Which has solution, (since $\\delta > - \\gamma$), \n",
"$$\n",
"Y(z) = \\begin{cases}\n",
" {\\I \\over \\gamma^2 }{ \\gamma - 1 \\over s+ \\I \\gamma} & \\Im z > \\delta\\\\\n",
" {\\I ( \\gamma - 1) \\over \\gamma^2 z} - {1 \\over \\gamma z^2} & \\Im z < \\delta \n",
" \\end{cases}\n",
"$$\n",
"We thus get\n",
"$$\n",
"\\Phi_-(z) = ({\\I ( \\gamma - 1) \\over \\gamma^2 z} - {1 \\over \\gamma z^2}) {z-\\I \\over z-\\I \\gamma} \n",
"$$\n",
"and Jordan's lemma gives us\n",
"$$\n",
"u(x) = {x \\over \\gamma^2} - \\E^{-x \\gamma}(\\gamma-1)/\\gamma^2\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 55,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-3.239075674343894e-14"
]
},
"execution_count": 55,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = Fun(0 .. 200)\n",
"λ = 0.1\n",
"γ = sqrt(1-2λ)\n",
"\n",
"u = t/γ^2 - exp(-t*γ)*(γ-1)/γ^2\n",
"\n",
"x = 0.1\n",
"\n",
"u(x) - λ*sum(exp(-abs(t-x))*u) - x"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Oddly, this is definitely a solution, but not in the form the question asked for. To get the other solution, consider now the bad winding number case of $-1 < \\delta < - \\gamma$. Motivated by 2.2, what if we allow $\\kappa$ to have different behaviour? Consider\n",
"$$\n",
"\\kappa(z) = \\begin{cases}\n",
" {1 \\over z + \\I} & \\Im z > \\delta\\\\\n",
" {(z-\\I) \\over (z-\\I \\gamma) (z+\\I \\gamma)} & \\Im z < \\delta \n",
" \\end{cases}\n",
"$$\n",
"Chosen so that both $\\kappa_+$ and $\\kappa_+^{-1}$ are analytic. \n",
"\n",
"Thus we want to solve\n",
"$$\n",
"Y_+(s) - Y_-(s) = h(s) \\kappa_+(s)^{-1} = { s + \\I \\over s^2} = {1 \\over s} + {\\I \\over s^2}\n",
"$$\n",
"but now we only need $Y_+(s) = O(1)$ and $Y_-(s) = O(1)$. Here is where the non-uniqueness comes in, as we can add an arbitrary constant:\n",
"$$\n",
"Y(z) = \\begin{cases}\n",
" A & \\Im z > 0 \\\\\n",
" A - {1 \\over z} - {\\I \\over z^2} & \\Im z < 0\n",
"\\end{cases}\n",
"$$\n",
"Thus we have\n",
"$$\n",
"\\Phi_-(z) = Y_-(z) \\kappa_-(z) = -( A + {1 \\over z} + {\\I \\over z^2} ){(z-\\I) \\over (z-\\I \\gamma) (z+\\I \\gamma)} \n",
"$$\n",
"Using Jordan's lemma, and now since $\\delta < - \\gamma$, we get\n",
"\\begin{align*}\n",
"u(x) &= \\I (\\Res_{z = 0} + \\Res_{z = \\I \\gamma} + \\Res_{z = - \\I \\gamma}) \\Phi_-(z) \\E^{\\I x z} \\\\\n",
"&= {x \\over \\gamma^2} - \\E^{-x \\gamma}( {\\gamma^2-1 \\over 2 \\gamma^3} + {\\gamma-1 \\over 2\\gamma^3} A) - \\E^{x \\gamma} ({1-\\gamma^2 \\over 2 \\gamma^3} + {\\gamma+1 \\over 2\\gamma^3} A)\\\\\n",
"&= {x \\over \\gamma^2} + {\\E^{x \\gamma} - \\E^{-x \\gamma} \\over 2} {\\gamma-\\gamma^{-1} \\over 2 \\gamma^2} - {A \\over \\gamma^3} ( {\\E^{x \\gamma} - \\E^{-x \\gamma} \\over 2} + \\gamma {\\E^{x \\gamma} + \\E^{-x \\gamma} \\over 2} )\n",
"\\end{align*}\n",
"Redefining $A$ and using the definition of $\\sinh$ and $\\cosh$ gives the form in the assignment.\n",
"\n",
"What's the moral of the story? \n",
"\n",
"1. Different choices of contours can give different solutions\n",
"2. When the winding number is non-trivial, the solution may not be unique"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"### 4.4\n",
"\n",
"\n",
"1. Integrating by parts we have\n",
"\\begin{align*}\n",
"\\widehat{u_\\rR'}(s) &= \\I s \\widehat{u_\\rR}(s) - u(0) = \\I s \\widehat{u_\\rR}(s) \\\\\n",
"\\widehat{u_\\rR''}(s) &= \\I s \\widehat{u_\\rR'}(s) - u'(0) = -s^2 \\widehat{u_\\rR}(s) - u'(0)\n",
"\\end{align*}\n",
"2. Our integral equation when cast on the whole real line is:\n",
"$$\n",
"u_\\rR''(x) - {72 \\over 5} \\int_{-\\infty}^\\infty \\E^{-5 |x-t|} u_\\rR(t)\\dt = 1_\\rR(x) + p_\\rL(x)\n",
"$$\n",
"where \n",
"$$\n",
"p(x) = {72 \\over 5} \\int_{-\\infty}^\\infty \\E^{-5 |x-t|} u_\\rR(t)\\dt = {72 \\over 5} \\int_0^\\infty \\E^{-5 |x-t|} u_\\rR(t)\\dt.\n",
"$$\n",
"Note that, for $-5 <\\Im s < 5$,\n",
"$$\n",
"\\hat K(s) = {10 \\over s^2 + 25}\n",
"$$\n",
"provided $s$ is in the lower half plane,\n",
"$$\n",
"\\widehat{1_\\rR}(s) = \\int_0^\\infty \\E^{-\\I s x} \\dx = {1 \\over \\I s}\n",
"$$\n",
"Thus our integral equation in frequency space is \n",
"\\begin{align*}\n",
"-\\alpha - s^2 \\widehat{u_\\rR}(s) - {72 \\over 5} \\widehat K(s) \\widehat{u_\\rR}(s) &=\\widehat{p_\\rL}(x) + \\widehat{1_\\rR}(s)\\\\\n",
"\\Phi_+(s) - (s^2 + {144 \\over s^2+25}) \\Phi_-(s) &= \\alpha + {1 \\over \\I s}\\\\\n",
"\\Phi_+(s) - {(s^2 + 9) (s^2+16) \\over s^2+25} \\Phi_-(s) &= \\alpha + {1 \\over \\I s}\n",
"\\end{align*}\n",
"where $s \\in \\R +\\I \\gamma$ for any $-5 < \\gamma < 0$.\n",
"3. We can factorize this to construct $g(s)$ as \n",
"$$\n",
"g(s) = \\kappa_+(s) \\kappa_-(s)^{-1} = {(s +3 \\I) (s+4 \\I) \\over s+5 \\I} {(s -3 \\I) (s-4 \\I) \\over s-5 \\I}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-5.551115123125783e-17 - 2.220446049250313e-16im"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"κ = z -> imag(z) > γ ? \n",
" (z+3im)*(z+4im)/(z+5im) :\n",
" (z-5im)/((z-3im)*(z-4im)) \n",
"\n",
"γ = -1.0\n",
"s = 0.1+γ*im\n",
"g = s -> (s^2+9)*(s^2+16)/(s^2+25)\n",
"\n",
"κ(s+eps()im) - g(s)κ(s-eps()im)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Writing $\\Phi(z) = \\kappa(z) Y(z)$ we get the subtractive RH problem\n",
"$$\n",
"Y_+(s) - Y_-(s) = {h(s) \\over \\kappa_+(s)} = (\\alpha + {1 \\over \\I s}) {s + 5 \\I \\over (s + 3 \\I)(s + 4 \\I)}\n",
"$$\n",
"We use partial fraction expansion to write\n",
"$$\n",
"{h(s) \\over \\kappa_+(s)} = -{\\alpha + 1/4 \\over s+4 \\I} + {2/3 + 2 \\alpha \\over s+3 \\I} - {5 \\over 12 s}\n",
"$$\n",
"Therefore we have\n",
"$$\n",
"Y(z) = \\begin{cases} \n",
"-{\\alpha + 1/4 \\over s+4 \\I} + {2/3 + 2 \\alpha \\over s+3 \\I} & 2 \\\\\n",
" {5 \\over 12 s} &1\n",
"\\end{cases}\n",
"$$\n",
"and hence\n",
"$$\n",
"\\Phi(z) = \\begin{cases}\n",
"{(z+3\\I)(z+4\\I) \\over z+5\\I} (-{α+1/4 \\over z+4\\I} + (2/3 + 2α)/(z+3\\I)) & \\Im z > \\gamma \\\\\n",
" {z-5\\I \\over (z-3\\I)(z-4\\I)} {5 \\over 12z}& \\Im z < \\gamma\n",
" \\end{cases}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can now invert the Fourier transform of \n",
"$$\n",
"\\Phi_-(s) = {s-5\\I \\over (s-3\\I)(s-4\\I)} {5 \\over 12s}\n",
"$$\n",
"This actually decays so fast that we don't need Jordan's lemma to justify here. This has three poles above our contour, so we sum over each residue to get\n",
"$$\n",
"u(x) = \\I (\\Res_{z = 0} +\\Res_{z = 3 \\I } +\\Res_{z = 4\\I} ) \\E^{\\I z x } {z-5\\I \\over (z-3\\I)(z-4\\I)} {5 \\over 12z} = -{25 \\over 144} - {5 \\E^{-4 x} \\over 48} + {5 \\E^{-3 x} \\over 18}\n",
"$$\n",
"\n",
"Here's we check the solution:"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.0000000000000018"
]
},
"execution_count": 42,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = Fun(0 .. 200)\n",
"\n",
"u = -25/144 - 5exp(-4t)/48 + 5exp(-3t)/18\n",
"\n",
"x = 1.1\n",
"u''(x) - 72/5*sum(exp(-5abs(x-t))*u)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here we check the jump of $Y$:"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.7755575615628914e-17 + 5.551115123125783e-17im"
]
},
"execution_count": 45,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"α = u'(0)\n",
"\n",
"h = s -> α + 1/(im*s)\n",
"\n",
"Y = z -> imag(z) > γ ? \n",
" (-(α+1/4)/(z+4im) + (2/3 + 2α)/(z+3im)) :\n",
" 5/(12z)\n",
"\n",
"Y(s+eps()im) - Y(s-eps()im) - h(s)/κ(s + eps()im)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here we check the jump of $\\Phi$:"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-1.1102230246251565e-16 + 1.3877787807814457e-17im"
]
},
"execution_count": 46,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"γ = -1.0\n",
"\n",
"\n",
"Φ = z -> imag(z) > γ ? \n",
" (z+3im)*(z+4im)/(z+5im) * (-(α+1/4)/(z+4im) + (2/3 + 2α)/(z+3im)) :\n",
" (z-5im)/((z-3im)*(z-4im)) * 5/(12z)\n",
"\n",
"\n",
"Φ(s + eps()*im) - g(s)*Φ(s - eps()*im) - h(s)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.0.0",
"language": "julia",
"name": "julia-1.0"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.0.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}