{
 "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",
    "<br>\n",
    "s.olver@imperial.ac.uk\n",
    "\n",
    "<br>\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
}