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

第2章 連續(xù)系統(tǒng)的計(jì)算機(jī)模擬 本章討論連續(xù)系統(tǒng)的模擬技術(shù),由于這類系統(tǒng)中狀態(tài)隨時(shí)間連續(xù)動(dòng)態(tài)地變化，常常具有一定的規(guī)律，故可用一些數(shù)學(xué)方程來描述，這些方程就是系統(tǒng)的數(shù)學(xué)模型，通常以微分方程、代數(shù)方程為多見。下面將介紹利用數(shù)值積分法對(duì)連續(xù)系統(tǒng)進(jìn)行數(shù)字模擬的基本原理和具體方法，并給出數(shù)值積分法中幾個(gè)常用的算法以及實(shí)現(xiàn)這些算法的計(jì)算程序，最后介紹兩個(gè)建模實(shí)例。 數(shù)值積分法不僅方法種類多，而且有較強(qiáng)的理論性，本章由淺入深地介紹幾種常見的數(shù)值解法。主要為單步法中的四階龍格-庫塔法與默森法和多步法中的亞當(dāng)斯法。 使用數(shù)字計(jì)算機(jī)對(duì)連續(xù)系統(tǒng)進(jìn)行模擬，首先必須將連續(xù)系統(tǒng)離散化，并將它轉(zhuǎn)化為差分方程，以建立所謂的模擬系統(tǒng)的數(shù)學(xué)模型。描述連續(xù)系統(tǒng)動(dò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)算精度高，對(duì)于不同的計(jì)算機(jī)，字長(zhǎng)一般在16位-72位之間，也可采用浮點(diǎn)運(yùn)算和雙精度運(yùn)算，其精度一般可達(dá)千萬分之一到百萬分之一。當(dāng)然結(jié)果的精度與所選的算法有關(guān)。這可以根據(jù)實(shí)際需要選擇機(jī)器、算法和模擬的步長(zhǎng)。數(shù)字計(jì)算機(jī)儲(chǔ)存容量大，可進(jìn)行各種運(yùn)算，在以前認(rèn)為是不可能解決的問題，利用數(shù)字計(jì)算機(jī)都可容易地或有可能得到解決。本章介紹的方法適應(yīng)性較強(qiáng)，應(yīng)用也十分廣泛。數(shù)字計(jì)算機(jī)上還有各種功能的軟件包（即一些子程序），用戶可以稍加修改或不經(jīng)修改就可以用于自己的模擬程序中，解決自己的實(shí)際問題，使用非常方便。 21 歐拉（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的分解速率正比與C的數(shù)量，即在任何時(shí)刻，如果a,b,和c是化學(xué)物質(zhì)A,B,和C的數(shù)量，即在任何時(shí)刻，如果a,b,和c是化學(xué)物質(zhì)A,B,C的數(shù)量，它們的增加和減少的速度服從下列微分方程。 （2.1） 其中K1和K2 是比例常數(shù)（一般而言這些比例常數(shù)會(huì)隨溫度和壓力發(fā)生變化，但在模擬過程中，為了簡(jiǎn)化模型，一般不允許其變化，故一律視為常數(shù)）。在給出常數(shù)K1 和K2 值以及A和B的數(shù)量（C=0）后，我們希望能確定有多少C物質(zhì)產(chǎn)在出來，這種化學(xué)反應(yīng)速率的決定在化學(xué)工業(yè)上是有意義的。 模擬該系統(tǒng)的一個(gè)直接的方法是在t0時(shí)開始，使t以t間隔增加。假定化學(xué)量在t時(shí)間步長(zhǎng)內(nèi)不變，而只能在t結(jié)束的瞬間發(fā)生變化，這樣在每個(gè)t結(jié)束時(shí)的A（或B或C）的數(shù)量就可以從t開始時(shí)的數(shù)值由下式求出 (22)同樣的方程b(t+t)和C(t+t),也可寫出。假定模擬周期為T,可將T分成N個(gè)小的時(shí)間步長(zhǎng)t,及 TNt 在時(shí)間為零時(shí)，我們知道a(0)、b(0)和c(0)的初始值,從這些初始值及常數(shù)K1和K2值出發(fā)，就可以計(jì)算出t時(shí)間內(nèi)的化學(xué)量的變化值。 a(t)=a(0)+K2C(0)-K1a(0)b(0)t b(t)=b(0)+K2C(0)-K1a(0)b(0)t c(t)=c(0)+2K1a(0)b(0)-2K2 c(0)t 使用這些值，又可計(jì)算系統(tǒng)的下一個(gè)狀態(tài)，即2t時(shí)的狀態(tài)。 a(2t)=a(t)+K2c (t)-K1a(t)b(t)t b(2t)=b(t)+K2c(t)-K1a(t )b(t )t c(2t)=c(t )+2K1a(t )b(t )-2K2 c(t)t 使用2t時(shí)的系統(tǒng)狀態(tài)，又可寫出3t時(shí)的狀態(tài)，依此類推，以t為間隔，進(jìn)行N步計(jì)算，就可寫出周期T的系統(tǒng)狀態(tài)得到所期望的結(jié)果，這個(gè)過程可用圖2.1表示。 輸入K1，K2，a(0),b(0),c(0) T, t和N I0 II1 保存a(I, t),b(I, t) c(I, t) Y N I=N? 打印a,b,c 模擬完畢 圖2.1 化學(xué)反應(yīng)例子模擬程序框圖 模擬程序使用C語言編寫，程序中的初值如下：K1=0.008/g.min K20.002/min, a(0)=100g, b(0)=50g C(0)=0, T=5mins, t=0.1min, N=50其程序清單如下：#include #include float k1,k2;static float A53,B53,C53,delta,t;void strut(int);main() int i; A1=100.0; B1=50.0; C1=0.0; t=0; delta=0.1; k1=0.008; k2=0.002; for(i=1;i53;i+) printf(%2d,i); printf(%10.2f,%10.2f,%10.2f,%10.2fn,t,Ai,Bi,Ci); strut(i); return; void strut(int i) 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é)果如下： 1 0.00, 100.00, 50.00, 0.00 2 0.10, 96.00, 46.00, 8.00 3 0.20, 92.47, 42.47, 15.06 4 0.30, 89.33, 39.33, 21.34 5 0.40, 86.52, 36.52, 26.95 6 0.50, 84.00, 34.00, 32.00 7 0.60, 81.72, 31.72, 36.55 8 0.70, 79.66, 29.66, 40.69 9 0.80, 77.77, 27.77, 44.4510 0.90, 76.05, 26.05, 47.8911 1.00, 74.48, 24.48, 51.0412 1.10, 73.03, 23.03, 53.9413 1.20, 71.70, 21.70, 56.6114 1.30, 70.46, 20.46, 59.0715 1.40, 69.32, 19.32, 61.3616 1.50, 68.26, 18.26, 63.4817 1.60, 67.28, 17.28, 65.4418 1.70, 66.36, 16.36, 67.2819 1.80, 65.51, 15.51, 68.9920 1.90, 64.71, 14.71, 70.5921 2.00, 63.96, 13.96, 72.0822 2.10, 63.26, 13.26, 73.4823 2.20, 62.60, 12.60, 74.7924 2.30, 61.99, 11.99, 76.0325 2.40, 61.41, 11.41, 77.1826 2.50, 60.86, 10.86, 78.2727 2.60, 60.35, 10.35, 79.3028 2.70, 59.87, 9.87, 80.2729 2.80, 59.41, 9.41, 81.1830 2.90, 58.98, 8.98, 82.0431 3.00, 58.57, 8.57, 82.8632 3.10, 58.19, 8.19, 83.6333 3.20, 57.82, 7.82, 84.3634 3.30, 57.48, 7.48, 85.0535 3.40, 57.15, 7.15, 85.7036 3.50, 56.84, 6.84, 86.3237 3.60, 56.55, 6.55, 86.9138 3.70, 56.27, 6.27, 87.4639 3.80, 56.00, 6.00, 87.9940 3.90, 55.75, 5.75, 88.5041 4.00, 55.51, 5.51, 88.9742 4.10, 55.29, 5.29, 89.4343 4.20, 55.07, 5.07, 89.8644 4.30, 54.86, 4.86, 90.2745 4.40, 54.67, 4.67, 90.6646 4.50, 54.48, 4.48, 91.0347 4.60, 54.31, 4.31, 91.3948 4.70, 54.14, 4.14, 91.7349 4.80, 53.98, 3.98, 92.0550 4.90, 53.82, 3.82, 92.3551 5.00, 53.68, 3.68, 92.6552 5.10, 53.54, 3.54, 92.93 以上例子的算法就是數(shù)值積分法中最簡(jiǎn)單的一種方法叫作歐拉法, 從這個(gè)例子中我們可以看出使用計(jì)算機(jī)解題的過程，由于歐拉法本身的特點(diǎn)，決定了其計(jì)算精度差, 現(xiàn)在幾乎無人在實(shí)際工作中使用，但它導(dǎo)出簡(jiǎn)單， 能說明構(gòu)造數(shù)值解法一般計(jì)算公式的基本思想, 模擬程序也清晰易懂，故人們樂意用它做為構(gòu)造數(shù)值解法的入門例子。其一般解法如下：設(shè)給定微分方程(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í)的積分值為 (2.5) 將上式寫成以下差分方程形式: （2.6）這就是歐拉公式。它是一個(gè)遞推的差分方程，任何一個(gè)新的數(shù)值解yn+1均是基于前一個(gè)數(shù)值解以及它的導(dǎo)數(shù)f(tn,yn)值求得的。只要給定初始條件y0及步長(zhǎng)h就可根據(jù)f(t0,y0)算出y1的值，再以y1算出y2，如此逐步算出y3, y4,，直到滿足所需計(jì)算的范圍才停止計(jì)算。歐拉法的基本思路是把連續(xù)的時(shí)間t分割成等間隔的t, 在這些離散的時(shí)刻算得函數(shù)值，根據(jù)這些值在函數(shù)圖上可得到一條折線，所以歐拉法又叫折線法，其特點(diǎn)是分析方法簡(jiǎn)單，計(jì)算量小，但計(jì)算精度低（后面將討論歐拉法與其它方法的比較）。下圖為歐拉折線法的幾何意義。 Y f(tn,yn) f(tn+1,yn+1) 歐拉折線 h y=f(t,y(t)dt+c 0 tn-1 tn tn+1 tn+2 圖2.2 歐拉法的幾何意義如果用梯形面積來代替一個(gè)小區(qū)間的曲線積分，就可克服以小矩形計(jì)算的缺點(diǎn)，提高精度，梯形法計(jì)算公式為(2.7)上式為隱式公式，因?yàn)楣接叶撕衴n+1, 這是未知的待求量，故梯形法不能自行啟動(dòng)運(yùn)算，而要依賴于其它算法的幫助，比如說用歐拉公式法求出初值，算出的近似值，然后將其帶入微分方程，計(jì)算的近似值，最后利用梯形公式反復(fù)迭代。如迭代一次后就認(rèn)為求得了近似解，這就是改進(jìn)的歐拉法，其公式為 預(yù)估 （2-8） 校正 （2-9）上式中第一個(gè)為預(yù)估公式，第二個(gè)為校正公式。通常這種方法稱為預(yù)估矯正法。在校正公式中計(jì)算了兩點(diǎn)的斜率，再求其平均值，故計(jì)算量比歐拉法要大些。23 數(shù)值積分的幾個(gè)基本概念 1. 單步法與多步法 數(shù)值積分法都用遞推公式求解, 而遞推公式是步進(jìn)式的, 當(dāng)由前一時(shí)刻的數(shù)值yn就能計(jì)算出后一時(shí)刻的數(shù)值yn+1時(shí), 這種方法稱為單步法, 它是一種能自啟動(dòng)的算法, 歐拉法是單步法. 反之, 如果求yn+1時(shí)要用到 yn, yn-1, yn-2 ,等多個(gè)值, 則稱為多步法,由于多步法在一開始時(shí)要用其它方法計(jì)算該時(shí)刻前面的函數(shù)值, 所以它不能自啟動(dòng),以上所講的改進(jìn)的歐拉法就是一個(gè)多步法的例子。 2 顯式與隱式 在遞推公式中, 計(jì)算yn+1時(shí)所用到的數(shù)據(jù)均已算出的計(jì)算公式稱為顯示公式. 相反,在算式中隱含有末知量yn+1 的則稱為隱含公式. 這需用另一個(gè)顯示公式估計(jì)一個(gè)值, 再用隱式公式進(jìn)行迭代運(yùn)算, 這就是預(yù)估-校正法, 這種方法不能自啟動(dòng). 用單步法與顯示公式在計(jì)算上比用多步法和隱式公式方便. 但有時(shí)要考慮精度與穩(wěn)定性時(shí), 則必須采用隱式公式運(yùn)算. 3 截?cái)嗾`差 一般使用臺(tái)勞級(jí)數(shù)來分析數(shù)值積分公式的精度. 為簡(jiǎn)化分析, 假定前一步得的結(jié)果yn是準(zhǔn)確的, 則用臺(tái)勞級(jí)數(shù)求得tn+1處的精確解為 (2.10)與前面的遞推公式比較, 歐拉法是從以上精確解中只取前兩 項(xiàng)之和來近似計(jì)算yn+1的 , 由這種方法單獨(dú)一步引進(jìn)的附加誤差通常稱著局部截?cái)嗾`差，它是該方法給出的值與微分方程解之間的差值，故又稱為局部離散誤差。歐拉法的局部截?cái)嗾`差為 不同的數(shù)值解法，其局部截?cái)嗾`差也不同。一般若截?cái)嗾`差為，則方法稱為r階的，所以方法的階數(shù)可以作為衡量方法精確度的一個(gè)重要標(biāo)志。歐拉法是一階精度的數(shù)值解法，而改進(jìn)的歐拉法(梯形法)相當(dāng)于取臺(tái)勞級(jí)數(shù)的前三項(xiàng)，故其解為二階精度。 4 舍入誤差 由于計(jì)算機(jī)的字長(zhǎng)有限, 數(shù)字不能表示得完全精確, 在計(jì)算過程中不可避免地會(huì)產(chǎn)生舍入誤差. 舍入誤差與步長(zhǎng) h 成反比, 如計(jì)算步長(zhǎng)小, 計(jì)算次數(shù)多則舍入誤差大. 產(chǎn)生舍入誤差的因素較多, 除了與計(jì)算機(jī)字長(zhǎng)有關(guān)外, 還與計(jì)算機(jī)所使用的數(shù)字系統(tǒng), 數(shù)的運(yùn)算次序及計(jì)算f(t,y)所用的程序的精確度等因素有關(guān). 5 數(shù)值解的穩(wěn)定性 以上對(duì)連續(xù)系統(tǒng)的模擬用到了差分方程, 把初始值代入遞推公式進(jìn)行迭代運(yùn)算, 如果原系統(tǒng)是穩(wěn)定的, 由數(shù)值積分法求得的解也應(yīng)該是穩(wěn)定的. 但是, 在計(jì)算機(jī)逐次計(jì)算時(shí), 初始數(shù)據(jù)的誤差及計(jì)算過程中的舍入誤差對(duì)后面的計(jì)算結(jié)果會(huì)產(chǎn)生影響. 而且計(jì)算步長(zhǎng)選擇得不合理時(shí), 有可能使模擬結(jié)果不穩(wěn)定. 以下我們簡(jiǎn)要地討論這一問題.差分方程的解與微分方程的解類似, 均可分為特解與通解兩部分. 與穩(wěn)定性有關(guān)的為通解部分, 它取決于該差分方程的特征根是否滿足穩(wěn)定性要求。例如, 為了考慮歐拉法的穩(wěn)定性，可用檢驗(yàn)方程y=y。其中為方程的特征根。對(duì)此方程有 要使該微分方程穩(wěn)定，須使下式成立 |1+h|1有此可得 |h|2或h8 精度階數(shù)234456n-2可見精度的階數(shù)與計(jì)算函數(shù)值f 的次數(shù)之間的關(guān)系并非等量增加。 其精度最低。 如果將h取得很小，能否用歐拉公式得到很高的精度呢？從理論上講是可以的，但是實(shí)際上由于計(jì)算機(jī)字長(zhǎng)有限，在計(jì)算中有舍入誤差。步長(zhǎng)h越小，計(jì)算的步數(shù)越多，舍入誤差的積累會(huì)越來越大，故用歐拉法時(shí)不能用很小的h。 2.5 四階龍格庫塔法模擬程序及應(yīng)用目前已有多種四階龍格庫塔法模擬（仿真）軟件包推出，雖然子程序稍有不同，但總的結(jié)構(gòu)還是相同的。對(duì)于連續(xù)系統(tǒng)的龍格庫塔法的計(jì)算機(jī)程序，其大致結(jié)構(gòu)如下圖所示： 主 函 數(shù)輸 入 及 運(yùn) 行 輸 出 初 始 化 微分方程 數(shù)值積分 顯 示 打印 繪圖 右函數(shù) 圖2.3 龍格庫塔法程序簡(jiǎn)要框圖 圖2.3中每一個(gè)程序模塊的功能為：1. 主函數(shù)：主函數(shù)實(shí)現(xiàn)對(duì)整個(gè)模擬系統(tǒng)的控制, 這是通過調(diào)用各個(gè)子程序?qū)崿F(xiàn)的。 2. 輸入及初始化函數(shù)：主要對(duì)系統(tǒng)參數(shù)輸入初值，例如模擬時(shí)間，積分步長(zhǎng)，方程階數(shù)等。這可以從鍵盤輸入，也可從數(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)在再將它寫成可編程的向量形式式中，i =1,2,3,.N, N為方程階數(shù)。程序求得，這些分別為變量的導(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.5DT RK2(I)=D(I)*0.5DT RK3(I)=D(I)*DT代入上式： Y(I) =YS(I) +0.166667*DT*D (I)+2.0*D (I)+2.0*D (I)+D (I) 即 yn+1 =yn +h (K1+2K2+2K3+K4)設(shè)有一個(gè)微分方程： y(t) +P1 y(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 積分步長(zhǎng)為0.1 打印點(diǎn)數(shù)為101程序清單如下：主程序：/*p23 example*/#include stdio.h#include conio.h#include math.h#define NORDER 3#define NPARAM 2float yNORDER,dNORDER,pNPARAM,t;void output(void)int j; printf(%9s,time); for(j=0;jNORDER;j+) printf( y%d,j); putchar(n);void difq(void) d0=-p0*y0*y1; d1=p0*y0*y1-p1*y1; d2=p1*y1; void run(float dt) int i; float rk3NORDER,ysNORDER; difq(); for(i=0;iNORDER;i+) rk0i=di*0.5*dt; ysi=yi; yi=ysi+rk0i; t+=0.5*dt; difq(); for(i=0;iNORDER;i+) rk1i=di*0.5*dt; yi=ysi+rk1i; difq(); for(i=0;iNORDER;i+) rk2i=di*dt; yi=ysi+rk2i; t+=0.5*dt; difq(); for(i=0;iNORDER;i+) yi=ysi+1.0/6*(2*(rk0i+rk2i)+4*rk1i+di*dt); void main(void)char c; do float dt,tmax,co,tnext,tol; int j; printf(This program will find the root of:n); printf( y0(t)=-p0*y0*y1n); printf(&y1(t)=p0*y0*y1-p1*y1n); printf(&y2(t)=p1*y1n); for(j=0;jNORDER;j+) printf(Now,please input the first y%d:,j); scanf(%f,&yj); for(j=0;jNPARAM;j+) printf(Please input p%d:,j); scanf(%f,&pj); printf(Please input the time of simulation:); scanf(%f,&tmax); printf(Please input the time of one step:); scanf(%f,&dt); printf(Now,this is the result:n); co=5*dt; tol=0.0001*co; tnext=0; output(); for (t=0;t=tmax+tol;) int k; if(fabs(tnext-t)tol) tnext+=co; printf(%10.4f,t); for(k=0;kNORDER;k+)printf(%10d,(int)(yk+0.5); putchar(n); run(dt); printf(nRun this program againY:); c=getche(); putchar(n); while(c=r|c=Y|c=y);程序中所用的變量和數(shù)組說明如下：Y(20) 狀態(tài)變量數(shù)組G(20) 狀態(tài)變量的一階導(dǎo)數(shù)P(20) 存數(shù)系統(tǒng)參數(shù)，本例中僅一個(gè)參數(shù)，P(1)=P1=0.1TMAX 模擬時(shí)間DT 積分步長(zhǎng)，即前面的hNP 打印點(diǎn)數(shù)NORDER 系統(tǒng)階數(shù)NPARAM 系統(tǒng)參數(shù)個(gè)數(shù)OUTPUT 打印輸出控制變量，.TRUE.為打印 FALSE 為不打印INIT 打印表頭控制變量A(8) 8個(gè)參數(shù)變量CO 打印時(shí)間間隔NOUT 已打印點(diǎn)數(shù)計(jì)數(shù)器TNEXT 打印點(diǎn)處的時(shí)間值SS 剩下的打印點(diǎn)數(shù)NLIST 輸出結(jié)果對(duì)以上程序編譯，運(yùn)行，得到如下結(jié)果 TIME Y（1） Y（2） 0.0000 0.0000 0.0000 0.5000 0.1204 0.4676 1.0000 0.4450 0.8008 1.5000 0.8863 0.9265 2.0000 1,3332 0.8247 2.5000 1.6788 0.5310 9.5000 1.6666 0.8333 10.0000 1.6666 0.8333 2.6 變步長(zhǎng)的龍格庫塔法以上提到，步長(zhǎng)h的選擇十分重要，h太大不能達(dá)到預(yù)期的精度要求，h太小則增加了模擬時(shí)間，且有可能影響計(jì)算精度，要克服這一問題，就要采用變步長(zhǎng)方法，使計(jì)算量盡可能的小，且精度也合乎要求。步長(zhǎng)的自動(dòng)控制是根據(jù)每一步的誤差的大小來實(shí)現(xiàn)的。為了得到每一步的局部誤差的估計(jì)值，可以采用兩種不同階次的遞推公式（一般采用低一階的RK公式，同時(shí)計(jì)算yn+1并取兩者之差），要使計(jì)算量最少，可以選擇系數(shù)，要求兩個(gè)公式中的Ki相同，使中間結(jié)果對(duì)兩種RK公式都適用，則這兩個(gè)公式的計(jì)算結(jié)果之差就作為誤差估計(jì)。1957年Merson（默森）給出了一個(gè)四階變步長(zhǎng)的方法，稱為龍格庫塔默森法。 其四階公式為 yn+1 =yn +h (K1+4K4+K5)/6 (2.15) K1 =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+ (h/2)(K1-3K2+4k4)計(jì)算估計(jì)誤差的三階公式如下 其絕對(duì)誤差為這一算法是四階精度，三階誤差，故稱為RKM34法，目前使用較多。其穩(wěn)定域和一般RK4法相近，缺點(diǎn)是計(jì)算量稍大，每步計(jì)算5次f值，除用絕對(duì)誤差外，也可用相對(duì)誤差，最大相對(duì)誤差RE為 在每一步計(jì)算后，先計(jì)算誤差，根據(jù)誤差的大小來控制步長(zhǎng)，控制策略如下： （1）當(dāng)誤差ERR大于預(yù)先設(shè)定的最大允許誤差Emax時(shí)，則縮小步長(zhǎng)，一般將步長(zhǎng)減半，并以新步長(zhǎng)重新計(jì)算以后再比較。 （2）當(dāng)誤差ERR小于預(yù)先定的最小允許誤差Emin時(shí)，步長(zhǎng)增加一倍，以新步長(zhǎng)計(jì)算下一步并計(jì)算誤差。 （3）如步長(zhǎng)小于某一下限hmin時(shí)不再減半，以免增加模擬時(shí)間，舍入誤差過大。 （4）如步長(zhǎng)大于某上限hmin（指打印間隔），則不再加大步長(zhǎng)。龍格庫塔默森法雖然增加了計(jì)算量，但在微機(jī)日益普及的今天，這一方法是值得廣泛采用的。采用這種方法的程序RKM清單如下：#include const int n=2;inline float ABS(float input)if (input=0.0f) return input;else return -input;void main()void rkm(float *t,float y2,float eps,float h,int *locp); int locp=1,i; float t=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);void rkm(float *t,float y2,float eps,float h,int *locp)void eql(float t,float y2,float f2); float w5n,hc,er; int loc=1,lk=1,i; hc=h/(float)(*locp);label_1: eql(*t,y,&w20); for (i=0;i=n-1;i+) w0i=w2i*hc/3.0f+yi; eql(hc/3.0f+*t,&w00,&w30); for (i=0;i=n-1;i+) w0i=(w2i+w3i)*hc/6.0f+yi; eql(hc/3.0f+*t,&w00,&w30); for (i=0;i=n-1;i+) w0i=(w2i+3.0f*w3i)*hc*0.125f+yi; eql(hc*0.5f+*t,&w00,&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;i1.0f) er=er/ABS(w1i); if (er=eps) hc=hc*0.5f; *locp=(*locp)*2; loc=2*loc; goto label_1; if (er*64.0feps) lk=0; *t=(*t)+hc; for (i=0;i=*locp) return; if (lk=0)|(loc%2=1) goto label_1; hc=2.0f*hc; loc=loc/2; *locp=*locp/2; goto label_1;void eql(float t,float y2,float f2)f0=1.0f/(y1-t); f1=1.0f-1.0f/y0;使用RKM程序計(jì)算微分方程組當(dāng)時(shí)， 步長(zhǎng)取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.2 t 精確解 RK(H=0.1) MS(H=0.1,EPS=10e-8) 0 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000 0.2 1.22140276 1.221402204 1.22140276 1.01873075 1.01873122 1.01873075 0.4 1.49182470 1.49182294 1.49182470 1.07032004 1.0732081 1.07032004 0.6 1.82211880 1.8211559 1.82211880 1.14881163 1.14881258 1.14881163 0.8 2.22554093 2.71827390 2.22554093 1,24932896 1.24932999 1.24932896 1.36787944 1.36788049 1.36787944 2.0 7.38905610 7.38