数学建模第六天：数学建模算法篇之插值及MATLAB实现

计算机鬼才～

已于 2023-04-27 14:32:33 修改

阅读量854

点赞数 3

分类专栏：数学建模两周速成文章标签： html 前端

于 2023-04-26 20:15:17 首次发布

本文链接：https://blog.csdn.net/m0_69299019/article/details/130380661

版权

数学建模两周速成专栏收录该内容

8 篇文章 18 订阅

订阅专栏

一、前言

1、引例

伍老师最近苦不堪言，最近胡吃海喝，管不住嘴，感觉自己最近长胖了，于是给自己制定了到2023年6月要减到45kg的目标（已知现在是2023年的3月）。于是伍老师兴致冲冲的找到了会数学建模的程序员小李，让小李帮她分析，这个月要减到多重才算完成目标。

小李一看，来活儿了，为了体现出自己的专业素养，于是，他向伍老师索要了近两年的伍老师的体重数据。伍老师反手一个大嘴巴子，并且说：“真有你的，竟然敢问女生的体重，不知道这很不礼貌吗？”
最后，伍老师还是妥协了，给了小李自己这两年的相关体重，如下表：(体重单位：kg)

年月	2022.3	2022.7	2022.8	2022.11	2023.1	2023.3	2023.6
体重	52.5	54	53	47	48.4	？	45

其余的月份，伍老师要么是没有称过体重，要么就是称过体重，但是数据没有留存。请大家思考一下，如果你是小李，你怎么帮伍老师预测体重？

2、插值与拟合模型

插值：求过已知有限个数据点的近似函数。

拟合：已知有限个数据点，求近似函数，不要求过已知数据点，只要求在某种意义下它在这些点上的总偏差最小。

插值和拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。而面对一个实际问题，究竟应该用插值还是拟合，有时容易确定，有时则并不明显。

二、插值

1、插值相关定义

其中已知有n+1个节点 $(x_{j},y_{j})$ （其中j=0,1,2,3……,n， $x_{j}$ 各不相同），求任意插值点 $x^{*}$ ( $\neq x_{j}$ )处的插值 $y^{*}$ 。

我们需要构造一个(相对简单的)函数y=f(x)，通过全部节点，即 $f(x_{j})=y_{j}$ ，（j=0,1,2……n），再用f(x)计算插值。即 $f(x^{*})=y^{*}$

2、拉格朗日插值

已知函数f(x)在n+1个点 $x_{0},x_{1},\cdot \cdot \cdot \cdot \cdot x_{n}$ 处的函数值为 $y_{0},y_{1},\cdot \cdot \cdot \cdot \cdot y_{n}$ 。求一n次多项式函数 $P_{n}(x)$ ，使其满足： $P_{n}(x_{i})=y_{i}$ ，i=1,2,3……n。

解决此问题的拉格朗日插值多项式公式为：

其中L（x）为n次多项式，它的表达式为：

它也被称为拉格朗插值基函数。我们以例题为例，算一算使用拉格朗日插值，得到伍老师五月的预测体重。我们以2022.3为起点，得到新的表格（设女性体重到18岁趋于稳定，而伍老师在2022年三月刚好经历成年后的第30个月）：

年月	30	34	35	38	40	42	45
体重	52.5	54	53	47	48.4	？	45

x=42
z1=(x-34)*(x-35)*(x-38)*(x-40)*(x-45)
t1=(0-4)*(0-5)*(0-8)*(0-10)*(0-15)
z2=(x-30)*(x-35)*(x-38)*(x-40)*(x-45)
t2=(4-0)*(4-5)*(4-8)*(4-10)*(4-15)
z3=(x-30)*(x-34)*(x-38)*(x-40)*(x-45)
t3=(5-0)*(5-4)*(5-8)*(5-10)*(5-15)
z4=(x-30)*(x-34)*(x-35)*(x-40)*(x-45)
t4=(8-0)*(8-4)*(8-5)*(8-10)*(8-15)
z5=(x-30)*(x-34)*(x-35)*(x-38)*(x-45)
t5=(10-0)*(10-4)*(10-5)*(10-8)*(10-15)
z6=(x-30)*(x-34)*(x-35)*(x-38)*(x-40)
t6=(15-0)*(15-4)*(15-5)*(15-8)*(15-10)
z=52.5*z1/t1+54*z2/t2+53*z3/t3+47*z4/t4+48.4*z5/t5+45*z6/t6

不建议写成一个完整的等式，因为这样算不出正确的结果，而且运气差的话还会得到一个相当大的负数或者正数。上述代码的结果为：

也就是说，伍老师可以胡吃海喝，随便长胖。但是当小李兴致冲冲的把结果交给伍老师时，伍老师反手又给了小李一个大嘴巴子：“当我憨呢？别人减肥前跑步锻炼，你这个程序员的建议却是让我胡吃海喝？”

所以说，用拉格朗日插值法，不一定能求出最好的结果，而且也有些许不合理。比如说，请大家用上述代码预测一下伍老师在2023年7月的体重大家得出的结果如下：

直接瘦成纸片人，很不合理，对吧？具体原因是因为，再研究插值问题的初期，所有人都想当然的认为，插值多项式的次数越高，插值精度越高。Runge通过对一个例子的研究发现，上述结论仅仅在插值多项式的次数不超过七次时成立；插值多项式的次数超过七时，插值多项式会出现严重的震荡现象，称为Runge现象。而我们的例子是五次的，虽不像七次那么多，但依然有这种震荡现象

3、分段线性插值

相比于拉格朗日插值，分段线性插值要容易许多，大家看下图就明白了：

其函数表达式为：

其特点为计算量与n无关，但是n越大，误差越小。

4、matlab实现

函数：

yi=interp1(x，y，xi，'method')

各个参数说明：

①x，y为插值点，通常为向量

②xi、yi为被插值点，通常为向量

③method为插值方法，有nearest最邻近插值、linear线性插值、spline三次样条插值、cubic立方插值这几种。

我们用例题用三次样条插值测试一下其合理性：

x=[30 34 35 38 40 45];
y=[52.5 54 53 47 48.4 45];
x1=42;
y1=interp1(x,y,x1,'spline')
xi=30:1/1000:45;
yi=interp1(x,y,xi,'spline');
plot(x,y,'o',xi,yi)

y1输出的结果为51.35kg，结果图如下：

这个故事告诉我们，体重反弹并不可怕，因为这很有可能是成功减肥的前兆！

5、二维插值及matlab实现

伍老师为了更好的减肥，决定去爬武当山，刚翻过了一个山头就已经开始气喘吁吁了，于是伍老师边喘气边问：“这山，爬的，好，好累，到底有，，多高啊？”已知武当山在某些位置的高度，请预测一下，这座山的最高处有多高？

已知，横坐标方向每隔400米，纵坐标方向每隔400米相应的海拔如下：

我们先画出大致的3D图，代码与图形如下：

x=0:400:4400;
y=0:400:2800;
z=[370 470 550 600 670 690 670 620 580 450 400 300 ;
  510 620 730 800 870 850 780 720 650 500 200 300 ;
  650 760 880 970 1020 1050 1020 830 900 700 300 500;
   740 880 1080 1130 1250 1280 1230 1040 900 500 700 780;
   830 980 1180 1320 1450 1420 1400 1300 700 900 850 840;
    950 1190 1370 1380 1410 1430 1450 1470 1080 940 780 620;
   930 930 960 1100 1250 1300 1200  980 850 750 550 500 ;
   870 890 910 950 1000 1110 1080 820 690 540 380 300];


mesh(x,y,z)
rotate3

相应的插值函数为：

z=interp2(x0,y0,z0,x,y,’method’)

各个参数说明：

①x0，y0、z0为插值点，通常为矩阵

②x、y为被插值点，通常为向量

③method为插值方法，有nearest最邻近插值、linear线性插值、cubic立方插值这几种。

则我们使用相应的代码段（精度：10米）来预测一下，最高峰有多少米：

x=0:400:4400;
y=0:400:2800;
z=[370 470 550 600 670 690 670 620 580 450 400 300 ;
  510 620 730 800 870 850 780 720 650 500 200 300 ;
  650 760 880 970 1020 1050 1020 830 900 700 300 500;
   740 880 1080 1130 1250 1280 1230 1040 900 500 700 780;
   830 980 1180 1320 1450 1420 1400 1300 700 900 850 840;
    950 1190 1370 1380 1410 1430 1450 1470 1080 940 780 620;
   930 930 960 1100 1250 1300 1200  980 850 750 550 500 ;
   870 890 910 950 1000 1110 1080 820 690 540 380 300];
mesh(x,y,z)
xi=linspace(0,10,4400);  %加密横坐标数据到450个
yi=linspace(0,10,2800);  %加密纵坐标数据到290个
[xii,yii]=meshgrid(xi,yi);  %生成网格数据
zii=interp2(x,y,z,xii,yii,'cubic');  %插值
mesh(xii,yii,zii)  %加密后的地貌图
hold on     % 保持图形
[xx,yy]=meshgrid(x,y);  %生成网格数据
plot3(xx,yy,z+0.1,'ob')  %原始数据用‘O’绘出

跑完该代码之后，观察一下峰值，就知道了最高海拔是多少了！

好啦，本期的数学建模课就到此结束啦，感兴趣的观众老爷们给小编一个小心心好吗？