微分方程模型

微分方程模型第1页，共43页，2023年，2月20日，星期六第3讲微分方程建模方法及案例在研究实际问题时，常常会联系到某些变量的变化率或导数，这样所得到变量之间的关系式就是微分方程模型。微分方程模型反映的是变量之间的间接关系，因此，要得到直接关系，就得求微分方程。第2页，共43页，2023年，2月20日，星期六求解微分方程的方法求解微分方程有三种方法1）求精确解；2）求数值解（近似解）；3）定性理论方法。第3页，共43页，2023年，2月20日，星期六建立微分方程模型的方法（1）根据规律列方程利用数学、力学、物理、化学等学科中的定理或经过实验检验的规律等来建立微分方程模型。（2）微元分析法利用已知的定理与规律寻找微元之间的关系式，与第一种方法不同的是对微元而不是直接对函数及其导数应用规律。第4页，共43页，2023年，2月20日，星期六（3）模拟近似法在生物、经济等学科的实际问题中，许多现象的规律性不很清楚，即使有所了解也是极其复杂的，建模时在不同的假设下去模拟实际的现象，建立能近似反映问题的微分方程，然后从数学上求解或分析所建方程及其解的性质，再去同实际情况对比，检验此模型能否刻画、模拟某些实际现象。第5页，共43页，2023年，2月20日，星期六3.1年代鉴定问题在巴基斯坦一个洞穴里，发现了具有古代某种动物特征的骨碎片，科学家把它带到实验室，作碳14年代测定，分析表明，与的比例仅仅是活组织内的6.24%，能否判断此动物生活在多少年前？第6页，共43页，2023年，2月20日，星期六

年代测定：活体中的碳有一小部分是放射性同位素，这种放射性碳是由于宇宙射线在高层大气中的撞击引起的，经过一系列交换过程进入活组织内，直到在生物体内达到平衡浓度，这意味着在活体中，的数量与稳定的的数量成定比，生物体死亡后，交换过程就停止了，放射性碳便以每年八千分之一的速度减少。背景第7页，共43页，2023年，2月20日，星期六设t为死后年数，第8页，共43页，2023年，2月20日，星期六用Matlab求解对数log(x)x的自然对数log2(x)x的以2为底的对数log10(x)x的以10为底的对数以a为底x的对数，根据换底公式可表示为：log(x)/log(a)第9页，共43页，2023年，2月20日，星期六Matlab中的求解方法-8000*log(0.0624)ans=2.2194e+004第10页，共43页，2023年，2月20日，星期六年代测定的修订

1966年，耶鲁实验室的MinzeStuiver和加利福尼亚大学圣地亚哥分校的HansE.Suess在一份报告中指出：在2500到10000年前这段时间中测得的结果有差异，其根本原因在于那个年代，宇宙射线的放射性强度减弱了，偏差的峰值发生在大约6000年以前。他们提出了一个很成功的误差公式，用来校正根据碳测定出的2300年到6000年前这期间的年代：

真正的年代＝第11页，共43页，2023年，2月20日，星期六3.2人口演化模型第12页，共43页，2023年，2月20日，星期六预备知识:最小二乘原理当由实验提供了大量数据时,不能要求拟合函数在数据点处的偏差,即(i=1,2,…,m)

严格为零,但为了使近似曲线尽量反映所给数据点的变化趋势,需对偏差有所要求.通常要求偏差平方和最小,此即称为最小二乘原理第13页，共43页，2023年，2月20日，星期六Matlab的最优化工具箱中提供了求非线性最小二乘拟合的函数：lsqcurvefit，调用格式为：

x=lsqcurvefit(‘fun’,x0,xdata,ydata,…);fun是一个事先建立的定义函数F(x,xdata)的M-文件,自变量为x和xdata迭代初值已知数据点用Matlab作非线性最小二乘拟和注意:1.如无合理初值，那就只能给出一个猜想初值2.拟合结果是初值敏感的，因为找到的不一定是全局最优而可能是初值附近的局部最优.第14页，共43页，2023年，2月20日，星期六

例

用下面一组数据拟合中的参数a，b，k该问题即解最优化问题：第15页，共43页，2023年，2月20日，星期六

1）编写M-文件curvefun1.m

functionf=curvefun1(x,tdata)f=x(1)+x(2)*exp(0.02*x(3)*tdata)%其中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)解法第16页，共43页，2023年，2月20日，星期六问题的提出人口的增长情况是当前世界引起普遍关注的问题，早在18世纪人们就开始进行人口预报工作了。20世纪90年代，我们经常可以在报刊上看见关于人口增长的预报，说到20世纪末，或21世纪中叶，全世界（或某地区）的人口将达到多少多少亿。这些人口预报的数值是从哪里来的？准确不准确？你能不能对某地人口数目的演化进行一下估算？第17页，共43页，2023年，2月20日，星期六年（公元）人口（百万）17903.918005.318107.218209.6183012.9184017.1185023.2年（公元）人口（百万）186031.4187038.6188050.2189062.9190076.0191092.01920106.5年（公元）人口（百万）1930123.21940131.71950150.71960179.31970204.01980226.51990251.4下表是近两百年的美国人口统计数据，试依此建立美国人口增长的数学模型，最后用它预报2000年、2010年美国人口.第18页，共43页，2023年，2月20日，星期六问题的分析从整体来说，人口的变化由两个因素决定：出生和死亡。出生使得人口增加，死亡使得人口减少。对于局部地区来说，除了出生和死亡外，影响人口的变化还有两个因素：迁入和迁出。迁入使得局部人口增加，迁出使得局部人口减少。在迁入、迁出人口的差别不大时，人口的变化主要由出生率和死亡率决定。根据上面的分析，不难建立起人口演化模型。第19页，共43页，2023年，2月20日，星期六求解步骤[1]假设[2]建立模型[3]模型求解[4]模型的参数估计第20页，共43页，2023年，2月20日，星期六[1]

假设人口增长率r(即出生率b-死亡率d)是常数（或单位时间内人口的增长量与当时的人口成正比）.第21页，共43页，2023年，2月20日，星期六[2]

建立模型记时刻t=0时人口数为x0,时刻t的人口为x(t),由于量大,可视为连续、可微函数.t到时间内人口的增量为：于是x(t)满足微分方程：

(1)第22页，共43页，2023年，2月20日，星期六[3]

模型求解解微分方程得

(2)

(2)被称为指数增长模型——马尔萨斯提出(1798)，表示随着时间增加，人口按指数规律无限增长被称为等比数列模型第23页，共43页，2023年，2月20日，星期六[4]

模型的参数估计要用模型的结果来预报人口，必须对其中的参数r进行估计，这可以用表中所给的数据通过拟合得到.

第24页，共43页，2023年，2月20日，星期六拟合的关键问题1.t的变换指数函数exp(t)当t很大时可能会溢出,为了减小数据误差,首先将时间域变换至[0,20],所用的变换为:t=1800+(t-1800)/102.x0和r初值的确定x0初始值自然应取t=0时的x的值3.9r初始值取增长率的平均值mean(diff(x)./diff(t)./x(1:20))第25页，共43页，2023年，2月20日，星期六Matlab求解过程首先编写M文件fun.m，其中参数c(1)表示x0,c(2)表示r。functionx=fun(c,t)x=c(1)*exp(c(2)*t);第26页，共43页，2023年，2月20日，星期六再编写一个M文件clearall;closeall;t=0:1:20;x=[3.95.37.29.612.917.123.231.438.650.262.976.092.0106.5123.2131.7150.7179.3204.0226.5251.4];plot(t,x,'o');holdon;c(1)=3.9;c(2)=mean(diff(x)./diff(t)./x(1:20));k=lsqcurvefit(@fun,c,t,x)tt=[2122];xx=fun(k,tt);plot(tt,xx,'r*');tt=0:0.1:22;xx=fun(k,tt);plot(tt,xx,'r');holdoff第27页，共43页，2023年，2月20日，星期六holdon和holdoffholdon作用:当前轴及图形保持而不被刷新，准备接受此后的绘制例如，在作图时，想在一张图上同时显示多组数据(便于观察等等原因)的时候用holdon(就是等一下,我还想在本图上画个东西)holdoff自然就是取消这个功能第28页，共43页，2023年，2月20日，星期六指数增长模型的应用及局限性与19世纪以前欧洲一些地区人口统计数据吻合适用于19世纪后迁往加拿大的欧洲移民后代可用于短期人口增长预测不符合19世纪后多数地区人口增长规律不能预测较长期的人口增长过程19世纪后人口数据人口增长率r不是常数(逐渐下降)第29页，共43页，2023年，2月20日，星期六美国实际人口与按指数增长模型计算的人口比较年实际人口（百万）指数增长模型预测人口（百万）误差（%）17903.918005.318107.27.31.418209.610.04.2183012.913.76.2184017.118.79.4185023.225.610.3186031.435.010.8187038.647.823.8188050.265.530.5189062.989.642.4190076.0122.561.2191092.0167.682.11920106.5229.3115.3第30页，共43页，2023年，2月20日，星期六阻滞增长模型(Logistic模型)马尔萨斯模型为什么不能预测未来的人口呢？这主要是地球上的各种资源只能供一定数量的人生活,随着人口的增加,自然资源、环境条件等因素对人口增长的阻滞作用越来越显著如果当人口较少时,人口的自然增长率可以看作常数的话,那么当人口增加到一定数量以后,这个增长率就要随人口的增加而减小.因此,应对马尔萨斯模型中关于净增长率为常数的假设进行修改.r是x的减函数第31页，共43页，2023年，2月20日，星期六[1]

假设最简单假定（线性函数），r为固有增长率(x很小时).1838年,荷兰生物数学家韦尔侯斯特(Verhulst)引入常数,用来表示自然环境条件所能容许的最大人口数（一般说来,一个国家工业化程度越高,它的生活空间就越大,食物就越多,从而就越大）

第32页，共43页，2023年，2月20日，星期六[2]建立模型当时，增长率应为0代入得：

(3)

第33页，共43页，2023年，2月20日，星期六将（3）式代入（1）得：模型为：（4）

第34页，共43页，2023年，2月20日，星期六[3]

模型求解解方程组(4)得(5)

根据方程（4）作出曲线图，见下图，由该图可看出人口增长率随人口数的变化规律.

ox第35页，共43页，2023年，2月20日，星期六根据结果（5）作出x~t曲线，见下图，由该图可看出人口数随时间的变化规律.

S形曲线,x增加先快后慢第36页，共43页，2023年，2月20日，星期六也可将方程（4）离散化，得

t=0,1,2,…(6)第37页，共43页，2023年，2月20日，星期六[4]

模型的参数估计利用表中数据对r和xm拟合作业1!第38页，共43页，2023年，2月20日，星期六模型检验用模型计算2000年美国人口，与实际数据比较实际为281.4(百万)阻滞增长模型(Logistic模型)第39页，共43页，2023年，2月20日，星期六年实际人口阻滞增长模型公式（5）公式（6）预测人口误差（%）预测人口误差（%）17903.918005.35.90250.11373.90000.264218107.27.26140.00856.50740.096218209.68.93320.06958.68100.0957183012.910.98990.148111.41530.1151184017.113.52010.209415.12320.1156185023.216.63280.283119.81970.1457186031.420.46210.348326.52280.1553187038.625.17310.347835.45280.0815188050.230.96870.383143.53290.1328189062.938.09860.394356.18840.1067190076.046.86990.383370.14590.0770191092.057.66070.373384.73050.07901920106.570.93590.3339102.46260.03791930123.287.26740.2917118.95090.03451940131.71