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

第2章 連續(xù)系統(tǒng)的計算機模擬 本章討論連續(xù)系統(tǒng)的模擬技術,由于這類系統(tǒng)中狀態(tài)隨時間連續(xù)動態(tài)地變化，常常具有一定的規(guī)律，故可用一些數(shù)學方程來描述，這些方程就是系統(tǒng)的數(shù)學模型，通常以微分方程、代數(shù)方程為多見。下面將介紹利用數(shù)值積分法對連續(xù)系統(tǒng)進行數(shù)字模擬的基本原理和具體方法，并給出數(shù)值積分法中幾個常用的算法以及實現(xiàn)這些算法的計算程序，最后介紹兩個建模實例。 數(shù)值積分法不僅方法種類多，而且有較強的理論性，本章由淺入深地介紹幾種常見的數(shù)值解法。主要為單步法中的四階龍格-庫塔法與默森法和多步法中的亞當斯法。 使用數(shù)字計算機對連續(xù)系統(tǒng)進行模擬，首先必須將連續(xù)系統(tǒng)離散化，并將它轉化為差分方程，以建立所謂的模擬系統(tǒng)的數(shù)學模型。描述連續(xù)系統(tǒng)動態(tài)特征的數(shù)學模型是多種多樣的，除微分方程外，還有傳遞函數(shù)、結構圖及狀態(tài)方程等，由于篇幅所限，本書不討論后兩種方法。 建立系統(tǒng)的模擬模型之后，就要選擇計算機語言（也叫算法語言）編寫系統(tǒng)模擬程序，在計算機上運行，將結果保留在數(shù)據(jù)文件中以待傳輸和處理。由于模擬的目的不同，可以選用不同的模擬模型和算法，其特點是運算精度高，對于不同的計算機，字長一般在16位-72位之間，也可采用浮點運算和雙精度運算，其精度一般可達千萬分之一到百萬分之一。當然結果的精度與所選的算法有關。這可以根據(jù)實際需要選擇機器、算法和模擬的步長。數(shù)字計算機儲存容量大，可進行各種運算，在以前認為是不可能解決的問題，利用數(shù)字計算機都可容易地或有可能得到解決。本章介紹的方法適應性較強，應用也十分廣泛。數(shù)字計算機上還有各種功能的軟件包（即一些子程序），用戶可以稍加修改或不經(jīng)修改就可以用于自己的模擬程序中，解決自己的實際問題，使用非常方便。 21 歐拉（Euler）法 在討論連續(xù)系統(tǒng)的計算機模擬之前，讓我們先看一個化學反應的例子，通過這個例子我們可以看到怎樣使用數(shù)字計算機模擬一個實際問題，雖然介紹的是歐拉法，但是分析問題的思路同樣適用與其它數(shù)值積分法。 當兩種物質A和B放到一起產生化學反應時，產生第三種物質C，一般一克A與一克B結合產生2克的C物質，形成C 的速率與A 和B的數(shù)量乘積成正比，同樣C也可分解為A和B，C的分解速率正比與C的數(shù)量，即在任何時刻，如果a,b,和c是化學物質A,B,和C的數(shù)量，即在任何時刻，如果a,b,和c是化學物質A,B,C的數(shù)量，它們的增加和減少的速度服從下列微分方程。 （2.1） 其中K1和K2 是比例常數(shù)（一般而言這些比例常數(shù)會隨溫度和壓力發(fā)生變化，但在模擬過程中，為了簡化模型，一般不允許其變化，故一律視為常數(shù)）。在給出常數(shù)K1 和K2 值以及A和B的數(shù)量（C=0）后，我們希望能確定有多少C物質產在出來，這種化學反應速率的決定在化學工業(yè)上是有意義的。 模擬該系統(tǒng)的一個直接的方法是在t0時開始，使t以t間隔增加。假定化學量在t時間步長內不變，而只能在t結束的瞬間發(fā)生變化，這樣在每個t結束時的A（或B或C）的數(shù)量就可以從t開始時的數(shù)值由下式求出 (22)同樣的方程b(t+t)和C(t+t),也可寫出。假定模擬周期為T,可將T分成N個小的時間步長t,及 TNt 在時間為零時，我們知道a(0)、b(0)和c(0)的初始值,從這些初始值及常數(shù)K1和K2值出發(fā)，就可以計算出t時間內的化學量的變化值。 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 使用這些值，又可計算系統(tǒng)的下一個狀態(tài)，即2t時的狀態(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時的系統(tǒng)狀態(tài)，又可寫出3t時的狀態(tài)，依此類推，以t為間隔，進行N步計算，就可寫出周期T的系統(tǒng)狀態(tài)得到所期望的結果，這個過程可用圖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 化學反應例子模擬程序框圖 模擬程序使用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; 運行此程序，輸出模擬結果如下： 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ù)值積分法中最簡單的一種方法叫作歐拉法, 從這個例子中我們可以看出使用計算機解題的過程，由于歐拉法本身的特點，決定了其計算精度差, 現(xiàn)在幾乎無人在實際工作中使用，但它導出簡單， 能說明構造數(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時的積分值為 (2.5) 將上式寫成以下差分方程形式: （2.6）這就是歐拉公式。它是一個遞推的差分方程，任何一個新的數(shù)值解yn+1均是基于前一個數(shù)值解以及它的導數(shù)f(tn,yn)值求得的。只要給定初始條件y0及步長h就可根據(jù)f(t0,y0)算出y1的值，再以y1算出y2，如此逐步算出y3, y4,，直到滿足所需計算的范圍才停止計算。歐拉法的基本思路是把連續(xù)的時間t分割成等間隔的t, 在這些離散的時刻算得函數(shù)值，根據(jù)這些值在函數(shù)圖上可得到一條折線，所以歐拉法又叫折線法，其特點是分析方法簡單，計算量小，但計算精度低（后面將討論歐拉法與其它方法的比較）。下圖為歐拉折線法的幾何意義。 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 歐拉法的幾何意義如果用梯形面積來代替一個小區(qū)間的曲線積分，就可克服以小矩形計算的缺點，提高精度，梯形法計算公式為(2.7)上式為隱式公式，因為公式右端含有yn+1, 這是未知的待求量，故梯形法不能自行啟動運算，而要依賴于其它算法的幫助，比如說用歐拉公式法求出初值，算出的近似值，然后將其帶入微分方程，計算的近似值，最后利用梯形公式反復迭代。如迭代一次后就認為求得了近似解，這就是改進的歐拉法，其公式為 預估 （2-8） 校正 （2-9）上式中第一個為預估公式，第二個為校正公式。通常這種方法稱為預估矯正法。在校正公式中計算了兩點的斜率，再求其平均值，故計算量比歐拉法要大些。23 數(shù)值積分的幾個基本概念 1. 單步法與多步法 數(shù)值積分法都用遞推公式求解, 而遞推公式是步進式的, 當由前一時刻的數(shù)值yn就能計算出后一時刻的數(shù)值yn+1時, 這種方法稱為單步法, 它是一種能自啟動的算法, 歐拉法是單步法. 反之, 如果求yn+1時要用到 yn, yn-1, yn-2 ,等多個值, 則稱為多步法,由于多步法在一開始時要用其它方法計算該時刻前面的函數(shù)值, 所以它不能自啟動,以上所講的改進的歐拉法就是一個多步法的例子。 2 顯式與隱式 在遞推公式中, 計算yn+1時所用到的數(shù)據(jù)均已算出的計算公式稱為顯示公式. 相反,在算式中隱含有末知量yn+1 的則稱為隱含公式. 這需用另一個顯示公式估計一個值, 再用隱式公式進行迭代運算, 這就是預估-校正法, 這種方法不能自啟動. 用單步法與顯示公式在計算上比用多步法和隱式公式方便. 但有時要考慮精度與穩(wěn)定性時, 則必須采用隱式公式運算. 3 截斷誤差 一般使用臺勞級數(shù)來分析數(shù)值積分公式的精度. 為簡化分析, 假定前一步得的結果yn是準確的, 則用臺勞級數(shù)求得tn+1處的精確解為 (2.10)與前面的遞推公式比較, 歐拉法是從以上精確解中只取前兩 項之和來近似計算yn+1的 , 由這種方法單獨一步引進的附加誤差通常稱著局部截斷誤差，它是該方法給出的值與微分方程解之間的差值，故又稱為局部離散誤差。歐拉法的局部截斷誤差為 不同的數(shù)值解法，其局部截斷誤差也不同。一般若截斷誤差為，則方法稱為r階的，所以方法的階數(shù)可以作為衡量方法精確度的一個重要標志。歐拉法是一階精度的數(shù)值解法，而改進的歐拉法(梯形法)相當于取臺勞級數(shù)的前三項，故其解為二階精度。 4 舍入誤差 由于計算機的字長有限, 數(shù)字不能表示得完全精確, 在計算過程中不可避免地會產生舍入誤差. 舍入誤差與步長 h 成反比, 如計算步長小, 計算次數(shù)多則舍入誤差大. 產生舍入誤差的因素較多, 除了與計算機字長有關外, 還與計算機所使用的數(shù)字系統(tǒng), 數(shù)的運算次序及計算f(t,y)所用的程序的精確度等因素有關. 5 數(shù)值解的穩(wěn)定性 以上對連續(xù)系統(tǒng)的模擬用到了差分方程, 把初始值代入遞推公式進行迭代運算, 如果原系統(tǒng)是穩(wěn)定的, 由數(shù)值積分法求得的解也應該是穩(wěn)定的. 但是, 在計算機逐次計算時, 初始數(shù)據(jù)的誤差及計算過程中的舍入誤差對后面的計算結果會產生影響. 而且計算步長選擇得不合理時, 有可能使模擬結果不穩(wěn)定. 以下我們簡要地討論這一問題.差分方程的解與微分方程的解類似, 均可分為特解與通解兩部分. 與穩(wěn)定性有關的為通解部分, 它取決于該差分方程的特征根是否滿足穩(wěn)定性要求。例如, 為了考慮歐拉法的穩(wěn)定性，可用檢驗方程y=y。其中為方程的特征根。對此方程有 要使該微分方程穩(wěn)定，須使下式成立 |1+h|1有此可得 |h|2或h8 精度階數(shù)234456n-2可見精度的階數(shù)與計算函數(shù)值f 的次數(shù)之間的關系并非等量增加。 其精度最低。 如果將h取得很小，能否用歐拉公式得到很高的精度呢？從理論上講是可以的，但是實際上由于計算機字長有限，在計算中有舍入誤差。步長h越小，計算的步數(shù)越多，舍入誤差的積累會越來越大，故用歐拉法時不能用很小的h。 2.5 四階龍格庫塔法模擬程序及應用目前已有多種四階龍格庫塔法模擬（仿真）軟件包推出，雖然子程序稍有不同，但總的結構還是相同的。對于連續(xù)系統(tǒng)的龍格庫塔法的計算機程序，其大致結構如下圖所示： 主 函 數(shù)輸 入 及 運 行 輸 出 初 始 化 微分方程 數(shù)值積分 顯 示 打印 繪圖 右函數(shù) 圖2.3 龍格庫塔法程序簡要框圖 圖2.3中每一個程序模塊的功能為：1. 主函數(shù)：主函數(shù)實現(xiàn)對整個模擬系統(tǒng)的控制, 這是通過調用各個子程序實現(xiàn)的。 2. 輸入及初始化函數(shù)：主要對系統(tǒng)參數(shù)輸入初值，例如模擬時間，積分步長，方程階數(shù)等。這可以從鍵盤輸入，也可從數(shù)據(jù)文件輸入。3運行模塊: 在運行子程序中, 將反復調用數(shù)值積分和方程右函數(shù)模塊,進行迭代計算，可以每計算一步便顯示一次結果，也可以計算數(shù)步顯示一次結果，屏幕的顯示使用戶可以隨時監(jiān)視系統(tǒng)運行的進程，以便人工干預。4輸出子程序：按要求輸出打印結果，可以打印，也可以繪圖，視需要而定。四階龍格庫塔法公式前面已經(jīng)給出, 現(xiàn)在再將它寫成可編程的向量形式式中，i =1,2,3,.N, N為方程階數(shù)。程序求得，這些分別為變量的導數(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)設有一個微分方程： y(t) +P1 y(t) +y(t) =1.0令y1(t), y2(t)為兩各個狀態(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ù)初值： 模擬時間為50.0 積分步長為0.1 打印點數(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)變量的一階導數(shù)P(20) 存數(shù)系統(tǒng)參數(shù)，本例中僅一個參數(shù)，P(1)=P1=0.1TMAX 模擬時間DT 積分步長，即前面的hNP 打印點數(shù)NORDER 系統(tǒng)階數(shù)NPARAM 系統(tǒng)參數(shù)個數(shù)OUTPUT 打印輸出控制變量，.TRUE.為打印 FALSE 為不打印INIT 打印表頭控制變量A(8) 8個參數(shù)變量CO 打印時間間隔NOUT 已打印點數(shù)計數(shù)器TNEXT 打印點處的時間值SS 剩下的打印點數(shù)NLIST 輸出結果對以上程序編譯，運行，得到如下結果 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 變步長的龍格庫塔法以上提到，步長h的選擇十分重要，h太大不能達到預期的精度要求，h太小則增加了模擬時間，且有可能影響計算精度，要克服這一問題，就要采用變步長方法，使計算量盡可能的小，且精度也合乎要求。步長的自動控制是根據(jù)每一步的誤差的大小來實現(xiàn)的。為了得到每一步的局部誤差的估計值，可以采用兩種不同階次的遞推公式（一般采用低一階的RK公式，同時計算yn+1并取兩者之差），要使計算量最少，可以選擇系數(shù)，要求兩個公式中的Ki相同，使中間結果對兩種RK公式都適用，則這兩個公式的計算結果之差就作為誤差估計。1957年Merson（默森）給出了一個四階變步長的方法，稱為龍格庫塔默森法。 其四階公式為 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)計算估計誤差的三階公式如下 其絕對誤差為這一算法是四階精度，三階誤差，故稱為RKM34法，目前使用較多。其穩(wěn)定域和一般RK4法相近，缺點是計算量稍大，每步計算5次f值，除用絕對誤差外，也可用相對誤差，最大相對誤差RE為 在每一步計算后，先計算誤差，根據(jù)誤差的大小來控制步長，控制策略如下： （1）當誤差ERR大于預先設定的最大允許誤差Emax時，則縮小步長，一般將步長減半，并以新步長重新計算以后再比較。 （2）當誤差ERR小于預先定的最小允許誤差Emin時，步長增加一倍，以新步長計算下一步并計算誤差。 （3）如步長小于某一下限hmin時不再減半，以免增加模擬時間，舍入誤差過大。 （4）如步長大于某上限hmin（指打印間隔），則不再加大步長。龍格庫塔默森法雖然增加了計算量，但在微機日益普及的今天，這一方法是值得廣泛采用的。采用這種方法的程序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程序計算微分方程組當時， 步長取0.1程序中的變量與數(shù)組為：N 方程個數(shù)T 模擬時間H 積分區(qū)間 Y(2) 存放因變量初值與積分結果EPS 允許的誤差LOCP 表示在為h的積分區(qū)間等分數(shù)W(N,5) 工作單元，二維數(shù)組運行以上程序得到如下，表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