Lagrange Polynomial Interpolation on Python.

It's a whole a lot easier than Newton's divided differences interpolation polynomial, because there is no divided difference part that need a recursive function.

I use these data points
(0,0)
(1,1)
(2,4)
(4,16)
(5,25)

from pylab import *

def L(i,xs,x,y):
    Ls = 1
    for j in arange (len(x)):
        if i != j:
            Ls *= (xs-x[j])/(x[i]-x[j])

    return Ls

def f(xs,x,y):
    fs =0 
    for i in arange(len(x)):
        fs += L(i,xs,x,y)*y[i]

    return fs

x = []
y = []
#initial value
x.append(0)
x.append(1)
x.append(2)
x.append(4)
x.append(5)

y.append(0)
y.append(1)
y.append(4)
y.append(16)
y.append(25)

#for testing
xc = 3

yc = f(xc,x,y)

print ''
print xc, yc
#plot
t = linspace(-6,6,100)
u = f(t,x,y)
plot(t,u)

grid(True)
show()

.

with other data points

(0,0)
(1,1)
(2,4)
(4,-16)
(5,25)