样条函数

目录

• 定义
  • 三次样条
  • 三次样条函数代码
  • 自然三次样条最优性定理
• 高次自然样条理论
• B样条
  • B 样条的性质
  • B样条的导数和积分
  • 附加性质

《数值分析》 David Kincaid, Ward Cheney 著, 王国荣 余耀明 徐兆亮 译.

因为看的论文里用到了样条函数，故查阅了一些资料并整理一下.

定义

样条函数是由一些具有某些连续性条件得子空间上得分段多项式构成. 给定n+1个点(t_0, t_1, ldots, t_n)并且满足(t_0 < t_1 < ldots t_n.) 这些点称为结点. 又假如指定一个整数(k ge0), 具有结点(t_0, t_1, ldots, t_n)的一个(k)次样条函数是指满足下列条件的函数(S):

1. 在每一个区间([t_{i-1}, t_i])上, (S)是一个次数(le k)的多项式.
2. 在([t_0, t_n])上(S)有((k-1))阶连续导数.

三次样条

假设有数据((t_0, y_0), (t_1, y_1), ldots, (t_n, y_n)). 我们来探讨如何唯一确定一个三次样条函数.

记([t_i, t_{i+1}])上表示(S)的多项式记为(S_i).
首先根据连续性条件有:

[S_{i-1}(t_i)=y_i=S_{i}(t_i) quad (1 le i le n-1), \ S_0(t_0)=y_0, S_n(t_n)=y_n. ]

这是(2n)个方程,
再利用2阶连续可导的条件:

[ S'_{i-1}(t_i)=S'_i(t_i) quad (1 le i le n-1) \ S''_{i-1}(t_i)=S''_i(t_i) quad (1 le i le n-1) \ ]

可知，这里有(2n-2)个方程，总共有(4n-2)个方程，但是唯一确定一个三次样条函数，需要(4n)个方程，所以还有2个自由度，这个怎么利用后面会提到.

定义(z_i=S''(t_i)), 因为容易证明(S_i''(t_i)=z_i, S_i''(t_{i+1})=z_{i+1}), 且因为(S_i)是三次多项式，所以其二阶导数为一阶多项式，所以:

[S_i''(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).

进行俩次积分，其结果是(S_i):

[S_i(x) = frac{z_i}{6h_i} (t_{i+1} - x)^3 + frac{z_{i+1}}{6 h_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})可得:

[egin{array}{ll} S_i(x) &= frac{z_i}{6h_i} (t_{i+1} - x)^3 + frac{z_{i+1}}{6 h_i}(x-t_i)^3 \ &+ (frac{y_{i+1}}{h_i}-frac{z_{i+1}h_i}{6})(x-t_i) + (frac{y_i}{h_i} - frac{z_ih_i}{6})(t_{i+1}-x). end{array} ]

注: 原书在此处有误.

求在结点处的一阶导数:

[S_i'(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}'(t_i) = frac{h_{i-1}}{6} z_{i-1} + frac{h_{i-1}}{3}z_i - frac{y_{i-1}}{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,ldots, n-1. ]

这个线性方程组有未知量(z_0, z_1, ldots, z_n), 方程个数为(n-1)个，所以需要另外增添条件，一个很好的选择是(z_0=z_n=0)(后面会有一个定理为其提供依据).

[left [ egin{array}{cccccc} u_1 & h_1 & & & & \ h_1 & u_2 & h_2 & & & \ & h_2 & u_3 & h_3 & & \ & & ddots & ddots & ddots & \ & & & h_{n-3} & u_{n-2} & h_{n-2} \ & & & & h_{n-2} & u_{n-1} end{array} ight] left [ egin{array}{c} z_1 \ z_2 \ z_3 \ vdots \ z_{n-2} \ z_{n-1} end{array} ight ] = left [ egin{array}{c} v_1 \ v_2 \ v_3 \ vdots \ v_{n-2} \ v_{n-1} end{array} ight ], ]

其中

[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 = b_i - b_{i-1}. ]

三次样条函数代码

"""
三次样条函数实现
"""

import numpy as np
import matplotlib.pyplot as plt

class Polynomial:

    def __init__(self, t1, t2, y1, y2, z1, z2):
        assert t2 > t1, "t2 > t1 needed..."
        self.z = (z1, z2)
        self.y = (y1, y2)
        self.t = (t1, t2)
        self.h = t2 - t1

    def __call__(self, x):
        if x < self.t[0] or x > self.t[1]:
            return 0.
        reduce1_t2_x = self.t[1] - x
        reduce1_x_t1 = x - self.t[0]

        return self.z[0] / (6 * self.h) * reduce1_t2_x ** 3 + 
                self.z[1] / (6 * self.h) * reduce1_x_t1 ** 3 + 
               (self.y[1] / self.h - self.z[1] * self.h / 6) * reduce1_x_t1 + 
               (self.y[0] / self.h - self.z[0] * self.h / 6) * reduce1_t2_x


class Spline3:

    def __init__(self, t, y):
        self.n = len(t) - 1
        assert len(y) == self.n + 1, "t and y not match"
        self.t = np.array(t, dtype=float)
        self.y = np.array(y, dtype=float)
        self.h = np.diff(t)
        assert np.all(self.h > 0), "Invalid t"
        self.b = np.diff(y) * 6 / self.h
        self.calc_z()  #计算z
        self.generate_splines()  #生成分段多项式函数

    def calc_z(self):
        u = []
        v = []
        z = [0.]
        u.append((self.h[0] + self.h[1]) / 2)
        v.append(self.b[1] - self.b[0])
        for i in range(1, self.n-1): #对角化
            u_i = 2 * (self.h[i] + self.h[i-1]) - 
                  self.h[i-1] ** 2 / u[-1]
            v_i = self.b[i] - self.b[i-1] - 
                  self.h[i-1] * v[i-1] / u[-1]
            u.append(u_i)
            v.append(v_i)
        z.append(v[-1] / u[-1])
        for i in reversed(range(self.n-2)): #计算解
            z.append((v[i] - self.h[i] * z[-1]) / u[i])
        z.append(0.)
        z.reverse()
        self.z = z
        return z

    def generate_splines(self):
        s = []
        for i in range(self.n):
            t1 = self.t[i]
            t2 = self.t[i+1]
            y1 = self.y[i]
            y2 = self.y[i+1]
            z1 = self.z[i]
            z2 = self.z[i+1]
            s.append(Polynomial(t1, t2, y1,
                                y2, z1, z2))
        self.s = s
        return s

    def __call__(self, x):
        def f(x):
            if x < self.t[0] or x > self.t[-1]:
                return 0.
            for item in self.s:
                result = item(x)
                if result:
                    return result
            raise ValueError("Something wrong happened...")

        if not hasattr(x, "__len__"):
            return f(x)
        else:
            y = []
            for item in x:
                y.append(f(item))
            return np.array(y)


    def plot(self):
        t = np.linspace(self.t[0], self.t[-1], 100)
        y = self(t)
        fig, ax = plt.subplots()
        ax.plot(t, y, "b--")
        ax.set(xlabel="x", ylabel="y")
        ax.scatter(self.t, self.y, color="red")
        plt.show()

自然三次样条最优性定理

这个“最优性”似乎用“最光滑”来理解更为恰当吧. 这个定理也解释了之前的为什么(z_0=z_n=0).

定理1（自然三次样条最优性定理) 设(f'')在([a, b])上连续且(a=t_0<t_1<cdots<t_n=b). 若(S)是(f)在节点(t_i)上的自然三次样条插值，(0le i le n), 则:

[int_a^b [S''(x)]^2 mathrm{d}x le int_a^b [f''(x)]^2 mathrm{d}x. ]

证明: 令(g equiv f-S), 从而对于(0 le i le n), (g(t_i)=0)并且

[int_a^b (f'')^2 mathrm{d}x = int_a^b (S'')^2 mathrm{d}x + int_a^b (g'')^2 mathrm{d}x + 2 int_a^b S'' g'' mathrm{d}x, ]

只需证明:

[int_a^b S''g'' mathrm{d}x ge 0. ]

注意到: (S''(t_0) = S''(t_n)=0), 以及在([t_{i-1}, t_i])上(S''')是常数(记为(c_i)), 我们有:

[egin{array}{ll} int_a^b S'' g'' mathrm{d}x &= sum_{i=1}^n int_{t_{i-1}}^{t_i} S'' g'' mathrm{d}x \ & = sum_{i=1}^n { (S''g')|_{t_{i-1}}^{t_{i}}-int_{t_{i-1}}^{t_{i}} S'''g'mathrm{d}x} \ & = -sum_{i=1}^n int_{t_{i-1}}^{t_i} S''' g' mathrm{d}x = - sum_{i=1}^n c_i int_{t_{i-1}}^{t_i} g' mathrm{d}x \ &= -sum_{i=1}^n c_i [g(t_i) - g(t_{i-1})] = 0. end{array} ]

实际上，条件可以放松，注意到在上述证明步骤中:

[sum_{i=1}^n (S''g')|_{t_{i-1}}^{t_i} = (S''g')(b) - (S''g')(a). ]

当这个表达式非负的时候，结论依然成立，即需要

[S''(b)[f'(b)-S'(b)] ge S''(a)[f'(a) - S'(a)]. ]

例如: (S'(a)=f'(a))和(S'(b) = f'(b)).

高次自然样条理论

自然样条： 一个(2m+1)次的自然样条是一个函数(S in C^{2m}(R)), 在每一个区间([t_0, t_1], cdots, [t_{n-1},t_n])内，它都化为一个次数(le 2m+1)的多项式，而在区间((-infty,t_0))和((t_n, +infty))内化为一个次数至多为(m)的多项式.

降在(n+1)个结点(t_0, t_1, ldots, t_n)上的全体(2m+1)次自然样条所构成的线性空间记为(mathcal{N}^{2m+1}(t_0,t_1,cdots,t_n)), 简记为(mathcal{N}_n^{2m+1}).

下面不加证明地给出定理:

定理2（截断幂函数定理）： (mathcal{N}_n^{2m+1})中地每个元有表达式

[S(x) = sum_{i=1}^m a_ix^i + sum_{i=0}^n b_i (x-t_i)_+^{2m+1}, ]

其中对于(0 le j le m, sum_{i=0}^n b_it_I^j=0.) ((x-t_i)_+=x-t_i, : x-t_ige0,) 否则为(0).

定理3（奇数次自然样条唯一性定理）： 给定结点(t_0 < t_1 < cdots < t_n,) 设(0 le m le n). 则存在唯一的(2m+1)次自然样条在这些结点上渠道给定的值.

证明很有趣，很在意奇数次的限制，这里给出证明.

证明:

[S(x) = sum_{i=1}^m a_ix^i + sum_{i=0}^n b_i (x-t_i)_+^{2m+1}, ]

假设(S(t_i)=lambda_i)，再结合定理2，只需下列方程组:

[left { egin{array}{ll} S(t_i)= lambda_i & 0 le i le n \ sum_{j=0}^n b_jt_j^i=0 & 0 le i le m end{array} ight. , ]

存在唯一解，这一个(m+n+2)个未知量(m+n+2)个方程的方程组，所以只需证明其对应的其次问题仅有零解即可. 假设(lambda_i=0, 0 le i le n). 下证:

[I equiv int_a^b [S^{(m+1)}(x)]^2 mathrm{d}x=0, ]

其中(a=t_0,b=t_n). 利用分部积分可以得到:

[I = S^{(m+1)}(x) S^{(m)}(x)|_a^b - int_a^b S^{(m)}(x)S^{(m+2)}(x)mathrm{d}x=-int_a^b S^{(m)}(x)S^{(m+2)}(x)mathrm{d}x, ]

其中最后一个等式成立的原因是(S^{(k)}(a)=S^{(k)}(b)=0, k>m). 重复上述操作可得:

[I = (-1)^m int_a^b S^{(1)}(x)S^{(2m+1)}(x) mathrm{d}x. ]

因为(S^{(2m+1)}(x))是分段常值函数，所以:

[I = (-1)^m sum_{i=1}^nint_{t_{i-1}}^{t_i} c_i S'(x) mathrm{d}x = (-1)^m sum_{i=1}^nc_i[S(t_i)-S(t_{i-1})]=0. ]

由此，我们可知(S^{(m+1)} equiv 0). 因此，(S)是一个次数至多是(m)次的多项式，但(S)有(n+1>m)个不同的零点，所以(S=0), 所以(a_i, b_i=0). 至此，唯一性得证.

考虑偶数次样条，我不知道是怎么定义的，也不想去查，假设是(2m, m)的类型，那么我们依然要证明(或者(m+2,m+3,ldots)):

[I equiv int_a^b [S^{(m+1)}(x)]^2 mathrm{d}x=0, ]

不然利用分部积分的时候，没法消项, 能顺利推导到:

[I = (-1)^{m-1} int_a^b S^{(2)}(x)S^{(2m)}(x) mathrm{d}x. ]

此时(S^{(2m)})已然是分段常值函数，则:

[I = (-1)^{m-1} sum_{i=1}^nint_{t_{i-1}}^{t_i} c_i S''(x) mathrm{d}x = (-1)^{m-1} sum_{i=1}^nc_i[S'(t_i)-S'(t_{i-1})]. ]

所以没办法证明其为0. (m+2, m+3, ldots)更不必证明，因为没法保证(n+1>m+1)等.

如果是(2m,m-1)类型的，则需要证明:

[I equiv int_a^b [S^{(m)}(x)]^2 mathrm{d}x=0, ]

可以顺利推导到:

[I = (-1)^{m-1} int_a^b S^{(1)}(x)S^{(2m-1)}(x) mathrm{d}x. ]

显然也没法往下推.

定理4（奇数次自然样条最优性定理）： 设(mle n)及(f in C^{m+1}[a,b].) 假设(S)是在结点(a=t_0<t_1<cdots<t_n=b)上插值(f)的(2m+1)次自然样条，则

[int_a^b [S^{(m+1)}(x)]^2 mathrm{d}x le int_a^b [f^{(m+1)}(x)]^2 mathrm{d}x. ]

B样条

这一节专门讨论样条函数系统，使得其他所有样条函数都可以由它的线性组合得到. 这些样条构成某些样条空间的基. 所以称为(B)样条. 一旦给定结点，(B)样条就很容易通过递归关系产生. 而且算法也比较简单. (B) 样条以其优美的理论和数值计算中的典型性质著称. 此外，(B)杨i套还可以得到进一步的推广.

我们在无限集上讨论:

[cdots < t_{-2} < t_{-1}<t_0 < t_1 <t_2 < cdots ]

但是，在后面，我们会发现，在实际使用中，我们都只会用到有限个结点，在无限集上只是便于讨论.

首先给出0次(B)样条的定义:

[B^0_i (x) = left { egin{array}{ll} 1 & x in [t_i, t_{i+1}) \ 0 & else end{array} ight. . ]

容易得到(sum_{i=-infty}^{infty} B_i^0(x)=1).

高次(B)样条由递归定义:

[B_i^k(x) = (frac{x-t_i}{t_{i+k}-t_i})B_i^{k-1}(x)+(frac{t_{i+k+1}-x}{t_{i+k+1}-t_{i+1}})B^{k-1}_{i+1}(x) quad (kge 1). ]

若令

[V_i^k (x) = frac{x-t_i}{t_{i+k}-t_i}, ]

则

[B_i^k = V_i^k B_i^{k-1} + (1-V_{i+1}^k)B^{k-1}_{i+1}. ]

B 样条的性质

引理1（B样条的支撑引理）： 若(kge1)并且(x ot in (t_i, t_{i+k+1})), 则(B_i^k(x)=0).

引理2（B样条正性引理）： 设(kge 0). 若(x in (t_i, t_{i+k+1})), 则(B_i^k(x)> 0).

引理3（B样条的递归关系引理）： 对一切(kge 0)， 我们有

[sum_{i=-infty}^{infty} c_iB_i^k = sum_{i=-infty}^{infty} [c_iV_i^k+c_{i-1}(1-V_i^k)]B_i^{k-1}. ]

引理4（B样条单位分解引理）：对一切(k ge 0), 我们有

[sum_{i=-infty}^{infty} B_i^k(x) = 1. ]

B样条的导数和积分

记(alpha_i^k = frac{1}{t_{i+k}-t_i}), 有

[frac{mathrm{d}}{mathrm{d}x} V_i^k(x) = alpha_i^k \ alpha_i^k V_i^{k+1} = alpha_i^{k+1}V_i^k \ alpha_{i+1}^k (1-V_i^{k+1})=alpha_i^{k+1} (1-V_{i+1}^k). ]

引理5（B样条导数引理）： 对于(kge 2), B样条函数的导数可如下计算:

[frac{mathrm{d}}{mathrm{d}x} B_i^k(x) = (frac{k}{t_{i+1}-t_i}) B_i^{k-1}(x) - (frac{k}{t_{i+k+1}-t_{i+1}})B_{i+1}^{k-1}(x). ]

引理6（B样条光滑性引理）： 对于(k ge 1), B 样条(B_i^k)属于连续函数类(C^{k-1}(R)).

从与B样条导数相关的引理，可得：

[frac{mathrm{d}}{mathrm{d}x} sum_{i=-infty}^{infty} c_i B_i^k(x) = k sum_{i=-infty}^{infty} (frac{c_i-c_{i-1}}{t_{i+1}-t_i})B_i^{k-1}(x) quad (k ge 2). ]

引理7（B样条积分引理） B样条函数的积分可如下计算：

[int_{-infty}^x B_i^k(s) mathrm{d}s = (frac{t_{i+k+1}-t_i}{k+1})sum_{j=i}^{infty}Bj^{k+1}(x). ]

附加性质

引理8（B样条线性无关性引理）： B 样条集合({B_j^k, B_{j+1}^k, cdots, B_{j+k}^k})在((t_{k+j},t_{k+j+1}))上线性无关.

引理9（B样条线性无关性引理）： B样条集合({B_{-k}^k, B_{-k+1}^k, cdots, B_{n-1}^k})在((t_0, t_n))上线性无关.

记由在区间([t_0, t_1], [t_1, t_2], cdots, [t_{n-1}, t_n])上的(k)次样条函数所构成线性空间为(mathcal{S}_n^k).

定理1（空间(mathcal{S}_n^k)的基定理）： 空间(mathcal{S}_n^k)的一个基是

[{ B_i^k| [t_0, t_n]:-kle i le n-1} ]

从而，(mathcal{S}_n^k)的维数是(k+n).

插值矩阵A的定义：

[A_{ij} = B_j^k(x_i) quad (i le i,j le n). ]

引理1 （插值矩阵引理1）： 若插值矩阵非奇异，则(A_{ii} ot = 0, : 1 le i le n).

引理2（插值矩阵引理2）： 若(k=1)并且对于(1 le i le n)有(t_i < x_i < t_{i+2}), 则(A)非奇异.

定理2 （Schoenberg-Whitney定理）： 设(x_1 < x_2 cdots < x_n). 已知(A_{ij}=B_j^k(x_i))，则矩阵(A)非奇异当且仅当其主对角线上不含零元.

定理3（B样条插值定理）： 若([t_0,t_n])中的结点(x_1 , x_2, cdots, x_{n+k})满足

[t_{i-k-1}< x_i < t_i quad (1 le i le n+k), ]

则能用样条空间(mathcal{S}_n^k)插值这组节点上的任意一组数据.

注意到，如果我们的插值结点恰好就是(t_0, cdots, t_n), 那么仅有(n+1)个方程，还有(k-1)个自由度，此时方程的解是不唯一的。