樣條函數

目錄

• 定義
  • 三次樣條
  • 三次樣條函數代碼
  • 自然三次樣條最優性定理
• 高次自然樣條理論
• 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}\)可得:

\[\begin{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 [ \begin{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} \right] \left [ \begin{array}{c} z_1 \\ z_2 \\ z_3 \\ \vdots \\ z_{n-2} \\ z_{n-1} \end{array} \right ] = \left [ \begin{array}{c} v_1 \\ v_2 \\ v_3 \\ \vdots \\ v_{n-2} \\ v_{n-1} \end{array} \right ], \]

其中

\[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\)上的自然三次樣條插值，\(0\le 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\)), 我們有:

\[\begin{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_i\ge0,\) 否則為\(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 \{ \begin{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} \right. , \]

存在唯一解，這一個\(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}^n\int_{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}^n\int_{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（奇數次自然樣條最優性定理）： 設\(m\le 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 \{ \begin{array}{ll} 1 & x \in [t_i, t_{i+1}) \\ 0 & else \end{array} \right. . \]

容易得到\(\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 (k\ge 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樣條的支撐引理）： 若\(k\ge1\)並且\(x \not \in (t_i, t_{i+k+1})\), 則\(B_i^k(x)=0\).

引理2（B樣條正性引理）： 設\(k\ge 0\). 若\(x \in (t_i, t_{i+k+1})\), 則\(B_i^k(x)> 0\).

引理3（B樣條的遞歸關系引理）： 對一切\(k\ge 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樣條導數引理）： 對於\(k\ge 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]:-k\le 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} \not = 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\)個自由度，此時方程的解是不唯一的。