三次样条插值原理+Matlab代码

三次样条插值原理+Matlab代码

1.1 分段插值

1.1.1 龙格现象

例1-1: 给定函数 $f(x)=\frac{1}{1+25x^2}$, $x\in [-1,1]$，构造10次拉格朗日插值多项式 $L_{10}(x)$。

解: 对 $[-1,1]$ 作等距分割，取 $h=\frac{2}{10}=0$, $x_i=-1+0.2i$，取插值点为 $(x_i,\frac{1}{1+25x^2})$, $i=0,1,\cdots ,10$。解得 $L_{10}(x)$ 如图所示。

龙格(Runge)现象揭示了高次插值多项式的缺陷.它说明高次多项式的插值效果不一定由于低次多项式的插值效果,同时表明等距插值不能保证有较好的插值效果.为此，我们引入了分段低次插值多项式来提高插值函数的逼近效果.

1.1.2分段线性插值

对给定区间 $[a,b]$ 作分割，$a=x_0<x_1<\cdots <x_n=b$。在每个小区间 $[x_i,x_{i+1}]$ 上作 $f(x)$ 以 { x i , x i + 2 } \{x_i,x_{i+2}\} {xi​,xi+2​} 为节点的线性插值，记插值函数为 $p_i(x)$，则 $p_i(x)=\frac{x-x_{i+1}}{x_i-x_{i+1}}f(x_i)+\frac{x-x_i}{x_{i+1}-x_i}f(x_{i+1})$，$x_i\le x\le x_{i+1}$，
则 $p(x)$ 为：
$$\{\begin{array}{ll}{p}_0(x)=\frac{x-x_1}{x_0-x_1}f(x_0)+\frac{x-x_0}{x_1-x_0}f(x_1),& x\in [x_0,x_1]\\ p_1(x)=\frac{x-x_2}{x_1-x_2}f(x_1)+\frac{x-x_1}{x_2-x_1}f(x_2),& x\in [x_1,x_2]\\ ⋮& ⋮\\ p_{n-1}(x)=\frac{x-x_n}{x_{n-1}-x_n}f(x_{n-1})+\frac{x-x_{n-1}}{x_n-x_{n-1}}f(x_n),& x\in [x_{n-1},x_n]\end{array}$$
对于例1-1,应用分段线性插值的结果如图所示。

分段线性插值算法简单，只要区间充分小，就能保证它的误差要求。它的一个显著特点是它的局部性，如果修改了某个节点 $(x_i,f(x_i))$ 的值，仅在相邻的两个区间 $[x_{i-1},x_i]$ , $[x_i,x_{i+1}]$ 受到影响，分段线性插值的缺点是在插值节点处不光滑。

1.2 三次样条插值

为了解决分段线性插值在插值节点处不光滑的问题，我们在每段上使用3次多项式进行分段，这样在节点处一阶和二阶连续可导。这样就可以使得插值函数 ( S(x) ) 的光滑程度高，而且保证二阶导函数的连续性。
定义1.2 给定区间 $[a,b]$ 上$n+1$ 个节点$a=x_0<x_1<\cdots <x_n=b$ 和这些节点上的函数值 { f ( x i ) = y i , i = 0 , 1 , ⋯   , n } \{ f(x_i) = y_i, i=0,1,\cdots,n \} {f(xi​)=yi​,i=0,1,⋯,n}，若 (S(x))满足： { S ( x i ) = y i , i = 0 , 1 , ⋯   , n } \{ S(x_i) = y_i, i=0,1,\cdots,n \} {S(xi​)=yi​,i=0,1,⋯,n}，$S(x)$ 在每个小区间 $[x_i,x_{i+1}]$ 上是一个三次多项式；$S(x)$在$[a,b]$ 上有连续的二阶导数，则称$S(x)$为$f(x)$关于在区间 $[a,b]$ 上的三次样条插值函数，称 { x 0 , x 1 , ⋯   , x n } \{ x_0,x_1,\cdots,x_n \} {x0​,x1​,⋯,xn​}为样条节点。

公式推导：
记 { M i = S ′ ′ ( x ) ,   m i = S ′ ( x ) ,   i = 0 , 1 , ⋯   , n } \{ M_i=S''(x),\ m_i=S'(x),\ i=0,1,\cdots,n \} {Mi​=S′′(x), mi​=S′(x), i=0,1,⋯,n}，$S_i(x)$为区间$[x_i,x_{i+1}]$上$S(x)$的表达式。由于$S_i(x)$是一个三次多项式，则$S_i^{\prime \prime }(x)$是一次多项式。用插值点$\{(x_i,S_i^{\prime \prime }(x_i)),(x_i,S_i^{\prime \prime }(x_{i+1}))\}$作线性插值。则：
$$S_i^{\prime \prime }(x)=\frac{x-x_{i+1}}{x_i-x_{i+1}}{M}_i+\frac{x-x_i}{x_{i+1}-x_i}{M}_{i+1},x\in [x_i,x_{i+1}].$$
对$S_i^{\prime \prime }(x)$进行两次积分，即：
$$S_i^{\prime }(x)=\frac{(x-x_{i+1}{)}^2}{2(x_i-x_{i+1})}{M}_i+\frac{(x-x_i{)}^2}{2(x_{i+1}-x_i)}{M}_{i+1}+c,x\in [x_i,x_{i+1}].$$
$$S_i(x)=\frac{(x-x_{i+1}{)}^3}{6(x_i-x_{i+1})}{M}_i+\frac{(x-x_i{)}^3}{6(x_{i+1}-x_i)}{M}_{i+1}+cx+d,x\in [x_i,x_{i+1}].$$
令$cx+d=C(x_{i+1}-x)+D(x-x_i)$，则：
$$S_i(x)=\frac{(x-x_{i+1}{)}^3}{6(x_i-x_{i+1})}{M}_i+\frac{(x-x_i{)}^3}{6(x_{i+1}-x_i)}{M}_{i+1}+C(x_{i+1}-x)+D(x-x_i),x\in [x_i,x_{i+1}].$$
记$h_i=x_{i+1}-x_i$，则：
$$S_i(x)=\frac{(x_{i+1}-x{)}^3}{6h_i}{M}_i+\frac{(x-x_i{)}^3}{6h_i}{M}_{i+1}+C(x_{i+1}-x)+D(x-x_i),x\in [x_i,x_{i+1}].$$
又因为$S_i(x_i)=y_i,S_i(x_{i+1})=y_{i+1}$，代入解得，$C=\frac{y_i}{h_i}-\frac{h_i{M}_i}{6}$，$D=\frac{y_{i+1}}{h_i}-\frac{h_i{M}_{i+1}}{6}$，则：
$$S_i(x)=\frac{(x_{i+1}-x{)}^3{M}_i+(x-x_i{)}^3{M}_{i+1}}{6h_i}+\frac{(x_{i+1}-x)y_i+(x-x_i)y_{i+1}}{h_i}-\frac{h_i}{6}\left[(x_{i+1}-x)M_i+(x-x_i)M_{i+1}\right],x\in [x_i,x_{i+1}].$$
所以，
$$S_i^{\prime }(x)=-\frac{(x-x_{i+1}{)}^2}{2h_i}{M}_i+\frac{(x-x_i{)}^2}{2h_i}{M}_{i+1}+\frac{y_{i+1}}{h_i}-\frac{h_i{M}_{i+1}}{6}-\frac{y_i}{h_i}+\frac{h_i{M}_i}{6}.$$

又因为 $x$ 在内节点 $x_i$ 处连续，即 $S_i^{\prime }(x_i)=S_{i-1}^{\prime }(x_i)$ ，可以得到：
$$-\frac{(x_i-x_{i+1}{)}^2}{2h_i}{M}_i+\frac{y_{i+1}}{h_i}-\frac{h_i{M}_{i+1}}{6}-\frac{y_i}{h_i}+\frac{h_i{M}_i}{6}=\frac{(x-x_{i-1}{)}^2}{2h_{i-1}}{M}_i+\frac{y_i}{h_{i-1}}-\frac{h_{i-1}{M}_i}{6}-\frac{y_{i-1}}{h_{i-1}}+\frac{h_{i-1}{M}_{i-1}}{6}$$
整理得：
$$f[x_i,x_{i+1}]-\frac{h_i{M}_i}{3}-\frac{h_i{M}_{i+1}}{6}=f[x_{i-1},x_i]+\frac{h_{i-1}{M}_i}{3}+\frac{h_{i-1}{M}_{i-1}}{6}$$
其中 $f[x_i,x_{i+1}]=\frac{y_{i+1}-y_i}{x_{i+1}-x_i}$。整理得：
6 { f [ x i , x i + 1 ] − f [ x i − 1 , x i ] } = 2 ( h i − 1 + h i ) M i + h i − 1 M i − 1 + h i M i + 1 6\{ f[x_i,x_{i+1}] - f[x_{i-1},x_i] \} = 2(h_{i-1}+h_i)M_i + h_{i-1}M_{i-1} + h_iM_{i+1} 6{f[xi​,xi+1​]−f[xi−1​,xi​]}=2(hi−1​+hi​)Mi​+hi−1​Mi−1​+hi​Mi+1​
两边同除 $(h_{i-1}-h_i)$ ，得：
6 { f [ x i , x i + 1 ] − f [ x i − 1 , x i ] } ( h i − 1 + h i ) = 2 M i + h i − 1 ( h i − 1 + h i ) M i − 1 + h i ( h i − 1 + h i ) M i + 1 \frac{6\{ f[x_i,x_{i+1}] - f[x_{i-1},x_i] \}}{(h_{i-1}+h_i)} = 2M_i + \frac{h_{i-1}}{(h_{i-1}+h_i)}M_{i-1} + \frac{h_i}{(h_{i-1}+h_i)}M_{i+1} (hi−1​+hi​)6{f[xi​,xi+1​]−f[xi−1​,xi​]}​=2Mi​+(hi−1​+hi​)hi−1​​Mi−1​+(hi−1​+hi​)hi​​Mi+1​
令 ${\lambda }_i=\frac{h_i}{(h_{i-1}+h_i)}$，${\mu }_i=1-{\lambda }_i=\frac{h_{i-1}}{(h_{i-1}+h_i)}$， d i = 6 { f [ x i , x i + 1 ] − f [ x i − 1 , x i ] } ( h i − 1 + h i ) = 6 f [ x i − 1 , x i , x i + 1 ] d_i = \frac{6\{ f[x_i,x_{i+1}] - f[x_{i-1},x_i] \}}{(h_{i-1}+h_i)} = 6f[x_{i-1},x_i,x_{i+1}] di​=(hi−1​+hi​)6{f[xi​,xi+1​]−f[xi−1​,xi​]}​=6f[xi−1​,xi​,xi+1​]，则
$${\mu }_i{M}_{i-1}+2M_i+{\lambda }_i{M}_{i+1}=d_i=6f[x_{i-1},x_i,x_{i+1}],i=1,2,\cdots ,n-1$$
上式是关于 $M_i$ （$i=0,1,\cdots ,n$）的线性方程组，共有 $n-1$ 个方程，比未知参数 $n+1$ 少2，因而还需要两个方程才可能唯一确定 $S(x)$。
一般通过对边界附加适当的条件，即用边界条件给出所需要的两方程。常用的有以下三种边界条件：

(1) 已知 $S^{\prime \prime }(x_0)=M_0$ 和 $S^{\prime \prime }(x_n)=M_n$，此时 $n-1$ 个方程组有 $n-1$ 个未知量 { M i , i = 1 , 2 , ⋯   , n − 1 } \{ M_i, i=1,2,\cdots,n-1 \} {Mi​,i=1,2,⋯,n−1}，
当 $M_0=0,M_n=0$ 时，称为自然边界条件。

$$\left[\begin{array}{ccccc}2& {\lambda }_1& & & \\ {\mu }_2& 2& {\lambda }_2& & \\ & \ddots & \ddots & \ddots & \\ & & {\mu }_{n-2}& 2& {\lambda }_{n-2}\\ & & & {\mu }_{n-1}& 2\end{array}\right]\left[\begin{array}{c}{M}_1\\ M_2\\ ⋮\\ M_{n-2}\\ M_{n-1}\end{array}\right]=\left[\begin{array}{c}6f[x_0,x_1,x_2]-{\mu }_1{M}_0\\ 6f[x_1,x_2,x_3]\\ ⋮\\ 6f[x_{n-3},x_{n-2},x_{n-1}]\\ 6f[x_{n-2},x_{n-1},x_n]-{\lambda }_{n-1}{M}_n\end{array}\right]=\left[\begin{array}{c}{d}_1-{\mu }_1{M}_0\\ d_2\\ ⋮\\ d_{n-2}\\ d_{n-1}-{\lambda }_{n-1}{M}_n\end{array}\right].$$

(2) 已知 $S^{\prime }(x_0)=m_0$ 和 $S^{\prime }(x_n)=m_n$，代入 $S^{\prime }(x)$ 在 $[x_0,x_1]$ 和 $[x_{n-1},x_n]$ 的表达式种，即

$$\{\begin{array}{l}2M_0+M_1=\frac{6}{h_0}[f[x_0,x_1]-m_0]=d_0\\ M_{n-1}+2M_n=\frac{6}{h_{n-1}}[m_n-f[x_{n-1},x_n]]=d_n\end{array},$$得到 $n+1$ 个方程, $n+1$ 个未知量。
$$\left[\begin{array}{ccccc}2& 1& & & \\ {\mu }_1& 2& {\lambda }_1& & \\ & \ddots & \ddots & \ddots & \\ & & {\mu }_{n-1}& 2& {\lambda }_{n-1}\\ & & & 1& 2\end{array}\right]\left[\begin{array}{c}{M}_0\\ M_1\\ ⋮\\ M_{n-1}\\ M_n\end{array}\right]=\left[\begin{array}{c}\frac{6}{h_0}(f[x_0,x_1]-m_0)\\ 6f[x_0,x_1,x_2]\\ ⋮\\ 6f[x_{n-2},x_{n-1},x_n]\\ \frac{6}{h_{n-1}}(m_n-f[x_{n-1},x_n])\end{array}\right]=\left[\begin{array}{c}{d}_0\\ d_1\\ ⋮\\ d_{n-1}\\ d_n\end{array}\right].$$

(3) 已知 $S(x)$ 是以 $b-a$ 为周期（或 $x_n-x_0$）的周期函数，即 $f(x_0)=f(x_n)$，$S^{(k)}(x_0)=S^{(k)}(x_n)$，$(k=0,1,2)$。

取 $m_0=m_n$，$M_0=M_n$，将 $S^{\prime }(x_0)=S^{\prime }(x_n)$ 加入方程组，此时化为 $n$ 个变量, $n$ 个方程的方程组。

令 ${\lambda }_0=\frac{h_0}{h_{-1}+h_0}=\frac{h_0}{h_0+h_{n-1}}$, ${\mu }_0=1-{\lambda }_0=\frac{h_{n-1}}{h_0+h_{n-1}}$, $d_0=\frac{6}{h_0+h_{n-1}}[f[x_0,x_1]-f[x_{-1},x_0]]=6f[x_{-1},x_0,x_1]$。
$$\left[\begin{array}{ccccc}2& {\lambda }_0& & & {\mu }_0\\ {\mu }_1& 2& {\lambda }_1& & \\ & \ddots & \ddots & \ddots & \\ & & {\mu }_{n-2}& 2& {\lambda }_{n-2}\\ {\lambda }_{n-1}& & & {\mu }_{n-1}& 2\end{array}\right]\left[\begin{array}{c}{M}_0\\ M_1\\ ⋮\\ M_{n-2}\\ M_{n-1}\end{array}\right]=\left[\begin{array}{c}6f[x_{-1},x_0,x_1]\\ 6f[x_0,x_1,x_2]\\ ⋮\\ 6f[x_{n-3},x_{n-2},x_{n-1}]\\ 6f[x_{n-2},x_{n-1},x_n]\end{array}\right]=\left[\begin{array}{c}{d}_0\\ d_1\\ ⋮\\ d_{n-2}\\ d_{n-1}\end{array}\right].$$

原稿截图

参考

[1] 数值计算方法.马东升,董宁.
[2] 数值计算方法与算法(第三版).张韵华.等.
[3] 高等工程数学(第4版).于寅.

后续

会增加一些例题，完善Matlab代码。(2024/11/24已完成)