插值算法

学习视频：【强烈推荐】清风：数学建模算法、编程和写作培训的视频课程以及Matlab

老师讲得很详细，很受用！！！

作用

数模比赛中，常常需要根据已知的函数点进行数据、模型的处理和分析，而有时候现有的数据是极少的，不足以支撑分析的进行，这时就需要使用一些数学的方法，“模拟产生”一些新的但又比较靠谱的值来满足需求，这就是插值的作用，另一个不常见的作用就是短期预测。

一维插值问题

定义

方法分类

本文重点介绍数学建模常用的两种方法：三次样条插值和分段三次埃尔米特插值

插值多项式

原理

拉格朗日插值法

方法主要有拉格朗日插值法，具体不介绍，它会出现龙波现象（在两端处波动极大,产生明显的震荡）。

但觉得可以出一个ACM题。

分段线性插值

• 插值多项式次数高精度未必显著提高
• 插值多项式次数越高摄入误差可能显著增大

如何提高插值精度呢
采用分段低次插值是一种办法

概念

牛顿插值法

评价

与拉格朗日插值法相比，牛顿插 值法的计算过程具有继承性。 （牛顿插值法每次插值只和前n项 的值有关，这样每次只要在原来 的函数上添加新的项，就能够产 生新的函数） 但是牛顿插值也存在龙格现象的 问题。

致命缺点

上面讲的两种插值仅仅要求插值多项式在插值节点处与被插函数有相等的函数值，而这种插值多项式却不能全面反映被插值函数的性态。
然而在许多实际问题中，不仅要求插值函数与被插值函数在所有节点处有相同的函数值，它也需要在一个或全部节点上插值多项式与被插函数有相同的低阶甚至高阶的导数值。
对于这些情况，拉格朗日插值和牛顿插值都不能满足

埃尔米特(Hermite)插值

概念

保持播值曲线在节点处有切线（光滑），使插值函数和被插函数的密和程度更好。

​ 不但要求在节点上的函数值相等，而且还要求对应的导数值也相等,甚至要求 高阶导数也相等，满足这种要求的插值多项式就是埃尔米特插值多项式。

缺点与改进

​ 直接使用Hermite插值得到的多项式次数较高，也存在着龙格现象， 因此在实际应用中，往往使用分段三次Hermite插值多项式(PCHIP)。

​ Matlab有内置的函数（实现过程已经帮我们封装好了，会调用就行了）： p = pchip(x,y,new_x)

​ x是已知的样本点的横坐标;y是已知的样本点的纵坐标;new_x是要插入处对应的横坐标

三次样条插值

Matlab有内置的函数： p = spline(x,y,new_x)

x是已知的样本点的横坐标;y是已知的样本点的纵坐标;new_x是要插入处对应的横坐标

代码

% 分段三次埃尔米特插值
x = -pi:pi; y = sin(x); 
new_x = -pi:0.1:pi;
p = pchip(x,y,new_x);
figure(1); % 在同一个脚本文件里面，要想画多个图，需要给每个图编号，否则只会显示最后一个图哦~
plot(x, y, 'o', new_x, p, 'r-')

% plot函数用法:
% plot(x1,y1,x2,y2) 
% 线方式： - 实线 :点线 -. 虚点线 - - 波折线 
% 点方式： . 圆点  +加号  * 星号  x x形  o 小圆
% 颜色： y黄； r红； g绿； b蓝； w白； k黑； m紫； c青


% 三次样条插值和分段三次埃尔米特插值的对比
x = -pi:pi; 
y = sin(x); 
new_x = -pi:0.1:pi;
p1 = pchip(x,y,new_x);   %分段三次埃尔米特插值
p2 = spline(x,y,new_x);  %三次样条插值
figure(2);
plot(x,y,'o',new_x,p1,'r-',new_x,p2,'b-')
legend('样本点','三次埃尔米特插值','三次样条插值','Location','SouthEast')   %标注显示在东南方向
% 说明：
% LEGEND(string1,string2,string3, …)
% 分别将字符串1、字符串2、字符串3……标注到图中，每个字符串对应的图标为画图时的图标。
% ‘Location’用来指定标注显示的位置


% n维数据的插值
% p = interpn(x1,x2,...,xn, y, new_x1,newx_2,...,new_xn, method) 
% x1,x2,...,xn是已知的样本点的横坐标 y是已知的样本点的纵坐标坐标 
% new_x1,newx_2,...,new_xn是要插入点的横坐标 method是要插值的方法 
%‘linear’：线性插值（默认算法）； 
%‘cubic’：三次插值； 
%‘spline’： 三次样条插值法； ( 最为精准 ) 
%‘nearest：最邻近插值算法。 
x = -pi:pi; y = sin(x); 
new_x = -pi:0.1:pi;
p = interpn (x, y, new_x, 'spline');
% 等价于 p = spline(x, y, new_x);
figure(3);
plot(x, y, 'o', new_x, p, 'r-')

小技巧：短期预测

根据过去10年的中国人口数据，预测接下来三年的人口数据

% 人口预测（注意：一般我们很少使用插值算法来预测数据，随着课程的深入，后面的章节会有更适合预测的算法供大家选择，例如灰色预测、拟合预测等）
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) %三次样条插值预测
figure(4);
plot(year, population,'o',2019:2021,p1,'r*-',2019:2021,p2,'bx-')
legend('样本点','三次埃尔米特插值预测','三次样条插值预测','Location','SouthEast')

建模实例

MathorCup第六届A题淡水养殖池塘水华发生及池水净化处理

%插值预测中间周的水体评价指标
load Z.mat
x=Z(1,:); %Z的第一行是星期Z: 1     3     5     7     9    11    13    15
[n,m]=size(Z);%n为Z的行数，m为Z的列数
% 注意Matlab的数组中不能保存字符串，如果要生成字符串数组，就需要使用元胞数组，其用大括号{}定义和引用
ylab={'周数','轮虫','溶氧','COD','水温','PH值','盐度','透明度','总碱度','氯离子','透明度','生物量'};  % 等会要画的图形的标签
disp(['共有' num2str(n-1) '个指标要进行插值。'])
disp('正在对一号池三次埃尔米特插值，请等待')%一号池共有十一组要插值的数据，算上星期所在的第一行，共十二行
P=zeros(11,15);%对要储存数据的矩阵P赋予初值
for i=2:n%从第二行开始都是要进行插值的指标
    y=Z(i,:);%将每一行依次赋值给y
    new_x=1:15;%要进行插值的x
    p1=pchip(x,y,new_x);%调用三次埃尔米特插值函数
    subplot(4,3,i-1);%将所有图依次变现在4*3的一幅大图上
    plot(x,y,'ro',new_x,p1,'-');%画出每次循环处理后的图像
    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的第一行加上周数