《数值计算插值与拟合(给药方案估计水塔的水流量)ppt课件.ppt》由会员分享，可在线阅读，更多相关《数值计算插值与拟合(给药方案估计水塔的水流量)ppt课件.ppt（60页珍藏版）》请在三一办公上搜索。
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,
2、y1, y2, ., yn，要求一个次数不超过n次的代数多项式 Pn(x)=a0+a1x+a2x2+anxn使之在节点上满足Pn(xi)=f(xi),几种常用的多项式插值,拉格朗日插值：,牛顿插值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 ( x )各节点的值
3、分别为： 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.拟合的基本原理,1.拟合问题引例,3.用MATLAB求解拟合问题,4.应用举例,5.插值与拟合的比较,拟合问题引例一 电阻问题,温度t(0C) 20.5 32.7 51.0 73.0 95.7电阻
4、R() 765 826 873 942 1032,已知热敏电阻电阻值与温度的数据：,求600C时的电阻R。,设 R=at+ba,b为待定系数,解答,拟合问题引例二 给药问题,求血药浓度随时间的变化规律c(t).,作半对数坐标系(semilogy)下的图形,MATLAB(aa1),解答,曲 线 拟 合 问 题 的 提 法,已知一组（二维）数据，即平面上 n个点（xi,yi) i=1,n, 寻求一个函数（曲线）y=f(x), 使 f(x) 在某种准则下与所有数据点最为接近，即曲线拟合得最好。,y=f(x),i 为点（xi,yi) 与曲线 y=f(x) 的距离,曲线拟合问题最常用的解法线性最小二乘法
5、的基本思路,第一步:先选定一组函数 r1(x), r2(x), rm(x), mn, 令 f(x)=a1r1(x)+a2r2(x)+ +amrm(x) （1） 其中 a1,a2, am 为待定系数。,第二步: 确定a1,a2, am 的准则（最小二乘准则）： 使n个点（xi,yi) 与曲线 y=f(x) 的距离i 的平方和最小 。,记,问题归结为，求 a1,a2, am 使 J(a1,a2, am) 最小。,线性最小二乘法的求解,定理：当RTR可逆时，超定方程组 存在最小二乘解，且其解可表示为下列形式： a=(RTR)-1RTy a=Ry,曲线拟合的最小二乘法要解决的问题，实际上就是求以下超定
6、方程组的最小二乘解的问题。,Ra=y,其中,Ra=y,线性最小二乘拟合 f(x)=a1r1(x)+ +amrm(x)中函数r1(x), rm(x)的选取,1. 通过机理分析建立数学模型来确定 f(x)；,2. 将数据 (xi,yi) i=1, n 作图，通过直观判断确定 f(x)：,用MATLAB解拟合问题,1、线性最小二乘拟合,2、非线性最小二乘拟合,用MATLAB作线性最小二乘拟合,1. 作多项式f(x)=a1xm+ +amx+am+1拟合,可利用已有程序:,a=polyfit(x,y,m),2. 对超定方程组,3.多项式在x处的值y可用以下命令计算： y=polyval（a，x）,1.
7、lsqcurvefit已知数据点： xdata=（xdata1，xdata2，xdatan）， ydata=（ydata1，ydata2，ydatan）,用MATLAB作非线性最小二乘拟合,Matlab的提供了两个求非线性最小二乘拟合的函数：lsqcurvefit和lsqnonlin。两个命令都要先建立M-文件fun.m，在其中定义函数f(x)，但两者定义f(x)的方式是不同的,可参考例题.,lsqcurvefit用以求含参量x（向量）的向量值函数F(x,xdata)=（F（x，xdata1），F（x，xdatan）T中的参变量x(向量),使得,输入格式为: (1) x = lsqcurvef
8、it (fun,x0,xdata,ydata); (2) x =lsqcurvefit (fun,x0,xdata,ydata,options); (3) x = lsqcurvefit (fun,x0,xdata,ydata,options,grad); (4) x, options = lsqcurvefit (fun,x0,xdata,ydata,); (5) x, options,funval = lsqcurvefit (fun,x0,xdata,ydata,); (6) x, options,funval, Jacob = lsqcurvefit (fun,x0,xdata,yda
9、ta,);,说明：x = lsqcurvefit (fun,x0,xdata,ydata,options);,lsqnonlin用以求含参量x（向量）的向量值函数 f(x)=(f1(x),f2(x),fn(x)T 中的参量x，使得 最小。 其中 fi（x）=f（x，xdatai，ydatai） =F(x,xdatai)-ydatai,2. lsqnonlin,已知数据点： xdata=（xdata1，xdata2，xdatan） ydata=（ydata1，ydata2，ydatan）,输入格式为： 1） x=lsqnonlin（fun，x0）； 2） x= lsqnonlin （fun，x0
10、，options）； 3） x= lsqnonlin （fun，x0，options，grad）； 4） x，options= lsqnonlin （fun，x0，）； 5） x，options，funval= lsqnonlin （fun，x0，）；,说明：x= lsqnonlin （fun，x0，options）；,例2 用下面一组数据拟合 中的参数a，b，k,该问题即解最优化问题：,MATLAB(fzxec1),1）编写M-文件 curvefun1.m function f=curvefun1(x,tdata) f=x(1)+x(2)*exp(-0.02*x(3)*tdata) %其中
11、x(1)=a; x(2)=b；x(3)=k;,2）输入命令tdata=100:100:1000cdata=1e-03*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59; x0=0.2,0.05,0.05; x=lsqcurvefit (curvefun1,x0,tdata,cdata) f= curvefun1(x,tdata),F(x，tdata)= ，x=(a，b，k),解法1. 用命令lsqcurvefit,3）运算结果为：f =0.0043 0.0051 0.0056 0.0059 0.0061 0.0062 0.0062 0.0063
12、 0.0063 0.0063 x = 0.0063 -0.0034 0.2542,4）结论：a=0.0063, b=-0.0034, k=0.2542,MATLAB(fzxec2),解法 2 用命令lsqnonlin f(x)=F(x,tdata,ctada)= x=（a，b，k）,1）编写M-文件 curvefun2.m function f=curvefun2(x) tdata=100:100:1000; cdata=1e-03*4.54,4.99,5.35,5.65,5.90, 6.10,6.26,6.39,6.50,6.59; f=x(1)+x(2)*exp(-0.02*x(3)*td
13、ata)- cdata,2）输入命令: x0=0.2,0.05,0.05;x=lsqnonlin(curvefun2,x0)f= curvefun2(x),函数curvefun2的自变量是x，cdata和tdata是已知参数，故应将cdata tdata的值写在curvefun2.m中,3）运算结果为 f =1.0e-003 *(0.2322 -0.1243 -0.2495 -0.2413 -0.1668 -0.0724 0.0241 0.1159 0.2030 0.2792 x =0.0063 -0.0034 0.2542,可以看出,两个命令的计算结果是相同的.,4）结论：即拟合得a=0.0
14、063 b=-0.0034 k=0.2542,MATLAB解拟合问题应用实例,1、电阻问题,2、给药方案问题,MATLAB(dianzu1),引例一 电阻问题 解答,得到 a1=3.3940, a2=702.4918,方法2. 直接用,结果相同。,MATLAB(dianzu2),一室模型：将整个机体看作一个房室，称中心室， 室内血药浓度是均匀的。快速静脉注射后，浓度立即上升；然后迅速下降。 当浓度太低时，达不到预期的治疗效果；当浓度太高，又可能 导致药物中毒或副作用太强。 临床上，每种药物有一个最小有效浓度c1和一个最大有效浓度c2。设计给药方案时，要使血药浓度 保持在c1c2之间。 本题设c
15、1=10，c2=25(ug/ml).,引例二 给药方案 解答,一种新药用于临床之前，必须设计给药方案.,药物进入机体后血液输送到全身，在这个过程中不断地被吸收、分布、代谢，最终排出体外，药物在血液中的浓度，即单位体积血液中的药物含量，称为血药浓度。,在实验方面,对某人用快速静脉注射方式一次注入该药物300mg后,在一定时刻t(小时)采集血药,测得血药浓度c(ug/ml)如下表:,要设计给药方案,必须知道给药后血药浓度随时间变化的规律。从实验和理论两方面着手：,给药方案,1. 在快速静脉注射的给药方式下，研究血药浓度 （单位体积血液中的药物含量）的变化规律。,t,问题,分析,理论：用一室模型研究
16、血药浓度变化规律,实验：对血药浓度数据作拟合，符合负指数变化规律,2. 给定药物的最小有效浓度和最大治疗浓度，设计给药方案：每次注射剂量多大；间隔时间多长。,3. 血液容积 v, t=0时注射剂量 d, 血药浓度即为 d/v.,2. 药物排除速率与血药浓度成正比，比例系数 k(0),模型假设,1. 机体看作一个房室，室内血药浓度均匀一室模型,模型建立,在此，d=300mg，t及c（t）在某些点处的值见前表，需经拟合求出参数k、v,用线性最小二乘拟合c(t),MATLAB(lihe1),计算结果：,用非线性最小二乘拟合c(t),给药方案设计,设每次注射剂量D, 间隔时间,血药浓度c(t) 应c1
17、 c(t) c2,初次剂量D0 应加大,给药方案记为：,2、,1、,计算结果：,给药方案：,c1=10,c2=25k=0.2347v=15.02,故可制定给药方案：,即: 首次注射 375 mg， 其余每次注射 225 mg， 注射的间隔时间为 4 小时。,用非线性最小二乘拟合c(t)-用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(c
18、urvefun3,x0,tdata,cdata);f=curvefun3(x,tdata) x,MATLAB(lihe2),1、用M-文件curvefun3.m定义函数function f=curvefun3(x,tdata)d=300f=(x(1)d)*exp(-x(2)*tdata) % x(1)=v; x(2)=k,MATLAB(FZXEC3),拟合与插值的比较,问题：给定一批离散的数据点，需确定满足特定要求的曲线或曲面, 从而获取整体的规律。即通过窥几斑来达到知全豹。,解决方案：,若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，这就是数据拟合，又称曲线拟合或曲面拟合
19、。,若要求所求曲线（面）通过所给所有数据点，就是插值问题；,从几何意义上看，拟合是给定了空间中的一些点，找到一个已知形式的连续曲面来最大限度地逼近这些点；而插值是找到一个(或几个分片光滑的)连续曲面来穿过这些点。,拟合与插值的区别,函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。,实例：下面数据是某次实验所得，希望得到X和 f之间的关系？,MATLAB(cn),最临近插值、线性插值、样条插值与曲线拟合结果：,美国某州的各公用水管理机构要求各社区提供各个时刻的用水率以及每天所用的总用水量但许多社区并没有测量流入或流出当地水塔的水量的设备
20、，他们只能代之以每小时测量水塔中的水位更为重要的是，无论什么时候，只要水塔中的水位下降到某一最低水位L时，水泵就启动向水塔重新充水直到某一最高水位，但也无法得到水泵的供水量的测量数据因此，在水泵正在工作时，人们不容易建立水塔中水位与水泵工作时的用水量之间的关系水泵每天向水塔充水一次或两次，每次大约二小时试估计在任何时候，甚至包括水泵正在工作的时间内，水从水塔流出的流量，并估计一天的总用水量表给出了某个小镇某一天的真实数据,估计水塔的水流量,表某小镇某天的水塔水位,表给出了从第一次测量开始的以秒为单位的时刻，以及该时刻的高度单位为百分之一英尺的水塔中水位的测量值，例如3316秒后，水塔中的水位达
21、到31.10英尺水塔是一个垂直圆形柱体，高为40英尺，直径为57英尺,二、问题分析,我们很容易想到应通过对所给的数据进行数值拟合来建模在讨论具体的建模方法以前，我们应先给出一些合理的假设(1)影响水从水塔中流出的流量的唯一因素是公众对水的传统要求因为表只给出了某一天（近26小时）水塔的水位数据，并没有对这些数据的产生有影响的因素作出具体的说明，我们只能假定所给数据反映了有代表性的一天，而不包括任何特殊情况，如自然灾害、火灾、水塔溢水、水塔漏水等对水的特殊要求(2)水塔中的水位不影响水流量的大小，气候条件、温度变化等也不影响水流量(3)水泵工作起止时间有它的水位决定，每次充水时间大约为两个小时,
22、二、问题分析,(4)水泵充水速度恒定，且水泵充水的水流量远大于水塔的水流量，以保证人们对水的需求水泵工作时不需要维修，也不中途停止工作(5)水塔的水流量与水泵状态独立，并不因水泵工作而增加或减少水流量的大小(6)水塔的水流量曲线可以用一条光滑的曲线了逼近这时，在每一个数据点，水流量的两阶导数是连续的因为水的消耗量是基于社区公众一天的活动，如洗澡、做饭、洗衣服等，每一个使用者的要求与整个社区的要求相比是微不足道的，而整个社区的需求是不可能同时增加或减少的，由于水的消耗的自然性，可以假设水流量曲线是一条连续光滑的曲线(7)表的数据是准确的,二、问题分析,对所给的问题，其建模方法是经典的，基本上是分
23、成三步：首先由所给数据得到在各数据点处的水流量，然后找出一个水从水塔流出的水流量的光滑拟合逼近，最后处理水泵工作时的充水水量以及一天该小镇公众的总用水量，同时也重建了水泵工作时所缺的数据所给数据的初步处理.我们把表所给的数据作为时间的函数画成图,二、问题分析,图1时间与水位的关系图,二、问题分析,从图可以看出，最大的困难是要解决如何描述水塔充水期间的水流量的行为，为此，我们先分析一下水泵充水期间的观察数据，要解决两个问题：一是两次充水准确的起始时间和停止时间，如果无法得到准确时间的话，以哪一时刻作为起止时间比较合理；二是充水期间的水流量如何描述从所给的数据自然无法知道水泵开始和停止的准确时间，
24、但是已知第一次充水前的最后一个数据为32284秒时水位为26.97英尺 ，充水中第二个数据为39332秒时而39332-32284=7048秒，即约为1.96小时，由水泵每次充水要大约小时可知，水泵是在32284秒时开始充水的停止时间在39332秒与39435秒之间，但这两个时刻的差距为103秒，约0.028小时，很短的时间，所以我们可以假定水泵停止工作时间为39332秒充水开始时水塔水位为26.97英尺，可以认为L大约为27.00英尺,二、问题分析,、水流量曲线的拟合 表给出的是水位与时间的关系，而题目要求我们求出的是水流量与时间的关系，因此，我们先将表的数据转化为水塔中水的体积与时间的关系
25、，然后再转化为水流量与时间的关系表、图代表水的体积与时间的关系（程序见实验解答中程序二）,二、问题分析,表2时间与体积的关系,二、问题分析,图2时间与体积的关系图,二、问题分析,我们用 从水的体积与时间的关系得到水流量与时间的关系（由于在充水时，没有水的体积与时间的关系，所以也没有水流量与时间的关系）我们采用差分法来解决这个问题由于水泵充水两次，数据被分割成三组，因而我们也分三组来处理数据对每一组数据，我们采用中心差分公式,二、问题分析,来计算每一组中间数据点的水流量而对每组前两个和最后两个数据点，采用如下的公式来计算,对于最后的倒数第二个数据，我们用下面的公式计算：,二、问题分析,(4)水泵
26、充水速度恒定，且水泵充水的水流量远大于水塔的水流量，以保证人们对水的需求水泵工作时不需要维修，也不中途停止工作(5)水塔的水流量与水泵状态独立，并不因水泵工作而增加或减少水流量的大小(6)水塔的水流量曲线可以用一条光滑的曲线了逼近这时，在每一个数据点，水流量的两阶导数是连续的因为水的消耗量是基于社区公众一天的活动，如洗澡、做饭、洗衣服等，每一个使用者的要求与整个社区的要求相比是微不足道的，而整个社区的需求是不可能同时增加或减少的，由于水的消耗的自然性，可以假设水流量曲线是一条连续光滑的曲线(7)表的数据是准确的,二、问题分析,计算结果见表和图,二、问题分析,图3时间与流量的关系,5、三次样条插值,6、对水泵两段充水时间的处理,（1）第一次平均水流量：,（2）第二次平均水流量：,（3）平均水流量：,7、一天总用水量,8、检验,以不同的时间为起点得到的一天总用水量相差多少,
