Polynomial interpolation

Contents

Polynomial interpolant

Given some points $(x_0,y_0),(x_1,y_1),\ldots,(x_n,y_n)$ we want to find a polynomial function that passes through all of them.

Consider the points:

$$C=\{(2,2),(3,6),(4,5),(5,5),(6,6)\}$$

Let us draw the points.

To insert the graphs in the notebook.

In [2]:
%matplotlib inline

We import the modules matplotlib.pyplot y numpy.

In [3]:
import numpy as np
import matplotlib.pyplot as plt

We store the points in two numpy arrays

In [4]:
x = np.array([2.,3.,4.,5.,6.])
y = np.array([2.,6.,5.,5.,6.])

Now we can draw the points

In [5]:
plt.plot(x,y,'ro')
plt.axis([1.5,6.5,1.5,6.5]);   # We draw in the rectangle [1.5,6.5]x[1.5,6.5]

Solution with Vandermonde's matrix

As we have $5$ points we can write $5$ equations and we need $5$ unknows that will be the coefficients of the polynomial. Thus:

$$P_4(x)=a_0+a_1x+a_2x^2+a_3x^3+a_4x^4$$

that is $a_0,a_1,a_2,a_3,a_4$, that are $5$ unknowns and the polynomial degree is $4$. In general, the polynomial that passes throught the points $(x_0,y_0),(x_1,y_1),\ldots,(x_n,y_n)$ has a degree of $\leq n$.

The equations are:

$$ \begin{array}{lll} P_{4}(x_{0})=y_{0} & \quad & a_{0}+a_{1}x_{0}+a_{2}x_{0}^{2}+a_{3}x_{0}^{3}+a_{4}x_{0}^{4}=y_{0}\\ P_{4}(x_{1})=y_{1} & \quad & a_{0}+a_{1}x_{1}+a_{2}x_{1}^{2}+a_{3}x_{1}^{3}+a_{4}x_{1}^{4}=y_{1}\\ P_{4}(x_{2})=y_{2} & \quad & a_{0}+a_{1}x_{2}+a_{2}x_{2}^{2}+a_{3}x_{2}^{3}+a_{4}x_{2}^{4}=y_{2}\\ P_{4}(x_{3})=y_{3} & \quad & a_{0}+a_{1}x_{3}+a_{2}x_{3}^{2}+a_{3}x_{3}^{3}+a_{4}x_{3}^{4}=y_{3}\\ P_{4}(x_{4})=y_{4} & \quad & a_{0}+a_{1}x_{4}+a_{2}x_{4}^{2}+a_{3}x_{4}^{3}+a_{4}x_{4}^{4}=y_{4} \end{array} $$

And the system to solve is:

$$ \left(\begin{array}{ccccc} 1 & x_{0} & x_{0}^{2} & x_{0}^{3} & x_{0}^{4}\\ 1 & x_{1} & x_{1}^{2} & x_{1}^{3} & x_{1}^{4}\\ 1 & x_{2} & x_{2}^{2} & x_{2}^{3} & x_{2}^{4}\\ 1 & x_{3} & x_{3}^{2} & x_{3}^{3} & x_{3}^{4}\\ 1 & x_{4} & x_{4}^{2} & x_{4}^{3} & x_{4}^{4} \end{array}\right)\left(\begin{array}{c} a_{0}\\ a_{1}\\ a_{2}\\ a_{3}\\ a_{4} \end{array}\right)=\left(\begin{array}{c} y_{0}\\ y_{1}\\ y_{2}\\ y_{3}\\ y_{4} \end{array}\right) $$

This system coefficient matrix is called Vandermonde's matrix.

Exercise 1

  • Write a function Vandermonde(x) whose input is the x vector that contains the $x$ coordinates of the points (nodes), and as output, the Vandermonde's matrix V.
  • Compute the interpolation polynomial coefficients solving the system Va=y. Use the command numpy.linalg.solve.
  • Plot the nodes and the polynomial. Use the command numpy.polyval to compute the polynomial in some points from the coefficients computed in the previous question.

Use the nodes:

In [6]:
x = np.array([2.,3.,4.,5.,6.])
y = np.array([2.,6.,5.,5.,6.])

Check that the script also works for the nodes:

In [7]:
x = np.array([0.,1.,2.,3.,4.,5.,6.])
y = np.array([3.,5.,6.,5.,4.,4.,5.])
In [8]:
%run Exercise1.py

We can also show the Vandermonde's matrix and the polynomial coefficients

In [9]:
# print options
np.set_printoptions(precision = 2)   # only 2 fractionary digits
np.set_printoptions(suppress = True) # do not use exponential notation

print 'V matrix'
print V

V matrix
[[    1.     2.     4.     8.    16.]
 [    1.     3.     9.    27.    81.]
 [    1.     4.    16.    64.   256.]
 [    1.     5.    25.   125.   625.]
 [    1.     6.    36.   216.  1296.]]
In [10]:
print 'Polynomial coefficients (from a0 to a4)'
print a

Polynomial coefficients (from a0 to a4)
[-75.    81.   -29.25   4.5   -0.25]

There is a Python function that computes Vandermonde's matrix. The element order is different.

In [11]:
va = np.vander(x, len(x))
print(va)

[[   16.     8.     4.     2.     1.]
 [   81.    27.     9.     3.     1.]
 [  256.    64.    16.     4.     1.]
 [  625.   125.    25.     5.     1.]
 [ 1296.   216.    36.     6.     1.]]

Nevertheless, it is not advisable to use Vandermonde's matriz to contruct the polynomial interpolant because, the matrix, in general, is badly conditioned and small errors in the data can lead to large errors in the solution of the system.

Interpolation with Lagrange fundamental polynomials

For each $k=0,1,\ldots,n$, there exists a unique polynomial $\ell_{k}$ of degree $\leq n$ such as $\ell_{k}\left(x_{j}\right)=\delta_{kj}$ (one if $k=j$, zero otherwise):

$$ \ell_{k}\left(z\right)= \frac{(z-x_{0})\ldots(z-x_{k-1})(z-x_{k+1})\ldots(z-x_{n})}{(x_{k}-x_{0}) \ldots(x_{k}-x_{k-1})(x_{k}-x_{k+1})\ldots(x_{k}-x_{n})}. $$

The polynomials $\ell_{0},\ell_{1},\ldots,\ell_{n}$ are the Lagrange fundamental polynomials.

The Lagrange polynomial for the points $x_{0,}x_{1},\ldots,x_{n}$ with values $y_{0},y_{1},\ldots,y_{n}$ is

$$ P_{n}\left(x\right)=y_{0}\ell_{0}\left(x\right)+ y_{1}\ell_{1}\left(x\right)+\cdots+y_{n}\ell_{n}\left(x\right) $$

Exercise 2

Write a script that computes the Lagrange's interpolant polynomial using the Lagrange fundamental polynomials. Create two functions: lagrange_fundamental(k,z,x) with

  • Input: k, index of the fundamental polynomial, z, point (or array of points) where we will evaluate the polynomial and x, the $x$ coordinates of the nodes (or points the polynomial passes through).
  • Output: the value of the $k-$th fundamental polynomial in the point (or array of points) z.

and a function lagrange_polynomial(z,x,y) with:

  • Input: z, point (or array of points) where we will evaluate the polynomial and the coordinates of the nodes x and y.
  • Output: the value of the Lagrange's polynomial in the point (or array or points) z.

Plot the Lagrange fundamental polynomial with index k=2 in the interval $[2,6]$. Draw the interpolant polynomial in $[2,6]$.

Use the nodes:

In [12]:
x = np.array([2.,3.,4.,5.,6.])
y = np.array([2.,6.,5.,5.,6.])

Check that it also works for the nodes:

In [13]:
x = np.array([0.,1.,2.,3.,4.,5.,6.])
y = np.array([3.,5.,6.,5.,4.,4.,5.])
In [14]:
%run Exercise2.py

A drawback of this method is that if we want to add a new node we have to do all the computations again.

Interpolation with divided differences

Formula of Newton

The Lagrange interpolant polynomial of degree n for the nodes $(x_{0},y_{0}),(x_{1},,y_{1})\ldots,(x_{n},y_{n})$ is given by

$$ P_{n}\left(x\right)= \left[y_{0}\right]+\left[y_{0},y_{1}\right]\left(x-x_{0}\right)+ \left[y_{0},y_{1},y_{2}\right]\left(x-x_{0}\right)\left(x-x_{1}\right)+\cdots +\left[y_{0},y_{1},\ldots,y_{n}\right]\left(x-x_{0}\right)\left(x-x_{1}\right)\cdots\left(x-x_{n-1}\right) $$

Table of divided differences

The coefficients of the latter polynomial are in the first row of the following table:

$$ \begin{array}{cccccc} x_{0} & y_{0} & \left[y_{0},y_{1}\right]=\frac{y_{0}-y_{1}}{x_{0}-x_{1}} & \left[y_{0},y_{1},y_{2}\right]=\frac{\left[y_{0},y_{1}\right]-\left[y_{1},y_{2}\right]}{x_{0}-x_{2}} & \ldots & \left[y_{0},y_{1},\ldots,y_{n}\right]=\frac{\left[y_{0},\ldots,y_{n-1}\right]-\left[y_{1},\ldots,y_{n}\right]}{x_{0}-x_{n}}\\ x_{1} & y_{1} & \left[y_{1},y_{2}\right]=\frac{y_{1}-y_{2}}{x_{1}-x_{2}} & \left[y_{1},y_{2},y_{3}\right]=\frac{\left[y_{1},y_{2}\right]-\left[y_{2},y_{3}\right]}{x_{1}-x_{3}} & \ldots\\ x_{2} & y_{2} & \left[y_{2},y_{3}\right]=\frac{y_{2}-y_{3}}{x_{2}-x_{3}} & \left[y_{2},y_{3},y_{4}\right]=\frac{\left[y_{2},y_{3}\right]-\left[y_{3},y_{4}\right]}{x_{2}-x_{4}} & \ldots\\ \ldots & \ldots & \ldots & \ldots & \ldots\\ x_{n-1} & y_{n-1} & \left[y_{n-1},y_{n}\right]=\frac{y_{n-1}-y_{n}}{x_{n-1}-x_{n}}\\ x_{n} & y_{n} \end{array} $$

Exercise 3

Write a script that computes the Lagrange's interpolant polynomial using the divided differences. Create two functions: div_dif(x,y) with

  • Input: the coordinate coordinates of the nodes x and y.
  • Output: a matrix that contains the divide differences.

and a function newton_polynomial(z,x,y) with:

  • Input: z, point (or array of points) where we will evaluate the polynomial and the coordinates of the nodes x and y.
  • Output: the value of the Lagrange's polynomial in the point (or array or points) z.

Print the divided differences matrix. Plot the Lagrange fundamental polynomial with index k=2 in the interval $[2,6]$. Draw the interpolant polynomial in $[2,6]$.

Use the nodes:

In [15]:
x = np.array([2.,3.,4.,5.,6.])
y = np.array([2.,6.,5.,5.,6.])

Comprobar que también funciona con los nodos:

In [16]:
x = np.array([0.,1.,2.,3.,4.,5.,6.])
y = np.array([3.,5.,6.,5.,4.,4.,5.])
In [17]:
%run Exercise3.py

[[ 2.    2.    4.   -2.5   1.   -0.25]
 [ 3.    6.   -1.    0.5   0.    0.  ]
 [ 4.    5.    0.    0.5   0.    0.  ]
 [ 5.    5.    1.    0.    0.    0.  ]
 [ 6.    6.    0.    0.    0.    0.  ]]

Interpolation with python functions

The interpolant polynomial can be computed with numpy function polyfit if we choose as polynomial degree the number of the nodes minus one. The output is the coefficients of the polynomial, from $a_n$ to $a_0$.

In [18]:
x = np.array([2.,3.,4.,5.,6.])
y = np.array([2.,6.,5.,5.,6.])
pol = np.polyfit(x,y,len(x)-1)  # polynomial coefficients

xx = np.linspace(min(x),max(x))
yy = np.polyval(pol,xx)         # polynomial value in the points contained in xx

We plot the polynomial and the nodes

In [19]:
plt.plot(xx, yy, '-',x, y, 'ro')
plt.axis([min(xx)-d, max(xx)+d, min(yy)-d, max(yy)+d]);

Exercise 4

Write a function chebyshev(f,a,b,n) that interpolates function f in the interval [a,b] using n nodes:

  1. Using equispaced nodes that include the interval boundaries.
  2. Using Chebyshev nodes, that are defined in the interval $[-1,1]$ as:
$$ x_{i}^{(n)}=\cos\frac{\left(2i-1\right)\pi}{2n}\quad i=1,2,\ldots,n $$

the number of nodes is n = 11, the interval is [a,b]=[-1,1] and the function f is

$$ f(x)=\frac{1}{1+25x^{2}} $$

To interpolate use the functions numpy.polyfit and numpy.polyval.

In both cases, plot the nodes, the function and the interpolant polynomial.

Check that it also works with 15 nodes and the function:

$$ f(x)=e^{-20x^2} $$
In [20]:
%run Exercise4.py

We can compute piecewise interpolation with scipy functions

In [21]:
from scipy.interpolate import interp1d

f1 = interp1d(x, y, kind='nearest')
f2 = interp1d(x, y, kind='linear')
f3 = interp1d(x, y, kind='cubic')

xx = np.linspace(min(x),max(x))

Zero degree

In [22]:
plt.plot(x, y, 'ro', xx, f1(xx), '-')
plt.axis([min(xx)-d, max(xx)+d, min(yy)-d, max(yy)+d]);

Degree one or linear

In [23]:
plt.plot(x, y, 'ro', xx, f2(xx), '-')
plt.axis([min(xx)-d, max(xx)+d, min(yy)-d, max(yy)+d]);

Degree three or cubic

In [24]:
plt.plot(x, y, 'ro', xx, f3(xx), '-')
plt.axis([min(xx)-d, max(xx)+d, min(yy)-d, max(yy)+d]);

And interpolation with splines (with option not-a-knot)

In [25]:
from scipy import interpolate
In [26]:
s = interpolate.InterpolatedUnivariateSpline(x, y)
xx = np.linspace(min(x),max(x))
yy = s(xx)

plt.figure()
plt.plot(x, y, 'ro', xx, yy)
plt.axis([min(xx)-d, max(xx)+d, min(yy)-d, max(yy)+d]);

Exercise 5

Write a script that, for the function $f(x) = \cos x$ in interval $[0,10]$ with nodes $x=0,1,\ldots,10$ plot the polynomial interpolants

  • Linear piecewise.
  • Cubic splines.

and calculates, for the points

xx = np.linspace(0,10,100)

the absolute error for each case. Use the funtion numpy.linalg.norm.

$$E_a = \left\Vert P(xx)-f(xx) \right\Vert $$
In [28]:
%run Exercise5.py

Piecewise linear interpolation error = 0.64506
Spline interpolation error           = 0.06776