登山第八梯：三次样条拟合——滑溜的很

一 简介

四 代码演示

        此处以固定边界为例，进行代码演示

function [x2,y2]=Spline_3(x,y)  
n=len(x);
if n<=3
    disp('数组长度过短!')
    return
end
if n~=len(y)
    disp('两个数组长度必须相同！')
    return
end
 
if nargin==2
    dx=zeros(1,2);
    dx(1)=(y(2)-y(1))/(x(2)-x(1));
    dx(2)=(y(n)-y(n-1))/(x(n)-x(n-1));
    %指定末端点的导数值
else
    return;
end
%========预分配内存空间===============================
h=zeros(n-1,1);  %步长
delta_y=zeros(n-1,1);   %增量
%解对应的系数
A=zeros(n-1,1);
B=zeros(n-1,1);
C=zeros(n-1,1);
D=zeros(n-1,1);
%G和f为需要解的n*n矩阵及其数值
G=zeros(n,n);  %是 n*n 矩阵,右端有n个二阶导数值
f=zeros(n,1);  % f是右边的矩阵
%=========================
 
for i=1:1:n-1
    h(i)=x(i+1)-x(i);
    delta_y(i)=y(i+1)-y(i); %每次y的增量
end
%   A是y(i)而B是导数
 
%if nargin==2     %这个留在后续可能要加三弯矩时使用
%采用三转角构造方式
    G(1,1)=2*h(1);
    G(1,2)=h(1);
    %n个数据共有n个导数值，这里取第n个来列方程（b(n)是倒数第一个）
    %参考资料的b(n-1)有误，应当是b(n)
    G(n,n-1)=h(n-1); 
    G(n,n)=2*h(n-1);
    a=dx(1);
    b=dx(2);
%=============构造矩阵=====================
    %填写G矩阵中其余的参数()
    for i=2:1:n-1  %从2行到n-1行的数据
        G(i,i-1)=h(i-1);
        G(i,i)=2*(h(i)+h(i-1));
        G(i,i+1)=h(i);
    end
    
    f(1)=6*( delta_y(1)/h(1) - a  ) ;    
    %在默认情况下，采用最前方的斜率代替导数值时，f(1)必然为0
    f(n)=6*(  b-delta_y(n-1)/h(n-1)  );
    %同理默认情况下  f（n-1）=0
    for i = 2:1:n-1
        f(i)=6*(delta_y(i)/h(i)-delta_y(i-1)/h(i));
    end
    %m为二阶导数值
    m=G\f;   
    %是G\f（实际上是inv(G)*f  ）解出m 的值（这个看错检查我好长时间QAQ）
    m=m';
    
%=====================================================
    %计算每段样条插值的系数：
    for i=1:1:n-1
        A(i)=y(i);
        B(i)=  (y(i+1)-y(i))/h(i) -h(i)*m(i)/2 - h(i)/6*(m(i+1)-m(i));
        C(i)=m(i)/2;
        D(i)=(m(i+1)-m(i))/(6*h(i));
    end
    
%*******注意:从这里设置插值后绘制函数的精度*****可以有不同的要求,
%这里设为0.001
    x_2=min(x):0.001:max(x);  %如果用x2会有问题
    x_2=x_2';   %行向量
    %绘制插值曲线
    y2=zeros(len(x_2),1);
    j=1;

        
    for i=1:1:len(x_2)
        while x_2(i)>x(j+1)  %变换插值函数
            j=j+1;
        end
        y2(i)=A(j)   +B(j)*(x_2(i)-x(j))  +C(j)*(x_2(i)-x(j))^2  +  D(j)*(x_2(i)-x(j))^3;
    end
 
%end    %三转角构造方式
x2=x_2;
 
 
end

clear all;
close all;
clc;
t_raw = 0:0.01:10;
y_raw = cos(t_raw);
t=0:1:10;
y=cos(t);
[t2,y2]=Spline_3(t,y)
plot(t2,y2,'r-')
hold on 
plot(t_raw,y_raw,'g-')