数学建模-插值算法（Matlab）

注意：代码文件仅供参考，一定不要直接用于自己的数模论文中
国赛对于论文的查重要求非常严格，代码雷同也算作抄袭
如何修改代码避免查重的方法：https://www.bilibili.com/video/av59423231   //清风数学建模

一、基础知识

简单来说是根据已知点进行线性数据预测，但数据太少需要通过数学方法生成一些数据，如根据十年人口数据预测未来五年人口增减情况。

1.两个常用的插值法

一般多用分段三次埃尔米特插值和三次样条插值

其内置函数分别为

p=pchip(x,y,new_x)

x,y为已知样本点横纵坐标，new_x为要插入对应点的横坐标（纵坐标会计算生成，即为p）

三次样条插值函数参数也同理：

 p=spline(x,y,new_x)

2.生成和打印图片

脚本文件中生成图片（在同一个脚本文件里面，要想画多个图，需要给每个图编号，否则只会显示最后一个图）

figure(1) 

plot函数，显示图像，至少要两种点的横纵坐标，多种会生成更多线

线方式：-实线   :点线    -.虚点线   --波折线

点方式： .圆点   +加号   *型号   o小圆

颜色： y黄  r红  g绿  b蓝  k黑

标注使用legend函数，跟在figure、plot之后，一组x，y（即一种点）对应一个名称，最后选择注释所在图中的位置（如东南是右下角）

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-')
legend('样本点','三次埃尔米特插值','Location','SouthEast') 

 3.n维数据的插值

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

 参数还是已知样本点的横纵坐标，new_x为要插入点的横坐标。最后method是调用方法，如下

'linear'：线性插值（默认） 'cubic'：三次插值

'spline'：三次样条插值（多用） 'nearest'：最邻近插值

 例：人口预测，给09-18年的人口，预测19-21年的

population=[133126,133770,134413,135069,135738,136427,137122,137866,138639, 139538];
year = 2009:2018;  %类似于等差数列，是09-18每个数值取整组成的矩阵
%可以加上 p1 = pchip(year, population, 2019:2021) 这样最后就是两条线了 
p2 = interpn(year, population, 2019:2021,'spline') %等价于p2 = spline(year, population, 2019:2021)
figure(4);
plot(year, population,'o',2019:2021,p2,'bx-')
legend('样本点','interpn三次样条插值预测','Location','SouthEast')

 二、实战演练

例题来源：第六届MathorCup大学生数学建模挑战A题 淡水养殖池塘水华发生及池水净化处理

 上图是附件二每个池子一些物质含量的变化，可以主要上沿，每个池子只有奇数周的数据，本演练是把1号池子的每周数据补充完整（即补上偶数周）

根据一号池的数据，把数据保存为Z（忽略叶绿体和汉字），如下图（可以另存为本地）

 多组数据图片的演示，可以for循环处理，同时要用到subplot（m,n,p）

m表示是图排成m行，n表示图排成n列.p表示图所在的位置，如p=1表示从左到右从上到下的第一个位置。

ylabel(ylab{i}) %给y轴标题,这里是直接引用元胞数组（先定义好的）中的字符串哦。

参考优秀代码（Mr.sandman同学提供）：

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的第一行加上周数

 矩阵P的数值即为要求的，可以粘贴回excel