欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页

matlab cubic spline interpolation 三次样条插值

程序员文章站 2023-12-27 18:16:39
...

在*上看到别人自己写出的三次样条插值函数MATLAB代码,先搬运下。

链接在这里:https://*.com/questions/7642921/cubic-spline-program

这里进行了验证,C语言程序运行结果同MATLAB端的SPLINE(x,y,xq)函数结果相同。

以下是MATLAB中的代码:

function [s0,s1,s2,s3]=cubic_spline(x,y) %x,y表示待插值点序列

if any(size(x) ~= size(y)) || size(x,2) ~= 1
   error('inputs x and y must be column vectors of equal length');
end

n = length(x);

h = x(2:n) - x(1:n-1);
d = (y(2:n) - y(1:n-1))./h;

lower = h(2:end);
main  = 2*(h(1:end-1) + h(2:end));
upper = h(1:end-1);

T = spdiags([lower main upper], [-1 0 1], n-2, n-2);
rhs = 6*(d(2:end)-d(1:end-1));

m = T\rhs;

% Use natural boundary conditions where second derivative
% is zero at the endpoints

m = [ 0; m; 0];

s0 = y;
s1 = d - h.*(2*m(1:end-1) + m(2:end))/6;
s2 = m/2;
s3 =(m(2:end)-m(1:end-1))./(6*h);

端点采用固定边界,s0,s1,s2,s3为三次函数中的系数,带入下个函数进行插值运算:

function plot_cubic_spline(x,s0,s1,s2,s3)

n = length(x);

inner_points = [28 29 30 20 21];

for i=1:n-1
   xx = linspace(x(i),x(i+1),inner_points(i));
   xi = repmat(x(i),1,inner_points(i));
   yy = s0(i) + s1(i)*(xx-xi) + ... 
     s2(i)*(xx-xi).^2 + s3(i)*(xx - xi).^3;
   plot(xx,yy,'b')
   plot(x(i),0,'r');
end

这里我的插值个数,由inner_points数组给出,每个间隔插值点数不同,yy为三次函数,即插值点对应的Y值。

至于转成C语言代码,已经有人给出了,可以点击:C语言版的三次样条插值,讲解的很详细,我就是这么看懂的。

上一篇:

下一篇: