04微分方程new

1、第四章 微分方程建模方法4.1 导言在许多实际问题中，例如物理中的速率问题，人口的增长问题，放射性衰变问题，经济学中的边际问题等，常常涉及到两个变量之间的变化规律。微分方程是研究上述问题的一种机理分析方法，它在科技、工程、生态、环境、人口以及经济管理等领域中有着十分广泛的应用。在应用微分方程解决实际问题时，必须经过两个阶段。一是微分方程的建立，建立一个微分方程的实质就是构建函数、自变量以及函数对自变量的导数之间的一种平衡关系。而正确地构建这种平衡关系，需要对实际问题的深入浅出的刻画，根据物理的和非物理的原理、定律或定理，作出合理的假设和简化并将它升华成数学问题。另一个是方程的求解和结果分析。对

2、一些常系数的或特殊函数形式的微分方程，往往能得到解析解，这对实际问题的分析和应用都是有利的,但是大多数变系数的、非线性函数形式的微分方程都是求不出解析解的,此时就需要应用求解微分方程的另一个重要方法数值解法。本章简要介绍有关微分方程模型的概念，微分方程的数值解法和图解法，主要介绍若干建模实例，通过它们展示微分方程模型的建模步骤及解决实际问题的全过程。4. 微分方程模型及其解法4.2.1微分方程模型由物理意义可知，速率就是变化率，而变化率即数学中的导数。如果一个实际问题中涉及到一个变量的变化速率，那么就有可能用微分方程模型来描述该问题。例如某国人口的增长率与其人口的总数成正比，现需要预测未来某年

3、该国的人口总数。这就是一个涉及到变量的变化速率（人口的增长率）的实际问题。要解决该问题，就必须将它转化为数学问题。如果设表示时间，表示时刻的人口总数，假设是的可导函数，则人口增长率就是，根据假设可建立微分方程一般地，包含变量、函数及其导数的等式称为微分方程。微分方程的一般形式是或在上述问题中，如果已知时刻该国的人口总数为，则这个人口预测问题就可完整地表示为 (4.1)关系式(4.1)就是上述人口预测问题的数学模型.可以求得（.1）的解析解。下面我们考察葫芦形状容器壁上的刻度问题.对于圆柱形状容器,由于体积与相对于容器底部的任意高度的函数关系明确:其中直径是常数,因此可以非常方便地在容器壁上标示

4、容积刻度.但是对于非简单规则几何形状的容器,在其壁上标示容积刻度则比较困难.现在假设有一葫芦形状的容器,人们通过测量得到高度与直径的离散数据如表4-1.表4-100.20.40.60.81.01.20.640.851.000.580.200.120.12 根据表4-1中的数据,可以大致画出该容器的状,如图4.1所示。0图4.1 葫芦形状容器 如何标出任意高度处容器体积的刻度？由微元法分析知 其中是任意高度，直径是高度的函数，记为.因此 （.）只要求解上述微分方程，就可以求出体积与高度之间的函数关系，从而标出任意高度处容器体积。（4.2）就是上述物理问题的数学模型。但是，由于无解析表达式，故无法

5、求得（.2）的解析解。这就需要利用微分方程的数值解法。下节我们主要简述关于微分方程求解的系列方法。4.2.2微分方程解法 在高等数学中，介绍了微分方程的基本知识以及一些比较特殊类型的微分方程的解析解法，但是大量的复杂微分方程，无法或难于求出其解析解。因此，在实际问题中经常需要用数值解法求微分方程和微分方程组。为了数学建模的需要，本节我们主要补充微分方程的数值解法。设有微分方程问题 （.）数值解法的基本方法是:根据需要给定自变量点列：定义步长为.常给定等间距自变量点列,此时步长为常数。数值解法一般只能求得微分方程的近似解。设微分方程(4.3)在点列上的精确解是,而其近似解为。根据一定的原理,用差

6、商代替导数,结合初始条件,推出计算的迭代公式。1 欧拉方法当较小时，在小区间上用差商代替导数，即用近似代替方程（4.3）左端导数，而方程右端函数中取,由方程(4.3)得近似表达式 (4.4)在(4.4)中近似地取,则由(4.4)得 (4.5)其中是已知的初始点。由(4.5)和已知的初始点,我们就可以一步一步地推算出:因此,我们称公式（4.5）为显式欧拉公式。如果在上述过程中，在中取，那末可得另一个欧拉公式 (4.6)大多数情况下是非线性函数,从而公式（4.6）是非线性方程,无法直接计算,因此,我们称公式（4.6）为隐式欧拉公式。一般而言,欧拉方法计算精度低,收敛速度慢。如果把显式欧拉公式（4.

7、5）和隐式欧拉公式（4.6）加以平均，得到 （4.7）我们称公式（4.7）为梯形公式。可以证明，梯形公式比欧拉公式精度高,收敛速度快。例4.1 设微分方程问题分别用显式欧拉法、隐式欧拉法、梯形法求解，取步长和0.001，并利用其解析解（精确解），对数值解误差进行分析比较。解 这是一个一阶线性非齐次微分方程，用解析解法得到的精确解是下面用数值解法：取h=0.1时，1） 显示欧拉法的迭代公式为： ， ，.2） 隐式欧拉法的迭代公式为： ， ， 将它变形为： ， 3） 梯形法的迭代公式为：， 将它变形为：， 计算结果列入表4-2。表4-2 解析解及时三种数值法计算结果解析解显式欧拉解隐式欧拉解梯形解

8、011110.11.004811.00911.00480.21.01871.01001.02641.01860.31.04081.02901.05131.04060.41.07031.05611.08301.07010.51.10651.09051.12091.10630.61.14881.13141.16451.14850.71.19661.17831.21321.19630.81.24931.23051.26651.24900.91.30661.28741.32411.306311.36791.34871.38551.3673取时，容易得到与上述类似的公式，其计算结果列入表4-3。表4-3

9、 解析解及时三种数值法计算结果解析解显式欧拉解隐式欧拉解梯形解011110.11.00481.00481.00491.00480.21.01871.01861.01881.01870.31.04081.04071.04091.04080.41.07031.07021.07051.07030.51.10651.10641.10671.10650.61.14881.14861.14901.14880.71.19661.19641.19681.19660.81.24931.24911.24951.24930.91.30661.30641.30681.306611.36791.36771.36811.

10、3679由表4-2和表4-3首先可见，取步长时，计算结果具有两位有效数字，而 时，计算结果具有四位有效数字。（设为一准确值，为的一个近似值，则称为的误差。若的误差不超过的某一位上的半个单位，而且从该位到的左边第一位非零数字共有位，则称具有位有效数字 ）一般说来，在迭代中，步长越小，计算结果越精确（当步长在数值方法的绝对稳定区域内时）。其次可以发现，计算点离开初始点越远，的误差越大。最后还可知，梯形法显然优于显式、隐式欧拉法，例如取步长时，梯形法的计算结果基本具有四位有效数字（仅在处，只有三位有效数字）。从上例看好象用梯形法计算非常容易，其实不然，梯形公式是一种隐式公式，一般情况下计算比较困难。

11、于是产生了一种改进的方法：第一步，由显式欧拉公式（4.5）算出的预测值，第二步，将代入梯形公式（4.7），进行校正，即 （4.8）我们称公式（4.8）为改进的欧拉公式。在程序编写中，它常常写为 （4.9）改进的欧拉公式是显式公式，易于计算，精度较高，收敛速度快，是人们最常用的方法之一。欧拉法和改进的欧拉法还可以推广到微分方程组的情形。例4.2 用改进欧拉公式求例4.1中的初值问题。解 由公式（4.8）知由于故化简可得计算结果列于表4-4。表4-400.10.20.30.40.511.005001.0190251.0412181.0708021.107076对比例4.1与例4.2可见改进欧拉公式

12、的精确性高于欧拉公式，与梯形公式相当。2.龙格-库塔方法利用泰勒公式可以得到龙格-库塔方法（简称方法）是由上式产生迭代公式进行计算。若则称以上迭代公式为阶公式，越大，截断误差越小，精度越高。要得到一个阶公式，关键是选择函数，使之满足截断误差为的要求。 经常使用的公式有1） 2阶公式 中点公式 （4.10）2）3阶公式 （4.11）3）4阶公式 经典龙格-库塔公式 （4.12）关于微分方程组也有类似的数值解法。在MATLAB软件中含有数值求解的系统函数，其实现原理就是龙格-库塔方法。例4.3 用经典龙格-库塔法求初值问题在处的解。 解 ，由公式（4.12）得 ， ，本题的精确解是，现将计算结果列

13、表如下：表4-500.20.40.60.81.011.183231.341671.483281.612511.7321411.183221.341641.483241.612451.73205 下面我们简要介绍图解法。 高等数学中介绍的解析解法，虽然可以得到精确解，但是应用范围狭窄。上面介绍的数值解法适应面广，虽然得到的是离散点上解的近似值，但可以控制其误差，因此得到了广泛地应用。下面简要介绍微分方程的图解法。该法的特点是可以将微分方程解的全局信息直观地、形象地展现出来。如果一个一阶微分方程可以写成如下形式其中是已知的二元函数，那么我们就能够确定积分曲线（即解函数）在任意点处的斜率。从图象上看

14、，给出平面上一个点，我们就可以画出一条通过该点的短直线，这条短直线的斜率等于积分曲线在该点处准确的斜率。适当设定在坐标平面上一定数量的点处画这种短直线以及这种短直线的长度，这种短直线布满整个坐标面所形成的图形称为斜率场。另一方面如果我们求出了微分方程的通解那么我们就得到一蔟积分曲线，不同的初始值对应于不同的解曲线。在这些曲线蔟中任意取定一条，进行如下处理：1） 曲线分段: 将曲线分割成若干弧线段。2） 以直代曲： 将每一弧线段用其左端点处的短切线段代替。3） 相邻不连： 短切线段足够地短，相邻短切线不相连。经过上述处理，由微分方程的解可以得到斜率场，反之，由微分方程直接得到斜率场后（不需求微分

15、方程的解析解），从斜率场可以隐约看出积分曲线的形状。可以将斜率场看作一簇路标，它指示短切线左端点处要走的方向。要得到微分方程初值问题的解曲线，就可以从坐标平面上的初始点出发沿该点斜率场所指示的方向前进，按斜率场所指示的走向一段一段地走下去，就会走出一条“路”。这条“路”就是微分方程的积分曲线在坐标平面上的图形。大家知道图形也是函数的一种表现形式，因此，能画出上述“路”就相当于确定了解函数。 4. 新产品的推销模型一种新产品上市后，经销者自然十分关心它的销售情况，尤其是销售一段时间后，经销者往往需要根据已有销售情况，预测该产品在本地区的总销售量，从而恰当地组织货源。下面我们根据两种不同的假设建立

16、两种产品销售速度模型。第一种模型 假设产品是以自然推销的方式卖出，即被卖出的产品起着实物性广告宣传作用，自然地吸引着暂未购买的潜在消费者。设在时刻已售出的产品总数为，并假设售出的每一产品在单位时间内平均可以吸引个顾客，使其购买个该产品，则满足微分方程 （4.13）设初始条件为 （4.14）上述微分方程的解是 （4.15）此即为著名的模型。关于这个模型，我们作如下几点说明：（1）经过与实际销售情况对比，发现（4.15）式的计算结果与实际销售量在初始阶段的增长情况比较符合。（2）在产品即将卖出、尚未卖出之时刻为，显然，这时由（4.15）式知，这一结果自然与事实不符。产生这一错误结果的原因在于我们假

17、定产品是自然推销的，然而，自然推销方式无法自动启动，因为在最初产品还未卖出之时，在消费者中还没有作为实物性广告的产品，从而无法吸引消费者来购买。为此，我们这样来解决上述矛盾。将第一个单位时间内看作是以非自然推销方式进行销售（商家做广告），而将第一个单位时间之后看作以自然推销方式进行销售，因此，我们将时刻的销售情况作为初始条，即 （4.16）微分方程初值问题（4.13）、（4.16）的解为 （4.17）(4) 由（4.15）或（4.17）得，一般对耐用消费品而言，这显然与事实不符。事实上，耐用消费品，比如冰箱、电视机、电脑和汽车，需求量总是有限的，因此往往是有上界的。从上面的分析可见，模型不宜用

18、于销售量的中、长期预测。下面我们来介绍一种有界模型。第二种模型 假设：1）需求量的上界为，2）经销商可以通过自然推销方式和其它方式推销产品，3）一个消费者仅购买一件该种产品。根据这三个假设，在某时刻时，产品销量的增长既与到时刻为止的已经购买该种产品消费者数目成正比，也与尚未购买该种产品的潜在消费者数目成正比，即设比例系数为，则 （4.18）注意到初始条件用分离变量法可求得上述微分方程的解 （4.19）0上式为模型。当时如果则由（4.19）式得到另外，在（4.19）中令（见图4.2）。可见上述模型是有界的。图4.2由（4.19）式知， 即是关于时刻的单调函数，这与实际情况吻合，卖出的产品累积量自

19、然越卖越多。此外，对（4.19）式两端关于求导，得故令，得到。如图4.2，当时，由函数图象为上凹弧；同理，当时，单调减小，函数图形为上凸弧。这说明，在销售量小于最大需求量的一半时，销售速度是不断增加的；销售量达到最大需求量的一半时，该产品正处于最畅销时刻，此后销售速度开始下降。实际情况表明，产品销售情形与模型十分相似，特别在销售后期更加吻合。如果时，可用初始条件此时可得解 （4.20）问题一 某公司生产的冰箱在某个城镇10年的销售记录为（单位：万台）年序号12345678910销量量1.01.02.25.57.58.59.08.57.85.7试预测该种冰箱在该城镇的总销售量。解 设在时刻已售出

20、的冰箱总数为，满足模型（4.20），即此处需要求出最大需求量。由的定义知， 一般地，可得，用差商近似代替导数，注意到，可将微分方程化为处（）的差分方程即写成矩阵形式为记,上式可写为这是一个超定线性代数方程组,在上式两端左乘得如果的逆存在，用左乘上式两端，就得到将具体数据代入得 ， ， ， ， 故。因此，我们预测该种冰箱在该城镇的总销售量为69.29万台。另外，我们还可利用模型预测以后某年的该种冰箱在该城镇的销售量。4.4抵押贷款购房问题问题 随着改革开放的不断深入，我国开始实行住房货币化分配制度，一般人家都无力用现金买下自己满意的住房，从而面临贷款购房问题。下面是1991年1月1日某大城市晚报

21、上刊登的一则商品房广告。 名流 名 流 花 园 之 高 尚 住 宅 公 寓 花园 首 期 隆 重 发 售 用薪金，买高品质住房 对于大多数工薪阶层的人士来说，想买房，简直是天方夜谈。现在有这样一栋：自备款只需七万元人民币，其余由银行贷款，分五年还清。相当于每月只需付1200元人民币。那么，这对于您还有什么问题呢？ 任何人看了这则广告都会产生许多疑惑，比如广告中没有说明住房面积、房型、设施以及所处地段等等，下面我们的数学建模仅针对人们更关心问题：如果一次付款买这栋房要多少钱呢？银行贷款的利息是多少呢？为什么每月要付1200元呢？是如何算出来的？因为人们都知道，如果知道了一次付款买房的价格，自己暂

22、时只能支付一部分，其余的房款就不得不找银行贷款，只要知道了利息，就应该可以算出五年内每月应付多少钱才能按时还清贷款了，从而就可以对是否要去买该广告中所说的房子作出决策。现在我们分四个部分讨论。一 明确变量、参数，显然下面的量是需要考虑的： 表示需要贷款量； 表示月利息（贷款通常按复利计）； 表示借期，按月数计；二 建立变量之间的数学关系。若用表示第个月时尚欠的款数，则一个月后，本钱加上利息为，假如每月还款元，则第个月时尚欠的款数为因此，我们得到如下数学模型 （4.21）三 数学模型（4.21）的求解。由 应用数学归纳法，易知 ， （4.22） 四 根据广告中的信息，我们可知，年个月，每月还款元

23、，但是需要贷款数（即一次性付款购买价减去首次付款70000万元后剩下的款数）和银行贷款利率，广告未予说明，这造成了我们决策的困难。然而，由（4.22）可知60个月还清，即，从而得从中将解出的得 （4.23）（4.23）表示给定时与之间的关系式，如果我们已知银行的贷款利率，就可以算出。例如，若，则由（4.23）可算得元。因此，我们就计算出了该房地产公司该种房的一次性付款的房价应该是70000+53946=123946元。如果银行对房地产公司和个人所发放的贷款，其贷款利率相同，都是，如果该房地产公司说一次性付款的房价小于123946元的话，就说明找房地产公司分期付款买房不划算，你就应该自己银行贷款

24、。利用（4.22）式，以下我们进一步讨论两个问题。例1 某单位一对年青夫妇为买房要用银行贷款60000元，月利率0.01，贷款期25年=300月，这对夫妇希望知道每月要还多少钱，25年就可还清。假设这对夫妇每月可有节余900元，是否可以去买房呢？解： 根据上面的讨论，现在的问题就是要求使的，由（4.22）式知将代入，算得元，这说明这对夫妇有能力贷款买房。例2 接着例1的背景讨论，这对夫妇又看到某借贷公司的一则广告：“若借款60000元，22年还清，只要（i）每半个月还316元；(ii)由于增加了文书工作量，你要预付三个月的款，即6323=1896元。”这对夫妇算了一下：每半个月跑一趟去交款31

25、6元，一个月正好交款632元，虽然预付三个月的款1896元，这使人不快，但提前三年还清可省下632123=22752元，这是1896元的12倍！面对这一计算结果，这对夫妇又产生了一个疑问：这家借贷公司是慈善机构呢还是仍然要赚我们的钱呢？现在我们来给他们一个满意的回答。解 我们先就（i）,(ii)两条作一个初步分析，先孤立地分析（i）和(ii)看看能否缩短还款期。分析（i），这时不变，月利率变为半月利率，可粗略地认为半月利率正好是月利率的一半，即，于是由（4.22）式可求出使的还款期，即由求得 （4.24）将，代入（4.24）式得（半月）年，即只能提前大约一个月还清贷款。由此可见，该借贷公司如果

26、在广告中只有条件（i）的话，它一定会亏本，而这是不可能的，除非它是一个慈善机构。分析（ii）,此时元，即你只借了58104元而不是60000元，按问题1中银行贷款的条件算一下，即令（月息），（每月还款），求使的，来看看能否提前还清。将，代入（4.24）得（月）（年），这表明实际上不需通过借贷公司，直接向银行借贷，就可以提前接近四年还清贷款。因此，该借贷公司只要去同样的银行借款，即使半个月收来的316元不动再过半个月合在一起去交给银行，它还可以坐收第22年的用户交款（22-21.09）（年）12（月）632（元）=6901.44（元），更何况它可以利用每上半个月收到的还款去做短期（半个月内）投资

27、赚取额外的钱。根据以上分析，我们现在可以告诉这对夫妇，这家借贷公司不是慈善机构，它还是要赚你们的钱的。所以你们应该直接找银行贷款。4.5 利用微迹元素反演膏盐的沉积环境在天然地质体系中，微迹元素在平衡共存的两相间的分配服从能斯特(Nernst)分配定律,即当两种相处于化学平衡时，某种元素i在两相间的浓度比为常数Kd。在一定的浓度范围内的Kd与i的浓度无关，当两种相均为聚集相时，Kd受压力影响较小，而与温度的关系较为显著。本文试图根据微迹元素在溶液中由于蒸发而结晶沉淀过程中服从能斯特分配定律，从而建立沉积盆地演化环境的数学模型。一、数学模型的建立蒸发岩层中盐岩（NaCl）总是含有微量的Br。这是

28、由于Br取代盐岩晶格中的Cl。Braitsch和Herrmann（1963）测定了Br在盐岩与海水之间的分配系数Kd，发现在25盐岩沉淀开始时，Kd0.15。这意味着在盐岩沉积过程海水中的Br含量就迅速增长，随着溴化物浓度增高，Kd值可近似地保持不变。 因为在溶液中，只有当溶质的含量超过它的溶解度时，才发生结晶沉淀。因此可以假设在沉淀过程中，海水中NaCl溶液总是饱和的，沉淀的原因是由于海水的蒸发而引起的。假设在沉淀发生过程中温度不发生变化，则1m3饱和的海水中含NaCl为（0.36kg），重量为11.36kg。、都可以看成是固定不变的。另一方面，假设海水的蒸发量为 (m/a)，海水的体积为

29、(m3)，海水中Br的浓度为 (ppm)，、是沉积开始时海水的体积和海水中Br的浓度，假设注入盆地的注水量为 (m/a)，其中1m3的注入水中含NaCl为 (kg)，注入水中Br的浓度为，为盆地的面积。 从时刻到时刻，由于海水的蒸发而晶出NaCl晶体，从而使海水中Br的含量减少了 (4.29)同时，由于注入了非饱和的含有NaCl的补充水，上式可以改写为 (4.30)式（4.30）满足的条件是。即海水的蒸发量要大于使注入水成为饱和溶液的量。 另一方面，从时刻到时刻，海水中Br的含量减少 (4.31)因此有 (4.32)令，有 (4.33)式（4.33）是的一阶微分方程令 (4.34) (4.35)则 (4.36) 如果盆地每年蒸发率与注水率均为常数，计为与，那么可以计算出盆地海水的体积为 ， (4.37)则（4.36）式可以简化为： 令 ， 当时 (4.38.1) 当时 (4.38.2)由此可以计算出NaCl晶体中Br的浓度 (4.39) 另一方面，由于海水的蒸发，NaCl结晶沉积，其沉积厚度可表示为 (4.40) (4.41)其中为晶体NaCl的密度(kg/m)。 从(4.38)、(4.39)、(4.41)可以消去时间t得出沉积盆地中晶体NaCl中Br的浓度（或含量）与其沉积厚度的关系式。 二、实例 假设某膏盐沉积盆地中晶体NaC