2.1 数值积分算法_控制系统计算机仿真-QQ阅读男生玄幻网

上QQ阅读APP看本书，新人免费读10天

设备和账号都新为新人

2.1 数值积分算法

2.1.1 数值积分算法的基本原理

连续系统仿真的数值积分算法是利用数值积分法将常微分方程（组）描述的连续系统变换成离散形式的仿真模型——差分方程（组），数值积分算法就是对一阶微分方程近似求解的公式。为了能在计算机上进行求解，首先要把被仿真系统的数学模型表示为一阶微分方程组或状态空间模型。

假设一阶向量微分方程及初值问题为

式（2.1）可以写成一阶微分方程组和初值问题

式（2.1）在t=t0，t1，…，tk，tk+1处的解析解为

用公式

来代替解析解，其中yk+1，yk，qk分别是y（tk+1），y（tk），的近似解。这样就从一阶向量微分方程出发建立了离散形式的仿真模型，即向量差分方程式（2.4）。因此，所谓数值解法就是寻求初值问题式（2.1）在一系列离散时间点t1，t2，…，tk，tk+1处的近似解y1，y2，…，yk，yk+1（即数值解）。相邻两个离散时间点的间距h=tk+1-tk称为计算步长或计算步距（简称为步长或步距）。本书中如果不特别说明，总是假定步长h为定值。

由式（2.4）可知，数值积分算法的主要问题归结为如何对向量函数f（t，y（t））进行数值积分，求出f（t，y（t））在区间（tk，tk+1）上定积分的近似解qk。求解初值问题的数值方法的共同特点是步进式的（即所谓的递推式）。

采用不同的方法求解qk，会推导出不同的数值积分算法。常用的基本算法可以分为单步法、多步法和预估-校正法3大类，而每类算法又可以分为显式算法和隐式算法两种类型。不同的算法对系统的求解精度、速度和数值稳定性等有着不同的影响。

为了讨论方便起见，先从标量形式开始讨论，然后再推广到向量形式。

考虑一阶微分方程及初值问题

与式（2.3）类似，式（2.5）的解析解为

用公式

来代替解析解，其中yk+1，yk，qk分别是y（tk+1），y（tk），的近似值。

2.1.2 欧拉法

欧拉（Euler）法是一种最简单的数值积分算法，由于其精度低，实际中很少使用，但常常用来说明数值积分的基本概念。对于方程式（2.5），在区间（tk，tk+1）上求积分，有

如果区间（tk，tk+1）足够小，则（tk，tk+1）上的f（t，y（t））变化很小，可以近似地看成常数f（tk，yk）。因此，可以用矩形面积hf（tk，yk）近似代替曲线积分，如图2.1所示，从而有

即有

qk=hf（tk，yk）

将式（2.9）写成差分方程形式

图2.1 欧拉积分

式中

fk=f（tk，yk）

式（2.10）就是著名的欧拉数值积分法公式（简称为欧拉法）。欧拉法计算简单，每步只需要计算1次函数f的值。但该算法的精度较低，原因在于将f（t，y）看成常数f（tk，yk），从而用矩形面积ɑbtk+1tk代替了准确的曲面面积ɑctk+1tk，产生了曲边三角形ɑcb这一块较大的面积误差，如图2.1所示。因此，欧拉法仅适用于步长h很小或f（t，y）变化较小的场合。

2.1.3 龙格-库塔法

龙格-库塔（Runge-Kutta）法是求解常微分方程初值问题式（2.5）的各种数值积分算法中应用得最广泛的一种。由于阶数的不同及数值积分公式中系数选择的不同，它可以包括许多不同的公式。众所周知，利用在t=tk处y（t）的函数值y（tk）及其各阶导数值，…，利用泰勒级数可以求出在t=tk+1处y（tk+1）的近似值，而且其计算精度随泰勒级数所取项数的增加而提高。但是，对于方程式（2.5），直接采用泰勒级数法计算函数f（t，y（t））在t=tk处的各阶导数值（即需要计算对应的各阶偏导数值），运用起来很不方便。德国数学家C.Runge和M. W.Kutta两人先后提出了间接利用泰勒级数法，即用若干个时间点处f（t，y（t））的函数值的线性组合来代替f（t，y（t））在t=tk处的各阶导数值，然后按泰勒公式展开确定其中的系数。这样既可以避免计算各阶偏导数值，又可以提高数值积分的精度。这就是龙格-库塔法（以下简称RK法）的基本思想。

RK法包含显式、隐式和半隐式等算法。根据控制系统计算机仿真中的应用情况，本书仅介绍显式RK法。

考虑方程式（2.5），将解析解在t=tk附近用泰勒级数展开，有

由于

所以，有

为了避免计算，等各阶偏导数值，可以令

式中，Wi为待定的权因子；r为使用的k值的个数（即计算函数f（t，y（t））的次数）；ki为不同时间点处函数f（t，y（t））的值；ci和ɑij为待定系数。

式（2.13）就是显式R K法的一般公式。将该式中的ki（i=1，2，…，r）用泰勒级数展开后与式（2.12）逐项对比确定待定权因子和系数，可以得到不同阶次的RK法。

当r=1时，只有一个

k1=f（tk，yk）

此时W1=1，即将式（2.12）中h2及h2以上项略去，得

当r=2时，式（2.13）可以表示为

将k1=fk代入k2中，并将k2在（tk，yk）附近用泰勒级数展开，可得

于是，有

将上式与式（2.12）逐项进行比较，令h，h2项的对应系数相等，得到以下关系

上述3个方程中有4个未知数W1，W2，c2，ɑ21，解不是唯一的。取c2为自由参数来确定W1，W2，ɑ21。如果取c2=1，有

相应的式（2.15）为

上式表示的数值解是基于y（t）的泰勒级数在二阶导数项以后截断求得的，故称为二阶RK法（简记为RK2法）。

当r=3时，仿照上述方法，令h，h2，h3项的系数与式（2.12）中对应项的系数相等，可得三阶RK法（简记为RK3法）。

当r=4时，可得到如下著名的四阶RK算法（亦称为经典RK算法，简记为RK4法）

平时提到的RK4法就是指上式而言的。这个算法是数字仿真中最常用的一种算法。该算法每步需要计算4次微分方程右端函数f（t，y（t））的值，计算量较大，但计算精度较高。由于四阶精度的计算结果对于一般工程问题而言已可满足精度要求，因而式（2.21）在仿真中得到了广泛的应用；并且在仿真中比较不同算法的计算精度时，也常常以RK4法的计算结果作为标准。

通过对RK法的讨论，可以将前面介绍的几种数值积分法统一起来。它们都可以看成是由y（t）在t=tk附近用泰勒级数展开而产生的，只不过是泰勒级数所取项数的多少不同而已。欧拉法只取前两项，RK2法取了前3项，RK4法取了前5项。从理论上讲，可以构造任意高阶的计算方法。但是，精度的阶数与计算微分方程右端函数f（t，y（t））值的次数之间的关系并非等量增加的，如表2.1所示。

表2.1 f（t,y（t））的计算次数与算法精度阶数的关系

由此可见，RK4法有其优越性，而四阶以上的RK法所需计算的f（t，y（t））值的次数要比阶数多，这将大大增加计算的工作量，从而限制了更高阶RK法的应用。

【例2.1】取h=0.2，试分别采用欧拉法、RK2法和RK4法求解微分方程初值问题

的数值解，并比较计算精度（解析解：。

【解】编制文件名为exam2_1.m的程序如下：

        % 这是例2.1的仿真程序
        clear
        t（1）=0；                                  % 置自变量初值
        y（1）=1；y-euler（1）=1；y-rk2（1）=1；y-rk4（1）=1；  % 置解析解和数值解的初值
        h=0.2；                                  % 置步长h=0.2
        % 求解析解
        fork=1：5
          t（k+1）=t（k）+h；
            y（k+1）=sqrt（1+2*t（k+1））；
        end
        % 利用欧拉法求解
        fork=1：5
            y-euler（k+1）=y-euler（k）+h*（y-euler（k）-2*t（k）/y-euler（k））；
        end
        % 利用RK2法求解
        fork=1：5
            k1=y-rk2（k）-2*t（k）/y-rk2（k）；
            k2=（y-rk2（k）+h*k1）-2*（t（k）+h）/（y-rk2（k）+h*k1）；
            y-rk2（k+1）=y-rk2（k）+h*（k1+k2）/2；
        end
        % 利用RK4法求解
        fork=1：5
            k1=y-rk4（k）-2*t（k）/y-rk4（k）；
            k2=（y-rk4（k）+h*k1/2）-2*（t（k）+h/2）/（y-rk4（k）+h*k1/2）；
            k3=（y-rk4（k）+h*k2/2）-2*（t（k）+h/2）/（y-rk4（k）+h*k2/2）；
            k4=（y-rk4（k）+h*k3）-2*（t（k）+h）/（y-rk4（k）+h*k3）；
            y-rk4（k+1）=y-rk4（k）+h*（k1+2*k2+2*k3+k4）/6；
        end
        % 输出结果
        disp（′   时间    解析解    欧拉法    RK2法    RK4法′）
        yt=[t′，y′，y-euler′，y-rk2′，y-rk4′]；
        disp（yt）

计算结果如表2.2所示。可以看出，计算精度从低到高的排列次序为欧拉法、RK2法和RK4法，而RK4法的计算精度明显高。

表2.2 例2.1的计算结果

2.1.4 微分方程数值积分的矩阵分析

前面介绍的各种数值积分法公式都是针对标量微分方程式（2.5）的，而实际系统中的大量仿真对象是用一阶向量微分方程式（2.1）或一阶微分方程组式（2.2）描述的。每个微分方程都可以用数值积分法来求解，这样就很容易将整个系统的动态响应全部解算出来。下面给出RK4法求解一阶向量微分方程式（2.1）的计算公式，其他算法的计算公式读者可以自己写出。

对于一阶向量微分方程及初值问题