from numpy import *
from matplotlib.pyplot import *

class NewtonInterpolation:
    def __init__(self, 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])
        self.x = x
        self.gamma = gamma
    def __call__(self, x0):
        r = 0
        for k in range(len(self.x)-1, -1, -1):
            r = r*(x0-self.x[k]) + self.gamma[k];
        return r

def f(x):
    return sin(x)

x5 = linspace(0, 2.*pi, 5)
fp5 = NewtonInterpolation(x5, f(x5))

x10 = linspace(0, 2.*pi, 10)
fp10 = NewtonInterpolation(x10, f(x10))

x15 = linspace(0, 2.*pi, 15)
fp15 = NewtonInterpolation(x15, f(x15))

x = linspace(0, 2.*pi, 100)
figure()
plot(x, f(x), label='f(x)')
plot(x, fp5(x), label='5th degree')
plot(x, fp10(x), label='10th degree')
plot(x, fp15(x), label='15th degree')
legend()

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

x5 = linspace(-5., 5., 5)
gp5 = NewtonInterpolation(x5, g(x5))

x10 = linspace(-5., 5., 10)
gp10 = NewtonInterpolation(x10, g(x10))

x15 = linspace(-5., 5., 15)
gp15 = NewtonInterpolation(x15, g(x15))

x = linspace(-5, 5, 1000)
figure()
plot(x, g(x), label='g(x)')
plot(x, gp5(x), label='5th degree')
plot(x, gp10(x), label='10th degree')
plot(x, gp15(x), label='15th degree')
legend()

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

x5 = linspace(1, 5., 5)
hp5 = NewtonInterpolation(x5, h(x5))

x10 = linspace(1, 5, 10)
hp10 = NewtonInterpolation(x10, h(x10))

x15 = linspace(1, 5., 15)
hp15 = NewtonInterpolation(x15, h(x15))

x = linspace(1, 5, 1000)
figure()
plot(x, h(x), label='h(x)')
plot(x, hp5(x), label='5th degree')
plot(x, hp10(x), label='10th degree')
plot(x, hp15(x), label='15th degree')
legend()

show()