微分方程建模实例

常微分方程在数学建模中的应用这里介绍几个典型的用微分方程建立数学模型的例子. 一、人口预测模型由于资源的有限性,当今世界各国都注意有计划地控制人口的增长,为了得到人口预测模型,必须首先搞清影响人口增长的因素,而影响人口增长的因素很多,如人口的自然出生率、人口的自然死亡率、人口的迁移、自然灾害、战争等诸多因素,如果一开始就把所有因素都考虑进去，则无从下手.因此,先把问题简化,建立比较粗糙的模型,再逐步修改,得到较完善的模型.例1（马尔萨斯（Malthus）模型）英国人口统计学家马尔萨斯（1766—1834）在担任牧师期间,查看了教堂100 多年人口出生统计资料,发现人口出生率是一个常数,于1789 年在《人口原理》一书中提出了闻名于世的马尔萨斯人口模型,他的基本假设是：在人口自然增长过程中,净相对增长（出生率与死亡率之差）是常数,即单位时间内人口的增长量与人口成正比,比例系数设为r ,在此假设下,推导并求解人口随时间变化的数学模型. 解设时刻t 的人口为N (t ) ,把N (t ) 当作连续、可微函数处理（因人口总数很大,可近似地这样处理,此乃离散变量连续化处理）,据马尔萨斯的假设,在t 到t + ∆t时间段内,人口的增长量为N (t + ∆t ) − N (t ) = rN (t )∆t ,并设t = t 0 时刻的人口为N 0 ,于是dN = rN， dt N (t 0 ) = N 0．这就是马尔萨斯人口模型,用分离变量法易求出其解为N (t ) = N 0 e r (t −t0 ) ,此式表明人口以指数规律随时间无限增长. 模型检验：据估计1961 年地球上的人口总数为3.06 × 10 ,而在以后7 年中,人口总数9以每年2%的速度增长,这样t 0 = 1961 , N 0 = 3.06 × 109, r = 0.02 ,于是N (t ) = 3.06 × 10 9 e 0.02(t −1961) .这个公式非常准确地反映了在1700—1961 年间世界人口总数.因为,这期间地球上的人口大约每35 年翻一番,而上式断定34.6 年增加一倍（请读者证明这一点）．但是,后来人们以美国人口为例,用马尔萨斯模型计算结果与人口资料比较,却发现有很大的差异,尤其是在用此模型预测较遥远的未来地球人口总数时,发现更令人不可思议的问题,如按此模型计算,到2670 年,地球上将有36 000 亿人口.如果地球表面全是陆地（事实上, 地球表面还有80%被水覆盖）,我们也只得互相踩着肩膀站成两层了,这是非常荒谬的,因此, 这一模型应该修改. 例2（逻辑Logistic 模型）马尔萨斯模型为什么不能预测未来的人口呢？这主要是地7球上的各种资源只能供一定数量的人生活,随着人口的增加,自然资源环境条件等因素对人口增长的限制作用越来越显著,如果当人口较少时,人口的自然增长率可以看作常数的话,那么当人口增加到一定数量以后,这个增长率就要随人口的增加而减小.因此,应对马尔萨斯模型中关于净增长率为常数的假设进行修改. 1838 年,荷兰生物数学家韦尔侯斯特(Verhulst)引入常数N m ,用来表示自然环境条件所能容许的最大人口数（一般说来,一个国家工业化程度越高,它的生活空间就越大,食物就越多,从而N m 就越大）,并假设将增长率等于r 1 −N (t ) ,即净增长率随着N (t ) 的增加而Nm减小,当N (t ) → N m 时,净增长率趋于零,按此假定建立人口预测模型. 解由韦尔侯斯特假定,马尔萨斯模型应改为dN N = r 1 − N N， dt 0 N (t ) = N ，0 0上式就是逻辑模型,该方程可分离变量,其解为,N (t ) =NmN 1 + m − 1 e −r (t −t0 ) N 0.下面,我们对模型作一简要分析. （1）当t → ∞ , N (t ) → N m ,即无论人口的初值如何,人口总数趋向于极限值N m ；（2）当0 < N < N m 时, 数；（3）由于dN N = r 1 − N N > 0 ,这说明N (t ) 是时间t 的单调递增函 dt mN d2N N 2 N d2N dN 1 − N ,所以当N < m 时, 2 > 0 , = r 2 1 −单增；2 N N 2 dt dt dt m m当N >Nm N d2N dN dN 时, <0, 单减,即人口增长率由增变减,在m 处最大,也就是说2 2 dt dt 2 dt在人口总数达到极限值一半以前是加速生长期,过这一点后,生长的速率逐渐变小,并且迟早会达到零,这是减速生长期；（4）用该模型检验美国从1790 年到1950 年的人口,发现模型计算的结果与实际人口在1930 年以前都非常吻合,自从1930 年以后,误差愈来愈大,一个明显的原因是在20 世纪60 年代美国的实际人口数已经突破了20 世纪初所设的极限人口.由此可见该模型的缺点之一是N m 不易确定,事实上,随着一个国家经济的腾飞,它所拥有的食物就越丰富, N m 的值也就越大；(5）用逻辑模型来预测世界未来人口总数.某生物学家估计, r = 0.029 ,又当人口总数8为3.06 × 10 时,人口每年以2%的速率增长,由逻辑模型得91 dN N = r 1 − N , N dt m即3.06 × 10 9 0.02 = 0.029 1 − NmN m = 9.86 × 10 9 ,,从而得即世界人口总数极限值近100 亿. 值得说明的是：人也是一种生物,因此,上面关于人口模型的讨论,原则上也可以用于在自然环境下单一物种生存着的其他生物,如森林中的树木、池塘中的鱼等,逻辑模型有着广泛的应用. 二、市场价格模型对于纯粹的市场经济来说,商品市场价格取决于市场供需之间的关系,市场价格能促使商品的供给与需求相等(这样的价格称为(静态)均衡价格).也就是说,如果不考虑商品价格形成的动态过程,那么商品的市场价格应能保证市场的供需平衡,但是,实际的市场价格不会恰好等于均衡价格,而且价格也不会是静态的,应是随时间不断变化的动态过程. 例 3 试建立描述市场价格形成的动态过程的数学模型解假设在某一时刻t ,商品的价格为p (t ) ,它与该商品的均衡价格间有差别,此时,存在供需差,此供需差促使价格变动.对新的价格,又有新的供需差,如此不断调节,就构成市场价格形成的动态过程,假设价格p (t ) 的变化率dp 与需求和供给之差成正比,并记f ( p, r ) dt为需求函数, g ( p ) 为供给函数（r 为参数）,于是dp = α [ f ( p, r ) − g ( p )] ， dt p ( 0) = p 0 ，其中p 0 为商品在t = 0 时刻的价格, α 为正常数. 若设f ( p, r ) = −ap + b , g ( p ) = cp + d ,则上式变为dp = −α (a + c) p + α (b −d ) ， dt p ( 0) = p 0 ，其中a, b, c, d 均为正常数,其解为①b − d −α ( a +c )t b −d . p (t ) = p 0 − +e a+ca+c9下面对所得结果进行讨论: （1）设p 为静态均衡价格,则其应满足f ( p, r ) −g ( p ) = 0 ,即于是得p =−ap +b = cp + d ,b−d ,从而价格函数p (t ) 可写为a+cp(t ) = ( p0 − p)e −α ( a +c )t + p ,令t → +∞ ,取极限得t → +∞lim p (t ) = p这说明,市场价格逐步趋于均衡价格.又若初始价格p 0 = p ,则动态价格就维持在均衡价格p 上,整个动态过程就化为静态过程；（2）由于dp = ( p − p 0 )α (a + c)e −α ( a + c ) t , dt所以,当p 0 > p 时,dp dp < 0 , p (t ) 单调下降向p 靠拢；当p0 < p 时, > 0 , p (t ) 单调增dt dt加向p 靠拢.这说明：初始价格高于均衡价格时,动态价格就要逐步降低,且逐步靠近均衡价格;否则,动态价格就要逐步升高.因此,式①在一定程度上反映了价格影响需求与供给,而需求与供给反过来又影响价格的动态过程,并指出了动态价格逐步向均衡价格靠拢的变化趋势. 三、混合溶液的数学模型例 4 设一容器内原有100L 盐,内含有盐10kg,现以3L/min 的速度注入质量浓度为0.01kg/L 的淡盐水,同时以2L/min 的速度抽出混合均匀的盐水,求容器内盐量变化的数学模型. 解设t 时刻容器内的盐量为x(t ) kg,考虑t 到t + dt 时间内容器中盐的变化情况,在dt 时间内容器中盐的改变量= 注入的盐水中所含盐量－抽出的盐水中所含盐量容器内盐的改变量为dx ,注入的盐水中所含盐量为0.01 × 3dt , t 时刻容器内溶液的质量浓度为x(t ) ,假设t 到t + dt 时间内容器内溶液的质量浓度不变（事实上,容器内100 + (3 − 2)t的溶液质量浓度时刻在变,由于dt 时间很短,可以这样看）.于是抽出的盐水中所含盐量为x(t ) 2dt ,这样即可列出方程100 + (3 − 2)t10dx = 0.03dt −即2x dt , 100 + tdx 2x = 0.03 −. dt 100 + t 又因为t = 0 时,容器内有盐10 kg,于是得该问题的数学模型为2x dx + dt 100 + t = 0.03， x(0) = 10 ,这是一阶非齐次线性方程的初值问题,其解为x(t ) = 0.01(100 + t ) +9 × 10 4 . (100 + t ) 2下面对该问题进行一下简单的讨论,由上式不难发现: t 时刻容器内溶液的质量浓度为x(t ) 9 × 10 4 p (t ) = = 0.01 + , 100 + t (100 + t ) 3且当t → +∞ 时, p (t ) → 0.01 ,即长时间地进行上述稀释过程,容器内盐水的质量浓度将趋于注入溶液的质量浓度. 溶液混合问题的更一般的提法是：设有一容器装有某种质量浓度的溶液,以流量V1 注入质量浓度为C1 的溶液（指同一种类溶液,只是质量浓度不同）,假定溶液立即被搅匀,并以V2 的流量流出这种混合溶液,试建立容器中质量浓度与时间的数学模型.首先设容器中溶质的质量为x (t ) ,原来的初始质量为x 0 , t =0 时溶液的体积为V2 ,在d t 时间内,容器内溶质的改变量等于流入溶质的数量减去流出溶质的数量,即dx = C1V1dt − C 2V2 dt ,其中C1 是流入溶液的质量浓度, 度, C 2 =C2 为t 时刻容器中溶液的质量浓x ，于是,有混合溶液的数学模型V0 + (V1 − V2 )tdx = C1V1 − C2V2 ， dt x(0) = x0 ．该模型不仅适用于液体的混合,而且还适用于讨论气体的混合. 四、振动模型11振动是生活与工程中的常见现象.研究振动规律有着极其重要的意义.在自然界中,许多振动现象都可以抽象为下述振动问题. 例 5 设有一个弹簧,它的上端固定,下端挂一个质量为m 的物体,试研究其振动规律. 假设（1）物体的平衡位置位于坐标原点,并取x 轴的正向铅直向下（见图4）.物解体的平衡位置指物体处于静止状态时的位置.此时,作用在物体上的重力与弹性力大小相等, 方向相反；（2）在一定的初始位移x0 及初始速度v0 下,物体离开平衡位置,并在平衡位置附近作没有摇摆的上下振动；（3）物体在t 时刻的位置坐标为x = x (t ) ,即t 时刻物体偏离平衡位置的位移；（4）在振动过程中,受阻力作用.阻力的大小与物体速度成正比,阻力的方向总是与速度方向相反,因此阻力为−hdx , h 为阻尼系数；（5）当质点有位移x (t ) 时,假设所受dt的弹簧恢复力是与位移成正比的,而恢复力的方向总是指向平衡位置,也就是总与偏离平衡（6）在振动过程中位置的位移方向相反,因此所受弹簧恢复力为− kx ,其中k 为劲度系数；受外力f (t ) 的作用.在上述假设下,根据牛顿第二定律得m这就是该物体的强迫振动方程.d2x dx = − h − kx + f ( x) , 2 dt dt①由于方程①中, f (t ) 的具体形式没有给出,所以,不能对式①直接求解.下面我们分四种情形对其进行讨论. 1. 无阻尼自由振动在这种情况下,假定物体在振动过程中,既无阻力、又不受外力作用.此时方程①变为图4O xd2x m 2 + kx = 0 , dt令k = ω 2 ,方程变为md2x +ω2x = 0, 2 dt特征方程为特征根为通解为或将其写为λ2 + ω 2 = 0 ,λ1, 2 = ±iω ,x = C 1 sin ω t + C 2 cos ω t ,C 12 + C 22 cos ω tx=C1 C 12 + C 22sin ω t +C2 C 12 + C 2212= A (cos ϕsin ω t + sin ϕcos ω t )= A sin( ω t + ϕ )其中，C2 C +C2 1 2 2A = C12 + C 22 , sin ϕ =, cos ϕ =C12 C + C2 2 1.2 这就是说,无阻尼自由振动的振幅A = C12 + C 2 ,频率ω =k 均为常数. m2.有阻尼自由振动2.有阻尼自由振动在该种情况下,考虑物体所受到的阻力,不考虑物体所受的外力.此时,方程①变为m令d2x dx + h + kx = 0 , 2 dt dtk h = ω 2 , = 2δ ,方程变为m m d2x dx + 2δ +ω2x = 0, 2 dt dt特征方程为λ + 2δλ + ω = 0 ,特征根2 2λ1, 2 = −δ±δ 2 −ω 2 .根据δ 与ω 的关系,又分为如下三种情形：(1)大阻尼情形, δ > ω .特征根为二不等实根,通解为x = C 1e ( −δ +δ 2 −ω 2 ) t+ C 2 e ( −δ +δ 2 −ω 2 ) t(2)临界阻尼情形, δ = ω .特征根为重根,通解为x = (C 1 + C 2t ) e−δt这两种情形,由于阻尼比较大,都不发生振动.当有一初始扰动以后,质点慢慢回到平衡位置,位移随时间 t 的变化规律分别如图 5 和图 6 所示.xxx0x0OtOt图5图6(3)小阻尼情形, δ < ω .特征根为共轭复根,通解为13x = e −δt ( C1 sin ω 2 −δ 2 t + C 2 sin ω 2 −δ 2 t ) 将其简化为x = Ae −δt sin( ω 2 −δ 2 t + ϕ )其中 A = C1 2 + C 2 2 , sin ϕ =C2 C1 + C 22 2, cos ϕC1 C1 + C 22 2−δt , 振幅 A e 随时间 t 的增加而减小.因此,这是一种衰减振动.位移随时间 t 的变化规律见图 7. 3.无阻尼强迫振动 3.无阻尼强迫振动在这种情形下,设物体不受阻力作用,其所受外力为简谐力 f (t ) = m sin pt ,此时,方程①化为xx0md2x + kx = m sin pt , dt 2d2x + ω 2 x = sin pt , dt 2根据 ip 是否等于特征根 iω,其通解分为如下两种情形：（1）当 p ≠ ω时,其通解为Ot图7x=1 sin pt + C 1 sin ω t + C2 cos ω t , ω− p22此时,特解的振幅1 为常数,但当 p 接近于ω时,将会导致振幅增大,发生类似共振的ω− p22现象；（2）当 p = ω时,其通解为x = −1 t cos pt + C 1 sin ω t + C2 cos ω t , 2p此时,特解的振幅1 t 随时间 t 的增加而增大,这种现象称为共振,即当外力的频率 p 等于 2p物体的固有频率ω时,将发生共振. 4.阻尼强迫振动 4.阻尼强迫振动在这种情形下,假定振动物体既受阻力作用,又受外力 f ( x ) = m sin pt 的作用,并设δ < ω ,方程①变为d2x dx + 2δ + ω 2 x = sin pt , 2 dt dt14特征根λ = −δ ± iω 2 −δ 2 , δ ≠ 0 ,则 ip 不可能为特征根,特解为x * = A sin pt + B cos pt ,其中 A =ω 2 − p2(ω− p ) + 4δ p2 2 2 2 2,B =(ω2− 2δ p − p 2 ) 2 + 4δ2p2,还可将其化为x* =1 [(w2 − p2 )sin pt − 2δ p cos pt ] , 2 2 (w − p ) + 4δ p2 2 2由此可见,在有阻尼的情况下,将不会发生共振现象,不过,当 p = ω时, x* = −1 cos pt , 2δp若δ很小,则仍会有较大的振幅；若δ比较大,则不会有较大的振幅.15。