实验四微分方程数值解

实验四 微分方程数值解及应用

一 引例

在18世纪末，英国人马尔萨斯在他出版的一本专著中，对人口数量的增长趋势进行了模拟，提出了人口的指数增长模式，导致最后会出现人口数量超过地球的承受容量，即人口爆炸问题。虽然马尔萨斯人口模拟忽略了一些人口增长的重要因素，但是他为以后人口模型的改进提供了基础，下面将介绍著名的马尔萨斯人口模型及改进。

马尔萨斯认为：单位时间内，人口的出生数量和死亡数量与人口总数成正比例，即人口出生率和死亡率都是常数，因此人口的净增长率为常数。

设时刻t 的人口数量为()p t ，人口出生率为b

，死亡率为d ，则有：

()()()()p t t p t bp t t dp t t +∆=+∆-∆，

从而

()

()()()dp t bp t dp t kp t dt

=-=， 其中：k

b d =- 称为净增长率（常数）。因此马尔萨斯人口模型如下：

()

()dp t kp t dt

=，()00p t p =。

该微积分方程初值问题的解析解：

()0()0k t t p t p e -=。

用此模型预测某国1850年人口数量的误差为1% ，1900年人口数量的误差为31%，2000年人口数量的误差达567%，2050年将达到80多个亿。 该模型对于长期预测不合理。

针对马尔萨斯人口模型的不足，1837年荷兰生物学家Verhulst 提出了如下改进： 由于资源的限制，人口存在最大值（极限）M 。因此，人口增长率不应该是常数，假设增长率k

是随着人口数量接近M 而线性递增：

()k r M p =-，

从而得到改进后的人口模型为

()

()dp t kp r M p p dt

==- ，()00p t p =， 称为逻辑斯谛增长模型。

刻画世界千变万化的规律，微分方程是最有力的工具。

二、微分方程的数值解法

讨论一阶微分方程的初值问题

00(,),()dy

f x y y x y dx

==， （1） 或一阶微分方程组的初值问题 00()(,),()x t f t x x t x ==

其中

()x t 是向量函数，t 为自变量。

微分方程的数值解法的基本原理： 引入自变量的取值点列

{}n x ，定义1n n n h x x -=-，称n h 为步长，常用等间距

的步长 （n h 与n

无关，记为h ），假设精确解为

{}()n y x 。为了寻求(n y x ）的近似

值n y ，根据一定的原理，结合当前得到的近似解，近似地表示该点或前一点的导数值，由此推出计算n y 的迭代公式。

1． 欧拉方法

欧拉方法是一种简单的求解初值问题的数值逼近方法，其基本思路为：

对方程（1），在小区间1[,]n n x x +上，用差商1()n n y x y x h

+-（） 代替导数'

y ，

1）用左端点n x 替换函数(,)f x y 中的x ，得到方程（1）的近似表达式

1()()(,())n n n n y x y x hf x y x +≈+

设(n n y x y ≈）

，则1()n y x +的近似值为 1

(,),0,1,

n n n n y y hf x y n +=+= （2）

其中00(,)x y 为初始点，（2）称为显示欧拉公式，也称向前欧拉公式。 向前欧拉法简单易于计算，但精度却不高，收敛速度慢。 2）用右端点1n x +替换函数(,)f x y 中的x ，

1

11(,),0,1,

n n n n y y hf x y n +++=+= （3）

其中00(,)x y 为初始点，（3）称为隐示欧拉公式，也称向后欧拉公式。 这是一个非线性方程，无法直接计算1n y +。在计算精度、收敛速度方面与向前

欧拉算法相同。

3）将（2）式与（3）式加以平均，得到

1111

[(,),(,)],0,1,2

n n n n n n y y h f x y f x y n +++=+

= （4）

称为梯形公式。与前两个方法相比，该方法的计算精度高、收敛速度快。但迭代计算与(3) 式一样繁！

4）改进的欧拉公式：先由公式（2）算出1n y +的预测值1n y

+，然后带入梯形公式（4）的右端，作为校正，即 1

(,)n n n n y y hf x y +=+

1111

[(,),(,)],0,1,

2

n n n n n n y y h f x y f x y n +++=+

= （5）

称为改进的欧拉公式。还可以写为

1121211

(),2(,),

(6)

(,)

n n n n n n h y y k k k f x y k f x y hk ++⎧

=++⎪⎪

=⎨⎪=+⎪⎩

人们常用的是向前欧拉公式和改进的欧拉公式。且欧拉法可以推广到求解微分方程组的情形。

练习1：用上述方法求解微分方程初值问题

'1y y x =-++

要求：取步长0.10.001h =和 ，分别用三种数值解法求解，并结合其精确解，对求解误差进行分析比较。

首先用解析法得到其精确解x

y x e -=+。

其次用数值解法(0.1)h =取 向前欧拉算法 迭代公式为

100.90.10.1,0,1,,1n n n y y x n y +=++==

向后欧拉算法 迭代公式为

111100.1(1),0.10.11)/1.1,0,1,

,1n n n n n n n y y y x y y x n y ++++=+-++=++==将它变形为（