在連續系統的仿真中，主要的計算工作是求解一階微分方程

y'=f(x,y)y（

x0）=

y0解析法只能用來求解一些特殊類型的方程，實際仿真問題中歸結出來的微分方程主要靠數值解法。由于實際運算只能完成有限項或有限步運算，因此要將有些需用極限或無窮過程進行的運算有限化，對無窮過程進行截斷，這樣產生的誤差成為截斷誤差。根據實際情況建立的數學模型往往難以求解。通常需要通過近似替代，將所求解的數學模型簡化為易求解的數值計算問題后再進行求解。數學模型的理論解與數值計算問題的精確解之間的誤差稱為截斷誤差。這是計算方法本身帶來的誤差，所以也成為方法誤差。2025/6/213得到高精度方法的一個直接想法是利用Taylor展開假設式

y'=f(x,y)(a≤x≤b)

中的

f(x,y)

充分光滑,將y(xi+1)在xi點作Taylor展開，若取右端不同的有限項作為y(xi+1)的近似值，就可得到計算y(xi+1)的各種不同截斷誤差的數值公式。例如：取前兩項可得到3.2龍格－庫塔方法2025/6/214其中P階泰勒方法若取前三項，可得到截斷誤差為O(h3)的公式類似地，若取前P+1項作為y(xi+1)的近似值，便得到2025/6/215顯然p=1時,

yi+1=yi+hf(xi,yi)它即為我們熟悉的Euler方法。當p≥2時,要利用泰勒方法就需要計算f(x,y)的高階微商。這個計算量是很大的，尤其當f(x,y)較復雜時，其高階導數會很復雜。因此，利用泰勒公式構造高階公式是不實用的。但是泰勒級數展開法的基本思想是許多數值方法的基礎。R-K方法不是直接使用Taylor級數,而是利用它的思想2025/6/216龍格-庫塔(R-K)法的基本思想Euler公式可改寫成

則yi+1的表達式與y(xi+1)的Taylor展開式的前兩項完全相同,即局部截斷誤差為O(h2)。Runge-Kutta

方法是一種高精度的單步法,簡稱R-K法2025/6/217同理，改進Euler公式可改寫成

上述兩組公式在形式上共同點:都是用f(x,y)在某些點上值的線性組合得出y(xi+1)的近似值yi+1,

且增加計算的次數f(x,y)的次數,可提高截斷誤差的階。如歐拉法:每步計算一次f(x,y)的值,為一階方法。改進歐拉法需計算兩次f(x,y)的值，為二階方法。局部截斷誤差為O(h3)2025/6/218

于是可考慮用函數f(x,y)在若干點上的函數值的線性組合來構造近似公式，構造時要求近似公式在(xi,yi)處的Taylor展開式與解y(x)在xi處的Taylor展開式的前面幾項重合，從而使近似公式達到所需要的階數。既避免求高階導數,又提高了計算方法精度的階數。或者說,在[xi,xi+1]這一步內多計算幾個點的斜率值，然后將其進行加權平均作為平均斜率，則可構造出更高精度的計算格式，這就是龍格—庫塔（Runge-Kutta）法的基本思想。一般龍格－庫塔方法的形式為2025/6/219其中ai,bij,ci為待定參數，要求上式yi+1在點(xi,yi)處作Tailor展開，通過相同項的系數確定參數。稱為P階龍格－庫塔方法。10Runge-Kutta方法的推導思想對于常微分方程的初值問題的解y=y(x)，在區間[xi,xi+1]上使用微分中值定理，有即2025/6/2111引入記號就可得到相應的Runge-Kutta方法2025/6/2112如下圖即則上式化為即Euler方法Euler方法也稱為一階Runge-Kutta方法2025/6/21二階龍格—庫塔法

在[xi,xi+1]上取兩點xi和xi+a2=xi+a2h,以該兩點處的斜率值K1和K2的加權平均(或稱為線性組合)來求取平均斜率k*的近似值K，即

式中:K1為xi點處的切線斜率值

K1=hf(xi,yi)=hy'(xi)

K2為xi+a2h點處的切線斜率值,比照改進的歐拉法,將xi+a2視為xi+1，即可得

2025/6/2113確定系數c1、c2、a2、b21

，可得到有2階精度的算法格式2025/6/2114因此

將y(xi+1)在x=xi處進行Taylor展開：

將在x=xi處進行Taylor展開：

2025/6/2115K1=hf(xi,yi)2025/6/2116這里有4個未知數，3個方程。存在無窮多個解。所有滿足上式的格式統稱為2階龍格-庫塔格式。令

對應項的系數相等，得到

2025/6/2117注意到，就是二階龍格-庫塔公式，也就是改進的歐拉法。

因此，凡滿足條件式有一簇形如上式的計算格式，這些格式統稱為二階龍格—庫塔格式。因此改進的歐拉格式是眾多的二階龍格—庫塔法中的一種特殊格式。若取，就是另一種形式的二階龍格-庫塔公式。2025/6/2118此計算公式稱為變形的二階龍格—庫塔法。式中為區間的中點。也稱中點公式。

Q:為獲得更高的精度，應該如何進一步推廣？2025/6/2119

二級R-K方法是顯式單步式,每前進一步需要計算兩個函數值。由上面的討論可知,適當選擇四個參數c1,c2,a2,

b21，可使每步計算兩次函數值的二階R-K方法達到二階精度。能否在計算函數值次數不變的情況下,通過選擇不同的參數值,使得二階R-K方法的精度再提高呢?

答案是否定的！無論四個參數怎樣選擇,都不能使公式的局部截斷誤差提高到三階。

這說明每一步計算兩個函數值的二階R-K方法最高階為二階。若要獲得更高階得數值方法,就必須增加計算函數值的次數。三階龍格—庫塔法2025/6/2120為進一步提高精度，在區間[xi,xi+1]上除兩點xi和xi+a2=xi+a2h,以外，再增加一點xi+a3=xi

+a3h

，用這三點處的斜率值K1、K2和K3的加權平均得出平均斜率K*的近似值K，這時計算格式具有形式:

2025/6/2121同理推導二階公式，將y(xi+1)和yi+1在x=xi處進行Taylor展開，使局部截斷誤差達到O(h4)，使對應項的系數相等，得到系數方程組：參數的選擇不唯一，從而構成一類不同的三階R-K公式，下面給出一種常用的三階R-K公式，形似simpson公式：2025/6/21222025/6/2123四階(經典)龍格—庫塔法

如果需要再提高精度，用類似上述的處理方法，只需在區間[xi,xi+1]上用四個點處的斜率加權平均作為平均斜率K*的近似值，構成一系列四階龍格—庫塔公式。具有四階精度，即局部截斷誤差是O(h5)。推導過程與前面類似，由于過程復雜，這里從略，只介紹最常用的一種四階經典龍格—庫塔公式。

2025/6/2124

K1=hf(xi,yi)

K2=hf(xi+a2h,yi+b21K1)

K3=hf(xi+a3h,yi+b31K1+b32K2)

K4=hf(xi+a4h,yi+b41K1+b42K2+b43K3)

其中c1、c2、c3、c4、a2、a3、a4、b21、b31、b32、b41、b42、b43均為待定系數。這里K1、K2、K3、K4為四個不同點上的函數值，分別設其為設yi+1=yi+c1K1+c2K2+c3K3+c4K42025/6/2125

類似于前面的討論，把K2、K3、K4分別在xi點展成h的冪級數，代入線性組合式中，將得到的公式與y(xi+1)在xi點上的泰勒展開式比較，使其兩式右端直到h4的系數相等，經過較復雜的解方程過程便可得到關于ci，ai，bij的一組特解

a2=a3=b21=b32=1/2

b31=b41=b42=0

a4=b43=1

c1=c4=1/6

c2=c3=1/326

四階(經典)Runge-Kutta方法2025/6/2127例1.使用高階R-K方法計算初值問題解:(1)使用三階R-K方法2025/6/2128其余結果如下:(2)如果使用四階R-K方法

ixik1k2k3yi1.00000.10000.10000.11030.12561.11112.00000.20000.12350.13760.15951.24993.00000.30000.15620.17640.20921.42844.00000.40000.20400.23420.28661.66645.00000.50000.27770.32590.41631.99932025/6/2129其余結果如下:

ixik1k2k3k4yi1.00000.10000.10000.11030.11130.12351.11112.00000.20000.12350.13760.13920.15631.25003.00000.30000.15620.17640.17910.20421.42864.00000.40000.20400.23420.23890.27811.66675.00000.50000.27770.32590.33480.40062.00002025/6/21【例】已知一階系統的微分方程為：，初始條件，取仿真步長h=0.1，分別用歐拉法、梯形法和龍格—庫塔法計算該系統仿真第一步的值。解：原方程可變為:即數值積分公式應用

（1）用歐拉法計算根據歐拉公式，將函數表達式及其初始值代入后，可得該系統仿真第一步的值：數值積分公式應用

（2）用梯形法計算：根據預報—校正公式，將函數表達式及其初始值代入后，可得仿真第一步的值。用預報公式求起始值：數值積分公式應用

再用校正公式得到系統仿真第一步的值：數值積分公式應用

二階龍格-庫塔公式（3）用二階龍格—庫塔法計算根據公式先計算出兩個系數，再計算仿真第一步的值：數值積分公式應用

則系統仿真第一步的值為：數值積分公式應用

四階龍格—庫塔（Runge—Kutta）法（4）用四階龍格—庫塔公式計算根據公式先計算出4個系數，再計算仿真第一步的值：數值積分公式應用

數值積分公式應用

則系統仿真第一步的值為：數值積分公式應用

從上述結果可以看出:

對于同一個系統進行仿真計算時，其值的精度是隨著數值積分公式的變化而改變的，其中歐拉法計算精度最低，其次為梯形法和二階龍格—庫塔法，四階龍格—庫塔法計算精度最高。數值積分公式應用

例2：用matlab演示龍格-庫塔法的誤差估計

一個高精度的仿真方法必須將步長控制作為手段。實現步長控制涉及局部誤差估計和步長控制策略兩方面的問題。

龍格-庫塔法的誤差估計RK方法的誤差估計通常是設法找一個低一階的RK公式，將兩個公式計算結果之差作為估計誤差。例如Runge-Kutta-Fehlberg法的計算公式是

RKF1-2公式用另一個一階公式來估計誤差龍格-庫塔法的誤差估計RKM3-4公式：誤差估計式用3階，計算公式為4階。RKM3-4公式誤差估計公式龍格-庫塔法的步長控制

龍格-庫塔法的誤差估計和步長控制策略的基本思想是：每積分一步都設法估計出本步的積分誤差en，然后判斷是否滿足允許誤差E，據此選擇相應的步長控制策略。每一步的局部誤差通常取以下形式en=En/(|yn|+1)其中|yn|是利用誤差估計式計算出的本步的估計誤差。當|yn|較大時，en是相對誤差，當|yn|較小時，en是絕對誤差。這樣作的目的是避免當y的值很小時，en變得過大。仿真模型的運行速度與實際系統運行速度一致，稱為實時仿真。一般方法難以滿足實時仿真的需要：所得模型的執行速度較慢；機理也不符合實時仿真的需要。假設對一般形式的系統進行仿真：以RK-2為例進行分析，其公式為

實時仿真

假定在h/2的時間內計算機剛好計算一次右端函數f,則計算分為兩步：

1在tk時刻利用當前的un、yn計算K1；

2在tn+h/2時刻計算K2,此時un+1無法得到，但實時仿真除了要滿足執行速度的要求外，還要求實時接收外部輸入，并實時得到輸出。?

此種情況下，解決的方法有兩個：對un+1進行預報（增大仿真誤差）或仿真延遲半個計算步距。后者的計算流程如下可見，后種方法的輸出也會延遲半個計算步距，為了克服這個缺陷，人們提出了如下形式的實時二階RK法圖RK-2的計算流程

實時RK-2公式：

其計算流程：假定在h/2的時間內計算機也剛好計算一次右端函數f,則計算也分為兩步：

1在tn時刻利用當前的un、yn計算K1；

2在tn+h/2時刻計算K2,此時un+1/2可以得到，不會引入新的誤差，可實