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）的計算公式，其他算法的計算公式讀者可以自己寫出。

對于一階向量微分方程及初值問題

RK4法的向量形式為

具體就是

式中，i=1，2，…，n，n為方程組階次。

顯然，對于n維向量y（t），每步需要計算4n個kij的函數值。

在控制系統的仿真中，最常見的向量微分方程是線性定常系統（假設為SISO系統）的狀態方程

即有

R K4法的4個Ki可表示為

對于n維向量x（t），取

和

則4個求Ki的公式可以合并為一個公式

注意：在式（2.29）中，方括號表示矩陣運算，圓括號表示求函數u（t）在tk+h i時刻的值。每步的計算工作量主要是4次矩陣與向量的乘法（即4n2次乘法）。