数学建模插值法与曲线拟合讲课

关于数学建模插值法与曲线拟合讲课第一页，共六十五页，2022年，8月28日

一、问题的提出

在生产和实验中，关于函数f(x)，经常存在两种情况：（1）其表达式不便于计算；（2）无表达式.

而只有函数在给定点的函数值，怎样预测其它点的函数值？xx0x1x2…xnyy0y1y2…yn第二页，共六十五页，2022年，8月28日飞机机翼制造

下表给出的x、y数据位于机翼端面的轮廓线上，Y1和Y2分别对应轮廓的上下线。假设需要得到x坐标每改变0.1时的y坐标，试完成加工所需数据，画出曲线.x035791112131415Y101.82.22.73.03.12.92.52.01.6Y201.21.72.02.02.01.81.21.01.6第三页，共六十五页，2022年，8月28日山体地貌要在某山区方圆大约27平方公里范围内修建一条公路，从山脚出发经过一个居民区，再到达一个矿区。横向纵向分别每隔400米测量一次，得到一些地点的高程：试做出该山区的地貌图.第四页，共六十五页，2022年，8月28日船在该海域会搁浅吗？---作业

在某海域测得一些点(x,y)处的水深z由下表给出，船的吃水深度为5英尺，在矩形区域（75，200）*（-50，150）里的哪些地方船要避免进入.第五页，共六十五页，2022年，8月28日

水深和流速的问题

在水文数据测量中，不同水深的流速是不同的.水文数据的测量时天天进行的，为了减少测量的工作，希望得到确定的水深和水流之间的关系.为此测量了一系列不同水深和流速值.下表给出了对某河流的测量数据，其中水深和流速根据适当的单位进行了规范化，共10个值.第六页，共六十五页，2022年，8月28日美国人口问题据美国人口普查局数据：从1790每隔10年至2000年的总人口（单位：百万）如下示

t=1790:10:2000;p=[3.9,5.3,7.2,9.6,12.9,17.1,23.1,31.4,38.6,50.2,62.9,76,92,105.7,122.8,131.7,150.7,179,205,226.5,251.4,281.422];

预测2001，2002年的美国人口数？并与调查数据285.318，288.369比较，选择拟合较好的模型。第七页，共六十五页，2022年，8月28日农作物施肥效果分析1992年A题

在农业生产试验研究中，对某地区土豆的产量与化肥的关系做了一实验，得到了氮肥、磷肥的施肥量与土豆产量的对应关系如下表：

1.根据上表数据分别给出土豆产量与氮、磷肥的关系式。

2.施肥问题优化策略氮肥量（公斤/公顷）03467101135202259336404471土豆产量（公斤）15.1821.3625.7232.293439.4543.1543.4640.8330.75磷肥量（公斤/公顷）024497398147196245294342土豆产量（公斤）33.4632.4736.0637.964140.141。342.240.442.7第八页，共六十五页，2022年，8月28日

配药方案---作业

一种新药用于临床之前,必须设计给药方案.在快速静脉注射的给药方式下,所谓给药方案是指,每次注射剂量多大,间隔时间多长.

药物进入机体后随血液输送到全身,在这个过程中不断地被吸收,分布,代谢,最终排出体外.药物在血液中的浓度,即单位体积血液中的药物含量,称血药浓度.在最简单的一室模型中,将整个机体看作一个房室,称中心室,室内的血药浓度是均匀的.快速静脉注射后,浓度立即上升;然后逐渐下降.当浓度太低时,达不到预期的治疗效果;血药浓度太高,又可能导致药物中毒或副作用太强.临床上,每种药物有一个最小有效浓度c1和一个最大治疗浓度c2.设计给药方案时,要使血药浓度保持在c1-c2之间.设本题所研究药物的最小有效浓度c1=10,最大治疗浓度c2=25().第九页，共六十五页，2022年，8月28日

显然,要设计给药方案,必须知道给药后血药浓度随时间变化的规律.为此,从实验和理论两方面着手.在实验方面,对某人用快速静脉注射方式一次注入该药物300mg后,在一定时刻t(小时)采集血样,测得血药浓度c.如表:血药浓度c(t)的测试数据

t0.250.511.523468c19.2118.1515.3614.1012.899.327.455.243.01问题：1.在快速静脉注射的给药方式下，研究血药浓度（单位体积血液中的药物含量）的变化规律；2.给定药物的最小有效浓度和最大治疗浓度，设计给药方案：每次注射剂量多大；间隔时间多长？

配药方案第十页，共六十五页，2022年，8月28日二、问题的解决（1）插值法；（2）曲线拟合法．

1、问题的抽象xx1x2…xmyy1y2…ym构造一个简单易于计算的近似函数

p(x)f(x)

（精确函数）。2、构造近似函数，p(x)的方法有两种：在实验中经常给出一组离散点，第十一页，共六十五页，2022年，8月28日插值法定义：当精确函数y=f(x)非常复杂或未知时，在一系列节点

x0…xn

处测得函数值y0

=f(x0),…，yn

=f(xn)，由此构造一个简单易算的近似函数p(x)

f(x)，满足条件p(xi)=f(xi)(i=0,…n)，（插值条件）这里的p(x)

称为f(x)的插值函数；构造插值函数的方法为插值法。第十二页，共六十五页，2022年，8月28日曲线拟合

但是不要求使p(xi)=yi,而只要p(xi)yi

总体上尽可能小。这种构造近似函数p(x)

的方法称为曲线拟合法，p(x)称为拟合函数。定义：当精确函数y=f(x)非常复杂或未知时，在一系列节点x0…xn

处,测得函数值y0

,…，yn

，由此构造一个简单易算的近似函数p(x)

f(x)，第十三页，共六十五页，2022年，8月28日插值与拟合的相同点都需要根据已知数据构造函数。可使用得到函数计算未知点的函数值。xx1x2…xmyy1y2…ym求一个简单易算的近似函数

p(x)f(x)

。第十四页，共六十五页，2022年，8月28日插值与拟合的不同点插值:过节点;；拟合:不过点,整体近似；第十五页，共六十五页，2022年，8月28日插值法拉格朗日插值牛顿插值三次埃尔米特插值法分段线性插值分段三次埃尔米特插值法三次样条插值第十六页，共六十五页，2022年，8月28日1、拉格朗日插值公式（１）定义对给定的n+1个节点x0,x1,x2,…,xn及对应的函数值y0,

y1,y2,…,yn，构造一个n次插值多项式：即为拉格朗日插值公式，其中插值基函数第十七页，共六十五页，2022年，8月28日拉格朗日插值的matlab实现functiony=lagrange(x0,y0,x)%x0插值节点，y0插值节点处的函数值，x要计算函数值的点；n=length(x0);%计算x0的长度m=length(x);%计算x的长度fori=1:ms=0;z=x(i);

fork=1:np=1.0;forj=1:n ifj~=kp=p*(z-x0(j))/(x0(k)-x0(j));%计算插值基函数

endends=p*y0(k)+s;endy(i)=s;%计算在x(i)处的函数值（拉格朗日）end第十八页，共六十五页，2022年，8月28日2、牛顿插值法牛顿插值公式：Nn(x)=f(x0)+f[x0,x1](x-x0)+…+f[x0,x1,x2,…,xn](x-x0)(x-x1)…(x-xn)其中：

f[x0,x1]一阶差商

f[x0,x1,x2,…,xn]n阶差商注：牛顿插值法与拉格朗日插值法，同一个多项式，不同的表达方式，但是计算量不一样，牛顿插值法的计算量小。第十九页，共六十五页，2022年，8月28日龙格现象Runge在上个世纪初发现：在[-5,5]上用n+1个等距节点作n次插值多项式Pn(x)，当在n→∞时，插值多项式Pn(x)在区间中部趋于f(x)=1/(1+x2),但对于3.63≤∣x∣≤1的x，Pn(x)严重发散。用图形分析问题。

第二十页，共六十五页，2022年，8月28日forn=10:2:20%从10等份到20等份x0=[-5:10/n:5];%插值节点y0=1./(1+x0.^2);%插值节点处的精确函数值x=[-5:0.1:5];%要进行计算函数值的点y=lagrange(x0,y0,x);%调用函数计算x点的函数值plot(x0,y0,‘*’,x,1./(1+x.^2),‘r’,x,y)%绘制图形pause%等待，按任意键end第二十一页，共六十五页，2022年，8月28日3、分段低次插值法（1）分段线性插值

定义：已知n+1个不同节点x0,x1,…,xn,构造分段多项式I(x)，使之满足l

I(x)在[a，b]上连续；l

I(xk)=yk；l

I(x)在[xi,xi+1]上是一次多项式；

I(x)=第二十二页，共六十五页，2022年，8月28日（2）分段三次埃尔米特插值法定义：已知n+1个不同节点x0,x1,…,xn,构造分段多项式I(x)，使之满足：l

I(x)在[a，b]上二阶连续导数；l

I(xk)=yk，I’(xk)=y’k，；l

I(x)在[xi,xi+1]上是三次次多项式。第二十三页，共六十五页，2022年，8月28日4、三次样条插值法

对于给定n+1个不同节点x0,x1,…,xn及函数值y0,y1,…,yn，其中a=x0<x1<…<xn=b，构造三次样条插值函数S(x)。

S(x)称为三次样条函数时需满足：l

S(x)在[a,b]上二阶导数连续；l

S(xk)=yk(k=0,1,…,n)；l

每个子区间[xk,xk+1]上S(x)是三次多项式(k=0,1,…,n)。第二十四页，共六十五页，2022年，8月28日插值法的matlab实现—一维插值

命令：interp1(x0,y0,x，’method’)

其中：x0：插值节点；

y0：插值节点处的函数值；

x：要计算函数值的点；

method：

linear

：分段线性插值；

cubic：分段三次埃尔米特插值；

spline：三次样条插值。第二十五页，共六十五页，2022年，8月28日插值法的应用一水库上游河段降暴雨，根据预报测算上游流入水库的流量为Q(t)(102立方米/秒)

：

t(时)81216243044485660Q（t）3654789210135251613

通过这个预报值，分别用不同的数值方法插值法来估计14和20时上游流入水库的流量。第二十六页，共六十五页，2022年，8月28日二维插值的MATLAB实现

在MATLAB中，二维插值命令常用的有两个，

1、一个是网格节点插值：

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

其中，

z：被插值点处的函数值；

x0,y0,z0：插值节点，x0,y0为向量，z0是矩阵，其列数等于x0的长度，行数等于y0的长度；

x,y：要计算函数值的点；

interp1(x0,y0,x，’method’)第二十七页，共六十五页，2022年，8月28日山体地貌要在某山区方圆大约27平方公里范围内修建一条公路，从山脚出发经过一个居民区，再到达一个矿区。横向纵向分别每隔400米测量一次，得到一些地点的高程：试做出该山区的地貌图，并对几种插值法进行比较.第二十八页，共六十五页，2022年，8月28日程序设计：clearx0=[1200:400:4000];y0=[1200:400:3600];z0=[11301250128012301040900500700;13201450142014001300700900850;139015001500140090011001060950;15001200110013501450120011501010;15001200110015501600155011801070;15001550160015501600160016001550;1480150015501510143013001200980];xi=1200:10:4000;%加密数据点yi=1200:10:3600;zil=interp2(x0,y0,z0,xi',yi,'linear');%线性插值zic=interp2(x0,y0,z0,xi',yi,'cubic');%三次插值zis=interp2(x0,y0,z0,xi',yi,'spline');%样条插值subplot(2,2,1)mesh(x0,y0,z0)subplot(2,2,2)mesh(xi,yi,zil)subplot(2,2,3)mesh(xi,yi,zic)subplot(2,2,4)mesh(xi,yi,zis)第二十九页，共六十五页，2022年，8月28日第三十页，共六十五页，2022年，8月28日二维插值的MATLAB实现2、另一个是离散数据节点的插值命令：

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

其中，

z：被插值点处的函数值；

x0,y0,z0：插值节点，x0,y0,z0均为向量；

x,y：被插值点；

method：插值方法，包括：

'linear'——线性插值；

'cubic'——三次插值；第三十一页，共六十五页，2022年，8月28日船在该海域会搁浅吗？

在某海域测得一些点(x,y)处的水深z由下表给出，船的吃水深度为5英尺，在矩形区域（75，200）*（-50，150）里的哪些地方船要避免进入.第三十二页，共六十五页，2022年，8月28日解决问题的步骤：1.作出测量点的分布图；2.求出矩形区域（75，200）*（-50，150）的细分网格节点之横、纵坐标向量；3.利用MATLAB中的散点插值函数求网格节点的水深；4.作出海底曲面图形和等高线图；5.作出水深小于5的海域范围.第三十三页，共六十五页，2022年，8月28日程序clearx0=[129140103.588185.5195105157.5107.57781162162117.5];y0=[7.5141.52314722.5137.585.5-6.5-81356.5-66.584-33.5];z0=[-4-8-6-8-6-8-8-9-9-8-8-9-4-9];subplot(2,2,1)plot(x0,y0,'+');%作出测量点的分布图；x=75:1:200;%加密y=-50:1:150;[x,y]=meshgrid(x,y);z=griddata(x0,y0,z0,x,y,'cubic');subplot(2,2,2)mesh(x,y,z),%用插值方法求出网格节点处的z坐标矩阵,绘制出三维图形subplot(2,2,3)meshc(x,y,z),%绘制等高线subplot(2,2,4)contour(x,y,z,[-5-5]);%水深5英尺处海底曲面的等高线gridon第三十四页，共六十五页，2022年，8月28日第三十五页，共六十五页，2022年，8月28日拟合的标准（1）用各点误差绝对值的和表示（2）用各点误差按绝对值的最大值表示（3）用各点误差的平方和表示第三十六页，共六十五页，2022年，8月28日最小二乘拟合式中R2称为均方误差。由于计算均方误差的最小值的原则容易实现而被广泛采用。按均方误差达到极小构造拟合曲线的方法称为最小二乘法。第三十七页，共六十五页，2022年，8月28日++++++++++++++++++++++++++++++p=a1+a2xp=a1+a2x+a3x2p=a1+a2x+a3x2p=a1+a2/xp=aebxp=ae-bx将数据(xi,yi)i=1,…,n

作图，通过直观判断确定p(x)：第三十八页，共六十五页，2022年，8月28日MATLAB---曲线拟合工具箱Matlab有一个功能强大的曲线拟合工具箱（CurveFittingToolbox

）cftool，使用方便，能实现多种类型的线性、非线性曲线拟合。调用：cftool界面如下所示第三十九页，共六十五页，2022年，8月28日第四十页，共六十五页，2022年，8月28日“Data”按钮数据的选取点击“Data”按钮，弹出“Data”窗口；利用Xdata和Ydata的下拉菜单读入数据x,y，可修改数据集名“Datasetname”，然后点击“Createdataset”按钮，退出“Data”窗口，返回工具箱界面，这时会自动画出数据集的曲线图；

第四十一页，共六十五页，2022年，8月28日第四十二页，共六十五页，2022年，8月28日“Fitting”按钮曲线拟合点击“Fitting”按钮，弹出“Fitting”窗口；点击“Newfit”按钮，可修改拟合项目名称“Fitname”，通过“Dataset”下拉菜单选择数据集，然后通过下拉菜单“Typeoffit”选择拟合曲线的类型。第四十三页，共六十五页，2022年，8月28日SSEThesumofsquaresduetoerror.Thisstatisticmeasuresthedeviationoftheresponsesfromthefittedvaluesoftheresponses.Avaluecloserto0indicatesabetterfit.偏差平方和，越接近0越好第四十四页，共六十五页，2022年，8月28日R-square

Thecoefficientofmultipledetermination.Thisstatisticmeasureshowsuccessfulthefitisinexplainingthevariationofthedata.Avaluecloserto1indicatesabetterfit.复相关系数平方(决定系数)，越接近1越好第四十五页，共六十五页，2022年，8月28日AdjustedR-squareThedegreeoffreedomadjustedR-square.Avaluecloserto1indicatesabetterfit.Itisgenerallythebestindicatorofthefitqualitywhenyouaddadditionalcoefficientstoyourmodel.修正的复相关系数平方，越接近1越好第四十六页，共六十五页，2022年，8月28日AdjustedR-square下列公式中的m为拟合函数中待估参数个数，如：对一元一次多项式拟合，f(x)=a+bx，此时m=2，n为数据点个数。该修正类似修正的样本方差使其为总体方差的无偏估计。第四十七页，共六十五页，2022年，8月28日RMSETherootmeansquarederror.Avaluecloserto0indicatesabetterfit.偏差平方的均值的算术平方根，越接近0越好第四十八页，共六十五页，2022年，8月28日曲线拟合好坏如何评价首要指标是目标函数误差最小（拟合度最大）；其次是应考虑关键点的吻合，这些关键点包括：初始点（有时是原点）、拐点、峰值点、极值点、中间点、渐近点、终值点等，在这些关键点上，数据观察值点与函数值点应尽可能一致；再次是拟合的模型应尽可能简单（模型的形式简单，参数数少）。第四十九页，共六十五页，2022年，8月28日实践中如何选择模型？在数据拟合实践中，理性模型毕竟是少数，大多数的情形是根据数据的趋势寻找合适的模型，有时好几个模型对数据都有较好的拟合，但通过对关键点的比较总会找到一种最合适的模型。在选择不同的模型时，合理性和可解释性是首要考虑的因素。第五十页，共六十五页，2022年，8月28日美国人口问题据美国人口普查局数据：从1790每隔10年至2000年的总人口（单位：百万）如下示

t=1790:10:2000;p=[3.9,5.3,7.2,9.6,12.9,17.1, 23.1,31.4,38.6,50.2,62.9,76,92, 105.7,122.8,131.7,150.7,179,205,226.5,251.4,281.422];第五十一页，共六十五页，2022年，8月28日美国人口问题t=1790:10:2000;p=[3.9,5.3,7.2,9.6,12.9,17.1, 23.1,31.4,38.6,50.2,62.9,76,92, 105.7,122.8,131.7,150.7,179,205,226.5,251.4,281.422];预测2001，2002年的美国人口数？并与调查数据285.318，288.369比较，选择拟合较好的模型。第五十二页，共六十五页，2022年，8月28日Matlab求解在命令窗口输入命令cftool回车，得拟合的图形用户界面第五十三页，共六十五页，2022年，8月28日第五十四页，共六十五页，2022年，8月28日第五十五页，共六十五页，2022年，8月28日结果分析LinearmodelPoly2:f(x)=p1*x^2+p2*x+p3Coefficients(with95%confidencebounds):p1=0.006757(0.006369,0.007144)p2=-24.32(-25.78,-22.85)p3=2.188e+004(2.049e+004,2.327e+004)Goodnessoffit:SSE: