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

1、关于数学建模插值法与曲线拟合讲课第一张，PPT共六十五页，创作于2022年6月一、问题的提出 在生产和实验中，关于函数f(x)，经常存在两种情况：（1）其表达式不便于计算；（2）无表达式. 而只有函数在给定点的函数值，怎样预测其它点的函数值？xx0 x1x2xnyy0y1y2yn第二张，PPT共六十五页，创作于2022年6月飞机机翼制造 下表给出的x、y数据位于机翼端面的轮廓线上，Y1和Y2分别对应轮廓的上下线。假设需要得到x坐标每改变0.1时的y坐标，试完成加工所需数据，画出曲线.x035791112131415Y101.82.22.73.03.12.92.52.01.6Y201.21.72

2、.02.02.01.81.21.01.6第三张，PPT共六十五页，创作于2022年6月山体地貌要在某山区方圆大约27平方公里范围内修建一条公路，从山脚出发经过一个居民区，再到达一个矿区。横向纵向分别每隔400米测量一次，得到一些地点的高程：试做出该山区的地貌图.第四张，PPT共六十五页，创作于2022年6月船在该海域会搁浅吗？-作业 在某海域测得一些点(x,y)处的水深z由下表给出，船的吃水深度为5英尺，在矩形区域（75，200）*（-50，150）里的哪些地方船要避免进入.第五张，PPT共六十五页，创作于2022年6月水深和流速的问题 在水文数据测量中，不同水深的流速是不同的. 水文数据的测

3、量时天天进行的，为了减少测量的工作，希望得到确定的水深和水流之间的关系. 为此测量了一系列不同水深和流速值. 下表给出了对某河流的测量数据，其中水深和流速根据适当的单位进行了规范化，共10个值.第六张，PPT共六十五页，创作于2022年6月美国人口问题据美国人口普查局数据： 从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

4、.7 , 179 , 205 , 226.5 , 251.4 , 281.422; 预测2001，2002年的美国人口数？并与调查数据285.318，288.369比较，选择拟合较好的模型。第七张，PPT共六十五页，创作于2022年6月农作物施肥效果分析1992年A题 在农业生产试验研究中，对某地区土豆的产量与化肥的关系做了一实验，得到了氮肥、磷肥的施肥量与土豆产量的对应关系如下表： 1.根据上表数据分别给出土豆产量与氮、磷肥的关系式。 2.施肥问题优化策略氮肥量（公斤/公顷）03467101135202259336404471土豆产量（公斤）15.1821.3625.7232.293439.

5、4543.1543.4640.8330.75磷肥量（公斤/公顷）024497398147196245294342土豆产量（公斤）33.4632.4736.0637.964140.141。342.240.442.7第八张，PPT共六十五页，创作于2022年6月配药方案-作业 一种新药用于临床之前, 必须设计给药方案. 在快速静脉注射的给药方式下, 所谓给药方案是指, 每次注射剂量多大, 间隔时间多长. 药物进入机体后随血液输送到全身, 在这个过程中不断地被吸收, 分布, 代谢, 最终排出体外. 药物在血液中的浓度, 即单位体积血液中的药物含量, 称血药浓度. 在最简单的一室模型中, 将整个机体看

6、作一个房室, 称中心室, 室内的血药浓度是均匀的. 快速静脉注射后, 浓度立即上升; 然后逐渐下降. 当浓度太低时, 达不到预期的治疗效果; 血药浓度太高, 又可能导致药物中毒或副作用太强. 临床上, 每种药物有一个最小有效浓度 c1 和一个最大治疗浓度 c2. 设计给药方案时, 要使血药浓度保持在 c1-c2 之间. 设本题所研究药物的最小有效浓度c1=10, 最大治疗浓度 c2=25( ). 第九张，PPT共六十五页，创作于2022年6月 显然, 要设计给药方案, 必须知道给药后血药浓度随时间变化的规律. 为此, 从实验和理论两方面着手. 在实验方面, 对某人用快速静脉注射方式一次注入该药

7、物300mg后, 在一定时刻 t (小时)采集血样, 测得血药浓度c. 如表: 血药浓度c(t) 的测试数据 t0.250.511.523468c19.2118.1515.3614.1012.899.327.455.243.01问题 ：1. 在快速静脉注射的给药方式下，研究血药浓度（单位体积血液中的药物含量）的变化规律；2. 给定药物的最小有效浓度和最大治疗浓度，设计给药方案：每次注射剂量多大；间隔时间多长？配药方案第十张，PPT共六十五页，创作于2022年6月二、问题的解决（1）插值法；（2）曲线拟合法 1、问题的抽象xx1x2xmyy1y2ym构造一个简单易于计算的近似函数 p(x) f

8、(x) （精确函数）。2、构造近似函数， p(x) 的方法有两种：在实验中经常给出一组离散点，第十一张，PPT共六十五页，创作于2022年6月插值法定义：当精确函数 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) 的插值函数； 构造插值函数的方法为插值法。第十二张，PPT共六十五页，创作于2022年6月曲线拟合 但是不要求使 p(xi) = yi ,而只要 p(xi)

9、 yi 总体上尽可能小。这种构造近似函数p(x) 的方法称为曲线拟合法， p(x) 称为拟合函数。定义： 当精确函数 y = f(x) 非常复杂或未知时，在一系列节点x0 xn 处,测得函数值 y0 , ，yn ，由此构造一个简单易 算的近似函数 p(x) f(x)，第十三张，PPT共六十五页，创作于2022年6月插值与拟合的相同点都需要根据已知数据构造函数。可使用得到函数计算未知点的函数值。xx1x2xmyy1y2ym求一个简单易算的近似函数 p(x) f (x) 。第十四张，PPT共六十五页，创作于2022年6月插值与拟合的不同点插值: 过节点; ；拟合: 不过点, 整体近似；第十五张，P

10、PT共六十五页，创作于2022年6月插值法拉格朗日插值 牛顿插值 三次埃尔米特插值法分段线性插值分段三次埃尔米特插值法三次样条插值第十六张，PPT共六十五页，创作于2022年6月1、 拉格朗日插值公式（）定义对给定的n+1个节点x0 , x1,x2,xn及对应的函数值y0 , y1,y2,yn，构造一个n次插值多项式：即为拉格朗日插值公式，其中插值基函数第十七张，PPT共六十五页，创作于2022年6月拉格朗日插值的matlab实现function y=lagrange(x0,y0,x) % x0插值节点， y0插值节点处的函数值，x要计算函数值的点；n=length(x0); %计算x0的长度

11、m=length(x); %计算x的长度for i=1:m s=0;z=x(i); for k=1:n p=1.0; for j=1:n if j=kp=p*(z-x0(j)/(x0(k)-x0(j); %计算插值基函数 end end s=p*y0(k)+s;endy(i)=s; %计算在x(i)处的函数值（拉格朗日）end第十八张，PPT共六十五页，创作于2022年6月2、牛顿插值法牛顿插值公式：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阶差商注：牛顿插值法与拉格

12、朗日插值法，同一个多项式，不同的表达方式，但是计算量不一样，牛顿插值法的计算量小。第十九张，PPT共六十五页，创作于2022年6月龙格现象 Runge在上个世纪初发现： 在-5,5上用n+1个等距节点作n次插值多项式Pn(x)，当在n时，插值多项式Pn(x)在区间中部趋于f(x)=1/(1+x2) ,但对于3.63x1的x，Pn(x)严重发散。用图形分析问题。第二十张，PPT共六十五页，创作于2022年6月for n=10:2:20 %从10等份到20等份x0=-5:10/n:5; %插值节点y0=1./(1+x0.2); %插值节点处的精确函数值x=-5:0.1:5; %要进行计算函数值的点

13、y=lagrange(x0,y0,x); %调用函数计算x点的函数值 plot(x0,y0,*,x,1./(1+x.2),r,x,y) %绘制图形pause %等待，按任意键end第二十一张，PPT共六十五页，创作于2022年6月3、分段低次插值法（1）分段线性插值 定义： 已知n+1个不同节点x0,x1,xn ,构造分段多项式I(x)，使之满足lI(x)在a，b上连续；l I(xk)=yk；lI(x)在xi,xi+1上是一次多项式； I(x)=第二十二张，PPT共六十五页，创作于2022年6月（2）分段三次埃尔米特插值法定义：已知n+1个不同节点x0,x1,xn ,构造分段多项式I(x)，使

14、之满足：lI(x)在a，b上二阶连续导数；l I(xk)=yk， I(xk)=yk， ；lI(x)在xi,xi+1上是三次次多项式。第二十三张，PPT共六十五页，创作于2022年6月4、三次样条插值法 对于给定n+1个不同节点x0,x1,xn及函数值y0,y1,yn， 其中a=x0 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)。第二十四张，PPT共六十五页，创作于2022年6月插值法的matlab实现一维插值 命令

15、：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 ：三次样条插值。第二十五张，PPT共六十五页，创作于2022年6月插值法的应用一水库上游河段降暴雨，根据预报测算上游流入水库的流量为Q(t) (102立方米/秒) ： t (时) 8 12 16 24 30 44 48 56 60 Q（t） 36 54 78 92 101 35 25 16 13 通过这个预报值，分别用不同的数值方法插值法来估计

16、14 和20时上游流入水库的流量。第二十六张，PPT共六十五页，创作于2022年6月二维插值的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)第二十七张，PPT共六十五页，创作于2022年6月山体地貌要在某山区方圆大约27平方公里范围内修建一条公路，从山脚出发经过一个居民区，再到达

17、一个矿区。横向纵向分别每隔400米测量一次，得到一些地点的高程：试做出该山区的地貌图，并对几种插值法进行比较.第二十八张，PPT共六十五页，创作于2022年6月程序设计：clearx0=1200:400:4000; y0=1200:400:3600; z0=1130 1250 1280 1230 1040 900 500 700; 1320 1450 1420 1400 1300 700 900 850; 1390 1500 1500 1400 900 1100 1060 950; 1500 1200 1100 1350 1450 1200 1150 1010; 1500 1200 1100

18、1550 1600 1550 1180 1070; 1500 1550 1600 1550 1600 1600 1600 1550; 1480 1500 1550 1510 1430 1300 1200 980;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)sub

19、plot(2,2,2) mesh(xi,yi,zil)subplot(2,2,3) mesh(xi,yi,zic)subplot(2,2,4) mesh(xi,yi,zis)第二十九张，PPT共六十五页，创作于2022年6月第三十张，PPT共六十五页，创作于2022年6月二维插值的MATLAB实现2、另一个是离散数据节点的插值命令： z=griddata(x0,y0,z0,x,y,method) 其中， z：被插值点处的函数值； x0,y0,z0：插值节点， x0,y0,z0均为向量； x,y：被插值点； method：插值方法，包括： linear线性插值； cubic 三次插值；第三十一张

20、，PPT共六十五页，创作于2022年6月船在该海域会搁浅吗？ 在某海域测得一些点(x,y)处的水深z由下表给出，船的吃水深度为5英尺，在矩形区域（75，200）*（-50，150）里的哪些地方船要避免进入.第三十二张，PPT共六十五页，创作于2022年6月解决问题的步骤：1.作出测量点的分布图；2.求出矩形区域（75，200）*（-50，150）的细分网格节点之横、纵坐标向量；3.利用MATLAB中的散点插值函数求网格节点的水深；4.作出海底曲面图形和等高线图；5.作出水深小于5的海域范围.第三十三张，PPT共六十五页，创作于2022年6月程序clearx0=129 140 103.5 88

21、185.5 195 105 157.5 107.5 77 81 162 162 117.5; y0=7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -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,

22、y,z), %用插值方法求出网格节点处的z坐标矩阵,绘制出三维图形subplot(2,2,3)meshc(x,y,z), %绘制等高线subplot(2,2,4)contour(x,y,z,-5 -5); %水深5英尺处海底曲面的等高线grid on第三十四张，PPT共六十五页，创作于2022年6月第三十五张，PPT共六十五页，创作于2022年6月拟合的标准（1）用各点误差绝对值的和表示（2）用各点误差按绝对值的最大值表示（3）用各点误差的平方和表示 第三十六张，PPT共六十五页，创作于2022年6月最小二乘拟合式中 R2 称为均方误差。由于计算均方误差的最小值的原则容易实现而被广泛采用。按均

23、方误差达到极小构造拟合曲线的方法称为最小二乘法。第三十七张，PPT共六十五页，创作于2022年6月+p=a1+a2xp=a1+a2x+a3x2p=a1+a2x+a3x2p=a1+a2/xp=aebxp=ae-bx将数据 (xi,yi) i=1, ,n 作图，通过直观判断确定 p(x)：第三十八张，PPT共六十五页，创作于2022年6月MATLAB-曲线拟合工具箱Matlab有一个功能强大的曲线拟合工具箱 （Curve Fitting Toolbox ） cftool， 使用方便，能实现多种类型的线性、非线性曲线拟合。调用：cftool界面如下所示第三十九张，PPT共六十五页，创作于2022年6

24、月第四十张，PPT共六十五页，创作于2022年6月“Data”按钮数据的选取点击“Data”按钮，弹出“Data”窗口；利用X data和Y data的下拉菜单读入数据x,y，可修改数据集名“Data set name”，然后点击“Create data set”按钮，退出“Data”窗口，返回工具箱界面，这时会自动画出数据集的曲线图；第四十一张，PPT共六十五页，创作于2022年6月第四十二张，PPT共六十五页，创作于2022年6月“Fitting”按钮 曲线拟合点击“Fitting”按钮，弹出“Fitting”窗口；点击“New fit”按钮，可修改拟合项目名称“Fit name”，通过“

25、Data set”下拉菜单选择数据集，然后通过下拉菜单“Type of fit”选择拟合曲线的类型。第四十三张，PPT共六十五页，创作于2022年6月SSEThe 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.偏差平方和，越接近0越好第四十四张，PPT共六十五页，创作于2022年6月R-square The

26、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.复相关系数平方(决定系数)，越接近1越好第四十五张，PPT共六十五页，创作于2022年6月Adjusted R-square The degree of freedom adjusted R-square. A value closer to 1 ind

27、icates a better fit. It is generally the best indicator of the fit quality when you add additional coefficients to your model.修正的复相关系数平方，越接近1越好第四十六张，PPT共六十五页，创作于2022年6月Adjusted R-square下列公式中的m为拟合函数中待估参数个数，如：对一元一次多项式拟合，f(x) = a + bx，此时m=2，n为数据点个数。该修正类似修正的样本方差使其为总体方差的无偏估计。第四十七张，PPT共六十五页，创作于2022年6月RMSE

28、The root mean squared error. A value closer to 0 indicates a better fit.偏差平方的均值的算术平方根，越接近0越好第四十八张，PPT共六十五页，创作于2022年6月曲线拟合好坏如何评价首要指标是目标函数误差最小（拟合度最大）；其次是应考虑关键点的吻合，这些关键点包括：初始点（有时是原点）、拐点、峰值点、极值点、中间点、渐近点、终值点等，在这些关键点上，数据观察值点与函数值点应尽可能一致；再次是拟合的模型应尽可能简单（模型的形式简单，参数数少）。第四十九张，PPT共六十五页，创作于2022年6月实践中如何选择模型？在数据拟合实

29、践中，理性模型毕竟是少数，大多数的情形是根据数据的趋势寻找合适的模型，有时好几个模型对数据都有较好的拟合，但通过对关键点的比较总会找到一种最合适的模型。在选择不同的模型时，合理性和可解释性是首要考虑的因素。第五十张，PPT共六十五页，创作于2022年6月美国人口问题据美国人口普查局数据： 从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.

30、7 , 150.7 , 179 , 205 , 226.5 , 251.4 , 281.422; 第五十一张，PPT共六十五页，创作于2022年6月美国人口问题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比较，选择拟合

31、较好的模型。第五十二张，PPT共六十五页，创作于2022年6月Matlab求解在命令窗口输入命令cftool回车，得拟合的图形用户界面第五十三张，PPT共六十五页，创作于2022年6月第五十四张，PPT共六十五页，创作于2022年6月第五十五张，PPT共六十五页，创作于2022年6月结果分析Linear model Poly2: f(x) = p1*x2 + p2*x + p3Coefficients (with 95% confidence bounds):p1 = 0.006757 (0.006369, 0.007144)p2 = -24.32 (-25.78, -22.85)p3 = 2.188e+004 (2.049e+004, 2.327e+004)Goodness of fit: SSE: 184.5 R-square: 0.9989 Adjusted R-square: 0.9987 RMSE: 3.116第五十六张，PPT共六十五页，创作于2022年6月预测第五十七张，PPT共六十五页，创