{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# The Airy function Riemann-Hilbert problem" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "using ApproxFun, SpecialFunctions, SingularIntegralEquations, RiemannHilbert, Plots" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Consider the Airy equation __[DLMF Sec. 9]__ $\\newcommand{\\Ai}{\\mathrm{Ai}}$$\\newcommand{\\mf}{\\mathfrak} \\newcommand{\\E}{\\mathrm{e}} \\newcommand{\\I}{\\mathrm{i}} \\newcommand{\\mf}{\\mathfrak} \\DeclareMathOperator{\\arg}{arg}\n", "$$\n", "Y''(z) = z Y(z).\n", "$$\n", "The canonical particular solution satisfies\n", "$$\n", "\\Ai(z) \\sim \\frac{1}{2 \\sqrt{\\pi}} z^{-1/4} \\E^{-\\frac 2 3 z^{3/2}}, \\quad - \\frac {2\\pi}{3} \\leq \\mathrm{arg}\\, z \\le \\frac{2\\pi}{3}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Define$\\omega = \\E^{- \\frac {2 \\I \\pi}{3}}and\n", "\\begin{align*}\n", " y(z) = \\begin{cases} \\begin{bmatrix} - \\omega \\Ai(\\omega z) & - \\I \\omega^2 \\Ai(\\omega^2 z) \\end{bmatrix} & - \\pi < \\arg z < - \\frac{2\\pi}{3}, \\\\\n", " \\begin{bmatrix} \\Ai(z) & - \\I \\omega^2 \\Ai(\\omega^2 z) \\end{bmatrix} & - \\frac{2 \\pi}{3} < \\arg z < 0,\\\\\n", " \\begin{bmatrix} \\Ai(z) & \\I \\omega \\Ai(\\omega z) \\end{bmatrix} & 0 < \\arg z < \\frac{2\\pi}{3},\\\\\n", " \\begin{bmatrix} - \\omega^2 \\Ai(\\omega^2 z) & \\I \\omega \\Ai(\\omega z) \\end{bmatrix} & \\frac{2\\pi}{3} < \\arg z < \\pi.\n", " \\end{cases}\n", "\\end{align*}" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Using the fact that\\Ai(z) + \\omega\\Ai(\\omega z) + \\omega^2\\Ai(\\omega^2 z) = 0$, it follows that$yhas the following jump conditions\n", "\\begin{align*}\n", "y^+(s) &= y^-(s) \\begin{bmatrix} 1 & -\\I \\\\ 0 & 1 \\end{bmatrix} \\quad s \\in (0,\\infty),\\\\\n", "y^+(s) &= y^-(s) \\begin{bmatrix} 1 & 0 \\\\ -\\I & 1 \\end{bmatrix} \\quad s \\in (0,\\omega^{-1}\\infty),\\\\\n", "y^+(s) &= y^-(s) \\begin{bmatrix} 1 & 0 \\\\ -\\I & 1 \\end{bmatrix} \\quad s \\in (0,\\omega\\infty),\\\\\n", "y^+(s) &= y^-(s) \\begin{bmatrix} 0 & -\\I \\\\ -\\I & 0 \\end{bmatrix} \\quad s \\in (-\\infty,0).\n", "\\end{align*}" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ω = exp(-2im*pi/3); ωi = exp(2im*pi/3);\n", "Γ₀ = Segment(0.,5) ∪ Segment(0,5ωi) ∪ Segment(0,5ω) ∪ Segment(-5,0)\n", "plot(Γ₀, arrow=:arrow, legend=false)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "This setup has two related problems:\n", "* The jump matrices do not tend to the identity matrix at infinity.\n", "* The asymptotics ofyat infinity is not, to leading order, constant." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "To fix some of these issues, define\n", "\\begin{align*}\n", "\\phi(z) = 2 \\sqrt \\pi z^{1/4} y(z) \\begin{bmatrix} \\E^{\\frac{2}{3} z^{3/2}} & 0 \\\\ 0 & \\E^{-\\frac{2}{3} z^{3/2} }\\end{bmatrix}.\n", "\\end{align*}" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Then\n", "\\begin{align*}\n", "\\phi^+(s) &= \\phi^-(s) \\begin{bmatrix} 1 & -\\I \\E^{- \\frac 4 3 s^{3/2}} \\\\ 0 & 1 \\end{bmatrix} \\quad s \\in (0,\\infty),\\\\\n", "\\phi^+(s) &= \\phi^-(s) \\begin{bmatrix} 1 & 0 \\\\ -\\I\\E^{ \\frac 4 3 s^{3/2}} & 1 \\end{bmatrix} \\quad s \\in (0,\\omega^{-1}\\infty),\\\\\n", "\\phi^+(s) &= \\phi^-(s) \\begin{bmatrix} 1 & 0 \\\\ -\\I\\E^{- \\frac 4 3 s^{3/2}} & 1 \\end{bmatrix} \\quad s \\in (0,\\omega\\infty),\\\\\n", "\\phi^+(s) &= \\phi^-(s) \\begin{bmatrix} 0 & -\\I \\\\ -\\I & 0 \\end{bmatrix} \\quad s \\in (-\\infty,0).\n", "\\end{align*}" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "And\n", "\\begin{align*}\n", "\\phi(z) \\sim \\begin{bmatrix} 1 & 1 \\end{bmatrix}, \\quad z \\to \\infty.\n", "\\end{align*}" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "We solved some of the problems with the jumps but we have an additional (major) problem:\n", "* Sincey(z)$is analytically extendible beyond each sector,$\\phi(z)$must be singular at$z = 0." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "To fix this, set\n", "\\begin{align*}\n", " J_1 = \\begin{bmatrix} 1 & 0 \\\\ -\\I & 1 \\end{bmatrix}, \\quad J_2 = \\begin{bmatrix} 1 & -\\I \\\\ 0 & 1 \\end{bmatrix}.\n", "\\end{align*}" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Instead of considering\\phi$, first consider$y_1$defined as follows\n", "\n", " ![title](y1.png) " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "$y_1 = youtside this circle." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "And set\n", "\\begin{align*}\n", " \\phi_1(z) = y_1(z) \\begin{cases} \n", " R(z) & |z| > 1 \\\\ \\\\\n", " I & |z| \\leq 1, \\end{cases}\n", "\\end{align*}\n", "where\n", "\\begin{align*}\n", "R(z) = z^{1/4} \\begin{bmatrix} \\E^{\\frac 2 3 z^{3/2}} & 0 \\\\ 0 & \\E^{- \\frac 2 3 z^{3/2}} \\end{bmatrix} \\begin{bmatrix} 1 & 1 - \\frac{1}{\\sqrt{z}} \\\\ 1 & 1 + \\frac{1}{\\sqrt{z}} \\end{bmatrix}.\n", "\\end{align*}" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ " ![title](airy_rhp.png) " ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r = 1.\n", "L = 15.\n", "Γ = Segment(r, L) ∪\n", " Segment(r*ωi, L*ωi) ∪\n", " Segment(r*ω, L*ω) ∪\n", " Segment(r, r*ωi) ∪\n", " Segment(r*ωi, r*ω) ∪\n", " Segment(r*ω, r)\n", "plot(Γ, arrow=:arrow, legend = false)" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "ϵ = 2.2e-16\n", "H(z) = 0.5*[1 (1 - 1/sqrt(z));\n", " 1 (1 + 1/sqrt(z))]\n", "H(z,s) = H(z + s*1im*ϵ)\n", "Hi(z) = inv(H(z))\n", "Hi(z,s) = inv(H(z,s))\n", "R(z) = 2*sqrt(π)*z^(1/4)*[exp(2/3*z^(3/2)) 0.;\n", " 0. exp(-2/3*z^(3/2))]*H(z)\n", "R(z,s) = R(z + s*1im*ϵ)\n", "Ri(z) = 1/(2*sqrt(π))*z^(-1/4)*Hi(z)*[exp(-2/3*z^(3/2)) 0.0;\n", " 0.0 exp(2/3*z^(3/2))]\n", "Ri(z,s) = Ri(z + s*1im*ϵ)\n", "\n", "J1 = [1 0; -1im 1]\n", "J2 = [1 -1im; 0 1]\n", "G0 = z -> if angle(z) ≈ 0\n", " Ri(z,1)*J2*R(z,1)\n", " elseif angle(z) ≈ 2π/3\n", " Ri(z,1)*J1*R(z,1)\n", " elseif angle(z) ≈ -2π/3\n", " Ri(z,1)*J1*R(z,1)\n", " end\n", "G1 = z -> Ri(z,1)*J1\n", "G2 = z -> imag(z) >= 0. ? Ri(z,1) : Ri(z,1)*J1*J2*J1\n", "G3 = z -> Ri(z,1)*J2*J1;" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Set up Riemann-Hilbert problem" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "rhp = Fun(function(z)\n", " z ∈ component(Γ,1) && return G0(z)\n", " z ∈ component(Γ,2) && return G0(z)\n", " z ∈ component(Γ,3) && return G0(z)\n", " z ∈ component(Γ,4) && return G1(z)\n", " z ∈ component(Γ,5) && return G2(z)\n", " z ∈ component(Γ,6) && return G3(z)\n", " error(\"Not in contour\") \n", " end, Γ);" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Check product condition atω r, \\omega^{-1} r, r,ω L, \\omega^{-1} L, L\$ and solve" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "true true true true true true " ] } ], "source": [ "for ℓ in [r,L]\n", " for z in [ω,ωi,1]*ℓ\n", " print(RiemannHilbert.productcondition(rhp,z) ≈ [1 0; 0 1],\" \")\n", " end\n", "end" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [], "source": [ "Φ = transpose(rhsolve(transpose(rhp),6*120));\n", "Φ = Φ[1,:] + Φ[2,:];" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "Q = z -> if 0 < angle(z) < 2π/3\n", " inv(J1)\n", " elseif 2π/3 < angle(z) <= π\n", " [1. 0im; 0im 1.]\n", " elseif -π < angle(z) < -2π/3\n", " inv(J1*J2*J1)\n", " elseif -2π/3 <= angle(z) <= 0\n", " inv(J2*J1)\n", " end\n", "P = z -> if real(z) >= -r*cos(π/3) && imag(z*exp(1im*π/3)) <= r*cos(pi/3) && imag(z*exp(-1im*π/3)) >= -r*cos(pi/3) \n", " Q(z)\n", " else \n", " Ri(z,1)\n", " end\n", "y = z -> transpose(Φ(z))*P(z);\n", "Airy = z -> if -2π/3 <= angle(z) <= 2π/3\n", " y(z + 1im*ϵ)\n", " elseif 2π/3 < angle(z) <= π\n", " yy = y(z + (1im-1)*ϵ)\n", " yy + 1im*yy\n", " else\n", " yy = y(z - (1im+1)*ϵ)\n", " yy - 1im*yy\n", " end;" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "text/plain": [ "8.55208292594855e-16" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "z = 1im; abs((Airy(z) - airyai(z))/airyai(z))" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2.070721306648171e-13" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "z = -100-1im; abs((Airy(z) - airyai(z))/airyai(z))" ] } ], "metadata": { "celltoolbar": "Slideshow", "kernelspec": { "display_name": "julia 1.3.0", "language": "julia", "name": "julia-1.3" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.3.0" } }, "nbformat": 4, "nbformat_minor": 2 }