插值与拟合给药方案估计水塔的水流量

1、一、插值1、插值问题： 不知道某一函数f(x)在待定范围a,b上 的具体表达式，而只能通过实验测量得到该函数在一系列点ax1, x2 , ., xn b上的值 y0, y1, y2, ., yn，需要找一个简单的函数P(x)来近似地代替f(x)，要求满足： P(xi)=yi (i=1,2,.,n)，此问题称为插值问题。 P(x)称为f(x)的插值函数， x1, x2 , ., xn 称为插值节点，f(xi)称为插值条件。几种常用的插值方法 1、多项式插值 2、样条插值1、多项式插值方法 设y=f(x)在n+1个互异点上的x0 , x1, x2 , ., xn 上的值 y0, y1, y2, .

2、, yn，要求一个次数不超过n次的代数多项式 Pn(x)=a0+a1x+a2x2+anxn使之在节点上满足Pn(xi)=f(xi)几种常用的多项式插值拉格朗日插值：00( )()nninjjijiijxxpxyxx 牛顿插值Hermite插值2、样条插值方法 设给定区间a,b的一个分化： a=x0 x1xn=b， 如果函数s(x)满足条件：在每个子区间xi-1,xi上是k次多项式，且具有直到k-1阶的连续导数，则称s(x)为一个k次多项式样条。广泛使用的样条函数（1）二次样条（2）三次样条（3）B样条。二次样条的定义 设a,b 的一个划分：a=x0 x1, x2 , ., xn= b，函数f

3、( x )各节点的值分别为： f ( xi )=yi (i=1,2,.,n) 如果二次样条函数：满足： S ( xi )=yi (i=1,2,.,n) 三次样条函数的定义 设a,b 的一个划分：a=x0 x1, x2 , ., xn= b，函数f ( x )各节点的值分别为： f ( xi )=yi (i=1,2,.,n) 如果三次样条函数：3满足： S ( xi )=yi (i=1,2,.,n) 数据的拟合数据的拟合2.2.拟合的基本原理拟合的基本原理1.拟合问题引例拟合问题引例3.3.用用MATLABMATLAB求解拟合问题求解拟合问题4.4.应用应用举例举例5.5.插值与拟合的比较插值与

4、拟合的比较拟合问题引例一拟合问题引例一 电阻问题电阻问题温度温度t(0C) 20.5 32.7 51.0 73.0 95.7电阻电阻R( ) 765 826 873 942 1032已知热敏电阻电阻值与温度的数据：已知热敏电阻电阻值与温度的数据：求求60600C时的电阻时的电阻R。2040608010070080090010001100 设设 R=at+ba,b为待定系数为待定系数解答解答拟合问题引例二拟合问题引例二 给药问题给药问题 t (h) 0.25 0.5 1 1.5 2 3 4 6 8c ( g/ml) 19.21 18.15 15.36 14.10 12.89 9.32 7.45

5、5.24 3.01已知一室模型快速静脉注射下的血药浓度数据已知一室模型快速静脉注射下的血药浓度数据(t=0注射注射300mg)求血药浓度随时间的变化规律求血药浓度随时间的变化规律c(t).作半对数坐标系作半对数坐标系(semilogy)下的图形下的图形为待定系数kcectckt,)(002468100101102MATLAB(aa1)解答解答曲曲 线线 拟拟 合合 问问 题题 的的 提提 法法已知一组（二维）数据，即平面上已知一组（二维）数据，即平面上 n个点个点（xi,yi) i=1,n, 寻求一个函数（曲线）寻求一个函数（曲线）y=f(x), 使使 f(x) 在某种准则下与所在某种准则下与

6、所有数据点最为接近，即曲线拟合得最好。有数据点最为接近，即曲线拟合得最好。 +xyy=f(x)(xi,yi)i i 为点为点（xi,yi) 与与曲线曲线 y=f(x) 的距离的距离曲线拟合问题最常用的解法曲线拟合问题最常用的解法线性最小二乘法线性最小二乘法的基本思路的基本思路第一步: :先选定一组函数先选定一组函数 r1(x), r2(x), rm(x), m0)k(0)模型假设模型假设1. 1. 机体看作一个房室，室内血药浓度均匀一室模型模型建立模型建立3 c(0)d/由 假 设 得 ：dc2 -kcdt由 假 设得 ：( )ktdc tev 在此，d=300mg，t及c（t）在某些点处的值

7、见前表，需经拟合求出参数k、v用线性最小二乘拟合用线性最小二乘拟合c(t)( )ktdc tev12ln ,ln(/ )ycakadv lnln(/ )cdvkt2121,/ayataka vd e MATLAB(lihe1)计算结果：0.2347(1/ ),15.02( )khvld=300;t=0.25 0.5 1 1.5 2 3 4 6 8;c=19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01;y=log(c);a=polyfit(t,y,1)k=-a(1)v=d/exp(a(2)程序：用非线性最小用非线性最小二乘拟合二乘拟合c(t)给药

8、方案设计给药方案设计cc2c10t 设每次注射剂量D, D, 间隔时间 血药浓度c(t) 应应c1 c(t) c2 初次剂量D D0 0 应加大,0DD给药方案记为：kecc2112ln1cck2、)( ,1220ccDcD1、计算结果：9 . 3, 3 .225, 5 .3750DD)(4),(225),(3750hmgDmgD给药方案：c1=10,c2=25k=0.2347v=15.02故可制定给药方案：)(4),(225),(3750hmgDmgD即: 首次注射 375 mg， 其余每次注射 225 mg， 注射的间隔时间为 4 小时。用非线性最小二乘拟合用非线性最小二乘拟合c(t)-用

9、用lsqcurvefit(lsqnonlin)2、主程序lihe2.m如下cleartdata=0.25 0.5 1 1.5 2 3 4 6 8;cdata=19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01; x0=10,0.5;x=lsqcurvefit(curvefun3,x0,tdata,cdata);f=curvefun3(x,tdata) xMATLAB(lihe2)1 1、用M-M-文件curvefun3.m定义函数function f=curvefun3(x,tdata)d=300f=(x(1)d)*exp(-x(2)*tdata

10、) % x(1)=v; x(2)=k ( )ktdc tevMATLAB(FZXEC3)拟合与插值的比较问题：给定一批离散的数据点，需确定满足特定要求的曲线或曲面, 从而获取整体的规律。即通过窥几斑来达到知全豹。 解决方案： 若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，这就是数据拟合，又称曲线拟合或曲面拟合。 若要求所求曲线（面）通过所给所有数据点，就是插值问题； 从几何意义上看，拟合拟合是给定了空间中的一些点，找到一个已知形式的连续曲面来最大限度地逼近这些点；而插值插值是找到一个(或几个分片光滑的)连续曲面来穿过这些点。 拟合与插值的区别 函数插值与曲线拟合都是要根据

11、一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。 实例：下面数据是某次实验所得，希望得到X和 f之间的关系？x1247912131517f1.53.96.611.715.618.8 19.620.621.1MATLAB(cn)最临近插值、线性插值、样条插值与曲线拟合结果：最临近插值、线性插值、样条插值与曲线拟合结果：0246810121416180510152025已知数据点spline三次多项式插值0246810121416180510152025已知数据点linest三次多项式插值0246810121416180510152025已知数据点nearest三次

12、多项式插值美国某州的各公用水管理机构要求各社区提供各个时刻的用水率以及每天所用的总用水量但许多社区并没有测量流入或流出当地水塔的水量的设备，他们只能代之以每小时测量水塔中的水位更为重要的是，无论什么时候，只要水塔中的水位下降到某一最低水位L时，水泵就启动向水塔重新充水直到某一最高水位，但也无法得到水泵的供水量的测量数据因此，在水泵正在工作时，人们不容易建立水塔中水位与水泵工作时的用水量之间的关系水泵每天向水塔充水一次或两次，每次大约二小时试估计在任何时候，甚至包括水泵正在工作的时间内，水从水塔流出的流量，并估计一天的总用水量表给出了某个小镇某一天的真实数据估计水塔的水流量表某小镇某天的水塔水位

13、 表给出了从第一次测量开始的以秒为单位的时刻，以及该时刻的高度单位为百分之一英尺的水塔中水位的测量值，例如3316秒后，水塔中的水位达到31.10英尺水塔是一个垂直圆形柱体，高为40英尺，直径为57英尺 二、问题分析l我们很容易想到应通过对所给的数据进行数值拟合来建模在讨论具体的建模方法以前，我们应先给出一些合理的假设l(1)影响水从水塔中流出的流量的唯一因素是公众对水的传统要求因为表只给出了某一天（近26小时）水塔的水位数据，并没有对这些数据的产生有影响的因素作出具体的说明，我们只能假定所给数据反映了有代表性的一天，而不包括任何特殊情况，如自然灾害、火灾、水塔溢水、水塔漏水等对水的特殊要求l

14、(2)水塔中的水位不影响水流量的大小，气候条件、温度变化等也不影响水流量l(3)水泵工作起止时间有它的水位决定，每次充水时间大约为两个小时二、问题分析(4)水泵充水速度恒定，且水泵充水的水流量远大于水塔的水流量，以保证人们对水的需求水泵工作时不需要维修，也不中途停止工作(5)水塔的水流量与水泵状态独立，并不因水泵工作而增加或减少水流量的大小(6)水塔的水流量曲线可以用一条光滑的曲线了逼近这时，在每一个数据点，水流量的两阶导数是连续的因为水的消耗量是基于社区公众一天的活动，如洗澡、做饭、洗衣服等，每一个使用者的要求与整个社区的要求相比是微不足道的，而整个社区的需求是不可能同时增加或减少的，由于水

15、的消耗的自然性，可以假设水流量曲线是一条连续光滑的曲线(7)表的数据是准确的二、问题分析对所给的问题，其建模方法是经典的，基本上是分成三步：首先由所给数据得到在各数据点处的水流量，然后找出一个水从水塔流出的水流量的光滑拟合逼近，最后处理水泵工作时的充水水量以及一天该小镇公众的总用水量，同时也重建了水泵工作时所缺的数据所给数据的初步处理.我们把表所给的数据作为时间的函数画成图二、问题分析图1时间与水位的关系图 二、问题分析从图可以看出，最大的困难是要解决如何描述水塔充水期间的水流量的行为，为此，我们先分析一下水泵充水期间的观察数据，要解决两个问题：一是两次充水准确的起始时间和停止时间，如果无法得

16、到准确时间的话，以哪一时刻作为起止时间比较合理；二是充水期间的水流量如何描述从所给的数据自然无法知道水泵开始和停止的准确时间，但是已知第一次充水前的最后一个数据为32284秒时水位为26.97英尺 ，充水中第二个数据为39332秒时而39332-32284=7048秒，即约为1.96小时，由水泵每次充水要大约小时可知，水泵是在32284秒时开始充水的停止时间在39332秒与39435秒之间，但这两个时刻的差距为103秒，约0.028小时，很短的时间，所以我们可以假定水泵停止工作时间为39332秒充水开始时水塔水位为26.97英尺，可以认为L大约为27.00英尺 二、问题分析、水流量曲线的拟合

17、表给出的是水位与时间的关系，而题目要求我们求出的是水流量与时间的关系，因此，我们先将表的数据转化为水塔中水的体积与时间的关系，然后再转化为水流量与时间的关系表、图代表水的体积与时间的关系（程序见实验解答中程序二） 二、问题分析表2时间与体积的关系二、问题分析图2时间与体积的关系图二、问题分析我们用 从水的体积与时间的关系得到水流量与时间的关系（由于在充水时，没有水的体积与时间的关系，所以也没有水流量与时间的关系）我们采用差分法来解决这个问题由于水泵充水两次，数据被分割成三组，因而我们也分三组来处理数据对每一组数据，我们采用中心差分公式 dtdVtf)( )dVf tdt)(128812112i

18、iiiiiittVVVVf二、问题分析来计算每一组中间数据点的水流量而对每组前两个和最后两个数据点，采用如下的公式来计算)(243121iiiiiittVVVf)(243121iiiiiittVVVf111122iiiiiiitttVVVf对于最后的倒数第二个数据，我们用下面的公式计算：二、问题分析(4)水泵充水速度恒定，且水泵充水的水流量远大于水塔的水流量，以保证人们对水的需求水泵工作时不需要维修，也不中途停止工作(5)水塔的水流量与水泵状态独立，并不因水泵工作而增加或减少水流量的大小(6)水塔的水流量曲线可以用一条光滑的曲线了逼近这时，在每一个数据点，水流量的两阶导数是连续的因为水的消耗量

19、是基于社区公众一天的活动，如洗澡、做饭、洗衣服等，每一个使用者的要求与整个社区的要求相比是微不足道的，而整个社区的需求是不可能同时增加或减少的，由于水的消耗的自然性，可以假设水流量曲线是一条连续光滑的曲线(7)表的数据是准确的二、问题分析计算结果见表和图二、问题分析图3时间与流量的关系 5、三次样条插值6、对水泵两段充水时间的处理（1）第一次平均水流量：（2）第二次平均水流量：（3）平均水流量：7、一天总用水量8、检验以不同的时间为起点得到的一天总用水量相差多少三次样条函数的定义 设a,b 的一个划分：a=x0 x1, x2 , ., xn= b，函数f ( x )各节点的值分别为： f (

20、xi )=yi (i=1,2,.,n) 如果三次样条函数：3满足： S ( xi )=yi (i=1,2,.,n) lsqnonlin用以求含参量x x（向量）的向量值函数 f(x)=(ff(x)=(f1 1(x),f(x),f2 2(x),(x),f,fn n(x)(x)T T 中的参量x x，使得 最小。 其中 fi（x）=f（x，xdatai，ydatai） =F(x,xdatai)-ydatai 22212( ) ( )( )( )( )Tnfx f xf xfxfx2. lsqnonlin已知数据点： xdata=xdata=（xdata1，xdata2，xdataxdatan n） ydata=ydata=（ydataydata1 1，ydataydata2 2，ydataydatan n）3 3）运算结果为 f =1.0e-003 f =1.0e-003 * *(0.2322 -0.1243 -0.2495 -0.2413 (0