{
"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",
"<br>\n",
"s.olver@imperial.ac.uk\n",
"\n",
"<br>\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<x f, g> = \\int_a^b x f(x) g(x) w(x) \\dx = \\ip<f, x g>\n",
"$$\n",
"Therefore, if $f_m$ is a degree $m < n-1$ polynomial, we have\n",
"$$\n",
"\\ip<x p_n, f_m> = \\ip<p_n, x f_m> = 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<x q_n, q_{n-1}> = \\ip<q_n, x q_{n-1}> = \\ip<x q_{n-1}, q_n> = 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<x p_n,p_{n+1}> \\over \\norm{p_{n+1}}^2} = {\\ip<x p_{n+1}, p_n> \\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<x q_n, q_n> \\\\\n",
"b_{n-1} &= \\ip<x q_n, q_{n-1}>\\\\\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": [
"<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
"<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 2400 1600\">\n",
"<defs>\n",
" <clipPath id=\"clip2200\">\n",
" <rect x=\"0\" y=\"0\" width=\"2000\" height=\"2000\"/>\n",
" </clipPath>\n",
"</defs>\n",
"<defs>\n",
" <clipPath id=\"clip2201\">\n",
" <rect x=\"0\" y=\"0\" width=\"2400\" height=\"1600\"/>\n",
" </clipPath>\n",
"</defs>\n",
"<polygon clip-path=\"url(#clip2201)\" points=\"\n",
"0,1600 2400,1600 2400,0 0,0 \n",
" \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
"<defs>\n",
" <clipPath id=\"clip2202\">\n",
" <rect x=\"480\" y=\"0\" width=\"1681\" height=\"1600\"/>\n",
" </clipPath>\n",
"</defs>\n",
"<polygon clip-path=\"url(#clip2201)\" points=\"\n",
"149.361,1503.47 2321.26,1503.47 2321.26,47.2441 149.361,47.2441 \n",
" \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
"<defs>\n",
" <clipPath id=\"clip2203\">\n",
" <rect x=\"149\" y=\"47\" width=\"2173\" height=\"1457\"/>\n",
" </clipPath>\n",
"</defs>\n",
"<polyline clip-path=\"url(#clip2203)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 210.632,1503.47 210.632,47.2441 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2203)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 722.971,1503.47 722.971,47.2441 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2203)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 1235.31,1503.47 1235.31,47.2441 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2203)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 1747.65,1503.47 1747.65,47.2441 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2203)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 2259.99,1503.47 2259.99,47.2441 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2203)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 149.361,1331.6 2321.26,1331.6 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2203)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 149.361,1043.57 2321.26,1043.57 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2203)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 149.361,755.532 2321.26,755.532 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2203)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 149.361,467.498 2321.26,467.498 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2203)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 149.361,179.464 2321.26,179.464 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2201)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 149.361,1503.47 2321.26,1503.47 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2201)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 149.361,1503.47 149.361,47.2441 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2201)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 210.632,1503.47 210.632,1481.63 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2201)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 722.971,1503.47 722.971,1481.63 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2201)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 1235.31,1503.47 1235.31,1481.63 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2201)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 1747.65,1503.47 1747.65,1481.63 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2201)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 2259.99,1503.47 2259.99,1481.63 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2201)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 149.361,1331.6 181.939,1331.6 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2201)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 149.361,1043.57 181.939,1043.57 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2201)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 149.361,755.532 181.939,755.532 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2201)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 149.361,467.498 181.939,467.498 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2201)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 149.361,179.464 181.939,179.464 \n",
" \"/>\n",
"<g clip-path=\"url(#clip2201)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 210.632, 1557.47)\" x=\"210.632\" y=\"1557.47\">-1.0</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip2201)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 722.971, 1557.47)\" x=\"722.971\" y=\"1557.47\">-0.5</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip2201)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1235.31, 1557.47)\" x=\"1235.31\" y=\"1557.47\">0.0</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip2201)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1747.65, 1557.47)\" x=\"1747.65\" y=\"1557.47\">0.5</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip2201)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 2259.99, 1557.47)\" x=\"2259.99\" y=\"1557.47\">1.0</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip2201)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 125.361, 1349.1)\" x=\"125.361\" y=\"1349.1\">-4</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip2201)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 125.361, 1061.07)\" x=\"125.361\" y=\"1061.07\">-2</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip2201)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 125.361, 773.032)\" x=\"125.361\" y=\"773.032\">0</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip2201)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 125.361, 484.998)\" x=\"125.361\" y=\"484.998\">2</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip2201)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 125.361, 196.964)\" x=\"125.361\" y=\"196.964\">4</text>\n",
"</g>\n",
"<polyline clip-path=\"url(#clip2203)\" style=\"stroke:#009af9; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 2259.54,661.594 2255.94,661.594 2248.76,661.594 2238.02,661.594 2223.75,661.594 2206.01,661.594 2184.87,661.594 2160.39,661.594 2132.66,661.594 2101.77,661.594 \n",
" 2067.85,661.594 2031,661.594 1991.35,661.594 1949.05,661.594 1904.24,661.594 1857.09,661.594 1807.74,661.594 1756.39,661.594 1703.21,661.594 1648.38,661.594 \n",
" 1592.1,661.594 1534.57,661.594 1475.99,661.594 1416.57,661.594 1356.5,661.594 1296.01,661.594 1235.31,661.594 1174.61,661.594 1114.12,661.594 1054.05,661.594 \n",
" 994.628,661.594 936.046,661.594 878.516,661.594 822.238,661.594 767.412,661.594 714.229,661.594 662.877,661.594 613.535,661.594 566.378,661.594 521.57,661.594 \n",
" 479.268,661.594 439.623,661.594 402.773,661.594 368.846,661.594 337.964,661.594 310.233,661.594 285.752,661.594 264.605,661.594 246.869,661.594 232.605,661.594 \n",
" 221.862,661.594 214.68,661.594 211.082,661.594 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2203)\" style=\"stroke:#e26f46; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 2259.59,632.753 2256.36,633.316 2249.93,634.439 2240.3,636.119 2227.51,638.351 2211.6,641.128 2192.61,644.441 2170.62,648.28 2145.68,652.632 2117.88,657.484 \n",
" 2087.3,662.821 2054.04,668.626 2018.2,674.88 1979.9,681.564 1939.26,688.657 1896.4,696.136 1851.47,703.979 1804.59,712.16 1755.92,720.653 1705.62,729.432 \n",
" 1653.83,738.47 1600.73,747.737 1546.48,757.205 1491.25,766.844 1435.22,776.624 1378.55,786.513 1321.43,796.481 1264.05,806.496 1206.57,816.527 1149.19,826.542 \n",
" 1092.07,836.51 1035.41,846.399 979.37,856.179 924.139,865.817 869.888,875.286 816.786,884.553 765.001,893.59 714.696,902.37 666.029,910.863 619.154,919.044 \n",
" 574.217,926.887 531.359,934.366 490.717,941.459 452.418,948.143 416.582,954.397 383.321,960.202 352.742,965.539 324.939,970.391 300.001,974.743 278.006,978.582 \n",
" 259.022,981.895 243.111,984.672 230.321,986.904 220.693,988.584 214.258,989.707 211.035,990.269 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2203)\" style=\"stroke:#3da44d; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 2259.63,604.59 2256.72,606.267 2250.92,609.599 2242.25,614.543 2230.71,621.039 2216.36,629.002 2199.23,638.332 2179.36,648.91 2156.82,660.601 2131.66,673.256 \n",
" 2103.97,686.711 2073.81,700.794 2041.27,715.321 2006.45,730.103 1969.45,744.946 1930.36,759.654 1889.31,774.03 1846.39,787.88 1801.75,801.015 1755.5,813.25 \n",
" 1707.78,824.413 1658.72,834.341 1608.46,842.883 1557.14,849.904 1504.91,855.286 1451.91,858.928 1398.3,860.748 1344.23,860.684 1289.85,858.696 1235.31,854.766 \n",
" 1180.77,848.895 1126.39,841.108 1072.32,831.454 1018.71,819.999 965.715,806.832 913.483,792.062 862.163,775.817 811.901,758.242 762.839,739.496 715.116,719.756 \n",
" 668.868,699.209 624.226,678.051 581.315,656.489 540.259,634.734 501.173,613.001 464.167,591.507 429.348,570.466 396.813,550.092 366.655,530.589 338.959,512.156 \n",
" 313.805,494.98 291.262,479.238 271.395,465.089 254.261,452.679 239.907,442.136 228.375,433.568 219.697,427.064 213.899,422.69 210.995,420.491 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2203)\" style=\"stroke:#c271d2; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 2259.66,581.283 2257.03,584.588 2251.78,591.118 2243.92,600.713 2233.47,613.136 2220.45,628.083 2204.91,645.188 2186.88,664.029 2166.41,684.144 2143.54,705.039 \n",
" 2118.35,726.199 2090.89,747.104 2061.23,767.238 2029.45,786.106 1995.64,803.242 1959.87,818.225 1922.24,830.685 1882.85,840.316 1841.79,846.883 1799.18,850.227 \n",
" 1755.13,850.273 1709.73,847.026 1663.12,840.581 1615.41,831.111 1566.73,818.874 1517.2,804.199 1466.94,787.482 1416.08,769.181 1364.77,749.8 1313.12,729.879 \n",
" 1261.27,709.984 1209.35,690.691 1157.5,672.574 1105.85,656.192 1054.54,642.074 1003.68,630.706 953.425,622.521 903.89,617.888 855.206,617.102 807.498,620.375 \n",
" 760.888,627.831 715.495,639.504 671.437,655.333 628.827,675.163 587.773,698.749 548.381,725.757 510.753,755.775 474.985,788.315 441.169,822.829 409.391,858.716 \n",
" 379.733,895.336 352.271,932.024 327.076,968.103 304.213,1002.9 283.74,1035.76 265.709,1066.05 250.167,1093.2 237.155,1116.68 226.704,1136.03 218.843,1150.89 \n",
" 213.59,1160.96 210.961,1166.04 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2203)\" style=\"stroke:#ac8d18; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 2259.69,560.903 2257.3,566.29 2252.52,576.861 2245.36,592.216 2235.85,611.776 2224,634.805 2209.84,660.436 2193.4,687.708 2174.73,715.603 2153.86,743.084 \n",
" 2130.85,769.138 2105.74,792.811 2078.6,813.253 2049.5,829.742 2018.49,841.722 1985.65,848.818 1951.06,850.853 1914.8,847.857 1876.95,840.065 1837.6,827.906 \n",
" 1796.85,811.99 1754.78,793.082 1711.5,772.072 1667.11,749.941 1621.71,727.721 1575.41,706.453 1528.31,687.146 1480.53,670.735 1432.18,658.042 1383.37,649.743 \n",
" 1334.21,646.335 1284.82,648.117 1235.31,655.171 1185.8,667.357 1136.41,684.313 1087.26,705.463 1038.44,730.04 990.089,757.106 942.308,785.588 895.211,814.316 \n",
" 848.909,842.061 803.51,867.587 759.118,889.693 715.839,907.258 673.774,919.288 633.019,924.953 593.672,923.621 555.823,914.889 519.561,898.595 484.97,874.84 \n",
" 452.133,843.98 421.124,806.624 392.017,763.618 364.879,716.019 339.774,665.06 316.761,612.118 295.893,558.663 277.219,506.211 260.783,456.279 246.623,410.326 \n",
" 234.772,369.713 225.258,335.651 218.103,309.161 213.324,291.043 210.931,281.845 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2203)\" style=\"stroke:#00a9ad; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 2259.72,542.566 2257.53,550.425 2253.16,565.735 2246.62,587.704 2237.92,615.196 2227.09,646.797 2214.13,680.889 2199.09,715.738 2181.99,749.59 2162.87,780.768 \n",
" 2141.77,807.758 2118.74,829.295 2093.82,844.429 2067.07,852.576 2038.54,853.551 2008.3,847.57 1976.41,835.241 1942.94,817.528 1907.96,795.689 1871.54,771.21 \n",
" 1833.76,745.718 1794.71,720.885 1754.47,698.332 1713.11,679.537 1670.74,665.749 1627.44,657.911 1583.3,656.605 1538.42,662.018 1492.89,673.921 1446.81,691.685 \n",
" 1400.28,714.308 1353.4,740.469 1306.26,768.604 1258.98,796.987 1211.64,823.834 1164.36,847.402 1117.22,866.093 1070.34,878.548 1023.81,883.736 977.733,881.017 \n",
" 932.204,870.197 887.322,851.547 843.183,825.805 799.88,794.154 757.507,758.167 716.153,719.738 675.908,680.991 636.856,644.173 599.081,611.541 562.664,585.244 \n",
" 527.682,567.21 494.211,559.038 462.321,561.906 432.081,576.501 403.554,602.97 376.803,640.896 351.884,689.307 328.85,746.702 307.751,811.119 288.631,880.212 \n",
" 271.531,951.357 256.488,1021.77 243.534,1088.63 232.697,1149.22 223.999,1201.04 217.459,1241.94 213.091,1270.2 210.905,1284.63 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2203)\" style=\"stroke:#ed5d92; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 2259.74,525.777 2257.73,536.439 2253.73,557.054 2247.73,586.259 2239.74,622.128 2229.79,662.309 2217.9,704.195 2204.08,745.099 2188.37,782.45 2170.78,813.968 \n",
" 2151.37,837.825 2130.17,852.768 2107.21,858.199 2082.55,854.217 2056.22,841.597 2028.3,821.727 1998.81,796.502 1967.84,768.173 1935.43,739.183 1901.65,711.977 \n",
" 1866.56,688.821 1830.24,671.632 1792.76,661.833 1754.18,660.243 1714.59,667.021 1674.06,681.65 1632.67,702.976 1590.5,729.293 1547.64,758.475 1504.16,788.139 \n",
" 1460.16,815.822 1415.72,839.18 1370.93,856.166 1325.87,865.197 1280.64,865.29 1235.31,856.143 1189.99,838.182 1144.75,812.541 1099.69,780.998 1054.9,745.851 \n",
" 1010.46,709.766 966.456,675.58 922.982,646.095 880.12,623.865 837.953,610.997 796.563,608.975 756.033,618.531 716.441,639.557 677.864,671.083 640.379,711.307 \n",
" 604.058,757.691 568.973,807.107 535.193,856.031 502.783,900.771 471.807,937.706 442.325,963.537 414.396,975.519 388.074,971.658 363.41,950.874 340.453,913.096 \n",
" 319.248,859.304 299.836,791.495 282.255,712.592 266.54,626.281 252.721,536.81 240.826,448.734 230.878,366.648 222.895,294.913 216.895,237.389 212.888,197.198 \n",
" 210.883,176.537 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2203)\" style=\"stroke:#c68125; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 2259.76,510.213 2257.91,523.955 2254.22,550.326 2248.7,587.195 2241.35,631.603 2232.18,680.016 2221.22,728.634 2208.48,773.716 2193.99,811.903 2177.77,840.506 \n",
" 2159.86,857.721 2140.27,862.782 2119.06,856.004 2096.25,838.741 2071.9,813.243 2046.03,782.436 2018.7,749.65 1989.97,718.3 1959.87,691.574 1928.46,672.134 \n",
" 1895.81,661.879 1861.97,661.768 1826.99,671.744 1790.96,690.74 1753.92,716.794 1715.94,747.234 1677.1,778.945 1637.46,808.665 1597.1,833.304 1556.09,850.248 \n",
" 1514.5,857.617 1472.4,854.459 1429.88,840.853 1387.01,817.922 1343.86,787.739 1300.52,753.147 1257.06,717.5 1213.56,684.349 1170.1,657.108 1126.76,638.72 \n",
" 1083.61,631.359 1040.74,636.206 998.219,653.298 956.125,681.49 914.533,718.525 873.519,761.202 833.158,805.646 793.521,847.64 754.68,883.006 716.705,907.992 \n",
" 679.665,919.634 643.626,916.062 608.654,896.717 574.81,862.467 542.157,815.589 510.753,759.644 480.655,699.22 451.917,639.594 424.59,586.311 398.725,544.732 \n",
" 374.367,519.582 351.56,514.538 330.346,531.881 310.763,572.267 292.846,634.61 276.627,716.109 262.136,812.409 249.399,917.889 238.438,1026.06 229.274,1130.04 \n",
" 221.923,1223.03 216.398,1298.88 212.709,1352.48 210.863,1380.23 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2203)\" style=\"stroke:#00a98d; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 2259.78,495.652 2258.07,512.705 2254.66,545.179 2249.56,589.983 2242.77,642.879 2234.3,698.927 2224.17,752.989 2212.39,800.252 2198.98,836.716 2183.98,859.588 \n",
" 2167.39,867.549 2149.25,860.855 2129.59,841.282 2108.45,811.895 2085.84,776.697 2061.83,740.169 2036.44,706.775 2009.71,680.469 1981.7,664.268 1952.44,659.93 \n",
" 1921.99,667.782 1890.4,686.698 1857.71,714.255 1823.99,747.019 1789.29,780.954 1753.67,811.894 1717.19,836.028 1679.9,850.346 1641.88,852.997 1603.17,843.505 \n",
" 1563.86,822.845 1524,793.335 1483.65,758.389 1442.9,722.127 1401.8,688.908 1360.42,662.826 1318.83,647.229 1277.11,644.309 1235.31,654.825 1193.52,677.977 \n",
" 1151.79,711.458 1110.2,751.678 1068.82,794.136 1027.72,833.912 986.967,866.215 946.625,886.932 906.763,893.12 867.448,883.384 828.745,858.091 790.719,819.397 \n",
" 753.433,771.082 716.949,718.19 681.327,666.528 646.628,622.056 612.908,590.238 580.225,575.424 548.632,580.319 518.181,605.616 488.925,649.81 460.91,709.255 \n",
" 434.184,778.424 408.792,850.387 384.776,917.445 362.175,971.867 341.027,1006.66 321.368,1016.27 303.23,997.229 286.643,948.533 271.636,871.859 258.232,771.502 \n",
" 246.454,654.068 236.323,527.941 227.854,402.573 221.062,287.654 215.958,192.258 212.55,124.021 210.845,88.4582 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip2203)\" style=\"stroke:#8e971d; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 2259.79,481.929 2258.21,502.483 2255.05,541.328 2250.33,594.206 2244.03,655.38 2236.18,718.309 2226.79,776.43 2215.87,823.925 2203.43,856.39 2189.51,871.329 \n",
" 2174.11,868.406 2157.26,849.428 2139,818.062 2119.34,779.338 2098.31,738.978 2075.96,702.658 2052.31,675.29 2027.4,660.412 2001.27,659.759 1973.95,673.087 \n",
" 1945.5,698.238 1915.96,731.468 1885.36,767.972 1853.76,802.542 1821.21,830.28 1787.75,847.255 1753.45,851.039 1718.34,841.031 1682.49,818.547 1645.95,786.644 \n",
" 1608.78,749.714 1571.03,712.897 1532.76,681.386 1494.03,659.72 1454.91,651.154 1415.44,657.2 1375.7,677.385 1335.75,709.278 1295.63,748.779 1255.43,790.634 \n",
" 1215.19,829.119 1174.99,858.796 1134.87,875.26 1094.92,875.754 1055.18,859.591 1015.71,828.315 976.588,785.562 937.862,736.655 899.595,687.952 861.845,646.048 \n",
" 824.672,616.912 788.131,605.077 752.28,612.988 717.174,640.596 682.867,685.258 649.412,741.962 616.86,803.866 585.261,863.095 554.665,911.696 525.119,942.644 \n",
" 496.667,950.785 469.355,933.576 443.223,891.545 418.313,828.398 394.663,750.74 372.309,667.443 351.285,588.713 331.625,524.969 313.358,485.642 296.512,478.049 \n",
" 281.115,506.467 267.188,571.515 254.754,669.933 243.832,794.784 234.439,936.065 226.59,1081.68 220.295,1218.64 215.566,1334.43 212.409,1418.26 210.829,1462.26 \n",
" \n",
" \"/>\n",
"</svg>\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
}