9.3 分段立方插值

主观认识

  分段立方是一种多项式插值，是将离散的原始数据还原为一个分段的立方函数（三次函数），英文为Cubic Spline Interpolation。其算法叫立方曲线插值算法Cubic spline interpolation algorithm。其基本思想是生成分段的三次多项式，在连接点处可导，也就是多个分段函数组成的整个函数是连续、可导的。
  我举个例子，以下是9个点组成的数据集，组成了9个小向量：
$$\begin{gathered} v_0=(1,1) \\ {v}_1=(2,2) \\ {v}_2=(3,1) \\ {v}_3=(4,-2) \\ {v}_4=(5,1) \\ {v}_5=(6,2) \\ {v}_6=(7,1) \\ {v}_7=(8,-2) \\ {v}_8=(9,1) \end{gathered}$$
  在坐标系上是这个样子：

  这个数据，凭直觉看，是一个周期函数，但是计算机不知道啊。如果我们用牛顿插值会是什么样子呢？会得到这样一个函数（近似值，用下面这个方程在绘图软件误差特别大）：
$$\begin{gathered} f(x)=0.0032x^8-0.130x^7+2.2444x^6 \\ -21.11x^5+117.378x^4-390.378x^3 \\ +747.3746x^2-740.381x+286.0 \end{gathered}$$
  算了，我给出精确方程吧：
$$\begin{gathered} f(x)=0.0031746031746031746x^8 \\ -0.13015873015873017x^7 \\ +2.244444444444444x^6 \\ -21.11111111111111x^5 \\ +117.3777777777778x^4 \\ -390.37777777777785x^3 \\ +747.3746031746031x^2 \\ -740.3809523809525x+286.0 \end{gathered}$$
  在绘图软件(Gnuplot)中得到的图形是这样的：

  可见偏离特别大，比如1.5，得到的是$-3.1679687500005684$，而我们预期的是一个在$(1,2)$区间内的数字。
  我们再试试分段立方插值。分段立方插值得到的是$1.6005154639175259$，与真实值比较接近。分段立方插值的图像是什么样子的呢？就是下图这个样子：

  与牛顿/拉格朗日插值对比可以更明显看出差异：
  主观上看，分段立方插值更接近原始数据，更精确。

样本方程

  分段立方插值获取的是分段立方（三次）函数，假设样本有$n$对数据，那么函数要分成$n-1$，这些样本点就是边界，首先要保证分段函数要经过样本点。
  为了简便，我们定义$h_i=x_{i+1}-x_i$。
  分段函数的形式是$y=a(x-x_i{)}^3+b(x-x_i{)}^2+c(x-x_i)+d$，定义域为$[x_i,x_{i+1}]$。要穿过样本点，那么意味着当$x=x_i$时，$y=y_i$，那么我们代进去算一下：
$$\begin{gathered} y_i=f(x_i) \\ {y}_i=a_i(x_i-x_i{)}^3+b_i(x_i-x_i{)}^2+c_i(x_i-x_i)+d_i \\ {y}_i=d_i \end{gathered}$$
  所以第一个未知参数$d$确定了，就等于样本的值$y_i$。那么还剩下三个未知数$a,b,c$。我们主观上看，插值函数形成的曲线是丝滑流畅的，那么在数学上，就是连续、可导的。因为每段函数经过两个点，左边界点是$(x_i,y_i)$，右边界点是$(x_{i+1},y_{i+1})$。
  知道这个条件后，我们得到了第二个方程：
$$\begin{gathered} y_{i+1}=a_i(x_{i+1}-x_i{)}^3+b_i(x_{i+1}-x_i{)}^2+c_i(x_{i+1}-x_i)+d_i \\ {y}_{i+1}=a_i(x_{i+1}-x_i{)}^3+b_i(x_{i+1}-x_i{)}^2+c_i(x_{i+1}-x_i)+y_i \\ {y}_{i+1}=a_i{h}_i^3+b_i{h}_i^2+c_i{h}_i+y_i \end{gathered}$$
  这只是连续，在连接处还需要可导。分段函数的形式是一个复合函数所以我们要用链式法则求导：
$$\begin{gathered} u(x)=x-x_i \\ f(u)=a{u}^3+b{u}^2+cu+d \\ {u}^{\prime }(x)=1 \\ {f}^{\prime }(x)=f^{\prime }(u)\cdot u^{\prime }(x)=3a(x-x_i{)}^2+2b(x-x_i)+c \end{gathered}$$
  对于中间的分段函数，左连接点要与前一段在连接处的导数相等，也要与后一段在连接处导数相等。假设前一段的分段函数已经解出来了，成为已知数，定义为g(x)，那么我们就得到了第三个方程，起点处的导数等于前一段函数在终点的导数：
$$\begin{gathered} f_i^{\prime }(x_i)=f_{i-1}^{\prime }(x_i) \\ 3a_i(x_i-x_i{)}^2+2b_i(x_i-x_i)+c_i=3a_{i-1}(x_i-x_{i-1}{)}^2+2b_{i-1}(x_i-x_{i-1})+c_{i-1} \\ {c}_i=3a_{i-1}(x_i-x_{i-1}{)}^2+2b_{i-1}(x_i-x_{i-1})+c_{i-1} \\ {c}_i=3a_{i-1}{h}_{i-1}^2+2b_{i-1}{h}_{i-1}+c_{i-1} \end{gathered}$$
  此外，在每个分段右边界，也是连续可导的，也就是后一段函数的一阶导相等，又可以得到一个方程：
$$\begin{gathered} f_i^{\prime }(x_{i+1})=f_{i+1}^{\prime }(x_{i+1}) \\ 3a_i(x_{i+1}-x_i{)}^2+2b_i(x_{i+1}-x_i)+c_i=3a_{i+1}(x_{i+1}-x_{i+1}{)}^2+2b_{i-1}(x_{i+1}-x_{i+1})+c_{i+1} \\ 3a_i{h}_i^2+2b_i{h}_i+c_i=c_{i+1} \end{gathered}$$
  这样，$c$这个未知数也出来了。分段立方插值还有个设定，分段函数的一阶导必须连续可导，丝滑流畅。其实就是二阶导要相等。为了直观展示，我把上面例子的一阶导用绿色线条画出来：

  设前一段函数的二阶导为$f_{i-1}^{\prime \prime }(x)$，因为二阶导的一般方程是$f^{\prime \prime }(x)=6a(x-x_i)+2b$，连接点处二阶导相等，我们又得到一个方程：
$$\begin{gathered} f_i^{\prime \prime }(x)=6a_i(x-x_i)+2b_i \\ {f}_i^{\prime \prime }(x_i)=f_{i-1}^{\prime \prime }(x_i) \\ 2b_i=6a_{i-1}(x_i-x_{i-1})+2b_{i-1} \\ {b}_i=3a_{i-1}{h}_{i-1}+b_{i-1} \end{gathered}$$
  也等同于：
$$a_{i-1}=\frac{b_i-b_{i-1}}{3h_{i-1}}$$
  与后一段在连接处二阶导相等，我们有：
$$\begin{gathered} f_i^{\prime \prime }(x_{i+1})=f_{i+1}^{\prime \prime }(x_{i+1}) \\ 6a_i(x_{i+1}-x_i)+2b_i=6a_{i+1}(x_{i+1}-x_{i+1})+2b_{i+1} \\ 6a_i{h}_i+2b_i=2b_{i+1} \\ 3a_i{h}_i+b_i=b_{i+1} \end{gathered}$$
 也等同于：
$$a_i=\frac{b_{i+1}-b_i}{3h_i}$$
  总结下就是四个方程：
$$\begin{gathered} d_i=y_i \\ {h}_i^3{a}_i+h_i^2{b}_i+h_i{c}_i=y_{i+1}-y_i \\ 3h_i^2{a}_i+2h_i{b}_i+c_i=c_{i+1} \\ 3h_i{a}_i+b_i=b_{i+1} \end{gathered}$$
  除去已经有解的$d_i$，还剩三个方程，却有5个未知数，所以还差一些方程。

首末分段

  这么依赖前一段函数，也就是要求出前段函数的$a,b,c$三个未知数，但是对于第一段函数，是没有前一段函数的，怎么办呢？
  这就需要讲讲分段立方插值的两个重要设定了，起点二阶导为0与终点二阶导为0。为什么要这两个设定呢？我们可以用物理学的概念来理解，二阶导相当于加速度，而加速度源于力，这其实假设起点和终点是受力平衡的状态。在有些文章里是写错了的，我举个例子，国外著名的网站，就把首末端二阶导为0写成了首末端一阶导为0，这就不够严谨了，哈哈：

  但是对于其他的英文网站，我推荐这篇加州伯克利写的文章。
  对于第一段函数，二阶导方程如下：
$$\begin{gathered} f^{\prime \prime }(x_0)=0 \\ 6a_0(x_0-x_0)+2b_0=0 \\ 2b_0=0 \end{gathered}$$
  所以第一段函数$b_0=0$，这就减少一个了一个未知数。对于终点的二阶导为0，假设终点为$(x_n,y_n)$我们计算一下：
$$\begin{gathered} f^{\prime \prime }(x_n)=6a_{n-1}(x_n-x_{n-1})+2b_{n-1}=0 \\ 3a_{n-1}(x_n-x_{n-1})+b_{n-1}=0 \end{gathered}$$
  哈哈，但是知道了这个方程，怎么与第一分段二阶导联系起来呢?
  因为n个数据点，那么有$n-1$个分段函数，除去$b_0=0$总共有$3n-4$个未知数。每一个段都有这个方程：
$$h_i^3{a}_i+h_i^2{b}_i+h_i{c}_i=y_{i+1}-y_i$$
  那么就现在就有了$n-1$个方程。每一段和后一段的关系，有以下两个方程：
$$\begin{gathered} 3h_i^2{a}_i+2h_i{b}_i+c_i=c_{i+1} \\ 3h_i{a}_i+b_i=b_{i+1} \end{gathered}$$
  那么就有$2n-4$个方程，加起来有$3n-5$个方程。再加上$3a_{n-1}(x_n-x_{n-1})+b_{n-1}=0$，这个方程，总共有$3n-4$个方程。未知数与方程的数量恰好相等，然后就可以解方程求出所有系数了。

泰勒级数

  大学时微积分最后一课就是泰勒级数和麦克劳林级数。大部分时候我们觉得这除了用来计算sin(x)和cos(x)的近似值外没什么用，而且真正到了计算机数值计算领域，已经是用切比雪夫近似值来计算了，不用泰勒级数了。原因很简单，泰勒级数一不精确，二性能低下。
  我们对任意一段函数进行泰勒级数展开：
$$f_i(x)=f_i(x_i)+f_i^{\prime }(x_i)(x-x_i)+\frac{f_i^{\prime \prime }(x_i)}{2!}(x-x_i{)}^2+\frac{f_i^{\prime \prime \prime }(x_i)}{3!}(x-x_i{)}^3$$
  再看看我们的公式：
$$f_i(x)=a_i(x-x_i{)}^3+b_i(x-x_i{)}^2+c_i(x-x_i)+d_i$$
  对比一下，换个顺序，就会发现：
$$\begin{gathered} a_i=\frac{f_i^{\prime \prime \prime }(x_i)}{3!} \\ {b}_i=\frac{f_i^{\prime \prime }(x_i)}{2!} \\ {c}_i=f_i^{\prime }(x_i) \\ {d}_i=f_i(x_i) \end{gathered}$$
  哈哈，数学是不是特别有意思？

$\mu$数列与z数列

  直接这样解方程，程序性能比较低，在实际生产中，我们使用 $\mu$数列与z数列法。我们换个思路，找到一种更简便的方法，而不是强解$3n-4$元一次方程组。
  $\mu$数列的定义是：
$${\mu }_i=\{\begin{array}{r}\frac{x_{i+1}-x_i}{2(x_{i+1}-x_{i-1})-{\mu }_{i-1}(x_i-x_{i-1})},i\in {\mathbb{Z}}^{+}\\ 0,i=0\end{array}$$
  $z$数列的定义是：
$$z_i=\{\begin{array}{r}\frac{\frac{3[y_{i+1}(x_i-x_{i-1})-y_i(x_{i+1}-x_{i-1})+y_{i-1}(x_{i+1}-x_i)]}{(x_i-x_{i-1})(x_{i+1}-x_i)}-(x_i-x_{i-1})z_{i-1}}{2(x_{i+1}-x_{i-1})-{\mu }_{i-1}(x_i-x_{i-1})},i\in {\mathbb{Z}}^{+}\\ 0,i=0\end{array}$$
  因为公式过于复杂，所以使用我们之前定义的$h_i=x_{i+1}-x_i$，将公式改写一下：
$$\begin{gathered} {\mu }_i=\{\begin{array}{r}\frac{h_i}{2(x_{i+1}-x_{i-1})-{\mu }_{i-1}{h}_{i-1}},i\in {\mathbb{Z}}^{+}\\ 0,i=0\end{array} \\ {z}_i=\{\begin{array}{r}\frac{\frac{3[y_{i+1}{h}_{i-1}-y_i(x_{i+1}-x_{i-1})+y_{i-1}{h}_i]}{h_{i-1}{h}_i}-h_{i-1}{z}_{i-1}}{2(x_{i+1}-x_{i-1})-{\mu }_{i-1}{h}_{i-1}},i\in {\mathbb{Z}}^{+}\\ 0,i=0\end{array} \end{gathered}$$
  搞这么复杂有什么用？哈哈，因为：
$$b_{n-1}=z_{n-1}$$
  并且有：
$$b_i=z_i-{\mu }_i{b}_{i+1},i<n-1$$
  得出了系数b，其他的a和c两个系数就可以很快求出来了。具体推导过程我就不去证明了，使用这两个数列，计算过程比直接解方程快多了。$\mu$数列与z数列的公式特别复杂，我强烈建议收藏本文！

Python实现($\mu$z数列法)

  首先简单实现一个霍纳多项式

# _*_ coding:utf-8 _*_

class HornerPolynomial:

    def __init__(self, items: list):
        self.__polynomial = items

    @property
    def polynomial(self):
        return self.__polynomial

    # 本身是个函数
    def __call__(self, *args, **kwargs):
        result = 0
        for coefficient in self.polynomial:
            result = result * args[0] + coefficient
        return result

    def __str__(self):
        result = ''
        first = True
        n= len(self.polynomial)
        for i, coefficient in enumerate(self.polynomial):
            if coefficient == 0:
                continue
            if first:
                first = False
            else:
                result += '+'
            if coefficient != 1:
                result += str(coefficient) + '*'
            order = n - i - 1
            if order > 0:
                result += 'x'
            else:
                result += '1'
            if order > 1:
                result += '**' + str(order)
        return result

  再实现立方曲线插值：

# _*_ coding:utf-8 _*_
from com.youngthing.mathalgorithm.interpolation.horner_polynominal import HornerPolynomial


class PeacewisePolynomials:

    def __init__(self, konts, polynomials):
        self.__konts = konts
        self.__polynomials = polynomials

    def __call__(self, *args, **kwargs):
        x = args[0]
        for i, e in enumerate(self.__polynomials):
            if x >= self.__konts[i] and x <= self.__konts[i + 1]:
                return e(x - self.__konts[i])
        return None

    def __str__(self):
        x = ''
        for i, e in enumerate(self.__polynomials):
            s = '' + str(e) + ',[' + str(self.__konts[i]) + ',' + str(self.__konts[i + 1]) + ']\n'
            x += s
        return x


def getx(value):
    return value[0]


def interpolate(values):
    values = sorted(values, key=getx)
    n = len(values) - 1

    h = [0] * n
    for i in range(0, n):
        h[i] = getx(values[i + 1]) - getx(values[i])

    mu = [0] * n
    z = [0] * (n + 1)

    x = [getx(e) for e in values]
    y = [e[1] for e in values]
    for i in range(1, n):
        g = 2 * (x[i + 1] - x[i - 1]) - h[i - 1] * mu[i - 1]
        mu[i] = h[i] / g
        temp = 3 * (y[i + 1] * h[i - 1] - y[i] * (x[i + 1] - x[i - 1]) + y[i - 1] * h[i])
        z[i] = (temp / (h[i - 1] * h[i]) - h[i - 1] * z[i - 1]) / g

    c = [0] * n
    b = [0] * (n + 1)
    a = [0] * n

    for j in range(n - 1, -1, -1):
        b[j] = z[j] - mu[j] * b[j + 1];
        c[j] = (y[j + 1] - y[j]) / h[j] - h[j] * (b[j + 1] + 2 * b[j]) / 3
        a[j] = (b[j + 1] - b[j]) / (3 * h[j])

    polys = [HornerPolynomial([a[i], b[i], c[i], y[i]]) for i in range(0, n)]
    return PeacewisePolynomials(x, polys)

代码测试

from com.youngthing.mathalgorithm.interpolation.cubic_spline import *


if __name__ == '__main__':
    values=[
        [1,1],
        [2,2],
        [3,1],
        [4, -2],
        [5, 1],
        [6, 2],
        [7, 1],
        [8, -2],
        [9, 1],
    ]
    #测试数据
    cubic = interpolate(values)

    print(cubic(1.5))
    print(cubic)

  测试结果完全正确：

1.6005154639175259
-0.26804123711340205*x**3+1.268041237113402*x+1,[1,2]
-0.6597938144329897*x**3+-0.8041237113402062*x**2+0.463917525773196*x+2*1,[2,3]
2.9072164948453607*x**3+-2.783505154639175*x**2+-3.1237113402061856*x+1,[3,4]
-2.9690721649484533*x**3+5.9381443298969065*x**2+0.030927835051546726*x+-2*1,[4,5]
0.9690721649484534*x**3+-2.9690721649484533*x**2+3.0*x+1,[5,6]
-0.9072164948453607*x**3+-0.06185567010309301*x**2+-0.030927835051546282*x+2*1,[6,7]
2.6597938144329896*x**3+-2.783505154639175*x**2+-2.8762886597938144*x+1,[7,8]
-1.731958762886598*x**3+5.195876288659794*x**2+-0.463917525773196*x+-2*1,[8,9]