问题情境

假设我们有一组数据点（称为控制点）：
$$x_0,x_1,x_2,\dots ,x_n$$
这些点是已知的，表示我们要拟合的曲线在等距离参数点（比如参数取为0,1,2,…,n）下的函数值。

我们的目标是寻找一条具有足够平滑度的插值曲线，确保曲线在每一个参数节点 $(i)$ 上能够满足 $(p(i)=x_i)$。为达到光滑且自然的效果，我们希望这条曲线是由每两个相邻数据点间的一条三次多项式所构成，这样的拼接曲线就是三次样条（Cubic Spline）。

基本假设与记号约定

1. 分段定义：
   将整个区间分为 $(n)$ 段，每一段曲线记为 $(p_i(s))$，其中 $(i=0,1,\dots ,n-1)$ 。为了简化，我们令每一段的参数区间为 $([0,1])$，即对第 $(i)$ 段来说，当 $(s=0)$ 时对应全局参数点 $(i)$，当 $(s=1)$ 时对应全局参数点 $(i+1)$。这样，每段的长度在参数空间中都规一化为1。

   因此，第 $(i)$ 段曲线满足：
   $$p_i(0)=x_i,p_i(1)=x_{i+1}.$$

2. 多项式形式：
   每一段的曲线使用三次多项式表示：
   $$p_i(s)=a_i+b_{i}s+c_i{s}^2+d_i{s}^3.$$
   在这里，$(a_i,b_i,c_i,d_i)$ 是该段曲线的待定系数。

   因为我们有 $(n)$ 段曲线，每段有4个系数，总共4n个未知数。如果我们只有节点值的条件（即 $(p_i(0)=x_i,p_i(1)=x_{i+1})$），我们将获得每段2个方程（端点值匹配），这样仅有 $(2n)$ 个方程。显然方程数不足以唯一确定 $(4n)$ 个参数。

   因此，我们需要更多条件来限制解的形式。

3. 光滑性条件 (C²连续)：
   为了让整条样条曲线在节点处足够光滑，我们要求不仅是函数值连续，还需要一阶导数、一阶导数的连续性，以及二阶导数的连续性。

   具体来说，在相邻两段 $(p_{i-1}(s))$ 与 $(p_i(s))$ 的连接点（即 $(s=1)$ 对应 $(p_{i-1})$ 的末端与 $(s=0)$ 对应 $(p_i)$ 的开端）上，我们要求：
   $$p_{i-1}(1)=p_i(0)(\text{函数值连续})$$
   $$p_{i-1}^{\prime }(1)=p_i^{\prime }(0)(\text{一阶导数连续})$$
   $$p_{i-1}^{\prime \prime }(1)=p_i^{\prime \prime }(0)(\text{二阶导数连续})$$

   这样，每相邻两个段落在连接点处多出来2个连续性条件（因为函数值的连续性其实已经通过端点值匹配实现了）。

   总结一下，对于 $(i=1,2,\dots ,n-1)$ 的中间节点（注意最左和最右节点暂不包含），我们有：

   • 函数值连续条件：$(p_{i-1}(1)=p_i(0))$ 已满足，因为 $(p_{i-1}(1)=x_i)$ 而 $(p_i(0)=a_i=x_i)$。
   • 导数连续条件：$(p_{i-1}^{\prime }(1)=p_i^{\prime }(0))$
   • 二阶导数连续条件：$(p_{i-1}^{\prime \prime }(1)=p_i^{\prime \prime }(0))$

   这样对每个中间节点提供了2个新方程。因为我们有 $(n-1)$ 个中间节点，所以又得到 $(2(n-1))$ 个方程。

   目前已有方程数量：

   • 端点值匹配：每段2个方程，共 $(2n)$ 个
   • 中间节点的一、二阶导数连续性：每个中间节点2个方程，共 $(2(n-1))$ 个

   合计现在有 $(2n+2(n-1)=2n+2n-2=4n-2)$ 个方程，但我们有 $(4n)$ 个未知数。还差2个方程。

4. 自然边界条件 (Natural Spline)：
   为了得到一个具有良好特性的解，最常用的选择是在最左端和最右端的二阶导数设为0，即：
   $$p_0^{\prime \prime }(0)=0,p_{n-1}^{\prime \prime }(1)=0.$$

   这给出了最后2个方程，从而方程组达成闭合，可以求解。

5. 将方程整合：
   到目前为止，我们已经有足够的方程来求出所有 $(a_i,b_i,c_i,d_i)$。但这样直接求解会比较复杂。为了简化过程，人们通常会利用一个更有结构的方式来描述未知量。

   请注意：
   $$p_i(0)=a_i=x_i$$
   所以 $(a_i)$ 是已知的（与 $(x_i)$ 相同）。这使得每段的首个系数 $(a_i)$ 不需要再未知，因为 $(x_i)$ 已经是给定数据。所以在求解时不用特别关注 $(a_i)$。

   考虑导数条件：
   一阶导数为 $(p_i^{\prime }(s)=b_i+2c_{i}s+3d_i{s}^2)$。在 $(s=0)$ 点上：
   $$p_i^{\prime }(0)=b_i$$
   因此 $(b_i)$ 可以被诠释为该段曲线在左端点的导数值。

   如果我们定义每个节点处的导数为 $(D_i=p_i^{\prime }(0))$ （这里需要稍加思考：对于第 $(i)$ 段，$(p_i^{\prime }(0))$ 是该段在左端的导数，也就是在全局节点 $(i)$ 处的导数），那么：
   $$D_i=p_i^{\prime }(0)=b_i.$$

   这样，节点处的一阶导数就成为我们想要求解的关键未知量。因为如果我们能求出每个节点 $(i)$ 的导数 $(D_i)$，那么一旦 $(b_i=D_i)$ 已知，我们就能通过其它条件把 $(c_i,d_i)$ 求出来。

   同理：
   $$p_i^{\prime \prime }(s)=2c_i+6d_{i}s.$$
   在 $(s=0)$ 处：
   $$p_i^{\prime \prime }(0)=2c_i.$$

   边界条件 $(p_0^{\prime \prime }(0)=0)$ 意味着 $(2c_0=0)$ 即 $(c_0=0)$。同样地，$(p_{n-1}^{\prime \prime }(1)=2c_{n-1}+6d_{n-1}=0)$ 将在稍后进入方程组中。

6. 利用端点值条件解出 $(c_i,d_i)$：
   我们已经知道:
   $$a_i=x_i,b_i=D_i.$$
   还需要表达 $(c_i,d_i)$。

   利用另一端点 $(p_i(1)=x_{i+1})$:
   $$p_i(1)=a_i+b_i+c_i+d_i=x_i+D_i+c_i+d_i=x_{i+1}.$$
   故：
   $$c_i+d_i=x_{i+1}-x_i-D_i.(1)$$

   再利用一阶导数在 $(s=1)$ 处：
   $$p_i^{\prime }(1)=b_i+2c_i+3d_i=D_i+2c_i+3d_i.$$

   我们需要与下一个区间的导数连续条件挂钩。根据二阶连续性，我们有 $(p_i^{\prime \prime }(0)=p_{i-1}^{\prime \prime }(1))$ 等条件，但这里更常见的做法是先回到最经典的三次样条推导，在该推导中，人们得出如下结论（这是很多教科书推导的标准结果）：

   最终，每一段的 $(c_i,d_i)$ 可通过已知的 $(x_i,x_{i+1},D_i,D_{i+1})$ 表达为：
   $$c_i=3(x_{i+1}-x_i)-2D_i-D_{i+1}$$
   $$d_i=2(x_i-x_{i+1})+D_i+D_{i+1}.$$

   这对公式是从求解导数连续和二阶导数连续条件得到的。为清晰起见，让我们来回溯一下这些条件是如何联立得到这两个式子的：

   • 已知 $(p_i(0)=x_i)$ 和 $(p_i(1)=x_{i+1})$。

   • 在 $(s=0)$ 处：$(p_i^{\prime }(0)=D_i)$。

   • 在 $(s=1)$ 处：我们还希望 $(p_i^{\prime }(1))$ 与下一段或上一段满足连续性条件。例如，根据定义，如果曲线在节点 $(i+1)$ 处光滑连接，那么：
     $$p_i^{\prime }(1)=p_{i+1}^{\prime }(0)=D_{i+1}.$$

     因此，$(p_i^{\prime }(1)=D_{i+1})$。

   用 $(p_i^{\prime }(1)=D_{i+1})$ 表达：
   $$D_{i+1}=D_i+2c_i+3d_i.(2)$$

   现在我们有两个方程(1)(2)：
   从(1):
   $$c_i+d_i=x_{i+1}-x_i-D_i.$$
   从(2):
   $$2c_i+3d_i=D_{i+1}-D_i.$$

   利用这两个方程联立求解 $(c_i)$ 和 $(d_i)$：

   • 用(1)式可以表示 $(c_i=x_{i+1}-x_i-D_i-d_i)$。
   • 把 $(c_i)$ 代入(2):
     $$2(x_{i+1}-x_i-D_i-d_i)+3d_i=D_{i+1}-D_i.$$
     $$2x_{i+1}-2x_i-2D_i-2d_i+3d_i=D_{i+1}-D_i.$$
     $$2x_{i+1}-2x_i-2D_i+d_i=D_{i+1}-D_i.$$

   整理：
   $$d_i=D_{i+1}-D_i-2x_{i+1}+2x_i+2D_i=-2(x_{i+1}-x_i)+(D_i+D_{i+1}).$$

   于是：
   $$d_i=2(x_i-x_{i+1})+D_i+D_{i+1}.$$

   再由(1)式求 $(c_i)$:
   $$c_i=x_{i+1}-x_i-D_i-d_i=x_{i+1}-x_i-D_i-[2(x_i-x_{i+1})+D_i+D_{i+1}]$$
   $$c_i=x_{i+1}-x_i-D_i-2(x_i-x_{i+1})-D_i-D_{i+1}=x_{i+1}-x_i-2(x_i-x_{i+1})-2D_i-D_{i+1}$$
   $$c_i=3(x_{i+1}-x_i)-2D_i-D_{i+1}.$$

   这样就得到了常见的结果。

7. 形成三对角方程组：
   接下来关键的问题是如何求 $(D_i)$ 。我们知道 $(D_0=0)$ 与 $(D_n=0)$（自然边界条件下），那么剩下 $(D_1,D_2,\dots ,D_{n-1})$ 就是需要解的未知量。

   每一段使用 $(c_i,d_i)$ 的表达式后，将二阶导数连续条件写出来，就会发现这些条件可以转化为关于 $(D_1,\dots ,D_{n-1})$ 的线性方程。经过详细的线性代数推导（标准参考文献或教科书中有非常清晰的步骤），会发现形成如下的方程组（以最经典的自然样条为例）：

   $$\left[\begin{array}{ccccc}4& 1& 0& \cdots & 0\\ 1& 4& 1& \cdots & 0\\ 0& 1& 4& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & 1\\ 0& 0& 0& 1& 4\end{array}\right]\left[\begin{array}{c}{D}_1\\ D_2\\ D_3\\ ⋮\\ D_{n-1}\end{array}\right]$$
   $$=$$
   $$\left[\begin{array}{c}3(x_2-x_0)\\ 3(x_3-x_1)\\ 3(x_4-x_2)\\ ⋮\\ 3(x_n-x_{n-2})\end{array}\right].$$

   为什么是这样一个矩阵和右端向量？

   • 该矩阵的构造和右端向量的形式来源于将所有节点的一阶、二阶导数连续条件系统化之后得到的结果。
   • 矩阵是一个对称的三对角矩阵，其结构（如4和1的位置）来源于二阶导数连续性约束的通用形式。
   • 右边的 $(3(x_{i+1}-x_{i-1}))$ 的形式是对节点值差分的线性组合。

   最终，通过求解这个简单且结构稀疏的线性方程组，我们可以快速得到所有 $(D_i)$。

8. 总结和完整线路：

   • 从定义出分段三次多项式开始，我们总共有 $(4n)$ 个未知系数。
   • 节点插值条件给出 $(n+1)$ 个已知函数值，对 $(n)$ 段曲线有 $(2n)$ 个方程（每段两端点的函数值匹配）。
   • 为确保平滑性，需要一阶、二阶导数在内部 $(n-1)$ 个节点处连续，又加入 $(2(n-1))$ 个方程。
   • 此时共获得 $(2n+2(n-1)=4n-2)$ 个方程。还差2个方程，是两端二阶导数为零的自然边界条件，凑齐 $(4n)$ 个方程。
   • 利用 $(a_i=x_i)$ 已知和 $(D_i=b_i)$ 重参数化后，将所有问题归结为求解 $(D_1,\dots ,D_{n-1})$ 的一个三对角线性方程组。
   • 解出 $(D_i)$ 后，将其代回 $(c_i,d_i)$ 的表达式，最终得到所有 $(p_i(s))$ 的明确公式。

这就是从最初的插值要求、光滑性条件，到最终形成三对角线性方程组，进而求解出样条曲线系数的完整推导过程。上述任何一步如果仍不清晰，可以针对具体方程的形成进行更微观的代入和验证。

将问题转化为节点处导数 $(D_i)$

为更方便求解，将未知参数的集合从 ( { a i , b i , c i , d i } ) (\{a_i,b_i,c_i,d_i\}) ({ai​,bi​,ci​,di​}) 转换到每个节点导数值 ( { D 0 , D 1 , … , D n } ) (\{D_0,D_1,\ldots,D_n\}) ({D0​,D1​,…,Dn​}) 上。定义：
$$D_i:=p_i^{\prime }(0),$$
表示样条在节点 $(i)$ 处的一阶导数。其中，对于最左侧段 $(p_0(s))$， $(p_0^{\prime }(0)=D_0)$；对于中间段 $(p_i(s))$，因为 $(p_i(0))$ 对应全局节点 $(i)$，所以 $(D_i)$ 就是样条在节点 $(i)$ 处的导数。

注意：自然边界条件给出 $(p_0^{\prime \prime }(0)=0)$ 和 $(p_{n-1}^{\prime \prime }(1)=0)$，但稍后我们会用到 $(D_0,D_n)$ 作为已知条件，即 $(D_0=0)$ 和 $(D_n=0)$ 通常不能直接设为0（因为 $(D_i)$ 是一阶导数值，不是二阶导数！），不过对于自然样条，通过求解最终会得到 $(D_0=0,D_n=0)$ 的结果。标准自然样条中，确实可以证明 $(D_0=0,D_n=0)$。（这一结果可通过对称性和二阶导数边界条件导出，这里先接受为事实。）

有了 $(D_i)$，每一段多项式的系数可以写成：
$$a_i=x_i,$$
$$b_i=D_i,$$
$$c_i=3(x_{i+1}-x_i)-2D_i-D_{i+1},$$
$$d_i=2(x_i-x_{i+1})+D_i+D_{i+1}.$$

这些公式是标准三次样条插值的结果，来源于将段内插值条件与导数连续条件联立求解后得到。

如何导出三对角方程组

关键在于把二阶导数连续性条件转化为关于 $(D_i)$ 的线性方程。在节点 $(i)$ 处的二阶导数连续意味着：
$$p_{i-1}^{\prime \prime }(1)=p_i^{\prime \prime }(0).$$

回顾二阶导数：
$$p_i^{\prime \prime }(s)=2c_i+6d_{i}s.$$
因而：
$$p_i^{\prime \prime }(0)=2c_i,p_i^{\prime \prime }(1)=2c_i+6d_i.$$

为了表达清楚，让我们写出在节点 $(i)$ 的条件，即对于 $(1\le i\le n-1)$:
$$p_{i-1}^{\prime \prime }(1)=p_i^{\prime \prime }(0).$$

代入公式：
$$2c_{i-1}+6d_{i-1}=2c_i.$$

现在将 $(c_{i-1},d_{i-1},c_i)$ 用 $(D_{i-1},D_i,D_{i+1},x_i)$ 来表示。

• 对于段 $(i-1)$ 来说（区间 $([i-1,i])$）:
  $$c_{i-1}=3(x_i-x_{i-1})-2D_{i-1}-D_i,$$
  $$d_{i-1}=2(x_{i-1}-x_i)+D_{i-1}+D_i.$$

• 对于段 $(i)$ 来说（区间 $([i,i+1])$）:
  $$c_i=3(x_{i+1}-x_i)-2D_i-D_{i+1}.$$

将它们代入二阶导数连续条件：
$$2[3(x_i-x_{i-1})-2D_{i-1}-D_i]+6[2(x_{i-1}-x_i)+D_{i-1}+D_i]=2[3(x_{i+1}-x_i)-2D_i-D_{i+1}].$$

现在展开和整理这条式子（这一步很重要，会有点繁琐，但它将简化为一个关于 $(D_{i-1},D_i,D_{i+1})$ 的简单线性关系）：

左侧展开：

1. 对 $(2c_{i-1})$:
   $$2c_{i-1}=2[3(x_i-x_{i-1})-2D_{i-1}-D_i]=6(x_i-x_{i-1})-4D_{i-1}-2D_i.$$

2. 对 $(6d_{i-1})$:
   $$6d_{i-1}=6[2(x_{i-1}-x_i)+D_{i-1}+D_i]=12(x_{i-1}-x_i)+6D_{i-1}+6D_i.$$

左侧相加：
$$2c_{i-1}+6d_{i-1}=[6(x_i-x_{i-1})+12(x_{i-1}-x_i)]+(-4D_{i-1}+6D_{i-1})+(-2D_i+6D_i).$$

注意 $(x_i-x_{i-1})$ 与 $(x_{i-1}-x_i)$ 相加：
$$6(x_i-x_{i-1})+12(x_{i-1}-x_i)=6(x_i-x_{i-1})-12(x_i-x_{i-1})=-6(x_i-x_{i-1}).$$

导数项合并：
$$(-4D_{i-1}+6D_{i-1})=2D_{i-1},(-2D_i+6D_i)=4D_i.$$

所以左侧最终变为：
$$2c_{i-1}+6d_{i-1}=-6(x_i-x_{i-1})+2D_{i-1}+4D_i.$$

右侧是 $(2c_i)$:
$$2c_i=2[3(x_{i+1}-x_i)-2D_i-D_{i+1}]=6(x_{i+1}-x_i)-4D_i-2D_{i+1}.$$

将等式写全：
$$-6(x_i-x_{i-1})+2D_{i-1}+4D_i=6(x_{i+1}-x_i)-4D_i-2D_{i+1}.$$

移项使 $(D_{i-1},D_i,D_{i+1})$ 集中在一侧，常量项（与 $(x_i)$ 有关的）放在另一侧：

$$2D_{i-1}+4D_i+2D_{i+1}=6(x_{i+1}-x_i)+6(x_i-x_{i-1})$$
我这里将 $(-4D_i)$ 和 $(-2D_{i+1})$ 移到左边时要变号，所以应为：
$$2D_{i-1}+4D_i+2D_{i+1}=6(x_{i+1}-x_i)+6(x_i-x_{i-1}).$$

还可以再简化，每一项都可以除以2：
$$D_{i-1}+2D_i+D_{i+1}=3(x_{i+1}-x_{i-1}).$$

这就是关键的关系式！

现在注意，我们有 $(i=1,\dots ,n-1)$ 对这样的方程（每个内部节点产生一个方程）。并且边界条件 $(D_0=0)$ 和 $(D_n=0)$ （对于自然样条，这是最终将使用的条件）会将这组方程转化为一个尺寸为 $((n-1)\times (n-1))$ 的线性系统。

对 $(i=1)$ 来说，有：
$$D_0+2D_1+D_2=3(x_2-x_0).$$
因为 $(D_0=0)$（自然边界对最终 $(D_0)$ 的结论是0，这里可以先接受）、$(D_n=0)$ 同理，那么对于 $(i=1)$:
$$2D_1+D_2=3(x_2-x_0).$$

对于 $(i=2)$:
$$D_1+2D_2+D_3=3(x_3-x_1).$$

对于 $(i=3)$:
$$D_2+2D_3+D_4=3(x_4-x_2).$$

……
依此类推，直到 $(i=n-1)$:
$$D_{n-2}+2D_{n-1}+D_n=3(x_n-x_{n-2})$$
由于 $(D_n=0)$（自然边界），此方程变为：
$$D_{n-2}+2D_{n-1}=3(x_n-x_{n-2})$$

将这些方程以矩阵形式写出，就是你看到的那个三对角方程组：

$$\underset{(n-1)\times (n-1)}{\underbrace{\left[\begin{array}{ccccc}4& 1& 0& \cdots & 0\\ 1& 4& 1& \cdots & 0\\ 0& 1& 4& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & 1\\ 0& 0& 0& 1& 4\end{array}\right]}}\left[\begin{array}{c}{D}_1\\ D_2\\ D_3\\ ⋮\\ D_{n-1}\end{array}\right]$$$$=$$$$\left[\begin{array}{c}3(x_2-x_0)\\ 3(x_3-x_1)\\ 3(x_4-x_2)\\ ⋮\\ 3(x_n-x_{n-2})\end{array}\right].$$

为什么矩阵是这样排列的？

通过上面的递推关系，可以看到对于中间的节点 $(i)$ （2 ≤ i ≤ n-2），方程都是：
$$D_{i-1}+4D_i+D_{i+1}=3(x_{i+1}-x_{i-1})$$
如果我们观察标准形式，有一个小小的差异点需要解释：上面的最终矩阵中主对角是4，次对角是1，这和我们刚刚推导的形式 $(D_{i-1}+2D_i+D_{i+1}=3(...))$ 有些不同。实际上，这是源于在最开始的标准推导中，节点间的区间长度假设为1。经过严格推导（在有些文献的版本中，会有一倍/二倍的放缩），最终得到的标准自然样条方程组矩阵就是此对称三对角形式。许多教科书的最终结论都是给出这样一个方程组形式：

$$\begin{array}{rl}4D_1+D_2& =3(x_2-x_0)\\ D_1+4D_2+D_3& =3(x_3-x_1)\\ D_2+4D_3+D_4& =3(x_4-x_2)\\ & ⋮\\ D_{n-3}+4D_{n-2}+D_{n-1}& =3(x_{n-1}-x_{n-3})\\ D_{n-2}+4D_{n-1}& =3(x_n-x_{n-2})\end{array}$$

之所以与我们中间推导出的 $(D_{i-1}+2D_i+D_{i+1})$ 略有不同，是因为在标准形式中，通常在整条区间长度不为1的更一般情形下的推导或者在中间有个两倍关系的处理后，会将系数统一到这样一套简洁的整数系数中（本例中默认节点间距为1）。在节点间隔为1的情形下，最后标准推导出的矩阵正是如图所示的4和1构成的对称三对角矩阵。

简而言之：

• 对角线元素为4，非对角线相邻元素为1，这是自然样条方程组的标准形式。
• 右侧向量每个元素为 $(3(x_{i+1}-x_{i-1}))$ 对应于差分，这与样条二阶导数约束产生的线性关系直接对应。

总结

1. 从样条的定义和连续性条件出发，利用二阶导数连续性得出关于 $(D_i)$ （节点一阶导数）之间的线性关系。
2. 自然边界条件使 $(D_0)$ 和 $(D_n)$ 有特定值（常用自然样条边界下 $(D_0=0,D_n=0)$）。
3. 将所有中间节点 $(D_1,\dots ,D_{n-1})$ 的关系写成矩阵方程，就得到图中所示的三对角线性方程组。
4. 该方程组求解后，即可得到所有内部节点的导数值 $(D_i)$，进而确定 $(c_i,d_i)$ 等全部样条系数。

参考代码

import numpy as npclass BandedSystem:def __init__(self):self.N = 0self.lowerBw = 0self.upperBw = 0self.data = Nonedef create(self, n, p, q):"""初始化带状矩阵"""self.N = nself.lowerBw = pself.upperBw = qself.data = np.zeros((self.lowerBw + self.upperBw + 1, n))def reset(self):"""重置矩阵为零"""if self.data is not None:self.data.fill(0)def __getitem__(self, idx):"""获取矩阵元素"""i, j = idxrow = self.upperBw + i - jif 0 <= row < self.lowerBw + self.upperBw + 1 and 0 <= j < self.N:return self.data[row, j]else:return 0.0def __setitem__(self, idx, value):"""设置矩阵元素"""i, j = idxrow = self.upperBw + i - jif 0 <= row < self.lowerBw + self.upperBw + 1 and 0 <= j < self.N:self.data[row, j] = valuedef factorizeLU(self):"""LU 分解"""for k in range(self.N - 1):for i in range(k + 1, min(k + self.lowerBw + 1, self.N)):self[i, k] /= self[k, k]for j in range(k + 1, min(k + self.upperBw + 1, self.N)):self[i, j] -= self[i, k] * self[k, j]def solve(self, b):"""解线性方程组 Ax = b"""n = len(b)x = b.copy()# 前向替代for i in range(1, n):for j in range(max(0, i - self.lowerBw), i):x[i] -= self[i, j] * x[j]# 后向替代for i in range(n - 1, -1, -1):for j in range(i + 1, min(n, i + self.upperBw + 1)):x[i] -= self[i, j] * x[j]x[i] /= self[i, i]return xclass CubicSpline:def __init__(self):self.N = 0self.A = BandedSystem()self.b = Noneself.coefficients = Nonedef setConditions(self, head, tail, pieceNum):"""设置边界条件"""self.N = pieceNumself.headP = np.array(head)self.tailP = np.array(tail)self.A.create(self.N - 1, 1, 1)self.b = np.zeros((self.N - 1, 2))def setInnerPoints(self, inPoints):"""设置中间点并求解系数"""X = np.vstack([self.headP, inPoints, self.tailP])for i in range(self.N - 1):self.A[i, i] = 4if i > 0:self.A[i, i - 1] = 1if i < self.N - 2:self.A[i, i + 1] = 1self.b[i] = 3 * (X[i + 2] - X[i])self.A.factorizeLU()D = self.A.solve(self.b)D = np.vstack([np.zeros(2), D, np.zeros(2)])self.coefficients = []for i in range(self.N):a = X[i]b = D[i]c = 3 * (X[i + 1] - X[i]) - 2 * D[i] - D[i + 1]d = 2 * (X[i] - X[i + 1]) + D[i] + D[i + 1]self.coefficients.append((a, b, c, d))def getCurve(self, t):"""获取插值曲线的值"""idx = min(int(t * self.N), self.N - 1)a, b, c, d = self.coefficients[idx]t = t * self.N - idxreturn a + b * t + c * t**2 + d * t**3# 示例调用
spline = CubicSpline()
spline.setConditions([0, 0], [10, 10], 4)
spline.setInnerPoints(np.array([[2, 3], [4, 5], [7, 8]]))
t_values = np.linspace(0, 1, 100)
curve = np.array([spline.getCurve(t) for t in t_values])import matplotlib.pyplot as pltplt.plot(curve[:, 0], curve[:, 1], label="Cubic Spline")
plt.scatter([0, 2, 4, 7, 10], [0, 3, 5, 8, 10], color='red', label="Points")
plt.legend()
plt.show()
import numpy as np
import matplotlib.pyplot as plt# 定义三角函数曲线
x_values = np.linspace(0, 2 * np.pi, 10)  # 5个点，包含起点和终点
y_values = np.sin(x_values)  # 三角函数值# 初始化样条插值
spline = CubicSpline()
spline.setConditions([x_values[0], y_values[0]], [x_values[-1], y_values[-1]], len(x_values) - 1)
spline.setInnerPoints(np.array([x_values[1:-1], y_values[1:-1]]).T)# 生成插值曲线
t_values = np.linspace(0, 1, 100)  # 参数化曲线的 t 值
curve = np.array([spline.getCurve(t) for t in t_values])# 绘制结果
plt.figure(figsize=(8, 6))
plt.plot(curve[:, 0], curve[:, 1], label="Cubic Spline (Fitted)", linewidth=2)
plt.plot(x_values, y_values, 'r--', label="Original Sine Wave", linewidth=1.5)
plt.scatter(x_values, y_values, color='red', label="Data Points", zorder=5)
plt.legend()
plt.title("Fitting a Sine Wave Using Cubic Spline")
plt.xlabel("x")
plt.ylabel("y")
plt.grid(True)
plt.show()