Matlab三次样条插值：如何获取超数据最大值点并可视化

已编写Matlab代码构建三次自然样条并绘制数据图，但不知道如何在图中显示原数据组外的点（例如f(2010)）。希望在2000年后的有效区间（如t=2010）实现，但不知从何入手。

原代码 #

clear; clc;

t= [1850, 1875, 1900, 1925, 1950, 1975, 2000];
y= [285.2, 288.6, 295.7, 305.3, 311.3, 331.36, 369.64];

N= length(t); %number of points I want 
n=N-1 ; % number of subintervals

h=(t(N)-t(1))/n; %step size

A=[1,1,1,0],B=[2,0,0,0,2],C=[0,1,1,1];

Trid=diag(4*ones(1,n-1))+diag(A,-1)+diag(B)+diag(C,1);

for i=1:n-1
    f(i)= 6/h^2*(y(i+2)-2*y(i+1)+y(i));
end

f=f';
w=inv(Trid)*f;%since sigma 1 and sigma n+1 are both 0, we need to add 0 in the beginning and also in the end of then matrix
sigma=[0;w;0];%it is a nx1 matrix, be careful.

for i=1:n
    d(i)=y(i);
    b(i)=sigma(i)/2;
    a(i)=(sigma(i+1)-sigma(i))/(6*h);
    c(i)=(y(i+1)-y(i))/h-h/6*(2*sigma(i)+sigma(i+1));
end

r= 25; %subsubintervals for t ex. between 1850 and 1875, here i seperate it into 1 years per slot
hh=h/r; %step size of subsubintervals

x=t(1):hh:t(N);

for i=1:n
    for j=r*(i-1)+1:r*i
        s(j)=a(i)*(x(j)-t(i))^3+b(i)*(x(j)-t(i))^2+c(i)*(x(j)-t(i))+d(i);
    end
end

s(r*n+1)=y(N);

plot(t,y,'o')
hold on
plot(x,s,'-x')
hold off

解决方案 #

要计算并显示样条在原数据范围外的点，核心是利用最后一段样条的多项式表达式进行短距离外推（注意：样条外推可靠性有限，仅建议在靠近原区间的范围内使用）。具体实现如下：

1. 计算目标外推点的值 #

以t=2010为例，它属于最后一个子区间（1975-2000）的延伸范围，直接用最后一段的样条系数计算：

target_t = 2010;
% 最后一段的系数为a(n), b(n), c(n), d(n)，对应区间t(n)到t(n+1)
target_s = a(n)*(target_t - t(n))^3 + b(n)*(target_t - t(n))^2 + c(n)*(target_t - t(n)) + d(n);

2. 扩展样条曲线的绘制范围 #

修改原代码中x的范围，让它覆盖到2010，同时用逻辑索引替代嵌套循环计算样条值（更简洁高效）：

% 扩展x到2010
x = t(1):hh:2010;
s = zeros(size(x)); % 预分配内存提升性能

% 计算原区间内的样条值
for i=1:n
    idx = x >= t(i) & x <= t(i+1);
    s(idx) = a(i)*(x(idx)-t(i)).^3 + b(i)*(x(idx)-t(i)).^2 + c(i)*(x(idx)-t(i)) + d(i);
end

% 计算2000-2010外推区间的样条值
idx_extra = x > t(N);
s(idx_extra) = a(n)*(x(idx_extra)-t(n)).^3 + b(n)*(x(idx_extra)-t(n)).^2 + c(n)*(x(idx_extra)-t(n)) + d(n);

完整修改后的代码 #

clear; clc;

t= [1850, 1875, 1900, 1925, 1950, 1975, 2000];
y= [285.2, 288.6, 295.7, 305.3, 311.3, 331.36, 369.64];

N= length(t); 
n=N-1 ; 

h=(t(N)-t(1))/n; 

A=[1,1,1,0]; B=[2,0,0,0,2]; C=[0,1,1,1];

Trid=diag(4*ones(1,n-1))+diag(A,-1)+diag(B)+diag(C,1);

for i=1:n-1
    f(i)= 6/h^2*(y(i+2)-2*y(i+1)+y(i));
end

f=f';
w=inv(Trid)*f;
sigma=[0;w;0];

for i=1:n
    d(i)=y(i);
    b(i)=sigma(i)/2;
    a(i)=(sigma(i+1)-sigma(i))/(6*h);
    c(i)=(y(i+1)-y(i))/h-h/6*(2*sigma(i)+sigma(i+1));
end

r= 25; 
hh=h/r; 

% 扩展x到2010
x=t(1):hh:2010;
s = zeros(size(x)); % 预分配内存

% 计算原区间内的样条值
for i=1:n
    idx = x >= t(i) & x <= t(i+1);
    s(idx) = a(i)*(x(idx)-t(i)).^3 + b(i)*(x(idx)-t(i)).^2 + c(i)*(x(idx)-t(i)) + d(i);
end

% 计算外推区间的样条值
idx_extra = x > t(N);
s(idx_extra) = a(n)*(x(idx_extra)-t(n)).^3 + b(n)*(x(idx_extra)-t(n)).^2 + c(n)*(x(idx_extra)-t(n)) + d(n);

% 计算并标记2010的点
target_t = 2010;
target_s = a(n)*(target_t - t(n))^3 + b(n)*(target_t - t(n))^2 + c(n)*(target_t - t(n)) + d(n);

plot(t,y,'o','DisplayName','原始数据')
hold on
plot(x,s,'-x','DisplayName','样条曲线')
plot(target_t, target_s, 'rs','MarkerSize',10,'DisplayName','f(2010)')
legend
hold off
xlabel('年份')
ylabel('数值')
title('三次自然样条曲线及外推点')

关键说明 #

• 样条外推是基于最后一段多项式的延伸，超出原数据区间越远，结果可靠性越低，仅适合短距离外推。
• 用逻辑索引替代嵌套循环，代码更简洁，同时预分配s的内存可提升运行效率。
• 新增图例、坐标轴标签和标题，让图表信息更清晰直观。

内容的提问来源于stack exchange，提问作者Teh Ais Kaw