Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import numpy as np
- import matplotlib.pyplot as plt
- def CubicNatural(x, y):
- m = x.size # m is the number of data points
- n = m-1
- a = np.zeros(m)
- b = np.zeros(n)
- c = np.zeros(m)
- d = np.zeros(n)
- for i in range(m):
- a[i] = y[i]
- h = np.zeros(n)
- for i in range(n):
- h[i] = x[i+1] - x[i]
- u = np.zeros(n)
- u[0] = 0
- for i in range(1, n):
- u[i] = 3*(a[i+1]-a[i])/h[i]-3*(a[i]-a[i-1])/h[i-1]
- s = np.zeros(m)
- z = np.zeros(m)
- t = np.zeros(n)
- s[0] = 1
- z[0] = 0
- t[0] = 0
- for i in range(1, n):
- s[i] = 2*(x[i+1]-x[i-1])-h[i-1]*t[i-1]
- t[i] = h[i]/s[i]
- z[i]=(u[i]-h[i-1]*z[i-1])/s[i]
- s[m-1] = 1
- z[m-1] = 0
- c[m-1] = 0
- for i in np.flip(np.arange(n)):
- c[i] = z[i]-t[i]*c[i+1]
- b[i] = (a[i+1]-a[i])/h[i]-h[i]*(c[i+1]+2*c[i])/3
- d[i] = (c[i+1]-c[i])/(3*h[i])
- return a, b, c, d
- def CubicNaturalEval(w, x, coeff):
- m = x.size
- if w<x[0] or w>x[m-1]:
- print('error: spline evaluated outside its domain')
- return
- n = m-1
- p = 0
- for i in range(n):
- if w <= x[i+1]:
- break
- else:
- p += 1
- # p is the number of the subinterval w falls into, i.e., p=i means
- # w falls into the ith subinterval $(x_i,x_{i+1}), and therefore
- # the value of the spline at w is
- # a_i+b_i*(w-x_i)+c_i*(w-x_i)^2+d_i*(w-x_i)^3.
- a = coeff[0]
- b = coeff[1]
- c = coeff[2]
- d = coeff[3]
- return a[p]+b[p]*(w-x[p])+c[p]*(w-x[p])**2+d[p]*(w-x[p])**3
- xaxis = np.linspace(0, 14.875, 100)
- 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
- for i in CubicNatural(xi,yi):
- print(i)
- coeff = CubicNatural(xi, yi)
- naturalspline = np.array(list(map(lambda x: CubicNaturalEval(x, xi, coeff), xaxis)))
- plt.plot(xi, yi, 'o', label='Data Points')
- plt.plot(xaxis, naturalspline, label='Natural cubic spline')
- plt.xlabel("X")
- plt.ylabel("Y")
- plt.title("Natural Cubic Spline")
- plt.legend()
- plt.grid(True)
- plt.show()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement