数学实验第七次讲稿课件

数学实验之七数学实验之七数学实验之七数学实验之七数据拟合数据拟合数据拟合数据拟合2024/6/231实验目的实验目的11 了解最小二乘拟合的原理，掌握用了解最小二乘拟合的原理，掌握用 MATLABMATLAB作线性最小二乘拟合的方法。作线性最小二乘拟合的方法。22 通过实例学习如何用拟合方法解决通过实例学习如何用拟合方法解决 实际问题；通过实例理解参数辨识问题实际问题；通过实例理解参数辨识问题的几种方法。的几种方法。2024/6/232主要内容主要内容范例范例2 2：薄膜渗透率的测定：薄膜渗透率的测定布置实验布置实验引例引例1,1,引例引例2 2拟合的基本原理拟合的基本原理拟合函数的选取拟合函数的选取用用MATLABMATLAB作拟合计算作拟合计算范例范例1 1：静脉注射的给药方案：静脉注射的给药方案2024/6/233求电阻求电阻R R随温度随温度t t的变化规律的变化规律。已知热敏电阻数据：已知热敏电阻数据：温度温度t(0C)20.5 32.7 51.0 73.0 95.7电阻电阻R()765 826 873 942 1032引例引例1 1：热敏电阻电阻值的变化规律：热敏电阻电阻值的变化规律2024/6/234 设设 R=at+b a,b为待定系数为待定系数引例引例1 1：热敏电阻电阻值的变化规律：热敏电阻电阻值的变化规律2024/6/235引例引例2 2：血药浓度的变化规律：血药浓度的变化规律 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.24 3.01 已知一室模型快速静脉注射下已知一室模型快速静脉注射下的血药浓度数据的血药浓度数据(t=0(t=0注射注射300300mg)mg)求血药浓度随时间的变化规律求血药浓度随时间的变化规律c(t).2024/6/236作图观察半对数坐标系半对数坐标系(semilogy)下的图形下的图形Log10c(t)=a t+b2024/6/237主要内容主要内容范例范例2 2：薄膜渗透率的测定：薄膜渗透率的测定布置实验布置实验引例引例1,1,引例引例2 2拟合的基本原理拟合的基本原理拟合函数的选取拟合函数的选取用用MATLABMATLAB作拟合计算作拟合计算范例范例1 1：静脉注射的给药方案：静脉注射的给药方案2024/6/238曲线拟合的基本原理曲线拟合的基本原理问题的提法问题的提法已知一组（二维）数据，即平面上已知一组（二维）数据，即平面上 n个点个点（xi,yi)i=1,n,寻求一个函数（曲线）寻求一个函数（曲线）y=f(x),使使 f(x)在某种准则下与所有数据点最为接近，在某种准则下与所有数据点最为接近，即曲线拟合得最好。即曲线拟合得最好。+xyy=f(x)(xi,yi)i i 为点为点（xi,yi)与与曲线曲线 y=f(x)的距离的距离2024/6/239问题的数学模型问题的数学模型 步骤：步骤：步骤：步骤：1)1)1)1)选定一类函数选定一类函数选定一类函数选定一类函数 f(xf(xf(xf(x，a a a a1 1 1 1，a a a a2 2 2 2，a a a am m m m)（1 1 1 1）其中其中其中其中 a a a a1 1 1 1,a,a,a,a2 2 2 2,a a a am m m m 为待定常数。为待定常数。为待定常数。为待定常数。2)2)2)2)确定参数确定参数确定参数确定参数a a a a1 1 1 1,a,a,a,a2 2 2 2,a a a am m m m，准则准则准则准则(最小二乘准则）：使最小二乘准则）：使最小二乘准则）：使最小二乘准则）：使n n n n个点（个点（个点（个点（x x x xi i i i,y,y,y,yi i i i)与曲线与曲线与曲线与曲线 y=f(x y=f(x y=f(x y=f(x，a a a a1 1 1 1，a a a a2 2 2 2，a a a am m m m)的距离的距离的距离的距离 i i i i 的平方和最小的平方和最小的平方和最小的平方和最小 。2024/6/2310记记问题归结为问题归结为:求求 a1,a2,am 使使 J(a1,a2,am)最小。最小。这样的拟合称为这样的拟合称为最小二乘拟合最小二乘拟合。问题的数学模型问题的数学模型 2024/6/2311 特别，特别，若选定一组函数若选定一组函数 r1(x),r2(x),rm(x),m0);2.2.血液容积血液容积v,t=0瞬时注射剂量瞬时注射剂量d,血药血药浓度立即为浓度立即为d/v.模型假设模型假设：由假设1，3.3.快速静脉注射下一室模型的血药浓度：快速静脉注射下一室模型的血药浓度：c(t)c(t)ctc00由假设2，2024/6/2328若若c1=10,c2=25(g/ml),给药方案设计归结为根据数据给药方案设计归结为根据数据(ti,ci)i=1,n(d 给定）给定）拟合曲线拟合曲线c(t),以确定系数以确定系数k,v.cc2c10t给药方案给药方案 设计：设计：D0:初次剂量；初次剂量；:注射时间间隔；注射时间间隔；D:重复注射剂量。重复注射剂量。范例范例1：静脉注射的给药方案：静脉注射的给药方案 2024/6/2329给药方案给药方案 血药浓度数据血药浓度数据(t=0注射注射300mg)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.24 3.01记记则则思考思考:取对数化为线性最小二乘取对数化为线性最小二乘,对结果有什么影响对结果有什么影响？2024/6/2330主要内容主要内容范例范例2 2：薄膜渗透率的测定：薄膜渗透率的测定布置实验布置实验引例引例1,1,引例引例2 2拟合的基本原理拟合的基本原理拟合函数的选取拟合函数的选取用用MATLABMATLAB作拟合计算作拟合计算范例范例1 1：静脉注射的给药方案：静脉注射的给药方案2024/6/2331范例范例2 2：薄膜渗透率的测定：薄膜渗透率的测定1问题背景；2假设；3问题的分析 4数学模型；5参数辨识方法：函数拟合法 非线性规划法；导函数拟合法；线性化迭代法。6结果分析VAVBS2024/6/2332 某种医用薄膜在试制时需测定其被物质分某种医用薄膜在试制时需测定其被物质分子穿透的能力。子穿透的能力。测定方法：用面积为测定方法：用面积为S S的薄膜将容器分成的薄膜将容器分成两部份，在两部分中分别注满该物质的两种不两部份，在两部分中分别注满该物质的两种不同浓度的溶液。此时该物质分子就会从一侧向同浓度的溶液。此时该物质分子就会从一侧向另一侧扩散。平均每单位时间通过单位面积薄另一侧扩散。平均每单位时间通过单位面积薄膜的物质分子量与膜两侧溶液的浓度差成正比，膜的物质分子量与膜两侧溶液的浓度差成正比，比例系数比例系数K K称为渗透率。定时测量容器中薄膜称为渗透率。定时测量容器中薄膜某一侧的溶液浓度，以此确定某一侧的溶液浓度，以此确定K K。VAVBS问题背景问题背景薄膜渗透率的测定薄膜渗透率的测定2024/6/23333 3）薄膜是双向同性的即物质从膜的任何一）薄膜是双向同性的即物质从膜的任何一侧向另一侧渗透的性能是相同的。侧向另一侧渗透的性能是相同的。VAVBS假设：假设：1）薄膜两侧的溶液始终是均匀的；薄膜两侧的溶液始终是均匀的；2）平均每单位时间通过单位面积薄膜的物平均每单位时间通过单位面积薄膜的物质分子量与膜两侧溶液的浓度差成正比。质分子量与膜两侧溶液的浓度差成正比。2024/6/2334VAVBS 第一步：通过机理分析确定容第一步：通过机理分析确定容器一侧的浓度随时间的变化规律：器一侧的浓度随时间的变化规律：C=CB(t)第二步：利用实验数据第二步：利用实验数据(tj,CB(tj)，和函数拟合的方法求出其中的未知，和函数拟合的方法求出其中的未知参数，包括渗透率参数，包括渗透率K.解决问题的思路：解决问题的思路：2024/6/2335 考考察察时时段段tt，t+tt+t薄薄膜膜的的一一侧侧容容器器中中该该物物质质量的变化。质质量的变化。1）以容器以容器A A侧为例，在时段侧为例，在时段tt，t+tt+t物质质量的增量为：物质质量的增量为：问题分析：问题分析：由假设由假设2,2,在时段在时段tt，t+tt+t，从，从B B侧渗透侧渗透至至A A侧的该物质的质量为：侧的该物质的质量为：VAVBS2024/6/2336于是有：于是有：两边除以两边除以tt，并令，并令t0t0取极限再稍加整理取极限再稍加整理即得：即得：（1）问题分析：问题分析：2024/6/2337分别表示在初始时刻两侧溶分别表示在初始时刻两侧溶其中其中2)注意到任意时刻，整个溶液中含有该注意到任意时刻，整个溶液中含有该 物质的质量物质的质量,与初始时刻该物质的含量相与初始时刻该物质的含量相同，因此同，因此 问题分析：问题分析：液的浓度液的浓度2024/6/2338从而：从而：加上初值条件：加上初值条件：代入式（代入式（1）得：）得：便可得出便可得出CB(t)的变化规律，从而根据实验数据进行的变化规律，从而根据实验数据进行拟合，估计出参数拟合，估计出参数K,。问题分析：问题分析：2024/6/2339基于假设和前面的分析，基于假设和前面的分析，B B侧的浓度侧的浓度CB(t)应满足如下微分方程和初始条件：应满足如下微分方程和初始条件：数学模型数学模型2024/6/2340模型求解方法模型求解方法 前面得到的模型是一个带初值的一阶前面得到的模型是一个带初值的一阶线性微分方程，解之得：线性微分方程，解之得：问题归结为利用问题归结为利用C CB B在时刻在时刻t tj j的测量数据的测量数据C Cj j(j=1,2,.,N)(j=1,2,.,N)来辨识来辨识 K K 和和 。1.1.函数拟合函数拟合法法2024/6/2341引入引入从而从而 用函数用函数CB(t)来拟合所给的实验数据，来拟合所给的实验数据，从而估计出其中的参数从而估计出其中的参数a,b,K。若若将其代入上式有：将其代入上式有：拟合函数化简拟合函数化简:2024/6/2342用用MATLABMATLAB软件进行计算软件进行计算.1 1）编写函数编写函数M-M-文件文件 nongdu.mnongdu.mfunction f=nongdu(x,tdata)function f=nongdu(x,tdata)f=x(1)+x(2)*exp(-0.02*x(3)*tdata);f=x(1)+x(2)*exp(-0.02*x(3)*tdata);其中其中 x(1)=a;x(2)=b;x(3)=k;x(1)=a;x(2)=b;x(3)=k;2)2)编写编写M M文件文件 (baomo.m)(baomo.m)tdata=linspace(100,1000,10);tdata=linspace(100,1000,10);cdata=4.54 4.99 5.35 5.65 5.90 6.10.cdata=4.54 4.99 5.35 5.65 5.90 6.10.6.26 6.39 6.50 6.59;6.26 6.39 6.50 6.59;x0=0.2,0.05,0.05;x0=0.2,0.05,0.05;x x=lsqcurvefit(lsqcurvefit(nongdunongdu,x0,tdata,cdata),x0,tdata,cdata)编写编写MATLABMATLAB程序程序:2024/6/23433)3)输出结果输出结果:x=0.007 -0.003 0.1012 即即 k=0.1012,a=0.007,b=-0.003,结果结果:进一步求得：进一步求得：To MATLAB(nongdu,baomo1)2024/6/23442.非线性规划法非线性规划法 利用利用C CB B在时刻在时刻t tj j的测量数据的测量数据C Cj j(j=1,2,.,N)(j=1,2,.,N)来来辨识辨识 K K 和和 。问题可转化为求函数问题可转化为求函数即求函数即求函数的最小值点（的最小值点（K K，a a，b b）。）。2024/6/23453.3.导函数拟合法导函数拟合法令令上式变为：上式变为：这可以看作这可以看作随随CB的变化规律的变化规律(j=1,2,.,N)若知道一组数据若知道一组数据则可用最小二乘拟合的方法来求出函数则可用最小二乘拟合的方法来求出函数中的未知参数中的未知参数K和和h。2024/6/2346即为求参数即为求参数K,a使下列误差函数达到最小：使下列误差函数达到最小：该问题等价于用函该问题等价于用函 数数 f(K,a,CB)=K(0.01a-0.02CB)来拟合数据来拟合数据(j=1,2,.,N)用上述方法求出参数K,a.3.3.导函数拟合法导函数拟合法2024/6/23474.4.线性化迭代法线性化迭代法前面带初始条件的一阶线性微分方程的解为前面带初始条件的一阶线性微分方程的解为其中：其中：如果得到了参数如果得到了参数K的一个较好的近似值的一个较好的近似值K*，则，则将将CB(t)关于关于K在在K*处展开，略去处展开，略去K的二次及以上的二次及以上的项得的项得CB(t)的一个近似式的一个近似式2024/6/23484.4.线性化迭代法线性化迭代法通过极小化确定a,b,d,再由K=d/0.02b得到K*的修正值K。K*K*-K,得到K的一个新的近似值，用同样的方法再求新的修正值K。这个过程可以不断重复，直到修正值足够小为止。2024/6/23491）当K的初值取为k=0.3时，出现奇异情况，迭代不收敛；2）当K的初值取为k=0.2时，经四次迭代，已经收敛到一个很好的解。迭代结果如下表。4.4.线性化迭代法线性化迭代法2024/6/23503）取K的初值为0.1009,只一次迭代就得到2）中的最后结果。4.4.线性化迭代法线性化迭代法导函数拟合法得出的参数值精度有限，线性化迭代法要求参数的初值比较接近精确值。可将导函数拟合法和线性化迭代法结合起来，把前者得到的参数K的值作为迭代法中K的初值，可使其收敛或收敛更快。2024/6/2351结果及误差分析结果及误差分析 几种方法得出的结果及相应的误差总结于几种方法得出的结果及相应的误差总结于下表，误差为计算数据与实验数据之差的平方下表，误差为计算数据与实验数据之差的平方和。和。2024/6/2352作出拟合函数曲线和数据点的图形作出拟合函数曲线和数据点的图形tdata=linspace(100,1000,10);cdata=1e-05.*454 499 535 565 590.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)函数拟合法的误差分析程序2024/6/2353函数拟合法的拟合效果函数拟合法的拟合效果返回返回2024/6/2354人口问题人口问题马尔萨斯人口模型马尔萨斯人口模型假设:人口出生率和死亡率是常数，因此人口的净增长率为常数.设时刻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）初始数据能不能也作为参数处理？这样会得到更好的结果吗？马尔萨斯人口模型马尔萨斯人口模型2024/6/2359处理方法方法1：直接使用原始数据（t和p），使用非线性拟合求解方程的系数，并和原始数据作比较。方法2：因为 ，也就是lnp(t)=k(t-t0)+lnp0。记q(t)=lnp(t),q0=lnp0。上面的问题变成一个线性回归问题，可以得到相应的系数，和上面的方法已经原始数据比较。方法3：还是使用方法1，仅仅对时间做一个变换，比如17900,18001，等等。重新使用方法1进行计算。1837 1837年荷兰生物数学家年荷兰生物数学家Verhulst Verhulst 提出改进提出改进.人口增长率不应该是常数，假人口增长率不应该是常数，假设增长率设增长率k k是随着人口数量接近最大数量是随着人口数量接近最大数量M M而而线性递减线性递减.逻辑斯谛人口模型逻辑斯谛人口模型 从而得到改进后的人口模型从而得到改进后的人口模型(逻辑斯谛增逻辑斯谛增长模型长模型)逻辑斯谛人口模型逻辑斯谛人口模型可以求得，该人口模型的解为：可以求得，该人口模型的解为：在这个方程中，将M，r，p0作为参数，使用前面类似的方法计算。并比较拟合数据和原始数据。1.1.掌握用掌握用MATLABMATLAB作最小二乘多项式拟作最小二乘多项式拟合和曲线拟合的方法；合和曲线拟合的方法；2.2.通过实例学习如何用拟合方法解决通过实例学习如何用拟合方法解决实际问题，注意与插值方法的区别。实际问题，注意与插值方法的区别。3.3.鼓励不囿于固定的模式或秩序，灵鼓励不囿于固定的模式或秩序，灵 活调整思路，突破思维的呆板性，活调整思路，突破思维的呆板性，找到打破常规的解法。并在文献检找到打破常规的解法。并在文献检索、动手和动脑等方面得到锻炼，索、动手和动脑等方面得到锻炼，树立创新意识。树立创新意识。布置实验布置实验实实验验目目的的2024/6/2362实验内容实验内容1.1.下表给出了近两个世纪某国人口的统计下表给出了近两个世纪某国人口的统计数据，试预测数据，试预测20102010年该国的人口年该国的人口。人口预测问题人口预测问题2024/6/23631)1)可先用以上数据拟合可先用以上数据拟合MulthusMulthus人口指数增长模人口指数增长模型，根据检验结果进一步讨论马尔萨斯人口模型型，根据检验结果进一步讨论马尔萨斯人口模型的改进。的改进。2)2)Malthus Malthus 模模型型的的基基本本假假设设是是：人人口口的的增增长长率率为为常常数数，记记为为r r。记记时时刻刻t t的的人人口口为为x(tx(t），且且初初始始时时刻的人口为刻的人口为x x0 0，于是得到如下微分方程：，于是得到如下微分方程：2024/6/23643)Malthus 3)Malthus 模型导致人口总数的剧烈上升，为此模型导致人口总数的剧烈上升，为此对该模型做一些改进。其中一种改进是认为增长对该模型做一些改进。其中一种改进是认为增长率应该是率应该是x x的函数，并且随着的函数，并且随着x x增大出现资源的竞增大出现资源的竞争导致增长率减小。因此一般可以假设争导致增长率减小。因此一般可以假设r(x)=r-r(x)=r-x/kx/k。得到如下微分方程：。得到如下微分方程：2024/6/2365实验内容实验内容2.2.某年美国旧车价格的调查资料如下表其中某年美国旧车价格的调查资料如下表其中x xi i表表示轿车的使用年数，示轿车的使用年数，y yi i表示相应的平均价格。试表示相应的平均价格。试分析用什么形式的曲线来拟合上述的数据，并预分析用什么形式的曲线来拟合上述的数据，并预测使用测使用4.54.5年后轿车的平均价格大致为多少？年后轿车的平均价格大致为多少？旧车价格预测旧车价格预测xi12345678910yi26151943149410877655384842902262042024/6/2366 经济增长模型经济增长模型3.增加生产、发展经济的主要因素有增加投资、劳动力以及技术革新等，在研究国民经济产值与这些因素的数量关系时，由于技术水平不像资金、劳动力那样容易定量化，作为初步的模型，可认为技术水平不变，只讨论产值和资金、劳动力之间的关系。用Q，K，L分别表示产值、资金、劳动力，要寻求Q(K,L)。经过简化与分析，在经济学中，推导出一个著名的Cobb-Douglas生产函数：Q(K,L)=aKL，0,1 （*）式中,，a要由经济统计数据确定。根据下表所给的统计数据，估计,，a的值。2024/6/2367 经济增长模型经济增长模型2024/6/2368 t Q K L t Q K L1900 1.05 1.04 1.051901 1.18 1.06 1.081902 1.29 1.16 1.181903 1.30 1.22 1.221904 1.30 1.27 1.171905 1.42 1.37 1.301906 1.50 1.44 1.391907 1.52 1.53 1.471908 1.46 1.57 1.311909 1.60 2.05 1.431910 1.69 2.51 1.581911 1.81 2.63 1.591912 1.93 2.74 1.661913 1.95 2.82 1.68 1914 2.01 3.24 1.65 1915 2.00 3.24 1.62 1916 2.09 3.61 1.86 1917 1.96 4.10 1.93 1918 2.20 4.36 1.96 1919 2.12 4.77 1.95 1920 2.16 4.75 1.90 1921 2.08 4.54 1.58 1922 2.24 4.54 1.67 1923 2.56 4.58 1.82 1924 2.34 4.58 1.60 1925 2.45 4.58 1.61 1926 2.58 4.54 1.64 经济增长模型经济增长模型4,小行星绕着太阳运动，其轨迹是一条椭圆曲线。现测量小行星的位置为（5.764，0.6480),(6.286,1.202),(6.759,1.823),(7.168,2.526),(7.408，3.360)。试利用这些数据确定曲线的方程和太阳的可能位置。（使用数据拟合方法）2024/6/2369