【数模学习笔记】插值算法

级数展开

已于 2024-01-28 14:59:04 修改

阅读量674

点赞数 7

分类专栏：数模学习笔记文章标签：学习笔记 matlab 数学建模

于 2024-01-28 14:53:54 首次发布

本文链接：https://blog.csdn.net/2301_79920063/article/details/135893942

版权

数模学习笔记专栏收录该内容

3 篇文章 0 订阅

订阅专栏

数模中主要使用三次埃尔米特插值和三次样条插值

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

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

三次埃尔米特插值

例如插值sinx的函数图像

x = -pi:pi; % 原始x, [-pi, pi]，步长为1
y = sin(x); % y和x必须对应
ins_x = -pi: 0.01 : pi; % 插值后的x，[-pi, pi]，步长为0.01
P1 = pchip(x, y, ins_x); % 插值得到的ins_x对应的y
plot(x, y, 'redo', ins_x, P1, 'black-') % 画图，原始数据画红色圆圈，插值后数据画黑色曲线
legend('原始数据', '三次埃尔米特插值', 'Location', 'southeast') % 添加图例到右下角

三次样条插值

就只是换了个函数，我们同样的对sinx进行插值

x = -pi:pi; % 原始x, [-pi, pi]，步长为1
y = sin(x); % y和x必须对应
ins_x = -pi: 0.01 : pi; % 插值后的x，[-pi, pi]，步长为0.01
P2 = spline(x, y, ins_x); % 插值得到的ins_x对应的y
plot(x, y, 'redo', ins_x, P2, 'blue-') % 画图，原始数据画红色圆圈，插值后数据画蓝色曲线
legend('原始数据', '三次样条插值', 'Location', 'southeast') % 添加图例到右下角

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

对人口进行短期预测

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

population = [133126,133770,134413,135069,135738,136427,137122,137866,138639, 139538]
year = [2009:2018]

population = [133126,133770,134413,135069,135738,136427,137122,137866,138639, 139538]; %初始y
year = [2009:2018]; % 初始x
ins_x = [2019:2023]; % 预测的年份
p1 = pchip(year, population, ins_x); % 三次埃尔米特预测
p2 = spline(year, population, ins_x); % 三次样条预测
real_p = [141008,141212,141260,141175,140967]; % 实际的人口数据
plot(year, population, 'o', ins_x, p1, 'r-*', ins_x, p2, 'b-x', ins_x, real_p, 'black-o');
legend('样本点', '三次埃尔米特插值预测', '三次样条插值预测', '实际人口数据', 'location', 'southeast');

可以看到用插值法做短期预测还是可以的，长期预测误差就太大了

建模实例

由于附件二缺少偶数周的数据，我们可以用插值法进行补充

我们的目标是得到1号池的插值图像和插值后的矩阵，由于这里有多个对象进行同样的操作，我们可以考虑用for循环来写代码，这是矩阵数据

Z=[1,3,5,7,9,11,13,15;
1913,1945,1920,2205,2260,2302,2385,2420;
5.12,3.2,6.72,3.36,2.4,4.14,6.43,4.6;
21.9,20,26.8,27.73,23.4,22.75,25.36,26.03;
24.8,25.7,26.8,28,30.4,30,27.6,30.8;
9.31,9.14,9.14,9.29,9.22,9.33,9.16,9.26;
1.8,2.3,1.9,2.1,2.1,1.1,1.5,1.5;
28,24,26,22,22,20,19,23;
425.11,457.99,492.48,492.08,501.93,598.48,604.44,623.89;
628.1,639.2,648.87,640.33,616.43,614.72,507.14,580;
28,24,26,22,22,20,19,23;
30.58,36.19,49.75,60.58,56.58,60.06,67.99,67.74]

clear;clc
Z=[1,3,5,7,9,11,13,15;
1913,1945,1920,2205,2260,2302,2385,2420;
5.12,3.2,6.72,3.36,2.4,4.14,6.43,4.6;
21.9,20,26.8,27.73,23.4,22.75,25.36,26.03;
24.8,25.7,26.8,28,30.4,30,27.6,30.8;
9.31,9.14,9.14,9.29,9.22,9.33,9.16,9.26;
1.8,2.3,1.9,2.1,2.1,1.1,1.5,1.5;
28,24,26,22,22,20,19,23;
425.11,457.99,492.48,492.08,501.93,598.48,604.44,623.89;
628.1,639.2,648.87,640.33,616.43,614.72,507.14,580;
28,24,26,22,22,20,19,23;
30.58,36.19,49.75,60.58,56.58,60.06,67.99,67.74];


init_x = Z(1, :); % 取得原始数据的x
ins_x = [1:15]; % 插值后的x
ylab={'轮虫','溶氧','COD','水温','PH值','盐度','透明度','总碱度','氯离子','透明度','生物量'}; % 画图用的标签
[row, col] = size(Z);
Z_ready = ones(row, 15); % 创建一个矩阵接收数据
Z_ready(1, :) = ins_x; % 将第一行置成[1:15]
for i = 2: row % 第二行开始就是对象了，所以从第二行开始
    p1 = pchip(init_x, Z(i, :), ins_x)
    subplot(4, 3, i-1); % 创建一个大窗口，里面有4*3个小图，对第i-1个图进行编辑
    plot(init_x, Z(i, :), 'ro', ins_x, p1, 'black-'); % 画图
    xlabel('周数'); % 设置x轴名称
    ylabel(ylab(i-1)); % 设置y轴名称
    axis([0, 15, -inf, inf]); % 设置x，y轴的范围
    Z_ready(i, :) = p1; % 将处理好的数据放进Z_ready矩阵中
end
legend('原始数据', '三次埃尔米特插值', 'Location', 'Southeast'); % 设置图例