三次样条插值算法及推导过程

1、定义

给定 $n + 1$ 个数据点，共有 $n$ 个区间，三次样条方程 $S (n)$ 满足以下条件：在每个分段区间内 $x_i,x_{i+1})$ ( $i = 0, 1, ..., n - 1, x$ 递增)， $S(x)=S_i(x)$ 是一个三次多项式；满足 $S(x_i)=y_i(i=0,1,...,n)$ ；导数 $S^{'} (x)$ 、二阶导数 $S^{''} (x)$ 在区间是连续的，即 $S (x)$ 曲线是光滑的。

那么 $n$ 个三次多项式分段可以写作：
$S_i(x)=a_i+b_i(x-x_i)+c_i(x-x_i)^2+d_i(x-x_i)^3\tag{1}$
$i = 0, 1, ..., n - 1$ 。
其中： $a_i、b_i、c_i、d_i$ 代表 $4 n$ 个未知数。

2、已知条件求解

在这里插入图片描述

3、具体推导

在这里插入图片描述

4、matlab案例

clc;
clear;


P = [4.0,4.2;
     4.3 5.7;
     4.6,6.6;
     5.3,4.8;
     5.9,4.6];

x = P(:,1);
y = P(:,2);

xx = linspace(min(x),max(x),20);
yy = csapi(x,y,xx); % 三次样条插值
plot(x,y,'b*',xx,yy,'r-',LineWidth=2)
title('Cubic Spline Interpolant')

5、案例结果

在这里插入图片描述

6、matlab仿真

4的案例代码是我自己实现的，做实验应付用足够。做研究还是要清楚具体实现过程的，于是乎，需要详细过程实现。仿真代码来自：三次样条Cubic Spline简介。其他博客翻阅了有数百篇，没一个好用的。

clc;
clear;
close;
 
P = [4.0,4.2;
     4.3 5.7;
     4.6,6.6;
     5.3,4.8;
     5.9,4.6];
 
x = P(:,1);
y = P(:,2);
 
% 计算矩阵A,B
n = size(x,1);
h = arr_diff(x);
A = calc_A(n,h);
B = calc_B(n,y,h);
 
% 求解abcd
a = y;
b = zeros(n-1,1);
c = A^-1*B;
d = zeros(n-1,1);
 
for i=1:n-1
    d(i) = (c(i+1)-c(i))/(3*h(i));
    b(i) = (a(i+1)-a(i))/h(i) - h(i)*(c(i+1)+2*c(i))/3;
end
 
% 绘制原始点
plot(x,y,'b*',LineWidth=2);hold on
 
% 三次样条曲线点采样
accuracy = 0.1;
num = int32((max(x) - min(x))/accuracy) + 1;
pnt = zeros(num,2);
i = 1;
for v = min(x):accuracy:max(x)
    idx = find_latest_index(v,x);
    ai = a(idx);
    bi = b(idx);
    ci = c(idx);
    di = d(idx);
    dx = v - x(idx);
    w = ((di*dx + ci)*dx + bi)*dx + ai;
    pnt(i,1) = v;
    pnt(i,2) = w;
    i = i+1;
end
plot(pnt(:,1),pnt(:,2),'r-',LineWidth=2);
 
% 求差值函数
function h = arr_diff(a)
    sz = size(a,1);
    h = zeros(sz-1,1);
    for i=1:sz-1
        h(i) = a(i+1)-a(i);
    end
end
 
% 计算A矩阵
function A = calc_A(n,h)
    A = zeros(n,n);
    A(1,1) = 1;
    for i=1:n-1
        if i ~= n-1
            A(i+1,i+1) = 2*(h(i) + h(i+1));
        end
        A(i,i+1) = h(i);
        A(i+1,i) = h(i);
    end
    A(1,2) = 0;
    A(n,n-1) = 0;
    A(n,n) = 1;
end
 
% 计算B矩阵
function B = calc_B(n,y,h)
    B = zeros(n,1);
    for i = 1:n-2
        B(i+1) = 3*((y(i+2)-y(i+1))/h(i+1) - (y(i+1)-y(i))/h(i));
    end
end
 
% 计算s中不大于si的值的最大索引
function idx = find_latest_index(si,s)
    [row,~] = size(s);
    if si <= 0 || row <= 1
       idx = 1;
       return;
    end
    for i=2:row
        if(s(i-1,1) <= si && si <= s(i,1))
             idx = i-1;
             return;
        end
    end
    idx = row;
end