微分方程建模学习.doc

1、1微分方程建模一般说来，微分方程建模的方法大致可以分为以下的几个步骤：1.根据实际问题的要求确定要研究的量，包括自变量、未知函数、必要的参数等以及它们各自的变化区间；2.列方程。可以在合理假设的前提下，利用导数表示斜率、速度、变化率的实际意义，根据一些基本定理（几何的、物理的、化学的或生物学的等等）或规律，找出未知函数的导数（或微分）与相关各量之间的等量关系式，建立微分方程并确定定解条件（注：如果没有现成的定理可供利用，也可以用微元分析法与模拟近似法列出微分方程） ；3.解微分方程；4.对模型的适用性作出评价，即用已知的数据检验微分方程的解是否与实际相符。若结果与实际存在一定的差距，则还要对方

2、程进行修正和调整，直到得出较满意的结果为止。下面，我们就通过一些实例说明微分方程建模的具体步骤。一.增长模型在自然界和社会的经济活动中，许多量的变化都遵循着一个基本的规律：任一单位时间内的增量都与该量自身当时的大小成正比。运用这一基本规律，就可以建立起各种各样的增长模型。1.马尔萨斯人口模型严格地讲，讨论人口问题所建立的模型应属于离散型模型。但在人口基数很大的情况下，突然增加或减少的只是单一的个体或少数几个个体，相对于全体数量而言，这种改变量是极其微小的，因此，我们可以近似地假设人口随时间连续变化甚至是可微的。这样，我们就可以采用微分方程的工具来研究这一问题。最早研究人口问题的是英国的经济系家

3、马尔萨斯(Malthus)（17661834） 。他根据百余年的人口资料，经过潜心研究，在 1798 年发表的人口论中首先提出了人口增长模型。他的基本假设是：任一单位时刻人口的增长量与当时的人口总数成正比，且比例系数为常数。于是，设 时刻的人口总数为 ，则单位时间内人口的增长量即为t)(tyty)(根据基本假设，有 （ 为比例系数）ty()(trr令 ，可得微分方程0t2（4.1）yrdt这就是著名的马尔萨斯人口方程。若假设 时的人口总数为 ，则不难求得该方程的00y特解为 (4.2)(0tre即任一时刻的人口总数都遵循指数规律向上增长。人们曾用这个公式对 17001961 达二百六十余年世界

4、的人口资料进行了检验，发现计算结果与人口的实际情况竟然是惊人的吻合！2.放射性元素衰变模型放射性元素的质量随时间的推移而逐渐减少（负增长） ，这种现象称为衰变。由物理学定律知，放射性元素任一时刻的衰变速度与该时刻放射性元素的质量成正比。根据这一原理，我们也可以通过微分方程研究放射性元素衰变的规律。设放射性元素 时刻的质量 ，则其衰变速度就是 ，于是可得t)(tmdtm(4.3)dt其中 是比例常数，可由该元素的半衰期（质量蜕变到一半所需要的时间）确定；0前置负号表明放射性元素的质量 是随时刻 递减的。t如果在初始时刻（ ）放射性元素的质量 ，则可求得该方程的特解为0t 0m(4.4)tet0)

5、(这说明放射性元素的质量也是随时刻按指数规律递减的。为了能将求得的放射性元素衰变规律应用于实际，还必须确定上式中的比例常数 。这时，我们可以假设放射性元素的半衰期为 ，从而有TTem02解之，得 ，于是反映放射性元素衰变规律的（4.4）式又可以表示为T2ln(4.5)tTet2ln0)(并由此可解得(4.6)(ln20tmt它所反映的是放射性元素由初始时刻的质量 衰减到 所需要的时间。0)(t放射性元素的衰变规律常被考古、地质方面的专家用于测定文物和地质的年代，其中最常用的是 （碳-12 的同位素）测定法。 这种方法的原理是：大气层在宇宙射线不断C143的轰击下所产生的中子与氮气作用生成了具有

6、放射性的 ，并进一步氧化为二氧化碳，C14二氧化碳被植物所吸收，而动物又以植物为食物，于是放射性碳就被带到了各种动植物的体内。对于具有放射性的 来说，不论是存在于空气中还是生物体内，它都在不断地蜕C14变。由于活着的生物通过新陈代谢不断地摄取 ，因而使得生物体内的 与空气中的14 14有相同的百分含量；一旦生物死亡之后，随着新陈代谢的停止，尸体内的 就会不C14 C断地蜕变而逐渐减少，因此根据 蜕变减少量的变化情况并利用（4.6）式，就可以判定14生物死亡的时间。下面，我们就来看一个运用 测定法确定年代的具体实例：C141972 年 8 月，湖南长沙出土了马王堆一号墓（注：出土时因墓中女尸历经

7、千年而未腐曾经轰动世界） 。经测定，出土的木炭标本中 的平均原子蜕变速度为 29.78 次/分，而新14砍伐烧成的木炭中 的平均原子蜕变速度为 38.37 次/分；如果 的半衰期取为 5568 年14 C14（注： 的半衰期在各种资料中说法不一，分别有 5568 年、5580 年和 5730 年不等） ，C14那么，怎样才能根据以上数据确定这座墓葬的大致年代呢？在确定衰变时间的公式（4.6）中，由于 和 表示的分别是该墓下葬时和出土时0m)(t木炭标本中 的含量，而测量到的是标本中 的平均原子蜕变速度，所有我们还要对14 C14（4.6）式作进一步的修改：对（4.4）式求导，得 )()(0te

8、tt从而有 0)()(mm上面两式相除，得 )(00tt代入（4.6） ，得(4.7)(ln2tmTt于是，衰变时间由(4.6)式根据 含量的变化情况确定就转化为由（ 4.7）式根据 衰变C14 C14速度的变化情况来确定，这就给实际操作带来了很大的方便。在本例中， 年， 次/分， 虽然表示的是下葬时所烧制的568T78.29)(t)0(木炭中 的衰变速度，但考虑到宇宙射线的强度在数千年内的变化不会很大，因而可以C14假设现代生物体中 的衰变速度与马王堆墓葬时代生物体中 的衰变速度相同，即可14 C14以用新砍伐烧成的木炭中 的平均原子蜕变速度 38.37 次/ 分替代 。代入（4.7）可)0

9、(m求得4（年）203678.9ln256)(0ln2tmTt若以 年或 年计算，则可分别算得 年或 年，即马王堆580T734t2095t一号墓大约是 2000 多年前我国汉代的墓葬。 ( 注：后经进一步考证，确定墓主人为汉代长沙国丞相利仓的夫人，名辛追。)3.固定资产折旧模型企业在进行成本核算的时候，经常要计算固定资产的折旧。一般说来，固定资产在任一单位时刻的折旧额与当时固定资产的价值都是成正比的。试研究固定资产价值 与时间P的函数关系。假定某固定资产五年前购买时的价格是 10000 元，而现在的价值为 6000t元，试估算再过 10 年后的价值。首先我们可以假设 时刻该固定资产的价值为

10、，则在 这段时间t )(tP,t内该固定资产单位时刻的折旧额可表示为 ，由题意可得t()()ktP)0(k令 ，即得0t kdt不难求得该方程的通解为 tkeCtP)(为方便计算，记五年前的时刻为 ，于是有初始条件01)(代入通解，可求得 ，故原方程的特解为10CtketP0)(为确定比例常数 ，可将另一个条件 代入上式，得k65k51解出 ，得k 3ln5k从而有 5ln)(1010)( ttetP这就是价值 与时间 之间的函数关系。于是，再过 10 年（即 ）该固定资产的价Pt 1t5值即为 （元）2160)35(10)(P二.阻滞增长模型与以上所讨论的增长模型不同，实际中存在着大量的另一

11、类增长模型：1.弗尔哈斯特人口模型在人口模型的研究中，马尔萨斯得出了“任一时刻的人口总数都遵循指数规律向上增长”的结论，并得到了一段时期内人口数据的验证，然而，随着人口基数的增大，公式（4.2）所暴露的不足之处也越来越明显了。根据公式（4.2）我们不难计算出，世界人口大约 35 年就要翻一番。事实上，设某时刻的世界人口数为 ，人口增长率为 2%，且经过 年就要翻一番，则有0yT即 Tey0.020.Te解之，即得 （年） 。6.342ln5T于是，我们以 1965 年的世界人口 33.4 亿为基数进行计算，可以得到如下的一系列人口数据：2515 年 200 万亿2625 年 1800 万亿26

12、60 年 3600 万亿若按人均地球表面积（包括水面、船上）计算，2625 年仅为 0.09 平方米/人，也就是说必须人挨着人站着才能挤得下；而 35 年后的 2660 年，人口又翻了一番，那就要人的肩上再站人了。而且随着时间的推移，我们有tlim)(00trey这显然不符合人口发展的实际。这说明，在人口基数不是很大的时候，马尔萨斯人口方程还能比较精确地反映人口增长的实际情况，但当人口数量变得很大时，其精确程度就大大降低了。究其根源，是随着人口的迅速膨胀，资源短缺、环境恶化等问题越来越突出，这些都将限制人口的增长。如果考虑到这些因素，就必须对马尔萨斯人口方程进行修改。1845 年，荷兰的数学、

13、生物学家弗尔哈斯特(Verhulst)提出了一个修改方案，即将方程修改为2byrdt)0(rb其中 、 称为“生命系数” 。由于 ，因此当 不太大时， 这一项相对于rbry2y6可以忽略不计；而当 很大时， 这一项所起的作用就不容忽视了，它降低了人yr y2by口的增长速度。于是，我们就有了下面的人口模型：（4.8）02ybrdtt这是一个可分离变量的一阶微分方程。解之，可得（4.9）)(000(trebyr这就是人口 随时间 的变化规律。下面，我们就对（4.9）作进一步的讨论，并根据它对yt人口的发展情况作一些预测。首先，由于brebyrytrtt )(000(limli即不论人口的基数如何

14、，随着时间的推移，人口总量最终将趋于一个确定的极限值 ；br其次，由 可得2byrdtybry )2(令 ，得 ，易知这正是函数（4.9）的图象（称为“人口增长曲线”或“ 型0ybr2 S曲线” ）拐点的纵坐标，它恰好位于人口总量极限值 brr一半的位置（如图所示） 。 2由于 时 ，故 是递增的，此时 bry20dty称为人口的“加速增长期” ；而当 时 ，故 是递减的，此时称为人口的br20dty“缓慢增长期” 。在利用（4.9）式对人口的发展情况进行预测之前，还必须确定恰当的 值，它可以按b以下方法来计算：由方程（4.8）可得7byrdt其中， 表示人口的理论增长率，而 则表示人口的实际

15、增长率。如果我们以 1965 年dtyyt的人口数 为初值，并把某些生态学家估计的 的自然值 0.029 及人口的实际增91034.r长率 0.02 代入上式，有)1034.(029 9b即可求得 。于是，世界人口的极限值120695.b（亿）br6.765.12若以 1965 年的人口数 为初值，则 2000 年的世界人口将达到934.（亿）.590234)160(9.20eyt2.单种群动物模型所谓单种群动物模型指的是某单一种群动物在一个相对封闭环境中生长衰亡规律的模型。这个模型可以运用微元分析法建立和求解。假设某动物种群 时刻的数量为 ，出生率与死亡率分别为 、 ，且任一t )(tNmn

16、时刻的出生数与死亡数都与当时的动物数量成正比，则在时间间隔 内，该动物,dt的出生量与死亡量就分别为 与 ，由于mdtnt增量 出生量 死亡量故有NdtntN)(即（4.10）mdt)(这也是一个可分离变量的微分方程，显然，它与马尔萨斯人口方程（4.1）在本质上是完全相同的。如果初始时刻该动物种群的数量为 ，即 ，则易求得该方程的特解0N0t为 tnme)(0这说明该动物种群的数量也是遵循指数规律增长的，其中 称为“自然增长率” 。n从以上的结果中可以看出，如果 ，则该动物种群的数量随着时间的推移将会无限制地增加，这显然是不可能的。之所以产生这样的谬误，是因为我们把 、 当成了常mn8量，即视

17、出生率与死亡率一成不变，而事实却并非如此。因为随着动物数量的增加，食品短缺、环境恶化、瘟疫流行等问题都将导致出生率的降低和死亡率的升高，为此，可令 qNpnbam,其中 、 、 、 均为正常数。这说明出生率与死亡率不再是常数，而是 的线性函数abpq（这是最简单的情形） ，前者随 均匀减少，后者随 均匀增加。于是方程（4.10）就被N修正为 Nqbpadt)()()(若记 、 ，则上式可简记为qbkpaA（4.11）)(AkdtN这就是描述动物种群数量 与时间 之间函数关系的微分方程。分离变量并积分，得 tkAeCN将初始条件 代入上式，可求得 ，故方程（4.11）的特解即为0Nt0tkAeA

18、0解出 并化简，可得（4.12）ttkAtkA eBeNeN1)(100其中 ， 。10AB由于 ，这说明该种群中动物的数量 不再随时间的推移AeBNttt limli N无限增加，而是逐渐接近于一个固定的上限 ，因此它还是比较符合实际情况的。事实上，自然界各种动物（包括人类）的生长规律，以及传染病的发病人数、商品销售量的变化趋势等问题，都可以用这个模型来描述。形如（4.11）的微分方程称为逻辑斯蒂(Logistic )方程，该方程的解（4.12）所表示的9曲线称为逻辑斯蒂曲线，即所谓的“ 型曲线” 。S三.传染病模型现代医学的发展已经能够有效地控制和预防许多传染病，如以前曾经肆虐一时的天花上

19、个世纪就已经在世界范围内被消灭，鼠疫、霍乱等恶性传染病也得到了控制，但是一些常见的传染病并没有完全绝迹，一遇到合适的气候就会不断地爆发和流行，至今仍然威胁着人类。而且随着时间的推移，还会出现一些新的传染病，如 AIDS、SARS 等。有些传染病传染得很快，导致很高的致残率和死亡率，危害极大，因而对传染病在人群中传播过程的定量研究并应用于传染病的防治就具有了重要的现实意义。传染病流行过程的研究与其它学科的研究有所不同，科学的数据只能通过以往传染病流行时的报告来获取而不能在人群中实验，而报告的数据往往难于满足研究的需要，因而通过建立数学模型与计算机仿真就成为研究传染病流行过程的重要途径之一。 1.

20、相对封闭环境中的传染病模型假设 个人共同生活在一个相对封闭的环境中，如果其中的一个人感染了某种传染病，N而这种传染病又有一定的潜伏期，因而未发病时人们（包括患者自己）是不知道的。一旦患者发病被隔离时，这种传染病实际上已经在传播了。在这种情况下，传染病在人群中的传播过程是怎样的呢？如果我们以 表示发现首例病人后 时刻被感染的人数，则 就表示此时刻)(tIt )(tIN未被感染的人数。在传染病流行的初期，由于 较小，能接触到感染者的人数少，单位)(I时间内被感染的人数也比较少，因而传播速度较慢；而在传染病流行的后期，由于大多数人已经被感染（ 较大） ，未被感染的人数 已经不多了，所以此时单位时间内

21、)(tI )(tIN被感染的人数也不多，因而传播速度也很慢。排除上述两种极端的情况，当有很多的感染者和很多的未感染者时，传染病的传播速度是很快的。因此，传染病的发病率一方面受感染人数的影响，另一方面也受未感染人数的制约。基于以上的分析，我们就可以建立微分方程如下：（4.13）)(INkdtI其中 是比例常数，可根据发病情况的统计数据来确定。 （注）k不难看出，这个方程与（4.11）完全相同，因而它的通解为 tkNeCI注意到 ，代入上式，可求得 ，故该方程的特解为1)0(I 110tkNeI1即（4.14）tkNeI)(这就是该传染病的传播规律。显然， 单调增，且当 时，有 ，这表明tAI如果

22、任其发展而不采取积极有效的措施，最终所有的人都将被感染。尽管这与实际情况有一定的距离，但在传染病流行的前期这个模型还是可用的，因而传染病学者曾用它来预报过传染病高潮到来的时刻。（注：假设发现首例病人 24 小时后被感染的人数为 3 人，则将 、 代入特24t3I解（4.14） ，即可求出比例常数 的值。 ）k2. 传染病的 SIR 模型在上一个模型中，我们实际上只是对人群进行了简单的划分：感染者与健康人。如果我们把病愈免疫和死亡者也考虑在内，情况就会有很大的不同。首先，我们把人群分成三类：类：称为易感类，该类成员没有染上传染病，但缺乏免疫力，可以被感染；S类：称为传染类，该类成员已经染上了传染

23、病，而且可以传播给 类成员；I S类：称为恢复类或排除类，该类成员具有免疫能力或者已经死亡。R记以上三类人员在 时刻的人数分别为 、 和 ，总人数为 ，显然在任t )(tSI)(tRN一时刻，都有 NtRItS)()(为了建立数学模型，我们作出以下的两点假设：（1）单位时间内，一个病人传染的人数与当时健康者的人数成正比，比例系数为（称为传染系数） ； k（2）单位时间内，病愈免疫（包括死亡）的人数与当时感染者的人数成正比，比例系数为 （称为恢复系数） 。由假设（1） ， 可知 类成员中在 时刻的单位时间内被感染的人数为 ；St )(tIkS由假设（2） ，可知 时刻的单位时间内由 类成员转化为

24、 类的人数为 。由于 时tIR)(tI刻单位时间内各类人员的人数就是该类人员总人数的变化率，故有 ISkdtI且 It11如果初始时刻的感染人数为 ，即 ，则 （此时 ） ，0I0)(I0)(INS0)(R于是我们就可以得到如下的传染病模型：（4.15）ISkdtI其初始条件为 ， （4.16）0)(I0)(N虽然微分方程组（4.15）的解析解难以求得，但是我们可以通过对 、 之间函数关IS系的研究讨论（4.15）解的性态。为此，将（4.15）中的两个方程相除，得（4.17）1SkIkdSI若将 记作 ，则称 为“特征指数” ，且对于同一地区、同一种传染病， 为常kk 数。于是方程（4.17）

25、可表示为（4.18）1SdI其初始条件为 00)(N方程（4.18）是一个可分离变量的微分方程，易求得其满足初始条件的解为(4.19)SSI0ln)(（4.19）表示的只是患者数量 与健康人数量 之间的函数关系，为了建立它们与时刻 之间的函数关系以便作进一步的分析，我们可以利用（4.19）式求 对 的导数：t )(SItdtSdtSIt )1(即 )()tStI由于 是单调减函数，即 ，故由上式可以看出，当 时，)(tS0 S，从而 单调增加；当 时， ，这时 单调减少。这也就是说，0tII)(tI)(tI仅当传染病开始流行时健康人数 超过 的情况下，传染病才会蔓延。因此 是一个0“阀值” （

26、俗称“门槛” ） 。注意到 ，而通常 是很小的，所以可近似地认为0NS0，于是在总人数 不变的情况下，提高门槛 的数值，对制止传染病的蔓延是十NS0分有利的。由于 ，因此一方面要及时采取隔离措施，想方设法降低传染系数 ，k k12另一方面要努力增大恢复系数 ，不断提高该地区的医疗水平、免疫水平和卫生保健水平，以避免传染病的大规模爆发和流行。 四.弱肉强食模型 在生态环境中，所有的动物都要通过食物来生存，不是吃其它的动物，就是吃植物。从这种意义上讲，没有一种动物是完全孤立生活的，而且彼此之间一定是相互影响的。为了搞清楚它们之间相互依存的关系，就必须把它们合在一起进行研究，这就产生了多种群模型。弱

27、肉强食模型就是其中较简单的一种。假设有一个包含两个群体的系统，其中的一个以另一个为食，比如狐狸与兔子、鲸鱼与磷虾、天敌与害虫、正常细胞与肿瘤细胞等等，我们把这种系统笼统地称为“捕食者被食者”系统。设 、 分别表示 时刻被食者与捕食者的数量。如果它们各自单独生活，则被)(txyt食者的增长速率正比与当时的数量，即(4.20)0(xdt而捕食者由于没有捕食对象，其数量减少的速率也正比与当时的数量，即(4.21)(yt现在两者生活在一起，被食者有一部分遭到了捕食者的捕食，增长速率要变缓，即(4.20)式中的 将减少，而且减少的量正比于捕食者的数量，所以(4.20)式应修改为)0()(xydtx类似地

28、，捕食者有了捕食对象，减少的速率也要变缓，即(4.21)式中的 也减少，而且减少的量正比于被食者的数量，所以(4.21)式应修改为)0()(yxdty将修改后的两个方程联立，可得方程组（4.22）yxdty)(称为 Volterra -Lotka 方程，其初始条件是 00)(,)(x虽然这个方程组是非线性方程组，不易直接求出它的解析解，但是我们可以采用传染13病模型中使用过的方法，通过对 、 之间函数关系的研究讨论（4.22）解的性态。为此，xy将（4.22）中的两个方程相除，得（4.23）xydx)(分离变量并积分，可得通解Cylnln利用 及 ，可将上式化为yelnxelxy这就是方程（4

29、.23）的隐式通解；当 、 时， ，于是0000xyeCxye00xye就是方程（4.23 ）的一个隐式特解，它的图象是 平面上的一条闭曲线（示意图如下） 。只要初始值 和 不为零，这条闭曲线就永远不 0xy通过零点。当被食者较多时，由于食物丰富使捕食者不断增加，于是被食者又开始减少；被食者减少使捕食者也随之减少，从而被食者又会增多。如此，两者的数量循环起伏，周而复始，维持着 Ox生态的平衡。下面我们再进一步讨论这个模型的另一个问题：假如由于某种外界的因素（比如天灾）使得捕食者与被食者同时按比例消亡，那么一旦这种因素的作用停止，是被食者恢复得快，还是捕食者恢复得快呢？这个问题可以从方程（4.2

30、2）中找到答案。在方程组（4.22）即 yxdtyx)(yxdtyx中，乘积项 对被食者来说是消亡项，对捕食者来说是生长项。当 、 都减小，x )(t譬如说 、 都减半时，则 、 也随之减半，但 和 却只减为原来的四分之yxy一。这说明数量减半以后，被食者增长速度的减少还不到一半，而捕食者生长的速度却开始变缓，因而被食者受到的损失比捕食者要小。所以灾难过后，被食者的恢复要比捕食者恢复得快。这个推断来自 Volterra 方程，因而也称之为 Volterra 原理。 这个原理在杀虫剂的施用方面得到了证实。1968 年，澳洲的吹棉蚧（一种昆虫）偶然14传入美国后，严重地威胁着美国的柑桔业，为此，人

31、们又引入了它的天敌澳洲瓢虫，才使吹棉蚧的数量得到了控制。后来，园艺家们使用了 DDT，以期能进一步减少吹棉蚧。然而，由于 DDT 在杀灭吹棉蚧的同时也杀死了它的天敌，一段时间过后，人们发现吹棉蚧的数量反而增多了。上述原理在临床医学上也是有意义的。癌症患者体内的正常细胞与癌细胞就构成了一个“捕食者被食者”系统，某些具有较高毒性的药物既能抑制癌细胞的生长，同时也能抑制正常细胞的生长。癌症患者在进行化疗的时候，药物的毒性使作为捕食者的正常细胞也受到了较大的伤害，一定停止用药，作为被食者的癌细胞可能会变本加厉地增长起来。这始终是癌症化疗中需要解决的一个难题。五. 溶液浓度模型溶液浓度问题是工农业生产和

32、治理环境污染中经常要碰到的问题，此类问题通常都可以描述为如下的实验室模型：一个容器有一个入口，一个出口，里面盛满了某种溶液。如果从入口以不变的速率向容器内注入一定浓度的相同溶液（或清水） ，搅拌均匀后以同样的速率从出口排出，假设搅拌是在瞬间完成的，那么容器内溶液浓度的变化规律是怎样的呢？为了讨论的方便，我们来看一个具体的例子。已知容器内盛有 1000 公斤的清水，如果以每分钟 5 公斤的速率注入浓度为 0.2 的盐水且不停地搅拌，并以同样的速率排出搅拌后的盐水，那么经过多少时间能使容器内的含盐量达到 100 公斤？这个问题我们可以采用微元分析的方法建立微分方程求解。为此，可以设 时刻容器t内的含盐量为 ，则此时溶液的浓度为)(ty即 ，10)(ty)(01.ty于是，在时间间隔 内，有,dt进盐量： t52.出盐量： ydty5.从而含盐量的微元即为或 tdt0.105.yt这是一个一阶线性非齐次微分方程，易求得该方程满足初始条件为 的特解为0)()1(205.tey这就是容器内的含盐量 随时刻 变化的规律。将 代入上式，可求得t1y15（分）62.13805.9.2lnt即经过约 2 小时 18 分 37 秒可使容器内的含盐量达到 100 公斤。