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

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-1d(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-1h(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-1if i ~= n-1A(i+1,i+1) = 2*(h(i) + h(i+1));endA(i,i+1) = h(i);A(i+1,i) = h(i);endA(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-2B(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 <= 1idx = 1;return;endfor i=2:rowif(s(i-1,1) <= si && si <= s(i,1))idx = i-1;return;endendidx = row;
end