# divideddifference.py

May 22nd, 2024 (edited)
658
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. import numpy as np
2. import matplotlib.pyplot as plt
3.
4. def diff(x, y):
5.     m = x.size # here m is the number of data points.
6.     # the degree of the polynomial is m-1
7.     a = np.zeros(m)
8.     for i in range(m):
9.         a[i] = y[i]
10.     for j in range(1, m):
11.         for i in np.flip(np.arange(j,m)):
12.             a[i] = (a[i]-a[i-1]) / (x[i]-x[i-(j)])
13.     return a
14.
15. def newton(x, y, z):
16.     m = x.size # here m is the number of data points, not the degree
17.     # of the polynomial
18.     a = diff(x, y)
19.     sum = a[0]
20.     pr = 1.0
21.     for j in range(m-1):
22.         pr *= (z-x[j])
23.         sum += a[j+1]*pr
24.     return sum
25.
26. 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
27. 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
28.
29. # print(diff(xi, yi))
30.
31. xaxis = np.linspace(0, 14.875, 100)
32. interp = newton(xi, yi, xaxis)
33. plt.plot(xaxis, interp, label='Newton Divided Difference')
34. plt.plot(xi, yi, 'o', label='data')
35. plt.xlabel("X")
36. plt.ylabel("Y")
37. plt.title("Newton Divided Difference")
38. plt.legend()
39. plt.grid(True)
40. plt.show()