{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# Polynomial interpolation \n", "\n", " \n", "**Anne Kværnø**\n", "\n", "Date: **Oct 16, 2018**\n", "\n", "# Introduction\n", "In the first part of Mathematics 4D/N you learned about the Fourier-series.\n", "By taking the sum of the series from 1 to $N$, you can find an *approximation* to\n", "a periodic function. \n", "Similar, polynomials can be used to approximate functions over some bounded\n", "interval $x \\in [a,b]$. Such polynomials can be used for different purposes.\n", "The function itself may be unknown, and only measured data are available. In\n", "this case, a polynomial may be used to find approximations to intermediate\n", "values of the function. Polynomials are\n", "easy to integrate, and can be used to find approximations of integrals of more complicated\n", "functions. This will be exploited later in the course. And there are plenty of other applications.\n", "\n", "In this part of the course, we will only discuss *interpolation polynomials*.\n", "\n", "**Interpolation problem.**\n", "\n", "Given $n+1$ points $(x_i,y_i)_{i=0}^n$. Find a polynomial $p(x)$ of lowest possible degree\n", "satisfying the *interpolation condition*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", " \\label{eq:intcond} \\tag{1}\n", " p(x_i) = y_i,\\qquad i=0,\\dotsc, n. \n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The solution $p(x)$ is called the *interpolation polynomial*, the $x_i$ values\n", "are called *nodes*, and the points $(x_i,y_i)$ *interpolation points*.\n", "\n", "\n", "\n", "\n", "**Example 1:** \n", "Given the points" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{array}{c|c|c|c}\n", "x_i & 0 & 2/3 & 1 \\\\ \\hline\n", "y_i & 1 & 1/2 & 0 \n", "\\end{array}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The corresponding interpolation polynomial is" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "p_2(x)=(-3x^2-x+4)/4\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The $y$-values of this example are chosen such that $y_i=\\cos{(\\pi x_i/2)}$. So\n", "$p_2(x)$ can be considered as an approximation to $\\cos{(\\pi x/2)}$ on the interval\n", "$[0,1]$. \n", "\n", "\n", "### Content of this note\n", "\n", "In this part, we will discuss the following: \n", "* Method: How to compute the polynomials?\n", "\n", "* Existence and uniqueness results. \n", "\n", "* Error analysis: If the polynomial is used to approximate a function, how good is the approximation?\n", "\n", "* Improvements: If the nodes $x_i$ can be chosen freely, how should we do it in order to reduce the error? \n", "\n", "## Notation etc.\n", "\n", "Let us start with some useful notation and facts about polynomials. \n", "* A polynomial of degree $n$ is given by" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \\label{eq:polynomial} \\tag{2}\n", " p_n(x) = c_{n}x^n + c_{n-1}x^{n-1} + \\cdots + c_1 x_1 + c_0, \\qquad c_i \\in\n", "\\mathbb{R}, \\quad i=0,1,\\dotsc,n.\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* $\\mathbb{P}_n$ is the set of all polynomials of degree $n$. \n", "\n", "* $C^m[a,b]$ is the set of all continuous functions that have continuous first $m$ derivatives.\n", "\n", "* The value $r$ is a root or a zero of a polynomial $p$ if $p(r)=0$.\n", "\n", "* A nonzero polynomial of degree $n$ can never have more than $n$ real roots (there may be less). \n", "\n", "* A polynomial of degree $n$ with $n$ real roots $r_1,r_2,\\dotsc,r_n$ can be written as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "p_n(x) = c(x-r_1)(x-r_2)\\dotsm(x-r_n) = c\\prod_{i=1}^n(x-r_i).\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Methods\n", "In this section, we present three techniques for finding the interpolation polynomial for a given set of data. \n", "\n", "## The direct approach\n", "For a polynomial of degree $n$ the interpolation condition ([eq:intcond](#eq:intcond)) is a linear systems of \n", "$n+1$ equations in $n+1$ unknowns:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\sum_{i=0}^n x_j^i c_i = y_j, \\qquad j=0,\\dotsc, n.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we are basically interested in the polynomials themself, given by the coefficients $c_i$, $i=0,1,\\dotsc, n$, this is a perfectly fine solution. It is for instance the strategy implemented in MATLAB's interpolation routines. However, in this course, polynomial interpolation will be used as a basic tool to construct other algorithms, in particular for integration. In that case, this is not the most convenient option, so we concentrate on a different strategy, which essentially makes it possible to just write up the polynomials. \n", "\n", "\n", "\n", "## Lagrange interpolation\n", "Given $n+1$ points $(x_i,y_i)_{i=0}^n$ with distinct $x_i$ values. \n", "The *cardinal functions* are defined by:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\ell_i(x) = \\prod_{j=0,j\\not=i}^n \\frac{x-x_j}{x_i-x_j}\n", " = \\frac{x-x_0}{x_i-x_0} \\dotsm \\frac{x-x_{i-1}}{x_i-x_{i-1}}\\cdot \\frac{x-x_{i+1}}{x_i-x_{i+1}} \\dotsm \\frac{x-x_n}{x_i-x_n} , \\qquad i=0,\\dotsc,n.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The cardinal functions have the following properties:\n", "* $\\ell_i \\in \\mathbb{P}_n$, $i=0,1,\\dotsc,n$.\n", "\n", "* $\\ell_i(x_j) = \\delta_{ij} = \\begin{cases} 1, & \\text{when } i=j \\\\ 0, & \\text{when }i\\not=j \\end{cases}$.\n", "\n", "* They are constructed solely from the nodes $x_i$'s.\n", "\n", "* They are linearly independent, and thus form a basis for $\\mathbb{P}_{n}$.\n", "\n", "The interpolation polynomial is now given by" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "p_n(x) = \\sum_{i=0}^n y_i \\ell_i(x)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "since" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "p_n(x_j) = \\sum_{i=0}^n y_i \\ell_i(x_j) = y_j, \\qquad j=0,\\dotsc,n.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Example 2:**\n", "Given the points:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{array}{c|ccc}\n", "x_i & 0 & 1 & 3 \\\\ \\hline y_i & 3 & 8 & 6\n", "\\end{array}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The corresponding cardinal functions are given by:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\\begin{align*}\n", " \\ell_0(x) & = \\frac{(x-1)(x-3)}{(0-1)(0-3)}\n", " = \\frac{1}{3}x^2-\\frac{4}{3}x+1 \\\\ \n", " \\ell_1(x) & = \\frac{(x-0)(x-3)}{(1-0)(1-3)}\n", " = -\\frac12 x^2 + \\frac32 x \\\\ \n", " \\ell_2(x) &= \\frac{(x-0)(x-1)}{(3-0)(3-1)} = \\frac16 x^2-\\frac16 x\n", "\\end{align*}\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and the interpolation polynomial is given by (check it yourself):" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "p_2(x) = 3 \\ell_0(x) + 8 \\ell_1(x) + 6 \\ell_2(x) = -2x^2 + 7x + 3.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Implementation\n", "The method above is implemented as two functions:\n", "* cardinal(xdata, x): Create a list of cardinal functions $\\ell_i(x)$ evaluated in $x$.\n", "\n", "* lagrange(ydata, l): Create the interpolation polynomial $p_n(x)$.\n", "\n", "Here, xdata and ydata are arrays with the interpolation points, and x is an \n", "array of values in which the polynomials are evaluated. \n", "\n", "You are not required to understand the implementation of these functions, but you should understand how to use them." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "from numpy import *\n", "from matplotlib.pyplot import *\n", "from math import factorial\n", "newparams = {'figure.figsize': (8.0, 4.0), 'axes.grid': True,\n", " 'lines.markersize': 8, 'lines.linewidth': 2,\n", " 'font.size': 14}\n", "rcParams.update(newparams)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "def cardinal(xdata, x):\n", " \"\"\"\n", " cardinal(xdata, x): \n", " In: xdata, array with the nodes x_i.\n", " x, array or a scalar of values in which the cardinal functions are evaluated.\n", " Return: l: a list of arrays of the cardinal functions evaluated in x. \n", " \"\"\"\n", " n = len(xdata) # Number of evaluation points x\n", " l = []\n", " for i in range(n): # Loop over the cardinal functions\n", " li = ones(len(x))\n", " for j in range(n): # Loop to make the product for l_i\n", " if i is not j:\n", " li = li*(x-xdata[j])/(xdata[i]-xdata[j])\n", " l.append(li) # Append the array to the list \n", " return l\n", "\n", "def lagrange(ydata, l):\n", " \"\"\"\n", " lagrange(ydata, l):\n", " In: ydata, array of the y-values of the interpolation points.\n", " l, a list of the cardinal functions, given by cardinal(xdata, x)\n", " Return: An array with the interpolation polynomial. \n", " \"\"\"\n", " poly = 0 \n", " for i in range(len(ydata)):\n", " poly = poly + ydata[i]*l[i] \n", " return poly" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Example 3:**\n", "Test the functions on the interpolation points of Example 2." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Example 3\n", "xdata = [0, 1, 3] # The interpolation points\n", "ydata = [3, 8, 6]\n", "x = linspace(0, 3, 101) # The x-values in which the polynomial is evaluated\n", "l = cardinal(xdata, x) # Find the cardinal functions evaluated in x\n", "p = lagrange(ydata, l) # Compute the polynomial evaluated in x\n", "plot(x, p) # Plot the polynomial\n", "plot(xdata, ydata, 'o') # Plot the interpolation points \n", "title('The interpolation polynomial p(x)')\n", "xlabel('x');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Numerical exercises:**\n", "1. Plot the cardinal functions for the nodes of Example 1. \n", "\n", "2. Plot the interpolation polynomials for some points of your own choice." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# Insert your code here (use \"+\" in the Toolbar menu for more cells)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Newton interpolation\n", "This is an alternative approach to find the interpolation polynomial.\n", "\n", "Let $x_0,x_1,\\ldots,x_n$ be $n+1$ distinct real numbers. The so-called Newton form of a polynomial of degree $n$ is an expansion of the form" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "p(x)=\\sum_{i=0}^{n-1} c_{n-i}\\prod_{j=0}^{n-1-i}(x-x_j) + c_0,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "or more explicitly" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "p(x)=c_n (x-x_0)(x-x_1)\\cdots(x-x_{n-1}) + c_{n-1}(x-x_0)(x-x_1)\\cdots(x-x_{n-2}) + \\cdots + c_1(x-x_0) + c_0.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the light of this form of writing a polynomial, the polynomial interpolation problem leads to the following observations. Let us start with a single node $x_0$, then $f(x_0)=p(x_0)=c_0$. Going one step further and consider two nodes $x_0,x_1$. Then we see that $f(x_0)=p(x_0)=c_0$ and $f(x_1)=p(x_1)=c_0 + c_1(x_1-x_0)$. The latter implies that the coefficient" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "c_1=\\frac{f(x_1)-f(x_0)}{x_1-x_0}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Given three nodes $x_0,x_1,x_2$ yields the coefficients $c_0,c_1$ as defined above, and from" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "f(x_2)=p(x_2)=c_0 + c_1(x_2-x_0) + c_2(x_2-x_0)(x_2-x_1)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "we deduce the coefficient" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "c_2=\\frac{f(x_2) - f(x_0) - \\frac{f(x_1)-f(x_0)}{x_1-x_0} (x_2-x_0)}{(x_2-x_0)(x_2-x_1)}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Playing with this quotient gives the much more structured expression" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "c_2=\\frac{\\frac{f(x_2)-f(x_1)}{x_2-x_1} - \\frac{f(x_1)-f(x_0)}{x_1-x_0}}{(x_2-x_0)}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This procedure can be continued and yields a so-called triangular systems that permits to define the remaining coefficients $c_3,\\ldots,c_n$. One sees quickly that the coefficient $c_k$ only depends on the interpolation points $(x_0,y_0),\\ldots,(x_k,y_k)$, where $y_i:=f(x_i)$, $i=0,\\ldots,n$.\n", "\n", "We introduce the folllwing so-called finite difference notation for a function $f$. The 0th order finite difference is defined to be $f[x_0]:=f(x_0)$. The 1st order finite difference is" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "f[x_0,x_1]:=\\frac{f(x_1)-f(x_0)}{x_1-x_0}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The second order finite difference is defined by" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "f[x_0,x_1,x_2]:= \\frac{f[x_1,x_2] - f[x_0,x_1]}{x_2-x_0}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In general, the nth order finite difference of the function $f$ is defined to be" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "f[x_0,\\ldots,x_n]:= \\frac{f[x_1,\\ldots,x_n] - f[x_0,\\ldots,x_{n-1}]}{x_n-x_0}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Newton's method to solve the polynomial interpolation problem can be summarized as follows. Given $n+1$ interpolation points $(x_0,y_0),\\ldots,(x_n,y_n)$, $y_i:=f(x_i)$. If the order $n$ interpolation polynomial is expressed in Newton's form" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "p_n(x)=c_n (x-x_0)(x-x_1)\\cdots(x-x_{n-1}) + c_{n-1}(x-x_0)(x-x_1)\\cdots(x-x_{n-2}) + \\cdots + c_1(x-x_0) + c_0,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "then the coefficients" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "c_k = f[x_0,\\ldots,x_k]\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "for $k=0,\\ldots,n$. In fact, a recursion is in place" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "p_n(x)=p_{n-1}(x) + f[x_0,\\ldots,x_n](x-x_0)(x-x_1)\\cdots(x-x_{n-1})\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is common to write the finite differences in a table, which for $n=3$ will\n", "look like:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{array}{c|cccc}\n", "x_0 & f[x_0] & \\\\ \n", " & & f[x_0,x_1] & \\\\ \n", "x_1 & f[x_1] & & f[x_0,x_1,x_2] \\\\ \n", " & & f[x_1,x_2] & & f[x_0,x_1,x_2, x_3] \\\\ \n", "x_2 & f[x_2] & & f[x_1,x_2,x_3] \\\\ \n", " & & f[x_2,x_3] & \\\\ \n", "x_3 & f[x_3] \\\\ \n", "\\end{array}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Example 1 again:**\n", "Given the points in Example 1. The corresponding table of divided differences\n", "becomes:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{array}{c|cccc}\n", "0 & 1 & \\\\ \n", " & & -3/4 & \\\\ \n", "2/3 & 1/2 & & -3/4 \\\\ \n", " & & -3/2 & \\\\ \n", "1 & 0 & \n", "\\end{array}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and the interpolation polynomial becomes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "p_2(x) = 1 - \\frac{3}{4}(x-0)-\\frac{3}{4}(x-0)(x-\\frac23) = 1 - \\frac{1}{4}x -\n", "\\frac{3}{4} x^2.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Implementation\n", "The method above is implemented as two functions:\n", "* divdiff(xdata, ydata): Create the table of divided differences\n", "\n", "* newtonInterpolation(F, xdata, x): Evaluate the interpolation polynomial.\n", "\n", "Here, xdata and ydata are arrays with the interpolation points, and x is an \n", "array of values in which the polynomial is evaluated." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def divdiff(xdata,ydata):\n", " # Create the table of divided differences based\n", " # on the data in the arrays x_data and y_data. \n", " n = len(xdata)\n", " F = zeros((n,n))\n", " F[:,0] = ydata # Array for the divided differences\n", " for j in range(n):\n", " for i in range(n-j-1):\n", " F[i,j+1] = (F[i+1,j]-F[i,j])/(xdata[i+j+1]-xdata[i])\n", " return F # Return all of F for inspection. \n", " # Only the first row is necessary for the\n", " # polynomial.\n", "\n", "def newton_interpolation(F, xdata, x):\n", " # The Newton interpolation polynomial evaluated in x. \n", " n, m = shape(F)\n", " xpoly = ones(len(x)) # (x-x)(x-x)...\n", " newton_poly = F[0,0]*ones(len(x)) # The Newton polynomial\n", " for j in range(n-1):\n", " xpoly = xpoly*(x-xdata[j])\n", " newton_poly = newton_poly + F[0,j+1]*xpoly\n", " return newton_poly" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Run the code on the example above:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# Example: Use of divided differences and the Newton interpolation\n", "# formula. \n", "xdata = [0, 2/3, 1]\n", "ydata = [1, 1/2, 0]\n", "F = divdiff(xdata, ydata) # The table of divided differences\n", "print('The table of divided differences:\\n',F)\n", "\n", "x = linspace(0, 1, 101) # The x-values in which the polynomial is evaluated\n", "p = newton_interpolation(F, xdata, x)\n", "plot(x, p) # Plot the polynomial\n", "plot(xdata, ydata, 'o') # Plot the interpolation points \n", "title('The interpolation polynomial p(x)')\n", "grid(True)\n", "xlabel('x');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Theory\n", "In this section we cover two theoretical aspects, and give the answer to one natural question:\n", "* Do interpolation polynomials always exist, and when they exist, are they unique?\n", "\n", "* If the polynomial is used to approximate a function, can we find an expression for the error?\n", "\n", "* How can the error be made as small as possible? \n", "\n", "## Existence and uniqueness of interpolation polynomials.\n", "We have already proved the existence of such polynomials, simply by constructing\n", "them. But are they unique? \n", "\n", "Suppose there exist two different interpolation polynomials $p_n$ and $q_n$ of\n", "degree $n$ interpolating the same $n+1$ points. The polynomial $r(x) = p_n(x)-q_n(x)$ is of degree $n$\n", "with zeros in all the nodes $x_i$, that is a total of $n+1$ zeros. But then\n", "$r\\equiv 0$, and the two polynomials $p_n$ and $q_n$ are identical. \n", "We have then proved the following result:\n", "\n", "**Theorem: Existence and uniqueness.**\n", "\n", "Given $n+1$ points $(x_i,y_i)_{i=0}^n$ with distinct $x$ values. Then there is\n", "one and only one polynomial $p_n(x) \\in \\mathbb{P}_n$ satisfying the\n", "interpolation condition" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "p_n(x_i) = y_i, \\qquad i=0,\\dotsc, n.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Error Analysis\n", "Given some function $f\\in C[a,b]$. Choose $n+1$ distinct nodes in $[a,b]$ and let $p_n(x) \\in \\mathbb{P}_n$ satisfy the interpolation condition" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "p_n(x_i) = f(x_i), \\qquad i=0,\\dots,n.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What can be said about the error $e(x)=f(x)-p_n(x)$? \n", "\n", "Let us start with an numerical experiment, to have a certain feeling of what to expect. \n", "\n", "**Example 4:**\n", "\n", "Let $f(x)=\\sin(x)$, $x\\in [0,2\\pi]$. Choose $n+1$ equidistributed nodes, that is $x_i=ih$, $i=0,\\dots,n$, and $h=2\\pi/n$. Calculate the interpolation polynomial by use of the functions cardinal and lagrange. Plot the error $e_n(x)=f(x)-p_n(x)$ for different values of $n$. Choose $n=4,8,16$ and $32$. Notice how the error is distributed over the interval, and find the maximum error $\\max_{x\\in[a,b]}|e_n(x)|$ for each $n$." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# Example 4\n", "\n", "# Define the function\n", "def f(x):\n", " return sin(x)\n", "\n", "# Set the interval \n", "a, b = 0, 2*pi # The interpolation interval\n", "x = linspace(a, b, 101) # The 'x-axis' \n", "\n", "# Set the interpolation points\n", "n = 8 # Interpolation points\n", "xdata = linspace(a, b, n+1) # Equidistributed nodes (can be changed)\n", "ydata = f(xdata) \n", "\n", "# Evaluate the interpolation polynomial in the x-values\n", "l = cardinal(xdata, x) \n", "p = lagrange(ydata, l)\n", "\n", "# Plot f(x) og p(x) and the interpolation points\n", "subplot(2,1,1) \n", "plot(x, f(x), x, p, xdata, ydata, 'o')\n", "legend(['f(x)','p(x)'])\n", "grid(True)\n", "\n", "# Plot the interpolation error\n", "subplot(2,1,2)\n", "plot(x, (f(x)-p))\n", "xlabel('x')\n", "ylabel('Error: f(x)-p(x)')\n", "grid(True)\n", "print(\"Max error is {:.2e}\".format(max(abs(p-f(x)))))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Numerical exercise:**\n", "* Repeat the experiment with Runge's function" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "f(x) = \\frac{1}{1+x^2}, \\qquad x\\in [-5,5].\n", "$$" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# Insert your code here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us now see if we can find an expression for the error $e(x)=f(x)-p_n(x)$.\n", "Assume that the nodes $x_i \\in [a,b]$, $i=0,\\dotsc,n$ are distinct, and that $f$\n", "is sufficiently differentiable, that is $f \\in C^r[a,b]$. How differentiable will be clear as we go. \n", "\n", "Define the function" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\omega(x) = \\prod_{i=0}^{n}(x-x_i) = x^{n+1} + \\dotsm.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Clearly, the error in the nodes, $e(x_i)=0$. \n", "Choose an *arbitrary* $x\\in [a,b]$, $x\\in [a,b]$, where $x\\not=x_i$,\n", "$i=0,1,\\dotsc,n$. For this fixed $x$, define a function in $t$ as:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\varphi(t) = e(t)\\omega(x) - e(x)\\omega(t).\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $e(t) = f(t)-p_n(t)$.\n", "Notice that $\\varphi(t)$ is as differentiable with respect to $t$ as $f(t)$. The\n", "function $\\varphi(t)$ has $n+2$ distinct zeros (the nodes and the fixed x). As a\n", "consequence of [Rolle's theorem](https://en.wikipedia.org/wiki/Rolle's_theorem), the derivative\n", "$\\varphi'(t)$ has at least $n+1$ distinct zeros, one between each of the zeros\n", "of $\\varphi(t)$. So $\\varphi''(t)$ has $n$ distinct\n", "zeros, etc. By repeating this argument, we can see that $\\varphi^{n+1}(t)$\n", "has at least one zero in $[a,b]$, let us call this $\\xi(x)$, as it do depend on the fixed $x$. \n", " Since\n", "$\\omega^{(n+1)}(t)=(n+1)!$ and $e^{(n+1)}(t)=f^{(n+1)}(t)$ we can conclude that" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\varphi^{(n+1)}(\\xi)= 0 = f^{(n+1)}(\\xi)\\omega(x) - e(x)(n+1)!,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The argument above is valid if $f$ is $n+1$ continuous\n", "differentiable on $[a,b]$. Altogether we have proved the following:\n", "\n", "**Theorem: Interpolation error.**\n", "\n", "Given $f \\in C^{(n+1)}[a,b]$. Let $p_{n} \\in \\mathbb{P}_n$ interpolate $f$ in\n", "$n+1$ distinct nodes $x_i \\in [a,b]$. For each $x\\in [a,b]$ there is at least\n", "one $\\xi(x) \\in (a,b)$ such that" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "f(x) - p_n(x) = \\frac{f^{(n+1)}(\\xi(x))}{(n+1)!}\\prod_{i=0}^n(x-x_i).\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The interpolation error consists of three elements: The derivative of the\n", "function $f$, the number of interpolation points $n+1$ and the distribution of\n", "the nodes $x_i$. We cannot do much with the first of these, but we can choose\n", "the two others. Let us first look at the most obvious choice of nodes.\n", "\n", "### Equidistributed nodes\n", "\n", "The nodes are *equidistributed* over the interval $[a,b]$ if $x_i=a+ih$, $h=(b-a)/n$. In this case it can\n", "be proved that:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "|\\omega(x)| \\leq \\frac{h^{n+1}}{4}n!\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "such that" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "|e(x)| \\leq \\frac{h^{n+1}}{4(n+1)}M, \\qquad M=\\max_{x\\in[a,b]}|f^{(n+1)}(x)|.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "for all $x\\in [a,b]$. \n", "\n", "Let us now see how good this error bound is by an example.\n", "\n", "**Example 5:**\n", "Let again $f(x)=\\sin(x)$ and $p_n(x)$ the polynomial interpolating $f(x)$ in\n", "$n+1$ equidistributed points. Find an upper bound for the error for different values of $n$. \n", "\n", "Clearly,\n", "$\\max_{x\\in[0,2\\pi]}|f^{(n+1)}(x)|=M=1$ for all $n$, so" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "|e_n(x)| = |f(x)-p_n(x)| \\leq\n", "\\frac{1}{4(n+1)}\\left(\\frac{2\\pi}{n}\\right)^{n+1}, \\qquad x\\in[a,b].\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Use the code in Example 4 to verify the result. How close is the bound to the real error? \n", "\n", "## Optimal choice of interpolation points\n", "So how can the error be reduced? For a given $n$ there is only one choice: to\n", "distribute the nodes in order to make\n", "$|\\omega(x)|= \\prod_{j=0}^{n}|x-x_i|$ as small as possible. We will first do this\n", "on a standard interval $[-1,1]$, and then transfer the results to some arbitrary\n", "interval $[a,b]$.\n", "\n", "Let us start taking a look at $\\omega(x)$ for equidistributed nodes on the\n", "interval $[-1,1]$, for\n", "different values of $n$:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "def omega(xdata, x):\n", " # compute omega(x) for the nodes in xdata\n", " n1 = len(xdata)\n", " omega_value = ones(len(x)) \n", " for j in range(n1):\n", " omega_value = omega_value*(x-xdata[j]) # (x-x_0)(x-x_1)...(x-x_n)\n", " return omega_value" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "# Plot omega(x) \n", "n = 8 # Number of interpolation points is n+1\n", "a, b = -1, 1 # The interval\n", "x = linspace(a, b, 501) \n", "xdata = linspace(a, b, n) \n", "plot(x, omega(xdata, x))\n", "grid(True)\n", "xlabel('x')\n", "ylabel('omega(x)')\n", "print(\"n = {:2d}, max|omega(x)| = {:.2e}\".format(n, max(abs(omega(xdata, x)))))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Run the code for different values of $n$. Notice the following: \n", "* $\\max_{x\\in[-1,1]} |\\omega(x)|$ becomes smaller with increasing $n$. \n", "\n", "* $|\\omega(x)|$ has its maximum values near the boundaries of $[-1, 1]$.\n", "\n", "A a consequence of the latter, it seems reasonable to move the nodes towards the boundaries. \n", "It can be proved that the optimal choice of nodes are the *Chebyshev-nodes*, given by" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\tilde{x}_i = \\cos \\left( \\frac{(2i+1)\\pi}{2(n+1)} \\right), \\qquad i=0,\\dotsc,n\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let $\\omega_{Cheb}(x) = \\prod_{j=1}^n(x-\\tilde{x}_i)$. It is then possible to prove that" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{1}{2^{n}} = \\max_{x\\in [-1, 1]} |\\omega_{Cheb}(x)| \\leq \\max_{x \\in [-1, 1]} |q(x)|\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "for all polynomials $q\\in \\mathbb{P}_n$ such that $q(x)=x^n + c_{n-1}x^{n-1}+\\dotsm+c_1x + c_0$. \n", "\n", "The distribution of nodes can be transferred to an interval $[a,b]$ by the linear transformation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "x = \\frac{b-a}{2}\\tilde{x} + \\frac{b+a}{2}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $x\\in[a,b]$ and $\\tilde{x} \\in [-1,1]$. By doing so we get" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\omega(x) = \\prod_{j=0}^n (x-x_i) =\n", " \\left(\\frac{b-a}{2}\\right)^{n+1} \\prod_{j=0}^n (\\tilde{x}-\\tilde{x}_i)\n", " = \\left(\\frac{b-a}{2}\\right)^{n+1} \\omega_{Cheb}(\\tilde{x}).\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From the theorem on interpolation errors we can conclude:\n", "\n", "**Theorem (interpolation error for Chebyshev interpolation).**\n", "\n", "Given $f \\in C^{(n+1)}[a,b]$, and let $M_{n+1} = \\max_{x\\in [a,b]}|f^{(n+1)}(x)|$. Let $p_{n} \\in \\mathbb{P}_n$ interpolate $f$ i $n+1$ Chebyshev-nodes $x_i \\in [a,b]$. Then" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\max_{x\\in[a,b]}|f(x) - p_n(x)| \\leq \\frac{(b-a)^{n+1}}{2^{2n+1}(n+1)!} M_{n+1}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Chebyshev nodes over an interval $[a,b]$ are evaluated in the following function:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "def chebyshev_nodes(a, b, n):\n", " # n Chebyshev nodes in the interval [a, b] \n", " i = array(range(n)) # i = [0,1,2,3, ....n-1]\n", " x = cos((2*i+1)*pi/(2*(n))) # nodes over the interval [-1,1]\n", " return 0.5*(b-a)*x+0.5*(b+a) # nodes over the interval [a,b]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Numerical exercises:**\n", "1. Plot $\\omega_{Cheb}(x)$ for $3, 5, 9, 17$ interpolation points.\n", "\n", "2. Repeat Example 3 using Chebyshev interpolation on the functions below. Compare with the results you got from equidistributed nodes." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\\begin{align*}\n", " f(x) &= \\sin(x), && x\\in[0,2\\pi] \\\\ \n", " f(x) &= \\frac{1}{1+x^2}, && x\\in[-5,5]. \n", "\\end{align*}\n", "" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "# Insert your code here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**For information**: \n", "[Chebfun](http://www.chebfun.org/) is software package which makes it possible to manipulate functions and to solve equations with accuracy close to machine accuracy. The algorithms are based on polynomial interpolation in Chebyshev nodes." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.3" } }, "nbformat": 4, "nbformat_minor": 2 }