{
"cells": [
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"using Plots, ComplexPhasePortrait, ApproxFun, SingularIntegralEquations, DifferentialEquations, LinearAlgebra\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\\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",
"Website: https://github.com/dlfivefifty/M3M6LectureNotes\n",
"\n",
"\n",
"# Lecture 15: Recurrence relationships\n",
"\n",
"\n",
"This lecture we do the following:\n",
"1. Jacobi operators and three-term recurences for general orthogonal polynomials\n",
" - Three-term recurrence relationship\n",
" - Jacobi operator and multiplication by $x$ \n",
" - Evaluating polynomials and Clenshaw's algorithm\n",
"2. Gram–Schmidt, revisited\n",
"\n",
"\n",
"A central theme: if you know the Jacobi operator / three-term recurrence, you know the polynomials. This is the __best__ way to evaluate expansions in orthogonal polynomials: even for Chebyshev $T_n(x) = \\cos n \\acos x$, using the recurrence avoids evaluating trigonometric functions.\n",
"\n",
"## Jacobi operators and three-term recurences for general orthogonal polynomials\n",
"### Three-term recurrence relationships\n",
"\n",
"Every family of orthogonal polynomials has a three-term recurrence relationship:\n",
"\n",
"**Theorem (three-term recurrence)** Suppose $\\{p_n(x)\\}$ are a family of orthogonal polynomials w.r.t. a weight $w(x)$. Then there exists constants $a_n \\neq 0$, $b_n$ and $c_n$ such that\n",
"\\begin{align*}\n",
"x p_0(x) = a_0 p_0(x) + b_0 p_1(x) \\\\\n",
"x p_n(x) = c_n p_{n-1}(x) + a_n p_n(x) + b_n p_{n+1}(x)\n",
"\\end{align*}\n",
"\n",
"**Proof**\n",
"The first part follows since $p_0(x)$ and $p_1(x)$ span all degree 1 polynomials.\n",
"\n",
"The second part follows essentially because multiplication by $x$ is \"self-adjoint\", that is,\n",
"$$\n",
"\\ip = \\int_a^b x f(x) g(x) w(x) \\dx = \\ip\n",
"$$\n",
"Therefore, if $f_m$ is a degree $m < n-1$ polynomial, we have\n",
"$$\n",
"\\ip = \\ip = 0.\n",
"$$\n",
"In particular, if we write\n",
"$$\n",
"x p_n(x) = \\sum_{k=0}^{n+1} C_k p_k(x)\n",
"$$\n",
"then\n",
"$$\n",
"C_k = {\\ip< x p_n, p_k> \\over \\norm{p_k}^2} = 0\n",
"$$\n",
"if $k < n-1$.\n",
"\n",
"⬛️\n",
"\n",
"\n",
"Monic polynomials clearly have $b_n = 1$. Orthonormal polynomials have an even simpler form:\n",
"\n",
"**Theorem (orthonormal three-term recurrence)** Suppose $\\{q_n(x)\\}$ are a family of orthonotms polynomials w.r.t. a weight $w(x)$. Then there exists constants $a_n$ and $b_n$ such that\n",
"\\begin{align*}\n",
"x q_0(x) = a_0 q_0(x) + b_0 q_1(x)\\\\\n",
"x q_n(x) = b_{n-1} q_{n-1}(x) + a_n q_n(x) + b_{n} q_{n+1}(x)\n",
"\\end{align*}\n",
"\n",
"**Proof**\n",
"Follows again by self-adjointness of multiplication by $x$:\n",
"$$\n",
"c_n = \\ip = \\ip = \\ip = b_{n-1}\n",
"$$\n",
"⬛️\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"*Example* Last lecture, we used the formula, derived via trigonometric manipulations, \n",
"$$\n",
"T_1(x) = x T_0(x) \\\\\n",
"T_{n+1}(x) = 2x T_n(x) - T_{n-1}(x)\n",
"$$\n",
"Rearranging, this becomes\n",
"$$\n",
" x T_0(x) = T_1(x) \\\\\n",
"x T_n(x) = {T_{n-1}(x) \\over 2} + {T_{n+1}(x) \\over 2}\n",
"$$\n",
"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$."
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x * T(0, x) - T(1, x) = 1.1102230246251565e-16\n",
"x * T(n, x) - (T(n - 1, x) + T(n + 1, x)) / 2 = 5.273559366969494e-16\n"
]
}
],
"source": [
"T = (n,x) -> cos(n*acos(x))\n",
"x = 0.5\n",
"n = 10\n",
"@show x*T(0,x) - (T(1,x))\n",
"@show x*T(n,x) - (T(n-1,x) + T(n+1,x))/2;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Corollary (symmetric three-term recurrence implies orthonormal)** Suppose $\\{p_n(x)\\}$ are a family of orthogonal polynomials w.r.t. a weight $w(x)$ such that\n",
"\\begin{align*}\n",
"x p_0(x) = a_0 p_0(x) + b_0 p_1(x)\\\\\n",
"x p_n(x) = b_{n-1} p_{n-1}(x) + a_n p_n(x) + b_{n} p_{n+1}(x)\n",
"\\end{align*}\n",
"with $b_n \\neq 0$. Suppose further that $\\norm{p_0} = 1$. Then $p_n$ must be orthonormal.\n",
"\n",
"**Proof** This follows from\n",
"$$\n",
"b_n = {\\ip \\over \\norm{p_{n+1}}^2} = {\\ip \\over \\norm{p_{n+1}}^2} = b_n {\\norm{p_n}^2 \\over \\norm{p_{n+1}}^2 }\n",
"$$\n",
"which implies\n",
"$$\n",
"\\norm{p_{n+1}}^2 = \\norm{p_n}^2 = \\cdots = \\norm{p_0}^2 = 1\n",
"$$\n",
"⬛️\n",
"\n",
"**Remark** We can scale $w(x)$ by a constant without changing the orthogonality properties, hence we can make $\\|p_0\\| = 1$ by changing the weight.\n",
"\n",
"**Remark** This is beyond the scope of this course, but satisfying a three-term recurrence like this such that coefficients are bounded with $p_0(x) = 1$ is sufficient to show that there exists a $w(x)$ (or more accurately, a Borel measure) such that $p_n(x)$ are orthogonal w.r.t. $w$. The relationship between the coefficients $a_n,b_n$ and the $w(x)$ is an object of study in spectral theory, particularly when the coefficients are periodic, quasi-periodic or random. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Jacobi operators and multiplication by $x$\n",
"\n",
"We can rewrite the three-term recurrence as\n",
"$$\n",
"x \\begin{pmatrix} p_0(x) \\cr p_1(x) \\cr p_2(x) \\cr \\vdots \\end{pmatrix} = J\\begin{pmatrix} p_0(x) \\cr p_1(x) \\cr p_2(x) \\cr \\vdots \\end{pmatrix}\n",
"$$\n",
"where $J$ is a Jacobi operator, an infinite-dimensional tridiagonal matrix:\n",
"$$\n",
"J = \\begin{pmatrix} \n",
"a_0 & b_0 \\cr\n",
"c_1 & a_1 & b_1 \\cr\n",
"& c_2 & a_2 & b_2 \\cr\n",
"&& c_3 & a_3 & \\ddots \\cr\n",
"&&&\\ddots & \\ddots\n",
"\\end{pmatrix} \n",
"$$\n",
"\n",
"When the polynomials are monic, we have $1$ on the superdiagonal. When we have an orthonormal basis, then $J$ is symmetric:\n",
"$$\n",
"J = \\begin{pmatrix} \n",
"a_0 & b_0 \\cr\n",
"b_0 & a_1 & b_1 \\cr\n",
"& b_1 & a_2 & b_2 \\cr\n",
"&& b_2 & a_3 & \\ddots \\cr\n",
"&&&\\ddots & \\ddots\n",
"\\end{pmatrix} \n",
"$$\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"Given a polynomial expansion, $J$ tells us how to multiply by $x$ in coefficient space, that is, if\n",
"$$\n",
"f(x) = \\sum_{k=0}^\\infty f_k p_k(x) = (p_0(x) , p_1(x) , \\ldots ) \\begin{pmatrix}f_0\\\\ f_1\\\\f_2\\\\\\vdots\\end{pmatrix}\n",
"$$\n",
"then (provided the relevant sums converge absolutely and uniformly)\n",
"$$\n",
"x f(x) = x (p_0(x) , p_1(x) , \\ldots ) \\begin{pmatrix}f_0\\\\ f_1\\\\f_2\\\\\\vdots\\end{pmatrix} =\n",
" \\left(J \\begin{pmatrix} p_0(x) \\cr p_1(x) \\cr p_2(x) \\cr \\vdots \\end{pmatrix}\\right)^\\top \\begin{pmatrix}f_0\\\\ f_1\\\\f_2\\\\\\vdots\\end{pmatrix} = (p_0(x) , p_1(x) , \\ldots ) J^\\top \\begin{pmatrix}f_0\\\\ f_1\\\\f_2\\\\\\vdots\\end{pmatrix}\n",
"$$\n",
"\n",
"\n",
"*Example* For the case of Chebyshev polynomials, we have \n",
"$$\n",
"J = \\begin{pmatrix} \n",
"0 & 1 \\cr\n",
"\\half & 0 & \\half \\cr\n",
"& \\half & 0 & \\half \\cr\n",
"&& \\half & 0 & \\ddots \\cr\n",
"&&&\\ddots & \\ddots\n",
"\\end{pmatrix} \n",
"$$\n",
"Therefore, the Chebyshev coefficients of $x f(x)$ are given by\n",
"$$\n",
"J^\\top \\vc f = \\begin{pmatrix} \n",
"0 & \\half \\cr\n",
"1 & 0 & \\half \\cr\n",
"& \\half & 0 & \\half \\cr\n",
"&& \\half & 0 & \\ddots \\cr\n",
"&&&\\ddots & \\ddots\n",
"\\end{pmatrix} \\begin{pmatrix} f_0\\\\ f_1\\\\f_2\\\\f_3\\\\\\vdots\\end{pmatrix}\n",
"$$\n",
"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$):"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"n = 14\n"
]
},
{
"data": {
"text/plain": [
"15×14 Adjoint{Float64,Array{Float64,2}}:\n",
" 0.0 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 1.0 0.0 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.5 0.0 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.5 0.0 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.5 0.0 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.5 0.0 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.5 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.5 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.5 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.5 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.5 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.5 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.5\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f = Fun(exp, Chebyshev())\n",
"n = ncoefficients(f) # number of coefficients\n",
"@show n\n",
"J = zeros(n,n+1)\n",
"J[1,2] = 1\n",
"for k=2:n\n",
" J[k,k-1] = J[k,k+1] = 1/2\n",
"end\n",
"J'"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4.163336342344337e-17"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cfs = J'*f.coefficients # coefficients of x*f\n",
"xf = Fun(Chebyshev(), cfs)\n",
"\n",
"xf(0.1) - 0.1*f(0.1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Evaluating polynomials\n",
"\n",
"\n",
"We can use the three-term recurrence to construct the polynomials. I think it's nicest to express this in terms of linear algebra. suppose we are given $p_0(x) = 1$ (which is pretty much always the case in practice). This can be written in matrix form as\n",
"$$\n",
"(1,0,0,0,0,\\ldots) \\begin{pmatrix} p_0(x) \\cr p_1(x) \\cr p_2(x) \\cr \\vdots \\end{pmatrix} = 1\n",
"$$\n",
"We can combine this with the Jacobi operator to get\n",
"$$\n",
"\\underbrace{\\begin{pmatrix}\n",
"1 \\\\\n",
"a_0-x & b_0 \\\\\n",
"c_1 & a_1-x & b_1 \\\\\n",
"& c_2 & a_2-x & b_2 \\cr\n",
"&& c_3 & a_3-x & b_3 \\cr\n",
"&&&\\ddots & \\ddots & \\ddots\n",
"\\end{pmatrix}}_{L_x} \\begin{pmatrix} p_0(x) \\cr p_1(x) \\cr p_2(x) \\cr \\vdots \\end{pmatrix} = \\begin{pmatrix} 1\\cr 0 \\cr 0 \\cr \\vdots \\end{pmatrix}\n",
"$$\n",
"For $x$ fixed, this is a lower triangular system, that is, the polynomials equal \n",
"$$\n",
"L_x^{-1} \\vc e_0\n",
"$$\n",
"This can be solved via forward recurrence:\n",
"\\begin{align*}\n",
" p_0(x) &= 1\\\\\n",
" p_1(x) &= {(x-a_0) p_0(x) \\over b_0}\\\\\n",
" p_2(x) &= {(x-a_1) p_0(x) - c_1 p_0(x) \\over b_1}\\\\ \n",
" p_3(x) &= {(x-a_2) p_1(x) - c_2 p_1(x) \\over b_2}\\\\ \n",
" &\\vdots\n",
"\\end{align*}\n",
"\n",
"\n",
"**Example** We can construct $T_0(x),\\ldots,T_{n-1}(x)$ via\n",
"\\begin{align*}\n",
" p_0(x) &= 1\\\\\n",
" p_1(x) &= x T_0(x) \\\\\n",
" T_2(x) &= 2x T_0(x) - T_0(x) \\\\ \n",
" T_3(x) &= 2x T_1(x) - T_1(x)\n",
" &\\vdots\n",
"\\end{align*}\n",
"Believe it or not, this is much faster than using $\\cos k \\acos x$:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"trig_Chebyshev (generic function with 1 method)"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function recurrence_Chebyshev(n,x)\n",
" T = zeros(n)\n",
" T[1] = 1.0\n",
" T[2] = x*T[1]\n",
" for k = 2:n-1\n",
" T[k+1] = 2x*T[k] - T[k-1]\n",
" end\n",
" T\n",
"end\n",
"\n",
"trig_Chebyshev(n,x) = [cos(k*acos(x)) for k=0:n-1]"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.1102230246251565e-16"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n = 10\n",
"recurrence_Chebyshev(n, 0.1) - trig_Chebyshev(n,0.1) |>norm"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.000031 seconds (6 allocations: 78.359 KiB)\n",
" 0.000269 seconds (6 allocations: 78.359 KiB)\n"
]
}
],
"source": [
"n = 10000\n",
"@time recurrence_Chebyshev(n, 0.1) \n",
"@time trig_Chebyshev(n,0.1);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can use this to evaluate functions as well: \n",
"$$\n",
"f(x) = (p_0(x),p_1(x),\\ldots) \\begin{pmatrix}f_0 \\\\ f_1\\\\ \\vdots \\end{pmatrix} = \\vc e_0^\\top L_x^{-\\top} \\begin{pmatrix}f_0 \\\\ f_1\\\\ \\vdots \\end{pmatrix}\n",
"$$\n",
"when $f$ is a degree $n-1$ polynomial, this becomes a problem of inverting an upper triangular matrix, that is, we want to solve the $n \\times n$ system\n",
"$$\n",
"\\underbrace{\\begin{pmatrix}\n",
"1 & a_0-x & c_1 \\\\\n",
"& b_0 & a_1-x & c_2 \\\\\n",
"& & b_1 & a_2-x & \\ddots \\\\\n",
"& & & b_2 & \\ddots & c_{n-2} \\\\\n",
"&&&&\\ddots & a_{n-2}-x \\\\\n",
"&&&&& b_{n-2}\n",
"\\end{pmatrix}}_{L_x^\\top} \\begin{pmatrix} \\gamma_0 \\\\\\vdots\\\\ \\gamma_{n-1} \\end{pmatrix}\n",
"$$\n",
"so that $f(x) = \\gamma_0$. We we can solve this via back-substitution:\n",
"\\begin{align*}\n",
"\\gamma_{n-1} &= {f_{n-1} \\over b_{n-2}} \\\\\n",
"\\gamma_{n-2} &= {f_{n-2} - (a_{n-2}-x) \\gamma_{n-1} \\over b_{n-3}} \\\\\n",
"\\gamma_{n-3} &= {f_{n-3} - (a_{n-3}-x) \\gamma_{n-2} - c_{n-2} \\gamma_{n-1} \\over b_{n-4}} \\\\\n",
"& \\vdots \\\\\n",
"\\gamma_1 &= {f_1 - (a_1-x) \\gamma_2 - c_2 \\gamma_3 \\over b_0} \\\\\n",
"\\gamma_0 &= f_0 - (a_0-x) \\gamma_1 - c_1 \\gamma_2\n",
"\\end{align*}\n",
"\n",
"*Example* For Chebyshev, we want to solve the system\n",
"$$\n",
"\\underbrace{\\begin{pmatrix}\n",
"1 & -x & \\half \\\\\n",
"& 1 & -x & \\half \\\\\n",
"& & \\half & -x & \\ddots \\\\\n",
"& & & \\half & \\ddots & \\half \\\\\n",
"&&&&\\ddots & -x \\\\\n",
"&&&&& \\half\n",
"\\end{pmatrix}}_{L_x^\\top} \\begin{pmatrix} \\gamma_0 \\\\\\vdots\\\\ \\gamma_{n-1} \\end{pmatrix}\n",
"$$\n",
"via\n",
"\n",
"\\begin{align*}\n",
"\\gamma_{n-1} &= 2f_{n-1} \\\\\n",
"\\gamma_{n-2} &= 2f_{n-2} + 2x \\gamma_{n-1} \\\\\n",
"\\gamma_{n-3} &= 2 f_{n-3} + 2x \\gamma_{n-2} - \\gamma_{n-1} \\\\\n",
"& \\vdots \\\\\n",
"\\gamma_1 &= f_1 + x \\gamma_2 - \\half \\gamma_3 \\\\\n",
"\\gamma_0 &= f_0 + x \\gamma_1 - \\half \\gamma_2\n",
"\\end{align*}\n",
"\n",
"then $f(x) = \\gamma_0$."
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"clenshaw_Chebyshev (generic function with 1 method)"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function clenshaw_Chebyshev(f,x)\n",
" n = length(f)\n",
" γ = zeros(n)\n",
" γ[n] = 2f[n]\n",
" γ[n-1] = 2f[n-1] +2x*f[n]\n",
" for k = n-2:-1:1\n",
" γ[k] = 2f[k] + 2x*γ[k+1] - γ[k+2]\n",
" end\n",
" γ[2] = f[2] + x*γ[3] - γ[4]/2\n",
" γ[1] = f[1] + x*γ[2] - γ[3]/2 \n",
" γ[1]\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-1.3322676295501878e-15"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f = Fun(exp, Chebyshev())\n",
"clenshaw_Chebyshev(f.coefficients, 0.1) - exp(0.1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With some high performance computing tweeks, this can be made more accurate: this is the algorithm used for evaluating functions in ApproxFun:"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f(0.1) - exp(0.1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Gram–Schmidt, revisited.\n",
"\n",
"Remember last lecture we introduced the Gram–Schmidt approach to constructing OPs. But the three-term recurrence means we can simplify it, and calculate the recurrence coefficients at the same time:\n",
"\n",
"\n",
"**Proposition (Gram–Schmidt)** Define\n",
"\\begin{align*}\n",
"p_0(x) &= 1 \\\\\n",
"q_0(x) &= {1 \\over \\norm{p_0}}\\\\\n",
"a_n &= \\ip \\\\\n",
"b_{n-1} &= \\ip\\\\\n",
"p_{n+1}(x) &= x q_n(x) - a_n q_n(x) - b_{n-1} q_{n-1}(x)\\\\\n",
"q_{n+1}(x) &= {p_{n+1}(x) \\over \\norm{p_n}}\n",
"\\end{align*}\n",
"Then $q_0(x), q_1(x), \\ldots$ are orthonormal w.r.t. $w$.\n",
"\n",
"**Remark** This can be made a bit more efficient by using $\\norm{p_n}$ to calculate $b_n$."
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [],
"source": [
"x = Fun()\n",
"w = exp(x)\n",
"ip = (f,g) -> sum(f*g*w)\n",
"nrm = f -> sqrt(ip(f,f))\n",
"n = 10\n",
"q = Array{Fun}(undef, n)\n",
"p = Array{Fun}(undef, n)\n",
"a = zeros(n)\n",
"b = zeros(n)\n",
"p[1] = Fun(1, -1 .. 1 )\n",
"q[1] = p[1]/nrm(p[1])\n",
"\n",
"p[2] = x*q[1]\n",
"a[1] = ip(p[2],q[1])\n",
"p[2] -= a[1]*q[1]\n",
"q[2] = p[2]/nrm(p[2])\n",
"\n",
"for k=2:n-1\n",
" p[k+1] = x*q[k] \n",
" b[k-1] =ip(p[k+1],q[k-1])\n",
" a[k] = ip(p[k+1],q[k])\n",
" p[k+1] = p[k+1] - a[k]q[k] - b[k-1]q[k-1]\n",
" q[k+1] = p[k+1]/nrm(p[k+1])\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3.2959746043559335e-16"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ip(q[5],q[2])"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p = plot(; legend=false)\n",
"for k=1:10\n",
" plot!(q[k])\n",
"end\n",
"p"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"9.45569169473238e-16"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"norm(x*q[4] - (b[3]q[3] + a[4]q[4] + b[4]q[5]))"
]
}
],
"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
}