2.2 数值积分算法的基本分析

2.2.1 单步法和多步法

上节中介绍的几种数值积分算法都有一个共同特点：在计算yk+1时只用到了yk，而不直接用yk-1，yk-2，… 项的值，即在本次计算中仅仅用到前一次的计算结果，而不需要利用更前面各步的结果。这类方法称为单步法。单步法有如下优点：

● 需要存储的数据量少，占用的存储空间少；

● 只需要知道初值，就可以启动递推公式进行运算，即可以自启动；

● 容易实现变步长运算。

单步法的缺点是随着精度的提高，计算工作量会增加。

与单步法相对应的还有一类数值积分算法，称为多步法。在它的计算公式中，本次计算不仅要用到前一次的计算结果，还要用到更前面的若干次结果。例如，四阶阿达姆斯（Adams）显式法（简记为AB4法）

就是多步法，式中

fk-i=f（tk-ih，yk-i），i=0，1，2，3

多步法与单步法相比，要达到相同的精度，计算工作量要少得多，这一点从式（2.30）与式（2.21）的比较中可以明显看出。在计算yk+1时，式（2.30）只需计算1次函数值fk，而fk-1，fk-2，fk-3已经由前面的计算得到；而式（2.31）则需要计算k1，k2，k3，k4，相当于计算4次函数f（t，y（t））的值。因此，在相同精度条件下，多步法比单步法快。

多步法的主要缺点是不能自启动，必须用其他方法求初始几步的值。例如，式（2.30）在计算y1时，需要用到y0，f0，f-1，f-2和f-3，从常微分方程初值问题式（2.5）仅能直接得到y0并计算出f0，而f-1，f-2和f-3必须用其他方法求取。另外，多步法不容易实现变步长运算。

2.2.2 显式算法和隐式算法

如果数值积分算法在计算yk+1 时所用到数据均已求出，则称其为显式算法。例如，式（2.10）、式（2.21）都是显式算法。

如果数值积分算法的右端中隐含有未知量yk+1，则称其为隐式算法。例如，四阶Adams隐式法（简记为AM4法）

就是隐式算法。

隐式算法在计算yk+1时，要用到的值，因此，必须用一个显式算法估计一个初值，然后再用隐式算法进行迭代运算，这就构成了预估-校正算法。当精度要求较高时，可广泛使用四阶Adams预估-校正法。它是用AB4法预估一个，然后将该预估值代入式（2.31）中进行校正，其计算公式为

式中，，。

由以上分析可见，单步法和显式算法在计算上比多步法和隐式算法更方便，但有时为了满足精度、稳定性等方面的要求，需要采用隐式算法。

2.2.3 截断误差和舍入误差

在分析数值积分算法的精度时，通常用泰勒级数作为工具。假设前一步求得结果是准确的，即有

则用泰勒级数求得在tk+1处的解析解为

不同的数值积分算法相当于在上式中取了不同项数之和而得到的近似解。欧拉法是从解析解中取前两项之和来近似计算yk+1的；R K4法则取前5项之和来近似计算yk+1。这种由数值积分算法单独一步引进的附加误差称为局部截断误差。它是数值积分算法给出的解与微分方程的解析解之间的差，故又称为局部离散误差。步长h越小，局部截断误差就越小。

不同的数值积分算法，其局部截断误差不同。若一种数值积分算法的局部截断误差为O（hr+1）（即相当于在式（2.34）中取了前r+1项之和），则该算法为r阶的。所以算法的阶数可以作为衡量算法精确度的一个重要标志。可见，欧拉法的局部截断误差为O（h2），具有一阶计算精度；RK2法的局部截断误差为O（h3），具有二阶计算精度；RK4法的局部截断误差为O（h5），具有四阶计算精度。

以上分析是在假设前一步所得结果是准确的前提下得出的，即式（2.33）成立时，有

但在求解过程中，实际上只有当k=0时，式（2.33）才成立。而当k=1，2，3，…时，式（2.33）并不成立，因而，采用r阶算法求得的解的实际误差要大一些。

设yk是在无舍入误差情况下由r阶算法算出的微分方程式（2.5）的近似解，y（t）为微分方程的解析解，εk=y（tk）-yk为算法的整体截断误差，则可以证明，εk=O（hr）。这说明整体截断误差比局部截断误差低一阶。

由于计算机的字长有限，可表示的数字个数也有限，在计算过程中不可避免地会产生舍入误差。舍入误差与步长h成反比。如果步长h小，计算次数多，则舍入误差就大。产生舍入误差的因素较多，除了与计算机字长有关以外，还与软件使用的数制系统、数的运算次序及计算函数f（t，y（t））的值所用子程序的精度等因素有关。

2.2.4 数值积分算法的计算稳定性

1.计算稳定性的概念

连续系统数值积分算法的仿真，实质上就是将微分方程（组）变换成差分方程（组），然后从初值开始，进行递推计算的过程。因此，采用数值积分算法对稳定系统进行仿真时，应该保持原系统的稳定性。然而，这一点能够确保吗？先看下面的例子。

【例2.2】 已知微分方程及初值

试比较在取不同步长h时，其解析解与数值积分算法解之间的差异。

【解】 该初值问题的解析解为

对应的结果曲线如图2.2（a）所示。

取h=0.1，0.075，0.05，分别用欧拉法和RK4法计算t=1.5处的y（t），所得结果曲线如图2.2（b）～（g）所示。

从图2.2可以看出，解析解单调下降并迅速收敛到0。

图2.2 例2.1的解析解和数值解

当h=0.1时，欧拉法和RK4法的解曲线均发散，数值积分算法的解是错误的。

当h=0.075时，欧拉法的解曲线仍然发散，对应的解是错误的；RK4法的解曲线单调下降并收敛到0，对应的解是正确的。

当h=0.05时，欧拉法和RK4法的解均收敛到0（虽然欧拉法的解曲线是振荡收敛的）。如果只要求得到t=1.5处y（t）的解，则两种数值积分算法的解都可以认为是正确的。事实上，此时有

解析解y（1.5）=9.5417286 × 10-2 1

欧拉法y（1.5）=3.1044085 × 10-1 0

RK4法y（1.5）=4.2522715 × 10-1 8

尽管解析解和数值积分算法的解在数量级上相差甚大，但从工程意义上讲，它们都是0。

为什么会出现上述情况呢？这是因为数值积分算法只是一种近似积分方法，它在反复的递推计算中会引进误差（这种误差通常是由初始数据的误差及计算过程中的舍入误差产生的）。如果误差的累积越来越大（通常是因为步长h的选择不合理，见例2.2），将使计算出现不稳定，从而得到错误的结果。此外，原系统的稳定性与计算稳定性是两个不同的概念。前者用原系统的微分方程、传递函数来讨论，后者则用逼近微分方程的差分方程来讨论。由于选用的数值积分算法不同，即使对于同一系统，差分方程也各不相同，计算稳定性也就各不一样。

如何分析数值积分算法的计算稳定性呢？对高阶微分方程或非线性微分方程做全面分析是非常困难的，通常用一个简单的一阶微分方程来考查数值积分算法的计算稳定性。

微分方程及初值问题

称为测试方程（TestEquation）。式中，λ=σ+jω为方程的特征根。根据稳定性理论，当特征根λ位于左半s平面，即Reλ=σ＜0时，原系统稳定。此时相应的数值积分算法的计算稳定性如何呢？

2.欧拉法的计算稳定性

用欧拉法对式（2.37）进行计算，有

按上式计算出来的yk+1并不是在t=（k+1）h处y（t）的真值，而是真值的一个包含了各种误差的近似值。随着递推次数的增加，这些误差是否会不断扩大，从而使yk+1完全不能表示此时的真值y（tk+1）呢？为了简化问题的讨论，假定仅仅在初始值y0处引入了初始误差ε0，而在递推计算过程没有引进任何其他误差。显然，此时yk+1的误差εk+1仅由yk的误差εk引起，所以有

用式（2.39）减去式（2.38），得到误差方程

依次类推，有

当|1+λh|＞1时，|εk+1|＞|εk|，表明若在递推计算过程中的某步引入了误差，则随着递推步数的增加，这个误差将逐渐扩大，最终导致差分方程的解完全失真。

反之，当|1+λh|≤1时，|εk+1|≤|εk|，表明随着递推步数的增加，引入的误差会被逐渐减小或保持有界。在这种情况下，称差分方程是计算稳定的。

显然，对于欧拉法而言，合理地选择步长h使其满足|1+λh|≤1是保持其计算稳定性的充要条件。

在例2.2中，λ=-30。当取h=0.1（或0.075）时，|1+λh|=2（或1.25）＞1，因而计算不稳定，相应的解完全失真。而当取h=0.05时，|1+λh|=0.5＜1，因而计算是稳定的，相应的解未完全失真，可以被接受。

令=λh，由|1+λh|≤1可知欧拉法在实轴上的稳定区域为（-2，0），即有

因此，|Reλ|=|σ|越大，步长h就应该取得越小。

为了保证计算稳定性而对步长h有所限制的数值积分算法称为条件稳定算法。在任何步长h＞0的情况下，计算都稳定的数值积分算法称为绝对稳定算法（亦称为无条件稳定算法或恒稳算法）。欧拉法是一种条件稳定的算法。

3.龙格-库塔法的计算稳定性

将RK法应用于测试方程式（2.37），容易得到下列误差方程

令=λh，代入上式，可得该式稳定的条件为

根据式（2.43），可得如表2.3所示的各阶RK法在实轴上的稳定区域。

表2.3 RK法在实轴上的稳定区域

从表2.3中可以看出，显式RK法都是条件稳定算法。对于条件稳定的数值积分算法，步长h的大小除了与所选用的算法有关外，还与方程本身的性质有关。从R K4法的实负稳定域（-2.785，0）可知，系统特征值的模|λ|越大，允许的步长h就越小。

在例2.2中，λ=-30。当取h=0.1时，λh=-3，超出R K 4法的稳定区域，因而计算不稳定，相应的解完全失真；而当取h=0.075（或0.05）时，λh=-2.25（或-1.5），落在稳定区域内，因而计算是稳定的，对应的解在相当大的程度上反映了真解，可以被接受。

对于实际系统，由于其特征值不一定为实数，因此，满足式（2.43）的也是复数。一般而言，步长h必须满足不等式

4.多步法的计算稳定性

同样，可以将多步法应用于测试方程式（2.37），得到类似的误差方程，只不过分析过程比较复杂。下面不加证明地给出如下结论：

● 在同阶的多步法中，隐式算法的稳定区域比显式算法的稳定区域大得多；

● 随着阶数的增加，多步法的稳定区域逐渐减小。

表2.4给出了各阶Adams法在实轴上的稳定区域，也就是h—=λh所允许的区域。

表2.4 Adams法在实轴上的稳定区域

由表2.4可知，除了AM1法和AM2法（隐式Adams法）为绝对稳定算法外，其他算法都是条件稳定算法。也就是说，步长h必须满足不等式

式中，M为由积分算法确定的常数。

2.2.5 数值积分算法的计算精度、速度、稳定性与步长的关系

由前面的分析可知，一般数值积分算法的计算稳定性与所选用的步长有关。在数字仿真中，步长h的选取是一个十分关键的因素。

各种数值积分算法对应的差分方程是对原微分方程的近似，存在着明显的截断误差；并且由于计算机及软件数制的字长有限，计算过程只能限制在有限的位数。因此，在递推计算过程中将引入舍入误差。步长h取得过大，会带来较大的截断误差，甚至会出现计算不稳定的现象；步长h取得过小，又会增加计算步数，而使舍入误差累积，总的误差也会变大。所以，仿真的总误差与步长h的关系不是单调函数关系，而是一个具有极值的函数关系，如图2.3所示。显然，并不是步长h越小，精度越高，而是存在一个最佳步长。为了加深对这个问题的认识，下面来看一个例子。

图2.3 步长和误差的关系曲线

【例2.3】 考虑如下二阶系统

在0≤t≤10上的数字仿真，并将不同步长下仿真结果与解析解进行精度比较。

【解】 该微分方程的解析解分别为

y（t）=100cost（当R=0）

采用RK4法进行计算，选择状态变量

则有如下状态空间模型及初值条件

采用RK4法进行计算，编制文件名为exam2-3.m的程序如下：

        % 这是例2.3的仿真程序
        clear
        h=input（′请输入步长h=′）；         % 输入步长h
        M=round（10/h）；                  % 置总计算步数
        t（1）=0；                        % 置自变量初值
        y-0（1）=100；y-05（1）=100；         % 置解析解的初值
        x1（1）=100；x2（1）=0；             % 置状态向量的初值
                                                      %（y-0和y-05分别对应于R=0和R=0.5）
        y-rk4-0（1）=x1（1）；y-rk4-05（1）=x1（1）；% 置数值解的初值
        % 求解析解
        fork=1：M
          t（k+1）=t（k）+h；
            y-0（k+1）=100*cos（t（k+1））；
            y-05（k+1）=100*exp（-t（k+1）/2）.*cos（sqrt（3）/2*t（k+1））+...
          100*sqrt（3）/3*exp（-t（k+1）/2）.*sin（sqrt（3）/2*t（k+1））；
        end
        % 采用RK4法求解
        %R=0
        fork=1：M
            k11=x2（k）；k12=-x1（k）；
            k21=x2（k）+h*k12/2；k22=-（x1（k）+h*k11/2）；
            k31=x2（k）+h*k22/2；k32=-（x1（k）+h*k21/2）；
            k41=x2（k）+h*k32；k42=-（x1（k）+h*k31）；
            x1（k+1）=x1（k）+h*（k11+2*k21+2*k31+k41）/6；
            x2（k+1）=x2（k）+h*（k12+2*k22+2*k32+k42）/6；
            y-rk4-0（k+1）=x1（k+1）；
        end
        %R=0.5
        fork=1：M
                k11=x2（k）；k12=-x1（k）-x2（k）；
                k21=x2（k）+h*k12/2；k22=-（x1（k）+h*k11/2）-（x2（k）+h*k12/2）；
                k31=x2（k）+h*k22/2；k32=-（x1（k）+h*k21/2）-（x2（k）+h*k22/2）；
                k41=x2（k）+h*k32；k42=-（x1（k）+h*k31）-（x2（k）+h*k32）；
                x1（k+1）=x1（k）+h*（k11+2*k21+2*k31+k41）/6；
                x2（k+1）=x2（k）+h*（k12+2*k22+2*k32+k42）/6；
                y-rk4-05（k+1）=x1（k+1）；
            end
            % 求出误差最大值
            err-0=max（abs（y-0-y-rk4-0））；
            err-05=max（abs（y-05-y-rk4-05））；
            % 输出结果
            disp（′最大误差（R=0）最大误差（R=0.5）′）
            err-max=[err-0，err-05]；
            disp（err-max）

步长h选取为0.0001到0.5之间变化，并且将数值解与解析解进行比较，求出最大误差数据列于表2.5中。

表2.5 RK4法的步长与误差关系

从表2.5中可以看出，当步长h=0.001时，总误差最小；当步长h＜0.001时，由于舍入误差变大而使总误差增加；当步长h＞0.001时，则由于截断误差的增加也使得总误差加大。另外，当系统的解变化剧烈时（如R=0），误差对步长h的变化较为敏感；当系统的解变化平稳时，步长h的变化对误差的影响要缓和得多。

此外，对于给定的仿真时间区间，步长h取得过小，会使得计算步数增加得太多，仿真过程太长。

因此，当数值积分算法确定以后，需要仔细地选择适当的步长h。在实际仿真中，考虑得比较多的一个重要因素是被仿真系统的动态响应特性。如果系统的动态响应较快，导数的变化较为剧烈，步长h应当取得小一些；反之，步长h可以取得大一些。为了保证计算的稳定性，步长h应限制在最小时间常数Tmin（即），假定系统是稳定的，n为系统阶数）的数量级。例如，对于RK4法，步长h的选择至少应满足

而为了保证足够的计算精度，在实际使用中，所选择的步长h还要更小一些。

对于一般工程系统的分析和设计而言，仿真结果的误差不超过0.5%就可以满足精度要求。经验表明，对于RK4法，通常可以根据系统方程中的最小时间常数来选择步长h，一般取

或者根据系统中反应最快的小闭环的开环剪切频率ωc来选择，一般取

值得一提的是，上述两种选择步长h的经验公式都是针对RK4法而言的。如果要用RK2法进行仿真，则为了达到相同的精度，步长h一般应取为RK4法的。虽然RK4法每步要计算4次微分方程右端函数f（t，y（t））的值，而RK2法只需要计算2次，但为了达到同样的精度，前者的步长h通常可以取为后者的10倍（在保证计算稳定性的前提下）。显然，前者的计算速度为后者的5倍。如果采用欧拉法，则为了达到相同精度，步长h还要取得更小。这就是控制系统数字仿真中通常选用RK4法的重要原因之一。

2.2.6 数值积分算法的选择原则

在对连续系统进行数字仿真时，存在一个数值积分算法选择问题。由于数值积分算法与仿真计算的精度、速度、误差积累、计算稳定性、自启动能力等因素都有密切的关系，因此，正确选择数值积分算法是一个十分重要而复杂的问题。到目前为止，还没有一个普遍适用的规律，这里仅给出一些选择时应考虑的原则。

1.精度要求

数值积分算法的仿真精度主要受截断误差、舍入误差和误差累积的影响，它们与积分算法、步长h、计算时间及所用的计算机精度等有关。在步长h相同的条件下，积分算法的阶数越高，截断误差越小；另外，多步法的精度比单步法高，而其中同阶的隐式算法的精度又高于显式算法。所以，当需要进行高精度仿真时，可以采用高阶的多步隐式算法和较小的步长h。然而，步长h的减小会增加递推次数，不仅增加计算工作量，而且会增大舍入误差和累积误差。因此，应当根据需要的精度，合理地选择积分算法和阶数。在具体运算时，并不是说算法的阶数越高，步长h越小，效果就越好。经验表明，低精度问题最好用低阶算法来处理。如果选用高阶算法，则应在保证计算稳定性的前提下，把步长h取得大一些。

2.计算速度

计算速度取决于在给定积分时间内的计算步数和每步计算所需的时间，而每步的计算时间又与积分算法有关，它主要取决于微分方程（组）右端函数f（t，y（t））（或f（t，y（t）））值的计算次数。例如，RK4法每步要计算4次右端函数f（t，y（t））（或f（t，y（t）））的值；而四阶Adams预估-校正法每步只要计算2次右端函数f（t，y（t））（或f（t，y（t）））的值，计算速度几乎快了一倍。因此，在右端函数f（t，y（t））（或f（t，y（t）））复杂、计算量大而精度要求又高的问题中，可以采用Adams预估-校正法。同时，为了加快计算速度，在积分算法选定后，应在保证精度的前提下，尽量选择较大的步长h，以减少计算步数。

3.计算稳定性

保证数值解的计算稳定性，是进行数字仿真的先决条件；否则，计算结果将失去真实意义。从计算稳定性的角度来看，同阶的RK法优于显式Adams法，但又不如隐式Adams法；而三阶以下的隐式Adams法具有较好的稳定性（参见表2.4）。所以，最好避免使用显式Adams法。

4.自启动能力

单步法具有自启动能力。多步法没有自启动能力，必须借助于单步法启动运算后才能开始递推计算，因此实现起来要困难一些。一般简单的仿真程序多采用单步法。

5.步长变化能力

单步法在整个仿真计算过程中，步长h可以在一定的范围内变化，而多步法对步长h的变化有严格的要求。如果要求仿真时步长h可变，最好选用单步法。

总之，数值积分算法的选择与多种因素有关，而各种因素之间又相互影响，究竟选择哪一种算法，要根据具体情况而定。经验表明，当微分方程（组）右端函数的计算量不大而精度要求又不很高时，RK4法是一种很好的选择。在控制系统的数字仿真中，最常见的是线性系统的状态方程，右端函数中绝大部分是矩阵与向量的乘积（即线性运算）。因此，对于一般控制系统的仿真，通常都采用RK4法。

2.2.7 误差估计与步长控制

2.2.5节中介绍了控制系统仿真中选择步长h的经验公式。但遗憾的是，高阶系统的Tmin和ωc有时是很难估计的，或者是由于系统的非线性，或者有时根本就无法估计出上述性能指标。另外，系统中最小时间常数对应的极点只影响到过渡过程起始段的形态，而系统过渡过程主要是由那些靠近虚轴的主导极点所决定的。固定步长积分算法的计算步长h是按照起始段来选取的，这就不可避免地造成后面阶段由于使用过小步长积分而引起计算量的增加和时间的浪费，从而使得计算速度下降。为此可以采用变步长策略。一个成熟的仿真程序或仿真语言通常将步长h的自动控制作为必要的手段。

实现步长自动控制的前提是要有一个好的局部截断误差估计公式。下面仅对RK法的误差估计与步长控制策略作一简单介绍。

为了估计一种RK法的误差，通常采用的方法是设法在原算法公式的基础上找到另一个低阶（一般是低一阶）的RK公式，要求这两个公式中的ki相同，则两个公式计算结果之差就可以作为局部截断误差估计。

第一个四阶变步长RK法是R.H.Merson于1957年给出的，称之为龙格-库塔-默森法（简称为RKM34法），四阶RKM公式为

三阶RKM公式为

式中，k1，k3，k4由式（2.50）定义。

用作为局部截断误差估计，即有

采用RKM34法可以得到四阶精度的计算结果和三阶精度的误差估计。该法是目前被广泛应用的一种数值积分算法，它在实轴上的稳定域为（-3.548，0）。其缺点是计算量较大，每步需要计算5次微分方程右端函数f（t，y（t））的值。

与RKM34法类似的还有E.Fehlberg于1969年提出的RKF45法（即龙格-库塔-费尔别格法），它每步计算6次微分方程右端函数f（t，y（t））的值，获得五阶精度的计算结果和四阶精度的误差估计，被公认为是对非病态系统进行仿真最为有效的算法之一。1978年，L.F.Shampine提出了RKS34法（即龙格-库塔-夏普勒法），它每步只计算4次微分方程右端函数f（t，y（t））的值，却能获得四阶精度的计算结果与三阶精度的误差估计。这两种算法的稳定性与RK4法相差不大。

有了误差估计以后，就可以实现步长的自动控制。通常采用的步长控制策略是“加倍-减半法”。

设定最小允许误差εmin和最大允许误差εmax，当估计的局部截断误差大于最大允许误差时，将步长减半，并重新计算本步；当误差在最大允许误差和最小允许误差之间时，本步计算结果有效，下步步长不变；当误差小于最小允许误差时，本步计算结果有效，下步步长加倍。

每步的局部误差通常取以下形式

式中，E k为利用前述的R K M 34法等计算出来的误差估计。由式（2.53）可知，当|yk|较大时，ek是相对误差；当|yk|较小时，ek是绝对误差。这样做的目的是避免当|yk|的值很小时，ek变得过大。

步长控制策略可以用下式来表示

在实际应用时，通常还设置步长下限hmin和步长上限hmax。在步长减半的过程中，如果步长小于步长下限hmin，则不再减半；否则将增加仿真时间，并增大舍入误差。同样，在步长加倍的过程中，如果步长大于步长上限hmax，则不再增大步长；否则将增加截断误差，减小计算稳定性。

这种步长控制策略可以很容易地推广到求解向量微分方程。由2.1.4节可知，只需要将RKM34法等应用到式（2.23）的求解中，计算出

并将ek修改为

即可。

上述误差估计和控制步长的策略，虽然增加了每步的计算量，但从总体上考虑往往是合算的。它可以很好地解决仿真中计算精度与计算量之间的矛盾，尤其在系统特征根分散程度较大的情况下效果尤为明显。

2.2.8 数值积分算法仿真实例

在MATLAB/Simulink环境下，利用数值积分算法对系统进行仿真的途径有两种。对于以微分方程（组）给出的数学模型，通常用MATLAB语言编程实现；而对于控制系统分析和设计中最常见的以系统结构图描述的数学模型，则采用Simulink实现更为方便。

1.数值积分算法仿真的编程实现

用MATLAB语言编程实现仿真的主要步骤是调用MATLAB中的ODE（Ordinary Differential Equation）解函数。MATLAB提供的常用ODE解函数如下：

●ode45此算法被推荐为首选算法；

●ode23这是一个比ode45低阶的算法；

●ode113用于更高阶或大的标量计算；

●ode23t用于解决难度适中的问题；

●ode23s用于解决难度较大的问题；

●ode15s与ode23相同，但要求的精度更高；

●ode23tb用于解决难度较大的问题。

这些ODE解函数的调用格式基本相同。例如，ode45的基本调用格式为

        [t，x]=ode45（′方程函数名′，tspan，x0，tol）；

其中，方程函数名为描述系统状态方程的M函数的名称；tspan一般为仿真时间范围（例如，取tspan=[t0，tf]，t0为起始计算时间，tf为终止计算时间）；x0为系统状态变量初值，应使该向量的元素个数等于系统状态变量的个数；tol用来指定精度，其默认值为10-3（即0.1%的相对误差），一般应用中没有必要修改其默认属性，可以直接采用默认值。函数返回两个结果：t向量和x阵（注意：MATLAB中认为所有变量都是矩阵，没有黑体标志。本书在介绍MATLAB函数时，用英文字母的正体配合文字说明表示向量和矩阵，下同）。由于计算中采用了步长自动控制策略，因而t向量不一定是等间隔的。但是，仿真结果可以用指令plot（t，x）绘制出来。

【例2.4】 某地区某病菌传染的系统动力学模型为

式中，x1（t）表示可能受到传染的人数，x2（t）表示已经被传染的病人数，x3（t）表示已治愈的人数。试调用ode45函数进行编程，对其进行仿真研究，并绘制出对应的时间相应曲线。

【解】 采用MATLAB语言进行编程，文件名为exam2_4.m。程序如下：

        % 这是例2.4的仿真程序
        clear
        x0=[6201070]；                         % 置状态变量初值
        tspan=[0，30]；                         % 置仿真时间区间
        [t，x]=ode45（′fun2-4′，tspan，x0）；            % 调用ode45求仿真解
        plot（t，x（：，1），′k′，t，x（：，2），′k--′，t，x（：，3），′k：′）；  % 用不同的线型绘制仿真结果曲线
        xlabel（′time（天）t0=0，tf=30′）；              % 对t-x轴进行标注
        ylabel（′x（人）：x1（0）=620，x2（0）=10；x3（0）=70′）；
        legend（′x1′，′x2′，′x3′）；
        grid；

其中，fun2-4是描述微分方程组的M函数文件，具体程序如下：

        functionxdot=fun2-4（t，x）
        xdot1（1）=-0.001*x（1）*x（2）；                % 第一个微分方程
        xdot1（2）=0.001*x（1）*x（2）-0.072*x（2）；       % 第二个微分方程
        xdot1（3）=0.072*x（2）；                      % 第三个微分方程
        xdot=xdot1′；

运行exam2_4.m后，得到如图2.4所示的曲线。

图2.4 病菌传染模型的仿真结果

从图2.4可以看出，随着时间的推移，受到病菌威胁的人数和被传染的病人数逐渐减少，而治愈的人数逐渐增加。这一点与病菌传染的医学统计结果是吻合的。

2.Simulink模型的构建与仿真

对于以系统结构图描述的数学模型，采用Simulink构模并仿真是十分方便的。在大多数情况下，用户甚至无须编制一条仿真程序就可以进行仿真。

（1）Simulink模型的构建

为了能用Simulink对系统进行仿真，首先要在Simulink环境下打开一个空白模型窗口；然后依据系统结构图给定的环节和信号流程，从Simulink模块库的各个子库中选择相应的模块，并用鼠标左键将它们拖入模型窗口；双击选择的模块，设置需要的参数；对各模块进行连接，这就构成了需要的Simulink模型（即仿真结构图，参见1.6.2节例1.2）。

（2）仿真参数的设置和ODE算法的选择

构建好系统的Simulink模型之后，接下来的事情就是运行模型，得出仿真结果。而设置仿真参数和选择ODE算法则是保证仿真有效进行的前提。Simulink模型实际上是一个计算机程序，它定义了描写被仿真系统的一组微分或差分方程。当选中模型窗口中的“仿真启动”图标后，Simulink就开始用选定的数值计算方法去求解方程。

在进行仿真前，如果用户不采用系统默认的仿真设置，就必须对各种仿真参数进行配置（Configuration），其中主要包括：仿真的起始时间和终止时间的设定；计算步长h的选择；各种仿真容差的选定；解算器的选择。

在模型窗口（参见图1.14）的Simulation菜单下选择其中的仿真参数配置子菜单（Config-uration Parameters），就会弹出一个仿真参数配置窗口，如图2.5所示。

图2.5 Simulation仿真参数配置界面

在图2.5的对话框中有若干个选项，常用的选项为仿真时间（Simulation time）和解算器选项（Solver options）。对于一般用户而言，设置这两个选项中的相应内容就可以完成求解算法及仿真参数的设置。

① 求解算法的选择

在解算器选项的解算器（Solver）下拉框中可以选择不同的解算器（即求解算法），定步长（Fixed-step）下支持的算法如图2.6（a）所示，变步长（Variable-step）下支持的算法如图2.6（b）所示。一般情况下，连续系统仿真应该选择定步长的ode4（即RK4）算法或者ode45变步长算法（默认设置），而离散系统一般默认地选择定步长的discrete（no continuousstates）算法。

图2.6 Simulink求解算法的选择

② 计算步长的选择

定步长算法的步长h可以通过在图2.6（a）中Fixed-step size（fundamental sample time）框输入需要的步长参数进行选择，也可以依赖计算机自动选择步长h（即采用默认设置）；而对于变步长算法，则建议最大步长（Max step size），最小步长（Min step size）和初始步长（Initial step size）都使用默认（即auto）设置（见图2.6（b））。

③ 仿真时间的设置

在图2.5的相关选项中还可以修改仿真的初始时间和终止时间。这里的时间概念与真实时间并不一样，只是计算机仿真中对时间的一种表示。比如，10s的仿真时间，如果步长h定为0.1，则需要执行100步，若把步长h减小，则计算点增加，实际的执行时间就会增加。一般仿真开始时间设为0，而结束时间视不同的因素而选择。

（3）Simulink仿真结果输出

完成了仿真参数的设置和ODE算法的选择后，就可以启动仿真。Simulink会自动将系统结构图转换成状态空间模型并调用所选择的算法进行计算。为了得到所需要的仿真结果，除了可以直接采用Scope模块显示仿真结果曲线外，还可以将仿真结果数据传送到MATLAB工作空间中，利用plot函数绘制相应的图线。

【例2.5】 考虑如图2.7所示的直流电机拖动系统，试研究外环PI控制器对系统阶跃响应的影响。

【解】 构建系统的Simulink模型（文件名：exam2-5.mdl），如图2.8所示。

启动仿真后，可以立即得出仿真结果，该结果将自动返回到MATLAB工作空间中，其中时间变量名默认为tout，输出信号的变量名默认为yout。在MATLAB指令窗中，运行指令：

图2.7 直流电机拖动系统

图2.8 直流电机拖动系统的Simulink模型

        ≫plot（tout，yout，′k′）；grid；

则可以得到如图2.9（a）所示的阶跃响应曲线。显然，该曲线不很理想，超调量较大。为此，可以将外环的PI控制器参数调整为，并分别选择ɑ=0.17，0.5，1，1.5，则可以得到如图2.9（b）所示的仿真结果。可以看出，如果选择PI控制器为的控制效果。，就能够得到较为满意

图2.9 直流电机拖动系统的阶跃响应

Simulink除了能将用系统结构图描述的数学模型进行建模仿真外，也可以把微分方程模型用图形表示。

【例2.6】 考虑著名的VandePol方程

试绘制其相轨迹（0≤t≤20）。

【解】 选择状态变量

则有如下非线性状态方程及初值

第1个方程可以看成将x2（t）信号作为一个积分器的输入，积分器的输出将成为x1（t）信号。同样，x2（t）信号本身也可以看成一个积分器的输出，而积分器的输入端信号应该为，于是，利用Simulink提供的各种模块可以得到如图2.10所示的仿真模型（文件名：exam2-6.mdl）。

图2.10 VanderPol方程的Simulink模型

可以看出，Simulink模型中，除了有各个模块及其连接之外，还给出了各个信号的文字描述。在Simulink模型中加入文字描述的方式很简单，在想加文字说明的位置双击鼠标，则将出现字符插入的标示，这时可以将任意的字符串写到该位置。文字描述写到模型中后，可以用鼠标单击并拖到指定位置。

此外，Simulink模型中有些模块需要将输入端和输出端（通常用于反馈路径）调换方向。为此，可以用鼠标单击需要调换方向的模块并选中它，则选中模块的4个角出现黑点，表明它处于选中的状态。然后，打开Simulink的Format菜单（见图2.11），选择其中的翻转项（Flip block）即可。

为了进一步绘制相轨迹工作的需要，本例的仿真结果数据将输出到MATLAB工作空间中，故使用了3个To Workspace模块。需要注意的是：在该模块的参数设置中，必须在Save format下拉框中选择Array，同时在Variable name框中输入所需变量名，如图2.12所示。

在图2.5中的Stop time框中输入仿真结束时间：20，并选择定步长的ode4算法（步长h=0.01）。启动仿真后，仿真结果将赋给MATLAB工作空间内的变量t，x1和x2。在MATLAB指

图2.11 Format子菜单

图2.12 To Workspace模块的对话框

令窗口中运行指令：

        ≫plot（x1，x2，′k′）；grid；

得到该方程的相轨迹曲线，如图2.13所示。

图2.13 Vanderpol方程的相轨迹

从这个例子可以看出，微分方程实际上是可以由Simulink用图示的方法完成的，因而可以将这样的思想应用于更复杂系统的建模。在本例的仿真中，还涉及两个积分器的初始值，它们可以在Integrator模块的参数Initial condition框中设置。

（4）仿真结构的参数化

Simulink模型中的参数可以是实际数值，也可以是用字母表示的变量名。用字母表示的变量可以在MATLAB的工作空间中赋值，或用M文件赋值，然后进行仿真计算。

【例2.7】 含有磁滞回环非线性环节的控制系统如图2.14（a）所示，其中，磁滞回环非线性环节的特性如图2.14（b）所示。试研究该非线性环节对系统性能的影响（0≤t≤3s，r（t）=2·1（t））。

图2.14 非线性控制系统

【解】 构建系统的Simulink模型（文件名：exam2_7.mdl），如图2.15所示。为了研究问题的方便起见，将Backlash模块的参数Deadband width设置为C1。

图2.15 磁滞回环控制系统的Simulink模型

在MATLAB指令窗口依次运行C1的不同值（C1=0.5，1，2）的指令后启动仿真，并在仿真结束后在MATLAB指令窗口中运行指令：

        ≫plot（tout，yout，′k′）；

可以得到如图2.16所示的仿真曲线。显然，不同磁滞宽度的阶跃响应曲线的形状并不相同。需要注意的是，为了将不同C1值对应的阶跃响应曲线绘制在同一坐标轴下，在第一次绘制图形后，应该在MATLAB指令窗口中运行指令：

        ≫holdon

图2.16 磁滞回环控制系统的阶跃响应

此外，如果输入信号的幅值增大或减小，则原来系统响应曲线的形状将可能出现不同。如图2.17所示为C1=1时，输入幅值等于3和0.6的仿真结果。显然，这一点与线性系统是不相同的。

图2.17 输入幅值改变后的仿真结果

（5）与M函数的组合仿真

如果Simulink模型中存在复杂的非线性环节或复杂的逻辑运算，而在MATLAB提供的所有工具箱中都找不到相应的模块，读者可以自己编制一个M函数，嵌入Simulink模型中。

【例2.8】 某非线性系统如图2.18所示，试求r（t）=2·1（t）时系统的动态响应。

图2.18 饱和非线性系统

【解】 构建系统的Simulink模型（文件名：exam2-8.mdl），如图2.19所示。为了研究问题的方便起见，不采用Discontinuities子库中的Saturation模块，而选择User-Defined-Functions子库中的MATLAB Fcn模块，并将参数MATLAB Function设置为：saturation_zone。

图2.19 带有M函数的非线性系统的Simulink模型

然后编制文件名为saturation_zone.m的M函数文件，程序如下：

        % saturation-zonefunction
        function[uo]=saturation-zone（ui）
        if ui＞=1
                uo=1；
          elseif ui＜=-1
                uo=-1；
          else
                uo=ui；
          end

启动仿真后，得到如图2.20所示的仿真结果。

图2.20 饱和非线性系统的阶跃响应