Interpolation

Newton interpolation :

from numpy import *
import scipy as sp
import matplotlib.pyplot as plt
from numpy.linalg import inv
from numpy import matmul

def nj(x_local,x_list,kmax):
	out = 1
	for i in range(kmax):
		out = out*(x_local-x_list[i])
	return out

x = sp.array([-2,-1,0,1,2,3])
y = x*x*x

npoints = len(x)


# create matrix to solve
ns = sp.zeros((npoints,npoints))
for i in range(npoints):
	for j in range(npoints): 
		ns[i][j] = nj(x[i],x,j)

# solve matrix
inv_ns = inv(ns)
coefs  = matmul(inv_ns,y)

newx = sp.linspace(-2.,3.,64)
newy = sp.zeros(64)

# get new points
for i in range(len(newx)):
	n = 0
	for j in range(npoints):
		n = n+coefs[j]*nj(newx[i],x,j)
	newy[i] = n

plt.plot(x,y)
plt.plot(newx,newy)
plt.show()

Lagrange interpolation:

from numpy import *
import scipy as sp
import matplotlib.pyplot as plt
from numpy.linalg import inv
from numpy import matmul

def lj(xi,x,j):
	out = 1
	for m in range(len(x)):
		if(m!=j):
			out = out*(xi-x[m])/(x[j]-x[m])
	return out

x = sp.array([-2,-1,0,1,2,3])
y = x*x*x

newx = sp.linspace(-2.,3.,64)
newy = sp.zeros(64)

for i in range(len(newx)):
	n = 0
	for j in range(len(x)):
		n = n+y[j]*lj(newx[i],x,j)
	newy[i] = n

plt.plot(x,y)
plt.plot(newx,newy)
plt.show()