数学建模讲稿插值,拟合,方程求根.ppt

插值与拟合 一、插值 在工程实践和科学实验中，常常需要从一组实验观测数据，揭 表示自变量x与因变量y之间的关系，通常可以采用两种方法：曲 线拟合和插值 插值在工程实践和科学实验中有着非常广泛而又十分重要的应 用，例如，信息技术中的图像重建、图像放大中为避免图像的 扭曲失真的插值补点、建筑工程的外观设计。化学工程实验数 据与模型的分析、天文观测数据、地理信息数据的处理如(天气 预报）以及社会经济现象的统计分析等等 1) 分段多项式插值 1. 插值方法 一维插值方法主要有分段多项式插值、三次样条插值等 2）分段三次埃尔米特插值 3）三次样条插值 上面介绍的分段线性插值，其总体光滑程度不够在数学上，光 滑程度的定量描述是：函数（曲线）的k阶导数存在且连续，则称 该曲线具有k阶光滑性自然，阶数越高光滑程度越好于是，分 段线性插值具有零阶光滑性，也就是不光滑；分段三次埃尔米特插 值具有一阶光滑性仅有这些光滑程度，在工程设计和机械加工等 实际中是不够的提高分段函数如多项式函数的次数，可望提高整 体曲线的光滑程度值得注意的是分段插值曲线的光滑性关键在于 段与段之间的衔接点（节点）处的光滑性 问：是否存在较低次多项式达到较高阶光滑性的方法？ 三次样条插值 例如飞机轮船等的外形曲线设计 如何确定三次样条函数在每一个小区间上的三次多项式函数 的系数 （略） 2. 利用MATLAB软件进行插值计算 1)一维插值 MATLAB中的插值函数为interp1，其调用格式为 yi=interp1(x, y, xi, method) 其中x，y为插值点，yi为在被插值点xi处的插值结果； x,y为向量，meathod表示采用的插值方法，MATLAB 提供的插值方法有几种:linear线性插值；spline三次样 条插值；cubic立方插值缺省时表示线性插值. 注意：所有的插值方法都要求x是单调的，并且xi不能够 超过x的范围。 例 在一 天24小时内，从零点开始每间隔2小时测得的环 境温度数据分别为 12，9，9，1 0 ，18 ，24，28，27，25，20， 18，15 ,13 推测 下午1点（即13点）时的温度 键人： x=0：2：24； y=12 9 9 10 18 24 28 27 25 20 18 15 13； x113 ； y1interp1（x，y，xl, spline） 输出的y1即为所求 若要得到一天24小时的温度曲线，只需继续键人： x20：13600：24； y2=interp1(x, y, x2, spline); plot(x,y,o,x2,y2) 2) 高维插值 N维插值函数interpN（）其中N可以为2，3， 如N2为二维插值， 调用格式为: zi=interp2(x,y,z,xi, yi, method) 其中x, y 为原始数据的第一和第二维坐标，z 为函数值， 函数返回在xi, yi所指定位置插值所得的函数值， method表示采用的插值方法，最邻近插值，线性插值， 双三次插值。缺省时表示线性插值所有的插值方法都 要求x和y是单调的网格，x和 y可以是等距的也可以是不 等距的 例 气旋变化情况的可视化 下表是气象学家测量得到的气象资料，它们分别表示 在南半球按不同纬度, 不同月份的平均气旋数字根据 这些数据，绘制出气旋分布曲面图形 y=5:10:85; x=1:12; z=2.4, 1.6, 2.4, 3.2, 1.0, 0.5, 0.4, 0.2, 0.5, 0.8, 2.4, 3.6;, 0.3, 0, 0, 0.3, 0, 0, 0.1, 0.2, 0.3, 0, 0.1, 0.3 xi,yi=meshgrid(1:0.1:12,5:85); zi=interp2(x, y, z, xi, yi,); mesh(xi, yi, zi) xlabel(月份) ylabel(纬度) zlabel(气旋) axis(0 12 0 90 0 50) title(南半球气旋可视化图形) x = 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 y = 5 5 5 5 5 5 5 5 5 5 5 5 15 15 15 15 15 15 15 15 15 15 15 15 25 25 25 25 25 25 25 25 25 25 25 25 35 35 35 35 35 35 35 35 35 35 35 35 45 45 45 45 45 45 45 45 45 45 45 45 55 55 55 55 55 55 55 55 55 55 55 55 65 65 65 65 65 65 65 65 65 65 65 65 75 75 75 75 75 75 75 75 75 75 75 75 85 85 85 85 85 85 85 85 85 85 85 85 x,y=meshgrid(1:12,5:10:85); 二、拟合 曲线拟合： 根据一组二维数据，即平面上的若干点，要求确定 一个一元函数yf(x),即曲线y=f(x)，使这些点与曲线总体来说 尽量接近。 在化学反应中，为研究某化合物的浓度随时间的变化规律， 测得一组数据如下表 1. 引例 t（分）12345678 y46.48.08.49.289.59.79.86 t（分）910111213141516 y1010.210.3210.4210.510.5510.5810.6 时间与浓度的散点图 根据散点图，可以假定y与t的函数为 如何具体地确定函数关系？ 2. 最小二乘法 最小二乘法的原理是求s(x),使 达到最小. 已知一组数据，用什么样的曲线拟合最好呢？可以根据散 点图进行直观判断，在此基础上，选择几种曲线分别拟合， 然后观察哪条曲线的最小二乘指标最小. 常用的拟合函数 多项式函数 双曲线 指数曲线等 3、拟合原理 由 可知 因此可假设 因此求最小二乘解转化为 二次函数 由多元函数取极值的必要条件 得 即 -(4) 即 引入记号 则由内积的概念可知 -(5) -(6) 显然内积满足交换律 方程组(4)便可化为 -(7) 将其表示成矩阵形式 -(8) 2. 曲线拟合的MATLAB实现 引例的求解 t=l:16； y= 4 6.4 8 8.4 9.28 9.5 9.7 9.86 10 10.2 10.32 10.42 10.5 10.55 10.58 10.6 ； p=polyfit(t,y,2) (二次多项式拟合） 计算结果： p=-0.0445 1.0711 4.3252 %二次多项式的系数 由此得到某化合物的浓度y与时间t的拟合函数 对函数的精度如何检测呢？仍然以图形来检测 3，非线性拟合命令：lsqcurvefit、lsqnonlin （1）lsqcurvefit: c = lsqcurvefit ( fun, x0, xdata, ydata) 其中fun为拟合函数的M - 函数文件名, x0为 初始向量, 也即拟合解是靠迭代求解得到的， 初始值的选取好坏直接影响最终的求解。 xdata, ydata为参与曲线拟合的实验数据。 函数返回值c为非线性函数fun的拟合系数。 例2: 2004年全国大学生数学建模竞赛C题( 酒 后驾车)中给出某人在短时间内喝下两瓶啤 酒后,间 隔一定的时间t测量他的血液中酒精含量y( 毫克/百 毫升) ,得到数据如表1所示。 根据微分方程模型得到短时间内喝酒后血液中酒精 浓度与时间的关系为 现通过已有数据来拟合C_1,C_2,C_3 先建立拟合函数examp21 function f = examp21(c, tdata) f=c(1).*(exp(-c(2).*tdata)-exp(-c(3).*tdata); 然后运行以下程序： tdata = 0.25 0.5 0.75 1 1.5 2 2.5 3 3.5 4 4.5 5 6 7 8 9 10 11 12 13 14 15 16 ; ydata = 30 68 75 82 82 77 68 68 58 51 50 41 38 35 28 25 18 15 12 10 7 7 4 ; c0=1 1 1; for i = 1: 50 c = lsqcurvefit(examp21,c0,tdata,ydata); c0 = c End 求解为：c = 114.2472 0.1852 2.0126 （2）lsqnonlin用法：lsqsnonlin ( fun, c0) 。 求含参量非线性函数fun中的参量c,使得各数据点 函数值fun的平方和最小。 例如用lsqsnonlin ( fun, c0)命令求解上例 的过程 如下: 先建立拟合函数，注意跟上面的不同 function ff=examp22(c,tdata,ydata) tdata = 0.25 0.5 0.75 1 1.5 2 2.5 3 3.5 4 4.5 5 6 7 8 9 10 11 12 13 14 15 16 ; ydata = 30 68 75 82 82 77 68 68 58 51 50 41 38 35 28 25 18 15 12 10 7 7 4 ; ff=c(1).*(exp(-c(2).*tdata)-exp(-c(3).*tdata)-ydata; 然后运行： c0=1 1 1; for i = 1: 50 c = lsqnonlin(examp22,c0); c0 = c end 求解为：c = 114.2472 0.1852 2.0126 例 估计水塔的水流量 某居民区的民用自来水是由一个圆柱形的水塔提供水塔高 12.2米，直径17.4米水塔是由水泵根据水塔内水位高低自动 加水，一般每天水泵工作两次现在需要了解该居民区用水规律 与水泵的工作功率按照设计，当水塔的水位降至最低水位，约 8.2米时，水泵自动启动加水；当水位升高到一个最高水位，约 10.8米时，水泵停止工作可以考虑用用水率（单位时间的用 水量）来反映用水规律，并通过间隔一段时间测量水塔里的水位 来估算用水率下表是某一天的测量记录数据，测量了28个时 刻，但是由于其中有4个时刻遇到水泵正在向水塔供水，而无水位 记录（表1中用符号表示）试建立合适的数学模型，推算任意 时刻的用水率、一天的总用水量和水泵工作功率 表1 水位测测量记录记录 时时刻(h)00.921.842.953.874.985.90 水位(cm)968948931913898881869 时时刻(h)7.017.938.979.9810.9210.9512.03 水位(cm)852839822/10821050 时时刻(h)12.9513.8814.9815.9016.8317.9319.04 水位(cm)1021994965941918892866 时时刻(h)19.9620.8420.0122.9623.8624.9925.91 水位(cm)843822/105910351018 问题分析与数据处理 由问题的要求，关键在于确定用水率函数，即单位时间内 用水体积，记为f(t)，如果能够通过测量数据，产生若干个 时刻的用水率，也就是f(t)在若干个点的函数值，则f(t)的计 算问题就可以转化为插值或拟合问题 1假设 1）水塔中水流量是时间的连续光滑函数，与水泵工作 与否无关，并忽略水位高度对水流速度的影响 2）水泵工作与否完全取决于水塔内水位的高度. 3）水塔为标准圆柱体 考虑到假设2）结合表1中具体数据，推断得出水泵第一次供 水时间段为9，11，第二次供水时间段为20.8 ，23 4）影响水箱流量的唯一因素是居民对水的普通要求. 5）由于水塔截面积是常数，先算出单位时间流出的水的高度. 分析 一天有两个供水段和三个水泵不工作的时间段（分别成为第一， 第二段，第三段） 一、二段数据多直接用多项式拟合得水位函数，求出流量。 第三段因记录少,利用前段用外推的方法来处理，水泵工作段因 没有记录，也得另想办法 记号 t-输入时刻 h-水位纪录（cm) ci-第i时刻段拟合多项式的系数 ai-第i时刻段拟合多项式的导数的系数 步骤1）拟合1，2，3段的水位，并求出流量 2）拟合供水段的流量 3）求出一天总用水量 4）流量及总用水量的检验 5）结果 6）分析与改进 程序： t=0 0.92 1.84 25.91; h=968 948 931 1018;%排除供水时刻的4个数据 c1=polyfit(t(1:10) ，h(1:10), 3); a1=polyder(c1); tp1=0:0.1:9 x1=-polyval(a1,tp1); c2=polyfit(t(11:21)， h(11:21), 3); a2=polyder(c2); tp2=11:0.1:20.8 x2=-polyval(a2,tp2); % 在第1供水段之前和之后各取几点，用来拟合第一供水段 %的流量 xx1=-polyval(a1,8,9); xx2=-polyval(a2,11,12); xx12=xx1,xx2; t12=8 9 11 12; c12=polyfit(t12,xx12,3); tp12=9:0.1:11 x12=polyval(c12,tp12); %第二段供水时段前取两点，该时段后，因记录少，用差分得到 %流量，然后取两点，用这4个数值拟合第2供水段的流量 dt3=diff(t(22:24); dh3=diff(h(22:24); dht3=-dh3./dt3; t3=20 20.8 t(22) t(23); xx3=- polyval(a2,t3(1:2),dht3); c3=polyfit(t3,xx3,3); t3=20.8:0.1:26; x3=polyval(c3,tp3) %用数值积分计算一天的总用水量 y1=0.1*trapz(x1); y2=0.1*trapz(x2); y12=0.1*trapz(x12); y3=0.1*trapz(x3); s=17.4*17.4*pi; % 圆柱的底面积 y=(y1+y2+y12+y3)*s*0.01 %m3 流水量及水量的检验： 1） 模型用水量为用水率函数f(t)的积分；实际用水量为水 塔内水的体积之差 2） 两次充水期间，水泵的注人水量因该相差不大 水泵充水量=充水后的水量+充水期间的流出量-充水前的水量 计算结果（略） 分析与改进（略） 三、Matlab中方程求解的基本 命令 1.roots(p) %求多项式的根，其中p是多项式向量。 例求 的根。 解：roots(1,-1,1,-1) 注： 1,-1,1,-1在matlab中表示多项式 2.solve(fun) %求方程fun=0的符号解，如果不 能求得精确的符号解，可以计算可变精度的数值解 例：用solve