第二章 人(虫)口-资源-环境动力模型 2007
2.1 地学建模的动力学机制分析 动力学建模理解系统变化的原因和机制,即在内部和外界因素作用下,系统有关的物理量(要素)将如何变化。对于变速运动: (2.1) 方程(2.1)是牛顿第二定律,描写质量为m的物体在n个合力的作用下其运动的变化情况。如果考虑物体只受到动力F和阻力f的作用,并假设物体的质量为1个单位,取m=1, (2.1)可写为: (2.2)
(2.3) 当用x表示任意物理量得到普适的一元(维)动力学方程: (2.3) x-广义物理量(如速度,温度,人口密度,生物量,经济总量等);F(/f)分别为导致(/阻碍)其变化的广义动力(∕阻力)因素; 建模对导致系统发生变化的有关广义动力和广义阻力进行数学表达。 ●广义动力是指有利于系统变化的所有因素(内部和外界的),它是正定的。 ●广义阻力则是指不利于系统变化的所有因素(内部和外界的),它是负定的。
1.1 f=0 考察任一不受阻力影响的一维(元)地学动力系统,如人口、经济总量、温度的变化等。如果广义动力为一正定的常数(c>0),则: (2.4) x=ct,t→∞,x→∞ 此时系统没有定态,x将随时间的增长而无限增大。不存在定态的模式系统不能真实描写客观系统。 模式(2.4)只能是一种极端近似下的模型。它代表了随时间线性变化的增长模式,如自增长模式,可再生模式,线性增值模式,线性扩张模式。
此系统有一个定态x=0,因a>0,定态不稳定。 如果广义动力为一正定的线性(正比)动力 (2.5) 则: (2.6) 此系统有一个定态x=0,因a>0,定态不稳定。 客观世界存在矛盾的两方面,即客观系统总是存在着两个以上的定态。只存在一个定态的模式系统也是不能真实描写客观系统的,充其量只能研究系统在某一阶段或短时间内的演化行为。由于真实的物理态都是稳定的,不稳定态是过渡态,只说明了系统演化的起点。
由于(2. 6)的解的形式为指数形,所以x将随时间的增长以指数形式快速地无限增大,所以模式(2
如果广义动力为一正定的2次非线性(正比)动力: (2.7) 则: (2.8) (2.8)是一简单的非线性动力系统。定态方程是一个一元二次方程,故系统有两个定态。非线性模式系统较线性模式系统的最大优点就是能够描写客观世界的双态问题。 凡是研究对象客观地存在双态,则所建立的模式系统至少必须是2次非线性的。如果有关参数是正定的,(2.8)的解析解是以反比形式快速地趋近某一定值。
以上的分析表明: ♥ 常数形动力系统不存在定态,系统的演化特性是随时间成正比增长的直线运动; ♥ 常数形动力系统不存在定态,系统的演化特性是随时间成正比增长的直线运动; ♥ 线性动力系统只能研究单态客观系统,系统的演化特性是随时间变化作指数增长的曲线运动; ♥ 2次非线性动力系统可研究双态或单态的客观系统,系统的演化特性是随时间变化作反比增长函数形的曲线运动; 只存在动力的模式系统是不稳定的。状态的变化意味着系统的突变,而环境的突变就意味着灾难。在多数情况下,客观状态都是稳定的。否则人类就无法生存。 所要建立的地学模式,必须具有稳定的定态的。即模式系统必须具有阻力项。
1.2 F=0 类似地,对只有阻力的模式动力系统,有: ♥ 常数形阻力系统不存在定态,系统的演化特性是随时间成正比衰减的直线运动; ♥ 线性阻力系统只能研究单态客观系统,系统演化特性是随时间变化作指数衰减的曲线运动; ♥ 2次非线性阻力系统可研究双态或单态的客观系统,系统的演化特性是随时间变化作反比衰减函数形的曲线运动。 结合1.1和1.2的的结论,在建模前可根据对研究对象的已知信息来决定模式的数学形式,如动力项(正的多项式或单项)和阻力项(负的多项式或单项)共存,而单态或多态的系统则分别对应于一次或高次的动力系统。
1.3 耦合系统 对于二元耦合系统,存在内部的相互作用和反馈,可分为单向或双向。 单向对应于主次分明或权重差异很大的耦合系统; 如鲨鱼-小鱼耦合系统里,鲨鱼对小鱼施加了一个单方向的捕食作用,而小鱼则无法向鲨鱼施加类似的捕食作用; 双向对应于权重相差不是太大的耦合系统。 如交战的两个国家都可以给对方施加一个破坏的作用。 故在二元(及以上)的耦合系统里: ①考虑不同物理要素各自所受到的动力和阻力; ②考虑两者之间的相互作用,这种相互作用可以是常数形、线性形,非线性形的;可以是单向或双向的。
例1 (2.9) 对于模式系统(2.9),其x的线性项代表线性动力项,自生长或线性正反馈,其结果是无限增长。而-x2则代表非线性阻力项,而(p-x)则代表环境或外界的线性约束(稳定因子),这两种表示都起稳定作用。 对于任意一个数学模式,可以赋以不同的物理内涵: 研究人口与环境关系时,x代表人口密度,x的线性项代表人口生产的动力项(父母越多后代越多),p代表环境资源所能供养的最大人口密度,(p-x)代表环境容量。 就是虫口模式,其中r=a-b为人口净增长率(a,b分别为出生率和死亡率);xm为环境资源所能供养的最大人口数。
考虑到各种不同的阻力和动力,可有下述形式: (2.10) 为线性动力, 为非线性环境约束, 为非线性阻力。 其中: (2.11) (2.12) (2.13) (2.14) (2.15) (2.16) 及:
(2.17) (2.17)已经是一个普适的虫(人)口模式或普适的可再生资源模式: ① n>1:虫(人)与环境的不协调,其发展受到了环境约束。 ② 0<n<1:虫(人)与环境协调,减小环境对其发展的约束影响; ③ n<0:虫(人)与环境的和谐关系。
§2.2 Malthus模型与Logistic模型 为了保持自然资源的合理开发与利用,人类必须保持并控制生态平衡,甚至必须控制人类自身的增长。 生物种群的增长,一方面取决于种群中生物个体的数目(尤其是雌性体),另一方面受环境的制约,如营养条件等。环境承载力!!! 种群的数量本应取离散值,但由于种群数量一般较大,为建立微分方程模型,可将种群数量看作连续变量,由此引起的误差将是十分微小的。
模型1 马尔萨斯(Malthus)模型 马尔萨斯在分析人口出生与死亡情况的资料后发现,人口净增长率r基本上是一常数,(r=b-d,b为出生率,d为死亡率),即: 或 (2.1) (2.1)的解为: (2.2) 其中N0=N(t0)为初始时刻t0时的种群数。 马尔萨斯模型的一个显著特点:种群数量翻一番所需的时间是固定的,且是一无界函数。 令种群数量翻一番所需的时间为T,则有: 故
模型预测 按马氏模型计算,人口数量每34.6年增加一倍,呈几何级数增长。 按马氏模型计算,人口数量每34.6年增加一倍,呈几何级数增长。 几何级数的增长 例如,到2510年人口达2×1014,即使海洋全部变成陆地,每人也只有9.3平方英尺的活动范围,而到2670年人口达36×1015,只好一个人站在另一人的肩上排成二层了。 故马尔萨斯模型是不完善的。
模型检验 Malthus模型只反映了单纯考虑生物个体繁殖能力,而未考虑到限制因素。 比较历年的人口统计资料,可发现人口增长的实际情况与马尔萨斯模型的预报结果基本相符,例如,1961年世界人口数为30.6 (即3.06×109),人口增长率约为2%,人口数大约每35年增加一倍。检查1700年至1961的260年人口实际数量,发现两者几乎完全一致,且按马氏模型计算,人口数量每34.6年增加一倍,两者也几乎相同。 Malthus模型只反映了单纯考虑生物个体繁殖能力,而未考虑到限制因素。 实际上只有在群体总数不太大时才合理,到总数增大时,生物群体的各成员之间由于有限的生存空间,有限的自然资源及食物等原因,就可能发生生存竞争等现象。 所以Malthus模型假设的人口净增长率不可能始终保持常数,它应当与人口数量有关。
模型2 Logistic模型 人口净增长率应当与人口数量有关,即: r=r(N) 从而有: (2.3)
对马尔萨斯模型引入一次项(竞争项) 令 r(N)=r-aN 此时得到微分方程: 或 (2.4) (2.4)被称为Logistic模型或生物总数增长的统计筹算律,是由荷兰数学生物学家弗赫斯特(Verhulst)首先提出。一次项系数是负的,因为当种群数量很大时,会对自身增大产生抑制性,故一次项又被称为竞争项。 (2.4)可改写成: (2.5)
(2.5) 由于空间和资源的约束,不可能供养无限增长的种群个体,当种群数量过多时,由于人均资源占有率的下降及环境恶化、疾病增多等原因,出生率将降低而死亡率却会提高。 设环境能供养的种群数量的上界为K(近似地将K看成常数),N表示当前的种群数量,K-N恰为环境还能供养的种群数量,种群增长率与两者的乘积成正比,正好符合统计规律,得到了实验结果的支持,(2.5)也被称为统计筹算律。
对(2.5)分离变量: 两边积分并整理得: 令N(0)=N0,求得: 故(2.5)的满足初始条件N(0)=N0的解为: (2.6) 易见: N(0)=N0 , N(t)的图形
模型检验 用Logistic模型来描述种群的增长效果较好。数学生物学家高斯(E·F·Gauss)把5只草履虫放进一个盛有0.5cm3营养液的小试管,发现开始时草履虫以每天230.9%的速率增长,此后增长速度不断减慢,到第五天达到最大量375个,实验数据与r=2.309,a=0.006157,N(0)=5的Logistic曲线: 几乎完全吻合
Malthus模型和Logistic模型的总结 (2.3) 两模型均为对微分方程(2.3)所作的模拟近似方程。前一模型假设了种群增长率r为一常数(r被称为该种群的内禀增长率)。后一模型则假设环境只能供养一定数量的种群,从而引入了一个竞争项。 用模拟近似法建立微分方程来研究实际问题,必须对解进行检验,是否与实际情况相符或基本相符。相符性越好则模拟得越好,否则就要找出原因,对模型进行修改。 两个模型虽然都是为了研究种群数量的增长情况而建立的,但它们也可用来研究其他实际问题,只要这些实际问题的数学模型有相同的微分方程即可。
稳定性问题 在研究许多实际问题时,人们最关心的也许并非系统与时间有关的变化状态,而是系统最终的发展趋势。 例如,在研究某频危种群时,虽然也想了解它当前或今后的数量,但我们更为关心的却是它最终是否会绝灭,用什么办法可以拯救这一种群,使之免于绝种等问题。 微分方程或微分方程组的稳定性理论。
一般的微分方程或微分方程组可以写成: 定义:称微分方程或微分方程组 (2.7) 为自治系统或动力系统。 若方程或方程组f(x)=0有解Xo,X=Xo显然满足(2.7)。称点Xo为微分方程或微分方程组(2.7)的平衡点或奇点。
例1 对Logistic模型 共有两个平衡点:N=0和N=K,分别对应微分方程的两个特殊解。前者为No=0时的解而后者为No=K时的解。 当No<K时,积分曲线N=N(t)位于N=K的下方; 当No>K时,则位于N=K的上方。 从图中不难看出,若No>0,积分曲线在N轴上的投影曲线(称为轨线)将趋于K。这说明,平衡点N=0和N=K有着极大的区别。
定义2 设x0是(2.7)的平衡点,称: (1)x0是稳定的,如果对于任意的ε>0,存在一个δ>0,只要|x(0)- x0|<δ,就有|x(t)- x0|<ε对所有的t都成立。 (2)x0是渐近稳定的,如果它是稳定的且 (3)x0是不稳定的,如果(1)不成立。 根据这一定义,Logistic方程的平衡点N=K是稳定的且为渐近稳定的,而平衡点N=0则是不稳定的。
§2.3 捕食系统的Volterra方程 意大利生物学家D’Ancona曾致力于鱼类种群相互制约关系研究,无意中发现了一些第一次世界大战期间地中海沿岸港口捕获的几种鱼类占捕获总量百分比的资料,发现各种软骨掠肉鱼,如鲨鱼、鳐鱼等鱼类(称之为捕食者或食肉鱼)占总渔获量的百分比有明显的增加。 1914-1923年意大利阜姆港收购鱼中食肉鱼所占比例 年代 1914 1915 1916 1917 1918 百分比 11.9 21.4 22.1 21.2 36.4 1919 1920 1921 1922 1923 27.3 16.0 15.9 14.8 10.7
捕获的各种鱼的比例近似地反映了地中海里各种鱼类的比例。战争期间捕鱼量大幅下降,但捕获量的下降为什么会导致鲨鱼、鳐鱼等食肉鱼比例的上升,即对捕食者有利而不是对食饵有利呢? 数学家V.Volterra(意大利)进行了数学建模。
模型假设 (Ⅰ)相互关系中仅有一种捕食者与一种猎物; (Ⅱ)如果捕食者数量下降到某一阈值以下,猎物数量就上升,而捕食者数量如果增多,猎物种数量就下降; (Ⅲ)如果猎物数量上升到某一阈值,捕食者数量增多,而猎物数量如果很少,捕食者数量就下降。 这一简单的模型做了一个有趣的预测:捕食者和猎物种群会发生动态循环,就象在自然的捕食者-猎物种群动态中所观察到的。
1、模型建立 Volterra将鱼划分为两类。一类为食用鱼(食饵),数量记为x1(t),另一类为食肉鱼(捕食者),数量记为x2(t),并建立双房室系统模型。 对于食饵(Prey)系统 : 大海中有食用鱼生存的足够资源,可假设食用鱼独立生存将按增长率为r1的指数律增长(Malthus模型),即设: 由于捕食者的存在,食用鱼数量减少,设减少的速率与两者数量的乘积成正比(竞争项的统计筹算律),即: λ1反映了捕食者掠取食饵的能力
对于捕食者(Predator)系统 : 捕食者设其离开食饵独立存在时的死亡率为r2,即: 但食饵提供了食物,使生命得以延续。这一结果也要通过竞争来实现,再次利用统计筹算律,得到: 综合以上分析,建立P-P模型(Volterra方程)的方程组: (2.10) 没有人工捕获的自然环境中食饵与捕食者间的相互制约关系
2、模型分析 方程组(2.10)是非线性的,不易直接求解。该方程组共有两个平衡点,即: 和 Po(0,0)是平凡平衡点且不稳定,故不研究 x1、x2轴是方程组的两条相轨线。 方程组还有两组平凡解: 和 当x1(0)、x2(0)均不为零时, ,应有x1(t)>0且x2(t)>0,相应的相轨线应保持在第一象限中。
求(2.10)的相轨线 将两方程相除消去时间t,得: 分离变量并两边积分得轨线方程: 令 用微积分知识容易证明: 有: 同理:对 有: (2.11) 令 两者应具有类似的性质 用微积分知识容易证明: 有: 同理:对 有:
与 的图形见图2-4 图2-4 (b) 图2-4 (a)
将Volterra方程中的第二个改写成: 等式左端为零,故可得: 平衡点P的两个坐标恰为食用鱼与食肉鱼在一个周期中的平均值 同理:
解释D’Ancona发现的现象 引入捕捞能力系数e,(0<e<1),e表示单位时间 内捕捞起来的鱼占总量的百分比。故Volterra方程应为: 平衡点P的位置移动到: A(0,0)
食用鱼的数量反而因捕捞它而增加?由于捕捞能力系数e的引入,食用鱼的平均量有了增加,而食肉鱼的平均量却有所下降,e越大,平衡点的移动也越大。 根据P-P模型,可以得出: (1)食用鱼的平均量取决于参数r1与λ1 (2)食用鱼繁殖率r1的减小将导致食肉鱼平均量的减小,食肉鱼捕食能力λ1的增大也会使自己的平均量减小;反之,食肉鱼死亡率r2的降低或食饵对食肉鱼供养效率λ2的提高都将导致食用鱼平均量的减少。 (3)捕鱼对食用鱼有利而对食肉鱼不利,多捕鱼(当然要在一定限度内,如e<r1)能使食用鱼的平均数量增加而使食肉鱼的平均数量减少。
Volterra的模型揭示了双种群之间内在的互相制约关系,成功解释了D’Ancona发现的现象。 ●P-P模型导出的结果在一定程度上是符合客观实际的,有着广泛的应用前景。 ●如当农作物发生病虫害时,使用杀虫剂在杀死害虫的同时也可能杀死害虫的天敌,(害虫与其天敌构成一个双种群捕食系统),导致使用杀虫剂的结果会适得其反。
§2.4 较一般的双种群生态系统 对捕食系统中存在周期性现象的结论,大多数生物学家并不完全赞同,因为更多的捕食系统并没有这种特征。 §2.4 较一般的双种群生态系统 对捕食系统中存在周期性现象的结论,大多数生物学家并不完全赞同,因为更多的捕食系统并没有这种特征。 一个捕食系统的数学模型未必适用于另一捕食系统,捕食系统除具有共性外,往往还具有本系统特有的个性,反映在数学模型上也应当有所区别。现考察较为一般的双种群系统。
一般的双种群系统 用x1(t)和x2(t)记t时刻的种群量(也可以是种群密度), 设 Ki为种群i的净相对增长率。 Ki随种群不同而不同,也随系统状态的不同而不同,即Ki为x1、x2的函数。Ki具体的函数形式,采用线性化方法,可得到下面的微分方程组: (2.12) (2.12)不仅可以用来描述捕食系统。也可以用来描述相互间存在其他关系的种群系统。
(2.12)式的一些说明: 式中a1、b2为本种群的亲疏系数,a2、b1为两种群间的交叉亲疏系数。a2b1≠0时,两种群间存在着相互影响,此时又可分为以下几类情况: (Ⅰ)a2>0,b1>0,共栖系统。 (Ⅱ)a2<0,b1>0( 或a2>0,b1<0 ),捕食系统。 (Ⅲ)a2<0,b1<0,竞争系统。 (Ⅰ)—(Ⅲ)构成了生态学中三个最基本的类型,种群间较为复杂的关系可以由这三种基本关系复合而成。 (2.12)
(2.12)是否具有周期解? 不同的系统具有不同的系数,在未得到这些系数之前先进行一般化的讨论。 首先,系统的平衡点为方程组: 的解。 (2.13) 的解。 均为平凡平衡点。 如果系统具有非平凡平衡点 则它应当对应于方程组 的根
定理 3 P存在时,P一般是稳定平衡点,此时平凡平衡点常为不稳定的鞍点。 解得: (无圈定理)若方程组(2.12)的系数满足 (i) A=a1b2-a2b1≠0 (ii)B= a1b0(a2-b2)-a0b2(a1-b1)≠0 则(2.12)不存在周期解 定理 3 证明: 记 作函数 ,并记 f(x1,x2)=x1(a0+a1x1+a2x2), g(x1,x2)=x2(b0+b1x1+b2x2),容易验证: 假设结论不真,则在x1~x2平面第一象限存在(2.13)的一个圈Γ,它围成的平面区域记为R。
对于Voltera方程,由a1=b2=0,得B=0;所以无圈定理不适用于Volterra方程。 于是由K(x1,x2)>0且连续以及AB≠0可知,函数 在第一象限中不变号且不为零,故二重积分: (2.14) 但另一方面,由格林公式 对于Voltera方程,由a1=b2=0,得B=0;所以无圈定理不适用于Volterra方程。 注意到 , ,又有: (3.36) (2.15) 其中T为周期。 (2.14)与(2.15)矛盾,说明圈Γ不可能存在。
对于一般的生态系统,如果通过求解的微分方程来讨论常常会遇到困难。 怎样讨论一般的生态系统? 对于一般的生态系统,如果通过求解的微分方程来讨论常常会遇到困难。 此时可以研究种群的变化率,搞清轨线的走向来了解各种群数量的最终趋势。
简化模型,设竞争系统的方程为: 定理4 :(竞争排斥原理) 其中:αβ不为0,否则为Logistic模型。 方便讨论取α=β=1,所用方法可适用一般情况。 定理4 :(竞争排斥原理) 若K1>K2,则对任一初态(x1(0),x2(0)),当t→+∞时,总有(x1(t),x2(t))→(K1,0),即物种2将绝灭,而物种1则趋于环境允许承担的最大总量。
作直线l1: x1+x2=K1及l2: x1+x2=K2, K1> K2,见图2-6。 有以下几个引理: 引理1 若初始点位于区域I中,则解(x1(t)、x2(t))从某一时刻起必开此区域而进入区域II 引理2 若初始点(x1(0)、x2(0))位于 区域II中,则(x1(t),x2(t))始终位于II中,且: dx1/dt<0 dx2/dt<0 图2-6 III II I k1 k2 dx1/dt>0 dx2/dt>0 引理3 若初始点位于区域III中,且对于任意t ,(x1(t),x2(t))仍位于 III中,则当t→+∞时,(x1(t), x2(t))必以(K1,0)为极限点。
易见只有上述三种可能,而在三种可能情况下(x1(t),x2(t))均以(K1,0)为极限,定理得证。 定理4的证明: 由引理1和引理2,初始点位于像限I和II的解必趋于平衡点(K1,0)。 由引理3,初始点位于III且(x1(t),x2(t))始终位于III中的解最终必趋于平衡点(K1,0),而在某时刻进入区域II的解由引理最终也必趋于(K1,0)。 易见只有上述三种可能,而在三种可能情况下(x1(t),x2(t))均以(K1,0)为极限,定理得证。
当无法求得解析解时,往往需要用数值解法。1999年美国大学生数学建模竞赛题A题(小行星撞击地球):如何描述南极地区的生态系统,如何定量研究小行星撞击地球对南级生态环境的影响? 将南极附近的生物划分成:海藻、鳞虾和其他海洋生物。鳞虾吃海藻,其他海洋动物吃鳞虾,可建立了一个三房室系统模型。 小行星撞击会影响大气层能见度,从而影响海藻的生长(光合作用),进而影响生物链中的其他生物。但此时无法得到模型中的参数值(小行星撞击南极的事件并未发生过),为此可取一系列不同的参数值,对不同参数值下模型的数值解进行了分析对比,研究了解对各参数变化的灵敏度。