常微分方程建模方法

在应用数学方法解决实际问题的过程中，很多时候，要直接导出变量之间的函数关系较为困难，但要导出包含未知函数的导数或微分的关系式却较为容易，在这种情况下，就需要我们建立微分方程模型来研究。事实上，微分方程是研究函数变化规律的有力工具，在物理、工程技术、经济管理、军事、社会、生态、环境、人口、交通等各个领域中有着广泛的应用.下面我们就介绍如何应用微分方程模型来解决实际问题.

利用微分方程解决的问题通常可以分为两类：一类问题要求把未知变量直接表示为已知量的函数，这时，有些问题可以求出未知函数的解析表达式，在很多情况下只能利用数值解法；另一类问题只要求知道未知函数的某些性质，或它的变化趋势，这时可以直接根据微分方程定性理论来研究.

2.1 微分方程的一般理论

2.1.1微分方程简介

所谓微分方程就是表示未知函数、未知函数的导数与自变量之间的关系的方程若未知函数是一元函数的微分方程 叫常微分方程而未知函数是多元函数的微分方程 叫偏微分方程 例如

y44y'''10y''12y'5ysin2x （2.1.1） x2y''12xy'5y0 (y')2xy0

（2.1.2） （2.1.3） （2.1.4） （2.1.5） （2.1.6）

2y''y'xy0 y(n)10

uta2uxx

其中，方程（2.1.6）是偏微分方程，其他都是常微分方程

微分方程中所出现的未知函数的最高阶导数的阶数 叫微分方程的阶例如，方程（2.1.1）是四阶微分方程，（2.1.3）是一阶微分方程.一般n阶微分方程具有形式

F(x y y     y(n) )0

或

y(n)f(x y y     y(n1) ) 

必须指出，y(n)是必须出现的，而x,y,y',,y(n1)等变量则可以不出现，如方程（2.1.5）.

若F(x y y     y(n) )是关于y及其各阶导数的线性函数，则称此方程是线性的，否则，称为非线性的.例如，方程（2.1.1）、（2.1.2）是线性微分方程，方程（2.1.3）是非线性微分方程.

线性微分方程可以分为常系数和变系数两大类，常系数线性微分方程中未知函数及其导数的系数均为常数，而变系数线性微分方程中未知函数及其导数的系数不完全是常数.例如，方程（2.1.1）、（2.1.5）是常系数线性微分方程，而方程（2.1.2）是变系数线性微分方程.

满足微分方程的函数(也就是，把函数代入微分方程能使该方程成为恒等式)叫做该微分方程的解 确切地说 设函数y(x)在区间I上有n阶连续导数 如果在区间I上

F[x (x) (x)    (n) (x)]0

那么函数y(x)就叫做微分方程F(x y y    y(n) )0在区间I上的解

如果微分方程的解中含有任意常数 且任意常数的个数与微分方程的阶数相同 这样的解叫做微分方程的通解若以显函数形式给出的解，称为显式解，以隐函数形式给出的解，就称为隐式解.

为了确定微分方程的一个特定的解，我们通常给出这个解所必须满足的条件，这就是所谓的定解条件。常见的定解条件是初始条件。用于确定通解中任意常数的条件 称为初始条件 如

xx0 时 yy0  y y0 

一般写成

 yxx0y0 yxx0y0求微分方程满足初始条件的特解这样一个问题，叫做微分方程的初值问题，例如求微分方程yf(x y)满足初始条件yxx0y0的解的问题 记为

yf(x,y) yyxx00微分方程的解的图形是一条曲线 叫做微分方程的积分曲线 2.1.2微分方程初值问题的适定性

在实际问题中，由于自然界本身就给出了问题唯一的答案，所以一个初值问题提得是否符合实际情况，从数学角度来看，可以从三个方面加以检验：

（1）解的存在性，即初值问题是否有解？ （2）解的唯一性，即初值问题的解是否只有一个？

（3）解的稳定性，即当初值条件有微小变动时，解相应的只有微小的变动.

一个初值问题的解如果满足存在性、唯一性和稳定性，则称此初值问题是适定的.微分方程的初值问题解的适定性具有重要的实际意义.微分方程模型通常是用来描述确定性的模型的.对于一个由实际问题所建立的微分方程模型，如果其初值问题的解不存在，或解不唯一，这样的模型本身就是不合理的，是没有实际意义的.因为在一定的条件下实际问题到最后总会有确定的结果，这反映在模型上，就是定解问题有唯一解.而解的稳定性更是具有重要的实际应用背景.由于由实

际问题导出初值问题时，总要经过一些简化、近似的过程以及一些附加的假设，并且在测量初始条件的值和测量方程中各项系数（或参数）等的值时，不可避免地会出现测量误差，从而致使我们得到的微分方程模型，通常只能是近似地描述所讨论的实际问题，难免存在误差.

当测量的数据出现微小的误差时，相应模型的“解”是否也只有微小的误差？如果回答是肯定的，我们就说这个模型的解（在某种意义下）是稳定的，否则，就说这个模型的解是不稳定的.显然，只有“稳定的”解才具有可靠性，只有“稳定的”解才会有使用价值.相反，“不稳定的”解是不会有任何使用价值的.因为初值、参数等的微小误差或干扰将导致“差之毫厘，谬以千里”的严重后果.同时，稳定性也是计算机利用数值方法求解的前提和保证.

2.2 微分方程的平衡点及稳定性

一个微分方程即使存在解，也有可能解不出.事实上，我们在学习高等数学的时候就知道，能用初等的方法求出解的微分方程只是极少数.更多的情况下，是没有初等解法的，这一事实为法国数学家刘维尔(Liouville)在1841年所证明.如果一个微分方程的解不是一个初等函数，由于我们不能将方程的解函数像初等函数一样地将它表示出来，也就可能出现方程解不出的情况.既然初等积分法有着不可克服的局限性，那么是否可以不求微分方程的解，而是从微分方程本身来推断其解的性质呢?定性理论和稳定性理论正是在这种背景下发展起来的.前者由法国数学家庞加莱(Poincaré，1854-1912)在19世纪80年代所创立，后者由数学家李雅普罗夫(Liapunov，1857-1918)在同年代所创立.它们共同的特点就是在不求出方程的解的情况下，直接根据微分方程本身的结构和特点，来研究其解的性质.由于这种方法的有效性，近一百多年以来它们已经成为常微分方程发展的主流，并且在实际中有大量的应用.比如，在研究许多实际问题时，其变量的变化率仅与平衡状态有关而与时间并无直接的联系，或者人们最为关心的并非系统与时间有关的变化状态，而是系统最终（时间充分大之后）的发展趋势.例如，在研究某频危种群时，虽然我们也想了解它当前或今后的数量，但我们更为关心的却是它最终是否会绝灭，用什么办法可以拯救这一种群，使之免于绝种等等问题.要解决这类问题，就需要用到微分方程或微分方程组的稳定性理论.本节对定性理论和稳定性理论的一些基本概念和基本方法作一简单介绍.

2.2.1一阶方程的平衡点及稳定性

函数的变化率只和函数本身有关而与自变量无关的微分方程或

微分方程组被称为自治系统，也称为动力系统.

通常，一阶微分方程可写成

dxf(t,x)，而自治系统则可写成dtdxf(x)，即右端不显含自变量t. dt方程f(x)0的实根xx0称为自治系统

dxf(x)的平衡点（或奇dt点）.显然，根据平衡点的定义，xx0也是自治系统的一个解（奇解），即微分方程不变化的解，也就是常数解.

如果对任意给定的0，存在0（一般与和t0有关），使得只要初始条件x(t0)满足x(t0)x0时，自治系统满足

x(t)x0(对所有的tt0)

dxf(x)的解x(t)均dt则称自治系统

dxf(x)的平衡点x0是（在李雅普洛夫意义下）稳定的. dtdx如果自治系统f(x)的平衡点x0稳定，且存在这样的00使当

dtx(t0)x00时，自治系统的解x(t)都满足

limx(t)x0

t则称平衡点x0是（在李雅普洛夫意义下）渐近稳定的；否则，称x0是不稳定的.特别的，如果从所有可能的初始条件出发，都是是渐近稳定的，则称平衡点x0是全局渐近稳定的.我们在这里讨论的稳定性都是指渐近稳定性.

判断平衡点x0是否稳定通常有两种方法： （1）间接法：求出自治系统的定义判断；

（2）直接法：不求自治系统

dxf(x)的解x(t)，按线性近似判定dtdxf(x)的解x(t)，利用上述稳定性dt稳定性，即利用f(x)在xx0处的泰勒展开式，只取一次项，

f(x)f'(x)(xx0),则方程

dxf(x)近似为： dtdx f'(x0)(xx0) （2.2.1）

dtdx方程（2.2.1）称为方程f(x)的近似线性方程.显然，xx0也

dt是近似线性方程（2.2.1）的平衡点.

因为方程（2.2.1）的通解为

x(t)cef'(x0)tx0 （2.2.2）

其中c是由初始条件决定的常数.由稳定性的定义很容易证明：

若f'(x0)0，则x0是方程（2.2.1）的稳定的平衡点； 若f'(x0)0，则x0是方程（2.2.1）的不稳定的平衡点. 同样，根据李雅普洛夫理论，对于自治系统

dx若f'(x0)0，f(x)，

dt则x0是稳定的平衡点；若f'(x0)0，则x0是不稳定的平衡点. 例2.2.1 讨论微分方程

dyy2y的平衡点稳定性. dt解 I 间接法：易知，方程f(y)y2y0有两个常数解

y1(t)0和y2(t)1

这也是原微分方程的两个平衡点. 当y0和y1时，原方程可写成

dydt y(y1)解得

ln|y|ln|y1|tc

即原方程的通解为

y1 1cet1， y0若有初始条件y(0)y0(y00,1)，求得

c1那么所给初值问题的解是

y1111ety0

易知

limy(t)limtt1111ety00

根据稳定性定义，x0是稳定的平衡点，x1是不稳定的平衡点.

II直接法：由于

f'(y)2y1

f'(0)10,f'(1)10，从而，x0是稳定的平衡点，x1是不稳定

的平衡点.

在微分方程模型中，微分方程解的这种特性对许多实际问题的讨论是非常重要的.例如，

研究对象为某温度控制系统.我们有一个理想温度x和一个实际温度y，x和y都是时间t的函数，而x,y满足某个微分方程，假如我们能够设定一个控制器，使得x和y的关系更接近我们的需求，那么保证这个控制器稳定就是一个非常重要的前提.我们以空调为例，假设室内温度为y，空调的设定温度为x，x和y都是时间t的函数，并且满足某个微分方程，现在我们要控制空调的制冷和加热系统，让y在更短的时间内更快的接近x或者空调最节能，首先就要保证这个控制系统稳定.特别是对于这种带时滞的系统，不稳定的情形往往是这样：假如室内温度是32度，设定温度是26度，模型不稳定的话有可能会过制冷一直到23度，然后又会加热到30度，接着又制冷到23度，再加热到30度，无限工作下去，这就是临界稳定，甚至在绝对不稳定的情况下，温度波动会离26度的平衡位置越来越远.

2.2.2二阶（平面）方程的平衡点和稳定性

二阶方程的一般形式可用两个一阶方程表示为

dx1(t)f(x1,x2)dt （2.2.3） dx(t)2g(x,x)12dt右端不显含t，这也是一个自治系统.

方程组

f(x1,x2)0 （2.2.4） g(x,x)01200)称为自治系统（2.2.4）的平衡点.记为P0(x10,x2). 的实根(x10,x2如果从所有可能的初始条件出发，自治系统（2.2.4）的解x1(t),x2(t)都满足

0limx1(t)x10, limx2(t)x2 （2.2.5） tt0)是稳定的（渐近稳定的）则称平衡点P0(x10,x2；否则，称P0是不稳

定的.

我们仍然用直接法讨论自治系统（2.2.4）的平衡点的稳定性，先做线性近似.考虑线性常系数方程

dx1(t)a11x1a12x2dt （2.2.6） dx(t)2a21x1a22x2dt系数矩阵记为

aaA1112 a21a22并假定det(A)0，则原点P0(0,0)是方程组（2.2.6）的唯一平衡点，它的稳定性由特征方程

det(IA)0

的根（A的特征根）决定，特征方程det(IA)0可以改写成下列形式：

2pq0p(a11a22) qdet(A)则特征根

1,2(pp24q)

12方程组（2.2.6）的解一般形式为c1etc2et（12）或(c1c2t)et121（12），其中c1,c2为任意实数.由平衡点稳定性的定义式（2.2.5）可知，当1,2全为负数或有负实部时，P0(0,0)是稳定的平衡点，反之，当1,2有一个为正数或有正实部时，P0(0,0)是不稳定的平衡点.

微分方程稳定性理论将平衡点分为结点、焦点、鞍点、中心等类型，完全由特征根1,2或相应的p,q取值决定，表1简明地给出了这些结果.

表1 平衡点的类型和稳定性 1,2 p,q 平衡点类型 稳定性 稳定 不稳定 不稳定 12 12 p0,q0,p2q 稳定结点 p0,q0,p2q 不稳定结点 10 12 12 q0 p0,q0,p2q 鞍点 稳定退化结点或稳定稳定 临界结点 不稳定退化结点或不p0,q0,p2q 不稳定 稳定临界结点 稳定 不稳定 不稳定 1,2i p0,q0,p2q 稳定焦点 1,2i p0,q0,p2q 不稳定焦点 1,2i p0,q0 中心 由表1可以看出，根据特征方程的系数p,q的正负很容易判断平衡点的稳定性，准则如下：

若p0,q0，则平衡点稳定； 若p0或q0，则平衡点不稳定.

以上是对线性方程组（2.2.6）的平衡点P0(0,0)稳定性的结论，根

据李雅普洛夫理论，对于一般的非线性自治系统（2.2.4），可以用近

0)的稳定性，即 似线性方程来判断其平衡点P0(x10,x2dx1(t)000000f(x,x)(xx)f(x,x)(xx)x1211x122212dt （2.2.7） dx(t)2g(x0,x0)(xx0)g(x0,x0)(xx0)x11211x21222dt系数矩阵

fx1Agx1fx2 gx2P0(x10,x20)特征方程系数为

,qdet(A) p(fx1gx2)P(x0,x0)0120)点对于方程（2.2.7）的稳定性可以由表1或特征显然，P0(x10,x20)是否自治系统方程的系数p,q的正负判定，同样方法可以判定P0(x10,x2（2.2.3）的稳定点.

2.3 微分方程模型

建立微分方程模型,一般有三种方法.

一是应用已知规律直接列方程建模 在数学、力学、物理、化学等学科中

已有许多经过实践检验的规律和定律，如牛顿运动定律、曲线的切线的性质等，这些都涉及

到某些函数的变化率.由于本身就是微分方程形式，我们就可以根据相应的规律直接列出方程，从而建立数学模型.

二是用微元法建模. 用微元法建立常微分方程模型，实际上是寻求微元之间的

关系式.在建立这些关系式时也要用到已知的规律和定理.与第一种方法不同之处在于，这里不是直接对未知函数及其导数应用规律和定理来求关系式，而是对某些微元来应用规律,从而建立相关模型.

三是用模拟近似法建模. 在社会科学、生物学、医学、经济学等学科的实践中，

常常要用模拟近似法来建立微分方程模型.在这些领域中的一些现象的规律性我们还不是很清楚，即使有所了解也通常极其复杂，因此，在实际应用中，总要经过一些简化、近似的过程，并在不同的假设下建立微分方程，从数学上求解或分析解的性质，再同实际情况作对比，观察这个模型能否模拟、近似某些实际的现象.

这三种方法中，我们在学习微分方程时做的应用题就属于第一种，这种方法相对比较简单，在这里，我们主要介绍用后两种方法来建微分方程模型.

在实际的微分方程建模过程中，往往都是上述三种方法的综合应用。不论应用哪种方法，通常都要根据实际情况，作出一定的假设与简化，并且要把模型的理论或计算结果与实际情况进行对照验证，以修改模型使之更准确地描述实际问题并进而达到预测预报的目的。

例2.3.1 体重问题

1、问题：研究人的体重随时间的变化规律. 人的体重变化的过程是一个非常复杂的过程.这里，我们进行了简化，只考虑饮食和活动这两个主要因素与体重的关系.比如，研究减肥者、运动员的体重等. 2、基本假定

（1）对于一个成年人来说体重主要由三部分组成：骨骼、水和脂肪.骨骼和水大体上可以认为是不变的，我们不妨以人体脂肪的重量作为体重的标志.而人体的脂肪是存储和提供能量的主要方式，体重的变化就是能量的摄取和消耗的过程.已知脂肪的能量转换率为100%，每千克脂肪可以转换为41868焦耳的能量.

（2）人体的体重仅仅看成是时间t的函数w(t)，而与其他因素无关，这意味着在研究体重变化的过程中，我们忽略了个体间的差异（年龄、性别、健康状况等）对体重的影响.

（3）体重的变化是一个渐变的过程，因此，可以认为w(t)随时间是连续变化的，即w(t)是连续函数且充分光滑，也就是说，我们认为能量的摄取和消耗是随时发生的.

（4）不同的活动对能量的消耗是不同的，例如：体重分别为50千克和100千克的人都跑1000米，所消耗的能量显然是不同的.可见，活动对能量的消耗也不是一个简单的问题，我们假设研究对象会为自

己制订一个合理且相对稳定的活动计划，我们可以假设在单位时间（1日）内人体活动所消耗的能量与其体重成正比.

（5）假设研究对象用于基本新陈代谢（即自动消耗）的能量是一定的.

（6）假设研究对象对自己的饮食有相对严格的控制，在本问题中，为简单计，我们可以假设人体每天摄入的能量（即食量）是一定的.

（7）根据能量的平衡原理，任何时间段内由于体重的改变所引起的人体内能量的变化应该等于这段时间内摄入的能量与消耗的能量的差.即体重的变化等于输入与输出之差，其中输入是指扣除了基本新陈代谢之后的净食量吸收；输出就是活动时的消耗.

为此，我们采集了某运动员的相关数据，食量是10467（焦/天），其中5038（焦/天）用于基本的新陈代谢（即自动消耗）.在健身训练中，他所消耗的热量大约是69（焦/公斤•天）乘以他的体重（公斤）.试研究此人的体重随时间变化的规律. 3、建模分析与量化

问题并没有直接给出有关“导数”的概念，但是，体重是时间的连续函数，就表示我们用“变化”观察来考察问题。 量化：w0: 为第一天开始时的体重，

t ：时间以天为单位

则，每天：体重量的变化＝输入－输出

输入＝总热量－基本新陈代谢热量

＝净热量吸收 ＝10467 焦－5038 焦 = 5429 焦

输出＝训练时消耗=69W 焦.

4、建立模型

WdW t0tdtlim即

(104675038)69WdW41868dt W|t0W05、模型求解：应用分离变量法，有

dWdt

129616W10,000解得

W1296(129616W0)exp(16t) 100001616这就是此人体重随时间的变化规律. 6、模型讨论

现在我们再来考虑一下：此人的体重会达到平衡吗?若能，那么，这个人的体重达到平衡时是多少公斤？事实上，从W(t)的表达式可知，当t时，W(t)81，因此，平衡时，体重为81公斤.

当然，如果我们需要知道的仅仅是这个平衡值，由稳定性理论可知，就是要求微分方程的平衡点.在平衡状态下，W(t)是不发生变化的，所以

dW(104675038)69W0，直接求得W平衡81. dt418687、模型改进

我们在这里只是简单的讨论了基本新陈代谢（即自动消耗）、饮食和活动取固定值时的规律，进一步，还可以考虑饮食量和活动量改变时的情况，还有基本新陈代谢随体重变化的情况，等等.在此，我们就不作讨论了.

例2.3.2 捕鱼业的持续收获模型

渔业资源是一种再生资源，要保持可持续发展，就要适度开发，在捕捞时，应在持续稳产的前提下再追求产量或最优经济效益.

考察一个渔场，其中的鱼量在天然环境下按一定规律增长，如果捕捞量恰好等于增长量，那么渔场鱼量将保持不变，这个捕捞量就可以持续，在这里，我们要建立在捕捞情况下渔场鱼量要遵从的数学模型（即微分方程），分析渔场鱼量稳定的条件，并且在稳定的条件下讨论如何控制捕捞强度，使持续产量或经济效益达到最大，以及研究捕捞过度的问题.

我们从产量和效益着手建立模型，鱼量的自然增长是无法控制的，分析如何控制捕捞强度实现产量最大或效益最高.由于渔场有其自身的，所以我们假设其有最大容量，并认为它与人口的增长模型相似．这样便于模型的假设和建立．

（一）产量模型 模型一

1、问题：研究在捕捞情况下渔场鱼量要遵从的规律，并讨论在可持续捕捞的条件下如何控制捕捞强度，使持续产量达到最大. 2、基本假定

记时刻t渔场中鱼量为x(t)，关于x(t)的自然增长和人工捕捞作如下假设：

（1）在无捕捞条件下x(t)的增长服从Gompertz规律，即：

x(t)f(x)rxlnN xr是固有增长率，N是环境容许的最大鱼量，用f(x)表示单位时间的

增长量.

（2）单位时间的捕捞量（即产量）与渔场鱼量x(t)成正比，比例

常数k表示单位时间的捕捞率，k可进一步分解成kqE，E称为捕捞强度，用可以控制的参数如出海渔船数量、渔船的规模、渔网的规格等来度量；q称为捕捞系数，表示单位强度下的捕捞率. 为方便起见，可以选择合适的捕捞强度单位，使q1.

3、建模分析

单位时间内渔场鱼量的变化=自然增长量-捕捞量 要使渔场鱼量将保持不变，就必须满足：

捕捞量=增长量

由假设可知，单位时间的捕捞量为

h(x)kxEx

单位时间内渔场鱼量的变化为

F(x)f(x)h(x)

4、建立模型 由上面的分析，可知

xdx =F(x) t0tdtlim即在捕捞情况下渔场鱼量满足的方程

x(t)FxfxhxrxlnNEx x5、模型求解

我们并不关心渔场鱼量随时间变化的详细过程，所以不需要解上述方程以得到x(t)的动态化过程，而只希望知道渔场的稳定鱼量与保持稳定的条件，即时间t足够长以后渔场鱼量x(t)的趋向，并由此确定最大持续产量.为此可以直接求方程的平衡点并分析其稳定性.

令

FxfxhxrxlnNEx0 x得到一个平衡点

x0NeEr

又F'(x0)r0，x0点稳定,即x0为稳定的平衡点.

由于模型的建立是在不破坏鱼群的基础上，所以易知Er，这表明，只要捕捞适度 (Er)，就可使渔场鱼量稳定在x0，从而获得持续产量h(x0)Ex0；而当捕捞过度时 (Er)，渔场鱼量将减至x1= 0，当然谈不上获得持续产量了.

进一步讨论渔场鱼量稳定在x0的前提下，如何控制捕捞强度E使持续产量最大的问题.

根据前面分析，当捕捞量等于自然增长量时，可得最大产量.要使产量最大，就要是自然增长量达到最大，由

f'(x)rlnN eNr eNr0 x得到获得最大持续产量时，稳定平衡点为

x0且单位时间的最大持续产量为

hmf(x0)不难算出保持渔场面鱼量稳定在x0的捕捞强度为

E*hmr x0由上面分析可知，此模型将捕捞强度控制在E*，或者说使捕捞率等于增长率时，可以获得最大持续产量.

模型二

1、问题：同模型一. 2、基本假定

模型的假设与模型一相似，只是增长规律不同.假设在无捕捞条件下x(t)的增长服从Logistic规律，即：

xx(t)f(x)rx1.

N3、建模分析 与模型一相同. 4、建立模型

此时，在捕捞情况下渔场鱼量满足的方程为

xx(t)F(x)rx1Ex

N5、模型求解

同样，我们并不关心渔场鱼量随时间变化的详细过程，而只关心渔场的稳定鱼量与保持稳定的条件，并由此确定最大持续产量.为此可以直接求方程的平衡点并分析其稳定性.

令

xF(x)rx1Ex0，

N得到两个平衡点：

Ex0N1,x10.

r从而

F(x0)Er,F(x1)rE，

所以若Er（即捕捞率小于固有增长率），则有F(x0)0,x0是稳

定点；F(x1)0,x1是不稳定点.

若Er（即捕捞率大于固有增长率），则结论正好相反.

由于E是捕捞率(此时kE)，r是固有增长率，上述分析表明，只要捕捞适度 (Er)，就可使渔场鱼量稳定在x0，从而获得持续产量h(x0)Ex0；而当捕捞过度时 (Er)，渔场鱼量将减至x1= 0，当然谈不上获得持续产量了.

进一步讨论渔场鱼量稳定在x0的前提下，如何控制捕捞强度E使持续产量达到最大的问题.

根据前面分析，当捕捞量等于自然增长量时，可得最大产量.要使产量最大，就要是自然增长量达到最大，由

f'(x)r2rx0 N得到获得最大持续产量时，稳定平衡点为

x0N 2rN 4且单位时间的最大持续产量为

hmf(x0)不难算出保持渔场面鱼量稳定在x0的捕捞强度为

E*hmr x02综合上述，产量模型的结论是将捕捞强度控制在E*，或者说使渔场鱼量保持在最大鱼量 N的一半时，可以获得最大持续产量.

6、模型讨论

事实上，我们讨论的模型还是非常简化的.比如，捕捞系数会随着捕捞强度的增大而减少，而我们认为是常数，这并不符合事实，特别是渔场较小的时候，更是如此.还有，通常各个国家每年都有禁渔期和捕捞期，并且对渔网的大小也规定，这是为了保护幼鱼，防止捕捞过度造成资源枯竭，等等.这些都会使捕捞系数难于估计.同时，由于渔场的边界难以控制，最大容量也就难以测得.而且在限定容量里鱼群的数量增多，则鱼的死亡率就会增大，这些都造成渔场鱼量的估计也很困难.通过对这些因素的讨论，可以得到更完善的一些模型，这也是进一步研究的方向.

（二）效益模型

1、问题：研究在捕捞情况下渔场鱼量要遵从的规律，并讨论在可持续捕捞的条件下如何控制捕捞强度，使持续经济效益达到最佳. 2、基本假定

如果经济效益用从捕捞所得的收入中扣除开支后的利润来衡量，并且简单地假设：鱼的销价单价为常数p，单位捕捞强度（如每条出海渔船）的费用为常数C. 其他假设如产量模型二的假设.

3、建模分析

单位时间的利润=收入-支出

那么单位时间的收入T和支出S分别为

Tph(x)pEx

SCE

4、建立模型 单位时间的利润为

RTSpExCE

在稳定条件下 (xx0)，

ER(E)T(E)S(E)pNE1CE

r5、模型求解 令

dR2pN1EC0 dEr容易求出使利润R达到最大的捕捞强度为：

rC ER12pN在最大利润下的渔场稳定鱼量xR及单位时间的持续产量hR为：

ENCxRN1R

r22prNC2hRERxR122 4pN,hm相比较，可以看出，在最大经将ER,xR,hR与产量模型中的E*,x0济效益原则下捕捞强度和持续产量均有所减少，而渔场稳定鱼量有所增加，并且减少或增加的比例随着捕捞成本C的增长而变大，随着销售价格p的增长而变小. 这显然是符合实际情况的.

（三）捕捞过度模型

过度开发造成资源枯竭是可再生资源面临的严重问题，人们盲目的追求自己眼前的利益，在很大程度上破坏了生态平衡，为了长久的

利益，以及生态环境的平衡，我们有必要研究一下捕鱼过度的问题，以及如何防止这种现象的发生.

1、问题：研究在捕捞情况下渔场鱼量要遵从的规律，并讨论盲目捕捞下捕鱼过度的问题，以及如何防止这种现象的发生. 2、基本假定

上面的效益模型是以计划捕捞（或称封闭式捕捞）为基础的，即渔场单独的经营者有计划地捕捞，可以追求最大利润，但是如果有众多盲目的捕捞者，即使有微薄的利润，经营者也会去捕捞，这种情况称为盲目捕捞. 在这种情况下，随着鱼的价格和捕捞成本的变动，鱼量将会迅速减少，出现捕捞过度，甚至导致鱼群种类的灭绝.同时，我们假设：鱼的销价单价为常数p，单位捕捞强度（如每条出海渔船）的费用为常数C.其他假设如产量模型二的假设.

3、建模分析

只要利润R(E)0，盲目的经营者会加大捕捞强度；一旦利润

他们当然就会减小捕捞强度.所以R(E)0是盲目捕捞下的临R(E)0，

界点.

类似效益模型，单位时间的收入T和支出S分别为

Tph(x)pEx

SCE

4、建立模型 单位时间的利润为

RTSpExCE

在稳定条件下 (xx0)，

ER(E)T(E)S(E)pNE1CE

r5、模型求解

令R(E)0 的解为ES，则

CESr1 pN当EES时，利润R(E)0，盲目的经营者会加大捕捞强度；当

EES时，利润 R(E)0，他们会减小捕捞强度.所以ES是盲目捕捞下

的临界捕捞强度.

ES也可以用图解法得到. 容易知道ES存在的必要条件是

pC N即销价大于（相对于总量而言）成本价. 由ESr1pN可知成本价

越低，售价越高，则ES越大.

在盲目捕捞下的渔场稳定鱼量为：

ECxSN1

rpxS完全由成本-价格比决定，随着价格的上升和成本的下降，xS将迅

C速减少，出现捕捞过度.

比较效益模型和捕捞过度模型，可知ES2ER. 即盲目捕捞强度比最大效益下捕捞强度大一倍.

从ES表达式还可以看出，当济学捕捞过度；当p2

2.4 建模案例及分析

我们下面以2003年高教社杯全国大学生数学建模竞赛A题SARS的传播为例，对问题2建立SARS传播模型进行了讨论.

SARS（Severe Acute Respiratory Syndrome，严重急性呼吸道综合症, 俗称：非典型肺炎）是21世纪第一个在世界范围内传播的传染病。SARS的爆发和蔓延给我国的经济发展和人民生活带来了很大影响，我们从中得到了许多重要的经验和教训，认识到定量地研究传染病的传播规律、为预测和控制传染病蔓延创造条件的重要性。请你们对SARS 的传播建立数学模型，具体要求如下：

（1）对附件1所提供的一个早期的模型，评价其合理性和实用性。 （2）建立你们自己的模型，说明为什么优于附件1中的模型；特别要说明怎样才能建立一个真正能够预测以及能为预防和控制提供可靠、足够的信息的模型，这样做的困难在哪里？对于卫生部门所采取

ccp2时，(ER)ESE*，称经NNc时，ESE*，称生态学捕捞过度. N的措施做出评论，如：提前或延后5天采取严格的隔离措施，对疫情传播所造成的影响做出估计。附件2提供的数据供参考。

（3）收集SARS对经济某个方面影响的数据，建立相应的数学模型并进行预测。附件3提供的数据供参考。

（4）给当地报刊写一篇通俗短文，说明建立传染病数学模型的重要性。

1问题重述（略） 2模型假设

1)假设所考查人群的总数恒定，忽略迁移的影响，忽略自然死亡。 2)将所考查人群分为患者、治愈者、死亡者、正常人四类。 3)假设患者治愈恢复后具有免疫能力，不考虑其再感染。

4)假设所有患者均为“他人感染型”患者，即不考虑人群个体自身发病。

5)假设各类人群在人群总体中分布均匀。 6)假设已被隔离的人群之间不会发生交叉感染。 7) 附件2提供的北京市疫情统计数据真实可信。 3符号说明

X(t)患者数

Y(t)累计病人数

R(t)累计治愈人数 D(t)累计死亡人数

T采取强制措施的时间

L1病人的死亡率 L2病人的治愈率

p采取控制措施后的隔离强度

r(t)未被隔离的病人平均每人每天感染的人数

4问题分析

该问题是一个比较典型的传染病模型问题。由于SARS的传播受社会、经济、文化、风俗习惯等因素的影响，而影响疫情发展趋势的最直接的因素是：感染者的数量、传播形式以及病毒本身的传播能力、隔离强度，入院时间等，我们在建立模型时不可能也没有必要考虑所有因素，只能抓住关键因素，进行合理的假设和建模。

在这里，我们把考查人群分为四类：正常人群、患病人群、治愈人类和死亡人群，分别用H(t)、X(t)、R(t)和D(t)表示。

在SARS爆发初期，由于整个社会对SARS病毒传播的速度和危害程度认识不够，和公众对之不予重视，没有采取任何有效的隔离控制措施。当疫情蔓延到4月20号，与社会开始采取强制措施，对SARS进行预防和控制。因此SARS的传播规律可分为“控前”和“控后”两个阶段。控前模型为近似于自然传播时的SIR模型，控后模型为介入隔离强度后的微分方程模型。

为了建立SIR和微分方程模型，在这里，我们先作一些数据上的准备。

SARS的死亡率和治愈率两个参数，一般只能通过医学界对治病机理的进一步研究加以控制，在短期内不会发生变化。根据附录2的

所给的累计病人数、累计死亡人数、累计治愈人数，我们可以对L1和

L2作最小平方误差估计。

死亡率L1累计死亡人数累计治愈人数，治愈率L2＝

累计病人数累计病人数用SAS对其作线性回归，得到

L10.0530，L20.0695

5模型的建立

基于以上的分析，我们对“控制前”和“控制后”分别进行建模。设T为实施强力控制的时间（以天为单位）。当tT时，适用于“控前模型”，tT时，适用于“控后模型”。 I 控前模型

假设某地区产生第一例SARS病人的时间为T0，在(T0,T)时段，是近乎于自由传播的时段，隔离强度为0，每个病人每天感染人数r为一常数。

我们现在来考虑在t到tt这段时间内几类人群的变化情况。并通过分析各类人群的状态转化关系，建立微分方程，得到病毒传播的动力学模型。经过分析，可知

患者数的变化＝新增病人数－(死亡人数＋治愈人数) 死亡累计人数的变化＝新增死亡人数 治愈累计人数的变化＝新增治愈人数

累计病人数＝患者数＋累计死亡人数＋累计治愈人数 于是有

X(tt)X(t)rX(t)(L1L2)X(t)tD(tt)D(t)L1X(t) tR(tt)R(t)L2X(t)tY(t)X(t)D(t)R(t)当t0时，我们得到了SARS传播的控前模型：

dX(t)dtrX(t)(L1L2)X(t)dD(t)L1X(t) dtdR(t)L2X(t)dtY(t)X(t)D(t)R(t)其中，初始值

X(0)1Y(0)1 D(0)0R(0)0II控后模型

控后隔离强度从控前的0变为p。未被隔离的病人平均每人每天感染的人数r随时间逐渐变化，它从初始的最大值ab逐渐减小至最

a、小值a。可从中科院遥感所SARS课题攻关组《SARSb的值客观存在，

传播时空模型研究简报》中查到。设每个未被隔离的病人每天感染的人数

r(t)abe(tT)

其中，用来反映r(t)的变化快慢，可以用附件2中的数据估计出它的大小。

类似于控前模型的分析，我们来考虑在t到tt时段内各类人群

的变化情况。类似分析，可得

X(tt)X(t)(1p)r(t)X(t)(L1L2)X(t)tD(tt)D(t)L1X(t) tR(tt)R(t)L2X(t)tY(t)X(t)D(t)R(t)当t0时，我们得到了SARS传播的控后模型：

dX(t)dt(1p)r(t)X(t)(L1L2)X(t)dD(t)L1X(t),tT dtdR(t)L2X(t)dtY(t)X(t)D(t)R(t)其中，

r(t)abet a0.245b0.6初始值X(T)取控前模型的最后一个值。 6模型的求解 I控前模型的求解

对于患者数X(t)，我们可以根据SARS传播的控前方程，求得它的解析解为

X(t)X(0)e(rL1L2)t,tT

其中，

r0.55L0.0531 L20.0695X(0)1再将X(t)分别代入SARS传播的控后方程，就可以给出D(t)、R(t)以及Y(t)的数值解。

II控后模型的求解

同理，我们求得患者数得解析解

X(t)X(T)e[(1p)aL1L2](tT)(1P)b(1e(tT)),tT

其中，

a0.245b0.6 T47我们已经分析过，a、b为一客观参数，可以从中科院遥感所SARS课题攻关组《SARS传播时空模型研究简报》中查到。由于3月5日第一例SARS进入北京，是我们记时的起点；4月20日即为T47的情况。

p和为待估计的参数，现在来估计p和。

根据附件2中的数据，将各时刻累计病人数减去累计治愈人数再减去死亡人数，可得到患者数，估计p和的值。估计时我们按均方最小误差原则，用SAS软件计算出其估计值分别为

P65%，0.02.

我们根据以上求出的解，作出了患者数、累计死亡人数、累计治愈人数、累计病人数的曲线图，与实际公布数据进行对照（略）。可以看出，方程的解与实际数据吻合的很好，说明我们的参数和模型都

是正确可靠的。

7模型结果分析与改进方向

(1) 灵敏度分析

根据我们所建的模型，卫生部门通常可以采取两种方案对疫情进行有效控制。一是改变控制时间点T；二是改变控制强度p。现在我们分别考察它们对模型的影响。

表2.4.1隔离强度p对的模型影响 隔离强度p 55％ 65％ 75％ 累计病人数 6996 2827 1339 由表2.4.1可以看出：隔离强度p，对疫情的传播具有极大的敏感度和相关性。

表2.4.2控制时间T对的模型影响 控制时间T 累计病人数 延后5天 5382 延后4天 4729 延后2天 3733 4月20日 2879 提前2天 27 提前4天 1576 提前5天 1621 由表2.4.2可以看出：控制时间的提前或延后，对累计病人影响显著。说明控制时间T，对疫情的传播具有极大的敏感度和相关性。

(2) 收敛性讨论

针对该模型，我们要判别控后模型方程组解的收敛性，X(t)的取值至关重要，D(t)、R(t)以及Y(t)的收敛性都直接依赖于X(t)是否收敛到0.

将控后模型中X(t)的解析解取极限得：

limX(t)X(T)limett[(1p)aL1L2](tT)(1P)b.

该式为t的指数函数，其收敛性取决于自变量的系数。

当(1p)aL1L20时，limX(t)0，模型收敛，疫情能够得到控

t制;

当(1p)aL1L20，limX(t)0，模型发散，疫情难以控制。

t分析发现，模型收敛得条件为

p1L1L2 a其中，

a0.245, L1 0.053,L20.0695，

所以，要使疫情得到控制，必须使隔离强度p50%。 (3) 计算机模拟检验

为了检验模型求解结果的正确性，我们进行了仿真模拟。（此处从略）从模拟结果来看：计算机模拟结果与模型计算结果有着良好的一致性。本模型是可以信赖的SARS传播模型。

(4) 对其他地区的疫情预测（略） (5)模型的评价（略）

参考书目：

[1] 朱建青，张国梁，数学建模方法，郑州大学出版社，2003 [2] 丁同仁，李承治，常微分方程教程，高等教育出版社，2004 [3] 于义良，数学建模，中国人民大学出版社，2004 [4] 参考了http://www.mcm.edu.cn/的一些资料