Solution Sheet 3.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "┌ Info: Recompiling stale cache file /Users/sheehanolver/.julia/compiled/v1.1/Plots/ld3vC.ji for Plots [91a5bcdd-55d7-5caf-9e0b-520d859cae80]\n",
      "└ @ Base loading.jl:1184\n",
      "┌ Info: Recompiling stale cache file /Users/sheehanolver/.julia/compiled/v1.1/ComplexPhasePortrait/tkJfA.ji for ComplexPhasePortrait [38ac1a67-8e16-5889-9a62-b2c9995eb50f]\n",
      "└ @ Base loading.jl:1184\n",
      "┌ Info: Recompiling stale cache file /Users/sheehanolver/.julia/compiled/v1.1/ApproxFun/jGqLz.ji for ApproxFun [28f2ccd6-bb30-5033-b560-165f7b14dc2f]\n",
      "└ @ Base loading.jl:1184\n",
      "┌ Info: Recompiling stale cache file /Users/sheehanolver/.julia/compiled/v1.1/SingularIntegralEquations/OCv8s.ji for SingularIntegralEquations [e094c991-5a90-5477-8896-c1e4c9552a1a]\n",
      "└ @ Base loading.jl:1184\n",
      "┌ Warning: Package SingularIntegralEquations does not have Random in its dependencies:\n",
      "│ - If you have SingularIntegralEquations checked out for development and have\n",
      "│   added Random as a dependency but haven't updated your primary\n",
      "│   environment's manifest file, try `Pkg.resolve()`.\n",
      "│ - Otherwise you may need to report an issue with SingularIntegralEquations\n",
      "└ Loading Random into SingularIntegralEquations from project dependency, future warnings for SingularIntegralEquations are suppressed.\n"
     ]
    }
   ],
   "source": [
    "using Plots, ComplexPhasePortrait, ApproxFun, SingularIntegralEquations\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\\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",
    "<br>\n",
    "s.olver@imperial.ac.uk\n",
    "\n",
    "<br>\n",
    "Website: https://github.com/dlfivefifty/M3M6LectureNotes\n",
    "\n",
    "# Solution Sheet 3\n",
    "\n",
    "\n",
    "## Problem 1.1 \n",
    "\n",
    "We actually start by showing the second properties of Problem 1.2, for all $\\alpha$:\n",
    "$$\n",
    "{\\D \\over \\dx} \\br[x^\\alpha \\E^{-x} L_n^{(\\alpha)}(x)] = {1 \\over n!}  {\\D^{n+1} \\over \\dx^{n+1}}\\br[x^{\\alpha+n}\\E^{-x}] = (n+1)x^{\\alpha-1}\\E^{-x}  {x^{1-\\alpha} \\E^{x} \\over (n+1)!}{ \\D^{n+1} \\over \\dx^{n+1}}\\br[x^{\\alpha+n}\\E^{-x}] = (n+1) x^{\\alpha-1}\\E^{-x} L_{n+1}^{(\\alpha-1)}(x).\n",
    "$$\n",
    "Expanding out the derivative we see\n",
    "$$\n",
    "x^{\\alpha-1} \\E^{-x}\\pr({(\\alpha -x)L_n^{(\\alpha)}(x) + (L_n^{(\\alpha)})'(x)}) = (n+1) x^{\\alpha-1}\\E^{-x} L_{n+1}^{(\\alpha-1)}(x)\n",
    "$$\n",
    "or in other words\n",
    "$$\n",
    "(\\alpha -x)L_n^{(\\alpha)}(x) + (L_n^{(\\alpha)})'(x) = (n+1) L_{n+1}^{(\\alpha-1)}(x)\n",
    "$$\n",
    "By induction with the fact $L_0^{(\\alpha)}(x) = 0$, we therefore get \n",
    "$$L_n^{(\\alpha)}(x) = {(\\alpha+1 -x)L_{n-1}^{(\\alpha+1)}(x) + (L_n^{(\\alpha+1)})'(x) \\over n}\n",
    "$$\n",
    "is a degree $n$ polynomial. We further have that the leading coefficient is\n",
    "\\begin{align*}\n",
    "L_n^{(\\alpha)}(x) &= -{x \\over n} L_{n-1}^{(\\alpha+1)}(x) +O(x^{n-1}) =   {x^2 \\over n(n-1)} L_{n-2}^{(\\alpha+2)}(x) +O(x^{n-1}) = \\cdots = {(-1)^n x^n \\over n!} L_0^{(\\alpha+n)}(x) +O(x^{n-1}) \\\\\n",
    "&= {(-1)^n x^n \\over n!} +O(x^{n-1})\n",
    "\\end{align*}\n",
    "\n",
    "We now show orthogonality with lower degree polynomials using integration by parts:\n",
    "$$\n",
    "\\int_0^\\infty L_n^{(\\alpha)}(x) p_m(x) x^\\alpha \\E^{-x} \\dx = \n",
    "\\int_0^\\infty  {1 \\over n!} {\\D^n \\over \\dx^n}\\br[x^{\\alpha+n}\\E^{-x}]  p_m(x)  \\dx = (-1)^n \\int_0^\\infty  {1 \\over n!} \\br[x^{\\alpha+n}\\E^{-x}]  p_m^{(n)}(x)  \\dx = 0\n",
    "$$\n",
    "since $p_m^{(n)}(x) = 0$.\n",
    "Note we use the fact that\n",
    "$$\n",
    "{\\D^k \\over \\dx^k}\\br[x^{\\alpha+n}\\E^{-x}]  = O(x^{\\alpha+n-k})\n",
    "$$\n",
    "hence vanishes at zero to ignore the boundary terms in integration by parts.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "\n",
    "## Problem 1.2 \n",
    "\n",
    "We showed the second property as part of 1.1. For the first part, it is clear that we have the correct constant. Now we show orthogonality with all degree $m < n-1$ polynomials (using the fact that $x^{\\alpha+1} \\E^{-x}$ is zero at $x = 0$):\n",
    "$$\n",
    "\\int_0^\\infty {\\D L_n^{(\\alpha)}(x) \\over \\dx} p_m(x) x^{\\alpha+1} \\E^{-x} \\dx = -\\int_0^\\infty L_n^{(\\alpha)}(x) (x p_m'(x) + (\\alpha+1) p_m - x p_m) x^{\\alpha} \\E^{-x} \\dx = 0\n",
    "$$\n",
    "since $(x p_m'(x) + (\\alpha+1) p_m - x p_m) $ is degree $m+1 < n$.\n",
    "\n",
    "For the third part, use the product rule on the last derivative:\n",
    "\\begin{align*}\n",
    "(n+1) L_{n+1}^{(\\alpha)}(x) &= {x^{-\\alpha}\\E^x \\over n!} {\\D^{n} \\over \\dx^{n}} {\\D \\over \\dx} \\br[x^{\\alpha+n+1} \\E^{-x}] \\\\\n",
    "  &= {x^{-\\alpha}\\E^x \\over n!} {\\D^{n} \\over \\dx^{n}}  \\br[(\\alpha+n+1)x^{\\alpha+n} \\E^{-x}-x^{\\alpha+n+1} \\E^{-x}]  \\\\\n",
    "    &= (\\alpha+n+1)L_n^{(\\alpha)}(x) -   xL_n^{(\\alpha+1)}(x)\n",
    "\\end{align*}\n",
    "For the last result, we apply the product rule $n$ times:\n",
    "\\begin{align*}\n",
    "L_{n}^{(\\alpha+1)}(x) &= {x^{-1-\\alpha}\\E^x \\over n!} {\\D^{n} \\over \\dx^{n}} {\\D \\over \\dx} \\br[x x^{\\alpha+n} \\E^{-x}] \\\\\n",
    "&= {x^{-1-\\alpha}\\E^x \\over n!} {\\D^{n-1} \\over \\dx^{n-1}} \\br[ x^{\\alpha+n} \\E^{-x}]  + {x^{-1-\\alpha}\\E^x \\over n!} {\\D^{n-1} \\over \\dx^{n-1}} x {\\D \\over \\dx} \\br[ x^{\\alpha+n} \\E^{-x}]  \\\\\n",
    "&= {2 \\over n}L_{n-1}^{(\\alpha+1)}(x)   + {x^{-1-\\alpha}\\E^x \\over n!} {\\D^{n-2} \\over \\dx^{n-2}} x {\\D^2 \\over \\dx^2} \\br[ x^{\\alpha+n} \\E^{-x}]  \\\\\n",
    "&= {3 \\over n}L_{n-1}^{(\\alpha+1)}(x)   + {x^{-1-\\alpha}\\E^x \\over n!} {\\D^{n-3} \\over \\dx^{n-3}} x {\\D^3 \\over \\dx^3} \\br[ x^{\\alpha+n} \\E^{-x}]  \\\\\n",
    "&\\vdots\\\\\n",
    "&={n \\over n}L_{n-1}^{(\\alpha+1)}(x)   + {x^{-\\alpha}\\E^x \\over n!} {\\D^{n} \\over \\dx^{n}}  \\br[ x^{\\alpha+n} \\E^{-x}] \\\\\n",
    "&=L_{n-1}^{(\\alpha+1)}(x) + L_n^{(\\alpha)}(x)\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "\n",
    "## Problem 1.3 \n",
    "\n",
    "Note that relationship 3 above did not depend on $\\alpha >-1$.  We therefore have from 1.2, comibing property (3) and (4),\n",
    "\\begin{align*}\n",
    "x L_n^{(\\alpha)}(x) &= -(n+1)L_{n+1}^{(\\alpha-1)}(x) +(n+\\alpha)L_n^{(\\alpha-1)}(x) \\\\\n",
    "  &= -(n+1)L_{n+1}^{(\\alpha)}(x) + (n+1)L_{n}^{(\\alpha)}(x) +(n+\\alpha)L_{n}^{(\\alpha)}(x) - (n+\\alpha)L_{n-1}^{(\\alpha)}(x) \\\\\n",
    "  &= - (n+\\alpha)L_{n-1}^{(\\alpha)}(x) + (2n+\\alpha+1) L_n^{(\\alpha)}(x) -(n+1)L_{n+1}^{(\\alpha)}(x)\n",
    "\\end{align*}\n",
    "\n",
    "The Jacobi operator therefore has the form\n",
    "$$\n",
    "x \\begin{pmatrix} L_0^{(\\alpha)}(x) \\\\\n",
    "            L_1^{(\\alpha)}(x) \\\\\n",
    "                \\vdots\n",
    "                \\end{pmatrix} = \\begin{pmatrix} \\alpha+1 &-1 \\\\\n",
    "                                                -1-\\alpha & \\alpha+3 & -2 \\\\\n",
    "                                                   & -2-\\alpha & \\alpha+5 &-3 \\\\                  &       & -3-\\alpha & \\alpha+7 &-4 \\\\\n",
    "                                                      &&   & -4-\\alpha & \\alpha+9 &\\ddots \\\\\n",
    "                                                      &&&&\\ddots & \\ddots\n",
    "                                                    \\end{pmatrix} \\begin{pmatrix} L_0^{(\\alpha)}(x) \\\\\n",
    "            L_1^{(\\alpha)}(x) \\\\\n",
    "                \\vdots\n",
    "                \\end{pmatrix}\n",
    "$$\n",
    "\n",
    "_Demonstration_ The following command creates the Jacobi operators for a Laguerre polynomial with `α = 1`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "TransposeOperator : Laguerre{Int64,Ray{false,Float64}}(1, 【0.0,∞❫) → Laguerre{Int64,Ray{false,Float64}}(1, 【0.0,∞❫)\n",
       "  2   -1     ⋅    ⋅    ⋅    ⋅    ⋅    ⋅     ⋅    ⋅  ⋅\n",
       " -2    4   -2     ⋅    ⋅    ⋅    ⋅    ⋅     ⋅    ⋅  ⋅\n",
       "   ⋅  -3    6   -3     ⋅    ⋅    ⋅    ⋅     ⋅    ⋅  ⋅\n",
       "   ⋅    ⋅  -4    8   -4     ⋅    ⋅    ⋅     ⋅    ⋅  ⋅\n",
       "   ⋅    ⋅    ⋅  -5   10   -5     ⋅    ⋅     ⋅    ⋅  ⋅\n",
       "   ⋅    ⋅    ⋅    ⋅  -6   12   -6     ⋅     ⋅    ⋅  ⋅\n",
       "   ⋅    ⋅    ⋅    ⋅    ⋅  -7   14   -7      ⋅    ⋅  ⋅\n",
       "   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅  -8   16    -8     ⋅  ⋅\n",
       "   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅  -9    18   -9   ⋅\n",
       "   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅  -10   20   ⋱\n",
       "   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅     ⋅    ⋱  ⋱"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ApproxFun.Recurrence(Laguerre(1))'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 2.1\n",
    "\n",
    "Note that\n",
    "$$\n",
    "{\\D \\over \\dx} \\E^{-x/2} u(x) =  \\E^{-x/2} (-{u(x) \\over 2} +u'(x)) \n",
    "$$\n",
    "Thus \n",
    "$$\n",
    "{\\D \\over \\dx} \\E^{-x/2} u(x)  = \\E^{-x/2} (u'(x) -{u(x) \\over 2} - x u(x))\n",
    "$$\n",
    "\n",
    "\n",
    "We have the derivative operator from $L_k(x)$ to $L_k^{(1)}(x)$ as:\n",
    "$$\n",
    "D = \\begin{pmatrix}\n",
    "0 & -1 \\\\\n",
    "  &&-1 \\\\\n",
    "  &&&\\ddots\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "and the Multiplication operator for $\\alpha = 0$ (from Problem 1.3)\n",
    "$$\n",
    "J^\\top = \\begin{pmatrix} 1 &-1\\\\\n",
    "                        -1 & 3 &-2\\\\\n",
    "                            &\\ddots & \\ddots & \\ddots\n",
    "                            \\end{pmatrix}\n",
    "$$\n",
    "and the conversion operator (from Problem 1.2 property 4)\n",
    "$$\n",
    "S = \\begin{pmatrix}\n",
    "        1 & -1 \\\\ & 1 & -1 \\\\&&\\ddots & \\ddots\n",
    "\\end{pmatrix}        \n",
    "$$\n",
    "We thus have multiplication by $x$ from basis to the other as\n",
    "$$\n",
    "S J^\\top = \\begin{pmatrix}\n",
    "2 & -4 & 2 \\\\\n",
    "-1 & 5 & -7 & 3 \\\\\n",
    "& -2 & 8 & -10 & 4\\\\\n",
    "&&\\ddots&\\ddots&\\ddots&\\ddots\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "Putting everything together, we get the operator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Recurrence : Laguerre{Int64,Ray{false,Float64}}(0, 【0.0,∞❫) → Laguerre{Int64,Ray{false,Float64}}(0, 【0.0,∞❫)\n",
       "  1   -1     ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅  ⋅\n",
       " -1    3   -2     ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅  ⋅\n",
       "   ⋅  -2    5   -3     ⋅    ⋅    ⋅    ⋅    ⋅    ⋅  ⋅\n",
       "   ⋅    ⋅  -3    7   -4     ⋅    ⋅    ⋅    ⋅    ⋅  ⋅\n",
       "   ⋅    ⋅    ⋅  -4    9   -5     ⋅    ⋅    ⋅    ⋅  ⋅\n",
       "   ⋅    ⋅    ⋅    ⋅  -5   11   -6     ⋅    ⋅    ⋅  ⋅\n",
       "   ⋅    ⋅    ⋅    ⋅    ⋅  -6   13   -7     ⋅    ⋅  ⋅\n",
       "   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅  -7   15   -8     ⋅  ⋅\n",
       "   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅  -8   17   -9   ⋅\n",
       "   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅  -9   19   ⋱\n",
       "   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋱  ⋱"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D = Derivative() :  Laguerre(0) → Laguerre(1)\n",
    "S = I : Laguerre(0) → Laguerre(1)\n",
    "\n",
    "Jt = ApproxFun.Recurrence(Laguerre(0))\n",
    "\n",
    "D - S/2 -S*Jt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Problem 2.2\n",
    "\n",
    "We have\n",
    "\\begin{align*}\n",
    "{\\E^x \\over x^\\alpha} {\\D \\over \\dx} \\br[x^{\\alpha+1} \\E^{-x} {\\D L_n^{(\\alpha)} \\over \\dx}] &= -{\\E^x \\over x^\\alpha} {\\D \\over \\dx} \\br[x^{\\alpha+1} \\E^{-x} { L_{n-1}^{(\\alpha+1)} \\over \\dx}] \\\\\n",
    "&= -n L_{n}^{(\\alpha)}(x) \n",
    "\\end{align*}\n",
    "Therefore $\\lambda_n = -n$. This can be expanded in the form:\n",
    "$$\n",
    "x {\\D^2 L_n^{(\\alpha)} \\over \\dx^2} + (\\alpha+1 - x) {\\D L_n^{(\\alpha)} \\over \\dx} = -n L_n^{(\\alpha)}(x)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 3.1\n",
    "\n",
    "Define $C_k^{(\\alpha)}(z) = \\CC[L_k^{(\\alpha)} \\diamond^\\alpha \\E^{-\\diamond}](z)$ and recall that \n",
    "\\begin{align*}\n",
    "C_1(z)  &= {{1 \\over 2 \\pi \\I} \\int_0^\\infty \\E^{-x} \\dx + (z-a_0) C_0(z) \\over b_0} = -{1 \\over 2 \\pi \\I}  - (z-1) C_0(z)    \\\\\n",
    "&= {(z-1) \\E^{-z} \\Ei z -1 \\over 2 \\pi \\I} \n",
    "\\end{align*}\n",
    "\n",
    "Here we double check the formula, noting that\n",
    "$L_1(x) = \\E^x {\\D \\over \\dx} x \\E^{-x} = 1 - x$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ei (generic function with 1 method)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "const ei₋₁ = let ζ = Fun(-100 .. -1)\n",
    "    sum(exp(ζ)/ζ)\n",
    "end\n",
    "function ei(z) \n",
    "    ζ = Fun(Segment(-1 , z))\n",
    "    ei₋₁ + sum(exp(ζ)/ζ)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.01868464429845791 + 0.04836133565334999im"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = Fun(0..10)\n",
    "w = exp(-x)\n",
    "z = 1+im\n",
    "cauchy((1-x)*w, z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.018683923002628996 + 0.048368487990893105im"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "((z-1)*exp(-z)*ei(z)-1)/(2π*im)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now use these to determine the results with $\\alpha = 1$. Note that:\n",
    "$$\n",
    "C_0^{(1)}(z) = \\CC[\\diamond \\E^{-\\diamond}](z) =C_0(z) - C_1(z) = { \\E^{-z} \\Ei z - (z-1) \\E^{-z} \\Ei z +1 \\over 2 \\pi \\I} \n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.09210173751684986 - 0.02967669135489207im"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cauchy(x*w, z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.09210253209837324 - 0.02968456498826411im"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(-exp(-z)*ei(z)-(z-1)*exp(-z)*ei(z)+1)/(2π*im)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore, we have\n",
    "\\begin{align*}\n",
    "C_1^{(1)}(z) &= {{1 \\over 2 \\pi \\I} \\int_0^\\infty x \\E^{-x} \\dx + (z-a_0^{(1)}) C_0^{(1)}(z) \\over b_0^{(1)}} \\\\\n",
    "&=  {{1 \\over 2 \\pi \\I}  + (z-2) C_0^{(1)}(z) \\over -1}\\\\\n",
    "&=  {1  + (z-2) (\\E^{-z} \\Ei z - (z-1) \\E^{-z} \\Ei z +1) \\over -2 \\pi \\I}\n",
    "\\end{align*}\n",
    "\n",
    "Let's check the result using\n",
    "$$\n",
    "L_1^{(1)}(x) = x^{-1} \\E^x {\\D \\over \\dx} x^2 \\E^{-x} = 2 - x\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.062425046161957945 + 0.03729703236453844im"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cauchy((2-x)*x*w, z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.06241796711010913 + 0.037367846005258im"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(1+(z-2)*(-exp(-z)*ei(z)-(z-1)*exp(-z)*ei(z)+1))/(-2π*im)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 3.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have \n",
    "$$\n",
    "\\int_x^\\infty L_2(x) \\E^{-x} \\dx = {1 \\over 2} x \\E^{-x} L_1^{(1)}(x) \\E^{-x}\n",
    "$$\n",
    "Thus from lectures we have\n",
    "$$\n",
    "{1 \\over 2 \\pi\\I} \\int_0^\\infty L_2(x) \\E^{-x} \\log(z-x) \\dx = {1 \\over 2}\\CC[ \\diamond \\E^{-\\diamond} L_1^{(1)}](z)\n",
    "$$\n",
    "and therefore\n",
    "$$\n",
    "{1 \\over \\pi} \\int_0^\\infty L_2(x) \\E^{-x} \\log|z-x| \\dx  = -\\Im \\CC[ \\diamond \\E^{-\\diamond} L_1^{(1)}](z) = \\Re{1  + (z-2) (\\E^{-z} \\Ei z - (z-1) \\E^{-z} \\Ei z +1) \\over -2 \\pi }\n",
    "$$\n",
    "Let's check the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2 + 1im"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = Fun(0 .. 100)\n",
    "w = exp(-x)\n",
    "z = 2+im"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.06972323453771472"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "-sum(1/2*(2 - 4x + x^2)*w*log(abs(z-x)))/π"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.06972323453771528"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "imag(sum(1/2*(2 - 4x + x^2)*w*log(z-x))/(π*im))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.06972323454062201"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "-imag(cauchy((2-x)*x*w,z))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.06972323454061685"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "real((1+(z-2)*(-exp(-z)*ei(z)-(z-1)*exp(-z)*ei(z)+1))/(-2π))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 4.1\n",
    "\n",
    "Consider integration contours $\\gamma_{+x}$ and $\\gamma_{-x}$ that avoid $0$ above and below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 2400 1600\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip7600\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2000\" height=\"2000\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7601\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2400\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip7601)\" points=\"\n",
       "0,1600 2400,1600 2400,0 0,0 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7602\">\n",
       "    <rect x=\"480\" y=\"0\" width=\"1681\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip7601)\" points=\"\n",
       "189.504,1503.47 2321.26,1503.47 2321.26,47.2441 189.504,47.2441 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7603\">\n",
       "    <rect x=\"189\" y=\"47\" width=\"2133\" height=\"1457\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip7603)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,1503.47 189.504,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7603)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  722.443,1503.47 722.443,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7603)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1255.38,1503.47 1255.38,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7603)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1788.32,1503.47 1788.32,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7603)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  2321.26,1503.47 2321.26,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7603)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,1503.47 2321.26,1503.47 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7603)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,1139.42 2321.26,1139.42 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7603)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,775.359 2321.26,775.359 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7603)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,411.302 2321.26,411.302 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7603)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,47.2441 2321.26,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7601)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,1503.47 2321.26,1503.47 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7601)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,1503.47 189.504,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7601)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,1503.47 189.504,1481.63 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7601)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  722.443,1503.47 722.443,1481.63 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7601)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1255.38,1503.47 1255.38,1481.63 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7601)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1788.32,1503.47 1788.32,1481.63 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7601)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  2321.26,1503.47 2321.26,1481.63 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7601)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,1503.47 221.48,1503.47 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7601)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,1139.42 221.48,1139.42 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7601)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,775.359 221.48,775.359 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7601)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,411.302 221.48,411.302 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7601)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,47.2441 221.48,47.2441 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7601)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 189.504, 1557.47)\" x=\"189.504\" y=\"1557.47\">-5.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7601)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 722.443, 1557.47)\" x=\"722.443\" y=\"1557.47\">-2.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7601)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1255.38, 1557.47)\" x=\"1255.38\" y=\"1557.47\">0.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7601)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1788.32, 1557.47)\" x=\"1788.32\" y=\"1557.47\">2.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7601)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 2321.26, 1557.47)\" x=\"2321.26\" y=\"1557.47\">5.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7601)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 1520.97)\" x=\"165.504\" y=\"1520.97\">-1.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7601)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 1156.92)\" x=\"165.504\" y=\"1156.92\">-0.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7601)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 792.859)\" x=\"165.504\" y=\"792.859\">0.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7601)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 428.802)\" x=\"165.504\" y=\"428.802\">0.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7601)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 64.7441)\" x=\"165.504\" y=\"64.7441\">1.0</text>\n",
       "</g>\n",
       "<circle clip-path=\"url(#clip7603)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"829.031\" cy=\"775.359\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip7603)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"829.031\" cy=\"775.359\" r=\"14\"/>\n",
       "<polyline clip-path=\"url(#clip7603)\" style=\"stroke:#e26f46; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1233.98,775.359 1233.35,775.359 1232.08,775.359 1230.17,775.359 1227.65,775.359 1224.5,775.359 1220.75,775.359 1216.4,775.359 1211.47,775.359 1205.98,775.359 \n",
       "  1199.93,775.359 1193.36,775.359 1186.28,775.359 1178.71,775.359 1170.68,775.359 1162.21,775.359 1153.32,775.359 1144.06,775.359 1134.44,775.359 1124.5,775.359 \n",
       "  1114.26,775.359 1103.77,775.359 1093.05,775.359 1082.13,775.359 1071.06,775.359 1059.86,775.359 1048.57,775.359 1037.23,775.359 1025.87,775.359 1014.53,775.359 \n",
       "  1003.24,775.359 992.039,775.359 980.964,775.359 970.048,775.359 959.326,775.359 948.831,775.359 938.596,775.359 928.654,775.359 919.035,775.359 909.771,775.359 \n",
       "  900.889,775.359 892.419,775.359 884.387,775.359 876.817,775.359 869.735,775.359 863.161,775.359 857.117,775.359 851.622,775.359 846.694,775.359 842.347,775.359 \n",
       "  838.595,775.359 835.45,775.359 832.922,775.359 831.019,775.359 829.747,775.359 829.111,775.359 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7603)\" style=\"stroke:#e26f46; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1276.7,775.356 1276.7,775.331 1276.7,775.282 1276.7,775.207 1276.7,775.107 1276.7,774.983 1276.7,774.833 1276.7,774.658 1276.7,774.457 1276.7,774.232 \n",
       "  1276.7,773.98 1276.69,773.703 1276.69,773.4 1276.69,773.07 1276.69,772.715 1276.68,772.332 1276.68,771.923 1276.67,771.487 1276.66,771.024 1276.65,770.533 \n",
       "  1276.64,770.014 1276.63,769.467 1276.62,768.892 1276.6,768.288 1276.58,767.656 1276.56,766.994 1276.53,766.302 1276.51,765.581 1276.48,764.831 1276.44,764.049 \n",
       "  1276.4,763.238 1276.36,762.396 1276.31,761.523 1276.26,760.619 1276.2,759.684 1276.14,758.718 1276.06,757.72 1275.99,756.692 1275.9,755.633 1275.81,754.542 \n",
       "  1275.71,753.421 1275.6,752.27 1275.48,751.089 1275.35,749.878 1275.21,748.639 1275.06,747.371 1274.9,746.076 1274.72,744.755 1274.54,743.408 1274.34,742.038 \n",
       "  1274.12,740.645 1273.89,739.231 1273.64,737.799 1273.38,736.349 1273.1,734.884 1272.81,733.408 1272.49,731.921 1272.16,730.428 1271.8,728.93 1271.43,727.433 \n",
       "  1271.04,725.938 1270.62,724.451 1270.19,722.974 1269.73,721.512 1269.25,720.07 1268.75,718.652 1268.23,717.262 1267.69,715.906 1267.12,714.588 1266.54,713.313 \n",
       "  1265.93,712.087 1265.3,710.915 1264.66,709.801 1263.99,708.751 1263.31,707.768 1262.61,706.859 1261.89,706.026 1261.16,705.275 1260.42,704.609 1259.66,704.032 \n",
       "  1258.9,703.546 1258.13,703.154 1257.35,702.858 1256.56,702.66 1255.78,702.56 1254.99,702.56 1254.2,702.66 1253.42,702.858 1252.64,703.154 1251.86,703.546 \n",
       "  1251.1,704.032 1250.35,704.609 1249.6,705.275 1248.87,706.026 1248.16,706.859 1247.46,707.768 1246.77,708.751 1246.11,709.801 1245.46,710.915 1244.83,712.087 \n",
       "  1244.23,713.313 1243.64,714.588 1243.08,715.906 1242.53,717.262 1242.01,718.652 1241.51,720.07 1241.03,721.512 1240.58,722.974 1240.14,724.451 1239.73,725.938 \n",
       "  1239.33,727.433 1238.96,728.93 1238.61,730.428 1238.27,731.921 1237.96,733.408 1237.66,734.884 1237.38,736.349 1237.12,737.799 1236.87,739.231 1236.64,740.645 \n",
       "  1236.43,742.038 1236.23,743.408 1236.04,744.755 1235.86,746.076 1235.7,747.371 1235.55,748.639 1235.41,749.878 1235.28,751.089 1235.16,752.27 1235.06,753.421 \n",
       "  1234.95,754.542 1234.86,755.633 1234.78,756.692 1234.7,757.72 1234.63,758.718 1234.56,759.684 1234.51,760.619 1234.45,761.523 1234.41,762.396 1234.36,763.238 \n",
       "  1234.32,764.049 1234.29,764.831 1234.26,765.581 1234.23,766.302 1234.21,766.994 1234.18,767.656 1234.17,768.288 1234.15,768.892 1234.13,769.467 1234.12,770.014 \n",
       "  1234.11,770.533 1234.1,771.024 1234.09,771.487 1234.09,771.923 1234.08,772.332 1234.08,772.715 1234.07,773.07 1234.07,773.4 1234.07,773.703 1234.07,773.98 \n",
       "  1234.07,774.232 1234.07,774.457 1234.07,774.658 1234.06,774.833 1234.06,774.983 1234.06,775.107 1234.06,775.207 1234.06,775.282 1234.06,775.331 1234.06,775.356 \n",
       "  \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7603)\" style=\"stroke:#e26f46; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  22568.8,775.359 22535.3,775.359 22468.4,775.359 22368.3,775.359 22235.4,775.359 22070.1,775.359 21872.8,775.359 21644.2,775.359 21385.1,775.359 21096.2,775.359 \n",
       "  20778.4,775.359 20432.8,775.359 20060.4,775.359 19662.4,775.359 19240,775.359 18794.7,775.359 18327.7,775.359 17840.6,775.359 17334.9,775.359 16812.1,775.359 \n",
       "  16274,775.359 15722.2,775.359 15158.4,775.359 14584.5,775.359 14002.2,775.359 13413.3,775.359 12819.8,775.359 12223.5,775.359 11626.2,775.359 11029.8,775.359 \n",
       "  10436.3,775.359 9847.47,775.359 9265.17,775.359 8691.23,775.359 8127.47,775.359 7575.66,775.359 7037.53,775.359 6514.77,775.359 6009.04,775.359 5521.92,775.359 \n",
       "  5054.95,775.359 4609.59,775.359 4187.25,775.359 3789.26,775.359 3416.86,775.359 3071.23,775.359 2753.46,775.359 2464.54,775.359 2205.39,775.359 1976.82,775.359 \n",
       "  1779.56,775.359 1614.21,775.359 1481.3,775.359 1381.25,775.359 1314.38,775.359 1280.89,775.359 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7603)\" style=\"stroke:#3da44d; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1233.98,775.359 1233.35,775.359 1232.08,775.359 1230.17,775.359 1227.65,775.359 1224.5,775.359 1220.75,775.359 1216.4,775.359 1211.47,775.359 1205.98,775.359 \n",
       "  1199.93,775.359 1193.36,775.359 1186.28,775.359 1178.71,775.359 1170.68,775.359 1162.21,775.359 1153.32,775.359 1144.06,775.359 1134.44,775.359 1124.5,775.359 \n",
       "  1114.26,775.359 1103.77,775.359 1093.05,775.359 1082.13,775.359 1071.06,775.359 1059.86,775.359 1048.57,775.359 1037.23,775.359 1025.87,775.359 1014.53,775.359 \n",
       "  1003.24,775.359 992.039,775.359 980.964,775.359 970.048,775.359 959.326,775.359 948.831,775.359 938.596,775.359 928.654,775.359 919.035,775.359 909.771,775.359 \n",
       "  900.889,775.359 892.419,775.359 884.387,775.359 876.817,775.359 869.735,775.359 863.161,775.359 857.117,775.359 851.622,775.359 846.694,775.359 842.347,775.359 \n",
       "  838.595,775.359 835.45,775.359 832.922,775.359 831.019,775.359 829.747,775.359 829.111,775.359 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7603)\" style=\"stroke:#3da44d; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1276.7,775.363 1276.7,775.387 1276.7,775.437 1276.7,775.512 1276.7,775.611 1276.7,775.736 1276.7,775.886 1276.7,776.061 1276.7,776.261 1276.7,776.487 \n",
       "  1276.7,776.739 1276.69,777.016 1276.69,777.319 1276.69,777.648 1276.69,778.004 1276.68,778.386 1276.68,778.796 1276.67,779.232 1276.66,779.695 1276.65,780.186 \n",
       "  1276.64,780.705 1276.63,781.251 1276.62,781.827 1276.6,782.43 1276.58,783.063 1276.56,783.725 1276.53,784.416 1276.51,785.137 1276.48,785.888 1276.44,786.669 \n",
       "  1276.4,787.481 1276.36,788.323 1276.31,789.196 1276.26,790.1 1276.2,791.035 1276.14,792.001 1276.06,792.998 1275.99,794.027 1275.9,795.086 1275.81,796.177 \n",
       "  1275.71,797.297 1275.6,798.449 1275.48,799.63 1275.35,800.841 1275.21,802.08 1275.06,803.348 1274.9,804.643 1274.72,805.964 1274.54,807.311 1274.34,808.681 \n",
       "  1274.12,810.074 1273.89,811.488 1273.64,812.92 1273.38,814.37 1273.1,815.834 1272.81,817.311 1272.49,818.798 1272.16,820.291 1271.8,821.788 1271.43,823.286 \n",
       "  1271.04,824.78 1270.62,826.268 1270.19,827.745 1269.73,829.207 1269.25,830.649 1268.75,832.067 1268.23,833.457 1267.69,834.813 1267.12,836.131 1266.54,837.406 \n",
       "  1265.93,838.631 1265.3,839.804 1264.66,840.918 1263.99,841.968 1263.31,842.951 1262.61,843.86 1261.89,844.692 1261.16,845.443 1260.42,846.109 1259.66,846.687 \n",
       "  1258.9,847.173 1258.13,847.565 1257.35,847.861 1256.56,848.059 1255.78,848.159 1254.99,848.159 1254.2,848.059 1253.42,847.861 1252.64,847.565 1251.86,847.173 \n",
       "  1251.1,846.687 1250.35,846.109 1249.6,845.443 1248.87,844.692 1248.16,843.86 1247.46,842.951 1246.77,841.968 1246.11,840.918 1245.46,839.804 1244.83,838.631 \n",
       "  1244.23,837.406 1243.64,836.131 1243.08,834.813 1242.53,833.457 1242.01,832.067 1241.51,830.649 1241.03,829.207 1240.58,827.745 1240.14,826.268 1239.73,824.78 \n",
       "  1239.33,823.286 1238.96,821.788 1238.61,820.291 1238.27,818.798 1237.96,817.311 1237.66,815.834 1237.38,814.37 1237.12,812.92 1236.87,811.488 1236.64,810.074 \n",
       "  1236.43,808.681 1236.23,807.311 1236.04,805.964 1235.86,804.643 1235.7,803.348 1235.55,802.08 1235.41,800.841 1235.28,799.63 1235.16,798.449 1235.06,797.297 \n",
       "  1234.95,796.177 1234.86,795.086 1234.78,794.027 1234.7,792.998 1234.63,792.001 1234.56,791.035 1234.51,790.1 1234.45,789.196 1234.41,788.323 1234.36,787.481 \n",
       "  1234.32,786.669 1234.29,785.888 1234.26,785.137 1234.23,784.416 1234.21,783.725 1234.18,783.063 1234.17,782.43 1234.15,781.827 1234.13,781.251 1234.12,780.705 \n",
       "  1234.11,780.186 1234.1,779.695 1234.09,779.232 1234.09,778.796 1234.08,778.386 1234.08,778.004 1234.07,777.648 1234.07,777.319 1234.07,777.016 1234.07,776.739 \n",
       "  1234.07,776.487 1234.07,776.261 1234.07,776.061 1234.06,775.886 1234.06,775.736 1234.06,775.611 1234.06,775.512 1234.06,775.437 1234.06,775.387 1234.06,775.363 \n",
       "  \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7603)\" style=\"stroke:#3da44d; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  22568.8,775.359 22535.3,775.359 22468.4,775.359 22368.3,775.359 22235.4,775.359 22070.1,775.359 21872.8,775.359 21644.2,775.359 21385.1,775.359 21096.2,775.359 \n",
       "  20778.4,775.359 20432.8,775.359 20060.4,775.359 19662.4,775.359 19240,775.359 18794.7,775.359 18327.7,775.359 17840.6,775.359 17334.9,775.359 16812.1,775.359 \n",
       "  16274,775.359 15722.2,775.359 15158.4,775.359 14584.5,775.359 14002.2,775.359 13413.3,775.359 12819.8,775.359 12223.5,775.359 11626.2,775.359 11029.8,775.359 \n",
       "  10436.3,775.359 9847.47,775.359 9265.17,775.359 8691.23,775.359 8127.47,775.359 7575.66,775.359 7037.53,775.359 6514.77,775.359 6009.04,775.359 5521.92,775.359 \n",
       "  5054.95,775.359 4609.59,775.359 4187.25,775.359 3789.26,775.359 3416.86,775.359 3071.23,775.359 2753.46,775.359 2464.54,775.359 2205.39,775.359 1976.82,775.359 \n",
       "  1779.56,775.359 1614.21,775.359 1481.3,775.359 1381.25,775.359 1314.38,775.359 1280.89,775.359 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip7601)\" points=\"\n",
       "1793.91,372.684 2249.26,372.684 2249.26,130.764 1793.91,130.764 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip7601)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1793.91,372.684 2249.26,372.684 2249.26,130.764 1793.91,130.764 1793.91,372.684 \n",
       "  \"/>\n",
       "<circle clip-path=\"url(#clip7601)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"1901.91\" cy=\"191.244\" r=\"25\"/>\n",
       "<circle clip-path=\"url(#clip7601)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"1901.91\" cy=\"191.244\" r=\"21\"/>\n",
       "<g clip-path=\"url(#clip7601)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 1985.91, 208.744)\" x=\"1985.91\" y=\"208.744\">x = -2.0</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip7601)\" style=\"stroke:#e26f46; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1817.91,251.724 1961.91,251.724 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7601)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 1985.91, 269.224)\" x=\"1985.91\" y=\"269.224\">gamma_+</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip7601)\" style=\"stroke:#3da44d; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1817.91,312.204 1961.91,312.204 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7601)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 1985.91, 329.704)\" x=\"1985.91\" y=\"329.704\">gamma_-</text>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = -2.0\n",
    "r = 0.1\n",
    "γ₊ₓ = Segment(-2.0 , -r) ∪ Arc(0.,r, (π,0)) ∪ Segment(r , 100)\n",
    "γ₋ₓ = Segment(-2.0 , -r) ∪ Arc(0.,r, (-π,0)) ∪ Segment(r , 100)\n",
    "scatter([x],[0.0];label=\"x = $x\")\n",
    "plot!(γ₊ₓ ; xlims=(-5,5), ylims=(-1,1), label=\"gamma_+\")\n",
    "plot!(γ₋ₓ; xlims=(-5,5), ylims=(-1,1), label=\"gamma_-\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So that\n",
    "$$\n",
    "\\Gamma_\\pm(\\alpha, x) = \\int_{\\gamma_{\\pm x}} \\zeta^{\\alpha-1} \\E^{-\\zeta} \\D \\zeta\n",
    "$$\n",
    "Note that\n",
    "$$\n",
    "\\int_x^{-r} (\\zeta_+^{\\alpha - 1} - \\zeta_-^{\\alpha -1}) \\E^{-\\zeta} \\D\\zeta = 0\n",
    "$$\n",
    "since $\\zeta_+^{\\alpha-1} = \\E^{\\pi \\I (\\alpha-1)}|\\zeta|^{\\alpha-1} = \\E^{2 \\I \\pi \\alpha} \\zeta_-^{\\alpha-1}$.\n",
    "Furthermore, the integrals over the arcs tend to zero as $r \\rightarrow 0$:\n",
    "$$\n",
    "|\\I r^\\alpha \\int_0^\\pi \\E^{- r \\E^{\\I \\theta}} \\E^{\\I \\theta \\alpha} \\D \\theta  | \\leq r^\\alpha \\pi \\E^r  \\rightarrow 0\n",
    "$$\n",
    "and similarly on the lower arc. Thus we have\n",
    "$$\n",
    "\\Gamma_+(\\alpha, x)-\\E^{2 \\I \\pi \\alpha} \\Gamma_-(\\alpha, x) = \\lim_{r \\rightarrow 0 } \\left(\\int_{\\gamma_{+x}} - \\E^{2 \\I \\pi \\alpha} \\int_{\\gamma_{-x}}\\right) \\zeta^{\\alpha-1} \\E^{-\\zeta} \\D \\zeta = (1 - \\E^{2 \\I \\pi \\alpha})\\int_0^\\infty x^{\\alpha-1} \\E^{-x} \\dx = (1 - \\E^{2 \\I \\pi \\alpha})\\Gamma(\\alpha)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 4.2\n",
    "\n",
    "Note that, for $0 < \\alpha < 1$,\n",
    "$$\n",
    "    \\psi(z) = z^{-\\alpha} \\E^z \\Gamma(\\alpha, z)\n",
    "$$\n",
    "has the following properties:\n",
    "1. $\\psi(z)$ decays as $z \\rightarrow \\infty$, via integration by parts:\n",
    "$$\n",
    "    z^{-\\alpha} \\E^z \\int_z^\\infty \\zeta^{\\alpha-1} \\E^{-\\zeta} \\D \\zeta =\n",
    "    z^{-1} \\E^z + z^{-\\alpha} \\int_z^\\infty \\zeta^{\\alpha-2} \\E^{z-\\zeta} \\D \\zeta\n",
    "$$\n",
    "and we have assuming $z$ is bounded away from the negative real axis:\n",
    "$$\n",
    "|\\int_z^\\infty \\zeta^{\\alpha-2} \\E^{z-\\zeta} \\D \\zeta| \\leq \\int_z^\\infty |\\zeta|^{\\alpha-2} \\D \\zeta = \\int_0^\\infty |x+z|^{\\alpha-2} \\D x  < \\infty\n",
    "$$\n",
    "(otherwise one would use a deformed contour). \n",
    "2. We have the subtractive jump:\n",
    "$$\\psi_+(x) - \\psi_-(x) = \\E^x(x_+^{-\\alpha} \\Gamma_+(\\alpha, x) -  \\Gamma_-(\\alpha, x)) = \\E^x |x|^\\alpha (\\E^{-\\I \\pi \\alpha} \\Gamma_+(\\alpha, z) - \\E^{\\I \\pi \\alpha} \\Gamma_-(\\alpha,x)) = \\E^x |x|^\\alpha \\E^{-\\I \\pi \\alpha} (1-\\E^{2 \\I \\pi \\alpha})\n",
    "$$\n",
    "\n",
    "We use these properties to verify that\n",
    "$$\n",
    "\\CC[\\diamond^\\alpha \\E^{-\\diamond}](z) = {1 \\over \\Gamma(-\\alpha)} {(-z)^\\alpha \\E^{-z} \\Gamma(-\\alpha, - z) \\over\n",
    "  \\E^{-\\I\\pi\\alpha} - \\E^{\\I\\pi\\alpha}}\n",
    "$$\n",
    "via Plemelj."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.07199876331437786 + 0.05850612397322342im"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = Fun(0 .. 20.0)\n",
    "α = -0.1\n",
    "z = 2.0+im\n",
    "cauchy(x^α*exp(-x), z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "#3 (generic function with 1 method)"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Γ = (α,z) -> let ζ = z + Fun(0 .. 500.0)\n",
    "    linesum(ζ^(α-1)*exp(-ζ))\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.07199876331505133 + 0.0585061239604885im"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "-(-z)^α*exp(-z)Γ(-α,-z)/(gamma(-α)*(exp(im*π*α)-exp(-im*π*α)))"
   ]
  }
 ],
 "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.1.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}