cube spline---三次样条插值

插值是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。与拟合不用经过每个已知点不同，插值需要经过每个已知点，另外并不是阶数越高越好，因为高阶插值容易出现龙格现象，即插值后在区间两端点处波动极大，产生明显的震荡。三次样条插值作为一种常见的插值方法，这里记录一下其基本概念及求解过程。

一、基本概念

设在区间\([a, b]\)上存在\(n+1\)个已知数据点如下，其把\([a, b]\)分成了\(n\)个子区间\([x_0, x_1], [x_1, x_2],\ ... ,\ [x_{n-1}, x_n]\)，

\[ \begin{align} x &: a \le x_0 < x_1 < ...< x_n \le b \nonumber\\ f(x) = y &: \quad\quad y_0\ \quad y_1\ \quad ... \quad y_n \nonumber \end{align} \]

如果函数\(S(x)\)满足以下三个条件，则称\(S(x)\)为\(f(x)\)关于节点\(x_0, x_1, ..., x_n\)的三次样条函数。

在每个子区间上都是一个不超过3次的多项式，即 \(S_i(x_i) = a_i + b_i(x-x_i) + c_i(x-x_i)^2 + d_i(x-x_i)^3\)；
在整个区间\([a,b]\)上连续且光滑，即一阶导数\(S^{'}(x_i)\)和二阶导数\(S^{''}(x_i)\)存在且连续；
满足插值条件，即\(S_i(x_i) = y_i,\ i = 0, 1, 2,..., n\)。

二、求解过程

2.1 获取方程组

\(n\)个子区间中每一个都有四个未知数\((a_i, b_i, c_i, d_i)\)，所以共有\(4n\)个未知参数需要求解。

插值条件

曲线在整个区间\([a, b]\)上所有点都满足插值条件，所以\(S_i(x_i) = y_i,\ i = 0, 1, 2,..., n-1,\ S_{n-1}(x_n) = y_n\)，依此可得到\(n+1\)个方程。

曲线连续

曲线在整个区间\([a, b]\)上连续，表明在端点\(i = 1, 2, ... , n-1\)处两边函数值相等，即：\(S_{i-1}(x_i) = S_{i}(x_{i})\)，其等价于\(S_i(x_{i+1}) = S_{i+1}(x_{i+1}) = y_{i+1}, i = 0, 1, 2, ..., n-2\)，这里可得\(n-1\)个方程。

曲线光滑

曲线在整个区间\([a, b]\)上光滑，表明在端点\(i = 1, 2, ... , n-1\)两边的一阶导数和二阶导数存在且相等，即

\[S^{'}_{i-1}(x_i) = S^{'}_{i}(x_i), i = 1, 2, ..., n-1.\ S^{''}_{i-1}(x_i) = S^{''}_{i}(x_i), i = 1, 2, ..., n-1 \]

其等价于

\[S^{'}_{i}(x_{i+1}) = S^{'}_{i+1}(x_{i+1}), i = 0, 1, ..., n-2.\ S^{''}_{i}(x_{i+1}) = S^{''}_{i+1}(x_{i+1}), i = 0, 1, ..., n-2 \]

共\(2(n-1)\)个方程。

区间\([a, b]\)左右两端点特性

由2.1---2.3一共可以得到\(n+1 + n-1 + 2(n-1) = 4n - 2\)个方程，要求解4n个未知数，还需要至少两个方程，所以这里可以考虑在两端点\(i = 0, n\)处的特性，加上边界处这两个方程可以对\(4n\)个参数进行求解。一般有三种边界条件：自然边界（Natural Spline），固定边界（Clamped Spline），非节点边界（Not-A-Knot Spline）。

自然边界
指定端点二阶导数为0，即\(S^{''}_0(x_0) = S^{''}_{n-1}(x_n)=0\)。
固定边界
人为指定端点一阶导数，这里分别定为\(A\)和\(B\)，即\(S^{'}_0(x_0) = A, S^{'}_{n-1}(x_n) = B\)。
非节点边界
强制第一个插值点的三阶导数值等于第二个点的三阶导数值，最后第一个点的三阶导数值等于倒数第二个点的三阶导数值. 即\(S^{'''}_0(x_0) = S^{'''}_1(x_1), S^{'''}_{n-2}(x_{n-1}) = S^{'''}_{n-1}(x_{n})\)。

2.2 方程推导

由\(S_i(x) = a_i + b_i(x-x_i) + c_i(x-x_i)^2 + d_i(x-x_i)^3\)可得其各阶导数如下：

\(S^{'}_i(x) = b_i + 2c_i(x-x_i) + 3d_i(x-x_i)^2\)

\(S^{''}_i(x) = 2c_i + 6d_i(x-x_i)\)

\(S^{'''}_i(x) = 6d_i\)

1、由\(S_i(x_i) = y_i,\ i = 0, 1, 2,..., n-1\)可得：

\(S_i(x_i) = a_i + b_i(x_i-x_i) + c_i(x_i-x_i)^2 + d_i(x_i-x_i)^3 = y_i\)，化简可得：

\(a_i = y_i \ \ \ \ \ (1)\)

2、由\(S_i(x_{i+1}) = y_{i+1},\ i = 0, 1, 2, ..., n - 2\)可得：

\(S_i(x_{i+1}) = a_i + b_i(x_{i+1}-x_i) + c_i(x_{i+1}-x_i)^2 + d_i(x_{i+1}-x_i)^3 = y_{i+1}\)，令\(h_i = x_{i+1}-x_i\)，则可简写为

\(a_i + b_ih_i + c_ih^2_i + d_ih^3_i = y_{i+1} \ \ \ \ \ (2)\)

3、由\(S^{'}_{i}(x_{i+1}) = S^{'}_{i+1}(x_{i+1}), i = 0, 1, ..., n-2\)可得：

\(S^{'}_i(x_{i+1}) = b_i + 2c_i(x_{i+1}-x_i) + 3d_i(x_{i+1}-x_i)^2 = b_i + 2c_ih_i + 3d_ih^2_i\)

\(S^{'}_{i+1}(x_{i+1}) = b_{i+1} + 2c_{i+1}(x_{i+1}-x_{i+1}) + 3d_{i+1}(x_{i+1}-x_{i+1})^2 = b_{i+1}\)

所以

\(b_i + 2c_ih_i + 3d_ih^2_i - b_{i+1} = 0 \ \ \ \ \ (3)\)

4、由\(S^{''}_{i}(x_{i+1}) = S^{''}_{i+1}(x_{i+1}), i = 0, 1, ..., n-2\)可得

\(2c_i + 6d_ih_i - 2C_{i+1} = 0 \ \ \ \ \ (4)\)

通过(2)(3)(4)式可知\(b_i, c_i, d_i\)存在某种关系，为了方便求解，这里作以下变换，将\(b_i, c_i, d_i\)三个参数变换为求解\(m_i\)一个参数。
由(4)式可知：

\(d_i = \frac{2C_{i+1} - 2C_{i}}{6h_i}\), 令\(2C_i = m_i\),则：

\(d_i = \frac{m_{i+1} - m_{i}}{6h_i}\)

由(2)式可知：

\(a_i + b_ih_i + c_ih^2_i + d_ih^3_i = y_{i+1} = y_i + h_ib_i + h^2_i\frac{m_i}{2} + h^3_i\frac{m_{i+1} - m_i}{6h_i}\)，化简可得：

\(b_i = \frac{y_{i+1} - y_i}{h_i} = \frac{h_i}{2}m_i - \frac{h_i}{6}(m_{i+1} - m_i)\)

将\(b_i, c_i, d_i\)代入(3)式可知：

\(h_im_i + 2(h_i + h_{i+1})m_{i+1} + h_{i+1}m_{i+2} = 6(\frac{y_{i+2} - y_{i+1}}{h_{i+1}} - \frac{y_{i+1} - y_i}{h_i}), i = 0, 1, ..., n-2.\)

通过以上推导可知关于\(m_i\)的一个方程组，维度为\(n-1\)，再考虑边界条件便可以得到\(n+1\)个方程，进而进行求解。
自然边界

由\(S^{''}_0(x_0) = S^{''}_{n-1}(x_n)=0\)可知\(m_0 = m_n = 0\)，因此方程组可以写成以下形式

\[\begin{bmatrix} 1 & 0 & 0 & 0 & 0 & \cdots & 0 \\ h_0 & 2(h_0 + h_1) & h_1 & 0 & 0 & \cdots & 0 \\ 0 & h_1 & 2(h_1 + h_2) & h_2 & 0 & \cdots & 0 \\ 0 & 0 & h_2 & 2(h_2 + h3) & h_3 & \cdots & \vdots \\ \vdots & \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & h_{n-2} & 2(h_{n-2} + h_{n-1}) & h_{n-1} \\ 0 & 0 & 0 & \cdots & 0 & 0 & 1 \\ \end{bmatrix} \begin{bmatrix} m_0 \\ m_1 \\ m_2 \\ m_3 \\ \vdots \\ m_{n-1} \\ m_n \\ \end{bmatrix} = 6\begin{bmatrix} 0 \\ \frac{y_2 - y_1}{h_1} - \frac{y_1 - y_0}{h_0} \\ \frac{y_3 - y_2}{h_2} - \frac{y_2 - y_1}{h_1} \\ \frac{y_4 - y_3}{h_3} - \frac{y_3 - y_2}{h_2} \\ \vdots \\ \frac{y_n - y_{n-1}}{h_{n-1}} - \frac{y_{n-1} - y_{n-2}}{h_{n-2}} \\ 0 \\ \end{bmatrix} \]

固定边界

由\(S^{'}_0(x_0) = A\)可知：

\[ \begin{align} A &= b_0 \nonumber\\ &= \frac{y_1 - y_0}{h_0} - \frac{h_0}{2}m_0 - \frac{h_0}{6}(m_1 - m_0) \nonumber \end{align} \]

变形可得：

\(2h_0m_0 + h_0m_1 = 6[\frac{y_1 - y_0}{h_0} - A]\)

由\(S^{'}_{n-1}(x_n) = B\)可知：

\[ \begin{align} B &= b_{n-1} \nonumber\\ &= \frac{y_n - y_{n-1}}{h_{n-1}} - \frac{h_{n-1}}{2}m_{n-1} - \frac{h_{n-1}}{6}(m_n - m_{n-1}) \nonumber \end{align} \]

所以：

\(h_{n-1}m_{n-1} + 2h_{n-1}m_n = 6[B - \frac{y_n - y_{n-1}}{h_{n-1}}]\)

因此方程组的系数矩阵可以改写为

\[\begin{bmatrix} 2h_0 & h_0 & 0 & 0 & 0 & \cdots & 0 \\ h_0 & 2(h_0 + h_1) & h_1 & 0 & 0 & \cdots & 0 \\ 0 & h_1 & 2(h_1 + h_2) & h_2 & 0 & \cdots & 0 \\ 0 & 0 & h_2 & 2(h_2 + h3) & h_3 & \cdots & \vdots \\ \vdots & \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & h_{n-2} & 2(h_{n-2} + h_{n-1}) & h_{n-1} \\ 0 & 0 & 0 & \cdots & 0 & h_{n-1} & 2h_{n-1} \\ \end{bmatrix} \begin{bmatrix} m_0 \\ m_1 \\ m_2 \\ m_3 \\ \vdots \\ m_{n-1} \\ m_n \\ \end{bmatrix} = 6\begin{bmatrix} \frac{y_1 - y_0}{h_0} - A \\ \frac{y_2 - y_1}{h_1} - \frac{y_1 - y_0}{h_0} \\ \frac{y_3 - y_2}{h_2} - \frac{y_2 - y_1}{h_1} \\ \frac{y_4 - y_3}{h_3} - \frac{y_3 - y_2}{h_2} \\ \vdots \\ \frac{y_n - y_{n-1}}{h_{n-1}} - \frac{y_{n-1} - y_{n-2}}{h_{n-2}} \\ B - \frac{y_n - y_{n-1}}{h_{n-1}} \\ \end{bmatrix} \]

非节点边界

由于\(S^{'''}_0(x_0) = S^{'''}_1(x_1), S^{'''}_{n-2}(x_{n-1}) = S^{'''}_{n-1}(x_{n})\)且\(S^{'''}_{i}(x_{i}) = 6d_i\)，所以：

\(d_0 = d_1,\ d_{n-2} = d_{n-1}\)，另外由\(d_i = \frac{m_{i+1} - m_i}{6h_i}\)，所以

\(h_1(m_1 - m_0) = h_0(m_2 - m_1), \ h_{n-1}(m_{n-1} - m_{n-2}) = h_{n-2}(m_{n} - m_{n-1})\)

变形可得：

\[ \begin{align} -h_1m_0 + (h_0 + h_1)m_1 - h_0m_1 &= 0 \nonumber\\ -h_{n-1}m_{n-2} + (h_{n-2} + h_{n-1})m_{n-1} - h_{n-2}m_n &= 0 \nonumber \end{align} \]

因此方程组的系数矩阵可以改写为

\[\begin{bmatrix} -h_1 & h_0 + h_1 & - h_0 & 0 & 0 & \cdots & 0 \\ h_0 & 2(h_0 + h_1) & h_1 & 0 & 0 & \cdots & 0 \\ 0 & h_1 & 2(h_1 + h_2) & h_2 & 0 & \cdots & 0 \\ 0 & 0 & h_2 & 2(h_2 + h3) & h_3 & \cdots & \vdots \\ \vdots & \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & h_{n-2} & 2(h_{n-2} + h_{n-1}) & h_{n-1} \\ 0 & 0 & 0 & \cdots & -h_{n-1} & h_{n-2} + h_{n-1} & - h_{n-2} \\ \end{bmatrix} \begin{bmatrix} m_0 \\ m_1 \\ m_2 \\ m_3 \\ \vdots \\ m_{n-1} \\ m_n \\ \end{bmatrix} = 6\begin{bmatrix} 0 \\ \frac{y_2 - y_1}{h_1} - \frac{y_1 - y_0}{h_0} \\ \frac{y_3 - y_2}{h_2} - \frac{y_2 - y_1}{h_1} \\ \frac{y_4 - y_3}{h_3} - \frac{y_3 - y_2}{h_2} \\ \vdots \\ \frac{y_n - y_{n-1}}{h_{n-1}} - \frac{y_{n-1} - y_{n-2}}{h_{n-2}} \\ 0 \\ \end{bmatrix} \]

三、算法总结

步骤1：分区间
将数据分成不同子区间\([x_0, x_1], [x_1, x_2],\ ... ,\ [x_{n-1}, x_n]\)
步骤2：计算步长
计算步长\(h_i = x_{i+1} - x_{i}\)
步骤3：求解方程组获得\(m_i\)
步骤4：计算每个子区间的参数\({a_i, b_i, c_i, d_i}\)
\(a_i = y_i\)
\(b_i = \frac{y_{i+1} - y_i}{h_i} - \frac{h_i}{2}m_i - \frac{h_i}{6}(m_{i+1} - m_i)\)
\(c_i = \frac{m_i}{2}\)
\(d_i = \frac{m_{i+1} - m_{i}}{6h_i}\)
步骤5：得到每个区间样条函数
\(S_i(x_i) = a_i + b_i(x-x_i) + c_i(x-x_i)^2 + d_i(x-x_i)^3\)

四、示例结果

样本\((4.0, 4.2), (4.3, 5.7), (4.6, 6.6), (5.3, 4.8), (5.9, 4.6)\)在不同边界情况下插值效果如下：

五、代码实现

参考三次样条插值(Cubic Spline Interpolation)及代码实现(C语言)

点击展开代码

#define S_FUNCTION_NAME  cubic
#define S_FUNCTION_LEVEL 2
#include "simstruc.h"
#include "malloc.h"  //方便使用变量定义数组大小

static void mdlInitializeSizes(SimStruct *S)
{
    /*参数只有一个，是n乘2的定点数组[xi, yi]:
     * [ x1，y1;
     *   x2, y2;
     *   ..., ...;
     *   xn, yn;
    */
    ssSetNumSFcnParams(S, 1); 
    if (ssGetNumSFcnParams(S) != ssGetSFcnParamsCount(S)) return;

    ssSetNumContStates(S, 0);
    ssSetNumDiscStates(S, 0);

    if (!ssSetNumInputPorts(S, 1)) return;  //输入是x
    ssSetInputPortWidth(S, 0, 1);
    ssSetInputPortRequiredContiguous(S, 0, true);
    ssSetInputPortDirectFeedThrough(S, 0, 1);

    if (!ssSetNumOutputPorts(S, 1)) return;  //输出是S(x)
    ssSetOutputPortWidth(S, 0, 1);

    ssSetNumSampleTimes(S, 1);
    ssSetNumRWork(S, 0);
    ssSetNumIWork(S, 0);
    ssSetNumPWork(S, 0);
    ssSetNumModes(S, 0);
    ssSetNumNonsampledZCs(S, 0);

    ssSetSimStateCompliance(S, USE_DEFAULT_SIM_STATE);

    ssSetOptions(S, 0);
}

static void mdlInitializeSampleTimes(SimStruct *S)
{
    ssSetSampleTime(S, 0, CONTINUOUS_SAMPLE_TIME);
    ssSetOffsetTime(S, 0, 0.0);
}



#define MDL_INITIALIZE_CONDITIONS
#if defined(MDL_INITIALIZE_CONDITIONS)
  static void mdlInitializeConditions(SimStruct *S)
  {
  }
#endif



#define MDL_START
#if defined(MDL_START) 
  static void mdlStart(SimStruct *S)
  {
  }
#endif /*  MDL_START */


static void mdlOutputs(SimStruct *S, int_T tid)
{
    const real_T *map = mxGetPr(ssGetSFcnParam(S,0));  //获取定点数据
    const int_T *mapSize = mxGetDimensions(ssGetSFcnParam(S,0));  //定点数组维数
    const real_T *x = (const real_T*) ssGetInputPortSignal(S,0);  //输入x
    real_T       *y = ssGetOutputPortSignal(S,0); //输出y
    int_T step = 0;  //输入x在定点数中的位置
    int_T i;
    real_T yval;

    for (i = 0; i < mapSize[0]; i++)
    {
        if (x[0] >= map[i] && x[0] < map[i + 1])
        {
            step = i;
            break;
        }
    }
    
    cubic_getval(&yval, mapSize, map, x[0], step);
    y[0] = yval;
    
}

//自然边界的三次样条曲线函数
void cubic_getval(real_T* y, const int_T* size, const real_T* map, const real_T x, const int_T step)
{
    int_T n = size[0];
    
    //曲线系数
    real_T* ai = (real_T*)malloc(sizeof(real_T) * (n-1));
    real_T* bi = (real_T*)malloc(sizeof(real_T) * (n-1));
    real_T* ci = (real_T*)malloc(sizeof(real_T) * (n-1));
    real_T* di = (real_T*)malloc(sizeof(real_T) * (n-1));
    
    real_T* h = (real_T*)malloc(sizeof(real_T) * (n-1));  //x的??
    
    /* M矩阵的系数
     *[B0, C0, ...
     *[A1, B1, C1, ...
     *[0,  A2, B2, C2, ...
     *[0, ...             An-1, Bn-1]
     */
    real_T* A = (real_T*)malloc(sizeof(real_T) * (n-2));
    real_T* B = (real_T*)malloc(sizeof(real_T) * (n-2));
    real_T* C = (real_T*)malloc(sizeof(real_T) * (n-2));
    real_T* D = (real_T*)malloc(sizeof(real_T) * (n-2)); //等号右边的常数矩阵
    real_T* E = (real_T*)malloc(sizeof(real_T) * (n-2)); //M矩阵
    
    real_T* M = (real_T*)malloc(sizeof(real_T) * (n));  //包含端点的M矩阵
    
    int_T i;
    
    //计算x的步长
    for ( i = 0; i < n -1; i++)
    {
        h[i] = map[i + 1] - map[i];
    }
    
    //指定系数
    for( i = 0; i< n - 3; i++)
    {
        A[i] = h[i]; //忽略A[0]
        B[i] = 2 * (h[i] + h[i+1]);
        C[i] = h[i+1]; //忽略C(n-1)
    }

    
    //指定常数D
    for (i = 0; i<n - 3; i++)
    {
        D[i] = 6 * ((map[n + i + 2] - map[n + i + 1]) / h[i + 1] - (map[n + i + 1] - map[n + i]) / h[i]);
    }
    
    
    //求解三对角矩阵，结果赋值给E
    TDMA(E, n-3, A, B, C, D);
    
    M[0] = 0; //自然边界的首端M为0
    M[n-1] = 0;  //自然边界的末端M为0
    for(i=1; i<n-1; i++)
    {
        M[i] = E[i-1]; //其它的M值
    }
    
    //?算三次?条曲?的系数
    for( i = 0; i < n-1; i++)
    {
        ai[i] = map[n + i];
        bi[i] = (map[n + i + 1] - map[n + i]) / h[i] - (2 * h[i] * M[i] + h[i] * M[i + 1]) / 6;
        ci[i] = M[i] / 2;
        di[i] = (M[i + 1] - M[i]) / (6 * h[i]);
    }
    
    *y = ai[step] + bi[step]*(x - map[step]) + ci[step] * (x - map[step]) * (x - map[step]) + di[step] * (x - map[step]) * (x - map[step]) * (x - map[step]);
    
    free(h);
    free(A);
    free(B);
    free(C);
    free(D);
    free(E);
    free(M);
    free(ai);
    free(bi);
    free(ci);
    free(di);

}

void TDMA(real_T* X, const int_T n, real_T* A, real_T* B, real_T* C, real_T* D)
{
    int_T i;
    real_T tmp;

    //上三角矩阵
    C[0] = C[0] / B[0];
    D[0] = D[0] / B[0];

    for(i = 1; i<n; i++)
    {
        tmp = (B[i] - A[i] * C[i-1]);
        C[i] = C[i] / tmp;
        D[i] = (D[i] - A[i] * D[i-1]) / tmp;
    }

    //直接求出X的最后一个值
    X[n-1] = D[n-1];

    //逆向迭代， 求出X
    for(i = n-2; i>=0; i--)
    {
        X[i] = D[i] - C[i] * X[i+1];
    }
}


#define MDL_UPDATE 
#if defined(MDL_UPDATE)
  static void mdlUpdate(SimStruct *S, int_T tid)
  {
  }
#endif

#define MDL_DERIVATIVES
#if defined(MDL_DERIVATIVES)
  static void mdlDerivatives(SimStruct *S)
  {
  }
#endif

static void mdlTerminate(SimStruct *S)
{
}

#ifdef  MATLAB_MEX_FILE  
#include "simulink.c"
#else
#include "cg_sfun.h"
#endif

参考链接

三次样条插值(Cubic Spline Interpolation)及代码实现(C语言)
三次样条（cubic spline）插值

posted @ 2022-11-05 21:39 半夜打老虎阅读(736) 评论(0) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

半夜打老虎