三次样条插值

2023年11月5日
#analysis

---

---

1. 样条函数

样条函数即满足一定光滑性的分段多项式。
对区间 $(-\infty ,+\infty )$ 的一个分割：
$$\Delta :-\infty <x_1<x_2<\cdots <x_n<+\infty$$
若分段函数 $s(x)$ 满足条件：

1. 在每个区间 $(-\infty ,x_1],[x_j,x_{j+1}](j=1,\cdots ,n1)$ 和 $[x_n,+\infty )$ 上， $s(x)$ 是一个次数不超过 $m$ 的实系数代数多项式。
2. 光滑性要求：$s(x)$ 在整个区间 $(-\infty ,+\infty )$ 上具有直至 $m-1$ 阶的连续微商（导数）

则称 $y=s(x)$ 为对应于分割 $\Delta$ 的 $m$ 次样条函数， $x_1,x_2,\cdots ,x_n$ 为样条节点。
以 $x_1,x_2,\cdots ,x_n$ 为样条节点的 $m$ 次样条函数的全体记为：
$$s_m(x_1,x_2,\cdots ,x_n)$$
$$s(x)=\{\begin{array}{ll}{p}_0(x)& ,x\le x_1\\ p_1(x)& ,x_1\le x\le x_2\\ & ⋮\\ p_j(x)& ,x_j\le x\le x_{j+1}\\ & ⋮\\ p_n(x)& ,x_n\le x\end{array},p_j(x)\in P_m(j=0,1,\cdots ,n)$$
对样条函数而言，若令相邻两段函数的差：
$$q_j(x)=p_j(x)-p_{j-1}(x)\in P_m$$
$$\Rightarrow q_j^{(i)}(x)=p_j^{(i)}(x)-p_{j-1}^{(i)}(x)=0,i=0,1,\cdots ,m-1$$
$$\Rightarrow q_j(x)=c_j(x-x_j{)}^m,\text{也叫做光滑因子}$$
$x_j$ 是 $q_j(x)$ 的 $m$ 重根。所以如果函数是样条函数，则相邻两段函数满足：
$$p_j(x)=p_{j-1}(x)+c_j(x-x_j{)}^m,j=0,1,\cdots ,n$$
所以s(x)是m次样条的充要条件是
$$\begin{array}{rl}{p}_0(x)=& a_0+a_{1}x+\cdots +a_m{x}^m\\ & \\ p_n(x)=& p_0(x)+\sum _{j=1}^n{c}_j(x-x_j{)}^m\end{array}$$

1.1 截断多项式

为了方便表示分段信息，引进截断多项式：
$$(x-a{)}_{+}^m=\{\begin{array}{ll}(x-a{)}^m& ,x\ge a\\ & \\ 0& ,x<a\end{array}$$
显然 $(x-a{)}_{+}^m$ 是 $(-\infty ,+\infty )$ 上 $m-1$ 次连续可微函数的集合（ ${\mathbb{C}}^{m-1}(-\infty ,+\infty )$ 类的分段 $m$ 次多项式）。
使用阶段多项式表示样条函数有：任意 $s(x)\in S_m(x_1,x_2,\cdots ,x_n)$， $c_j\in \mathbb{R}$ 均可唯一地表示为：
$$s(x)=p_m(x)+\sum _{j=1}^n{c}_j(x-x_j{)}_{+}^m,-\infty <x<+\infty$$
样条函数的基底：
S m = span { 1 , x , ⋯ , x m , ( x − x 1 ) + m , ( x − x 2 ) + m , ⋯ , ( x − x n ) + m } S_m= \text{span} \lbrace 1,x,\cdots ,x^m,(x-x_1)_+^m,(x-x_2)_+^m,\cdots ,(x-x_n)_+^m \rbrace Sm​=span{1,x,⋯,xm,(x−x1​)+m​,(x−x2​)+m​,⋯,(x−xn​)+m​}
$$\text{dim}(S_m)=m+n+1$$

[!example]-
验证分片多项式是三次样条函数
$$s(x)=\{\begin{array}{ll}1-2x& ,x<-3\\ 28+25x+9x^2+x^3& ,-3\le x<-1\\ 26+19x+3x^2-x^3& ,-1\le x<0\\ 26+19x+3x^2& ,0\le x\end{array}$$
解：利用光滑因子验证。
$$(28x+25x+9x^2+x^3)-(1-2x)=(x+3{)}^3$$
$$(26+19x+3x^2-x^3)-(28+25x+9x^2+x^3)=-2(x+1{)}^3$$
$$(26+19x+3x^2)-(26+19x+3x^2-x^3)=x^3$$
该函数为三次样条函数。
光滑因子的零点是已知的，即为边界，光滑因子的常数项只需看等式左边最高次多项式的系数就能得到。之后把右边的多项式展开，看是不是和左边的相等就行。

---

2. 三次样条插值

设给定节点 $a=x_0<x_1<\cdots <x_n=b$ ，及节点上的函数值
$$f(x_i)=y_i,i=0,1,\cdots ,n$$
节点是中间的点！头尾点不是！三次样条问题就是构造 $s(x)\in S_3(x_1,x_2,\cdots ,x_{n-1})$ 满足插值条件
$$s(x_i)=y_i,i=0,1,\cdots ,n$$
且有一定光滑性。
边界条件分类

1. $s^{\prime \prime }(x_0)=y_0^{\prime \prime },s^{\prime \prime }(x_n)=y_n^{\prime \prime }$ ，当 $y_0^{\prime \prime }=y_n^{\prime \prime }=0$ ，为自然样条/自然边界
2. $s^{\prime }(x_0)=y_0^{\prime },s^{\prime }(x_n)=y_n^{\prime }$
3. $s^{\prime }(x_0^{+})=s^{\prime }(x_n^{-}),s^{\prime \prime }(x_0^{+})=s^{\prime \prime }(x_n^{-})$ 起始点和终止点导数相等，适用于周期函数

2.1 B样条为基底的三次样条插值函数

对 $[a,b]$ 进行 $n$ 等分 时候的情况
$$s(x)=\sum _{j=0}^{n+2}{c}_j{\Omega }_3(\frac{x-x_{j-1}}{h}),a\le x\le b,h=\frac{b-a}{n}$$
其中B样条函数
$$
\Omega_3 (y)= \begin{cases}
0 &, |y|\ge2 \\
\frac{1}{2}|y|^3-y2+ \frac{2}{3}&,|y|\le1 \ \

• \frac{1}{6} |y|^3+y2-2|y|+ \frac{4}{3} &, 1<|y|<2
  \end{cases}
  $$
  关键在于求 $c_j$

2.1.1 第一种边界条件

对第一种边界条件，有三对角矩阵方程组：
$$\left[\begin{array}{ccccc}4& 1& 0& \cdots & 0\\ 1& 4& 1& \cdots & 0\\ ⋮& \ddots & \ddots & \ddots & ⋮\\ 0& \cdots & 1& 4& 1\\ 0& \cdots & 0& 1& 4\end{array}\right]\left[\begin{array}{c}{c}_2\\ c_3\\ ⋮\\ c_{n-1}\\ c_n\end{array}\right]=\left[\begin{array}{c}6y_1-y_0+\frac{h^2}{6}{y}_0^{\prime \prime }\\ 6y_2\\ 6y_3\\ ⋮\\ 6y_{n-1}\\ 6y_{n-1}-y_n+\frac{h^2}{6}{y}_n^{\prime \prime }\end{array}\right]$$
$$\{\begin{array}{l}{c}_0=2c_1-c_2+h^2{y}_0^{\prime \prime }\\ \\ c_1=y_0-\frac{h^2}{6}{y}_0^{\prime \prime }\\ \\ c_{n+1}=y_n-\frac{h^2}{6}{y}_n^{\prime \prime }\\ \\ c_{n+2}=6c_{n-1}-c_n+h^2{y}_n^{\prime \prime }\end{array}$$
求出所有式子后带入B样条为基底的样条函数就得到了样条插值函数。

2.1.2 第二种边界条件

$$\left[\begin{array}{ccccc}4& 2& 0& \cdots & 0\\ 1& 4& 1& \cdots & 0\\ ⋮& \ddots & \ddots & \ddots & ⋮\\ 0& \cdots & 1& 4& 1\\ 0& \cdots & 0& 2& 4\end{array}\right]\left[\begin{array}{c}{c}_1\\ c_2\\ ⋮\\ c_n\\ c_{n+1}\end{array}\right]=\left[\begin{array}{c}6y_0+2h{y}_0^{\prime }\\ 6y_1\\ 6y_2\\ ⋮\\ 6y_{n-1}\\ 6y_n-2h{y}_n\end{array}\right]$$
$$\{\begin{array}{l}{c}_0=c_2-2h{y}_0^{\prime }\\ \\ c_{n+2}=c_n+2h{y}_n^{\prime }\end{array}$$

2.1.3 第三种边界条件

$$\left[\begin{array}{ccccc}4& 1& 0& \cdots & 1\\ 1& 4& 1& \cdots & 0\\ ⋮& \ddots & \ddots & \ddots & ⋮\\ 0& \cdots & 1& 4& 1\\ 1& \cdots & 0& 1& 4\end{array}\right]\left[\begin{array}{c}{c}_2\\ c_3\\ ⋮\\ c_n\\ c_{n+1}\end{array}\right]=\left[\begin{array}{c}6y_1\\ 6y_2\\ 6y_3\\ ⋮\\ 6y_{n-1}\\ 6y_n\end{array}\right]$$
$$\{\begin{array}{l}{c}_{n+2}=c_2\\ \\ c_1=c_{n+1}\\ \\ c_0=c_n\end{array}$$

[!example]-
$$\begin{array}{cccccccc}x& ||& x_0=1& |& x_1=2& |& x_2=3& |\\ =& =& =& =& =& =& =& =\\ f(x)& ||& 2& |& 4& |& 8& |\\ -& -& -& -& -& -& -& -\\ f^{\prime }(x)& ||& 1.3863& |& & |& 5.5452& |\end{array}$$
求三次样条插值函数。使其满足 $s(x_i)=f(x_i),s^{\prime }(x_0)=f^{\prime }(x_0),s^{\prime }(x_2)=f^{\prime }(x_2)$ 。
解：边界条件二，区间上被分成两段， $n=2$
$$\left[\begin{array}{ccc}4& 2& 0\\ 1& 4& 1\\ 0& 2& 4\end{array}\right]\left[\begin{array}{c}{c}_1\\ c_2\\ c_3\end{array}\right]=\left[\begin{array}{c}6y_0+2h{y}_0^{\prime }\\ 6y_1\\ 6y_2-2h{y}_2^{\prime }\end{array}\right]=\left[\begin{array}{c}14.7726\\ 24\\ 36.9096\end{array}\right]$$
$$\left[\begin{array}{c}{c}_1\\ c_2\\ c_3\end{array}\right]=\left[\begin{array}{c}1.84657\\ 3.69317\\ 7.38081\end{array}\right]$$

$$c_0=c_2-2h{y}_0^{\prime }=0.92057$$
$$c_4=c_2+2h{y}_2^{\prime }=14.7836$$
$$\begin{array}{rl}\therefore s(x)=& 0.92057{\Omega }_3(x-0)+1.84657{\Omega }_3(x-1)\\ & +3.69317{\Omega }_3(x-2)+7.38081{\Omega }_3(x-3)+14.7836{\Omega }_3(x-4)\\ & 1\le x\le 3\end{array}$$

2.2 三弯矩法求三次样条插值函数

用于 区间长度 $h$ 不一致的情况
$$\begin{array}{rl}s(x)=& \frac{M_{i-1}}{6h_i}(x_i-x{)}^3+\frac{M_i}{6h_i}(x-x_{i-1}{)}^3\\ & +(\frac{y_{i-1}}{h_i}-\frac{M_{i-1}}{6}{h}_i)(x_i-x)+(\frac{y_i}{h_i}-\frac{M_i}{6}{h}_i)(x-x_i)\\ & x_{i-1}\le x\le x_i,i=1,2,\cdots ,n\end{array}$$
关键在求M。三弯矩方程
$$r_i{M}_{i-1}+2M_i+{\alpha }_i{M}_{i+1}={\beta }_i,i=1,2,\cdots ,n-1$$
$${\alpha }_i=\frac{h_{i+1}}{h_i+h_{i+1}},r_i=1-{\alpha }_i$$
$${\beta }_i=\frac{6}{h_i+h_{i+1}}(\frac{y_{i+1}-y_i}{h_{i+1}}-\frac{y_i-y_{i-1}}{h_i})$$
可以求得
$${\alpha }_1\cdots {\alpha }_{n-1},r_1\cdots r_{n-1},{\beta }_1\cdots {\beta }_{n-1}$$

2.2.1 第一种边界条件

$${\alpha }_0=0,{\beta }_0=2y_0^{\prime \prime },r_n=0,{\beta }_n=2y_n^{\prime \prime }$$

$$\left[\begin{array}{ccccc}2& {\alpha }_0& & & \\ r_1& 2& {\alpha }_1& & \\ & & \ddots & & \\ & & r_{n-1}& 2& {\alpha }_{n-1}\\ & & & r_n& 2\end{array}\right]\left[\begin{array}{c}{M}_0\\ M_1\\ ⋮\\ M_{n-1}\\ M_n\end{array}\right]=\left[\begin{array}{c}{\beta }_0\\ {\beta }_1\\ ⋮\\ {\beta }_{n-1}\\ {\beta }_n\end{array}\right]$$
可得 $M_0,M_1,\cdots ,M_n$ 。

2.2.2 第二种边界条件

$${\alpha }_0=1,r_n=1$$
$${\beta }_0=\frac{6}{h_1}(\frac{y_1-y_0}{h_1}-y_0^{\prime }),{\beta }_n=\frac{6}{h_n}(y_n^{\prime }-\frac{y_n-y_{n-1}}{h_n})$$
代入的矩阵式子和第一种边界条件的式子相同。求得 $M_0$ 到 $M_n$ 。

2.2.3 第三种边界条件

$$M_0=M_n,{\alpha }_n=\frac{h_1}{h_1+h_n},r_n=1-{\alpha }_n$$
$${\beta }_n=\frac{6}{h_1+h_n}(\frac{y_1-y_0}{h_1}-\frac{y_n-y_{n-1}}{h_n})$$
$$\left[\begin{array}{ccccc}2& {\alpha }_1& & & r_1\\ r_2& 2& {\alpha }_2& & \\ & & \ddots & & \\ & & r_{n-1}& 2& {\alpha }_{n-1}\\ {\alpha }_n& & & r_n& 2\end{array}\right]\left[\begin{array}{c}{M}_1\\ M_1\\ ⋮\\ M_{n-1}\\ M_n\end{array}\right]=\left[\begin{array}{c}{\beta }_1\\ {\beta }_1\\ ⋮\\ {\beta }_{n-1}\\ {\beta }_n\end{array}\right]$$

---

下链

---