数模笔记——三、插值算法

（一）适用问题

要根据已知的函数点进行数据、模型的处理和分析时，数据过少不足以支撑分析的进行，需要用一些数学的方法，模拟产生一些新的但又比较靠谱的值来满足需求。

注意：高次插值会产生龙格现象，即在两端处波动极大，产生明显的震荡。在不熟悉运动趋势的前提下，不要轻易使用高次插值。

为提高插值精度可以选择使用分段低次插值的方法。

（二）运用Matlab中的内置函数实现插值

1.Matlab中绘图函数的用法

Matlab中的绘图函数：plot(x1,y1,x2,y2)

线方法：- 实线 : 点线 -. 虚点线 – 波折线

点方法：. 圆点 + 加号 * 星号 x 交叉 o 圆圈

颜色：y 黄 r 红 g 绿 b 蓝 w 白 k 黑 m 紫 c 青

2.Matlab中求插值的两个内置函数

分段三次埃尔米特插值：p=pchip(x,y,new_x)

三次样条插值：p=spline(x,y,new_x)

x：已知样本点的横坐标

y：已知样本点的纵坐标

new_x：要插入处对应的横坐标

3.两种插值方法的比较

例如：

分别用两种方法求出$sin(x)$函数中在区间$[-\pi ,\pi ]$​中间隔为0.1的插值

x = -pi : pi;%从-pi到pi，默认步长为1
y = sin(x);
new_x = -pi : 0.1 : pi;%从-pi到pi，步长为0.1
p1 = pchip(x , y , new_x);%分段三次埃尔米特插值
p2 = spline(x , y , new_x);%三次样条插值
figure(1) %命名防止被覆盖
plot(x , y , 'o' , new_x , p1 , 'r-' , new_x , p2 , 'b-')
legend('样本点','三次埃尔米特插值','三次样条插值','Location','SouthEast')%标注显示在东南方向

结果如下：

可以看到，三次样条生成的曲线更光滑，在实际建模中，由于我们不知道数据的生成过程，因此这两种插值均可以使用

补充：n维数据插值

p = interpn(x1, x2, ... ,xn, y, new_x1, new_x2, ... , new_xn, method)

xi：已知样本点的横坐标

yi：已知样本点的纵坐标

new_xi：要插入处对应的横坐标

method：要插值的方法

如：

“linear” : 线性插值（默认算法）

“cubic” : 三次插值

“spline” : 三次样条插值（最精确）

“nearest” : 最邻近插值算法

（三）插值算法的运用

1. 预测

例如：根据过去10年的中国人口数据，预测接下来3年的数据

Matlab代码

population = [133126, 133770, 134413, 135069, 135738, 136427, 137122, 137866, 138639, 139538];
year = 2009 : 2018;
p1 = pchip(year, population, 2019 : 2021); %分段三次埃尔米特插值预测
p2 = spline(year, population, 2019 : 2021); %三次样条插值预测
plot(year, population, 'o', 2019 : 2021, p1, 'r*-', 2019 : 2021, p2, 'bx-')
legend('样本点','三次埃尔米特插值预测','三次样条插值预测','Location','SouthEast');

预测结果：

(但在实际建模中一般不用插值算法预测)

2.补齐数据

例如：如图为某一池塘各水体评价指标的单周数据，请补齐第1到15周的数据

Matlab代码：

%用插值算法预测中间周的水体评价指标
load Z.mat
x = Z(1, :); %Z的第一行（第1到15的奇数周）
[n, m] = size(Z); %Z的行数n，列数m
ylab = {'周数','轮虫','溶氧','COD','水温','PH值','盐度','透明度','总碱度','氯离子','透明度','生物量'};  % 每个图对应的y
disp(['共有' num2str(n - 1) '个指标要进行插值']) %Z的第2到第12行为要插值的指标
disp(['正在进行三次埃尔米特插值……'])

P = zeros(11, 15); %用一个11行15列的矩阵存储插值后的数据并初始化矩阵为0
for i = 2 : n %Z的第2到第12行为要插值的指标
    y = Z(i, :); %将每个指标的数据赋给y（y为一个行向量）
    new_x = 1: 15; %要进行插值的x
    p1 = pchip(x, y, new_x); %调用三次埃尔米特插值函数
    subplot(4, 3, i - 1); %将当前图窗分成4*3的网格，并在第i-1个位置创建坐标轴
    plot(x, y, 'ro', new_x, p1, 'b-'); %画出图像
    axis([0 15, -inf, inf]); %横坐标为0到15，纵坐标不变
    xlabel('星期') %x轴标题
    ylabel(ylab{i}) %y轴标题
    P(i - 1, :) = p1; %将每次插值的结果保存在P矩阵中
end
legend('原始数据','三次埃尔米特插值数据', 'Location', 'SouthEast') %标注
P = [1:15; P] %在P的第一行加上周数

补齐结果：