Advertisement
Hustinyano

divideddifference.py

May 22nd, 2024 (edited)
658
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 1.35 KB | None | 0 0
  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()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement