{
"cells": [
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"using Plots, ComplexPhasePortrait, ApproxFun, SingularIntegralEquations, OscillatoryIntegrals\n",
"import ApproxFun: UnionDomain\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\\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\\qqand{\\qquad\\hbox{for}\\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\\vc#1{{\\mathbf #1}}\n",
"\\def\\ip<#1,#2>{\\left\\langle#1,#2\\right\\rangle}\n",
"\\def\\norm#1{\\left\\|#1\\right\\|}\n",
"\\def\\half{{1 \\over 2}}\n",
"$$\n",
"\n",
"Dr. Sheehan Olver\n",
"
\n",
"s.olver@imperial.ac.uk\n",
"\n",
"\n",
"
\n",
"Website: https://github.com/dlfivefifty/M3M6LectureNotes\n",
"\n",
"# Chapter 4: The Wiener–Hopf method\n",
"\n",
"The Wiener–Hopf method is an approach for solving integral equations of the form\n",
"$$\n",
"\\lambda u(x) + \\int_0^\\infty K(x-t)u(t) \\dt = f(x)\\qqfor 0 \\leq x < \\infty.\n",
"$$\n",
"Here $f$ and $K$ are given functions, $\\lambda$ is a given constant, and we want to find $u$. The approach similarly extends to integro-differential equations of the form\n",
"$$\n",
"c_1 u''(x) + c_2 u'(x) + \\lambda u(x) + \\int_0^\\infty K(x-t)u(t) \\dt = f(x)\\qqfor 0 \\leq x < \\infty.\n",
"$$\n",
"Here $c_1, c_2, \\lambda$ are given constants and again we want to find $u$.\n",
"\n",
"\n",
"In this course, we will only consider $K(x) = \\E^{-\\gamma |x|}$, though the methodology translates to other choices of $K$. \n",
"\n",
"\n",
"# Lecture 20: Integral equations on the real line\n",
"\n",
"Before discussing problems on the half line, we consider problems on the whole line:\n",
"$$\n",
"\\lambda u(x) + \\int_{-\\infty}^\\infty K(x-t)u(t) \\dt = f(x)\\qqfor -\\infty < x < \\infty.\n",
"$$\n",
"Outline:\n",
"\n",
"1. Fourier transforms and convolutions\n",
"2. Some simple examples\n",
"\n",
"## Fourier transforms and convolutions\n",
"\n",
"\n",
"Recall the Fourier transform of a function (using $s$ so I don't get confused with $K$):\n",
"$$\n",
" {\\hat{K}}(s) = {\\cal F} K (s) = \\int_{-\\infty}^\\infty K(t) \\E^{-\\I s t} \\dt\n",
"$$\n",
"so that\n",
"$$\n",
"K(x) = {\\cal F}^{-1} \\hat K(x) = {1 \\over 2\\pi}\\int_{-\\infty}^\\infty \\hat K(s) \\E^{ \\I s x} \\D s\n",
"$$\n",
"What effect does translating a function have on it's Fourier transform? \n",
"$$\n",
"\\int_{-\\infty}^\\infty K(t+x_0) \\E^{-\\I s t} \\dt = \\int_{-\\infty}^\\infty K(\\tau) \\E^{-\\I s (\\tau - x_0)} \\dt = \\E^{\\I s x_0} \\hat K(s)\n",
"$$\n",
"\n",
"We are interested in the Fourier transform of the _convolution_ term \n",
"$$\n",
"\\int_{-\\infty}^\\infty K(x-t)u(t) \\dt \n",
"$$\n",
"But we have\n",
"$$\n",
"\\int_{-\\infty}^\\infty \\E^{-\\I s x} \\int_{-\\infty}^\\infty K(x-t)u(t) \\dt \\dx = \\int_{-\\infty}^\\infty \\int_{-\\infty}^\\infty \\E^{-\\I s x} K(x-t)\\dx u(t) \\dt =\n",
"\\int_{-\\infty}^\\infty \\E^{-\\I s t} \\hat K(s) u(t) \\dt = \\hat K(s) \\hat u(s)\n",
"$$\n",
"That is, convolution becomes multiplication. Thus our integral equation in Fourier space becomes:\n",
"$$\n",
"(\\lambda + \\hat K(s)) \\hat u(s) = \\hat f(s) \\qqfor -\\infty < s < \\infty\n",
"$$\n",
"Or in other words,\n",
"$$\n",
"u(x) = {\\cal F}^{-1} \\left[{\\hat f \\over \\lambda + \\hat K }\\right](x)\n",
"$$\n",
"\n",
"### Example: \n",
"\n",
"Consider the kernel $K(x) = \\E^{-\\gamma |x|}$. Then we have\n",
"\\begin{align*}\n",
"\\hat K(s) &= \\int_{-\\infty}^\\infty K(t) \\E^{-\\I s t} \\dt = \\int_{-\\infty}^0 \\E^{( \\gamma - \\I s ) t} \\dt + \\int_0^\\infty \\E^{(-\\gamma - \\I s ) t} = {1 \\over \\gamma - \\I s} - {1 \\over -\\gamma - \\I s } \\\\\n",
"&= {2\\gamma \\over \\gamma^2 + s^2}\n",
"\\end{align*}"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-3.553823901825126e-13 - 8.326672684688674e-17im"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"γ = 2.0\n",
"Γ = ApproxFun.UnionDomain((-Inf .. 0) , (0 .. Inf)) # Line split in two\n",
"K = Fun(x -> exp(-γ*abs(x)), Γ)\n",
"K̂ = s -> fourier(K, -s) # magic routine for Fourier transforms\n",
"K̂(2.0) - 2γ/(γ^2+2.0^2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider a simple RHS like\n",
"$$\n",
"f(x) = {1 \\over x^2 + a^2}\n",
"$$\n",
"Because of the decay, we can use Residue calculus to determine, for $s < 0$ deforming to the upper half plane,\n",
"$$\n",
"\\hat f(s) = \\int_{-\\infty}^\\infty {\\E^{-\\I s t} \\over t^2 + a^2}\\dt = 2 \\pi \\I \\Res_{z = a \\I}{\\E^{-\\I s z} \\over z^2 + a^2} = \\pi {\\E^{s a} \\over a}\n",
"$$\n",
"and for $s > 0$ deforming to the lower half plane\n",
"$$\n",
"\\hat f(s) =-2 \\pi \\I \\Res_{z =- a \\I}{\\E^{-\\I s z} \\over z^2 + a^2} = \\pi {\\E^{-s a} \\over a}\n",
"$$\n",
"and for $s = 0$ the results match:\n",
"$$\n",
"\\hat f(0) = {\\pi \\over a}\n",
"$$\n",
"\n",
"**Remark** The Fourier transform of smooth functions that exponentially decay are smooth functions that exponentially decay. When we only have algebraic decay, that indicates that the Fourier transform is smooth apart from at $0$. "
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"9.594547378810603e-13"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a = 2.0\n",
"f = Fun(x -> 1/(x^2 + a^2), Γ)\n",
"f̂ = s -> fourier(f, -s)\n",
"f̂(0) - π/a"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.5959455978986625e-16 + 0.0im"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f̂(2.0) - (π*exp(-2.0*a)/a)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.5959455978986625e-16 + 0.0im"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f̂(-2.0) - (π*exp(-2.0*a)/a)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So we can now solve the integral equation\n",
"$$\n",
"\\lambda u(x) + \\int_{-\\infty}^\\infty K(x-t)u(t) \\dt = f(x)\\qqfor -\\infty < x < \\infty.\n",
"$$\n",
"by\n",
"converting to Fourier space\n",
"$$\n",
"(\\lambda + \\hat K(s)) \\hat u(s) = \\hat f(s)\\qqfor -\\infty < x < \\infty.\n",
"$$\n",
"to get\n",
"$$\n",
"\\hat u(s) = {\\pi \\over a} { \\E^{-a |s|} (\\gamma^2+s^2) \\over \\lambda (\\gamma^2+s^2) + 2\\gamma} \n",
"$$\n",
"Therefore, $u(x)$ is the inverse Fourier transform of this. We calculate it numerically first:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"λ = 1.0\n",
"û = π/a*Fun(s -> exp(-a*abs(s))*(γ^2 + s^2)/(λ*(γ^2 + s^2) + 2γ), Γ)\n",
"\n",
"u = Fun(x -> fourier(û, x), Γ)/(2π)\n",
"\n",
"plot(-10.0:0.1:10.0, u)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.894040480010517e-13 + 6.805472279452433e-19im"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"λ*u(0.0) + sum(u*K) - f(0.0)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.3935020632137594e-12 + 1.5832361733373612e-15im"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x =2.0\n",
"\n",
"λ*u(x) + sum(Fun(t -> K(x-t)*u(t), UnionDomain(-Inf .. x, x .. Inf))) - f(x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Working out the analytic formula for this example takes a lot more work and involves other special functions, and is outside the scope of the course."
]
}
],
"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
}