import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

def evalPoly(a,xData,x):
    n = len(xData)-1  # order of polynomial
    p = a[n]
    for i in range(1,n+1):
        #print "a[", (n-i),"]=", a[n-i]
        p = a[n-i] + (x-xData[n-i])*p
    return p

def coefficients(xData,yData):
    n = len(xData)
    a = np.zeros((n, n));
    for i in range(0, n):
        a[i, 0] = yData[i];
    for j in range(1, n):
        for k in range(0, n-j):
            a[k, j] = (a[k+1, j-1] - a[k, j-1])/(xData[k+j]-xData[k])
    return a[0, :] # return the zeroth row


# plot the data points
numDataPoints = 8
xData = np.linspace(-1, 1, numDataPoints)
print xData
yData = 1.0/(1+25*xData**2)
print yData
plt.plot(xData,yData,'bo',markersize=10) 

# plot the exact function
x = np.linspace(-1, 1, 50)
y = 1.0/(1+25*x**2)
plt.plot(x,y)

# plot the interpolation polynomial
a = coefficients(xData, yData)
xs=np.linspace(-1, 1, 50)
ys=np.zeros(50)
i=0
for x in xs:
    ys[i] = evalPoly(a,xData,x)
    i=i+1
plt.plot(xs,ys,'r--')

# plot the data points
numDataPoints = 12
xData = np.linspace(-1, 1, numDataPoints)
print xData
yData = 1.0/(1+25*xData**2)
print yData
plt.plot(xData,yData,'bo',markersize=10) 

# plot the exact function
x = np.linspace(-1, 1, 50)
y = 1.0/(1+25*x**2)
plt.plot(x,y)

# plot the interpolation polynomial
a = coefficients(xData, yData)
xs=np.linspace(-1, 1, 50)
ys=np.zeros(50)
i=0
for x in xs:
    ys[i] = evalPoly(a,xData,x)
    i=i+1
plt.plot(xs,ys,'r--')

# plot the data points
numDataPoints = 20
xData = np.linspace(-1, 1, numDataPoints)
print xData
yData = 1.0/(1+25*xData**2)
print yData
plt.plot(xData,yData,'bo',markersize=10) 

# plot the exact function
x = np.linspace(-1, 1, 50)
y = 1.0/(1+25*x**2)
plt.plot(x,y)

# plot the interpolation polynomial
a = coefficients(xData, yData)
xs=np.linspace(-1, 1, 50)
ys=np.zeros(50)
i=0
for x in xs:
    ys[i] = evalPoly(a,xData,x)
    i=i+1
plt.plot(xs,ys,'r--')