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 飽和非線性系統的階躍響應