分段三次样条

一、插值和样条？

给定若干已知点位，连接、逼近这些点位的方式有很多种，我们既可以给定一个足够阶数的多项式一次性穿过、逼近这些点位，这些方法可以是：牛顿插值法、拉格朗日插值法；也可以按段给出每个区间的表达式，比如说三次样条。

上面展示的就是按段给出的表达式，按段一条条给出的，所以叫做样条。当然你也可以用一个多项式穿过这些点：

红色的线就是用一个多项式连接的方法。用一个多项式拟合所有的点位不都是好的，尤其是均匀分布的插值点两端因为它会产生龙格现象：

从上图可以看出，左右两端的拟合对于原函数完全不匹配。缓解方法有很多种，可以使用切比雪夫或者分段三次样条进行处理。也就是文章的重点，分段三次多项式样条。

二、分段三次样条（自然边界条件）

样条的表达式是通过多个分段函数给出的，约束条件有：

• 满足节点自变量递增关系；
• 相邻多项式之间一二阶导数连续；
• 多项式的次数是3次；

如果每相邻两个点之间的连线是三次多项式，那么对应的整段组合而成的曲线就叫做三次样条函数。假设点数是 $n+1$，那么段数为：$n$

• 三次多项式有四个未知数，一共有 $n$ 个，所以待求的方程组个数为 $4n$；
• 一共有 $n+1$ 个点，我已有 $(n+1)$ 个方程组，这样一来我就剩下 $4n-(n+1)=3n-1$；
• 一共有 $n-1$ 个衔接点，衔接点处函数值连续、一阶导连续，二阶导连续，提供 $3(n-1)$个方程，剩余$3n-1-3(n-1)=2$

根据不同的约束可以将三次样条分为三类：

• 端点处二阶导函数值约束，如果都是0，则为自然样条；
• 端点处一阶导数约束；
• 周期约束，略；

给定以下插值点：

$x$	$t_0$	$t_1$	$\cdots$	$t_n$
$y$	$y_0$	$y_1$	$\cdots$	$y_n$

根据这个表格结合三次样条多项式的定义：
$$S(x)=\{\begin{array}{rl}& S_0(x)x\in [t_0,t_1]\\ & S_1(x)x\in [t_1,t_2]\\ & ⋮\\ & S_{n-1}(x)x\in [t_{n-1},t_n]\end{array}$$

在此之前我们先来推导区间 $[t_i,t_{i+1}]$ 上 $S_i(x)$ 的表达式。首先我们定义一组数 $z_i=S^{\prime \prime }(t_i)$，因为 $S^{\prime \prime }$ 在每个内结点上连续，所以显然，对于 $0\le i\le n$ 存在 $z_i$ 满足：
$$\underset{x\to t_i^{+}}{\lim }{S}^{\prime \prime }(x)=\underset{x\to t_i^{-}}{\lim }{S}^{\prime \prime }(x)$$
因为 $S_i$ 是在 $[t_i,t_{i+1}]$ 上的三次多项式，因此是满足 $S_i^{\prime \prime }(t_i)=z_i$、 $S_i^{\prime \prime }(t_{i+1})=z_{i+1}$ 的线性函数，根据根据拉格朗日插值法有：
$$S_i^{\prime \prime }(x)=\frac{z_i}{h_i}(t_{i+1}-x)+\frac{z_{i+1}}{h_i}(x-t_i)$$

其中 $h_i=t_{i+1}-t_i$ 。为了利用插值点，我们对$(3)$ 对变量 $x$积分两次有：
$$S_i(x)=\frac{z_i}{6h_i}(t_{i+1}-x{)}^3+\frac{z_{i+1}}{6h_i}(x-t_i{)}^3+C(x-t_i)+D(t_{i+1}-x)$$

这里直接用整体法进行积分不要直接拆开会简单一些

其中 $C$ $D$ 是积分常数。

代入插值条件 $S_i(t_i)=y_i$ 和 $S_i(t_{i+1})=y_{i+1}$ 作用在 $S_i$ 就可以去确定 $C$ $D$ 两个常数，其结果为：

$$\begin{array}{r}{S}_i(x)=\frac{z_i}{6h_i}(t_{i+1}-x{)}^3+\frac{z_{i+1}}{6h_i}(x-t_i{)}^3+\left(\frac{y_{i+1}}{h_i}-\frac{z_{i+1}{h}_i}{6}\right)(x-t_i)+\left(\frac{y_i}{h_i}-\frac{z_i{h}_i}{6}\right)(t_{i+1}-x)\end{array}$$

此时整个方程只有二阶导数是未知的，为了利用一阶导数的连续性，我们对 $(5)$ 微分一次代入 $S_{i-1}^{\prime }(t_i)=S_i^{\prime }(t_i)$ 条件：
$$S_i^{\prime }(t_i)=-\frac{h_i}{3}{z}_i-\frac{h_i}{6}{z}_{i+1}-\frac{y_i}{h_i}+\frac{y_{i+1}}{h_i}$$
同理，可以求得 $S_{i-1}^{\prime }(t_i)$：
$$S_{i-1}^{\prime }(t_i)=\frac{h_{i-1}}{6}{z}_{i-1}+\frac{h_{i-1}}{3}{z}_i-\frac{y_i}{h_{i-1}}+\frac{y_i}{h_{i-1}}$$
因为一阶导数连续，所以：
$$h_{i-1}{z}_{i-1}+2(h_i+h_{i-1})z_i+h_i{z}_{i+1}=\frac{6}{h_i}(y_{i+1}-y_i)-\frac{6}{h_{i-1}}(y_i-y_{i-1})$$
上面的方程对于 $i=1,2,\cdots ,n-1$ 成立。注意这个式子没有对两个端点成立，这是以为他们前面和后面都没有对应的三次多项式，所以我们给出了一个含有 $n+1$ 个未知量 $z_0,z_1,\cdots ,z_n$ 的一个 $n-1$ 个方程的线性方程组。只要任选 $z_0$ $z_1$ 就可以计算出确定这个分段多项式方程。最常见的方法就是取 $z_0=z_n=0$，由此得到的样条函数称为自然三次样条。对于 $1\le i\le n-1$ 和 $z_0=z_n=0$ ：
$$\begin{gathered} h_0{z}_0+2(h_1+h_0)z_1+h_1{z}_2=\frac{6}{h_1}(y_2-y_1)-\frac{6}{h_0}(y_1-y_0) \\ {h}_1{z}_1+2(h_2+h_1)z_2+h_2{z}_3=\frac{6}{h_2}(y_3-y_2)-\frac{6}{h_1}(y_2-y_1) \\ {h}_2{z}_2+2(h_3+h_2)z_3+h_3{z}_4=\frac{6}{h_3}(y_4-y_3)-\frac{6}{h_2}(y_3-y_2) \\ \cdots \\ {h}_{n-2}{z}_{n-2}+2(h_{n-1}+h_{n-2})z_{n-1}+h_{n-1}{z}_n=\frac{6}{h_{n-1}}(y_n-y_{n-1})-\frac{6}{h_{n-2}}(y_{n-1}-y_{n-2}) \end{gathered}$$
令 $u_i=2(h_i+h_{i-1})$ $b_i=\frac{6}{h_i}(y_{i+1}-y_i)$ $v_i=b_i-b_{i-1}$
$$\begin{gathered} h_0{z}_0+u_1{z}_1+h_1{z}_2=v_1 \\ {h}_1{z}_1+u_2{z}_2+h_2{z}_3=v_2 \\ {h}_2{z}_2+u_3{z}_3+h_3{z}_4=v_3 \\ \cdots \\ {h}_{n-2}{z}_{n-2}+u_{n-1}{z}_{n-1}+h_{n-1}{z}_n=v_{n-1} \end{gathered}$$
写成矩阵形式

$$\left[\begin{array}{ccccccc}{u}_1& h_1& & & & & \\ h_1& u_2& h_2& & & & \\ & h_2& u_3& h_3& & & \\ ⋮& & & & & & \\ & & & h_{n-3}& u_{n-2}& h_{n-2}& \\ & & & & h_{n-2}& u_{n-1}& \end{array}\right]\left[\begin{array}{c}{z}_1\\ z_2\\ z_3\\ ⋮\\ z_{n-2}\\ z_{n-1}\end{array}\right]=\left[\begin{array}{c}{v}_1\\ v_2\\ v_3\\ ⋮\\ v_{n-2}\\ v_{n-1}\end{array}\right]$$
其中
$$\begin{array}{rl}{h}_i& =t_{i+1}-t_i\\ u_i& =2(h_i+h_{i-1})\\ b_i& =\frac{6}{h_i}(y_{i+1}-y_i)\\ v_i& =b_i-b_{i-1}\end{array}$$
$(11)$ 的系数矩阵是一个三对角矩阵，求解出对应的二阶导数值 $z_i$ 结合 $(5)$ 便可得到任意区间对应的插补函数值。

三、实现

为了求得每段的表达式，我们需要求解 $(11)$ 这个线性方程组，简单用python实现了以下并将其与标准库对比：

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import CubicSpline
from scipy.linalg import solve_banded

# 给定数据点
P = np.array(
    [[2, 3], [3, 5], [7, 11], [10, 20], [14, 10], [28, 34], [33, 23]],
)

x = P[:, 0]
y = P[:, 1]


def getUHV(x,y):
    h=np.zeros(len(x)-1)
    u=np.zeros(len(h)-1)
    b=np.zeros(len(y)-1)
    v=np.zeros(len(b)-1)
    
    h=np.diff(x)
    
    for i in range(len(u)):
        u[i]=2*(h[i+1]+h[i])
    
    for i in range(len(b)):
        b[i]=6*(y[i+1]-y[i])/h[i]
    
    v=np.diff(b)
    
    return u,h,v

def getZ(u,h,v):
    
    # | u_1     h_1     0       0      ...      0        0     | | z_1     |   = | v_1     |
    # | h_1     u_2     h_2     0      ...      0        0     | | z_2     |   = | v_2     |
    # | 0       h_2     u_3     h_3    ...      0        0     | | z_3     |   = | v_3     |
    # | .       .       .       .      ...      .        .     | | .       |   = | .       |
    # | .       .       .       .      ...      .        .     | | .       |   = | .       |
    # | .       .       .       .      ...      .        .     | | .       |   = | .       |
    # | 0       0       0       ...    h_{n-3} u_{n-2} h_{n-2} | | z_{n-2} |   = | v_{n-2} |
    # | 0       0       0       ...    0       h_{n-2} u_{n-1} | | z_{n-1} |   = | v_{n-1} |
    M=np.zeros((len(v),len(v)))
    np.fill_diagonal(M, u)
    np.fill_diagonal(M[1:],h[1:])
    np.fill_diagonal(M[:,1:], h[1:])
    z=np.linalg.solve(M,v)
    z=np.insert(z,0,0)
    z=np.append(z,0)
    
    return z
    
def cubicInterpolation(xx, z):
    n = len(x) - 1
    value = np.zeros_like(xx)
    count = 0

    for xi in xx:
        for k in range(n):
            if x[k] <= xi <= x[k + 1]: # chech which interval x falls into
                s1 = ((x[k + 1] - xi) ** 3) * (z[k] / (6 * h[k]))
                s2 = ((xi - x[k]) ** 3) * (z[k + 1] / (6 * h[k]))
                s3 = (y[k + 1] / h[k] - z[k + 1] * h[k] / 6) * (xi - x[k])
                s4 = (y[k] / h[k] - z[k] * h[k] / 6) * (x[k + 1] - xi)
                value[count] = s1 + s2 + s3 + s4
                count += 1
                break

    return value


# 自定义的三次样条插值
u, h, v = getUHV(x, y)
z = getZ(u, h, v)
xx = np.linspace(min(x), max(x), 50000)
yy_custom = cubicInterpolation(xx, z)

# 使用标准库进行三次样条插值
cs = CubicSpline(x, y, bc_type="natural")
yy_scipy = cs(xx)

# 绘图
plt.figure(figsize=(10, 6))
plt.plot(x, y, "o", label="Data points")
plt.plot(xx, yy_custom, "-", label="Custom cubic spline interpolation")
plt.plot(xx, yy_scipy, "--", label="Scipy cubic spline interpolation")
plt.legend()
plt.xlabel("x")
plt.ylabel("y")
plt.title("Comparison of Custom and Scipy Cubic Spline Interpolation")
plt.show()

在实际工程实现中，通常直接通过求解线性方程组的方式进行求解，利用三对角方程的性质递归求解系数比较合适。当然对于空间中的点需要进行一些参数化处理，对C++ 代码实现感兴趣，可以参考我上传的在github的代码：https://github.com/FISHFISHFISHLOVECAT/cublic_spline_curve

四、优缺点分析

分段三次多项式的缺点

• 三次只能保证段和段之间的一阶和二阶之间的连续性；
• 对于非常密集的点，为了准确通过这些点可能导致的过拟合（也就是龙格现象）；
• 是一种局部插补的方法，可能从整体上看不是我们想要的趋势；

可以想象，如果你有一条曲线，通过几个数据点在这条曲线附近取样，然后尝试用一段段的三次多项式将这些点连起来。如果这些点非常密集且分布不均匀，插值曲线可能会在不同点之间“剧烈地”上下摆动以适应所有点，而这种摆动可能并不代表数据的真实趋势，而只是为了准确地通过所有给定的点。这就是过拟合的核心问题——在追求高精度拟合时，失去了对整体趋势的把握。

---

[1]李庆扬 王能超.数值分析(第4版)[M].华中科技大学出版社(原华中理工出版社)出版社,2006.
[2]金凯德,WardCheney,金凯德,等.数值分析[M].机械工业出版社,2005.

原文地址：https://blog.csdn.net/weixin_39258979/article/details/143014207

免责声明：本站文章内容转载自网络资源，如本站内容侵犯了原著者的合法权益，可联系本站删除。更多内容请关注自学内容网（zxcms.com）！