1、微分方程模 型,数学建模,1.了解问题的实际背景,明确建模目的,收集掌握必要的数据资料。2.在明确建模目的,掌握必要资料的基础上,通过对资料的分析计 算, 找出起主要作用的因素,经必要的精炼、简化,提出若干符合客观实际的假设。3.在所作假设的基础上,利用适当的数学工具去刻划各变量之间的关系,建立相应的数学结构 即建立数学模型。 4.模型求解。 5.模型的分析与检验。,3.0 数学建模的一般步骤,3.1 微分方程的几个简单实例,在许多实际问题中,涉及变化率时,可导出包含未知函数的导数或微分的关系式,从而用建立微分方程模型的方法来研究该问题,,本节将通过一些最简单的实例来说明微分方程建模的一般方法
2、。在连续变量问题的研究中,微分方程是十分常用的数学工具之一。,例1 (理想单摆运动)建立理想单摆运动满足的微分方程,并得出理想单摆运动的周期公式。,从图3-1中不难看出,小球所受的合力为mgsin,根据牛顿第二定律可得:,这是理想单摆应满足的运动方程,(3.1)是一个两阶非线性方程,不易求解。当很小时,sin,此时,可考察(3.1)的近似线性方程:,由此即可得出,(3.1)的近似方程,为了保持自然资料的合理开发与利用,人类必须保持并控制生态平衡,甚至必须控制人类自身的增长。本节将建立几个简单的单种群增长模型,以简略分析一下这方面的问题。,种群的数量本应取离散值,但由于种群数量一般较大,为建立微
3、分方程模型,可将种群数量看作连续变量,由此引起的误差将是十分微小的。,3.2 Malthus模型与Logistic模型,模型1 马尔萨斯(Malthus)模型,马尔萨斯在分析人口出生与死亡情况的资料后发现,人口净增长率r基本上是一常数,(r=b-d,b为出生率,d为死亡率), 既:,马尔萨斯模型的一个显著特点:种群数量翻一番所需的时间是固定的。,Malthus模型实际上只有在群体总数不太大时才合理,到总数增大时,生物群体的各成员之间由于有限的生存空间,有限的自然资源及食物等原因,就可能发生生存竞争等现象。,所以Malthus模型假设的人口净增长率不可能始终保持常数,它应当与人口数量有关。,模型
4、2 Logistic模型,人口净增长率应当与人口数量有关,即: r=r(N),r(N)是未知函数,但根据实际背景,它无法用拟合方法来求 。,为了得出一个有实际意义的模型,我们不妨采用一下工程师原则。工程师们在建立实际问题的数学模型时,总是采用尽可能简单的方法。,r(N)最简单的形式是常数,此时得到的就是马尔萨斯模型。对马尔萨斯模型的最简单的改进就是引进一次项(竞争项),(3.9)式还有另一解释,由于空间和资源都是有限的,不可能供养无限增长的种群个体,当种群数量过多时,由于人均资源占有率的下降及环境恶化、疾病增多等原因,出生率将降低而死亡率却会提高。设环境能供养的种群数量的上界为K(近似地将K看
5、成常数),N表示当前的种群数量,K-N恰为环境还能供养的种群数量,(3.9)指出,种群增长率与两者的乘积成正比,正好符合统计规律,得到了实验结果的支持,这就是(3.9)也被称为统计筹算律的原因。,图3-5,对(3.9)分离变量:,两边积分并整理得:,令N(0)=N0,求得:,N(t)的图形请看图3.5,大量实验资料表明用Logistic模型来描述种群的增长,效果还是相当不错的。例如,高斯把5只草履虫放进一个盛有0.5cm3营养液的小试管,他发现,开始时草履虫以每天230.9%的速率增长,此后增长速度不断减慢,到第五天达到最大量375个,实验数据与r=2.309,a=0.006157,N(0)=
6、5的Logistic曲线:几乎完全吻合,见图3.6。,图3-6,Malthus模型和Logistic模型的总结,Malthus模型和Logistic模型均为对微分方程(3.7)所作的模拟近似方程。前一模型假设了种群增长率r为一常数,(r被称为该种群的内禀增长率)。后一模型则假设环境只能供养一定数量的种群,从而引入了一个竞争项。,用模拟近似法建立微分方程来研究实际问题时必须对求得的解进行检验,看其是否与实际情况相符或基本相符。相符性越好则模拟得越好,否则就得找出不相符的主要原因,对模型进行修改。,Malthus模型与Logistic模型虽然都是为了研究种群数量的增长情况而建立的,但它们也可用来研
7、究其他实际问题,只要这些实际问题的数学模型有相同的微分方程即可。,例2 新产品的推广,经济学家和社会学家一直很关心新产品的推销速度问题。怎样建立一个数学模型来描述它,并由此析出一些有用的结果以指导生产呢?以下是第二次世界大战后日本家电业界建立的电饭包销售模型。,设需求量有一个上界,并记此上界为K,记t时刻已销售出的电饭包数量为x(t),则尚未使用的人数大致为Kx(t),于是由统计筹算律:,记比例系数为k,则x(t)满足:,此方程即Logistic模型,解为:,对x(t)求一阶、两阶导数:,x(t)0,即x(t)单调增加。,令x(t0)=0,有,当tt0时,x(t)单调减小。,在销出量小于最大需
8、求量的一半时,销售速度是不断增大的,销出量达到最大需求量的一半时,该产品最为畅销,接着销售速度将开始下降。,所以初期应采取小批量生产并加以广告宣传;从有20%用户到有80%用户这段时期,应该大批量生产;后期则应适时转产,这样做可以取得较高的经济效果。,3.3 药物在体内的分布,何为房室系统?,在用微分方程研究实际问题时,人们常常采用一种叫“房室系统”的观点来考察问题。根据研究对象的特征或研究的不同精度要求,我们把研究对象看成一个整体(单房室系统)或将其剖分成若干个相互存在着某种联系的部分(多房室系统)。,房室具有以下特征:它由考察对象均匀分布而成,房室中考察对象的数量或浓度(密度)的变化率与外
9、部环境有关,这种关系被称为“交换”且交换满足着总量守衡。在本节中,我们将用房室系统的方法来研究药物在体内的分布。在下一节中,我们将用多房室系统的方法来研究另一问题。,药物的分解与排泄(输出)速率通常被认为是与药物当前的浓度成正比的,即:,药物分布的单房室模型,单房室模型是最简单的模型,它假设:体内药物在任一时刻都是均匀分布的,设t时刻体内药物的总量为x(t);系统处于一种动态平衡中,即成立着关系式:,药物的输入规律与给药的方式有关。下面,我们来研究一下在几种常见的给药方式下体内药体的变化规律。,情况1 快速静脉注射,与放射性物质类似,医学上将血浆药物浓度衰减一半所需的时间称为药物的血浆半衰期:
10、,情况2 恒速静脉点滴,易见:,对于多次点滴,设点滴时间为T1,两次点滴之间的间隔时间设为T2,则在第一次点滴结束时病人体内的药物浓度可由上式得出。其后T2时间内为情况1。故:,类似可讨论以后各次点滴时的情况,区别只在初值上的不同。第二次点滴起,患者 体内的初始药物浓度不为零。,情况3 口服药或肌注,口服药或肌肉注射时,药物的吸收方式与点滴时不同,药物虽然瞬间进入了体内,但它一般都集中与身体的某一部位,靠其表面与肌体接触而逐步被吸收。设药物被吸收的速率与存量药物的数量成正比,记比例系数为K1,即若记t时刻残留药物量为y(t),则y满足:,因而:,所以:,解得:,从而药物浓度:,图3-9给出了上
11、述三种情况下体内血药浓度的变化曲线。容易看出,快速静脉注射能使血药浓度立即达到峰值,常用于急救等紧急情况;口服、肌注与点滴也有一定的差异,主要表现在血药浓度的峰值出现在不同的时刻,血药的有效浓度保持时间也不尽相同。,图3-9,我们已求得三种常见给药方式下的血药浓度C(t),当然也容易求得血药浓度的峰值及出现峰值的时间,因而,也不难根据不同疾病的治疗要求找出最佳治疗方案。,上述研究是将机体看成一个均匀分布的同质单元,故被称单房室模型,但机体事实上并不是这样。药物进入血液,通过血液循环药物被带到身体的各个部位,又通过交换进入各个器官。因此,要建立更接近实际情况的数学模型就必须正视机体部位之间的差异
12、及相互之间的关联关系,这就需要多房室系统模型。,图3-10表示的是一种常见的两房室模型,其间的k12表示由室I渗透到室II的变化率前的系数,而k21则表示由室II返回室I的变化率前的系数,它们刻划了两室间的内在联系,其值应当用实验测定,使之尽可能地接近实际情况。,当差异较大的部分较多时,可以类似建立多房室系统,即N房室系统,hy3.5 传染病模型,传染病是人类的大敌,通过疾病传播过程中若干重要因素之间的联系建立微分方程加以讨论,研究传染病流行的规律并找出控制疾病流行的方法显然是一件十分有意义的工作。在本节中,我们将主要用多房室系统的观点来看待传染病的流行,并建立起相应的多房室模型。,医生们发现
13、,在一个民族或地区,当某种传染病流传时,波及到的总人数大体上保持为一个常数。即既非所有人都会得病也非毫无规律,两次流行(同种疾病)的波及人数不会相差太大。如何解释这一现象呢?试用建模方法来加以证明。,问题的提出:,设某地区共有n+1人,最初时刻共有i人得病,t时刻已感染(infective)的病人数为i(t),假定每一已感染者在单位时间内将疾病传播给k个人(k称为该疾病的传染强度),且设此疾病既不导致死亡也不会康复,模型1,此模型即Malthus模型,它大体上反映了传染病流行初期的病人增长情况,在医学上有一定的参考价值,但随着时间的推移,将越来越偏离实际情况。,已感染者与尚未感染者之间存在着明
14、显的区别,有必要将人群划分成已感染者与尚未感染的易感染,对每一类中的个体则不加任何区分,来建立两房室系统。,模型2,记t时刻的病人数与易感染人数(susceptible)分别为i(t)与s(t),初始时刻的病人数为 i。根据病人不死也不会康复的假设及(竞争项)统计筹算律,,其中:,统计结果显示,(3.17)预报结果比(3.15)更接近实际情况。医学上称曲线 为传染病曲线,并称 最大值时刻t1为此传染病的流行高峰。,模型2仍有不足之处,它无法解释医生们发现的现象,且当时间趋与无穷时,模型预测最终所有人都得病,与实际情况不符。,为了使模型更精确,有必要再将人群细分,建立多房室系统,(3.18),求
15、解过程如下:,对(3)式求导,由(1)、(2)得:,解得:,将人群划分为三类(见右图):易感染者、已感染者和已恢复者(recovered)。分别记t时刻的三类人数为s(t)、i(t)和r(t),则可建立下面的三房室模型:,模型3,由(1)式可得:,从而解得:,为揭示产生上述现象的原因(3.18)中的第(1)式改写成:,其中 通常是一个与疾病种类有关的 较大的常数。,下面对 进行讨论,请参见右图,如果 ,则开始时 ,i(t)单增。但在i(t)增加的同时, 伴随地有s(t)单减。当s(t)减少到小于等于 时, i(t)开始减小,直至此疾病在该地区消失。,鉴于在本模型中的作用, 被医生们称为此疾病在
16、该地区的阀值。 的引入解释了为什么此疾病没有波及到该地区的所有人。,图3-14,综上所述,模型3指出了传染病的以下特征:,(1)当人群中有人得了某种传染病时,此疾病并不一定流传,仅当易受感染的人数与超过阀值时,疾病才会流传起来。,(2)疾病并非因缺少易感染者而停止传播,相反,是因为缺少传播者才停止传播的,否则将导致所有人得病。,(3)种群不可能因为某种传染病而绝灭。,模型检验:,医疗机构一般依据r(t)来统计疾病的波及人数 ,从广义上理解,r(t)为t时刻已就医而被隔离的人数,是康复还是死亡对模型并无影响。,通常情况下,传染病波及的人数占总人数的百分比不会太大,故 一般是小量。利用泰勒公式展开
17、取前三项,有:,代入(3.20)得近似方程:,积分得:,其中:,这里双曲正切函数 :,而:,图3-14给出了(3.21)式曲线的图形,可用医疗单位每天实际登录数进行比较拟合得最优曲线。,图3-14(a),3.7 稳定性问题,在研究许多实际问题时,人们最为关心的也许并非系统与时间有关的变化状态,而是系统最终的发展趋势。例如,在研究某频危种群时,虽然我们也想了解它当前或今后的数量,但我们更为关心的却是它最终是否会绝灭,用什么办法可以拯救这一种群,使之免于绝种等等问题。要解决这类问题,需要用到微分方程或微分方程组的稳定性理论。在下两节,我们将研究几个与稳定性有关的问题。,一般的微分方程或微分方程组可
18、以写成:,若方程或方程组f(x)=0有解Xo,X=Xo显然满足(3.28)。称点Xo为微分方程或微分方程组(3.28)的平衡点或奇点。,例7 本章第2节中的Logistic模型,共有两个平衡点:N=0和N=K,分别对应微分方程的两两个特殊解。前者为No=0时的解而后者为No=K时的解。,当NoK时,则位于N=K的上方。从图3-17中不难看出,若No0,积分曲线在N轴上的投影曲线(称为轨线)将趋于K。这说明,平衡点N=0和N=K有着极大的区别。,图3-17,定义1 自治系统 的相空间是指以(x1,xn)为坐标 的空间Rn。,特别,当n=2时,称相空间为相平面。,空间Rn的点集(x1,xn)|xi
19、=xi(t)满足(3.28),i=1,n称为系统的轨线,所有轨线在相空间的分布图称为相图。,定义2 设x0是(3.28)的平衡点,称:,(1)x0是稳定的,如果对于任意的0,存在一个0,只要|x(0)- x0|,就有|x(t)- x0|对所有的t都成立。,(2)x0是渐近稳定的,如果它是稳定的且 。,微分方程平衡点的稳定性除了几何方法,还可以通过解析方法来讨论,所用工具为以下一些定理。,(3)x0是不稳定的,如果(1)不成立。,根据这一定义,Logistic方程的平衡点N=K是稳定的且为渐近稳定的,而平衡点N=0则是不稳定的。,解析方法,证 由泰勒公式,当x与xo充分接近时,有:,由于xo是平
20、衡点,故f(xo)=0。若 ,则当x0,从而x单增;当xxo时,又有f(x)0,从而x单减。无论在哪种情况下都有xxo,故xo是渐进稳定的。,的情况可类似加以讨论。,高阶微分方程与高阶微分方程组平衡点的稳定性讨论较为复杂,大家有兴趣可参阅微分方程定性理论。为了下两节的需要,我们简单介绍一下两阶微分方程组平衡点的稳定性判别方法。,考察两阶微分方程组:,(3.29),令 ,作一坐标平移,不妨仍用x记x,则平衡点xo的稳定性讨论转化为原点的稳定性讨论了。将f(x1,x2)、g(x1,x2)在原点展开,(3.29)又可写成:,其中:,令p=a+d, q=ad-bc=|A|,则 ,记 。,讨论特征值与零
21、点稳定的关系, 如果只有一个特征向量当p0时,零点不 稳定当p0时,零点稳定,综上所述:仅当p0时, (3.30)零点才是渐近稳定的;当p=0且q0时(3.30)有周期解,零点是稳定的中心(非渐近稳定);在其他情况下,零点均为不稳定的。,非线性方程组(3.29)平衡点稳定性讨论可以证明有下面定理成立:,定理2 若(3.30)的零点是渐近稳定的,则(3.29)的平衡点 也是渐近稳定的;若(3.30)的零点是不稳定的,则(3.29) 的平衡点也是不稳定的。,3.8 捕食系统的Volterra方程,问题背景:,意大利生物学家DAncona曾致力于鱼类种群相互制约关系的研究,在研究过程中他无意中发现了
22、一些第一次世界大战期间地中海沿岸港口捕获的几种鱼类占捕获总量百分比的资料,从这些资料中他发现各种软骨掠肉鱼,如鲨鱼、鳐鱼等我们称之为捕食者(或食肉鱼)的一些不是很理想的鱼类占总渔获量的百分比。在 19141923年期间,意大利阜姆港收购的鱼中食肉鱼所占的比例有明显的增加:,他知道,捕获的各种鱼的比例近似地反映了地中海里各种鱼类的比例。战争期间捕鱼量大幅下降,但捕获量的下降为什么会导致鲨鱼、鳐鱼等食肉鱼比例的上升,即对捕食者有利而不是对食饵有利呢?他百思不得其解,无法解释这一现象,就去求教当时著名的意大利数学家V.Volterra,希望他能建立一个数学模型研究这一问题。,Volterra将鱼划分
23、为两类。一类为食用鱼(食饵),数量记为x1(t),另一类为食肉鱼(捕食者),数量记为x2(t),并建立双房室系统模型。,1、模型建立,大海中有食用鱼生存的足够资源,可假设食用鱼独立生存将按增长率为r1的指数律增长(Malthus模型),既设:,由于捕食者的存在,食用鱼数量因而减少,设减少的速率与两者数量的乘积成正比(竞争项的统计筹算律),即:,对于食饵(Prey)系统 :,对于捕食者(Predator)系统 :,捕食者设其离开食饵独立存在时的死亡率为r2,即:,但食饵提供了食物,使生命得以延续。这一结果也要通过竞争来实现,再次利用统计筹算律,得到:,方程组(3.31)反映了在没有人工捕获的自然
24、环境中食饵与捕食者之间的相互制约关系。下面我们来分析该方程组。,2、模型分析,方程组(3.31)是非线性的,不易直接求解。容易看出,该方程组共有两个平衡点,即:,方程组还有两组平凡解:,和,和,当x1(0)、x2(0)均不为零时, ,应有x1(t)0且x2(t)0,相应的相轨线应保持在第一象限中。,求(3.31)的相轨线,将两方程相除消去时间t,得:,令,用微积分知识容易证明:,有:,与 的图形见图3-20,易知仅当 时(3.32)才有解,当 时,轨线退化为平衡点。,当 时,轨线为一封闭曲线(图3-21),即周期解。,证明具有周期解。,只需证明:存在两点 及 , 时,方 程无解。,由 的性质,
25、 , 而 ,使得:,。同样根据的性质知,当 x1 时,。此时:,由 的性质, ,使 成立。,当x1= 或 时, ,,仅当 时才能成立。,而当x1 时,由于 ,,故 无解。,得证。,确定闭曲线的走向,在每一子区域, 与 不变号,据此确定轨线的走向(图3-22),将Volterra方程中的第二个改写成:,将其在一个周期长度为T的区间上积分,得,等式左端为零,故可得:,同理:,解释DAncona发现的现象,引入捕捞能力系数,(01),表示单位时间 内捕捞起来的鱼占总量的百分比。故Volterra方程应为:,平衡点P的位置移动到了:,由于捕捞能力系数的引入,食用鱼的平均量有了增加,而食肉鱼的平均量却有
26、所下降,越大,平衡点的移动也越大。,食用鱼的数量反而因捕捞它而增加, 真的是这样?!,P-P模型导出的结果虽非绝对直理,但在一定程度上是附合客观实际的,有着广泛的应用前景。例如,当农作物发生病虫害时,不要随随便便地使用杀虫剂,因为杀虫剂在杀死害虫的同时也可能杀死这些害虫的天敌,(害虫与其天敌构成一个双种群捕食系统),这样一来,使用杀虫剂的结果会适得其反,害虫更加猖獗了。,(3)捕鱼对食用鱼有利而对食肉鱼不利,多捕鱼(当然要在一定限度内,如r1)能使食用鱼的平均数量增加而使食肉鱼的平均数量减少。,根据P-P模型,我们可以导出以下结论:,(1)食用鱼的平均量取决于参数r1与1,(2)食用鱼繁殖率r
27、1的减小将导致食肉鱼平均量的减小,食肉鱼捕食能力1的增大也会使自己的平均量减小;反之,食肉鱼死亡率r2的降低或食饵对食肉鱼供养效率2的提高都将导致食用鱼平均量的减少。,3.9 较一般的双种群生态系统,Volterra的模型揭示了双种群之间内在的互相制约关系,成功解释了DAncona发现的现象。然而,对捕食系统中存在周期性现象的结论,大多数生物学家并不完全赞同,因为更多的捕食系统并没有这种特征。,一个捕食系统的数学模型未必适用于另一捕食系统,捕食系统除具有共性外,往往还具有本系统特有的个性,反映在数学模型上也应当有所区别。现考察较为一般的双种群系统。,一般的双种群系统,Ki随种群不同而不同,同时
28、也随系统状态的不同而不同,即Ki应为x1、x2的函数。Ki究竟是一个怎样的函数,我们没有更多的信息。不妨再次采用一下工程师们的原则,采用线性化方法。这样,得到下面的微分方程组:,(3.33)不仅可以用来描述捕食系统。也可以用来描述相互间存在其他关系的种群系统。,(3.33),(3.33)式的一些说明,式中a1、b2为本种群的亲疏系数,a2、b1为两种群间的交叉亲疏系数。a2b10时,两种群间存在着相互影响,此时又可分为以下几类情况:,(i)a20,b10,共栖系统。,(ii)a20( 或a20,b10 ),捕食系统。,(iii)a20,b10,竞争系统。,(i)(iii)构成了生态学中三个最基
29、本的类型,种群间较为复杂的关系可以由这三种基本关系复合而成。,(3.33)是否具有周期解,不同的系统具有不同的系数,在未得到这些系数之前先来作一个一般化的讨论。,首先,系统的平衡点为方程组:,(3.34),的解。,如果系统具有非平凡平衡点 则它应当对应于方程组,均为平凡平衡点。,的根,解得:,证明:,记,假设结论不真,则在x1x2平面第一象限存在(3.33)的一个圈,它围成的平面区域记为R。,于是由K(x1,x2)0且连续以及AB0可知,函数 在第一象限中不变号且不为零,故二重积分:,(3.35),但另一方面,由格林公式,注意到 , ,又有:,(3.36),其中T为周期。,(3.35)与(3.
30、36)矛盾,说明圈不可能存在。,对于Voltera方程,由a1=b2=0,得B=0;所以无圈定理不适用于Volterra方程。,对于一般的生态系统,如果通过求解的微分方程来讨论常常会遇到困难。,怎样来讨论一般的生态系统,如果困难的话可以研究种群的变化率,搞清轨线的走向来了解各种群数量的最终趋势。,简化模型,设竞争系统的方程为:,其中不为0,否则为Logistic模型 。,方便讨论取=1,但所用方法可适用一般情况。,作直线l1: x1+x2=K1及l2: x1+x2=K2, K1 K2,见图3-26。,有以下几个引理:,引理1 若初始点位于区域I中,则解 (x1(t)、x2(t))从某一时刻起
31、必开此区域而进入区域II,引理2 若初始点(x1(0)、x2(0))位于 区域II中,则(x1(t),x2(t))始 终位于II中,且:,引理3 若初始点位于区域III中,且对于 任意t ,(x1(t),x2(t))仍位于 III中,则当t+时,(x1(t),x2(t))必以(K1,0)为极限点。,由引理1和引理2,初始点位于像限I和II的解必趋于平衡点(K1,0)。由引理3,初始点位于III且(x1(t),x2(t))始终位于III中的解最终必趋于平衡点(K1,0),而在某时刻进入区域II的解由引理最终也必趋于(K1,0)。易见只有上述三种可能,而在三种可能情况下(x1(t),x2(t))均
32、以(K1,0)为极限,定理得证。,定理4的证明:,在研究实际课题时,数值解方法也许会用得更多。当解析解无法求得时,计算机作为强大的辅助工具发挥了它应起的作用。我校学生在研究1999年美国大学生数学建模竞赛题A(小行星撞击地球)时就遇到了一个棘手的问题:如何描述南极地区的生态系统,如何定量化地研究小行星撞击地球对南级生态环境的影响?在上网查阅了南极附近的海洋生态状况后,他们将南极附近的生物划分成三个部分:海藻、鳞虾和其他海洋生物。鳞虾吃海藻,其他海洋动物吃鳞虾,运用基本建模技巧建立了一个三房室系统模型。小行星的撞击会影响大气层的能见度,从而影响到海藻的生长(光合作用),进而影响到生物链中的其他生物。他们无法得到模型中的参数值(事实上,小行星撞击南极的事件并未发生过),就取了一系列不同的参数值,对不同参数值下模型的数值解进行了分析对比,研究了解对各参数变化的灵敏度,取得了十分有意义的结果并获得了当年国际竞赛的一等奖。,