Notebook

In [2]:

using Plots, ComplexPhasePortrait, ApproxFun, SingularIntegralEquations, DifferentialEquations
gr();

M3M6: Methods of Mathematical Physics¶

$\def\dashint{{\int\!\!\!\!\!\!-\,}} \def\infdashint{\dashint_{\!\!\!-\infty}^{\,\infty}} \def\D{\,{\rm d}} \def\E{{\rm e}} \def\dx{\D x} \def\dt{\D t} \def\dz{\D z} \def\C{{\mathbb C}} \def\R{{\mathbb R}} \def\CC{{\cal C}} \def\HH{{\cal H}} \def\I{{\rm i}} \def\qqqquad{\qquad\qquad} \def\qqand{\qquad\hbox{for}\qquad} \def\qqfor{\qquad\hbox{for}\qquad} \def\qqwhere{\qquad\hbox{where}\qquad} \def\Res_#1{\underset{#1}{\rm Res}}\, \def\sech{{\rm sech}\,} \def\acos{\,{\rm acos}\,} \def\vc#1{{\mathbf #1}} \def\ip<#1,#2>{\left\langle#1,#2\right\rangle} \def\norm#1{\left\|#1\right\|} \def\half{{1 \over 2}}$

Dr. Sheehan Olver

s.olver@imperial.ac.uk
Website: https://github.com/dlfivefifty/M3M6LectureNotes

Lecture 18: Orthogonal polynomials and singular integrals¶

This lecture we do the following:

Cauchy transforms of weighted orthogonal polynomials
- Three-term recurrence and calculation
- Hilbert transform of weighted orthogonal polynomials
- Hilbert transform of weighted Chebyshev polynomials

Cauchy transforms of orthogonal polynomials¶

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:

$\begin{align*} x p_0(x) &= a_0 p_0(x) + b_0 p_1(x) \\ x p_k(x) &= c_k p_{k-1}(x) + a_k p_k(x) + b_k p_{k+1}(x) \end{align*}$

Consider now the Cauchy transform of the weighted orthogonal polynomial:

$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$

Theorem (Three-term recurrence Cauchy transform of weighted OPs) $C_k(z)$ satisfies the same recurrence relationship as $p_k(x)$ for $k=1,2,\ldots$ :

$\begin{align*} 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 \\ z C_k(z) &= c_k C_{k-1}(z) + a_k C_k(z) + b_k C_{k+1}(z) \end{align*}$

Proof

$\begin{align*} 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 \\ &= -{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 \\ &= -{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) \end{align*}$ when

$k > 0$ , the integral term disappears.
⬛️

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:

$\begin{pmatrix} 1 \\ a_0-z & b_0 \\ c_1 & a_1-z & b_1 \\ &c_2 & a_2-z & b_2 \\ && c_3 & a_3-z &\ddots\\ &&&\ddots & \ddots \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}$

Example (Chebyshev Cauchy transform)

Consider the Chebyshev case $w(x) = {1 \over \sqrt{1-x^2}}$ , which satisfies $\int_{-1}^1 w(x) \dx = {\pi}$ . Recall that

$C_0(z) ={\cal C}w(z) = { \I \over 2\sqrt{z-1}\sqrt{z+1}}$ Further, we have

$\begin{align*} x T_0(x) = T_1(x) \\ x T_k(x) = {T_{k-1}(x) \over 2} + { T_{k+1}(x) \over 2} \end{align*}$ hence

$\begin{align*} z C_0(z) = C_1(z) - {1 \over 2 \I} \\ z C_k(z) = {C_{k-1}(z) \over 2} +{C_{k+1}(z) \over 2} . \end{align*}$ In other words, we want to solve

$\begin{pmatrix} 1 \\ -z & 1 \\ 1/2 & -z & 1/2 \\ &1/2 & -z & 1/2 \\ && 1/2 & -z &\ddots\\ &&&\ddots & \ddots \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}$ with forward substitution.

In [5]:

x = Fun()
w = 1/sqrt(1-x^2)
z = 0.1+0.1im

n = 10
L = zeros(ComplexF64,n,n)
L[1,1] = 1
L[2,1] = -z
L[2,2] = 1
for k=3:n
    L[k,k-1] = -z
    L[k,k-2] = L[k,k] = 1/2
end



C = L \ [ im/(2sqrt(z-1)sqrt(z+1)); 1/(2im); zeros(n-2)]

Out[5]:

10-element Array{Complex{Float64},1}:
   0.4999250218677838 + 0.004998750393615982im
 0.049492627147416784 - 0.44950762277386im    
 -0.40012497188352847 - 0.08500174951890464im 
 -0.11251727162034156 + 0.35248227849337344im 
   0.3071250618607855 + 0.13299475089351104im 
  0.14734333381379644 - 0.2644583159425141im  
 -0.22476473190952334 - 0.15641774731925456im 
 -0.16101273073185018 + 0.18822182009675853im 
   0.1549178217438016 + 0.16185956519223624im 
  0.15962438204216325 - 0.12486634270955096im

In [6]:

T₅ = Fun(Chebyshev(), [zeros(5);1])
cauchy(T₅*w,z) - C[6]

Out[6]:

-5.551115123125783e-17 + 5.551115123125783e-17im

Warning This fails for large $n$ or large $z$ :

In [11]:

x = Fun()
w = 1/sqrt(1-x^2)
z = 5+6im

n = 100
L = zeros(ComplexF64,n,n)
L[1,1] = 1
L[2,1] = -z
L[2,2] = 1
for k=3:n
    L[k,k-1] = -z
    L[k,k-2] = L[k,k] = 1/2
end

C = L \ [ im/(2sqrt(z-1)sqrt(z+1)); 1/(2im); zeros(n-2)]

T₂₀ = Fun(Chebyshev(), [zeros(20);1])

C[21], cauchy(T₂₀*w, z)

Out[11]:

(-2.2332625421658917e6 + 521568.0612707699im, 0.0 + 8.834874115176436e-18im)

Get around it by dropping the first row:

In [12]:

L[2:end,1:end-1]

Out[12]:

99×99 Array{Complex{Float64},2}:
 -5.0-6.0im   1.0+0.0im   0.0+0.0im  …   0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.5+0.0im  -5.0-6.0im   0.5+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.5+0.0im  -5.0-6.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.5+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im  …   0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im  …   0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
     ⋮                               ⋱                                    
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im  …   0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im  …   0.5+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im     -5.0-6.0im   0.5+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.5+0.0im  -5.0-6.0im   0.5+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.5+0.0im  -5.0-6.0im

In [13]:

C = L[2:end,1:end-1]\ [1/(2im); zeros(n-2)]

C[6]- cauchy(T₅*w, z)

Out[13]:

-3.072062106493749e-17 + 4.410020129853371e-17im

Hilbert transform of weighted orthogonal polynomials¶

Now consider the Hilbert transform of weighted orthogonal polynomials:

$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$

Just like Cauchy transforms, the Hilbert transforms have

Corollary (Hilbert transform recurrence)

$\begin{align*} x H_0(x) &= a_0 H_0(x) + b_0 H_1(x) - {1 \over \pi} \int_a^b w(x) \dx\\ x H_k(x) &= c_k H_{k-1}(x) + a_k H_k(x) + b_k H_{k+1}(x) \end{align*}$

Proof Recall

$\CC^+ f(x) + \CC^- f(x) = -\I \HH f(x)$ Therefore, we have

$C_k^+(x) + C_k^-(x) = -\I \HH[w p_k](x)$ hence we have

$\begin{align*} 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]\\ &= a_0 H_0(x) + b_0 H_1(x) - {1 \over \pi} \int_a^b w(x) \dx \end{align*}$ Other

$k$ follows by a similar argument.

⬛️

Example: weighted Chebyshev¶

For

$H_k(x) = \int_{-1}^1 {T_k(t) \over (t-x)\sqrt{1-t^2}} \dt$ The recurrence gives us

$\begin{align*} x H_0(x) &= H_1(x) -1 \\ x H_k(x) &= {H_{k-1}(x) \over 2} + {H_k(x) \over 2} \\ \end{align*}$ In this case, we have

$H_0(x) = \HH[w](x) = 0$ . Therefore, we can rewrite this recurrence as

$\begin{align*} H_1(x)&= 1 \\ x H_1(x) &= {H_2(x) \over 2} \\ x H_k(x) &= {H_{k-1}(x) \over 2} + {H_k(x) \over 2} \\ \end{align*}$ This is precisely the three-term recurrence satisfied by

$U_{k-1}$ ! We therefore have

$H_k(x) = U_{k-1}(x)$

In [14]:

x = 0.1

T = Fun(Chebyshev(),[zeros(n);1])
hilbert(w*T,x) - Fun(Ultraspherical(1), [zeros(n-1);1])(x)

Out[14]:

0.0

This gives a very easy way to compute Hilbert transforms: if

$f(x) = \sum_{k=0}^\infty f_k T_k(x)$ then

$\HH\left[{f \over \sqrt{1-\diamond^2}}\right](x) = \sum_{k=0}^\infty f_{k+1} U_k(x)$