{
"cells": [
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"using Plots, ComplexPhasePortrait, ApproxFun, SingularIntegralEquations,SpecialFunctions, OscillatoryIntegrals\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\\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\\questionequals{= \\!\\!\\!\\!\\!\\!{\\scriptstyle ? \\atop }\\,\\,\\,}\n",
"$$\n",
"\n",
"Dr Sheehan Olver\n",
"
\n",
"s.olver@imperial.ac.uk\n",
"\n",
"
\n",
"Website: https://github.com/dlfivefifty/M3M6LectureNotes\n",
"\n",
"\n",
"\n",
"# Lecture 25: The Wiener–Hopf method\n",
"\n",
"We can now employ the Wiener–Hopf method to solve\n",
"$$\n",
"\\lambda u(x) + \\int_{0}^\\infty K(x-t)u(t) \\dt = f(x)\\qqfor 0 < x < \\infty.\n",
"$$\n",
"For simplicity, take $\\lambda = 1$ (note we can divide through). \n",
"\n",
"The Wiener–Hopf method consists of the following steps:\n",
"\n",
"1. Calculate Fourier transforms $\\widehat K(s)$ and $\\widehat{f_{\\rm R}}(s)$ to reduce to RH problem $\\Phi_+(s) -g(s) \\Phi_-(s) = h(s)$ with $\\lim \\Phi(z) = 0$\n",
"2. Find homogenous solution $\\kappa_+(s) = g(s) \\kappa_-(s)$ with $\\lim \\kappa(z) = 1$\n",
"3. Solve Cauchy transform problem $Y_+(s) - Y_-(s) = {h(s) \\over \\kappa_+(s)}$\n",
"4. Find $u(x)$ by inverting Fourier transform $\\widehat{u_{\\rm R}}(s) = \\Phi_-(s) = \\kappa_-(s) Y_-(s)$\n",
"\n",
"We will demonstrate these four steps using the example $K(x) = \\E^{-|x|}$ and $f(x) = \\E^{-x}$. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Step 1: Calculate Fourier transforms\n",
"\n",
"We have \n",
"$$\n",
"\\widehat{K}(s) = \\int_{-\\infty}^0 \\E^{t}\\E^{-\\I s t} \\dt + \\int_0^\\infty \\E^{-t} \\E^{-\\I s t} \\dt = {1 \\over 1-\\I s} - {1 \\over -1-\\I s} = {2 \\over 1+s^2}\n",
"$$\n",
"Because we have exponential decay at $\\pm \\infty$, $\\widehat{K}$ is analytic in a strip.\n",
"Furthermore, \n",
"$$\n",
"\\widehat{f_{\\rm R}}(s) = \\int_0^\\infty \\E^{-t} \\E^{-\\I t s} \\dt = -{1 \\over -1-\\I s} = {\\I \\over \\I-s}\n",
"$$\n",
"As predicted, $\\widehat{\\fR}(s)$ is analytic in the lower half-plane.\n",
"\n",
"We thus have the RH problem:\n",
"$$\n",
"\\underbrace{\\Phi_+(s)}_{\\widehat{p_{\\rm L}}(s)} - \\underbrace{g(s)}_{1 + \\widehat K(s)} \\underbrace{\\Phi_-(s)}_{\\widehat{u_{\\rm R}}(s)} = \\underbrace{h(s)}_{-\\widehat{\\fR}(s)} \\qqand \\lim\\Phi(z) = 0\n",
"$$\n",
"where \n",
"$$\n",
"p(x) = \\int_{-\\infty}^\\infty K(x-t)u_{\\rm R}(t) \\dt\n",
"$$\n",
"is unknown, or in our case\n",
"$$\n",
"\\Phi_+(s) - {3 + s^2 \\over 1 + s^2} \\Phi_-(s) = {\\I \\over s - \\I} \\qqand \\lim\\Phi(z) = 0\n",
"$$\n",
"\n",
"Here we confirm our calculations are correct for the Fourier transforms:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.19999999999999996 - 0.4000000000000001im, 0.2 - 0.4im)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = Fun(0.0 .. 40)\n",
"f = exp(-t)\n",
"fourier(f, -2.0), im/(im-2.0)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(1.4 - 4.996003610813204e-16im, 1.4)"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"1-(fourier(-exp(-t), -2.0) + fourier(-exp(Fun(-40 .. 0)), -2.0)), (2.0^2 + 3)/(2.0^2+1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Step 2: Find homogeneous solution\n",
"\n",
"Last lecture we employed the \"guess and check\" kernel factorization method to find\n",
"$$\n",
"g(s) = {3 + s^2 \\over 1 + s^2} = \\underbrace{s + \\I \\sqrt{3} \\over s + \\I}_{\\kappa_+(s)} \\underbrace{s - \\I \\sqrt{3} \\over s-\\I }_{\\kappa_-(s)^{-1}}\n",
"$$\n",
"That is, we have the solution\n",
"$$\n",
"\\kappa(z) = \\begin{cases} {z + \\I \\sqrt{3} \\over z + \\I} & \\Im z > 0 \\\\\n",
" {z - \\I \\over z - \\I \\sqrt{3}} & \\Im z < 0\n",
" \\end{cases}\n",
" $$\n",
"Always verify this solution satisfies the right conditions:\n",
"\n",
"1. $\\kappa(z)$ is analytic off $\\R$, since the pole of $ {z + \\I \\sqrt{3} \\over z + \\I}$ is in the lower-half plane and the pole of $ {z - \\I \\over z - \\I \\sqrt{3}} $ is in the upper-half plane\n",
"2. $\\lim \\kappa(z) = 1$ by L'Hopital's rule or similar.\n",
"3. It has the right jump:\n",
"$$\n",
"\\kappa_+(s) = {s + \\I \\sqrt{3} \\over s + \\I} = {3 + s^2 \\over 1 + s^2}{s - \\I \\over s - \\I \\sqrt{3}} = g(s)\\kappa_-(s) \n",
"$$\n",
"\n",
"Here we confirm that we have the right jump:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-2.220446049250313e-16 - 4.163336342344337e-17im"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"κ = z -> imag(z) > 0 ? (z+im*sqrt(3))/(z+im) :\n",
" (z-im)/(z-im*sqrt(3))\n",
"\n",
"κ(0.1+eps()im) - κ(0.1-eps()im)*(3+0.1^2)/(1+0.1^2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## Step 3: Solve Cauchy transform problem\n",
"\n",
"Writing $\\Phi(z) = \\kappa(z) Y(z)$ we have\n",
"$$\n",
"\\Phi_+(s) - g(s) \\Phi_-(s) = \\kappa_+(s) (Y_+(s) - Y_-(s))\n",
"$$\n",
"hence we need to solve\n",
"$$\n",
"Y_+(s) - Y_-(s) = {h(s) \\over \\kappa_+(s)} = {\\I \\over s- \\I} {s + \\I \\over s+ \\I \\sqrt{3}}\n",
"$$\n",
"As in last lecture, we found:\n",
"$$\n",
"Y(z) = \\begin{cases} \n",
"{-\\I \\over 1+\\sqrt 3} {1-\\sqrt 3 \\over z+ \\I \\sqrt{3}} & \\Im z > 0 \\\\\n",
"{-2 \\I \\over z-\\I}{1 \\over 1+ \\sqrt{3}} & \\Im z < 0\n",
"\\end{cases}\n",
"$$\n",
"Let's double check that $Y$ has the right jump:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-0.5706161966689216 + 0.08138224450306746im"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Y = z -> imag(z) > 0 ? -im*(1-sqrt(3))/(1+sqrt(3))/(z+im*sqrt(3)) :\n",
" -2im/((z-im)*(1+sqrt(3)))\n",
"\n",
"s = 0.1\n",
"\n",
"Y(s+eps()im) - Y(s-eps()im)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-0.5706161966689216 + 0.08138224450306748im"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"im/(s-im)*(s+im)/(s + im*sqrt(3))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Thus we get\n",
"$$\n",
"\\begin{align*}\n",
"\\Phi(z) &= \\kappa(z) \\CC_\\R\\br[{h \\over \\kappa_+}](z) = \\begin{cases} \n",
"{-\\I \\over 1+\\sqrt 3} {1-\\sqrt 3 \\over z+\\I} & \\Im z > 0 \\\\\n",
"{-2 \\I \\over 1+ \\sqrt{3}} {1 \\over z-\\sqrt{3}\\I} & \\Im z < 0\n",
"\\end{cases}\n",
"\\end{align*}\n",
"$$\n",
"\n",
"We can confirm it satisfies the right RHP:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.1102230246251565e-16 - 2.7755575615628914e-17im"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"φ = z -> imag(z) > 0 ? -im*(1-sqrt(3))/(1+sqrt(3))/(z+im) :\n",
" -2im/((z-sqrt(3)*im)*(1+sqrt(3)))\n",
"\n",
"g = s -> (s^2+3)/(s^2+1)\n",
"\n",
"s = 0.1\n",
"φ(s+eps()im) - φ(s-eps()im)*g(s) - im/(s-im)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Step 4: Invert the Fourier transform\n",
"\n",
"In particular, we have\n",
"$$\n",
"\\widehat{u_{\\rm R}}(s) = \\Phi_-(s) = \\kappa_-(s) Y_-(s) = {-2 \\I \\over (1+\\sqrt{3})(s-\\I \\sqrt 3)}\n",
"$$\n",
"This is analytic in the lower half plane, and for $x > 0$ we can use Jordan's lemma hence use Residue calculus in the upper-half plane:\n",
"$$\n",
"u(x) = {1 \\over 2 \\pi} \\int_{-\\infty}^\\infty \\widehat{u_{\\rm R}}(s) \\E^{\\I s x} \\ds = {-\\I \\over \\pi(1+\\sqrt{3})} \\int_{-\\infty}^\\infty {1 \\over s-\\I \\sqrt 3} \\E^{\\I s x} \\ds = {2 \\over (1+\\sqrt{3})} \\Res_{z=\\I \\sqrt 3} {1 \\over z-\\I \\sqrt 3} \\E^{\\I z x} = {2 \\E^{-\\sqrt{3}x} \\over 1+\\sqrt{3}} \n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## Step 5: Check I didn't make a mistake\n",
"\n",
"Did this work? Amazingly, yes (and only took me four tries to get the sign right! 🕺):\n",
"\\begin{align*}\n",
"u(x) + \\int_0^\\infty K(t -x) u(t) \\dt &= {2 \\over 1+\\sqrt{3}} \\left( \\E^{-\\sqrt{3}x} + \\int_0^x \\E^{t-x} \\E^{-\\sqrt{3}t} \\dt + \\int_x^\\infty \\E^{x-t} \\E^{-\\sqrt{3}t}\\dt \\right) \\\\\n",
"&={2 \\over 1+\\sqrt{3}}\\left( \\E^{-\\sqrt{3}x} + \\E^{-x} \\int_0^x \\E^{(1-\\sqrt{3})t} \\dt + \\E^x \\int_x^\\infty \\E^{-(1+\\sqrt{3})t}\\dt \\right) \\\\\n",
"&={2 \\over 1+\\sqrt{3}}\\left( \\E^{-\\sqrt{3}x} + {1 \\over 1-\\sqrt 3} (\\E^{-\\sqrt{3}x} - \\E^{-x}) + { \\E^{-\\sqrt 3 x } \\over 1+\\sqrt{3}} \\right) \\\\\n",
"&=\\E^{-x} + {2 \\E^{-\\sqrt{3}x} \\over 1+\\sqrt{3}}\\left({-2 \\over -2} + {1+ \\sqrt 3 \\over -2} + {1 - \\sqrt 3 \\over -2} \\right) \\\\\n",
"&= \\E^{-x}\n",
"\\end{align*}\n",
"\n",
"We can verify it numerically:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-4.036770917537069e-13"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = Fun(0 .. 10)\n",
"\n",
"u = 2exp(-sqrt(3)*t)/(1+sqrt(3))\n",
"x = 0.1\n",
"\n",
"u(x) + sum(exp(-abs(t-x))*u) - exp(-x)"
]
}
],
"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
}