自然三次样条

三次样条：

三次样条是使用满足给定m 个控制点的三次多项式的样条。为了推导出三次样条的解，我们假设端点处的第二个推导为 0，这反过来提供了一个边界条件，将两个方程添加到m-2方程以使其可解。一维三次样条的方程组可由下式给出：

$S_i (x) = a_i + b_i *x + c_i * x^{2} + d_i x^{3}$

我们为函数y = f(x)取一组点 [x _i , y _i ] for i = 0, 1, ..., n 。三次样条插值是一条分段连续曲线，通过表中的每个值。

以下是 K=3 次样条的条件：
- s 的域在 [a, b] 的区间内。
- S, S', S" 都是 [a,b] 上的连续函数。

$S(x) =\begin{bmatrix} S_0 (x), x \epsilon [x_0, x_1] \\ S_1 (x), x \epsilon [x_1, x_2] \\ ... \\ ... \\ ... \\ S_{n-1} (x), x \epsilon [x_1, x_2] \end{bmatrix}$

这里S _i (x)是将在子区间[x _i , x _i+1 ]上使用的三次多项式。

由于每个方程有 n 个区间和 4 个系数，因此我们需要4n 个参数来求解样条，我们可以从每个三次样条方程都满足两端值的事实得到 2n 个方程：

$S_i(x_i) = y_i \\ S_i(x_{i+1}) = y_{i+1}$

上述三次样条方程不仅应该是连续的和可微的，而且还定义了在控制点上也是连续的一阶和二阶导数。

$S_{i-1}^{'}(x_i) = S_{i}^{'}(x_i)$

和

$S_{i-1}^{''}(x_i) = S_{i}^{''}(x_i)$

对于1, 2, 3…n-1提供2n -2方程约束。因此，我们还需要 2 个方程来求解上述三次样条。为此，我们将使用一些自然边界条件。

自然三次样条：

在自然三次样条中，我们假设样条在边界点处的二阶导数为 0：

$S^{''}(x_0) = 0 \\ S^{''}(x_n) = 0$

现在，由于 S(x) 是一个三阶多项式，我们知道 S”(x) 是一个插值的线性样条。因此，首先，我们构造 S”(x)，然后对其进行两次积分以获得 S(x)。

现在，让我们假设 t_i = x_i for i= 0, 1,…n，并且 $z_i = S^{''}(x_i), i=0,1,2..n$ 并且从自然边界条件 $z_0= z_n =0$ .二次微分三次样条给出线性样条，可以写成：

$S_i^{''}(x) = z_i \frac{x- t_{i+1}}{t_i - t_{i+1}} + z_{i+1} \frac{x - t_i}{t_{i+1} - t_i}$

在哪里，

$h_i = t_{i+1} - t_i ; t=0,1, 2,...,n$

现在，等式变为：

$S''(x) = z_{i+1} \frac{x- t_i}{h_i} + z_{i} \frac{t_{i+1} - x}{h_i}$

将该方程积分两次以获得三次样条：

$S(x) =\frac{z_{i+1}}{6h_i}(x - t_i)^{3} + \frac{z_{i}}{6h_i}( t_{i+1} - x)^{3} + C_i(x_i -t) + D_i(t_{i+1} -x)$

在哪里，

$C_i = \frac{y_{i+1}}{h_i} - \frac{z_{i+1} * h_i}{6} \\ D_i = \frac{y_{i}}{h_i} - \frac{z_{i} * h_i}{6}$

现在，检查 t_i 处导数的连续性，即 S^{'}_{i}(t_i) = S^{'}_{i-1}(t_{i})。我们首先需要找到导数并设置条件：

$S'(x) = \frac{z_{i+1}}{2h_i}(x - t_i)^{2} - \frac{z_{i}}{2h_i}( t_{i+1} - x)^{2} + \frac{1}{h_i}(y_{i+1}-y_i)- \frac{h_i}{6}(z_{i+1} -z_i)$

在哪里，

$b_i = \frac{1}{h_i}(y_{i+1}-y_i)$

将上述方程用于连续性并求解它给出以下方程：

$6(b_i -b_{i-1}) = h_{i-1}z_{i-1} + 2(h_{i-1} + h_i)z_i + h_{i}z_{i+1}$

取，v_i = 6(b_i – b_{i-1}) 上面的等式可以写成矩阵的形式：

$\begin{bmatrix} 1& 0 & 0& .& .& .& .& .& .& .& .& 0& \\ h_0& 2(h_0 + h_1)& h_1& & & & & & & & & .& \\ .& h_1 & 2(h_1 + h_2)& h_2& & & & & & & & .& \\ .& .& .& .& & & & & & & & .& \\ .& & & .& .& .& & & & & & .& \\ .& & & & .& .& .& & & & & .& \\ .& & & & & .& .& .& & & & .& \\ .& & & & & & .& .& .& & & .& \\ .& & & & & & & & .& h_{n-2}& 2(h_{n-1} + h_{n-2})& h_{n-1}& \\ 0& .& .& .& .& .& .& .& .& 0& 0& 1& \\ \end{bmatrix} \begin{bmatrix} z_0 \\ z_1\\ z_2\\ .\\ .\\ .\\ .\\ .\\ z_{n-1}\\ z_{n} \end{bmatrix} = \begin{bmatrix} 0\\ v_1 \\ .\\ .\\ .\\ .\\ .\\ .\\ v_{n-1}\\ 0 \end{bmatrix}$

执行

在这个实现中，我们将使用满足自然边界条件的三次样条对函数f(x) = 1/x的点 b/w 2-10 执行样条插值。

Python3

#imports
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import CubicSpline, interp1d
plt.rcParams['figure.figsize'] =(12,8)
  
x = np.arange(2,10)
y = 1/(x)
# apply cubic spline interpolation
cs = CubicSpline(x, y, extrapolate=True)
# apply natural cubic spline interpolation
ns = CubicSpline(x, y,bc_type='natural', extrapolate=True)
  
# Apply Linear interpolation
linear_int = interp1d(x,y)
  
xs = np.arange(2, 9, 0.1)
ys = linear_int(xs)
  
# plot linear interpolation
plt.plot(x, y,'o', label='data')
plt.plot(xs,ys,  label="interpolation", color='green')
plt.legend(loc='upper right', ncol=2)
plt.title('Linear Interpolation')
plt.show()
  
# define a new xs
xs = np.arange(1,15)
#plot cubic spline and natural cubic spline
plt.plot(x, y, 'o', label='data')
plt.plot(xs, 1/(xs), label='true')
plt.plot(xs, cs(xs), label="S")
plt.plot(xs, ns(xs), label="NS")
plt.plot(xs, ns(xs,2), label="NS''")
plt.plot(xs, ns(xs,1), label="NS'")
  
plt.legend(loc='upper right', ncol=2)
plt.title('Cubic/Natural Cubic Spline Interpolation')
plt.show()
  
# check for boundary condition
print("Value of double differentiation at boundary conditions are %s and %s"
      %(ns(2,2),ns(10,2)))

输出：

Value of double differentiation at boundary conditions are -1.1102230246251565e-16 and -0.00450262550915248

线性插值

三次/自然样条插值

参考：

自然三次样条 NTNU 大学