Lecture 21

README.md

@@ -48,4 +48,5 @@ Examples of previous projects:
 17. [Logarithmic singular integrals](notes/Lecture17.pdf)
 18. [Logarithmic singular integral examples](notes/Lecture18.pdf)
 19. [Inverting logarithmic singular integrals](notes/Lecture19.pdf)
-20. [Orthogonal polynomials](notes/Lecture20.pdf)
+20. [Orthogonal polynomials](notes/Lecture20.pdf)
+21. [Classical orthogonal polynomials](notes/Lecture21.pdf)

src/Lecture21.jmd

@@ -9,37 +9,34 @@ [email protected]
 # Lecture 21: Classical orthogonal polynomials
 
 
-We will also investigate the properties of _classical_ OPs:
-
-1. Explicit formulae
-3. Rodriguez formulae
-2. Ordinary differential equations
-2. Derivatives
-
+We will also investigate the properties of _classical_ OPs.
 A good reference is  [Digital Library of Mathematical Functions, Chapter 18](http://dlmf.nist.gov/18).
 
-2. Definition of classical orthogonal polynomials
-    - Hermite, Laguerre, and Jacobi polynomials
-    - Legendre, Chebyshev, and ultraspherical polynomials
-    - Explicit construction for Chebyshev polynomials
+This lecture we discuss
+1. Hermite, Laguerre, and Jacobi polynomials
+2. Legendre, Chebyshev, and ultraspherical polynomials
+3. Explicit construction for Chebyshev polynomials
 
 
 ## Definition of classical orthogonal polynomials
 
 Classical orthogonal polynomials are orthogonal with respect to the following three weights:
 
 
-| Name        | Interval $(a,b)$ |Weight  function  $w(x)$      | Standard polynomial   | highest order coefficient $k_n$ |
+| Name        |  $(a,b)$ |  $w(x)$      | Notation   |  $k_n$ |
 |:-------------|:------------- |:----------------------|:-----|:-----|
 | Hermite     |$(-\infty,\infty)$ | $\E^{-x^2}$             | $H_n(x)$ | $2^n$ |
-| Laguerre    | $(0,\infty)$ | $x^\alpha \E^{-x}$      |   $L_n^{(\alpha)}(x)$ | See [Table 18.3.1](http://dlmf.nist.gov/18.3) |
-| Jacobi      | $(-1,1)$ | $(1-x)^{\alpha} (1+x)^\beta$      |    $P_n^{(\alpha,\beta)}(x)$ | See [Table 18.3.1](http://dlmf.nist.gov/18.3) |
+| Laguerre    | $(0,\infty)$ | $x^\alpha \E^{-x}$      |   $L_n^{(\alpha)}(x)$ | [Table 18.3.1](http://dlmf.nist.gov/18.3) |
+| Jacobi      | $(-1,1)$ | $(1-x)^{\alpha} (1+x)^\beta$      |    $P_n^{(\alpha,\beta)}(x)$ | [Table 18.3.1](http://dlmf.nist.gov/18.3) |
+
 
 
-Note out of convention the parameters for Jacobi polynomials are in the "wrong" order.
 
-We can actually construct these polynomials in Julia: first consider Hermite:
+Note out of convention the parameters for Jacobi polynomials are right-to-left order.
 
+We can actually construct these polynomials in Julia, first consider Hermite:
+```julia
+using ApproxFun, Plots, LinearAlgebra, ComplexPhasePortrait
 H₀ = Fun(Hermite(), [1])
 H₁ = Fun(Hermite(), [0,1])
 H₂ = Fun(Hermite(), [0,0,1])
@@ -54,17 +51,17 @@ plot!(xx, H₂.(xx); label="H_2", ylims=(-400,400))
 plot!(xx, H₃.(xx); label="H_3", ylims=(-400,400))
 plot!(xx, H₄.(xx); label="H_4", ylims=(-400,400))
 plot!(xx, H₅.(xx); label="H_5")
-
+```
 
 We verify their orthogonality:
-
+```julia
 w = Fun(GaussWeight(), [1.0])
 
 @show sum(H₂*H₅*w)  # means integrate
 @show sum(H₅*H₅*w);
-
+```
 Now Jacobi:
-
+```julia
 α,β = 0.1,0.2
 P₀ = Fun(Jacobi(β,α), [1])
 P₁ = Fun(Jacobi(β,α), [0,1])
@@ -84,28 +81,29 @@ plot!(xx, P₅.(xx); label="P_5^($α,$β)")
 w = Fun(JacobiWeight(β,α), [1.0])
 @show sum(P₂*P₅*w)  # means integrate
 @show sum(P₅*P₅*w);
+```
 
 ## Legendre, Chebyshev, and ultraspherical polynomials
 
 There are special families of Jacobi weights with their own name. 
 
-| Name        | Jacobi parameters  |Weight  function  $w(x)$      | Standard polynomial   | highest order coefficient $k_n$ |
+| Name        | Jacobi parameters  |  $w(x)$      | Notation   |  $k_n$ |
 |:-------------|:------------- |:----------------------|:-----|:------|
-| Jacobi      | $\alpha,\beta$ | $(1-x)^{\alpha} (1+x)^\beta$      |    $P_n^{(\alpha,\beta)}(x)$ | See [Table 18.3.1](http://dlmf.nist.gov/18.3) | 
+| Jacobi      | $\alpha,\beta$ | $(1-x)^{\alpha} (1+x)^\beta$      |    $P_n^{(\alpha,\beta)}(x)$ | [Table 18.3.1](http://dlmf.nist.gov/18.3) | 
 | Legendre      | $0,0$ | $1$      |    $P_n(x)$ | $2^n(1/2)_n/n!$ | 
-| Chebyshev (first kind)      | $-{1 \over 2},-{1 \over 2}$ | $1 \over \sqrt{1-x^2}$ | $T_n(x)$ | $1 (n=0), 2^{n-1} (n \neq 0)$ |
-| Chebyshev (second kind)      | ${1 \over 2},{1 \over 2}$ | $\sqrt{1-x^2}$      |    $U_n(x)$ | $2^n$
+| Chebyshev (1st)      | $-{1 \over 2},-{1 \over 2}$ | $1 \over \sqrt{1-x^2}$ | $T_n(x)$ | $1 (n=0), 2^{n-1} (n \neq 0)$ |
+| Chebyshev (2nd)      | ${1 \over 2},{1 \over 2}$ | $\sqrt{1-x^2}$      |    $U_n(x)$ | $2^n$
 | Ultraspherical      | $\lambda-{1 \over 2},\lambda-{1 \over 2}$ | $(1-x^2)^{\lambda - 1/2}, \lambda \neq 0$      |    $C_n^{(\lambda)}(x)$ | $2^n(\lambda)_n/n!$ |
 
 Note that other than Legendre, these polynomials have a different normalization than $P_n^{(\alpha,\beta)}$:
-
-
+```julia
 T₂ = Fun(Chebyshev(), [0.0,0,1])
 P₂ = Fun(Jacobi(-1/2,-1/2), [0.0,0,1])
 plot(T₂; label="T_2", title="T_2 is C*P_2 for some C")
 plot!(P₂; label="P_2") 
-
-But because they are orthogonal w.r.t. the same weight, they must be a constant multiple of each-other.
+```
+But because they are orthogonal w.r.t. the same weight, they must be a constant multiple of each-other,
+as discussed last lecture.
 
 ### Explicit construction of Chebyshev polynomials (first kind and second kind)
 
@@ -120,12 +118,16 @@ $$
 
 **Proof** We first show that they are orthogonal w.r.t. $1/\sqrt{1-x^2}$. Too easy: do $x = \cos \theta$, $\dx = -\sin \theta$ to get (for $n \neq m$)
 $$
-    \int_{-1}^1 {\cos n \acos x \cos m \acos x \dx \over \sqrt{1-x^2}} = -\int_\pi^0  \cos n \theta \cos m \theta \D \theta =  \int_0^\pi  {\E^{\I (-n-m)\theta} + \E^{\I (n-m)\theta} + \E^{\I (m-n)\theta} + \E^{\I (n+m)\theta}    \over 4} \D \theta =0
+\begin{align*}
+    \int_{-1}^1 {\cos n \acos x \cos m \acos x \dx \over \sqrt{1-x^2}} &= -\int_\pi^0  \cos n \theta \cos m \theta \D \theta \ccr
+    =  \int_0^\pi  {\E^{\I (-n-m)\theta} + \E^{\I (n-m)\theta} + \E^{\I (m-n)\theta} + \E^{\I (n+m)\theta}    \over 4} \D \theta =0
+\end{align*}
 $$
 
 We then need to show it has the right highest order term $k_n$. Note that $k_0 = k_1 = 1$.  Using $z = \E^{\I \theta}$ we see that $\cos n \theta$ has a simple recurrence for $n=2,3,\ldots$:
 $$
-\cos n \theta = {z^n + z^{-n} \over 2} = 2 {z + z^{-1} \over 2} {z^{n-1} + z^{1-n} \over 2}- {z^{n-2} + z^{2-n} \over 2} = 2 \cos \theta \cos (n-1)\theta - \cos(n-2)\theta 
+\cos n \theta = {z^n + z^{-n} \over 2} = 2 {z + z^{-1} \over 2} {z^{n-1} + z^{1-n} \over 2}- {z^{n-2} + z^{2-n} \over 2} = 
+2 \cos \theta \cos (n-1)\theta - \cos(n-2)\theta 
 $$
 thus 
 $$
@@ -135,10 +137,9 @@ It follows that
 $$
 k_n = 2 x k_{n-1} = 2^{n-1} k_1 = 2^{n-1}
 $$
-By uniqueness we have $T_n(x) \cos n \acos x$.
-
-⬛️
+By uniqueness we have $T_n(x) = \cos n \acos x$.
 
+■
 
 
 **Proposition (Chebyshev second kind formula)** 
@@ -170,8 +171,10 @@ J^\top \vc f = \begin{pmatrix}
 &&&\ddots & \ddots
 \end{pmatrix} \begin{pmatrix} f_0\\ f_1\\f_2\\f_3\\\vdots\end{pmatrix}
 $$
-In the case where $f$ is a degree $n-1$  polynomial, we can represent $J^\top$ as an $n+1 \times n$ matrix (this makes sense as $x f(x)$ is one more degree than $f$):
 
+### Demonstration
+In the case where $f$ is a degree $n-1$  polynomial, we can represent $J^\top$ as an $n+1 \times n$ matrix (this makes sense as $x f(x)$ is one more degree than $f$):
+```julia
 f = Fun(exp, Chebyshev())
 n = ncoefficients(f) # number of coefficients
 @show n
@@ -181,25 +184,25 @@ for k=2:n
     J[k,k-1] = J[k,k+1] = 1/2
 end
 J'
-
+```
+```julia
 cfs = J'*f.coefficients # coefficients of x*f
 xf = Fun(Chebyshev(), cfs)
 
 xf(0.1) - 0.1*f(0.1)
-
-
-**Example** We can construct $T_0(x),\ldots,T_{n-1}(x)$ via
+```
+We can construct $T_0(x),\ldots,T_{n-1}(x)$ via
+$$
 \begin{align*}
-    p_0(x) &= 1\\
-    p_1(x) &= x T_0(x) \\
+    T_0(x) &= 1\\
+    T_1(x) &= x T_0(x) \\
     T_2(x) &= 2x  T_0(x) -  T_0(x) \\    
     T_3(x) &= 2x T_1(x) - T_1(x)
     &\vdots
 \end{align*}
+$$
 Believe it or not, this is much faster than using $\cos k \acos x$:
-
-
-
+```julia
 function recurrence_Chebyshev(n,x)
     T = zeros(n)
     T[1] = 1.0
@@ -214,14 +217,14 @@ trig_Chebyshev(n,x) = [cos(k*acos(x)) for k=0:n-1]
 
 n = 10
 recurrence_Chebyshev(n, 0.1) - trig_Chebyshev(n,0.1) |>norm
-
+```
+```julia
 n = 10000
 @time recurrence_Chebyshev(n, 0.1) 
 @time trig_Chebyshev(n,0.1);
+```
 
-
-
-*Example* For Chebyshev, we want to solve the system
+We can also demonstrate Clenshaw's algorithm for evaluating polynomials. For Chebyshev, we want to solve the system
 $$
 \underbrace{\begin{pmatrix}
 1 & -x & \half \\
@@ -233,7 +236,7 @@ $$
 \end{pmatrix}}_{L_x^\top} \begin{pmatrix} \gamma_0 \\\vdots\\ \gamma_{n-1} \end{pmatrix}
 $$
 via
-
+$$
 \begin{align*}
 \gamma_{n-1} &= 2f_{n-1} \\
 \gamma_{n-2} &= 2f_{n-2} + 2x \gamma_{n-1} \\
@@ -242,9 +245,9 @@ via
 \gamma_1 &= f_1 + x \gamma_2 - \half \gamma_3 \\
 \gamma_0 &= f_0 + x \gamma_1 - \half \gamma_2
 \end{align*}
-
+$$
 then $f(x) = \gamma_0$.
-
+```julia
 function clenshaw_Chebyshev(f,x)
     n = length(f)
     γ = zeros(n)
@@ -260,19 +263,18 @@ end
 
 f = Fun(exp, Chebyshev())
 clenshaw_Chebyshev(f.coefficients, 0.1) - exp(0.1)
+```
 
-With some high performance computing tweeks, this can be made more accurate: this is the algorithm used for evaluating functions in ApproxFun:
-
+With some high performance computing tweeks, this can be made more accurate.
+This is the algorithm used for evaluating functions in ApproxFun:
+```julia
 f(0.1) - exp(0.1)
-
+```
 
 
 ## Approximation with Chebyshev polynomials
 
-
-
-
-*Example* Last lecture, we used the formula, derived via trigonometric manipulations, 
+Last lecture, we used the formula, derived via trigonometric manipulations, 
 $$
 T_1(x) = x T_0(x) \\
 T_{n+1}(x) = 2x T_n(x) - T_{n-1}(x)
@@ -284,82 +286,57 @@ x T_n(x)  =  {T_{n-1}(x) \over 2} + {T_{n+1}(x) \over 2}
 $$
 This tells us that we have the three-term recurrence with $a_n = 0$, $b_0 = 1$, $c_n = b_n = {1 \over 2}$ for $n > 0$.
 
-T = (n,x) -> cos(n*acos(x))
-x = 0.5
-n = 10
-@show x*T(0,x) - (T(1,x))
-@show x*T(n,x) - (T(n-1,x) + T(n+1,x))/2;
-
-
-Here, we demonstrate this with Chebyshev polynomials:
-```julia
-f = Fun(x -> 1 + x + x^2 + x^3, Chebyshev())
-f₀, f₁, f₂, f₃ = f.coefficients
-```
-```julia
-x = 0.1
-@show f₀*1 + f₁*x + f₂*cos(2acos(x)) + f₃*cos(3acos(x))
-@show 1 + x + x^2 + x^3;
-```
-
-
-plot(Fun(exp))
-plot!(Fun(t-> exp(cos(t)), -pi .. pi))
-
-x = Fun()
-ip = (f,g) -> sum(f*g/sqrt(1-x^2))
-
-T₂ = cos.(2 .* acos.(x))
-ip(T₂,f)/ip(T₂,T₂) # gives back f₂
-
-Of course, if $p_k$ are othernormal than we don't need the denominator.
-
-This can be extended to function approximation Provided the sum converges absolutely and uniformly in $x$, we can write
+This can be extended to function approximation. Provided the sum converges absolutely and uniformly in $x$, we can write
 $$
-f(x) = \sum_{k=0}^\infty f_k p_k(x).
+f(x) = \sum_{k=0}^\infty f_k T_k(x).
 $$
 In practice, we can approximate smooth functions by a finite truncation:
 $$
-f(x) \approx f_n(x) = \sum_{k=0}^{n-1} f_k p_k(x)
+f(x) \approx f_n(x) = \sum_{k=0}^{n-1} T_k p_k(x)
 $$
 
 Here we see that $\E^x$ can be approximated by a Chebyshev approximation using 14 coefficients and is accurate to 16 digits:
-
+```julia
 f = Fun(x -> exp(x), Chebyshev())
 @show f.coefficients
 @show ncoefficients(f)
 
 @show f(0.1) # equivalent to f.coefficients'*[cos(k*acos(x)) for k=0:ncoefficients(f)-1]
 @show exp(0.1);
+```
 
-The accuracy of this approximation is typically dictated by the smoothness of $f$: the more times we can differentiate, the faster it converges. For analytic functions, it's dictated by the domain of analyticity, just like Laurent/Fourier series. In the case above, $\E^x$ is entire hence we get faster than exponential convergence.
-
+The accuracy of this approximation is typically dictated by the smoothness of $f$: the more times we can differentiate, the faster it converges. 
+For analytic functions, it's dictated by the domain of analyticity, just like Laurent/Fourier series. In the case above, $\E^x$ is entire hence we get faster than exponential convergence.
 
-Chebyshev expansions work even when Taylor series do not. For example, the following function has poles at $\pm {\I \over 5}$, which means the radius of convergence for the Taylor series is $ |x| < {1 \over 5}$, but Chebyshev polynomials continue to work on $[-1,1]$:
 
+Chebyshev expansions work even when Taylor series do not. For example, the following function has poles at $\pm {\I \over 5}$, which means the radius of convergence for the Taylor series is $|x| < {1 \over 5}$, 
+but Chebyshev polynomials continue to work on $[-1,1]$:
+```julia
 f = Fun( x -> 1/(25x^2 + 1), Chebyshev())
 @show ncoefficients(f)
 plot(f)
+```
 
 This can be explained for Chebyshev expansion by noting that it is cosine expansion / Fourier expansion of an even function:
 $$
 f(x) = \sum_{k=0}^\infty f_k T_k(x) \Leftrightarrow f(\cos \theta) = \sum_{k=0}^\infty f_k \cos k \theta
 $$
-From Lecture 4 we saw that Fourier coefficients decay exponentially (thence the approximation converges exponentially fast) uf  $f\left({z + z^{-1} \over 2}\right)$ is analytic in an ellipse. In the case of $f(x) = {1 \over 25 x^2 + 1}$, we find that 
+From Lecture 4 we saw that Fourier coefficients decay exponentially (thence the approximation converges exponentially fast) of  
+$f\left({z + z^{-1} \over 2}\right)$ is analytic in an ellipse. In the case of $f(x) = {1 \over 25 x^2 + 1}$, we find that 
 $$
 f(z) = {4 z^2 \over 25 + 54 z^2 + 25 z^4}
 $$
 which has poles at 
 $
 \pm 0.8198040\I,\pm1.2198\I:
 $
-
+```julia
 f = x -> 1/(25x^2 + 1)
 portrait(-3..3, -3..3, z -> f((z+1/z)/2))
-
+```
 Hence we predict a rate of decay of about $1.2198^{-k}$:
-
+```julia
 f = Fun( x -> 1/(25x^2 + 1), Chebyshev())
 plot(abs.(f.coefficients) .+ 1E-40; yscale=:log10, label="Chebyshev coefficients")
 plot!( 1.2198.^(-(0:ncoefficients(f))); label="R^(-k)")
-
+```

src/M3M6AppliedComplexAnalysis.jl

@@ -20,6 +20,7 @@ weave("src/Lecture17.jmd", doctype="md2tex", informat="markdown", out_path=pwd()
 weave("src/Lecture18.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
 weave("src/Lecture19.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
 weave("src/Lecture20.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
+weave("src/Lecture21.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
 
 weave("src/Solutions1.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")
 weave("src/Solutions2.jmd", doctype="md2tex", informat="markdown", out_path=pwd()*"/output/", template="src/template.tpl")

src/template.tpl

@@ -37,6 +37,7 @@
 \def\qqand{\qquad\hbox{and}\qquad}
 \def\qqfor{\qquad\hbox{for}\qquad}
 \def\qqas{\qquad\hbox{as}\qquad}
+\def\half{ {1 \over 2} }
 \def\D{ {\rm d} }
 \def\I{ {\rm i} }
 \def\E{ {\rm e} }