{
"cells": [
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"using Plots, ComplexPhasePortrait, ApproxFun, SingularIntegralEquations, DifferentialEquations\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",
"<br>\n",
"s.olver@imperial.ac.uk\n",
"\n",
"<br>\n",
"Website: https://github.com/dlfivefifty/M3M6LectureNotes\n",
"\n",
"\n",
"# Lecture 18: Orthogonal polynomials and singular integrals\n",
"\n",
"\n",
"This lecture we do the following:\n",
"\n",
"1. Cauchy transforms of weighted orthogonal polynomials\n",
" - Three-term recurrence and calculation\n",
" - Hilbert transform of weighted orthogonal polynomials\n",
" - Hilbert transform of weighted Chebyshev polynomials\n",
" \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Cauchy transforms of orthogonal polynomials\n",
"\n",
"Given a family of orthogonal polynomials $p_k(x)$ with respect to the weight $w(x)$ on $(a,b)$, we always know it satisfies a three-term recurrence:\n",
"\\begin{align*}\n",
"x p_0(x) &= a_0 p_0(x) + b_0 p_1(x) \\\\\n",
"x p_k(x) &= c_k p_{k-1}(x) + a_k p_k(x) + b_k p_{k+1}(x)\n",
"\\end{align*}\n",
"\n",
"Consider now the Cauchy transform of the weighted orthogonal polynomial:\n",
"$$\n",
"Q_k(z) := {\\cal C}_{[a,b]}[p_k w](z) = {1 \\over 2 \\pi \\I} \\int_a^b {p_k(x) w(x) \\over x -z} \\dx\n",
"$$\n",
"\n",
"\n",
"**Theorem (Three-term recurrence Cauchy transform of weighted OPs)** \n",
"$C_k(z)$ satisfies the same recurrence relationship as $p_k(x)$ for $k=1,2,\\ldots$:\n",
"\\begin{align*}\n",
"z C_0(z) &= a_0 C_0(z) + b_0 C_1(z) - {1 \\over 2 \\pi \\I} \\int_a^b w(x) \\dx \\\\\n",
"z C_k(z) &= c_k C_{k-1}(z) + a_k C_k(z) + b_k C_{k+1}(z)\n",
"\\end{align*}\n",
"\n",
"**Proof**\n",
"\\begin{align*}\n",
"z C_k(z) &= {1 \\over 2 \\pi \\I} \\int_a^b {z p_k(x) w(x) \\over x -z} \\dx = {1 \\over 2 \\pi \\I} \\int_a^b {(z -x) p_k(x) w(x) \\over x -z} \\dx + \\int_a^b {x p_k(x) w(x) \\over x -z} \\dx \\\\\n",
" &= -{1 \\over 2 \\pi \\I} \\int_a^b p_k(x) w(x) \\dx + \\int_a^b {(c_k p_{k-1}(x) + a_k p_k(x) + b_k p_{k+1}(x) w(x) \\over x -z} \\dx \\\\\n",
" &= -{1 \\over 2 \\pi \\I} \\int_a^b p_k(x) w(x) \\dx + c_k C_{k-1}(z) + a_k C_k(z) + b_k C_{k+1}(z)\n",
"\\end{align*}\n",
"when $k > 0$, the integral term disappears. \n",
"⬛️\n",
"\n",
"This gives a convenient way to calculate the Cauchy transforms: if we know $C_0(z) ={\\cal C}w(z)$ and $\\int_a^b w(x) \\dx$, solve the lower triangular system:\n",
"$$\n",
"\\begin{pmatrix}\n",
"1 \\\\\n",
"a_0-z & b_0 \\\\\n",
"c_1 & a_1-z & b_1 \\\\\n",
"&c_2 & a_2-z & b_2 \\\\\n",
"&& c_3 & a_3-z &\\ddots\\\\\n",
"&&&\\ddots & \\ddots\n",
"\\end{pmatrix}\\begin{pmatrix}C_0(z) \\\\C_1(z) \\\\C_2(z) \\\\C_3(z) \\\\ \\vdots \\end{pmatrix} = \\begin{pmatrix}C_0(z) \\\\{1 \\over 2 \\pi \\I} \\int_a^b w(x) \\dx \\\\0 \\\\0 \\\\ \\vdots \\end{pmatrix} \n",
"$$\n",
"\n",
"**Example (Chebyshev Cauchy transform)** \n",
"\n",
"Consider the Chebyshev case $w(x) = {1 \\over \\sqrt{1-x^2}}$, which satisfies $\\int_{-1}^1 w(x) \\dx = {\\pi}$. Recall that\n",
"$$\n",
" C_0(z) ={\\cal C}w(z) = { \\I \\over 2\\sqrt{z-1}\\sqrt{z+1}}\n",
"$$\n",
"Further, we have\n",
"\\begin{align*}\n",
"x T_0(x) = T_1(x) \\\\\n",
"x T_k(x) = {T_{k-1}(x) \\over 2} + { T_{k+1}(x) \\over 2} \n",
"\\end{align*}\n",
"hence\n",
"\\begin{align*}\n",
"z C_0(z) = C_1(z) - {1 \\over 2 \\I} \\\\\n",
"z C_k(z) = {C_{k-1}(z) \\over 2} +{C_{k+1}(z) \\over 2} .\n",
"\\end{align*}\n",
"In other words, we want to solve\n",
"$$\n",
"\\begin{pmatrix}\n",
"1 \\\\\n",
"-z & 1 \\\\\n",
"1/2 & -z & 1/2 \\\\\n",
"&1/2 & -z & 1/2 \\\\\n",
"&& 1/2 & -z &\\ddots\\\\\n",
"&&&\\ddots & \\ddots\n",
"\\end{pmatrix}\\begin{pmatrix}C_0(z) \\\\C_1(z) \\\\C_2(z) \\\\C_3(z) \\\\ \\vdots \\end{pmatrix} = \\begin{pmatrix} { \\I \\over 2\\sqrt{z-1}\\sqrt{z+1}} \\\\{1 \\over 2 \\I} \\\\0 \\\\0 \\\\ \\vdots \\end{pmatrix} \n",
"$$\n",
"with forward substitution.\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"10-element Array{Complex{Float64},1}:\n",
" 0.4999250218677838 + 0.004998750393615982im\n",
" 0.049492627147416784 - 0.44950762277386im \n",
" -0.40012497188352847 - 0.08500174951890464im \n",
" -0.11251727162034156 + 0.35248227849337344im \n",
" 0.3071250618607855 + 0.13299475089351104im \n",
" 0.14734333381379644 - 0.2644583159425141im \n",
" -0.22476473190952334 - 0.15641774731925456im \n",
" -0.16101273073185018 + 0.18822182009675853im \n",
" 0.1549178217438016 + 0.16185956519223624im \n",
" 0.15962438204216325 - 0.12486634270955096im "
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = Fun()\n",
"w = 1/sqrt(1-x^2)\n",
"z = 0.1+0.1im\n",
"\n",
"n = 10\n",
"L = zeros(ComplexF64,n,n)\n",
"L[1,1] = 1\n",
"L[2,1] = -z\n",
"L[2,2] = 1\n",
"for k=3:n\n",
" L[k,k-1] = -z\n",
" L[k,k-2] = L[k,k] = 1/2\n",
"end\n",
"\n",
"\n",
"\n",
"C = L \\ [ im/(2sqrt(z-1)sqrt(z+1)); 1/(2im); zeros(n-2)]"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-5.551115123125783e-17 + 5.551115123125783e-17im"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"T₅ = Fun(Chebyshev(), [zeros(5);1])\n",
"cauchy(T₅*w,z) - C[6]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Warning** This fails for large $n$ or large $z$:"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(-2.2332625421658917e6 + 521568.0612707699im, 0.0 + 8.834874115176436e-18im)"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = Fun()\n",
"w = 1/sqrt(1-x^2)\n",
"z = 5+6im\n",
"\n",
"n = 100\n",
"L = zeros(ComplexF64,n,n)\n",
"L[1,1] = 1\n",
"L[2,1] = -z\n",
"L[2,2] = 1\n",
"for k=3:n\n",
" L[k,k-1] = -z\n",
" L[k,k-2] = L[k,k] = 1/2\n",
"end\n",
"\n",
"C = L \\ [ im/(2sqrt(z-1)sqrt(z+1)); 1/(2im); zeros(n-2)]\n",
"\n",
"T₂₀ = Fun(Chebyshev(), [zeros(20);1])\n",
"\n",
"C[21], cauchy(T₂₀*w, z)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Get around it by dropping the first row:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"99×99 Array{Complex{Float64},2}:\n",
" -5.0-6.0im 1.0+0.0im 0.0+0.0im … 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.5+0.0im -5.0-6.0im 0.5+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.5+0.0im -5.0-6.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.5+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im … 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im … 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" ⋮ ⋱ \n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im … 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im … 0.5+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im -5.0-6.0im 0.5+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.5+0.0im -5.0-6.0im 0.5+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.5+0.0im -5.0-6.0im"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"L[2:end,1:end-1]"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-3.072062106493749e-17 + 4.410020129853371e-17im"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"C = L[2:end,1:end-1]\\ [1/(2im); zeros(n-2)]\n",
"\n",
"C[6]- cauchy(T₅*w, z)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Hilbert transform of weighted orthogonal polynomials\n",
"\n",
"Now consider the Hilbert transform of weighted orthogonal polynomials:\n",
"$$\n",
"H_k(x) = \\HH_{[a,b]}[p_k w](x) = {1 \\over \\pi} \\int_a^b {p_k(t) w(t) \\over t-x} \\dt\n",
"$$\n",
"\n",
"Just like Cauchy transforms, the Hilbert transforms have \n",
"\n",
"**Corollary (Hilbert transform recurrence)**\n",
"\\begin{align*}\n",
"x H_0(x) &= a_0 H_0(x) + b_0 H_1(x) - {1 \\over \\pi} \\int_a^b w(x) \\dx\\\\\n",
"x H_k(x) &= c_k H_{k-1}(x) + a_k H_k(x) + b_k H_{k+1}(x) \n",
"\\end{align*}\n",
"\n",
"**Proof**\n",
"Recall\n",
"$$\n",
"\\CC^+ f(x) + \\CC^- f(x) = -\\I \\HH f(x)\n",
"$$\n",
"Therefore, we have\n",
"$$\n",
"C_k^+(x) + C_k^-(x) = -\\I \\HH[w p_k](x)\n",
"$$\n",
"hence we have\n",
"\\begin{align*}\n",
"x H_0(x) &= \\I x (C_0^+(x) + C_0^-(x)) = \\I \\left[a_0 (C_0^+(x) + C_0^-(x)) + b_0 (C_1^+(x) + C_1^-(x)) -{1 \\over \\pi \\I} \\int_a^b w(x) \\dx \\right]\\\\\n",
" &= a_0 H_0(x) + b_0 H_1(x) - {1 \\over \\pi} \\int_a^b w(x) \\dx\n",
"\\end{align*}\n",
"Other $k$ follows by a similar argument.\n",
"\n",
"⬛️\n",
"\n",
"\n",
"\n",
"### Example: weighted Chebyshev\n",
"For\n",
"$$\n",
"H_k(x) = \\int_{-1}^1 {T_k(t) \\over (t-x)\\sqrt{1-t^2}} \\dt\n",
"$$\n",
"The recurrence gives us\n",
"\\begin{align*}\n",
"x H_0(x) &= H_1(x) -1 \\\\\n",
"x H_k(x) &= {H_{k-1}(x) \\over 2} + {H_k(x) \\over 2} \\\\\n",
"\\end{align*}\n",
"In this case, we have $H_0(x) = \\HH[w](x) = 0 $. Therefore, we can rewrite this recurrence as\n",
"\\begin{align*}\n",
" H_1(x)&= 1 \\\\\n",
" x H_1(x) &= {H_2(x) \\over 2} \\\\\n",
"x H_k(x) &= {H_{k-1}(x) \\over 2} + {H_k(x) \\over 2} \\\\\n",
"\\end{align*}\n",
"This is precisely the three-term recurrence satisfied by $U_{k-1}$! We therefore have\n",
"$$\n",
"H_k(x) = U_{k-1}(x)\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = 0.1\n",
"\n",
"T = Fun(Chebyshev(),[zeros(n);1])\n",
"hilbert(w*T,x) - Fun(Ultraspherical(1), [zeros(n-1);1])(x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This gives a very easy way to compute Hilbert transforms: if\n",
"$$\n",
"f(x) = \\sum_{k=0}^\\infty f_k T_k(x)\n",
"$$\n",
"then\n",
"$$\n",
"\\HH\\left[{f \\over \\sqrt{1-\\diamond^2}}\\right](x) = \\sum_{k=0}^\\infty f_{k+1} U_k(x)\n",
"$$"
]
}
],
"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.0.2"
}
},
"nbformat": 4,
"nbformat_minor": 2
}