import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
x = [1, 2, 3]
y = [6, 11, 18]
plt.plot(x,y,'bo',markersize=10)
plt.xlim(0, 4)
plt.ylim(5, 19)
x=np.linspace(0, 4, 20)
y=x**2+2*x+3
plt.plot(x, y, 'r--')

import numpy as np

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

# y = x**2 + 2*x + 3
xData = [1, 2, 3]
yData = [6, 11, 18]
a = coefficients(xData, yData)
print a
print " x    yInterp   yExact"
print "-----------------------"
for x in np.arange(0.0,4,0.2):
    y = evalPoly(a,xData,x)
    yExact = y = x**2 + 2*x + 3
    print "%3.1f %9.5f %9.5f"% (x,y,yExact)

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

xAll = [0, 10, 15, 20, 22.5, 30]
yAll = [0, 227.04, 362.78, 517.35, 602.97, 901.67]
plt.plot(xAll,yAll,'bo',markersize=10)

xData = xAll[2:4]
print "xData=",xData
yData = yAll[2:4]
print "yData=",yData
a = coefficients(xData, yData)
print "a=",a
y = evalPoly(a, xData, 15)
print "f(15)=", y
y = evalPoly(a, xData, 16)
print "f(16)=", y

xs=np.linspace(0, 30, 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--')

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

xAll = [0, 10, 15, 20, 22.5, 30]
yAll = [0, 227.04, 362.78, 517.35, 602.97, 901.67]
plt.plot(xAll,yAll,'bo',markersize=10)

xData = xAll[1:4]
print "xData=",xData
yData = yAll[1:4]
print "yData=",yData
a = coefficients(xData, yData)
print "a=",a
y = evalPoly(a, xData, 10)
print "f(10)=", y
y = evalPoly(a, xData, 15)
print "f(15)=", y
y = evalPoly(a, xData, 16)
print "f(16)=", y

xs=np.linspace(0, 30, 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--')

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

xAll = [0, 10, 15, 20, 22.5, 30]
yAll = [0, 227.04, 362.78, 517.35, 602.97, 901.67]
plt.plot(xAll,yAll,'bo',markersize=10)

xData = xAll
print "xData=",xData
yData = yAll
print "yData=",yData
a = coefficients(xData, yData)
print "a=",a

y = evalPoly(a, xData, 16)
print "f(16)=", y

xs=np.linspace(0, 30, 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--')

# add the data points (0,8) (8,0) (3,2) (4,5) (1,12)
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

xAll = [8, 4, 0, 1, 3]
yAll = [0, 5, 8, 12, 2]
plt.plot(xAll,yAll,'bo',markersize=10)

xData = xAll
print "xData=",xData
yData = yAll
print "yData=",yData
a = coefficients(xData, yData)
print "a=",a


xs=np.linspace(-0.2, 8.2, 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--')