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

File Size: | 0 kb |

File Type: | 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.

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.

## Leave a Reply.

## Archives

February 2015

January 2015

November 2014

September 2014

August 2014

July 2014

June 2014

May 2014

April 2014

March 2014

February 2014

January 2014

December 2013

November 2013

October 2013

September 2013

August 2013

July 2013

June 2013

April 2013

March 2013

February 2013

January 2013

December 2012

November 2012

September 2012

August 2012

July 2012

June 2012

May 2012

April 2012

## Categories

All

Basis Functions

Biharmonic

Boundary Conditions

Boundary Value Problem

Boundary Value Problems

Chebyshev

Differentiation

Discontinuous Galerkin

Domain Decomposition

Eigenvalue Problems

Fast Transform

Finite Difference

Finite Differences

Finite Element

Fourth Order

Galerkin

Galerkin Method

Initial Value Problems

Integral Equation

Integration

Iterative Solver

Jacobi Polynomials

Kronecker Product

Krylov Method

Lagrange Interpolation

Lagrange Polynomial

Legendre

Manufactured Solution

Modal

Newton

Nodal Galerkin

Nonuniform

Numerical Differentiation

Numerical Integration

Optimization

Plotting

Poisson Equation

Preconditioning

Pseudospectral

Pseudospectral Method

Quadrature

Quantum

Radau Quadrature

Radial Equation

Sipg

Spectral Methods

Symmetric

Symplectic Integration

Time-dependent

Triangles

Two Dimensions

Weak Form