数学实验第七次讲稿

数学实验之七数据拟合2023/2/21实验目的

[1]

了解最小二乘拟合的原理，掌握用MATLAB作线性最小二乘拟合的方法。[2]

通过实例学习如何用拟合方法解决实际问题；通过实例理解参数辨识问题的几种方法。2023/2/22主要内容范例2：薄膜渗透率的测定布置实验引例1,引例2拟合的基本原理拟合函数的选取用MATLAB作拟合计算范例1：静脉注射的给药方案2023/2/23求电阻R随温度t的变化规律。已知热敏电阻数据：温度t(0C)20.532.751.073.095.7电阻R()7658268739421032引例1：热敏电阻电阻值的变化规律2023/2/24

设

R=at+ba,b为待定系数引例1：热敏电阻电阻值的变化规律2023/2/25引例2：血药浓度的变化规律

t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01

已知一室模型快速静脉注射下的血药浓度数据(t=0注射300mg)求血药浓度随时间的变化规律c(t).2023/2/26作图观察半对数坐标系(semilogy)下的图形Log10c(t)=at+b2023/2/27主要内容范例2：薄膜渗透率的测定布置实验引例1,引例2拟合的基本原理拟合函数的选取用MATLAB作拟合计算范例1：静脉注射的给药方案2023/2/28曲线拟合的基本原理——问题的提法已知一组（二维）数据，即平面上n个点（xi,yi)i=1,…n,

寻求一个函数（曲线）y=f(x),

使f(x)

在某种准则下与所有数据点最为接近，即曲线拟合得最好。

+++++++++xyy=f(x)(xi,yi)ii为点（xi,yi)与曲线y=f(x)的距离2023/2/29问题的数学模型步骤：1)选定一类函数

f(x，a1，a2，

…，am)（1）其中

a1,a2,…am

为待定常数。2)确定参数a1,a2,…am，准则(最小二乘准则）：使n个点（xi,yi)

与曲线y=f(x，a1，a2，

…，am)的距离

i的平方和最小

。2023/2/210记问题归结为:求

a1,a2,…am

使

J(a1,a2,…am)

最小。这样的拟合称为最小二乘拟合。问题的数学模型2023/2/211

特别，若选定一组函数

r1(x),r2(x),…rm(x),m<n,令

f(x,a1,a2,…am

)=a1r1(x)+a2r2(x)+…+amrm(x)其中

a1,a2,…am

为待定系数。问题的解法——线性最小二乘法

当f具有上述形式，且按最小二乘准则确定待定系数时，称这样的拟合为线性最小二乘拟合，称对应的函数f(x,a1,a2,…am

)

为已知数据点的线性最小二乘拟合函数。2023/2/2121)除了最小二乘准则（即各点误差的平方和最小），你认为还可以用怎样的拟合准则？

2)比较起来，最小二乘准则有什么优点？2023/2/213主要内容范例2：薄膜渗透率的测定布置实验引例1,引例2拟合的基本原理拟合函数的选取用MATLAB作拟合计算范例1：静脉注射的给药方案2023/2/214最小二乘拟合函数的选取+++++++++++++++f=a1+a2x将数据(xi,yi)i=1,…n作图，通

过直观判断确定f：f=a1+a2x+a3x2f=a1+a2x+a3x22023/2/215最小二乘拟合函数的选取2.通过机理分析建立数学模型来确定f。+++++++++++++++f=a1+a2/xf=a1exp(a2x)f=a1exp(a2x)2023/2/216主要内容范例2：薄膜渗透率的测定布置实验引例1,引例2拟合的基本原理拟合函数的选取用MATLAB作拟合计算范例1：静脉注射的给药方案2023/2/217多项式拟合：作多项式f(x)=a1xm+…+amx+am+1函数拟合,可利用已有程序polyfit,其调用格式为:a=polyfit(x,y,m)用MATLAB作最小二乘拟合数据点拟合多项式次数系数2023/2/218例1.由数据温度t(0C)20.532.751.073.095.7电阻R()7658268739421032拟合R=f(t)用MATLAB作多项式拟合2023/2/2192.用命令polyfit(x,y,m)作最小二乘拟合3.编写MATLAB程序dianzu1.m，并运行得到：a1=3.3940,a2=702.49181.选取拟合函数R=

a1t+a2用MATLAB作多项式拟合2023/2/220t=[20.532.5517395.7];r=[7658268739421032];aa=polyfit(t,r,1);a=aa(1)b=aa(2)y=polyval(aa,t);plot(t,r,'k+',t,y,'r')用MATLAB作多项式拟合2023/2/2212.曲线拟合：

作一般的最小二乘曲线拟合，可利用已有程序lsqcurvefit,其调用格式为：

[a,resnorm,residual]=lsqcurvefit(‘f’,a0,x,y)

用MATLAB作曲线拟合数据点待定常数a的初值函数M文件拟合函数y=f(a,x)的函数M—文件误差平方和各数据点处的误差2023/2/222用MATLAB作最小二乘曲线拟合xdata=[0:.1:2]ydata=[5.89553.56392.51731.97901.89901.39381.13591.00961.03430.84350.68560.61000.53920.39460.39030.54740.34590.13700.22110.17040.2636]例2：用函数

y(x)=a1*exp(-a2*x)+a3*exp(-a4*x)拟合下列数据点：1.用命令lsqcurvefit(‘fun’,a0,x,y)2023/2/223a0=[1,1,1,0];xdata=[0:.1:2];%ydata省略;[a,resnorm,residual,flag,output]=lsqcurvefit('fitf',a0,xdata,ydata)xi=linspace(0,2,200);yi=fitf(a,xi);plot(xdata,ydata,'ro',xi,yi)xlabel('x'),ylabel('y=f(x)'),title('nonlinearcurvefitting')用MATLAB作最小二乘曲线拟合2023/2/224functiony=fitf(x,xdata)y=x(1)*exp(-x(2)*xdata)+x(3)*exp(-x(4)*xdata);主要内容范例2：薄膜渗透率的测定布置实验引例1,引例2拟合的基本原理拟合函数的选取用MATLAB作拟合计算范例1：静脉注射的给药方案2023/2/226范例1：静脉注射的给药方案问题背景

一种新药用于临床之前，必须设计给药方案。在快速静脉注射下，所谓给药方案是指，每次注射量多大，间隔时间多长。药物进入肌体后随血液输送到全身，在这过程中不断被吸收、分解、代谢，最终排出体外。药物向体外排出的速率与血药浓度成正比。单位体积血液中的药物含量，称血药浓度。临床上，每种药物有一个最小有效浓度c1和最大治疗浓度c2。2023/2/227范例1：静脉注射的给药方案1.药物排除速率与血药浓度成正比，比例系数

k(>0);2.血液容积v,t=0瞬时注射剂量d,血药浓度立即为d/v.模型假设：由假设1，3.快速静脉注射下一室模型的血药浓度：c(t)ctc00由假设2，2023/2/228若c1=10,c2=25(g/ml),

给药方案设计归结为根据数据(ti,ci)i=1,…n(d给定）拟合曲线c(t),以确定系数k,v.cc2c10t给药方案设计：D0:初次剂量；:注射时间间隔；D:重复注射剂量。范例1：静脉注射的给药方案2023/2/229给药方案血药浓度数据(t=0注射300mg)

t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01记则思考:取对数化为线性最小二乘,对结果有什么影响？2023/2/230主要内容范例2：薄膜渗透率的测定布置实验引例1,引例2拟合的基本原理拟合函数的选取用MATLAB作拟合计算范例1：静脉注射的给药方案2023/2/231范例2：薄膜渗透率的测定[1]问题背景；[2]假设；[3]问题的分析[4]数学模型；[5]参数辨识方法：

函数拟合法

非线性规划法；

导函数拟合法；线性化迭代法。[6]结果分析VAVBS2023/2/232

某种医用薄膜在试制时需测定其被物质分子穿透的能力。测定方法：用面积为S的薄膜将容器分成两部份，在两部分中分别注满该物质的两种不同浓度的溶液。此时该物质分子就会从一侧向另一侧扩散。平均每单位时间通过单位面积薄膜的物质分子量与膜两侧溶液的浓度差成正比，比例系数K称为渗透率。定时测量容器中薄膜某一侧的溶液浓度，以此确定K。VAVBS问题背景薄膜渗透率的测定2023/2/2333）薄膜是双向同性的即物质从膜的任何一侧向另一侧渗透的性能是相同的。VAVBS假设：1）薄膜两侧的溶液始终是均匀的；2）平均每单位时间通过单位面积薄膜的物质分子量与膜两侧溶液的浓度差成正比。2023/2/234VAVBS

第一步：通过机理分析确定容器一侧的浓度随时间的变化规律：

C=CB(t)

第二步：利用实验数据(tj,CB(tj))，和函数拟合的方法求出其中的未知参数，包括渗透率K.解决问题的思路：2023/2/235

考察时段[t，t+Δt]薄膜的一侧容器中该物质质量的变化。

1）以容器A侧为例，在时段[t，t+Δt]物质质量的增量为：

问题分析：

由假设2,在时段[t，t+Δt]，从B侧渗透至A侧的该物质的质量为：VAVBS2023/2/236于是有：两边除以Δt，并令Δt→0取极限再稍加整理即得：（1）问题分析：2023/2/237分别表示在初始时刻两侧溶其中2)注意到任意时刻，整个溶液中含有该物质的质量,与初始时刻该物质的含量相同，因此问题分析：液的浓度2023/2/238从而：加上初值条件：代入式（1）得：便可得出CB(t)的变化规律，从而根据实验数据进行拟合，估计出参数K,。问题分析：2023/2/239基于假设和前面的分析，B侧的浓度CB(t)应满足如下微分方程和初始条件：数学模型2023/2/240模型求解方法

前面得到的模型是一个带初值的一阶线性微分方程，解之得：

问题归结为利用CB在时刻tj的测量数据Cj(j=1,2,...,N)来辨识K和。1.函数拟合法2023/2/241引入从而

用函数CB(t)来拟合所给的实验数据，从而估计出其中的参数a,b,K。若将其代入上式有：拟合函数化简:2023/2/242用MATLAB软件进行计算.1）编写函数M-文件nongdu.mfunctionf=nongdu(x,tdata)f=x(1)+x(2)*exp(-0.02*x(3)*tdata);其中x(1)=a;x(2)=b;x(3)=k;2)编写M文件(baomo.m)tdata=linspace(100,1000,10);cdata=[4.544.995.355.655.906.10...6.266.396.506.59];x0=[0.2,0.05,0.05];x=lsqcurvefit(‘nongdu’,x0,tdata,cdata)编写MATLAB程序:2023/2/2433)输出结果:x=0.007-0.0030.1012

即k=0.1012,a=0.007,b=-0.003,结果:进一步求得：ToMATLAB(nongdu,baomo1)2023/2/2442.非线性规划法

利用CB在时刻tj的测量数据Cj(j=1,2,...,N)来辨识K和。问题可转化为求函数即求函数的最小值点（K，a，b）。2023/2/2453.导函数拟合法令上式变为：这可以看作随CB的变化规律(j=1,2,...,N)若知道一组数据则可用最小二乘拟合的方法来求出函数中的未知参数K和h。2023/2/246即为求参数K,a使下列误差函数达到最小：

该问题等价于用函数f(K,a,CB)=K(0.01a-0.02CB)来拟合数据(j=1,2,...,N)用上述方法求出参数K,a.‡3.导函数拟合法2023/2/2474.线性化迭代法前面带初始条件的一阶线性微分方程的解为其中：

如果得到了参数K的一个较好的近似值K*，则将CB(t)关于K在K*处展开，略去K的二次及以上的项得CB(t)的一个近似式2023/2/2484.线性化迭代法通过极小化确定a,b,d,再由K=d/0.02b得到K*的修正值K。K*K*-K,得到K的一个新的近似值，用同样的方法再求新的修正值K。这个过程可以不断重复，直到修正值足够小为止。2023/2/2491）当K的初值取为k=0.3时，出现奇异情况，迭代不收敛；2）当K的初值取为k=0.2时，经四次迭代，已经收敛到一个很好的解。迭代结果如下表。4.线性化迭代法2023/2/2503）取K的初值为0.1009,只一次迭代就得到2）中的最后结果。4.线性化迭代法导函数拟合法得出的参数值精度有限，线性化迭代法要求参数的初值比较接近精确值。可将导函数拟合法和线性化迭代法结合起来，把前者得到的参数K的值作为迭代法中K的初值，可使其收敛或收敛更快。2023/2/251结果及误差分析

几种方法得出的结果及相应的误差总结于下表，误差为计算数据与实验数据之差的平方和。2023/2/252％作出拟合函数曲线和数据点的图形tdata=linspace(100,1000,10);cdata=1e-05.*[454499535565590...

610626639650659];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)函数拟合法的误差分析程序2023/2/253函数拟合法的拟合效果返回2023/2/254人口问题马尔萨斯人口模型假设:人口出生率和死亡率是常数，因此人口的净增长率为常数.设时刻t的人口数量为p(t)，人口出生率为b，死亡率为d,(k=b-d)因此马尔萨斯人口模型如下：该微分方程初值问题的解析解：马尔萨斯人口模型为了验证，采用某国家的人口历史数据,见下页表.年份17901800181018201830184018501860187018801890人口/103392953087240963812866170692319231433385585015662948年份19001910192019301940195019601970198019902000人口/1037599591972105711122755131669150697179323203212226505248710281416首先确定参数k：因为所以利用1790年和1840年的数据计算得出：=0.0294.马尔萨斯人口模型由此预测，1850年的人口数量为22,898,000，误差为1%，1900年的人口数量为99,476,000，误差为31%，2000年的人口数量为1,877,463,000，误差为567%，2050年将达到80多个亿？问题：（1）能不能用尽可能多的数据拟合微分方程的参数？（2）初始数据能不能也作为参数处理？这样会得到更好的结果吗？马尔萨斯人口模型2023/2/259处理方法方法1：直接使用原始数据（t和p），使用非线性拟合求解方程的系数，并和原始数据作比较。方法2：因为，也就是lnp(t)=k(t-t0)+lnp0。记q(t)=lnp(t),q0=lnp0。上面的问题变成一个线性回归问题，可以得到相应的系数，和上面的方法已经原始数据比较。方法3：还是使用方法1，仅仅对时间做一个变换，比如1790—0,1800—1，等等。重新使用方法1进行计算。1837年荷兰生物数学家Verhulst提出改进.人口增长率不应该是常数，假设增长率k是随着人口数量接近最大数量M而线性递减.逻辑斯谛人口模型

从而得到改进后的人口模型(逻辑斯谛增长模型)逻辑斯谛人口模型可以求得，该人口模型的解为：在这个方程中，将M，r，p0作为参数，使用前面类似的方法计算。并比较拟合数据和原始数据。掌握用MATLAB作最小二乘多项式拟合和曲线拟合的方法；2.通过实例学习如何用拟合方法解决实际问题，注意与插值方法的区别。3.鼓励不囿于固定的模式或秩序，灵活调整思路，突破思维的呆板性，找到打破常规的解法。并在文献检索、动手和动脑等方面得到锻炼，树立创新意识。‡布置实验实验目的2023/2/262实验内容‡1.下表给出了近两个世纪某国人口的统计数据，试预测2010年该国的人口。人口预测问题2023/2/2631)可先用以上数据拟合Multhus人口指数增长模型，根据检验结果进一步讨论马尔萨斯人口模型的改进。2)Malthus模型的基本假设是：人口的增长率为常数，记为r。记时刻t的人口为x(t），且初始时刻的人口为x0，于是得到如下微分方程：2023/2/2643)Malthus模型导致人口总数的剧烈上升，为此对该模型做一些改进。其中一种改进是认为增长率应该是x的函数，并且随着x增大出现资源的竞争导致增长率减小。因此一般可以假设r(x)=r-x/k。得到如下微分方程：2023/2/265实验内容‡2.某年美国旧车价格的调查资料如下表其中xi表示轿车的使用年数，yi表示相应的平均价格。试分析用什么形式的曲线来拟合上述的数据，并预测使用4.5年后轿车的平均价格大致为多少？旧车价格预测xi12345678910yi26151943149410877655384842902262042023/2/266‡

经济增长模型3.增加生产、发展经济的主要因素有增加投资、劳动力以及技术革新等，在研究国民经济产值与这些因素的数量关系时，由于技术水平不像资金、劳动力那样容易定量化，作为初步的模型，可认为技术水平不变，只讨论产值和资金、劳动力之间的关系。用Q，K，L分别表示产值、资金、劳动力，要寻求Q(K,L)。经过简化与分析，在经济学中，推导出一个著名的Cobb-Douglas生产函数：Q(K,L)=aKαLβ，0<α,β<1（*）式中α,β，a要由经济统计数据确定。根据下表所给的统计数据，估计α,β，a的值。

2023/2/267‡

经济增长模型2023/2/268

tQK