divideddifference.py

import numpy as np
import matplotlib.pyplot as plt

def diff(x, y):
    m = x.size # here m is the number of data points.
    # the degree of the polynomial is m-1
    a = np.zeros(m)
    for i in range(m):
        a[i] = y[i]
    for j in range(1, m):
        for i in np.flip(np.arange(j,m)):
            a[i] = (a[i]-a[i-1]) / (x[i]-x[i-(j)])
    return a

def newton(x, y, z):
    m = x.size # here m is the number of data points, not the degree
    # of the polynomial
    a = diff(x, y)
    sum = a[0]
    pr = 1.0
    for j in range(m-1):
        pr *= (z-x[j])
        sum += a[j+1]*pr
    return sum

xi = np.array([0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9, 9.5, 10, 10.5, 11, 11.5, 12, 12.5, 13, 13.5, 14, 14.5, 14.875])  # Replace with your data points
yi = np.array([0.44, 1.5, 2.875, 4.125, 5.5, 6.875, 8.25, 9.25, 9.375, 9.25, 8.75, 8, 7.375, 6.875, 6.375, 6.375, 6.625, 7.25, 8, 8.25, 7.5, 5.625, 5.25, 5, 5, 5.125, 5.625, 6.375, 7, 7.5, 8.25])  # Replace with your function values at x_data points

# print(diff(xi, yi))

xaxis = np.linspace(0, 14.875, 100)
interp = newton(xi, yi, xaxis)
plt.plot(xaxis, interp, label='Newton Divided Difference')
plt.plot(xi, yi, 'o', label='data')
plt.xlabel("X")
plt.ylabel("Y")
plt.title("Newton Divided Difference")
plt.legend()
plt.grid(True)
plt.show()