求解常微分方程組初值問題的龍格庫塔法分析及其C代碼

1、求解常微分方程組初值問題的龍格庫塔法分析及其C代碼1、概 述由高等數(shù)學(xué)的知識可知，一些特殊類型的常微分方程（組）能夠求出給定初始值的解析解，而在科學(xué)與工程問題中遇到的常微分方程（組）往往是極其復(fù)雜的，要想求得其給定初始值的解析解就變得極其困難，甚至是得不到解析解。盡管如此，在研究實際問題時，往往只需要獲得若干點上的近似值就行了，這就說明研究常微分方程（組）的數(shù)值解法是很有必要的。求解常微分方程（組）的數(shù)值解法有多種，比如歐拉法、龍格庫塔法、線性多步法等等，其中龍格庫塔法是這幾種方法中比較常用的，因為其計算精度較高且具有自啟動的特點。對于用龍格庫塔法求解單個常微分方程和求解常微分方程組的思路基本

2、相似（注意一點一個微分方程組是常微分方程組即表明微分方程中的各階導(dǎo)數(shù)都是對同一個變量求導(dǎo)，例如可以把各個量對時間求導(dǎo)得到一個常微分方程組，如果一個微分方程組中的有對不同變量的導(dǎo)數(shù)那么這個方程組就成了偏微分方程組），都是根據(jù)泰勒展開得到其迭代計算形式，其基本思想都是按照一定原則求取當(dāng)前點附近一些點的斜率，通過這些斜率的線性組合作為當(dāng)前點處的斜率，進(jìn)行遞推求解。2、數(shù)學(xué)模型2.1 常微分方程初值問題的數(shù)學(xué)模型 （1）式中為的已知函數(shù)，為給定的初始值。常微分方程的數(shù)值解法的任務(wù)就是要求出函數(shù)值在微分自變量取如下序列：時的值的近似值，一般情況下都采取等距結(jié)點的方式，即：，其中為相鄰兩結(jié)點的距離，稱為步

3、長。2.2 常微分方程組初值問題的數(shù)學(xué)模型 （2）式中為的已知函數(shù)，為初始值，利用數(shù)值解法求解常微分方程組的任務(wù)就是求出函數(shù)值在微分自變量取如下序列：時其值的近似值3 龍格庫塔法的迭代形式龍格庫塔法的迭代形式推導(dǎo)在相關(guān)的數(shù)值分析書籍中均有詳細(xì)的講解，這里就不進(jìn)行推導(dǎo)直接給出常用的四階龍格庫塔法的形式作為例子。3.1 單個常微分方程的迭代形式對應(yīng)于2.1節(jié)的數(shù)學(xué)模型有： （3）通過上面的迭代形式可知實際上就是如下四個點 （點1） （點2） （點3） （點4）處函數(shù)關(guān)于的斜率。這四個點的確定又是一環(huán)扣一環(huán)的，比如要求得點2的位置，就必須先求得在點1處的斜率3.2 常微分方程組的迭代形式對應(yīng)于2.2

4、節(jié)的數(shù)學(xué)模型有： （4）對比于單個微分方程的龍格庫塔迭代形式，可以看出整個求解的過程完全一樣，只不過在每步算斜率的時候需要一次性的解算完（方程的個數(shù)）個，即是要把的斜率一次性算完，才能繼續(xù)往下遞推運算。在求解常微分方程組初值問題的情況下，可將常微分方程組寫成向量的形式，這樣更加便于理解，下面用圖形來展示整個求解過程：圖1 龍格庫塔法的求解過程沿著箭頭的方向依次解出各個向量，即可依次遞推下去 ，下面的表示為點處的斜率，處情形類似。4 算 例4.1 問題描述利用四階龍格庫塔法求解初值問題取步長為。4.2 問題分析 為了和第3節(jié)的數(shù)學(xué)模型對應(yīng)，我們令，則上面的問題就變?yōu)楹颓懊娴哪Ｐ鸵粯拥男问?4.3

5、 C程序代碼#include #include #define M 2/ 維數(shù),方程組中方程的個數(shù)/ 算斜率的右端函數(shù)void RightSlopeFun( double *pRightK,/ 算出來的右端斜率項out double iX, / 計算點的微分自變量值 double *iY / 計算點的函數(shù)初始值 ) pRightK0 =iX*iY0-iY1;pRightK1 =(iX+iY0)/iY1;/ 龍格庫塔法解算程序double LgktSoulution( double *pOY, / 返回函數(shù)值out double iX, / 積分自變量輸入值 double *iY, / 初始值

6、double h, / 單步步長 int Count=1 / 解算步數(shù)，默認(rèn)形參為1，調(diào)用時不給此/ 參數(shù)賦值，則只解算一步 ) double K1M; / 第一個斜率double K2M; / 第二個斜率double K3M; / 第三個斜率double K4M; / 第四個斜率double tempYM; / 保存Y的中間值的臨時空間int i; / 計數(shù)器 double tempX; / 臨時保存x的空間tempX=iX; / 載入初始的x值 for (i=0;i0)RightSlopeFun(K1,tempX,tempY); / 解算斜率K1tempX += h/2;for (i=0;

7、iM;i+)pOYi=tempYi+K1i*h/2; RightSlopeFun(K2,tempX,pOY); / 解算斜率K2for (i=0;iM;i+)pOYi=tempYi+K2i*h/2;RightSlopeFun(K3,tempX,pOY); / 解算斜率K3tempX += h/2;for (i=0;iM;i+)pOYi=tempYi+K3i*h;RightSlopeFun(K4,tempX,pOY); / 解算斜率K4for (i=0;iM;i+)/ 得到函數(shù)值pOYi=tempYi+h*(K1i+2*K2i+2*K3i+K4i)/6;/ 推算了一步的函數(shù)值，此處可以對該數(shù)據(jù)進(jìn)

8、行相應(yīng)的操作/ 遞推數(shù)據(jù)更新Count-;if (Count=0)break;for (i=0;iM;i+) tempYi=pOYi; return tempX;void main() FILE *fp=fopen(d:lgktData.txt,a); double iyM; double ix; int i; iy0=1; iy1=2; ix=0; for(i=0;i500;i+) ix=LgktSoulution(iy,ix,iy,0.001); fprintf(fp,%-8.3f%-8.3f%-8.3fn,ix,iy0,iy1); fclose(fp);4.4 結(jié)果分析在運行完畢以上的C

9、程序后，在D盤的根目錄下就會生成名為lgktData.txt的文件，這個文件中就保存了龍格庫塔數(shù)值解法的數(shù)據(jù)，我們利用這個數(shù)據(jù)在MATLAB中繪圖，在MATLAB的命令窗口或者M(jìn)文件中輸入如下代碼clear clcx,y1,y2=textread(d:lgktData.txt,%.3f %.3f %.3f);plot(x,y1,r,x,y2,g);便可得到如下的圖形：從以上運用的四階龍格庫塔法求解常微分方程組的程序中，我們可以看出，只需要根據(jù)實際需要求解的模型，設(shè)定好方程組的個數(shù)M，將計算斜率的函數(shù)RightSlopeFun實現(xiàn)即可，求解函數(shù)LgktSoulution的內(nèi)容完全是固定不變的，除非你要使用其他形式龍格庫塔迭代形式；用龍格庫塔法求解單個微分方程和求解微分方程組的思路和計算流程完全一樣，只是在問題規(guī)模上有差異。5 小 結(jié)本文將單個常微分方程初值問題的龍格庫塔求解方法和常微分方