1、数学实验之七数据拟合,2018/4/14,1,实验目的1 了解最小二乘拟合的原理，掌握用 MATLAB作线性最小二乘拟合的方法。2 通过实例学习如何用拟合方法解决 实际问题；通过实例理解参数辨识问题的几种方法。,2018/4/14,2,主要内容,范例2：薄膜渗透率的测定,布置实验,引例1,引例2,拟合的基本原理,拟合函数的选取,用MATLAB作拟合计算,范例1：静脉注射的给药方案,2018/4/14,3,求电阻R随温度t的变化规律。,已知热敏电阻数据：,引例1：热敏电阻电阻值的变化规律,2018/4/14,4,设 R=at+b a,b为待定系数,引例1：热敏电阻电阻值的变化规律,2018/4/
2、14,5,引例2：血药浓度的变化规律,已知一室模型快速静脉注射下的血药浓度数据(t=0注射300mg),求血药浓度随时间的变化规律c(t).,2018/4/14,6,作图观察,半对数坐标系(semilogy)下的图形,Log10c(t)=a t + b,2018/4/14,7,主要内容,范例2：薄膜渗透率的测定,布置实验,引例1,引例2,拟合的基本原理,拟合函数的选取,用MATLAB作拟合计算,范例1：静脉注射的给药方案,2018/4/14,8,曲线拟合的基本原理问题的提法,已知一组（二维）数据，即平面上 n个点（xi,yi) i=1,n, 寻求一个函数（曲线）y=f(x), 使 f(x) 在
3、某种准则下与所有数据点最为接近，即曲线拟合得最好。,y=f(x),i 为点（xi,yi) 与曲线 y=f(x) 的距离,2018/4/14,9,问题的数学模型,步骤：1) 选定一类函数 f(x，a1，a2， ，am) （1） 其中 a1,a2, am 为待定常数。,2) 确定参数a1,a2, am， 准则(最小二乘准则）：使n个点（xi,yi) 与曲线 y=f(x ，a1，a2， ，am) 的距离 i 的平方和最小 。,2018/4/14,10,记,问题归结为: 求 a1,a2, am 使 J(a1,a2, am) 最小。 这样的拟合称为最小二乘拟合。,问题的数学模型,2018/4/14,11
4、,特别，若选定一组函数 r1(x), r2(x), rm(x), m0);,2. 血液容积v, t=0瞬时注射剂量d, 血药浓度立即为d/v.,模型假设：,由假设1，,3.快速静脉注射下一室模型的血药浓度：c(t),c,由假设2，,2018/4/14,28,若c1=10, c2=25(g/ml), 给药方案设计归结为根据数据(ti,ci) i=1,n (d 给定）拟合曲线c(t), 以确定系数k, v.,给药方案 设计：D0:初次剂量； :注射时间间隔；D:重复注射剂量。,范例1：静脉注射的给药方案,2018/4/14,29,给药方案,记,则,思考:取对数化为线性最小二乘,对结果有什么影响？,
5、2018/4/14,30,主要内容,范例2：薄膜渗透率的测定,布置实验,引例1,引例2,拟合的基本原理,拟合函数的选取,用MATLAB作拟合计算,范例1：静脉注射的给药方案,2018/4/14,31,范例2：薄膜渗透率的测定,1问题背景； 2假设；3问题的分析 4数学模型；5参数辨识方法： 函数拟合法 非线性规划法； 导函数拟合法； 线性化迭代法。6结果分析,2018/4/14,32,某种医用薄膜在试制时需测定其被物质分子穿透的能力。 测定方法：用面积为S的薄膜将容器分成两部份，在两部分中分别注满该物质的两种不同浓度的溶液。此时该物质分子就会从一侧向另一侧扩散。平均每单位时间通过单位面积薄膜的
6、物质分子量与膜两侧溶液的浓度差成正比，比例系数K称为渗透率。定时测量容器中薄膜某一侧的溶液浓度，以此确定K。,问题背景,薄膜渗透率的测定,2018/4/14,33,3）薄膜是双向同性的即物质从膜的任何一侧向另一侧渗透的性能是相同的。,假设：,1）薄膜两侧的溶液始终是均匀的；,2）平均每单位时间通过单位面积薄膜的物质分子量与膜两侧溶液的浓度差成正比。,2018/4/14,34,第一步：通过机理分析确定容器一侧的浓度随时间的变化规律： C=CB(t),第二步：利用实验数据(tj,CB(tj)，和函数拟合的方法求出其中的未知参数，包括渗透率K.,解决问题的思路：,2018/4/14,35,考察时段t
7、，t+t薄膜的一侧容器中该物质质量的变化。,1）以容器A侧为例，在时段t，t+t物质质量的增量为：,问题分析：,由假设2,在时段t，t+t，从B侧渗透至A侧的该物质的质量为：,2018/4/14,36,于是有：,两边除以t，并令t0取极限再稍加整理即得：,（1）,问题分析：,2018/4/14,37,分别表示在初始时刻两侧溶,其中,2) 注意到任意时刻，整个溶液中含有该 物质的质量,与初始时刻该物质的含量相同，因此,问题分析：,液的浓度,2018/4/14,38,便可得出CB(t)的变化规律，从而根据实验数据进行拟合，估计出参数K, 。,问题分析：,2018/4/14,39,基于假设和前面的分
8、析，B侧的浓度CB(t)应满足如下微分方程和初始条件：,数学模型,2018/4/14,40,模型求解方法,前面得到的模型是一个带初值的一阶线性微分方程，解之得：,问题归结为利用CB在时刻tj的测量数据Cj(j=1,2,.,N)来辨识 K 和 。,1. 函数拟合法,2018/4/14,41,引入,从而,用函数CB(t)来拟合所给的实验数据，从而估计出其中的参数a,b,K。,若,将其代入上式有：,拟合函数化简:,2018/4/14,42,用MATLAB软件进行计算.1）编写函数M-文件 nongdu.mfunction f = nongdu(x,tdata)f = x(1)+x(2)*exp(-0
9、.02*x(3)*tdata);其中 x(1) = a;x(2) = b;x(3) = k;2) 编写M文件 (baomo.m) tdata = linspace(100,1000,10); cdata =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(nongdu,x0,tdata,cdata),编写MATLAB程序:,2018/4/14,43,3) 输出结果: x = 0.007 -0.003 0.1012 即 k = 0.1012, a = 0.007, b
10、= -0.003,结果:,进一步求得：,To MATLAB (nongdu, baomo1),2018/4/14,44,2. 非线性规划法,利用CB在时刻tj的测量数据Cj(j=1,2,.,N)来辨识 K 和 。,问题可转化为求函数,即求函数,的最小值点（K，a，b）。,2018/4/14,45,3. 导函数拟合法,2018/4/14,46,即为求参数K, a使下列误差函数达到最小：,该问题等价于用函 数 f(K,a,CB)=K(0.01a-0.02CB)来拟合数据,(j=1, 2,.,N),用上述方法求出参数K, a.,3. 导函数拟合法,2018/4/14,47,4. 线性化迭代法,前面带
11、初始条件的一阶线性微分方程的解为,其中：,如果得到了参数K的一个较好的近似值K*，则将CB(t)关于K在K*处展开，略去K的二次及以上的项得CB(t)的一个近似式,2018/4/14,48,4. 线性化迭代法,通过极小化,确定a, b, d, 再由K=d/0.02b得到K*的修正值K。K*K*- K, 得到K的一个新的近似值，用同样的方法再求新的修正值K。这个过程可以不断重复，直到修正值足够小为止。,2018/4/14,49,1）当K的初值取为k=0.3时，出现奇异情况，迭代不收敛；,2）当K的初值取为k=0.2时，经四次迭代，已经收敛到一个很好的解。迭代结果如下表。,4. 线性化迭代法,20
12、18/4/14,50,3）取K的初值为0.1009,只一次迭代就得到2）中的最后结果。,4. 线性化迭代法,导函数拟合法得出的参数值精度有限，线性化迭代法要求参数的初值比较接近精确值。可将导函数拟合法和线性化迭代法结合起来，把前者得到的参数K的值作为迭代法中K的初值，可使其收敛或收敛更快。,2018/4/14,51,结果及误差分析,几种方法得出的结果及相应的误差总结于下表，误差为计算数据与实验数据之差的平方和。,2018/4/14,52,作出拟合函数曲线和数据点的图形tdata = linspace(100,1000,10);cdata = 1e-05.*454 499 535 565 590
13、. 610 626 639 650 659;x0 = 0.2,0.05,0.05;x,resnorm,residual = lsqcurvefit(nongdu,x0,tdata,cdata)t=linspace(100,1000,100);c=nongdu(x,t);plot(tdata,cdata,o,t,c)求拟合在各节点的误差平方和c1=nongdu(x,tdata);e=c1-cdata;e1=sum(e.*e),函数拟合法的误差分析程序,2018/4/14,53,函数拟合法的拟合效果,返回,2018/4/14,54,人口问题,马尔萨斯人口模型假设:人口出生率和死亡率是常数，因此人口
14、的净增长率为常数.设时刻t的人口数量为p(t)，人口出生率为b，死亡率为d ,(k=b-d),因此马尔萨斯人口模型如下：,该微分方程初值问题的解析解：,马尔萨斯人口模型,为了验证，采用某国家的人口历史数据,见下页表.,首先确定参数k：因为,所以,利用1790年和1840年的数据计算得出：,= 0.0294.,马尔萨斯人口模型,由此预测，1850年的人口数量为22,898,000，误差为1%，1900年的人口数量为 99,476,000，误差为31%，2000年的人口数量为 1,877,463,000，误差为567%，2050年将达到80多个亿？问题：（1）能不能用尽可能多的数据拟合微分方程的参
15、数？（2）初始数据能不能也作为参数处理？这样会得到更好的结果吗？,马尔萨斯人口模型,2018/4/14,59,处理方法,方法1：直接使用原始数据（t和p），使用非线性拟合求解方程的系数，并和原始数据作比较。方法2：因为 ，也就是lnp(t)=k(t-t0)+lnp0 。记q(t)=lnp(t),q0= lnp0 。上面的问题变成一个线性回归问题，可以得到相应的系数，和上面的方法已经原始数据比较。方法3：还是使用方法1，仅仅对时间做一个变换，比如17900,18001，等等。重新使用方法1进行计算。,1837年荷兰生物数学家Verhulst 提出改进. 人口增长率不应该是常数，假设增长率k是随着
16、人口数量接近最大数量M而线性递减.,逻辑斯谛人口模型,从而得到改进后的人口模型(逻辑斯谛增长模型),逻辑斯谛人口模型,可以求得，该人口模型的解为：,在这个方程中，将M，r，p0作为参数，使用前面类似的方法计算。并比较拟合数据和原始数据。,掌握用MATLAB作最小二乘多项式拟合和曲线拟合的方法；2.通过实例学习如何用拟合方法解决实际问题，注意与插值方法的区别。3.鼓励不囿于固定的模式或秩序，灵 活调整思路，突破思维的呆板性，找到打破常规的解法。并在文献检索、动手和动脑等方面得到锻炼，树立创新意识。,布置实验,实验目的,2018/4/14,62,实验内容,1. 下表给出了近两个世纪某国人口的统计数
17、据，试预测2010年该国的人口。,人口预测问题,2018/4/14,63,1) 可先用以上数据拟合Multhus人口指数增长模型，根据检验结果进一步讨论马尔萨斯人口模型的改进。2) Malthus 模型的基本假设是：人口的增长率为常数，记为r。记时刻t的人口为x(t），且初始时刻的人口为x0，于是得到如下微分方程：,2018/4/14,64,3) Malthus 模型导致人口总数的剧烈上升，为此对该模型做一些改进。其中一种改进是认为增长率应该是x的函数，并且随着x增大出现资源的竞争导致增长率减小。因此一般可以假设r(x)=r-x/k。得到如下微分方程：,2018/4/14,65,实验内容,2.
18、某年美国旧车价格的调查资料如下表其中xi表示轿车的使用年数，yi表示相应的平均价格。试分析用什么形式的曲线来拟合上述的数据，并预测使用4.5年后轿车的平均价格大致为多少？,旧车价格预测,2018/4/14,66,经济增长模型,3. 增加生产、发展经济的主要因素有增加投资、劳动力以及技术革新等，在研究国民经济产值与这些因素的数量关系时，由于技术水平不像资金、劳动力那样容易定量化，作为初步的模型，可认为技术水平不变，只讨论产值和资金、劳动力之间的关系。用Q，K，L分别表示产值、资金、劳动力，要寻求Q(K,L)。经过简化与分析，在经济学中，推导出一个著名的Cobb-Douglas生产函数：Q(K,L) = aKL， 0,1 （*）式中,，a要由经济统计数据确定。根据下表所给的统计数据，估计,，a的值。,2018/4/14,67,经济增长模型,2018/4/14,68,经济增长模型,4, 小行星绕着太阳运动，其轨迹是一条椭圆曲线。现测量小行星的位置为（5.764，0.6480),(6.286, 1.202),(6.759, 1.823),(7.168, 2.526),(7.408，3.360)。试利用这些数据确定曲线的方程和太阳的可能位置。（使用数据拟合方法）,2018/4/14,69,
