2.7 常微分方程（组）的数值解

非线性常微分方程通常无法求解析解，即使是线性常微分方程，如果不能获得变系数的积分，也无法给出解的表达式。因此在实际应用中，研究人员经常需要计算常微分方程的数值解。下面介绍几种求常微分方程（组）数值解的方法。

求常微分方程（组）的数值解，就是选取自变量的若干个离散节点，求常微分方程（组）在这些节点处的近似值。函数的导数是以差商的极限的形式定义的，在求常微分方程（组）的数值解时，一个很容易想到的思路就是反其道而行之，将以差商作为导数的近似值。

首先介绍Euler法。对于带有初值条件的常微分方程

y′=f（t, y）, y（t0）=y0

首先要选取适当的节点t0<t1<t2<…<tn<tf，比如采用固定步长的方式，令tk+1-tk≡h。然后在各个节点处用差商代替导数：

由此可得y（tk+1）≈y（tk）+h·f（tk, y（tk））。根据这个递推公式，我们可以从y0推出所有的y（tk），这就是Euler法。

在使用Euler法求常微分方程的数值解时，也可以有多种不同的处理方法。比如可以采用向后Euler法：

y（tk+1）≈y（tk）+h·f（tk+1, y（tk+1））

以及改进Euler法：

在实际应用中，人们发现Euler法的精度无法满足要求，而综合考虑多个节点处的斜率，可以有效提高精度。龙格-库塔法就是一种考虑多个节点的方法，适用面比较广，精度也比较高。具体思路如下：

其中，K1=f（tk, y（tk））, Ki=, h为步长。

有多种算法可以用来求解常微分方程的数值解，因此MATLAB中可供使用的函数也有多个。常用的函数如表2-2所示。

需要注意的是，上述函数仅适用于非刚性（Nonstiff）方程（组）。所谓刚性方程（组），就是其数值解只有在步长很小时才会稳定，步长较大时就会很不稳定。在具体应用中，如果使用常用函数长时间没有结果，则可以考虑换用表2-3所示函数。

求常微分方程（组）数值解的命令格式为：

[t, y]=solver（odefun, tspan, y0, options）

其中，solver选择ode45等函数名；odefun为根据待解方程或方程组编写的M文件名；tspan为自变量的区间 [t0, tf]，即数值解的计算范围，t0表示初始点；y0表示初始值；options用于设定误差限制，其命令格式为：

options=odeset（reltol, rt, abstol, at）

其中，rt输入相对误差，at输入绝对误差。

例2-28 假设任意时刻，人口的增速与人口总数成正比。据此建立常微分方程模型：x′（t）=kx（t）, x（0）=x0。利用分离变量法，可以迅速求出该初值问题的解为x（t）=x0·ekt。要求其数值解，需要先设定初值问题中的各项参数。比如设k=0.05, x0=1000，建立描述方程的M文件：

    function dx=human（t, x）
    dx=0.05∗x

并保存为human. m．接着使用命令：

    >> [ t, x]=ode45（ human , [0, 100] , 1000）；    %计算 [0, 100] 范围内的数值解
    >> plot（t, x）
    >> hold on
    >> x=1000∗exp（0.05∗t）；
    >> plot（t, x）

上述命令输入后，可以得到数值解各点连线与真实解各点连线的对比图，如图2-17所示。

从图2-17可以看到，两条曲线重合度非常高。

例2-29 假设甲、乙两个种群竞争同一种资源。t时刻甲、乙各自的数量分别是x（t）, y（t），则其变化规律可以用如下常微分方程组来描述：

其中，a、b分别表示两个种群的自然增长率；N、M分别表示两物种单独生存时的数量上限；p、q分别表示乙对甲和甲对乙的竞争强度或影响能力。

要求数值解，首先编写描述方程组的M文件：

    function dx=compete（t, x）
    dx=zeros（2, 1）；
    dx（1）=0.01∗x（1）∗（1-x（1）/50000-0.1∗x（2）/60000）；
    dx（2）=0.02∗x（2）∗（1-0.2∗x（1）/50000-x（2）/60000）；

并将文件命名为compete. m．其中，x（1）、x（2）表示方程组中的x、y。在输入方程组时，设定的乙的增长率和上限都相对占优。接着输入命令：

    >> [ t, x]=ode45（compete , [0, 500] , [10, 10]）；
    >> plot（t, x（ :, 1）, t, x（ :, 2））  %在一张图上绘制t-x和t-y曲线
    >> plot（x（:, 1）, x（:, 2））          %绘制x-y曲线图

对比图2-18可以看出，在乙增长率和上限都占优的情况下，甲的数量增长会受到压制。

对于来自现实问题的常微分方程（组）模型而言，求解往往不是最重要的，更重要的是模型中函数关系变化、参数变化等会对解的性态产生怎样的影响。因此，在求解数值解之前，最好能利用常微分方程定性理论对参数的影响做一界定，这样才能更加全面地了解事物内在的运行规律。