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

1、 插值与曲线拟合插值与曲线拟合-5-4-3-2-1012345-0.500.511.52一、问题的提出 在生产和实验中，关于函数f(x)，经常存在两种情况：（1）其表达式不便于计算；（2）无表达式. 而只有函数在给定点的函数值，怎样预测其它点的函数值？xx0 x1x2xnyy0y1y2yn飞机机翼制造 下表给出的下表给出的x x、y y数据位于机翼端面的轮廓线上，数据位于机翼端面的轮廓线上，Y1Y1和和Y2Y2分分别对应轮廓的上下线。假设需要得到别对应轮廓的上下线。假设需要得到x x坐标每改变坐标每改变0.10.1时的时的y y坐标，试完成加工所需数据，画出曲线坐标，试完成加工所需数据，画出曲

2、线. .x035791112131415Y101.82.22.73.03.12.92.52.01.6Y201.21.72.02.02.01.81.21.01.6山体地貌n要在某山区方圆大约27平方公里范围内修建一条公路，从山脚出发经过一个居民区，再到达一个矿区。横向纵向分别每隔400米测量一次，得到一些地点的高程： 山区地貌：山区地貌： 在某山区测得一些地点的高程如下表。平面区域为在某山区测得一些地点的高程如下表。平面区域为 1200=x=4000,1200=y=3600) 试作出该山区的地貌图和等高线图，并对几种插值方法进行比较。试作出该山区的地貌图和等高线图，并对几种插值方法进行比较。 X

3、 Y 1200 1600 2000 2400 2800 3200 3600 4000 1200 1130 1250 1280 1230 1040 900 500 700 1600 1320 1450 1420 1400 1300 700 900 850 2000 1390 1500 1500 1400 900 1100 1060 950 2400 1500 1200 1100 1350 1450 1200 1150 1010 2800 1500 1200 1100 1550 1600 1550 1380 1070 3200 1500 1550 1600 1550 1600 1600 1600

4、1550 3600 1480 1500 1550 1510 1430 1300 1200 980 试做出该山区的地貌图.船在该海域会搁浅吗？船在该海域会搁浅吗？-作业作业 在某海域测得一些点(x,y)处的水深z由下表给出，船的吃水深度为5英尺，在矩形区域（75，200）*（-50，150）里的哪些地方船要避免进入.xyz129 140 103.5 88 185.5 195 1057.5 141.5 23 147 22.5 137.5 85.54 8 6 8 6 8 8xyz157.5 107.5 77 81 162 162 117.5-6.5 -81 3 56.5 -66.5 84 -33.5

5、9 9 8 8 9 4 9水深和流速的问题 在水文数据测量中，不同水深的流速是不同的. 水文数据的测量时天天进行的，为了减少测量的工作，希望得到确定的水深和水流之间的关系. 为此测量了一系列不同水深和流速值. 下表给出了对某河流的测量数据，其中水深和流速根据适当的单位进行了规范化，共10个值.2.983.063.133.203.193.233.253.263.223.19流速0.90.80.70.60.50.40.30.20.10水深美国人口问题n据美国人口普查局数据： 从1790每隔10年至2000年的总人口（单位：百万）如下示n t = 1790:10:2000; np = 3.9, 5.

6、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; n预测2001，2002年的美国人口数？并与调查数据285.318，288.369比较，选择拟合较好的模型。农作物施肥效果分析1992年A题 在农业生产试验研究中，对某地区土豆的产量与化肥的关系做了一实验，得到了氮肥、磷肥的施肥量与土豆产量的对应关系如下表： 1.根据上表数据分别给出土豆产量与氮、磷肥的关系式。

7、 2.施肥问题优化策略氮肥量（公斤/公顷）03467101135202259336404471土豆产量（公斤）15.1821.3625.7232.293439.4543.1543.4640.8330.75磷肥量（公斤/公顷）024497398147196245294342土豆产量（公斤）33.4632.4736.0637.964140.141。342.240.442.7配药方案-作业 一种新药用于临床之前, 必须设计给药方案. 在快速静脉注射的给药方式下, 所谓给药方案是指, 每次注射剂量多大, 间隔时间多长. 药物进入机体后随血液输送到全身, 在这个过程中不断地被吸收, 分布, 代谢, 最终

8、排出体外. 药物在血液中的浓度, 即单位体积血液中的药物含量, 称血药浓度. 在最简单的一室模型中, 将整个机体看作一个房室, 称中心室, 室内的血药浓度是均匀的. 快速静脉注射后, 浓度立即上升; 然后逐渐下降. 当浓度太低时, 达不到预期的治疗效果; 血药浓度太高, 又可能导致药物中毒或副作用太强. 临床上, 每种药物有一个最小有效浓度 c1 和一个最大治疗浓度 c2. 设计给药方案时, 要使血药浓度保持在 c1-c2 之间. 设本题所研究药物的最小有效浓度c1=10, 最大治疗浓度 c2=25( ). /g ml 显然, 要设计给药方案, 必须知道给药后血药浓度随时间变化的规律. 为此,

9、 从实验和理论两方面着手. 在实验方面, 对某人用快速静脉注射方式一次注入该药物300mg后, 在一定时刻 t (小时)采集血样, 测得血药浓度c. 如表: 血药浓度c(t) 的测试数据 t0.250.511.523468c19.2118.1515.3614.1012.899.327.455.243.01问题问题 ：1. 1. 在快速静脉注射的给药方式下，研究血药浓度（单位体积血液中的药物含在快速静脉注射的给药方式下，研究血药浓度（单位体积血液中的药物含量）的变化规律；量）的变化规律；2. 2. 给定药物的最小有效浓度和最大治疗浓度，设计给药方案：每次注射剂量给定药物的最小有效浓度和最大治疗浓

10、度，设计给药方案：每次注射剂量多大；间隔时间多长？多大；间隔时间多长？配药方案二、问题的解决（1）插值法；（2）曲线拟合法 1、问题的抽象、问题的抽象xx1x2xmyy1y2ym构造一个构造一个简单易于计算简单易于计算的近似函数的近似函数 p(x) f (x) （精确函数）。（精确函数）。2、构造近似函数，、构造近似函数， p(x) 的方法有两种：的方法有两种：在实验中经常给出一组离散点，在实验中经常给出一组离散点，插值法定义：定义：当精确函数当精确函数 y = f(x) 非常复杂非常复杂或或未知时未知时，在一系列节点，在一系列节点 x0 xn 处测得函数值处测得函数值 y0 = f(x0),

11、 ，yn = f(xn)，由此构造一个由此构造一个简单易算简单易算的近似函数的近似函数 p(x) f(x)，满足条件，满足条件p(xi) = f(xi) (i = 0, n)，（插值条件），（插值条件）这里的这里的 p(x) 称为称为f(x) 的的插值函数；插值函数； 构造插值函数的方法为构造插值函数的方法为插值法插值法。曲线拟合 但是不要求使但是不要求使 p(xi) = yi ,而只要而只要 p(xi) yi 总体上总体上尽可能小。这种构尽可能小。这种构造近似函数造近似函数p(x) 的方法称为的方法称为曲线拟合法曲线拟合法， p(x) 称为称为拟合函数。拟合函数。定义：定义： 当精确函数当精

12、确函数 y = f(x) 非常复杂非常复杂或或未知时未知时，在一系列节点，在一系列节点x0 xn 处处,测得函数值测得函数值 y0 , ，yn ，由此构造一个，由此构造一个简单易简单易 算算的近似函的近似函数数 p(x) f(x)，插值与拟合的相同点n都需要根据已知数据构造函数。n可使用得到函数计算未知点的函数值。xx1x2xmyy1y2ym求一个求一个简单易算简单易算的近似函数的近似函数 p(x) f (x) 。插值与拟合的不同点n插值: 过节点; ；n拟合: 不过点, 整体近似；插值法n拉格朗日插值 n牛顿插值 n三次埃尔米特插值法n分段线性插值n分段三次埃尔米特插值法n三次样条插值1、

13、拉格朗日插值公式（）定义对给定的n+1个节点x0 , x1,x2,xn及对应的函数值y0 , y1,y2,yn，构造一个n次插值多项式：即为拉格朗日插值公式拉格朗日插值公式，其中)(0 xlyyknkk)()()()()()()(110110nkkkkkknkkkxxxxxxxxxxxxxxxxxl插值基函插值基函数数拉格朗日插值的matlab实现function y=lagrange(x0,y0,x) % x0插值节点， y0插值节点处的函数值，x要计算函数值的点；n=length(x0); %计算x0的长度m=length(x); %计算x的长度for i=1:m s=0;z=x(i);

14、nfor k=1:nn p=1.0;n for j=1:nn if j=knp=p*(z-x0(j)/(x0(k)-x0(j); n%计算插值基函数n endn endn s=p*y0(k)+s;nendny(i)=s; %计算在x(i)处的函数值（拉格朗日）nend2、牛顿插值法牛顿插值公式：Nn(x)=f(x0)+fx0,x1(x-x0)+fx0,x1,x2,xn(x-x0)(x-x1)(x-xn)其中： fx0,x1 一阶差商 fx0,x1,x2,xn n阶差商n注：牛顿插值法与拉格朗日插值法，同一个多项式，不同的表达方式，但是计算量不一样，牛顿插值法的计算量小。龙格现象 Runge在上

15、个世纪初发现： 在-5,5上用n+1个等距节点作n次插值多项式Pn(x)，当在n时，插值多项式Pn(x)在区间中部趋于f(x)=1/(1+x2) ,但对于3.63x1的x，Pn(x)严重发散。用图形分析问题。 for n=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 %等待，按任意键end3、分段低次插值法（

16、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)=,11111kkkkkkkkkkxxxxxxxyxxxxy（2）分段三次埃尔米特插值法定义：已知n+1个不同节点x0,x1,xn ,构造分段多项式I(x)，使之满足：l I(x)在a，b上二阶连续导数；l I(xk)=yk， I(xk)=yk， ；l I(x)在xi,xi+1上是三次次多项式。4、三次样条插值法 对于给定n+1个不同节点x0,x1,xn及函数值y0,y1,yn， 其中a=x0

17、x1xn=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)。插值法的matlab实现一维插值 命令：interp1(x0,y0,x，method) 其中：x0：插值节点； y0：插值节点处的函数值； x：要计算函数值的点； method： l i n e a r ：分段线性插值；：分段线性插值； c u b i c ：分段三次埃尔米特插值； s p l i n e ：三次样条插值。插值法的应用n一水库上游河段降暴雨，根据预报测算

18、上游流入水库的流量为一水库上游河段降暴雨，根据预报测算上游流入水库的流量为Q(tQ(t) (10) (102 2立方米立方米/ /秒秒) ) ： t (时时) 8 12 16 24 30 44 48 56 60 Q（t） 36 54 78 92 101 35 25 16 13 通过这个预报值，分别用不同的数值方法插值法来估计通过这个预报值，分别用不同的数值方法插值法来估计14 14 和和2020时上游流入水库的流量。时上游流入水库的流量。二维插值的MATLAB实现 在MATLAB中，二维插值命令常用的有两个， 1、一个是网格节点插值：、一个是网格节点插值： z=interp2(x0,y0,z0

19、,x,y,method) 其中，n z：被插值点处的函数值；n x0,y0,z0：插值节点， x0,y0为向量，z0是矩阵，其列数等于x0的长度，行数等于y0的长度；n x,y： 要计算函数值的点； interp1(x0,y0,x，method)山体地貌n要在某山区方圆大约27平方公里范围内修建一条公路，从山脚出发经过一个居民区，再到达一个矿区。横向纵向分别每隔400米测量一次，得到一些地点的高程： 山区地貌：山区地貌： 在某山区测得一些地点的高程如下表。平面区域为在某山区测得一些地点的高程如下表。平面区域为 1200=x=4000,1200=y=3600) 试作出该山区的地貌图和等高线图，并

20、对几种插值方法进行比较。试作出该山区的地貌图和等高线图，并对几种插值方法进行比较。 X Y 1200 1600 2000 2400 2800 3200 3600 4000 1200 1130 1250 1280 1230 1040 900 500 700 1600 1320 1450 1420 1400 1300 700 900 850 2000 1390 1500 1500 1400 900 1100 1060 950 2400 1500 1200 1100 1350 1450 1200 1150 1010 2800 1500 1200 1100 1550 1600 1550 1380 10

21、70 3200 1500 1550 1600 1550 1600 1600 1600 1550 3600 1480 1500 1550 1510 1430 1300 1200 980 试做出该山区的地貌图，并对几种插值法进行比较.n程序设计：程序设计：nclearnx0=1200:400:4000; ny0=1200:400:3600; nz0=1130 1250 1280 1230 1040 900 500 700;n 1320 1450 1420 1400 1300 700 900 850;n 1390 1500 1500 1400 900 1100 1060 950;n 1500 120

22、0 1100 1350 1450 1200 1150 1010;n 1500 1200 1100 1550 1600 1550 1180 1070;n 1500 1550 1600 1550 1600 1600 1600 1550;n 1480 1500 1550 1510 1430 1300 1200 980;nxi=1200:10:4000; %加密数据点nyi=1200:10:3600; nzil=interp2(x0,y0,z0,xi,yi,linear); %线性插值nzic=interp2(x0,y0,z0,xi,yi,cubic); %三次插值nzis=interp2(x0,y0

23、,z0,xi,yi,spline); %样条插值nsubplot(2,2,1) nmesh(x0,y0,z0)nsubplot(2,2,2) nmesh(xi,yi,zil)nsubplot(2,2,3) nmesh(xi,yi,zic)nsubplot(2,2,4) nmesh(xi,yi,zis)二维插值的MATLAB实现2、另一个是离散数据节点的插值命令： z=griddata(x0,y0,z0,x,y,method) 其中，n z：被插值点处的函数值；n x0,y0,z0：插值节点， x0,y0,z0均为向量；n x,y：被插值点； method：插值方法，包括： n linear线性

24、插值；n cubic 三次插值；船在该海域会搁浅吗？船在该海域会搁浅吗？ 在某海域测得一些点(x,y)处的水深z由下表给出，船的吃水深度为5英尺，在矩形区域（75，200）*（-50，150）里的哪些地方船要避免进入.xyz129 140 103.5 88 185.5 195 1057.5 141.5 23 147 22.5 137.5 85.54 8 6 8 6 8 8xyz157.5 107.5 77 81 162 162 117.5-6.5 -81 3 56.5 -66.5 84 -33.59 9 8 8 9 4 9解决问题的步骤：1.作出测量点的分布图；2.求出矩形区域（75，200）

25、*（-50，150）的细分网格节点之横、纵坐标向量；3.利用MATLAB中的散点插值函数求网格节点的水深；4.作出海底曲面图形和等高线图；5.作出水深小于5的海域范围.程序nclearnx0=129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162 162 117.5; ny0=7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5; nz0=-4 -8 -6 -8 -6 -8 -8 -9 -9 -8 -8 -9 -4 -9; nsubplot(2,2,1)nplot(x0,

26、y0,+); %作出测量点的分布图；作出测量点的分布图； nx=75:1:200; %加密加密ny=-50:1:150;nx,y=meshgrid(x,y); nz=griddata(x0,y0,z0,x,y,cubic);nsubplot(2,2,2)nmesh(x,y,z), %用插值方法求出网格节点处的用插值方法求出网格节点处的z坐标矩阵坐标矩阵,绘制出三绘制出三维图形维图形nsubplot(2,2,3)nmeshc(x,y,z), %绘制等高线绘制等高线nsubplot(2,2,4)ncontour(x,y,z,-5 -5); %水深水深5英尺处海底曲面的等高线英尺处海底曲面的等高线n

27、grid on拟合的标准n（1）用各点误差绝对值的和表示n（2）用各点误差按绝对值的最大值表示n（3）用各点误差的平方和表示 11()miiiRp xy1( )maxiii mRp xy 221()miiiRp xy最小二乘拟合n式中 R2 称为均方误差。由于计算均方误差的最小值的原则容易实现而被广泛采用。n按均方误差达到极小构造拟合曲线的方法称为最小二乘法。221( ()miiiRp xy+p=a1+a2xp=a1+a2x+a3x2p=a1+a2x+a3x2p=a1+a2/xp=aebxp=ae-bx将数据将数据 (xi,yi) i=1, ,n 作图，通过直观判断确定作图，通过直观判断确定

28、p(x)：MATLAB-曲线拟合工具箱nMatlab有一个功能强大的曲线拟合工具箱 （Curve Fitting Toolbox ） cftool， 使用方便，能实现多种类型的线性、非线性曲线拟合。n调用：cftooln界面如下所示“Data”按钮数据的选取n点击“Data”按钮，弹出“Data”窗口；n利用X data和Y data的下拉菜单读入数据x,y，可修改数据集名“Data set name”，然后点击“Create data set”按钮，退出“Data”窗口，返回工具箱界面，这时会自动画出数据集的曲线图；“Fitting”按钮 曲线拟合n点击“Fitting”按钮，弹出“Fitt

29、ing”窗口；n点击“New fit”按钮，可修改拟合项目名称“Fit name”，通过“Data set”下拉菜单选择数据集，然后通过下拉菜单“Type of fit”选择拟合曲线的类型。SSEnThe sum of squares due to error. This statistic measures the deviation of the responses from the fitted values of the responses. A value closer to 0 indicates a better fit.n偏差平方和，越接近0越好niiyxpSSE1i2)(R-

30、squaren The coefficient of multiple determination. This statistic measures how successful the fit is in explaining the variation of the data. A value closer to 1 indicates a better fit.n复相关系数平方(决定系数)，越接近1越好Adjusted R-squaren The degree of freedom adjusted R-square. A value closer to 1 indicates a be

31、tter fit. It is generally the best indicator of the fit quality when you add additional coefficients to your model.n修正的复相关系数平方，越接近1越好Adjusted R-squaren下列公式中的m为拟合函数中待估参数个数，如：对一元一次多项式拟合，f(x) = a + bx，此时m=2，n为数据点个数。该修正类似修正的样本方差使其为总体方差的无偏估计。RMSEnThe root mean squared error. A value closer to 0 indicates

32、 a better fit.n偏差平方的均值的算术平方根，越接近0越好曲线拟合好坏如何评价n首要指标是目标函数误差最小（拟合度最大）；n其次是应考虑关键点的吻合，这些关键点包括：初始点（有时是原点）、拐点、峰值点、极值点、中间点、渐近点、终值点等，在这些关键点上，数据观察值点与函数值点应尽可能一致；n再次是拟合的模型应尽可能简单（模型的形式简单，参数数少）。实践中如何选择模型？n在数据拟合实践中，理性模型毕竟是少数，大多数的情形是根据数据的趋势寻找合适的模型，有时好几个模型对数据都有较好的拟合，但通过对关键点的比较总会找到一种最合适的模型。n在选择不同的模型时，合理性和可解释性是首要考虑的因素

33、。美国人口问题n据美国人口普查局数据： 从1790每隔10年至2000年的总人口（单位：百万）如下示n t = 1790:10:2000; np = 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; 美国人口问题nt = 1790:10:2000; np = 3.9, 5.3 , 7.2 , 9.6 , 12.9 , 17.1 , 23.1 ,

34、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; n预测2001，2002年的美国人口数？并与调查数据285.318，288.369比较，选择拟合较好的模型。Matlab求解n在命令窗口输入命令ncftooln回车，得拟合的图形用户界面结果分析nLinear model Poly2: f(x) = p1*x2 + p2*x + p3nCoefficients (with 95% confidence bounds):np1 = 0.006757 (0.006369, 0.007144)np2 = -24.32 (-25.78, -22.85)np3 = 2.188e+004 (2.049e+004, 2.327e+004)nGoodness of fit:n SSE: 184.5n R-square: