The Lagrange and Hermite interpolation polynomials already appeared in the previous post, but here I'd like to look at a couple of efficient implementations for polynomial interpolation. First, I am going to assume that the function we want to interpolate is sufficiently smooth and that we want to do the interpolation on the interval \([-1,1]\). Any other finite interval can be mapped to this standard interval with an affine transformation. |

chebinterp.py |

Barycentric Lagrange interpolation

\[

\ell(x) =(x-x_0)(x-x_1)\cdots(x-x_n)

\]

and the

*barycentric weights*as

\[

w_j = \left(\prod\limits_{\substack{j=0\\j\neq k}^n}(x_j-x_k) \right)^{-1}

\]

so that we can write the \(j\)th Lagrange polynomial as

\[

\ell_j(x)=\ell(x)\frac{w_j}{x-x_j}

\]

so that the interpolated function is

\[

p(x)=\ell(x)\sum\limits_{j=0}^n \frac{w_j f(x_j)}{x-x_j}

\]

The barycentric formula is

\[p(x) = \frac{\sum\limits_{j=0}^n \frac{w_j f(x_j)}{x-x_j}}{\sum\limits_{j=0}^n \frac{w_j }{x-x_j}

}\]

I was thinking about writing my own Barycentric Lagrange interpolation code, however, there is already one included in scipy: scipy.interpolate.BarycentricInterpolator

Here's a sample usage. Let's say we want to interpolate \(e^x\) on \([-1,1]\). First we need a grid on which to evaluate the function. I choose the five point Chebyshev-Gauss-Lobatto grid

\[ x_j=\cos\left(\frac{\pi j}{4}\right),\quad j=0,1,2,3,4\]

so that with exact arithmetic, we have \(x=\{-1,-\tfrac{1}{\sqrt{2}},0,\tfrac{1}{\sqrt{2}},1\}\)

import numpy as np

x = np.cos(np.pi*np.arange(5)/float(4))

Now we need to make an interpolator object

from scipy.interpolate import BarycentricInterpolator as bi

L = bi(x) # initialize the object given the weights

the function we want to interpolate is

y = np.exp(x)

L.set_yi(y)

Now we can evaluate it on a fine grid of 20 uniformly distributed points

xf = np.linspace(-1,1,20)

yf = np.exp(xf) # Exact function on fine grid

yint = L(xf) # Interpolated function on fine grid

Finally we can see how close to the exact function the interpolated function is

# Let's save space by using only four digits after the decimal

np.set_printoptions(precision=4)

yint-yf

the result is

array([ 0.0000e+00, -4.0254e-04, -2.7953e-04, 8.9807e-05, 4.9947e-04, 8.1042e-04, 9.4373e-04, 8.7294e-04, 6.1552e-04, 2.2343e-04, -2.2735e-04, -6.4848e-04, -9.5226e-04, -1.0659e-03, -9.4782e-04, -6.0488e-04, -1.1263e-04, 3.6303e-04, 5.4141e-04, 0.0000e+00])

We could also get the barycentric weights with L.wi

array([ 1., -2., 2., -2., 1.])

This class also supports adding additional interpolation points with L.add_xi.

## Comments

What about choosing the Chebyshev nodes in order to minimize the error between f(x) and p(x)? From what I see here you are providing these nodes in the nparray x.

I saw this on the book by K. Atkinson...

I think you are referring to the Chebyshev-Gauss nodes which have a very uniform Lebesgue function. We can do the following:

import numpy as np

n = 10 # for example

x = -np.cos((2*np.arange(n+1)-1)*np.pi/(2*n))

I think you will find that provided you use a nodal distribution which is quadratically clustered at the boundaries you get a small Lebesgue constant. The roots and extrema of both Chebyshev and Legendre polynomials are particularly good. I'd recommend taking a look at Walter Gautschi's "Numerical Analysis: An Introduction," if you are interested.

I was talking about this: http://www.math.uiowa.edu/~atkinson/ftp/ENA_Materials/Overheads/sec_4-6.pdf

I don´t know if it corresponds to what you´ve answered me. Could u answer in my email:?

How about if we apply Chebyshev interpolation to the problem where we employed Lagrange polynomials.

y" + 9y = 0,y(0)=0,y(1)=1

I know you said that the Chebyshev polynomials don't hold at the BC's but we can check and see how the code behaves.

Again I am facing the same problem of addressing 'f' since in my case it is 0.