火電廠仿真課件第5章 面向微分方程的數(shù)值積分法仿真

第五章面向微分方程的數(shù)值積分法仿真數(shù)值積分是數(shù)值分析的一個基本問題。也是復(fù)雜計算問題中的一個基本組成部分。數(shù)值積分往往用極簡單的方法就能較好地得出對所求解的具體數(shù)值問題的解答。但數(shù)值積分的難點在于計算時間有時會過長，有時會出現(xiàn)數(shù)值不穩(wěn)定現(xiàn)象。另外，數(shù)值積分的理論性較強(qiáng)。其理論和方法都已經(jīng)比較成熟，計算精度也比較高。5.1仿真中研究數(shù)值積分法的意義5.2數(shù)值積分法仿真的基本原理5.3歐拉(Euler)法5.4龍格-庫塔（Runge-Kutta）法5.5計算方法的選擇5.6計算步長的選擇5.7面向微分方程的仿真程序構(gòu)成5.1仿真中研究數(shù)值積分法的意義由模擬仿真可知，一個n階連續(xù)系統(tǒng)可以由n個積分器組成的模擬仿真模型予以描述。因此，利用數(shù)字計算機(jī)進(jìn)行連續(xù)系統(tǒng)仿真時，從本質(zhì)上講，就是在數(shù)字計算機(jī)上構(gòu)造出幾個數(shù)字積分器，進(jìn)行幾次數(shù)值積分運算。在系統(tǒng)的仿真中，連續(xù)系統(tǒng)的數(shù)學(xué)模型一般都可用微分方程來表示，但它們通常是高階的微分方程，由前所述，高階的微分方程可化為若干個一階微分方程組成的一個方程組。（5-2）或用狀態(tài)方程表示：（5-1）所以，連續(xù)系統(tǒng)的仿真就是從給定的初始條件出發(fā)，對描述系統(tǒng)動態(tài)特性的常微分方程或常微分方程組進(jìn)行求解，從而得到系統(tǒng)在一定輸入作用下的變化過程。在求解這些微分方程(組)時，最常用、也是最有效的一種方法就是數(shù)值積分法。數(shù)值解的一種近似方法。數(shù)值積分法是求解一系列微分方程：（5-3）5.2數(shù)值積分法仿真的基本原理若對微分方程（5-3）兩端同時取積分，可得當(dāng)時，上式變?yōu)椋海?-4）（5-5）將積分項拆成兩項（5-6）故上式可寫為：（5-7）——此式是方程（5-3）在tn+1時刻的精確解。

在數(shù)值解法中，希望用近似解：代替準(zhǔn)確解，其中：（5-8）為為為的近似值是從常微分方程（5-3）出發(fā)建立的離散數(shù)學(xué)模型——差分方程。上述（5-8）式由此可見，數(shù)值積分法就是在已知微分方程初值的情況下，求解該方程在一系列離散點處的近似值，其特點是步進(jìn)式——根據(jù)初始值逐步遞推地計算出以后各時刻的值。從式（5-8）可知，數(shù)值積分法的主要問題歸結(jié)是如何對f（t，y）進(jìn)行數(shù)值積分——求出f（t，y）在區(qū)間[tn，tn+1]上定積分的近似值Qn。采用不同的方法求Qn，就出現(xiàn)了各種各樣的數(shù)值積分方法。不同的數(shù)值積分將對求解的精度、速度和數(shù)值穩(wěn)定性會產(chǎn)生不同的影響，這將在下述內(nèi)容中具體介紹。在種類繁多的數(shù)值積分法中，在此從實用角度介紹幾種基本的方法5.3歐拉(Euler)法5.3.1簡單歐拉法

歐拉法是一種最簡單的數(shù)值積分法，對于方程：在區(qū)間[t0，

tn+1]上求積分，得到：f(t,y)在[tn,tn+1]區(qū)間上包圍的面積Qn若區(qū)間[tn，tn+1]足夠小，則[tn，tn+1]上的f（t，y）可近似地看成常數(shù)f（tn，yn）

。故可用矩形面積Qn近似代替即：tntn+1f（tn，yn）y(tn)Qn于是有：（5-9）其中：——稱為計算步長或步距將此式寫成差分方程為：（5-10）——著名的歐拉公式（n=1，2，3……）5.3.2改進(jìn)的歐拉法

如果用梯形面積而不是矩形面積來代替每一個小區(qū)間上的曲線積分，就可以提高計算精度，梯形法的計算公式為：（5-11）式中的右端含有待求量yn+1，因而它是隱函數(shù)形式。這種方法不能自行啟動運算，需要依賴其它算法的幫助。每次計算都用歐拉法算出y(t

n+1

)的近似值，以此計算近似值，然后利用梯形公式（5-11）求出修正后的。即有：幫助方法：預(yù)估式校正式（5-12）——改進(jìn)的歐拉公式5.3.3幾個基本概念

簡單的歐拉法是用前一時刻tn的yn求出后一時刻的yn+1，這種方法稱為單步法，它是一種自行啟動的算法。如果求yn+1時需要tn,

tn-1,

tn-2……時刻yn,

yn-1,

yn-2……的值，這種方法為多步法（改進(jìn)的歐拉法為兩步法），它是一種不能自行啟動的算法。1、單步法與多步法簡單的歐拉法表達(dá)式的右端，計算所用的數(shù)據(jù)均已求出，這種公式稱為顯式公式。改進(jìn)的歐拉法表達(dá)式的右端，有待求量，這種公式稱為隱式公式。隱式公式不能自行啟動，需要用預(yù)估-校正法。單步法和顯式在計算上比多步法和隱式方便，但有時為了滿足精度、穩(wěn)定性等方面的要求，需要采用隱式算法。2、顯式與隱式3、截斷誤差這里用泰勒級數(shù)為工具來分析數(shù)值積分公式的精度。假定yn是精確的，用泰勒級數(shù)表示處的精確解，即：顯然，簡單的歐拉法是從以上精確解中取前兩項之和來近似計算，每一步由這種方法引入的誤差稱為局部截斷誤差，簡稱截斷誤差。簡單的歐拉法的截斷誤差為：不同的數(shù)值方法有不同的截斷誤差。一般若截斷誤差為，則方法為r階的。所以方法的階數(shù)可以作為衡量方法精確度的一個重要標(biāo)志。（5-13）由此可見，簡單的歐拉法是1階精度；改進(jìn)的歐拉法由于采用了平均斜率，相當(dāng)于取了泰勒級數(shù)的前3項，因此為2階精度。

分析歐拉法截斷誤差的思想，同樣也適用于其它數(shù)值積分方法。4、舍入誤差由于計算機(jī)的字長有限，數(shù)字不能表示得完全精確，在計算過程中，不可避免地會產(chǎn)生舍入誤差。舍入誤差與計算步長成反比。如果計算步長小，計算次數(shù)多，則舍入誤差就大。舍入誤差除了與計算機(jī)字長有關(guān)以外，還與計算機(jī)所使用的數(shù)字系統(tǒng)、數(shù)的運算次序以及計算所用的子程序的精度等因素有關(guān)。5、數(shù)值解的穩(wěn)定性問題采用數(shù)值積分法求解穩(wěn)定的常微分方程，應(yīng)該保持原系統(tǒng)的穩(wěn)定特征。但是：（1）在計算機(jī)逐步計算時，初始數(shù)據(jù)的誤差及計算過程的舍入誤差對后面的計算結(jié)果將產(chǎn)生影響。（2）如果計算步長取的不合格，有可能使仿真出現(xiàn)不穩(wěn)定的結(jié)果。下面我們簡單討論一下這個問題。差分方程的解與微分方程的解類似，可分為特解和通解兩部分。與穩(wěn)定性有關(guān)的是方程的通解，它取決于差分方程的特征根是否滿足穩(wěn)定性條件。顯然，要使該差分方程是穩(wěn)定的，必須使下式成立。即：表明：為使數(shù)字仿真穩(wěn)定，對計算步長應(yīng)有所限制。另外，穩(wěn)定性還與系統(tǒng)的特性以及數(shù)值積分方法有關(guān)。上述分析歐拉法穩(wěn)定性的思想，同樣適用于其它數(shù)值積分方法。（5-14）（5-15）例如：為了考查歐拉法的穩(wěn)定性，我們研究檢驗方程（TestEquation）：其中，為方程的特征根，對此有：5.4龍格-庫塔（Runge-Kutta）法由前面的分析可知，將泰勒展開式多取幾項以后截斷，就能提高截斷誤差的階數(shù)和計算精度。然而，直接采用泰勒展開方法要計算函數(shù)的高階導(dǎo)數(shù)，運用起來不便。龍格-庫塔方法的基本思想是：用幾個點上的函數(shù)值的線性組合代替函數(shù)的各階導(dǎo)數(shù)，然后按泰勒級數(shù)展開確定其中的系數(shù)，這樣既可避免計算高階導(dǎo)數(shù)，又可提高積分的精度。龍格-庫塔法有多種形式，以下從實用角度直接給出公式的形式和相應(yīng)的精度。5.4.1龍格-庫塔方法的基本思想5.4.2二階龍格-庫塔方法2階龍格-庫塔方法的公式為：（5-16）上式表示的數(shù)值解是用的泰勒級數(shù)在2階導(dǎo)數(shù)以后截斷所求得的，因此稱為2階方法。故2階龍格-庫塔法與式改進(jìn)的歐拉法相比，實質(zhì)完全相同。所以改進(jìn)的歐拉法實質(zhì)上是2階龍格-庫塔法。截斷誤差為：實時仿真的2階龍格-庫塔方法:（5-17）其截斷誤差為：5.4.3四階龍格-庫塔方法4階龍格-庫塔法是一種最常用的方法。其經(jīng)典表達(dá)式為：（5-18）其截斷誤差為：顯然，4階龍格-庫塔法的計算量較大，但計算精度較高，在比較不同算法的計算精度時，常以它的計算結(jié)果作為標(biāo)準(zhǔn)。實時仿真的4階龍格-庫塔方法（5-19）以上公式都是標(biāo)量形式，如果要換成向量形式，只要把式中的標(biāo)量y，f，k換成向量Y，F(xiàn)，K即可。從理論上講，可以構(gòu)造任意階數(shù)的龍格-庫塔方法，但是，精度的階數(shù)與計算函數(shù)值f

的次數(shù)之間并非等量增加的關(guān)系，見下表所列：由此可見，4階經(jīng)典龍格-庫塔方法有其優(yōu)越性，而4階以上的龍格-庫塔方法計算f所需的次數(shù)比階數(shù)多，增加了計算量，從而限制了更高龍格-庫塔方法的應(yīng)用。對于大量的實際問題，4階方法已可滿足精度要求，所以得到了廣泛的應(yīng)用。我們?nèi)圆捎脵z驗方程進(jìn)行討論，對它利用泰勒級數(shù)展開得：5.4.4龍格-庫塔公式的穩(wěn)定區(qū)域?qū)τ冢?，將它代入上式得：?-20）（5-21）令：——將其代入上式得到該式的穩(wěn)定條件為：（5-22）由此穩(wěn)定條件，下表給出了各階龍格-庫塔公式的穩(wěn)定區(qū)域。表5.2龍格-庫塔公式的穩(wěn)定區(qū)域在使用龍格-庫塔公式時，選取的步長應(yīng)使落在穩(wěn)定區(qū)域內(nèi)。否則，在計算時會產(chǎn)生很大的誤差，從而得不到穩(wěn)定的數(shù)值解。例如：用4階龍格-庫塔方法解：（1）用解析法求解：本例是穩(wěn)定的；（2）用數(shù)值法求解：當(dāng)h=0.1時，數(shù)值解是穩(wěn)定的;當(dāng)h=0.2時，數(shù)值解就不穩(wěn)定了，這時因為：此數(shù)值在穩(wěn)定區(qū)間以外，所以數(shù)值解不能收斂。這種對步長有限制的數(shù)值積分法稱為條件穩(wěn)定積分法。從4階龍格-庫塔方法的穩(wěn)定條件：中可以看出，系統(tǒng)的特征根越大，需要的積分步長就越小，這一點可作為選擇步長的依據(jù)。步長的大小，除了與數(shù)值積分方法的階數(shù)有關(guān)外，還與方程本身的性質(zhì)有關(guān)。除以上介紹的歐拉法、龍格-庫塔法外，數(shù)值積分方法還有許多種，如亞當(dāng)姆斯方法、吉爾方法等等。此處不一一介紹。5.4.5四階龍格-庫塔法仿真程序設(shè)計目前，應(yīng)用于系統(tǒng)仿真的商用軟件包諸多，但無論是開發(fā)者還是應(yīng)用者，了解程序的設(shè)計思想、分析程序的結(jié)構(gòu)與功能，完善和編寫仿真程序都是必要的。因此，此處通過對經(jīng)典龍格-庫塔方法（即式（5-18））仿真程序的介紹，使讀者了解和掌握編寫仿真程序的基本技術(shù)。為了便于說明程序的實質(zhì)，以下用具體的例子來討論。假設(shè)有一個2階系統(tǒng)：令：則：以下是用C語言編寫的原理性程序：

程序中變量說明：i，n: 分別為循環(huán)控制和積分器個數(shù)變量。h，t，tmax：分別為步長、時間和仿真總時間變量。yn[10]:

存放tn時刻y值的臨時變量數(shù)組。y[10]:

存放y值的變量數(shù)組。doty[10]:

存放的數(shù)組。k1[10]，k2[10]，k3[10]，k4[10]:分別為存放K1，K2，K3，K4的數(shù)組。程序清單：（“/*”和“*

/”之間的文字是非執(zhí)行語句，作為注釋使用。）#include<stdio.h>voiddiffeq();inti,n;floath,tmax,t,yn[10],y[10],doty[10],k1[10],k2[10],k3[10],k4[10];main() /*主函數(shù)*/{n=2; /*設(shè)置積分器個數(shù)*/y[1]=0;y[2]=0 /*設(shè)置積分器初值*/

h=0.01;tmax=20;t=0/*設(shè)置步長、仿真總時間和初始時間*/

whilt(t<=tmax) /*判斷是否仿真時間到*/{printf(“t=%f”,t)； /*打印時間*/

printf(“y[1]=%fy[2]=%f\r\n”,y[1],y[2]);/*打印積分器輸出*/for(i=1,i<=n,i++)

yn[i]=y[i]; /*保留tn時刻的初值*/dif(); /*計算*/for(i=1;i<=n;i++){k1[i]=doty[i]; /*計算K1*/y[i]=yn[i]+(h/2)*k1[i];}t=t+h/2dif(); /*計算*/for(i=1;i<=n;i++){k2[i]=doty[i]; /*計算K2*/y[i]=yn[i]+(h/2)*k2[i];}dif(); /*計算*/for(i=1;i<=n;i++){k3[i]=doty[i]; /*計算K3*/y[i]=yn[i]+(h/2)*k3[i];}t=t+h/2dif(); /*計算*/for(i=1;i<=n;i++){k4[i]=doty[i]; /*計算K4*/y[i]=yn[i]+(h/6.0)*(k1[i]+2*k2[i]+2*k3[i]+k4[i]);

/*計算tn+1時刻的y值*/}}}voiddif() /*計算變量導(dǎo)數(shù)的函數(shù)*/{floatu;u=1; /*計算u=1(t)*/

doty[1]=y[2];

doty[2]=4*u-2.828*y[2]-4*y[1];}下面是程序的主函數(shù)main()的處理流程圖

編寫以下實時仿真的4階龍格-庫塔法的仿真程序,對以下2階系統(tǒng)進(jìn)行數(shù)字仿真:練習(xí):并通過仿真實驗,驗證使該算法穩(wěn)定的步長h取值區(qū)間.5.5計算方法的選擇數(shù)值積分方法很多，在實際使用時存在一個選擇問題。對于一個具體問題，如何選擇具有一定的難度，至今尚無一種確定的方法。一般來說，數(shù)值積分方法的選擇應(yīng)考慮的因素有：1、精度要求數(shù)值積分法的精度受截斷誤差舍入誤差積累誤差的影響，而這些誤差與積分方法、步長、計算時間、計算機(jī)精度等有關(guān)。一般：（1）積分方法階數(shù)越高，截斷誤差越小，精度越高；（2）步長越小，精度越高；（3）多步法的精度高于單步法；（4）隱式算法的精度高于顯式算法。因此，當(dāng)需要高精度時，可采用高階、多步、較小步長、隱式算法。但是，步長的減小往往會增加迭代次數(shù)并增大舍入誤差和積累誤差。與此相反，在精度要求不高時，最好使用低階方法。計算速度取決于計算步數(shù)、每步積分所需的時間。而每步的計算時間又與積分方法有關(guān)，它主要取決于計算導(dǎo)數(shù)的次數(shù)（4階龍格-庫塔方法每步要計算4次導(dǎo)數(shù)），導(dǎo)數(shù)的計算是最費時的；為了加快計算速度，在積分方法已定的條件下，應(yīng)要保證精度的前提