Polynomial interpolation

We import the modules matplotlib.pyplot y numpy and the numpy functions for polynomials.

In [1]:
import numpy as np
import matplotlib.pyplot as plt
import numpy.polynomial.polynomial as pol

Lagrange polynomial

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=\{(-1,1),(0,3),(2,4),(3,3),(5,1)\}$$

Let us plot them.

We store the points in two numpy arrays

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

And we plot the points and the Lagrange polynomial

In [3]:
xp = np.linspace(min(x),max(x))
p  = pol.polyfit(x,y,len(x)-1)
yp = pol.polyval(xp,p)

plt.figure()
plt.plot(xp,yp,'b-', label = 'interpolant polynomial')
plt.plot( x, y,'ro', label = 'points')
plt.legend()
plt.show()

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})}{(x_{k}-x_{0})}\cdots \frac{(z-x_{k-1})}{(x_{k}-x_{k-1})} \cdot \frac{(z-x_{k+1})}{(x_{k}-x_{k+1})} \cdots \frac{(z-x_{n})}{(x_{k}-x_{n})} $$

For instance, with our 5 nodes, these polynomials are

$$ \ell_{0}\left(z\right)= \Bigg\rvert^{x_0} \frac{(z-x_1)}{(x_0-x_1)} \cdot \frac{(z-x_2)}{(x_0-x_2)} \cdot \frac{(z-x_3)}{(x_0-x_3)} \cdot \frac{(z-x_4)}{(x_0-x_4)} $$$$ \ell_{1}\left(z\right)= \frac{(z-x_0)}{(x_1-x_0)} \Bigg\rvert^{x_1} \frac{(z-x_2)}{(x_1-x_2)} \cdot \frac{(z-x_3)}{(x_1-x_3)} \cdot \frac{(z-x_4)}{(x_1-x_4)} $$$$ \ell_{2}\left(z\right)= \frac{(z-x_0)}{(x_2-x_0)} \cdot \frac{(z-x_1)}{(x_2-x_1)} \Bigg\rvert^{x_2} \frac{(z-x_3)}{(x_2-x_3)} \cdot \frac{(z-x_4)}{(x_2-x_4)} $$$$ \ell_{3}\left(z\right)= \frac{(z-x_0)}{(x_3-x_0)} \cdot \frac{(z-x_1)}{(x_3-x_1)} \cdot \frac{(z-x_2)}{(x_3-x_2)} \Bigg\rvert^{x_3} \frac{(z-x_4)}{(x_3-x_4)} $$$$ \ell_{4}\left(z\right)= \frac{(z-x_0)}{(x_4-x_0)} \cdot \frac{(z-x_1)}{(x_4-x_1)} \cdot \frac{(z-x_2)}{(x_4-x_2)} \cdot \frac{(z-x_3)}{(x_4-x_3)} \Bigg\rvert^{x_4} $$

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(z\right)=y_{0}\ell_{0}\left(z\right)+ y_{1}\ell_{1}\left(z\right)+\cdots+y_{n}\ell_{n}\left(z\right) $$

Exercise 1

Write a script that computes Lagrange's fundamental polynomials. Create a function: lagrange_fundamental(k,x,z) with

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

Plot the Lagrange's fundamental polynomials.

Use the nodes:

x = np.array([-1., 0, 2, 3, 5])

Hints:

  • You can start building a function that works for a float $z$ and then modify it so you can use a vector as parameter. For example, the first version would work
In [6]:
k = 2
z = 1.3
yz = lagrange_fundamental(k,x,z)
print(yz)

1.0448388888888889

so we know that $\ell_{2}\left(1.3\right)= 1.04484$. The second version would do

In [7]:
k = 2
z = np.array([1.3, 2.1, 3.2])
yz = lagrange_fundamental(k,x,z)
print(yz)

[ 1.04483889  0.94395    -0.2688    ]

so we know that $\ell_{2}\left([1.3,\, 2.1,\, 3.2]\right)= [ 1.04483889,\, 0.94395,\, -0.2688 ]$

  • To plot the dots, you can take into account that the $x$s are in the x vector of the nodes and the corresponding $y$s can be obtained as the rows of the identity matrix:
In [8]:
print(np.eye(len(x)))

[[1. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0.]
 [0. 0. 1. 0. 0.]
 [0. 0. 0. 1. 0.]
 [0. 0. 0. 0. 1.]]
In [9]:
%run Exercise1.py

Let us remember that 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(z\right)=y_{0}\ell_{0}\left(z\right)+ y_{1}\ell_{1}\left(z\right)+\cdots+y_{n}\ell_{n}\left(z\right) $$

Exercise 2

Write a script that computes the Lagrange's interpolant polynomial using the Lagrange fundamental polynomials. Use the function lagrange_fundamental(k,x,z) in the function lagrange_polynomial(x,y,z):

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

Plot the Lagrange's interpolant polynomial.

Use the nodes:

x = np.array([-1., 0, 2, 3, 5])
y = np.array([ 1., 3, 4, 3, 1])

Check that the script also works for the nodes:

x1 = np.array([-1., 0, 2, 3, 5, 6, 7])
y1 = np.array([ 1., 3, 4, 3, 2, 2, 1])
In [10]:
%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.

Chebyshev nodes

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 [11]:
x = np.array([-1., 0, 2, 3, 5])
y = np.array([ 1., 3, 4, 3, 1])
p = pol.polyfit(x,y,len(x)-1)  # polynomial coefficients

xp = np.linspace(min(x),max(x))
yp = pol.polyval(xp,p)         # polynomial value in the points contained in xp

We plot the polynomial and the nodes

In [12]:
plt.figure()
plt.plot(xp, yp, 'b-',x, y, 'ro')
plt.show()

Exercise 3

Write a function chebyshev(f,a,b,n) that interpolates the 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_1(x)=\frac{1}{1+25x^{2}} $$

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

$$ f_2(x)=e^{-20x^2} $$

Follow the steps:

  • Build the nodes: the n nodes will be stored in an n element vector x in both cases. For the equispaced nodes the functionse np.linspace or np.arange may be used. The corresponding y are calculated using f1 and f2.

  • Build the interpolant polynomials: To interpolate the functions, use pol.polyfit and pol.polyval.

  • Plot the graphs: In both cases, plot the nodes, the function and the interpolant polynomial. Use plt.axis([-1.05, 1.05, -0.3, 2.3]) so all the graphs are plotted in the rectangle $[-1.05, 1.05]\times[ -0.3, 2.3]$. Use more than 50 points for np.linespace as the interpolant polynomials will oscillate a lot.

In [13]:
%run Exercise3.py

------------  Function f1  ------------


------------  Function f2  ------------

More exercises

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(z\right)= \left[y_{0}\right]+\left[y_{0},y_{1}\right]\left(z-x_{0}\right)+ \left[y_{0},y_{1},y_{2}\right]\left(z-x_{0}\right)\left(z-x_{1}\right)+\cdots +\left[y_{0},y_{1},\ldots,y_{n}\right]\left(z-x_{0}\right)\left(z-x_{1}\right)\cdots\left(z-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 4

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 difference matrix. Draw the interpolant polynomial.

Use the nodes:

x = np.array([-1., 0, 2, 3, 5])
y = np.array([ 1., 3, 4, 3, 1])

Check that the script also works for the nodes:

x1 = np.array([-1., 0, 2, 3, 5, 6, 7])
y1 = np.array([ 1., 3, 4, 3, 2, 2, 1])
In [14]:
%run Exercise4.py

DIVIDED DIFFERENCES TABLE
[[-1.    1.    2.   -0.5   0.    0.02]
 [ 0.    3.    0.5  -0.5   0.1   0.  ]
 [ 2.    4.   -1.    0.    0.    0.  ]
 [ 3.    3.   -1.    0.    0.    0.  ]
 [ 5.    1.    0.    0.    0.    0.  ]]

We can compute piecewise interpolation with scipy functions

In [15]:
x = np.array([-1., 0, 2, 3, 5])
y = np.array([ 1., 3, 4, 3, 1])
In [16]:
from scipy.interpolate import interp1d, CubicSpline

f1 = interp1d(x, y, kind='nearest')
f2 = interp1d(x, y, kind='linear')
f3 = CubicSpline(x, y, bc_type='natural')

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

Zero degree

In [17]:
plt.figure()
plt.plot(x, y, 'ro', xp, f1(xp), 'b-')
plt.show()

Degree one or linear

In [18]:
plt.figure()
plt.plot(x, y, 'ro', xp, f2(xp), 'b-')
plt.show()

Degree three, a natural cubic spline

In [19]:
plt.plot(x, y, 'ro', xp, f3(xp), 'b-')
plt.show()

Exercise 5

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

  • Linear piecewise.
  • Cubic splines.

and calculates, for the points

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

the absolute error for each case. Use the funtion Ea = np.linalg.norm(yp-f(xp)).

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

Piecewise linear interpolation error = 
Spline interpolation error           =

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) $$

that is, $Vp=y$, where $V$ is the system coefficient matrix, called Vandermonde's matrix, $y$, that contains the independent terms of the system and the system solution, the polynomial coefficients, contained in $p$.

Exercise 6

  1. Write a function, Vandermonde(x), whose input argument is the x vector that contains the $x$ coordinates of the points (nodes), and as output argument, the Vandermonde's matrix V.
  2. Write a function polVandermonde(x,y), whose input argumets are the coordinates x and y of the nodes and whose output argument is p, that contains the coefficients of the interpolation polynomial, that is, $p = (a_0,a_1,a_2,a_3,a_4)$..
    • Compute V calling the Vandermonde function.
    • Compute the interpolation polynomial coefficients solving the system $Vp = y$. Use the command p = np.linalg.solve(V,y).
  3. Plot the nodes and the polynomial.
    • xp = np.linspace(min(x),max(x)) creates a vector with 50 points in the interpolation interval.
    • Use the command yp = pol.polyval(xp,p) to compute the polynomial in the points xp using the coefficients p computed in the previous question.

Use the nodes:

x = np.array([-1., 0, 2, 3, 5])
y = np.array([ 1., 3, 4, 3, 1])

Check that the script also works for the nodes:

x1 = np.array([-1., 0, 2, 3, 5, 6, 7])
y1 = np.array([ 1., 3, 4, 3, 2, 2, 1])
In [21]:
%run Exercise6.py

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

In [22]:
print('V matrix')
print(V)

V matrix
[[  1.  -1.   1.  -1.   1.]
 [  1.   0.   0.   0.   0.]
 [  1.   2.   4.   8.  16.]
 [  1.   3.   9.  27.  81.]
 [  1.   5.  25. 125. 625.]]
In [23]:
print('Polynomial coefficients')
print(p)

Polynomial coefficients
[ 3.    1.6  -0.48 -0.07  0.02]

There is a Python function that computes Vandermonde's matrix

In [24]:
va = pol.polyvander(x, len(x)-1)
print(va)

[[  1.  -1.   1.  -1.   1.]
 [  1.   0.   0.   0.   0.]
 [  1.   2.   4.   8.  16.]
 [  1.   3.   9.  27.  81.]
 [  1.   5.  25. 125. 625.]]

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