微分方程建模学习.doc

微分方程建模一般说来，微分方程建模的方法大致可以分为以下的几个步骤：1.根据实际问题的要求确定要研究的量，包括自变量、未知函数、必要的参数等以及它们各自的变化区间；2.列方程。可以在合理假设的前提下，利用导数表示斜率、速度、变化率的实际意义，根据一些基本定理（几何的、物理的、化学的或生物学的等等）或规律，找出未知函数的导数（或微分）与相关各量之间的等量关系式，建立微分方程并确定定解条件（注：如果没有现成的定理可供利用，也可以用微元分析法与模拟近似法列出微分方程）；3.解微分方程；4.对模型的适用性作出评价，即用已知的数据检验微分方程的解是否与实际相符。若结果与实际存在一定的差距，则还要对方程进行修正和调整，直到得出较满意的结果为止。下面，我们就通过一些实例说明微分方程建模的具体步骤。一.增长模型在自然界和社会的经济活动中，许多量的变化都遵循着一个基本的规律：任一单位时间内的增量都与该量自身当时的大小成正比。运用这一基本规律，就可以建立起各种各样的增长模型。1.马尔萨斯人口模型严格地讲，讨论人口问题所建立的模型应属于离散型模型。但在人口基数很大的情况下，突然增加或减少的只是单一的个体或少数几个个体，相对于全体数量而言，这种改变量是极其微小的，因此，我们可以近似地假设人口随时间连续变化甚至是可微的。这样，我们就可以采用微分方程的工具来研究这一问题。最早研究人口问题的是英国的经济系家马尔萨斯(Malthus)（17661834）。他根据百余年的人口资料，经过潜心研究，在1798年发表的人口论中首先提出了人口增长模型。他的基本假设是：任一单位时刻人口的增长量与当时的人口总数成正比，且比例系数为常数。于是，设时刻的人口总数为，则单位时间内人口的增长量即为根据基本假设，有 （为比例系数）令，可得微分方程 （4.1）这就是著名的马尔萨斯人口方程。若假设时的人口总数为，则不难求得该方程的特解为 (4.2)即任一时刻的人口总数都遵循指数规律向上增长。人们曾用这个公式对17001961达二百六十余年世界的人口资料进行了检验，发现计算结果与人口的实际情况竟然是惊人的吻合！2.放射性元素衰变模型放射性元素的质量随时间的推移而逐渐减少（负增长），这种现象称为衰变。由物理学定律知，放射性元素任一时刻的衰变速度与该时刻放射性元素的质量成正比。根据这一原理，我们也可以通过微分方程研究放射性元素衰变的规律。设放射性元素时刻的质量，则其衰变速度就是，于是可得 (4.3)其中是比例常数，可由该元素的半衰期（质量蜕变到一半所需要的时间）确定；前置负号表明放射性元素的质量是随时刻递减的。如果在初始时刻（）放射性元素的质量，则可求得该方程的特解为 (4.4)这说明放射性元素的质量也是随时刻按指数规律递减的。为了能将求得的放射性元素衰变规律应用于实际，还必须确定上式中的比例常数。这时，我们可以假设放射性元素的半衰期为，从而有解之，得，于是反映放射性元素衰变规律的（4.4）式又可以表示为 (4.5)并由此可解得 (4.6)它所反映的是放射性元素由初始时刻的质量衰减到所需要的时间。放射性元素的衰变规律常被考古、地质方面的专家用于测定文物和地质的年代，其中最常用的是（碳-12的同位素）测定法。 这种方法的原理是：大气层在宇宙射线不断的轰击下所产生的中子与氮气作用生成了具有放射性的，并进一步氧化为二氧化碳，二氧化碳被植物所吸收，而动物又以植物为食物，于是放射性碳就被带到了各种动植物的体内。对于具有放射性的来说，不论是存在于空气中还是生物体内，它都在不断地蜕变。由于活着的生物通过新陈代谢不断地摄取，因而使得生物体内的与空气中的有相同的百分含量；一旦生物死亡之后，随着新陈代谢的停止，尸体内的就会不断地蜕变而逐渐减少，因此根据蜕变减少量的变化情况并利用（4.6）式，就可以判定生物死亡的时间。下面，我们就来看一个运用测定法确定年代的具体实例：1972年8月，湖南长沙出土了马王堆一号墓（注：出土时因墓中女尸历经千年而未腐曾经轰动世界）。经测定，出土的木炭标本中的平均原子蜕变速度为29.78次/分，而新砍伐烧成的木炭中的平均原子蜕变速度为38.37次/分；如果的半衰期取为5568年（注：的半衰期在各种资料中说法不一，分别有5568年、5580年和5730年不等），那么，怎样才能根据以上数据确定这座墓葬的大致年代呢？在确定衰变时间的公式（4.6）中，由于和表示的分别是该墓下葬时和出土时木炭标本中的含量，而测量到的是标本中的平均原子蜕变速度，所有我们还要对（4.6）式作进一步的修改：对（4.4）式求导，得 从而有上面两式相除，得代入（4.6），得 (4.7)于是，衰变时间由(4.6)式根据含量的变化情况确定就转化为由（4.7）式根据衰变速度的变化情况来确定，这就给实际操作带来了很大的方便。在本例中，年，次/分，虽然表示的是下葬时所烧制的木炭中的衰变速度，但考虑到宇宙射线的强度在数千年内的变化不会很大，因而可以假设现代生物体中的衰变速度与马王堆墓葬时代生物体中的衰变速度相同，即可以用新砍伐烧成的木炭中的平均原子蜕变速度38.37次/分替代。代入（4.7）可求得（年）若以年或年计算，则可分别算得年或年，即马王堆一号墓大约是2000多年前我国汉代的墓葬。 ( 注：后经进一步考证，确定墓主人为汉代长沙国丞相利仓的夫人，名辛追。)3.固定资产折旧模型企业在进行成本核算的时候，经常要计算固定资产的折旧。一般说来，固定资产在任一单位时刻的折旧额与当时固定资产的价值都是成正比的。试研究固定资产价值与时间的函数关系。假定某固定资产五年前购买时的价格是10000元，而现在的价值为6000元，试估算再过10年后的价值。首先我们可以假设时刻该固定资产的价值为，则在这段时间内该固定资产单位时刻的折旧额可表示为，由题意可得 令，即得不难求得该方程的通解为为方便计算，记五年前的时刻为，于是有初始条件代入通解，可求得，故原方程的特解为为确定比例常数，可将另一个条件代入上式，得解出，得从而有这就是价值与时间之间的函数关系。于是，再过10年（即）该固定资产的价值即为 （元）二.阻滞增长模型与以上所讨论的增长模型不同，实际中存在着大量的另一类增长模型：1.弗尔哈斯特人口模型 在人口模型的研究中，马尔萨斯得出了“任一时刻的人口总数都遵循指数规律向上增长”的结论，并得到了一段时期内人口数据的验证，然而，随着人口基数的增大，公式（4.2）所暴露的不足之处也越来越明显了。根据公式（4.2）我们不难计算出，世界人口大约35年就要翻一番。事实上，设某时刻的世界人口数为，人口增长率为2%，且经过年就要翻一番，则有 即 解之，即得 （年）。于是，我们以1965年的世界人口33.4亿为基数进行计算，可以得到如下的一系列人口数据：2515年 200万亿2625年 1800万亿2660年 3600万亿 若按人均地球表面积（包括水面、船上）计算，2625年仅为0.09平方米/人，也就是说必须人挨着人站着才能挤得下；而35年后的2660年，人口又翻了一番，那就要人的肩上再站人了。而且随着时间的推移，我们有 这显然不符合人口发展的实际。这说明，在人口基数不是很大的时候，马尔萨斯人口方程还能比较精确地反映人口增长的实际情况，但当人口数量变得很大时，其精确程度就大大降低了。究其根源，是随着人口的迅速膨胀，资源短缺、环境恶化等问题越来越突出，这些都将限制人口的增长。如果考虑到这些因素，就必须对马尔萨斯人口方程进行修改。1845年，荷兰的数学、生物学家弗尔哈斯特(Verhulst)提出了一个修改方案，即将方程修改为 其中、称为“生命系数”。由于，因此当不太大时，这一项相对于可以忽略不计；而当很大时，这一项所起的作用就不容忽视了，它降低了人口的增长速度。于是，我们就有了下面的人口模型： （4.8）这是一个可分离变量的一阶微分方程。解之，可得 （4.9）这就是人口随时间的变化规律。下面，我们就对（4.9）作进一步的讨论，并根据它对人口的发展情况作一些预测。首先，由于 即不论人口的基数如何，随着时间的推移，人口总量最终将趋于一个确定的极限值；其次，由 可得 令，得，易知这正是函数（4.9）的图象（称为“人口增长曲线”或“型曲线”）拐点的纵坐标，它恰好位于人口总量极限值 一半的位置（如图所示）。 由于时，故是递增的，此时 称为人口的“加速增长期”；而当时，故是递减的，此时称为人口的“缓慢增长期”。在利用（4.9）式对人口的发展情况进行预测之前，还必须确定恰当的值，它可以按以下方法来计算：由方程（4.8）可得 其中，表示人口的理论增长率，而则表示人口的实际增长率。如果我们以1965年的人口数为初值，并把某些生态学家估计的的自然值0.029及人口的实际增长率0.02代入上式，有 即可求得 。于是，世界人口的极限值 （亿）若以1965年的人口数为初值，则2000年的世界人口将达到 （亿）2.单种群动物模型所谓单种群动物模型指的是某单一种群动物在一个相对封闭环境中生长衰亡规律的模型。这个模型可以运用微元分析法建立和求解。假设某动物种群时刻的数量为，出生率与死亡率分别为、，且任一时刻的出生数与死亡数都与当时的动物数量成正比，则在时间间隔内，该动物的出生量与死亡量就分别为与，由于增量 出生量 死亡量故有 即 （4.10）这也是一个可分离变量的微分方程，显然，它与马尔萨斯人口方程（4.1）在本质上是完全相同的。如果初始时刻该动物种群的数量为，即，则易求得该方程的特解为这说明该动物种群的数量也是遵循指数规律增长的，其中称为“自然增长率”。从以上的结果中可以看出，如果，则该动物种群的数量随着时间的推移将会无限制地增加，这显然是不可能的。之所以产生这样的谬误，是因为我们把、当成了常量，即视出生率与死亡率一成不变，而事实却并非如此。因为随着动物数量的增加，食品短缺、环境恶化、瘟疫流行等问题都将导致出生率的降低和死亡率的升高，为此，可令其中、均为正常数。这说明出生率与死亡率不再是常数，而是的线性函数（这是最简单的情形），前者随均匀减少，后者随均匀增加。于是方程（4.10）就被修正为若记、，则上式可简记为 （4.11）这就是描述动物种群数量与时间之间函数关系的微分方程。分离变量并积分，得将初始条件代入上式，可求得，故方程（4.11）的特解即为解出并化简，可得 （4.12）其中，。由于，这说明该种群中动物的数量不再随时间的推移无限增加，而是逐渐接近于一个固定的上限，因此它还是比较符合实际情况的。事实上，自然界各种动物（包括人类）的生长规律，以及传染病的发病人数、商品销售量的变化趋势等问题，都可以用这个模型来描述。形如（4.11）的微分方程称为逻辑斯蒂(Logistic)方程，该方程的解（4.12）所表示的曲线称为逻辑斯蒂曲线，即所谓的“型曲线”。三.传染病模型现代医学的发展已经能够有效地控制和预防许多传染病，如以前曾经肆虐一时的天花上个世纪就已经在世界范围内被消灭，鼠疫、霍乱等恶性传染病也得到了控制，但是一些常见的传染病并没有完全绝迹，一遇到合适的气候就会不断地爆发和流行，至今仍然威胁着人类。而且随着时间的推移，还会出现一些新的传染病，如AIDS、SARS等。有些传染病传染得很快，导致很高的致残率和死亡率，危害极大，因而对传染病在人群中传播过程的定量研究并应用于传染病的防治就具有了重要的现实意义。传染病流行过程的研究与其它学科的研究有所不同，科学的数据只能通过以往传染病流行时的报告来获取而不能在人群中实验，而报告的数据往往难于满足研究的需要，因而通过建立数学模型与计算机仿真就成为研究传染病流行过程的重要途径之一。 1.相对封闭环境中的传染病模型假设个人共同生活在一个相对封闭的环境中，如果其中的一个人感染了某种传染病，而这种传染病又有一定的潜伏期，因而未发病时人们（包括患者自己）是不知道的。一旦患者发病被隔离时，这种传染病实际上已经在传播了。在这种情况下，传染病在人群中的传播过程是怎样的呢？如果我们以表示发现首例病人后时刻被感染的人数，则就表示此时刻未被感染的人数。在传染病流行的初期，由于较小，能接触到感染者的人数少，单位时间内被感染的人数也比较少，因而传播速度较慢；而在传染病流行的后期，由于大多数人已经被感染（较大），未被感染的人数已经不多了，所以此时单位时间内被感染的人数也不多，因而传播速度也很慢。排除上述两种极端的情况，当有很多的感染者和很多的未感染者时，传染病的传播速度是很快的。因此，传染病的发病率一方面受感染人数的影响，另一方面也受未感染人数的制约。基于以上的分析，我们就可以建立微分方程如下： （4.13）其中是比例常数，可根据发病情况的统计数据来确定。（注）不难看出，这个方程与（4.11）完全相同，因而它的通解为注意到，代入上式，可求得，故该方程的特解为即 （4.14）这就是该传染病的传播规律。显然，单调增，且当时，有，这表明如果任其发展而不采取积极有效的措施，最终所有的人都将被感染。尽管这与实际情况有一定的距离，但在传染病流行的前期这个模型还是可用的，因而传染病学者曾用它来预报过传染病高潮到来的时刻。（注：假设发现首例病人24小时后被感染的人数为3人，则将、代入特解（4.14），即可求出比例常数的值。）2. 传染病的SIR模型在上一个模型中，我们实际上只是对人群进行了简单的划分：感染者与健康人。如果我们把病愈免疫和死亡者也考虑在内，情况就会有很大的不同。首先，我们把人群分成三类：类：称为易感类，该类成员没有染上传染病，但缺乏免疫力，可以被感染；类：称为传染类，该类成员已经染上了传染病，而且可以传播给类成员；类：称为恢复类或排除类，该类成员具有免疫能力或者已经死亡。记以上三类人员在时刻的人数分别为、和，总人数为，显然在任一时刻，都有为了建立数学模型，我们作出以下的两点假设：（1）单位时间内，一个病人传染的人数与当时健康者的人数成正比，比例系数为（称为传染系数）； （2）单位时间内，病愈免疫（包括死亡）的人数与当时感染者的人数成正比，比例系数为（称为恢复系数）。由假设（1）， 可知类成员中在时刻的单位时间内被感染的人数为；由假设（2），可知时刻的单位时间内由类成员转化为类的人数为。由于时刻单位时间内各类人员的人数就是该类人员总人数的变化率，故有且如果初始时刻的感染人数为，即，则（此时），于是我们就可以得到如下的传染病模型： （4.15）其初始条件为 ， （4.16）虽然微分方程组（4.15）的解析解难以求得，但是我们可以通过对、之间函数关系的研究讨论（4.15）解的性态。为此，将（4.15）中的两个方程相除，得 （4.17）若将记作，则称为“特征指数”，且对于同一地区、同一种传染病，为常数。于是方程（4.17）可表示为 （4.18）其初始条件为 方程（4.18）是一个可分离变量的微分方程，易求得其满足初始条件的解为 (4.19)（4.19）表示的只是患者数量与健康人数量之间的函数关系，为了建立它们与时刻之间的函数关系以便作进一步的分析，我们可以利用（4.19）式求对的导数： 即 由于是单调减函数，即，故由上式可以看出，当时，从而单调增加；当时，这时单调减少。这也就是说，仅当传染病开始流行时健康人数超过的情况下，传染病才会蔓延。因此是一个“阀值”（俗称“门槛”）。注意到，而通常是很小的，所以可近似地认为，于是在总人数不变的情况下，提高门槛的数值，对制止传染病的蔓延是十分有利的。由于，因此一方面要及时采取隔离措施，想方设法降低传染系数，另一方面要努力增大恢复系数，不断提高该地区的医疗水平、免疫水平和卫生保健水平，以避免传染病的大规模爆发和流行。 四.弱肉强食模型 在生态环境中，所有的动物都要通过食物来生存，不是吃其它的动物，就是吃植物。从这种意义上讲，没有一种动物是完全孤立生活的，而且彼此之间一定是相互影响的。为了搞清楚它们之间相互依存的关系，就必须把它们合在一起进行研究，这就产生了多种群模型。弱肉强食模型就是其中较简单的一种。假设有一个包含两个群体的系统，其中的一个以另一个为食，比如狐狸与兔子、鲸鱼与磷虾、天敌与害虫、正常细胞与肿瘤细胞等等，我们把这种系统笼统地称为“捕食者被食者”系统。设、分别表示时刻被食者与捕食者的数量。如果它们各自单独生活，则被食者的增长速率正比与当时的数量，即 (4.20)而捕食者由于没有捕食对象，其数量减少的速率也正比与当时的数量，即 (4.21)现在两者生活在一起，被食者有一部分遭到了捕食者的捕食，增长速率要变缓，即(4.20)式中的将减少，而且减少的量正比于捕食者的数量，所以(4.20)式应修改为 类似地，捕食者有了捕食对象，减少的速率也要变缓，即(4.21)式中的也减少，而且减少的量正比于被食者的数量，所以(4.21)式应修改为 将修改后的两个方程联立，可得方程组 （4.22）称为Volterra -Lotka方程，其初始条件是 虽然这个方程组是非线性方程组，不易直接求出它的解析解，但是我们可以采用传染病模型中使用过的方法，通过对、之间函数关系的研究讨论（4.22）解的性态。为此，将（4.22）中的两个方程相除，得 （4.23）分离变量并积分，可得通解 利用 及 ，可将上式化为这就是方程（4.23）的隐式通解；当、时，于是就是方程（4.23）的一个隐式特解，它的图象是平面上的一条闭曲线（示意图如下）。只要初始值和不为零，这条闭曲线就永远不 通过零点。当被食者较多时，由于食物丰富使捕食者不断增加，于是被食者又开始减少；被食者减少使捕食者也随之减少，从而被食者又会增多。如此，两者的数量循环起伏，周而复始，维持着 生态的平衡。下面我们再进一步讨论这个模型的另一个问题：假如由于某种外界的因素（比如天灾）使得捕食者与被食者同时按比例消亡，那么一旦这种因素的作用停止，是被食者恢复得快，还是捕食者恢复得快呢？这个问题可以从方程（4.22）中找到答案。在方程组（4.22） 即 中，乘积项对被食者来说是消亡项，对捕食者来说是生长项。当、都减小，譬如说、都减半时，则、也随之减半，但和却只减为原来的四分之一。这说明数量减半以后，被食者增长速度的减少还不到一半，而捕食者生长的速度却开始变缓，因而被食者受到的损失比捕食者要小。所以灾难过后，被食者的恢复要比捕