import numpy
import matplotlib.pyplot as plt
import math
import scipy.interpolate

def neville(x, y):
    n = len(x)
    gamma = y.copy()
    for k in range(1, n):
        for i in range(n-k-1, -1, -1):
            gamma[i+k] = (gamma[i+k] - gamma[i+k-1])/(x[i+k] - x[i])
    return gamma

def horner(x0, x, gamma):
    r = 0
    for k in range(len(x)-1, -1, -1):
        r = r*(x0-x[k]) + gamma[k];
    return r

def f(x):
    return numpy.sin(x)

x5 = numpy.linspace(0, 2.*math.pi, 5)
fgamma5 = neville(x5, f(x5))

x10 = numpy.linspace(0, 2.*math.pi, 10)
fgamma10 = neville(x10, f(x10))

x15 = numpy.linspace(0, 2.*math.pi, 15)
fgamma15 = neville(x15, f(x15))

x = numpy.linspace(0, 2.*math.pi, 100)
plt.figure()
plt.plot(
    x, f(x),
    x, horner(x, x5, fgamma5), 'r',
    x, horner(x, x10, fgamma10), 'b',
    x, horner(x, x15, fgamma15), 'g',
    x5, f(x5), 'ro',
    x10, f(x10), 'bx',
    x15, f(x15), 'g*',
    )
plt.legend(('f(x)', '5', '10', '15'))         

def g(x):
    return 1/(1+x*x)

x5 = numpy.linspace(-5, 5., 5)
ggamma5 = neville(x5, g(x5))

x10 = numpy.linspace(-5, 5, 10)
ggamma10 = neville(x10, g(x10))

x15 = numpy.linspace(-5, 5., 15)
ggamma15 = neville(x15, g(x15))

x = numpy.linspace(-5, 5, 100)
plt.figure()
plt.plot(
    x, g(x),
    x, horner(x, x5, ggamma5), 'r',
    x, horner(x, x10, ggamma10), 'b', 
    x, horner(x, x15, ggamma15), 'g',
    x5, g(x5), 'ro',
    x10, g(x10), 'bx',
    x15, g(x15), 'g*',
    )
plt.legend(('g(x)', '5', '10', '15'))         

def h(x):
    return x**-12 - x**-6

x5 = numpy.linspace(1, 5., 5)
hgamma5 = neville(x5, h(x5))

x10 = numpy.linspace(1, 5, 10)
hgamma10 = neville(x10, h(x10))

x15 = numpy.linspace(1, 5., 15)
hgamma15 = neville(x15, h(x15))

x = numpy.linspace(1, 5, 100)
plt.figure()
plt.plot(
    x, h(x),
    x, horner(x, x5, hgamma5), 'r',
    x, horner(x, x10, hgamma10), 'b', 
    x, horner(x, x15, hgamma15), 'g',
    x5, h(x5), 'ro',
    x10, h(x10), 'bx',
    x15, h(x15), 'g*',
    )
plt.legend(('h(x)', '5', '10', '15'))         

plt.show()