[Python] Why does cubic spline algorithm produce spiky and incorrect graph on a parametric...

I recently implemented my implementation of the cubic spline algorithm described on wikipedia. Here is my current implementation which should calculate the coefficients:

def get_spline_natural_fifth_degree(x_vals, y_vals):
# This is taken straight from wikipedia: https://en.wikipedia.org/wiki/Spline_(mathematics)#Algorithm_for_computing_natural_cubic_splines
#x_vals = [p[0] for p in points]
#y_vals = [p[1] for p in points]
x_vals = list(x_vals)
y_vals = list(y_vals)
# 1. Create new array a of size n + 1 and for i = 0, …, n set a_i = y_i
n = len(x_vals)-1
a = [y_vals for i in range(len(y_vals))]# + [0.0] # Initialize the thing.
assert len(a) == n + 1
# a[-1] = 0.0 # Because the index
# 2. Create new arrays b and d, each of size n.
assert len(a) == n + 1
b = [0.0 for _ in range(n)]
d = [0.0 for _ in range(n)]
# 3. Create new array h of size n and for i = 0, …, n – 1 set h_i = x_(i+1) - x_i
h = [x_vals[i+1] - x_vals for i in range(n)]
# 4. Create new array α of size n and for i = 1, …, n – 1 set alpha_1 = (3/h_i)*(a_(i+1) - a_i) - (3/h_(i-1))*(a_i-a_(i-1))
alpha = [(3.0/h)*(a[i+1]-a)-(3.0/h[i-1])*(a-a[i-1]) for i in range(1,n)] # Actually n-1, but because python ranges are dumb, we need to do this.
#alpha.append(0.0)
alpha = [0.0] + alpha
assert len(alpha) == n
# 5. Create new arrays c, l, μ, z, each of size n + 1.
c = [0.0 for _ in range(n+1)]
assert len(c) == len(x_vals)
l = [0.0 for _ in range(n+1)]
mu = [0.0 for _ in range(n+1)]
z = [0.0 for _ in range(n+1)]
# 6. Set l_0 = 1 , mu_0 = z_0 = 0
l[0] = 1.0
mu[0] = 0.0
z[0] = 0.0
# 7. For i = 1 .. n-1 set the following: l_i = 2*(x_(i+1)-x_(i-1))-(h_(i-1))*(mu_(i-1)) mu_i = h_i/l_i z_i = (alpha_i-h_(i-1)*z_(i-1))/l_i
for i in range(1, n):
l = 2*(x_vals[i+1]-x_vals[i-1])-(h[i-1])*(mu[i-1]) # Stuff.
mu = h/l
z = (alpha-h[i-1]*z[i-1])/l
assert l[0] == 1.0
# 8. Set l_n = 1; z_n = c_n = 0
l[n] = 1.0
assert c[n] == 0.0 # Should be zero...
z[n] = 0.0
# 9. For j = n – 1, n – 2, …, 0, set the following: c_j = z_j - mu_j*c_(j+1) b_j = (a_(j+1)-a_j)/h_j - (h_j*(c_(j+1)+2*c_j))/3 and d_j = (c_(j+1)-c_j)/(3*h_j)
for j in range(n - 1, -1, -1):
c[j] = z[j] - mu[j]*c[j+1]
for j in range(n - 1, -1, -1):
b[j] = (a[j+1]-a[j])/h[j] + (h[j]*(2*c[j+1]+c[j]))/3.0
for j in range(n - 1, -1, -1):
d[j] = (c[j+1]-c[j])/(3.0*h[j])
a.pop(0)
c.pop(0)
splines = [[] for _ in range(4)]
for i in range(n):
splines[0].append(a)
splines[1].append(b)
splines[2].append(c)
splines[3].append(d)
return splines # Return the output....


and it works in some cases. I found this implementation on stackoverflow and I cross referenced it and it produces the exact same coefficients as my function.

Here is a script which demonstrates this:

#!/bin/python3

import matplotlib.pyplot as plt
import numpy as np
from spline import *
from scipy.interpolate import CubicSpline

def evaluate_spline(x, x_i, a, b, c, d):
return a + (b + (c + d * (x - x_i)) * (x - x_i)) * (x - x_i)

def correct_interval(x, intervals): # This returns the correct interval index and the x0 start value.
# Now loop over each element and check the thing.
for i in range(len(intervals)-1):
if x >= intervals and x <= intervals[i+1]:
return i, intervals # Return the index thing...
assert False

def calculate_spline(x, t_knots, splines):
i, t0 = correct_interval(x, t_knots)
#print("Correct interval: "+str(i))
if len(splines[0]) == 5:
a, b, c, d, _ = splines[0], splines[1], splines[2][i], splines[3][i], splines[4][i]
else:
a, b, c, d = splines[0][i], splines[1][i], splines[2][i], splines[3][i]# , splines[4][i]
return evaluate_spline(x, t0, a, b, c, d)


def compute_t_knots(points):
""" Computes t values based on cumulative distance along the curve. """
t_knots = [0]
for i in range(1, len(points)):
dist = np.sqrt((points[i][0] - points[i - 1][0])**2 + (points[i][1] - points[i - 1][1])**2)
t_knots.append(t_knots[-1] + dist) # Accumulate distances
return t_knots

def render_result(x_splines: list, y_splines: list, points: list, t_knots: list) -> None:
x_knots = [p[0] for p in points] # Something like this????
y_knots = [p[1] for p in points]
n = len(points)
t_values = np.linspace(min(t_knots), max(t_knots), 20000)
# Compute interpolated x and y values
x_things = [calculate_spline(t, t_knots, x_splines) for t in t_values]
y_things = [calculate_spline(t, t_knots, y_splines) for t in t_values]
plt.plot(x_things, y_things, label="Cubic Spline Curve")
plt.scatter(x_knots, y_knots, color="red", label="Control Points")
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.title('Spline Graph Using Manually Calculated Coefficients')
plt.grid(True)
plt.show()
return

from math import sqrt

def compute_t_range(points):
t_knots = [0]
for i in range(1, len(points)):
dist = sqrt((points[i][0] - points[i - 1][0])**2 + (points[i][1] - points[i - 1][1])**2)
t_knots.append(t_knots[-1] + dist) # Accumulate distances
return t_knots

from plagiarized import * # Import the implementation which I copied from stackoverflow. It is stored in plagiarized.py

USE_MODE = 0

if __name__=="__main__":
x_vals=[-0.83,0.14,-1.09,1.09,-0.54,2.03,3.0]
y_vals=[-2.03,-2.06,0.71,1.49,2.06,2.43,3.0]

assert len(x_vals) == len(y_vals)
assert all([isinstance(x, float) for x in x_vals])
assert all([isinstance(x, float) for x in y_vals])
points = [[x_vals[i], y_vals[i]] for i in range(len(x_vals))]
x_knots = [p[0] for p in points]
y_knots = [p[1] for p in points]
t_i = compute_t_knots(points)
assert len(t_i) == len(points) == len(x_knots) == len(y_knots)
points_x_t = [(t_i[i], x_knots[i]) for i in range(len(t_i))]
points_y_t = [(t_i[i], y_knots[i]) for i in range(len(t_i))]
spline_x_reference = CubicSpline(t_i, x_knots).c
spline_y_reference = CubicSpline(t_i, y_knots).c
spline_y = get_spline_natural_fifth_degree(t_i, y_knots)
spline_x = get_spline_natural_fifth_degree(t_i, x_knots)
spline_x_reference = list(spline_x_reference)
spline_y_reference = list(spline_y_reference)
spline_x_reference.reverse()
spline_y_reference.reverse()
spline_x_stuff = calc_spline_params(np.array(t_i), np.array(x_knots))
spline_y_stuff = calc_spline_params(np.array(t_i), np.array(y_knots))
if USE_MODE == 0: # Use the implementation which I took from stackoverflow
render_result(spline_x_stuff, spline_y_stuff, points, t_i)
elif USE_MODE == 1: # Use Scipy. This one actually displays the correct result.
render_result(spline_x_reference, spline_y_reference, points, t_i)
else: # Use my own implementation
render_result(spline_x, spline_y, points, t_i)
exit(0)



the implementation from the other post on stackoverflow is in a file called plagiarized.py and we call the coefficient calculation function from this script. The USE_MODE variable describes which method to use in calculation of the coefficients as explained in the code.

The graph displays correctly like so:

[URL='https://i.sstatic.net/8MJfNCuT.png'][IMG]https://i.sstatic.net/8MJfNCuT.png[/IMG][/URL]

when running with USE_MODE == 1 meaning that the scipy implementation should be used. However, when running in either USE_MODE == 0 meaning the copied implementation from stackoverflow or when USE_MODE is something else meaning to use my own implementation, the following graph gets produced:

[URL='https://i.sstatic.net/Tp2ctvFJ.png'][IMG]https://i.sstatic.net/Tp2ctvFJ.png[/IMG][/URL]

I tried looking at the scipy implementation of cubic splines [URL='https://github.com/scipy/scipy/blob/main/scipy/interpolate/_cubic.py#L738']here on github[/URL], but I wasn't able to really understand it.

I have uploaded the relevant files to [URL='https://github.com/personnumber3377/spline_debugging']my github[/URL]. Just run python spline_plot.py to observe the wrong graph. Modify the value of USE_MODE to 1 and you can see what it is supposed to look like.

So, my question is why doesn't my code (or the other implementation for that matter) produce the same graph and why is the graph so spiky?

Thanks in advance for your answers!

[url="https://stackoverflow.com/questions/79410158/why-does-cubic-spline-algorithm-produce-spiky-and-incorrect-graph-on-a-parametri"]Continue reading...[/url][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i]