Lecture 15.ipynb

{
 "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
}