#331 Interpolation with Newton’s polynomial using python

interpolation.py

def interpolation(c,x,x0):
    # Evaluate Newton's polynomial at x0.
    # Degree of polynomial
    n = len(x) - 1
    y0 = c[n]
    for k in range(1,n+1):
        y0 = c[n-k] + (x0 - x[n-k]) * y0
    return y0

def coef(x,y):
    # Computes the coefficients of Newton's polynomial.
    # Number of data points.
    m = len(x)
    c = y.copy()
    for k in range(1,m):
        c[k:m] = (c[k:m] - c[k-1])/(x[k:m] - x[k-1])
    return c 

ex331.py

import numpy as np
import matplotlib.pyplot as plt
from interpolation import interpolation,coef

def f(x):
    return 1.0/(1.0+25.0*x**2)

a = -1.0
b = 1.0
x1 = np.linspace(a,b,200)
y1 = f(x1)
plt.plot(x1,y1)
nList = [4,6,10]
sglist = ['--','-.',':']
for k in range(len(nList)):
    n = nList[k]
    x = np.linspace(a,b,n+1)
    y = f(x)
    plt.scatter(x,y,marker='o',color='black')
    c = coef(x,y)
    y1 = interpolation(c,x,x1)
    sl = 'n=' + str(n)
    sg = sglist[k]
    plt.plot(x1,y1,sg,label=sl)
plt.xlabel('x')
plt.grid(True)    
plt.legend(loc=0)
plt.savefig('ex331.png', dpi=72)
plt.show()