Advertisement
Hustinyano

naturalcubicspline.py

May 22nd, 2024
698
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 2.41 KB | None | 0 0
  1. import numpy as np
  2. import matplotlib.pyplot as plt
  3.  
  4. def CubicNatural(x, y):
  5.     m = x.size # m is the number of data points
  6.     n = m-1
  7.     a = np.zeros(m)
  8.     b = np.zeros(n)
  9.     c = np.zeros(m)
  10.     d = np.zeros(n)
  11.     for i in range(m):
  12.         a[i] = y[i]
  13.     h = np.zeros(n)
  14.     for i in range(n):
  15.         h[i] = x[i+1] - x[i]
  16.     u = np.zeros(n)
  17.     u[0] = 0
  18.     for i in range(1, n):
  19.         u[i] = 3*(a[i+1]-a[i])/h[i]-3*(a[i]-a[i-1])/h[i-1]
  20.     s = np.zeros(m)
  21.     z = np.zeros(m)
  22.     t = np.zeros(n)
  23.     s[0] = 1
  24.     z[0] = 0
  25.     t[0] = 0
  26.     for i in range(1, n):
  27.         s[i] = 2*(x[i+1]-x[i-1])-h[i-1]*t[i-1]
  28.         t[i] = h[i]/s[i]
  29.         z[i]=(u[i]-h[i-1]*z[i-1])/s[i]
  30.     s[m-1] = 1
  31.     z[m-1] = 0
  32.     c[m-1] = 0
  33.     for i in np.flip(np.arange(n)):
  34.         c[i] = z[i]-t[i]*c[i+1]
  35.         b[i] = (a[i+1]-a[i])/h[i]-h[i]*(c[i+1]+2*c[i])/3
  36.         d[i] = (c[i+1]-c[i])/(3*h[i])
  37.     return a, b, c, d
  38.  
  39. def CubicNaturalEval(w, x, coeff):
  40.     m = x.size
  41.     if w<x[0] or w>x[m-1]:
  42.         print('error: spline evaluated outside its domain')
  43.         return
  44.     n = m-1
  45.     p = 0
  46.     for i in range(n):
  47.         if w <= x[i+1]:
  48.             break
  49.         else:
  50.             p += 1
  51.     # p is the number of the subinterval w falls into, i.e., p=i means
  52.     # w falls into the ith subinterval $(x_i,x_{i+1}), and therefore
  53.     # the value of the spline at w is
  54.     # a_i+b_i*(w-x_i)+c_i*(w-x_i)^2+d_i*(w-x_i)^3.
  55.     a = coeff[0]
  56.     b = coeff[1]
  57.     c = coeff[2]
  58.     d = coeff[3]
  59.     return a[p]+b[p]*(w-x[p])+c[p]*(w-x[p])**2+d[p]*(w-x[p])**3
  60.  
  61. xaxis = np.linspace(0, 14.875, 100)
  62.  
  63. 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
  64. 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
  65.  
  66. for i in CubicNatural(xi,yi):
  67.     print(i)
  68.  
  69. coeff = CubicNatural(xi, yi)
  70. naturalspline = np.array(list(map(lambda x: CubicNaturalEval(x, xi, coeff), xaxis)))
  71.  
  72. plt.plot(xi, yi, 'o', label='Data Points')
  73. plt.plot(xaxis, naturalspline, label='Natural cubic spline')
  74. plt.xlabel("X")
  75. plt.ylabel("Y")
  76. plt.title("Natural Cubic Spline")
  77. plt.legend()
  78. plt.grid(True)
  79. plt.show()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement