連續(xù)系統(tǒng)的計(jì)算機(jī)模擬

1、第2章連續(xù)系統(tǒng)的計(jì)算機(jī)模擬本章討論連續(xù)系統(tǒng)的模擬技術(shù)，由于這類系統(tǒng)中狀態(tài)隨時(shí)間連續(xù)動態(tài)地變化，常常具有一定的規(guī)律，故可用一些數(shù)學(xué)方程來描述，這些方程就是系統(tǒng)的數(shù)學(xué)模型，通常以微分方程、代數(shù)方程為多見。下面將介紹利用數(shù)值積分法對連續(xù)系統(tǒng)進(jìn)行數(shù)字模擬的基本原理和具體方法，并給出數(shù)值積分法中幾個(gè)常用的算法以及實(shí)現(xiàn)這些算法的計(jì)算程序，最后介紹兩個(gè)建模實(shí)例。數(shù)值積分法不僅方法種類多，而且有較強(qiáng)的理論性，本章由淺入深地介紹幾種常見的數(shù)值解法。主要為單步法中的四階龍格-庫塔法與默森法和多步法中的亞當(dāng)斯法。使用數(shù)字計(jì)算機(jī)對連續(xù)系統(tǒng)進(jìn)行模擬，首先必須將連續(xù)系統(tǒng)離散化，并將它轉(zhuǎn)化為差分方程，以建立所謂的模擬系統(tǒng)的

2、數(shù)學(xué)模型。描述連續(xù)系統(tǒng)動態(tài)特征的數(shù)學(xué)模型是多種多樣的，除微分方程外，還有傳遞函數(shù)、結(jié)構(gòu)圖及狀態(tài)方程等，由于篇幅所限，本書不討論后兩種方法。建立系統(tǒng)的模擬模型之后，就要選擇計(jì)算機(jī)語言（也叫算法語言）編寫系統(tǒng)模擬程序，在計(jì)算機(jī)上運(yùn)行，將結(jié)果保留在數(shù)據(jù)文件中以待傳輸和處理。由于模擬的目的不同，可以選用不同的模擬模型和算法，其特點(diǎn)是運(yùn)算精度高，對于不同的計(jì)算機(jī)，字長一般在16位-72位之間，也可采用浮點(diǎn)運(yùn)算和雙精度運(yùn)算，其精度一般可達(dá)千萬分之一到百萬分之一。當(dāng)然結(jié)果的精度與所選的算法有關(guān)。這可以根據(jù)實(shí)際需要選擇機(jī)器、算法和模擬的步長。數(shù)字計(jì)算機(jī)儲存容量大，可進(jìn)行各種運(yùn)算，在以前認(rèn)為是不可能解決的問題，

3、利用數(shù)字計(jì)算機(jī)都可容易地或有可能得到解決。本章介紹的方法適應(yīng)性較強(qiáng)，應(yīng)用也十分廣泛。數(shù)字計(jì)算機(jī)上還有各種功能的軟件包（即一些子程序），用戶可以稍加修改或不經(jīng)修改就可以用于自己的模擬程序中，解決自己的實(shí)際問題，使用非常方便。2.1歐拉（Euler）法在討論連續(xù)系統(tǒng)的計(jì)算機(jī)模擬之前，讓我們先看一個(gè)化學(xué)反應(yīng)的例子，通過這個(gè)例子我們可以看到怎樣使用數(shù)字計(jì)算機(jī)模擬一個(gè)實(shí)際問題，雖然介紹的是歐拉法，但是分析問題的思路同樣適用與其它數(shù)值積分法。當(dāng)兩種物質(zhì)A和B放到一起產(chǎn)生化學(xué)反應(yīng)時(shí)，產(chǎn)生第三種物質(zhì)C,一般一克A與一克B結(jié)合產(chǎn)生2克的C物質(zhì)，形成C的速率與A和B的數(shù)量乘積成正比，同樣C也可分解為A和B,C的分

4、解速率正比與C的數(shù)量，即在任何時(shí)刻，如果a,b,和c是化學(xué)物質(zhì)A,B,和C的數(shù)量，即在任何時(shí)刻，如果a,b,和c是化學(xué)物質(zhì)A,B,C的數(shù)量，它們的增加和減少的速度服從下列微分方程。dada=KC-Kabdt21db=K2c-Kabdt21dc_(2.1 )dc=2Kab-2KCdt12其中K1和K2是比例常數(shù)（一般而言這些比例常數(shù)會隨溫度和壓力發(fā)生變化，但在模擬過程中，為了簡化模型，一般不允許其變化，故一律視為常數(shù)）。在給出常數(shù)K1和K2值以及A和B的數(shù)量（C=0）后，我們希望能確定有多少C物質(zhì)產(chǎn)在出來，這種化學(xué)反應(yīng)速率的決定在化學(xué)工業(yè)上是有意義的。模擬該系統(tǒng)的一個(gè)直接的方法是在t=0時(shí)開始，

5、使t以At間隔增加。假定化學(xué)量在At時(shí)間步長內(nèi)不變，而只能在At結(jié)束的瞬間發(fā)生變化，這樣在每個(gè)At結(jié)束日的A(或B或C)的數(shù)量就可以從At開始時(shí)的數(shù)值由下式求出a(t：t)=a(t)da.t(2同樣的方程b(t+At)和C(t+At),也可寫出。假定模擬周期為T,可將T分成N個(gè)小的時(shí)間步長At,及T=N-At.2)K1和K2值在時(shí)間為零時(shí)，我們知道a(0)、b(0)和c(0)的初始值，從這些初始值及常數(shù)出發(fā)，就可以計(jì)算出a(b(c(使用這些值,a(2b(2c(2At時(shí)間內(nèi)的化學(xué)量的變化值。At)=a(0)+K2?C(0)-Kia(0)-b(0)-At)=b(0)+K2C(0)-Kia(0)-b

6、(0)-At)=c(0)+2Ki-a(0)-b(0)-2K2c(0)又可計(jì)算系統(tǒng)的下一個(gè)狀態(tài)，即At)=a(At)=b(At)=c(AtAtAt2At時(shí)的狀態(tài)。At)+K2?c(At)-Kia(At)b(At)At)+K2c(At)-Kia(At)-b(At)At)+2Kia(At)-b(At)-2K2c(AtAtAt)At2.1表不圖2.i化學(xué)反應(yīng)例子模擬程序框圖C語言編寫，程序中的初值如下：使用2At時(shí)的系統(tǒng)狀態(tài)，又可寫出3At時(shí)的狀態(tài)，依此類推，以At為間隔，進(jìn)彳TN步計(jì)算，就可寫出周期T的系統(tǒng)狀態(tài)得到所期望的結(jié)果，這個(gè)過程可用圖模擬程序使用Ki=0.008/g.minC(0)=0,T=

7、5mins,K20.002/min,a(0)=100g,/t=0.1min,b(0)=50gN=50其程序清單如下：#include<stdio.h>#include<string.h>floatk1,k2;staticfloatA53,B53,C53,delta,t;voidstrut(int);main()inti;A1=100.0;B1=50.0;C1=0.0;t=0;delta=0.1;k1=0.008;k2=0.002;for(i=1;i<53;i+)printf("%2d",i);printf("%10.2f,%10.2f

8、,%10.2f,%10.2fn",t,Ai,Bi,Ci);strut(i);return;voidstrut(inti)Ai+1=Ai+(k2*Ci-k1*Ai*Bi)*delta;Bi+1=Bi+(k2*Ci-k1*Ai*Bi)*delta;Ci+1=Ci+2.0*(k1*Ai*Bi-k2*Ci)*delta;t=t+delta;運(yùn)行此程序，輸出模擬結(jié)果如下：10.00,100.00,50.00,0.0020.10,96.00,46.00,8.0030.20,92.47,42.47,15.0640.30,89.33,39.33,21.3450.40,86.52,36.52,26.9

9、560.50,84.00,34.00,32.0070.60,81.72,31.72,36.5580.70,79.66,29.66,40.6990.80,77.77,27.77,44.45100.90,76.05,26.05,47.89111.00,74.48,24.48,51.04121.10,73.03,23.03,53.94131.20,71.70,21.70,56.61141.30,70.46,20.46,59.07151.40,69.32,19.32,61.36161.50,68.26,18.26,63.48171.60,67.28,17.28,65.44181.70,66.36,16

10、.36,67.28191.80,65.51,15.51,68.99201.90,64.71,14.71,70.59212.00,63.96,13.96,72.08222.10,63.26,13.26,73.48232.20,62.60,12.60,74.79242.30,61.99,11.99,76.03252.40,61.41,11.41,77.18262.50,60.86,10.86,78.27272.60,60.35,10.35,79.30282.70,59.87,9.87,80.27292.80,59.41,9.41,81.18302.90,58.98,8.98,82.04313.00

11、,58.57,8.57,82.86323.10,58.19,8.19,83.63333.20,57.82,7.82,84.36343.30,57.48,7.48,85.05353.40,57.15,7.15,85.70363.50,56.84,6.84,86.32373.60,56.55,6.55,86.91383.70,56.27,6.27,87.46393.80,56.00,6.00,87.99403.90,55.75,5.75,88.50414.00,55.51,5.51,88.97424.10,55.29,5.29,89.43434.20,55.07,5.07,89.86444.30,

12、54.86,4.86,90.27454.40,54.67,4.67,90.66464.50,54.48,4.48,91.03474.60,54.31,4.31,91.39484.70,54.14,4.14,91.73494.80,53.98,3.98,92.05504.90,53.82,3.82,92.35515.00,53.68,3.68,92.65525.10,53.54,3.54,92.93從這個(gè)例子中我們以上例子的算法就是數(shù)值積分法中最簡單的一種方法叫作歐拉法可以看出使用計(jì)算機(jī)解題的過程，由于歐拉法本身的特點(diǎn)，決定了其計(jì)算精度差，現(xiàn)在幾乎無人在實(shí)際工作中使用，但它導(dǎo)出簡單，能說明構(gòu)造數(shù)

13、值解法一般計(jì)算公式的基本思想，模擬程序也清晰易懂，故人們樂意用它做為構(gòu)造數(shù)值解法的入門例子。其一般解法如下：設(shè)給定微分方程：v=f(t,y)y=y0(2.3)在區(qū)間(tn,tn+1)上求積分,得y(tn+1)=y(tn)+/f(t,y)dt(2.4)如積分間隔足夠小，使得tn與tn+1之間的f(t,y)可近似的看成常數(shù)f(tn,yn),就可以用矩形面積近似地代替在該區(qū)間上的曲線積分，于是在tn+1時(shí)的積分值為y(tn1)ynf(tn，yn)h=yn.i(2.5)將上式寫成以下差分方程形式：yn+=yn+hfnn=1,2,3，(2.6)這就是歐拉公式。它是一個(gè)遞推的差分方程，任何一個(gè)新的數(shù)值解y

14、n+1均是基于前一個(gè)數(shù)值解以及它的導(dǎo)數(shù)f(tn,yn)值求得的。只要名定初始條件yo及步長h就可根據(jù)f(t0,y0)算出yi的值，再以yi算出y2,如此逐步算出y3,y4,，直到滿足所需計(jì)算的范圍才停止計(jì)算。歐拉法的基本思路是把連續(xù)的時(shí)間t分割成等間隔的At,在這些離散的時(shí)刻算得函數(shù)值，根據(jù)這些值在函數(shù)圖上可得到一條折線，所以歐拉法又叫折線法，其特點(diǎn)是分析方法簡單，計(jì)算量小，但計(jì)算精度低(后面將討論歐拉法與其它方法的比較)。下圖為歐拉折線法的幾何意義。Y f(t圖2.2歐拉法的幾何意義如果用梯形面積來代替一個(gè)小區(qū)間的曲線積分，就可克服以小矩形計(jì)算的缺點(diǎn)，提高精度，梯形法計(jì)算公式為hyn1-yn

15、2f(tn，yn)f(tn1,yn1)=yn2(fnfni)-,汴(2(2.7)上式為隱式公式，因?yàn)楣接叶撕衴n+1,這是未知的待求量，故梯形法不能自行啟動運(yùn)算,而要依賴于其它算法的幫助，比如說用歐拉公式法求出初值，算出y(yn中)的近似值yF書，然后將其帶入微分方程，計(jì)算fn書的近似值fnl=f(tn書，yP書)，最后利用梯形公式反復(fù)迭代。如迭代一次后就認(rèn)為求得了近似解，這就是改進(jìn)的歐拉法，其公式為預(yù)估yP 書=yn +hf(tn,yn)(2-8)校正yC書=yn+2f(tn.Yn)+fP(tn+.yP+)(2-9)h上式中第一個(gè)為預(yù)估公式，第二個(gè)為校正公式。通常這種方法稱為預(yù)估矯正法。

16、在校正公式中計(jì)算了兩點(diǎn)的斜率，再求其平均值，故計(jì)算量比歐拉法要大些。2. 3數(shù)值積分的幾個(gè)基本概念1 .單步法與多步法數(shù)值積分法都用遞推公式求解，而遞推公式是步進(jìn)式的，當(dāng)由前一時(shí)刻的數(shù)值yn就能計(jì)算出后一時(shí)刻的數(shù)值yn+1時(shí)，這種方法稱為單步法，它是一種能自啟動的算法，歐拉法是單步法.反之，如果求yn+1時(shí)要用到y(tǒng)n,yn-1,yn-2,等多個(gè)值，則稱為多步法，由于多步法在一開始時(shí)要用其它方法計(jì)算該時(shí)刻前面的函數(shù)值，所以它不能自啟動，以上所講的改進(jìn)的歐拉法就是一個(gè)多步法的例子。2 顯式與隱式在遞推公式中，計(jì)算yn+1時(shí)所用到的數(shù)據(jù)均已算出的計(jì)算公式稱為顯示公式.相反，在算式中隱含有末知量yn+

17、1的則稱為隱含公式.這需用另一個(gè)顯示公式估計(jì)一個(gè)值，再用隱式公式進(jìn)行迭代運(yùn)算，這就是預(yù)估-校正法，這種方法不能自啟動.用單步法與顯示公式在計(jì)算上比用多步法和隱式公式方便.但有時(shí)要考慮精度與穩(wěn)定性時(shí)，則必須采用隱式公式運(yùn)算.3 截?cái)嗾`差一般使用臺勞級數(shù)來分析數(shù)值積分公式的精度.為簡化分析，假定前一步彳#的結(jié)果yn是準(zhǔn)確的，則用臺勞級數(shù)求得tn+1處的精確解為1一1(2.10)ylnkyghg)”，(«)Phn與前面的遞推公式比較，歐拉法是從以上精確解中只取前兩項(xiàng)之和來近似計(jì)算yn+1的，由這種方法單獨(dú)一步引進(jìn)的附加誤差通常稱著局部截?cái)嗾`差，它是該方法給出的值與微分方程解之間的差值，故又

18、稱為局部離散誤差。歐拉法的局部截?cái)嗾`差為O(h2),則方法稱為r階的,y(tn1)-y(tn)=O(h2)不同的數(shù)值解法，其局部截?cái)嗾`差也不同。一般若截?cái)嗾`差為所以方法的階數(shù)可以作為衡量方法精確度的一個(gè)重要標(biāo)志。歐拉法是一階精度的數(shù)值解法，而改進(jìn)的歐拉法(梯形法)相當(dāng)于取臺勞級數(shù)的前三項(xiàng)，故其解為二階精度。4 舍入誤差由于計(jì)算機(jī)的字長有限，數(shù)字不能表示彳#完全精確，在計(jì)算過程中不可避免地會產(chǎn)生舍入誤差.舍入誤差與步長h成反比，如計(jì)算步長小，計(jì)算次數(shù)多則舍入誤差大.產(chǎn)生舍入誤差的因素較多，除了與計(jì)算機(jī)字長有關(guān)外，還與計(jì)算機(jī)所使用的數(shù)字系統(tǒng)，數(shù)的運(yùn)算次序及計(jì)算f(t,y)所用的程序的精確度等因素有

19、關(guān).5 數(shù)值解的穩(wěn)定性以上對連續(xù)系統(tǒng)的模擬用到了差分方程，把初始值代入遞推公式進(jìn)行迭代運(yùn)算，如果原系統(tǒng)是穩(wěn)定的，由數(shù)值積分法求得的解也應(yīng)該是穩(wěn)定的.但是，在計(jì)算機(jī)逐次計(jì)算時(shí)，初始數(shù)據(jù)的誤差及計(jì)算過程中的舍入誤差對后面的計(jì)算結(jié)果會產(chǎn)生影響.而且計(jì)算步長選擇得不合理時(shí)，有可能使模擬結(jié)果不穩(wěn)定.以下我們簡要地討論這一問題.差分方程的解與微分方程的解類似，均可分為特解與通解兩部分.與穩(wěn)定性有關(guān)的為通解部分，它取決于該差分方程的特征根是否滿足穩(wěn)定性要求。例如，為了考慮歐拉法的穩(wěn)定性，可用檢驗(yàn)方程y=Xyo其中入為方程的特征根。對此方程有yn1=yn'ynh=(11h)yn要使該微分方程穩(wěn)定，須使

20、下式成立|1+入h|<1有此可得|入h|<2或h<2T要保證用歐拉法進(jìn)行的計(jì)算是穩(wěn)定的，積分步長h必須小于系統(tǒng)時(shí)間常數(shù)的兩倍.所以積分步長的選擇就要慎重，不能隨意取任何值。2.4龍格貝塔(Runge-Kutta)法仍以(2.3)方程為例.已知y(t0)=y0,假設(shè)我t0開始以h增長，ti=t0+h,t1時(shí)刻為2一一一»yi=y(M+h),在b附近展開成臺老級數(shù)，保留h項(xiàng)，則有.h3fdy；:f2yi=y()f(t0,y0)h()h(2-ii)2二ydt二t上式括號內(nèi)下標(biāo)為0,表示括號中的函數(shù)用t=to,y=yo代入，以下均同。假設(shè)上式的解可寫成如下形式y(tǒng)i=yo(b

21、iKib2K2)h(2.i2)其中,Ki=f(10,yo)K2=f(t0C2h,y°aiKih)對K2式右端得函數(shù)在t=t0,y=y。處展開成臺勞級數(shù)，保留h項(xiàng)，可得到：-汗汗K2：f(t。y。)(C2aiKi-)0h二t二y將Ki與K2代入(2-i2)式得:f;:fy1=y°,b)hf(to,y°),bzhf(to,y°),(C2aiKi)ohftFy將上式與(2-12)式比較，可得系數(shù)ai,bi,C2,b2方程如下:bi+b2=1,b2c2=二2ib2ai二二、2以上三個(gè)方程，四個(gè)未知數(shù)，因此有無窮個(gè)解，若限定bi=b2則可得一個(gè)解：ai=C2=i將

22、它們代入(2-i2)式可得出一組計(jì)算公式hyi=y。(KJ&)-其中Ki=f(t0,y0),k2=f(t0+h,y0+kih),寫成一般的遞推形式如下：hy(tni)：yni二yn2(KiK2)(2.i3)其中，Ki=f(tn,yn),K2=f(tn+h,yn十K1h),式中只取h,h2兩項(xiàng)，而將h3以上的高階項(xiàng)略去，所以此遞推公式的截?cái)嗾`差正比于h的三次方，計(jì)算過程中只取h和h2兩項(xiàng)，故此方法被稱為二階龍格-庫塔法。一般使用的為四階龍格-庫塔公式，在展開臺勞級數(shù)時(shí)保留h,h2,h3,和h4項(xiàng),其截?cái)嗾`差為h5,四階龍格庫塔公式為y(tni)：yni=ynh(Ki2K22K3K4)(2

23、.i4)6式中：Ki=f(tn,yn)K2=f(tn+、+Ki)rhhK3=f(tn-,yn-K2)K4=f(tnh.ynhK3)四階龍格一庫塔法公式精度較高，故其應(yīng)用較為普遍。下面再給出幾組系數(shù)值設(shè)a1=1C2C2C21212b1 =0bibi1412b2b2b2二1二1讀者可以寫出其相應(yīng)的四階龍格一庫塔法公式。龍格一庫塔法的特點(diǎn)如下:1 .龍格庫塔法有多組a1,C2b1b2值，故可有多種龍格庫塔公式，常使用的有:yn書=yn+K?hK=f(tn,yn)K2=f(tn+,yn+K1h)所有龍格-庫塔公式都有以下幾個(gè)特點(diǎn)：(1)在計(jì)算yn+1時(shí)只用到y(tǒng)n,而不直接用yn-1,yn-2等項(xiàng)，也就

24、是計(jì)算后一步時(shí)只用到前一步的計(jì)算結(jié)果，故稱為單步法，單步法可以自啟動，且使用的存儲量也小。(2)在計(jì)算過程中，步長h可以改變，這要由計(jì)算精度來定，但是在一步中，為計(jì)算若干個(gè)系數(shù)Ki(稱成為龍格-庫塔系數(shù))，必須使用同一個(gè)步長ho(3).不論是幾階龍格-庫塔公式，計(jì)算式中均為兩部分組成。第一部分為上一步的結(jié)果yn,第二部分為h=tn-tn-1中對他y)的積分。它是步長h乘上各點(diǎn)斜率的加權(quán)平均。例如上面計(jì)算的四階公式中，取四點(diǎn)的斜率：K1、K2、K3和K4,然后對K2和K3取兩份，K1和K4各取一份，進(jìn)行加權(quán)平均。3 .龍格一庫塔法的精度取決于步長h的大小及求解的方法。實(shí)踐證明：為達(dá)到同樣的精度，

25、四階方法的步長h可以比二階方法的h大10倍，而四階方法的每步計(jì)算量僅比二階方法大一倍，所以總的計(jì)算量仍比二階方法小，使用的CPU機(jī)時(shí)也少。一般工程中四階龍格庫塔公式已能達(dá)到要求的精度，故不再使用更高階的公式。下表示出了f的計(jì)算次數(shù)與精度階數(shù)的關(guān)系。每一步計(jì)算f的次數(shù)234567n=>8精度階數(shù)234456n-2可見精度的階數(shù)與計(jì)算函數(shù)值f的次數(shù)之間的關(guān)系并非等量增加。4 .若在展開成臺勞級數(shù)時(shí)，只取h這一項(xiàng)，而將h2以及h2以上的高次項(xiàng)略去，可得yn1=yn,hf(tn,Yn)這是歐拉公式，因此歐拉公式也可看成時(shí)一階龍格-庫塔公式，她的截?cái)嗾`差正比h2,其精度最低。如果將h取得很小，能否

26、用歐拉公式得到很高的精度呢？從理論上講是可以的，但是實(shí)際上由于計(jì)算機(jī)字長有限，在計(jì)算中有舍入誤差。步長h越小，計(jì)算的步數(shù)越多，舍入誤差的積累會越來越大，故用歐拉法時(shí)不能用很小的ho2.5四階龍格一庫塔法模擬程序及應(yīng)用目前已有多種四階龍格一庫塔法模擬(仿真)軟件包推出，雖然子程序稍有不同，但總的結(jié)構(gòu)還是相同的。對于連續(xù)系統(tǒng)的龍格一庫塔法的計(jì)算機(jī)程序，其大致結(jié)構(gòu)如下圖所示：圖2.3龍格一庫塔法程序簡要框圖圖2.3中每一個(gè)程序模塊的功能為：1 .主函數(shù)：主函數(shù)實(shí)現(xiàn)對整個(gè)模擬系統(tǒng)的控制，這是通過調(diào)用各個(gè)子程序?qū)崿F(xiàn)的。2 .輸入及初始化函數(shù)：主要對系統(tǒng)參數(shù)輸入初值，例如模擬時(shí)間，積分步長，方程階數(shù)等。這

27、可以從鍵盤輸入，也可從數(shù)據(jù)文件輸入。3 .運(yùn)行模塊：在運(yùn)行子程序中，將反復(fù)調(diào)用數(shù)值積分和方程右函數(shù)模塊，進(jìn)行迭代計(jì)算,可以每計(jì)算一步便顯示一次結(jié)果，也可以計(jì)算數(shù)步顯示一次結(jié)果，屏幕的顯示使用戶可以隨時(shí)監(jiān)視系統(tǒng)運(yùn)行的進(jìn)程，以便人工干預(yù)。4 .輸出子程序：按要求輸出打印結(jié)果，可以打印，也可以繪圖，視需要而定。四階龍格一庫塔法公式前面已經(jīng)給出，現(xiàn)在再將它寫成可編程的向量形式1yy(Kfko2K°KJyi,n1yi,ni1i2i3i46Kilhfi(tn,y1n,y2n,yNn)Ki2= hfi(tnh2，y1n112K22, ,yNn 2 Kn2)Ki3rh111=hfi(tn2,yn2K

28、12,丫2n萬心,，yNn2KN2)Ki4=hfi(tnh)1n*13/2n七3,yNn*N3)式中，i=1,2,3,.N,N為方程階數(shù)。令DT=h,Ki1，Ki2，Ki3和Ki4分別換成D,D(I),D(3)(I)和口(4),D(I)為調(diào)微分子程序求得，這些分別為變量的導(dǎo)數(shù)，這樣上式變?yōu)椋篩(I)=YS(I)+0.166667*2.0*(RK1(I)+RK3(I)+4.0*RK2(I)+D(I)*DT其中RK1(I)=D(I)*0.5DTRK2(I)=D(I)*0.5DTRK3(I)=D(I)*DT代入上式：Y(I)=YS(I)+0.166667*DT*D(I)+2.0*D(I)+2.0*D

29、(I)+D(I)即yn+1=yn+h(K1+2K2+2K3+K4)設(shè)有一個(gè)微分方程：y(t)+P1y(t)+y(t)=1.0令y1(t),y2(t)為兩各個(gè)狀態(tài)變量，且y1(t)=y(t),y1(t)=y2(t)代入原方程得：y2(t)=-y1(t)-P1y1(t)+1.0其中P1為常數(shù)，且P1=0.1初始條件：y1(0)=0,y2(0)=0模擬參數(shù)初值：模擬時(shí)間為50.0積分步長為0.1打印點(diǎn)數(shù)為101程序清單如下：主程序：/*p23example*/#include"stdio.h"#include"conio.h"#include"mat

30、h.h"#defineNORDER3#defineNPARAM2floatyNORDER,dNORDER,pNPARAM,t;voidoutput(void)intj;printf("%9s","time");for(j=0;j<NORDER;j+)printf("y%d",j);putchar('n');voiddifq(void)d0=-p0*y0*y1;d1=p0*y0*y1-p1*y1;d2=p1*y1;voidrun(floatdt)inti;floatrk3NORDER,ysNORDER;d

31、ifq();for(i=0;i<NORDER;i+)rk0i=di*0.5*dt;ysi=yi;yi=ysi+rk0i;t+=0.5*dt;difq();for(i=0;i<NORDER;i+)rk1i=di*0.5*dt;yi=ysi+rk1i;difq();for(i=0;i<NORDER;i+)rk2i=di*dt;yi=ysi+rk2i;t+=0.5*dt;difq();for(i=0;i<NORDER;i+)yi=ysi+1.0/6*(2*(rk0i+rk2i)+4*rk1i+di*dt);voidmain(void)charc;dofloatdt,tmax,

32、co,tnext,tol;intj;printf("Thisprogramwillfindtherootof:n");printf("y0'(t)=-p0*y0*y1n");printf("&&y1'(t)=p0*y0*y1-p1*y1n");printf("&&y2'(t)=p1*y1n");for(j=0;j<NORDER;j+)printf("Now,pleaseinputthefirsty%d:",j);scanf("

33、;%f",&yj);for(j=0;j<NPARAM;j+)printf("Pleaseinputp%d:",j);scanf("%f",&pj);printf("Pleaseinputthetimeofsimulation:");scanf("%f",&tmax);printf("Pleaseinputthetimeofonestep:");scanf("%f",&dt);printf("Now,thisisther

34、esult:n");co=5*dt;tol=0.0001*co;tnext=0;output();for(t=0;t<=tmax+tol;)intk;if(fabs(tnext-t)<tol)tnext+=co;printf("%10.4f",t);for(k=0;k<NORDER;k+)printf("%10d",(int)(yk+0.5);putchar('n');run(dt);printf("nRunthisprogramagainY:");c=getche();putchar(&#

35、39;n');while(c='r'|c='Y'|c='y');程序中所用的變量和數(shù)組說明如下：Y(20)狀態(tài)變量數(shù)組G(20)狀態(tài)變量的一階導(dǎo)數(shù)P(1)=P1=0.1FALSEP(20)存數(shù)系統(tǒng)參數(shù)，本例中僅一個(gè)參數(shù)，TMAX模擬時(shí)間DT積分步長，即前面的hNP打印點(diǎn)數(shù)NORDER系統(tǒng)階數(shù)NPARAM系統(tǒng)參數(shù)個(gè)數(shù)OUTPUT打印輸出控制變量，.TRUE.為打印為不打印INIT打印表頭控制變量A(8)8個(gè)參數(shù)變量CO打印時(shí)間間隔NOUT已打印點(diǎn)數(shù)計(jì)數(shù)器TNEXT打印點(diǎn)處的時(shí)間值SS剩下的打印點(diǎn)數(shù)NLIST輸出結(jié)果對以上程序編譯，運(yùn)行，得到

36、如下結(jié)果TIMEY(1)Y(2)0.00000.00000.00000.50000.12040.46761.00000.44500.80081.50000.88630.92652.00001,33320.82472.50001.67880.53109.50001.66660.833310.00001.66660.83332.6 變步長的龍格一庫塔法以上提到，步長h的選擇十分重要，h太大不能達(dá)到預(yù)期的精度要求，h太小則增加了模擬時(shí)間，且有可能影響計(jì)算精度，要克服這一問題，就要采用變步長方法，使計(jì)算量盡可能的小，且精度也合乎要求。步長的自動控制是根據(jù)每一步的誤差的大小來實(shí)現(xiàn)的。為了得到每一步的局部

37、誤差的估計(jì)值，可以采用兩種不同階次的遞推公式（一般采用低一階的RK公式，同時(shí)計(jì)算yn+1并取Ki相同，使中間結(jié)果1957 年 Merson （默(2.15)兩者之差），要使計(jì)算量最少，可以選擇RK系數(shù)，要求兩個(gè)公式中的對兩種RK公式都適用，則這兩個(gè)公式的計(jì)算結(jié)果之差就作為誤差估計(jì)。森）給出了一個(gè)四階變步長的方法，稱為龍格一庫塔默森法。其四階公式為yn+1=yn+h(K1+4K4+K5)/6K1=f(tn,yn)K2=f(tn+h/3,yn+(h/3)*K1)K3=f(tn+h/3,yn+(h/6)(K1+K2)K4=f(tn+h/3,yn+(h/8)(K1+3K3)K5=f(tn+h,yn+(

38、h/2)(K1-3K2+4k4)計(jì)算估計(jì)誤差的三階公式如下n+=yn+?(3K1-9K3+12K4)6其絕對誤差為一hE=yn1-yn1(2K1-9K38K4-七)6這一算法是四階精度，三階誤差，故稱為RKM34法，目前使用較多。其穩(wěn)定域和一般RK4法相近，缺點(diǎn)是計(jì)算量稍大，每步計(jì)算5次f值，除用絕對誤差外，也可用相對誤差，最大相對誤差RE為RE = max(Eyn 1 - yn在每一步計(jì)算后，先計(jì)算誤差，根據(jù)誤差的大小來控制步長，控制策略如下：(1)當(dāng)誤差ERR大于預(yù)先設(shè)定的最大允許誤差Emax時(shí)，則縮小步長，一般將步長減半,并以新步長重新計(jì)算以后再比較。(2)當(dāng)誤差ERR小于預(yù)先定的最小允

39、許誤差Emm時(shí)，步長增加一倍，以新步長計(jì)算下(3)如步長小于某一下限hmin時(shí)不再減半，以免增加模擬時(shí)間，舍入誤差過大。(4)如步長大于某上限hmin(指打印間隔)，則不再加大步長。龍格一庫塔默森法雖然增加了計(jì)算量，但在微機(jī)日益普及的今天，這一方法是值得廣泛采用的。采用這種方法的程序RKM清單如下：#include<stdio.h>constintn=2;inlinefloatABS(floatinput)if(input>=0.0f)returninput;elsereturn-input;voidmain()voidrkm(float*t,floaty2,floateps

40、,floath,int*locp);intlocp=1,i;floatt=0.0f,eps=1.0E-8f,h=0.1f,y2;y0=1.0f;y1=1.0f;printf("%-1.8f,%-1.8f,%-1.8fn",t,y0,y1);for(i=1;i<=20;i+)rkm(&t,y,eps,h,&locp);printf("%-1.8f,%-1.8f,%-1.8fn",t,y0,y1);voidrkm(float*t,floaty2,floateps,floath,int*locp)voideql(floatt,floaty

41、2,floatf2);floatw5n,hc,er;intloc=1,lk=1,i;hc=h/(float)(*locp);label1:eql(*t,y,&w20);for(i=0;i<=n-1;i+)w0i=w2i*hc/3.0f+yi;eql(hc/3.0f+*t,&w00,&w3。);for(i=0;i<=n-1;i+)w0i=(w2i+w3i)*hc/6.0f+yi;eql(hc/3.0f+*t,&w00,&w3。);for(i=0;i<=n-1;i+)w0i=(w2i+3.0f*w3i)*hc*0.125f+yi;eql(h

42、c*0.5f+*t,&w0,&w40);for(i=0;i<=n-1;i+)w0i=w2i*hc*0.5f-w3i*hc*1.5f+w4i*hc*2.0f+yi;eql(hc+*t,&w00,&w30);for(i=0;i<=n-1;i+)w1i=(w2i+w3i+4.0f*w4i)*hc/6.0f+yi;lk=1;for(i=0;i<=n-1;i+)er=ABS(w0i-w1i)*0.2f;if(ABS(w1i)>1.0f)er=er/ABS(w1i);if(er>=eps)hc=hc*0.5f;*locp=(*locp)*2;l

43、oc=2*loc;gotolabel_1;if(er*64.0f>eps)lk=0;*t=(*t)+hc;for(i=0;i<=n-1;i+)yi=w1i;loc+;if(loc>=*locp)return;if(lk=0)|(loc%2=1)gotolabel_1;hc=2.0f*hc;10c=loc/2;*locp=*locp/2;gotolabel_1;voideql(floatt,floaty2,floatf2)f0=1.0f/(y1-t);f1=1.0f-1.0f/y0;使用RKM程序計(jì)算微分方程組dy_1dty2-t蛆二1dty1當(dāng)t=02時(shí)，y1It=1y21t

44、_0=1步長取0.1程序中的變量與數(shù)組為:N方程個(gè)數(shù)T模擬時(shí)間H積分區(qū)間Y(2)存放因變量初值與積分結(jié)果EPS允許的誤差LOCP表示在為h的積分區(qū)間等分?jǐn)?shù)W(N,5)工作單元，二維數(shù)組運(yùn)行以上程序得到如下，表2.2引出了精確解與R法與RKM法的比較表2.2t精確解RK(H=0.1)MS(H=0.1,EPS=10e-8)01.000000001.000000001.000000001.000000001.000000001.000000000.21.221402761.2214022041.221402761.018730751.018731221.018730750.41.491824701.

45、491822941.491824701.070320041.07320811.070320040.61.822118801.82115591.822118801.148811631.148812581.148811630.82.225540932.718273902.225540931,249328961.249329991.249328961.367879441.367880491.367879442.07.389056107.389013487.389056112.135335282.135336042.135335282.7 線性多步法上幾節(jié)我們討論了單步法，即在計(jì)算yn+1值時(shí)，只要知

46、道前一步的yn和fn的值就可以進(jìn)行計(jì)算，而每往后計(jì)算一步都可以使用前面已經(jīng)算出的結(jié)果。多步法就不同了，在計(jì)算yn+1時(shí)，不僅要用到前一步y(tǒng)n和fn的值，還要使用yn+1,yn+2及fn-1,fn-2的值，一般而言，充分利用前面多步信息來計(jì)算yn+1,可期望加快模擬速度，也會使精度得到提高，這顯然比單步法要好一些，這就是構(gòu)成多步法計(jì)算的基本思想。線性多步法中以亞當(dāng)姆斯(Adams)法使用得最為普遍，下面介紹這種方法。歐拉法在計(jì)算函數(shù)y時(shí)，在tn-tn-1中積分，用矩形公式求得。tn1y(tn1)=yntf(t,y)dttnuyn.f(tn,yn)(tn4tn)=ynf(tn,Yn*其過程和示意圖

47、在前面已經(jīng)給出。歐拉法計(jì)算誤差比較大，為了提高計(jì)算精度，將矩形改為梯形，即在求積分時(shí)使用下式tn1htf(t,y)dt：-f(tn1,yn1)f(tn,yn)tn2=h(fn1fn)(2.16)2則整個(gè)數(shù)值積分公式變成h(2-17)y(tn1)yn1=yn-f(tn1,yn1)f(tn,yn)1在(2-17)式中，為求yn+1值時(shí)，要用到y(tǒng)n+1本身，而這是一個(gè)未知數(shù)，故此式稱為二階隱式亞當(dāng)姆斯公式，也稱為梯形公式。(1)h(0).yn=yn2f(tn1，yn1)f(tn，yn)Adams法的幾何意義可用下圖表示用梯形公式進(jìn)行數(shù)值積分圖2-4f,時(shí)由于要用到 yn*本身，那就要用迭代法求解，即

48、先猜出一個(gè)初值(0)y n+ ，y：) = yn。£ (tn 1,yn0)1) f (tn,yn)然后用下式求出yn1)1再用y，1求出yn21。)h (八yn1 7n , 2 f(tnd，yn1) , f(tn，yn4)(2.18)如此重復(fù)下去，直到y(tǒng)nn；與yn、1)的差很小為止，此時(shí)可以認(rèn)為yn? = yn+，其截?cái)嗾`差為h2還有二階亞當(dāng)姆斯顯式公式，其可表示如下h _y(tn 1) ： yn 1 -3f(tn,yn) - f(tn4,yn)(2-19)yn1 =yn hB 4f (tn1,yn1B0 f (tn,yn)B-f(tnA1,yn*1)(2 - 20)表2-1到出了

49、各階亞當(dāng)姆斯法計(jì)算公式中的系數(shù)值，這是一個(gè)多步公式，為了計(jì)算yn+i,不僅要知道yn,還要知道yn-1,yn-2,，因此這些值都要保存，階數(shù)越搞，要用到的值越多，存儲量也就大了。多步法的缺點(diǎn)是不能自啟動，即開始時(shí)要用單步法，然后才用多步法。名稱B1B。B1B2B3一階顯式01000二階顯式03/2-1/200三階顯式023/12-16/125/120四階顯式055/24-59/2437/24-9/24:一階隱式10000二階隱式1/21/2000三階隱式5/128/12-1/1200四階隱式19/2419/24-5/241/240表2.1亞當(dāng)姆斯系數(shù)表在隱式多步法中，一般先用顯式多步法計(jì)算處值

50、，然后用隱式多步法作一次校正計(jì)算，這樣使計(jì)算過程得到簡化。例如，用一階顯式計(jì)算初值，即yn01=ynhf(tn,yn)(2-21)yn1=ynhf(tn,yn)f(tnh),y:0:(2-22)2以上計(jì)算公式與二階龍格一庫塔法公式相同。這種在多步法公式中先用顯式公式估計(jì)一個(gè)值的方法叫預(yù)估一校正法，它的優(yōu)點(diǎn)是在精度相同(四階)，步長h相同時(shí)，預(yù)估一校正法，每步計(jì)算兩次右函數(shù)，而龍格一庫塔法要計(jì)算四次右函數(shù)，因此這種方法應(yīng)用較為普遍。四階亞當(dāng)姆斯預(yù)估一校正法公式為：預(yù)估yF1=ynT-(55fn-59fnd37fn_2-9fnJ3)24校正YC+=yn+?(9爐蟲+19fn5fn+fnj)24亞當(dāng)

51、姆斯的顯示公式與隱式公式的特點(diǎn)如下：(1)相同階數(shù)的隱式公式的系數(shù)值比顯式公式小(一階例外)，因而隱式公式比顯式公式精確得多.(2)隱式公式的穩(wěn)定性亦比顯式的好。下表給出了Adams法的穩(wěn)定區(qū)。對于同樣階次，隱式的穩(wěn)定域比顯式的要大.隨著階數(shù)的增加，Adams法的穩(wěn)定域逐步減小.階數(shù)1234顯式-216/11-3/10-8-8-6-3(3)從計(jì)算上看，顯式比隱示計(jì)算量小。隱示公式需要計(jì)算fn+1,通常需要用Adams顯式公式對它提供一個(gè)首次近似值yn+1。這種將顯式公式和隱示公式聯(lián)合起來使用以改進(jìn)穩(wěn)定性和精度的方法就是上面介紹的預(yù)估一校正法。亞當(dāng)姆斯法每步只需要計(jì)算一次f函數(shù)值，因而計(jì)算量較小

52、，這是線性多步法的共同優(yōu)點(diǎn)。但它需要知道多個(gè)出發(fā)值才能計(jì)算，而微分方程只能提供一個(gè)初值，顯然多步法不象單步法那樣可以自啟動，它必須用其他方法先獲得所求時(shí)刻以前多步的解。為了保證出發(fā)值的計(jì)算精度，一般采用比較高階的單步法，以較小步長計(jì)算。此外，多步法改變步長不能象單步法那樣方便。下面介紹一個(gè)使用亞當(dāng)姆斯法進(jìn)行模擬的程序，此程序可用于單輸入單輸出和多輸入多輸出線性系統(tǒng)的模擬，在模擬時(shí)面向狀態(tài)方程。線性定常系統(tǒng)的狀態(tài)空間表達(dá)式包括兩個(gè)矩陣方程，X(t) = AX(t)Bu(T)(Y(t)=CX(t)+Du(t),X(t0)=X(0)第一個(gè)方程由n個(gè)一階微分方程組成，稱為狀態(tài)方程，第二個(gè)方程是一個(gè)線性代數(shù)方程，稱為輸出方程，式中X為n*l的狀態(tài)向量，u為m*1的控制向量。丫為l*l的輸出向量，A為n*n階狀態(tài)矩陣，由控制對象的參數(shù)決定，B為n*m的控制矩陣，C為l*n的輸出矩陣，D為l*m和直接傳輸矩陣。該程序用的計(jì)算方法為預(yù)估一校正線性多步法，首先用亞當(dāng)姆斯顯示公式計(jì)算預(yù)估值，再用亞當(dāng)姆斯隱式公式計(jì)算校正值，其計(jì)算公式為預(yù)估：Xn1=Xn'(55Xn-59*37Xn_2-9Xnj3)24校正：XC書=Xn+；h(9XnF+19Xn-5XnJXn/)24輸出：Yni=CX；.1DVAdambm函數(shù)流程圖如下:計(jì)算預(yù)報(bào)值