为何三次样条算法在参数曲线上生成尖刺状错误图形？

我最近实现了Wikipedia上描述的三次样条算法，下面是我当前计算系数的实现：

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[i] 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[i] 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[i])*(a[i+1]-a[i])-(3.0/h[i-1])*(a[i]-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[i] = 2*(x_vals[i+1]-x_vals[i-1])-(h[i-1])*(mu[i-1]) # Stuff.
        mu[i] = h[i]/l[i]
        z[i] = (alpha[i]-h[i-1]*z[i-1])/l[i]
    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[i])
        splines[1].append(b[i])
        splines[2].append(c[i])
        splines[3].append(d[i])
    return splines # Return the output....

这个实现在某些情况下可以正常工作，我找到了Stack Overflow上的一个实现，交叉验证后发现它和我的函数生成的系数完全一致。

下面是一个演示脚本：

#!/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[i] and x <= intervals[i+1]:
            return i, intervals[i] # 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][i], splines[1][i], 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)

Stack Overflow上的实现放在名为plagiarized.py的文件中，脚本会调用其中的系数计算函数。USE_MODE变量用于选择系数计算的方法，代码中已有说明。

当设置USE_MODE == 1使用SciPy实现时，图形显示正常：

但当设置USE_MODE == 0使用Stack Overflow上的实现，或设置为其他值使用我自己的实现时，生成的图形会出现尖刺：

我尝试查看了GitHub上SciPy的三次样条实现，但没能完全理解。我已将相关文件上传到我的GitHub仓库，运行python spline_plot.py即可看到错误图形，修改USE_MODE为1就能看到预期的正确效果。

我的问题是：为什么我的代码（以及那个Stack Overflow上的实现）无法生成正确的图形？为什么图形会出现这么多尖刺？

提前感谢各位的解答！

备注：内容来源于stack exchange，提问作者Some nerd who does not have a