Visual C++常微分方程初值問題求解

　　摘要：本文講述了以計算機為輔助計算工具，分別使用歐拉算法、改進歐拉算法以及經典龍格-庫塔算法對常微分方程的初值問題進行數值求解的實現算法。

　　一、 引言

　　在工程計算中我們經常要去解一些常微分方程，雖然在高等數學和其他一些涉及微分方程的專業書籍中介紹了不少類型的常微分方程，及各自的解法。但工程技術人員在工程和科學研究中所關心的往往只是常微分方程的近似數值解，而非從事數學研究的技術人員所注重的"過程"。采用常規的人工推導、求解無疑是效率非常低下的，而且工程上的常微分方程往往結構非常復雜，要給出一般方程解的表達式也是非常困難的。實際上到目前為止，我們只能對有限的幾種特殊類型的方程求精確解，這遠不能滿足工程需要，對那些不能用初等函數來表達的方程就只能去求其近似的數值解，而且這樣還可以借助于運算速度快的計算機來進行輔助求解，大大提高求解的速度和精度，修改也比較靈活。

　　二、 使用歐拉算法及其改進算法進行一般求解

　　所謂的數值求解，就是求問題的解y(x)在一系列點上的值y(xi)的近似值yi。歐拉(Euler)算法是其中最基本、最簡單的算法，但其求解精度較低，一般不在工程中單獨進行計算。其實現的依據是用向前差商來近似代替導數。對于常微分方程：

　　dy/dx=f(x,y)，x∈[a,b]
　　y(a)=y0

　　可以將區間[a,b]分成n段，那么方程在第xI點有y'(xI)=f(xI,y(xI))，再用向前差商近似代替導數則為：(y(xI+1)-y(xI))/h= f(xI,y(xI))，因此可以根據xI點和yI點的數值計算出yI+1來：

　　yI+1= yI+h*f(xI ,yI)

　　下面就在Visual C++ 6.0編程環境下對一個簡單的常微分方程

　　y'=x-y+1,x∈[0,0.5]
　　y(0)=1

　　求近似數值解，由于該簡單方程可以用數學方法求得其精確描述式y(x)=x+e-x，所以可以據此檢驗近似數值解同真實解的誤差情況。對于其他一些結構復雜的常微分方程的數值解實現方法也是一樣的。下面就通過代碼來實現上述算法，并對計算結果作了比較：

float y[6]; file://用于存放計算出的常微分方程數值解 float r; file://同真實解的誤差情況 memset(y,0,sizeof(float)*6);//清零 y[0]=1; file://y(0)=1 …… for(float x=0;x<0.6;x+=0.1) file://區間分5段，步長為0.1 { r=x+expf(-x); file://真實解y(x)=x+e-x y[i+1]=y[i]+0.1*(x-y[i]+1); file://數值解（近似） r=fabs(r-y[i]); file://誤差 str.Format("y[%d]=%f r=%f\r\n",i,y[i],r); i++; msg+=str; } AfxMessageBox(msg); ……

　　經過程序計算，得出y(xi)在各點的近似數值解及各自同真實解的誤差，現列表如下，以茲對照：

xI（各分點）	yI （數值解）	y(xi) （真實值）	| y(xi)- yI | (誤差)
0.0	1.000000	1.000000	0.000000
0.1	1.000000	1.004837	0.004837
0.2	1.010000	1.018731	0.008731
0.3	1.029000	1.040818	0.011818
0.4	1.056100	1.070320	0.014220
0.5	1.090490	1.106531	0.016041

　　雖然從實驗結果看誤差不算太大，但這僅僅是一個用于實驗的非常簡單的常微分方程，對于實際工程中應用的結構復雜的方程其求解結果的誤差要遠比此大的多，由于還存在著局部截斷誤差和整體截斷誤差，有必要采取措施來抑制、減少誤差，盡量使結果精確。在構造歐拉公式時采取的一個重要步驟--用向前差商來代替導數，如將其改為向后差商也是行的通的。此時的歐拉公式就變成了：yI+1= yI+h*f(xI+1,yI+1)，由于該式是一個隱式公式，所以可用迭代法進行計算，直至獲取到滿足精度要求的yI+1。從數學上可以證明，該式的局部截斷誤差和前面的歐拉公式的截斷誤差在主部上之相差正負號，所以只要將顯示和隱式的兩個歐拉公式相加后再行求解會大大減少誤差?？梢越獾酶倪M后的歐拉公式的表達式為：

　　yI+1= yI+h*(f(xI, yI)+f(xI+1, yI+hf(xI,yI)))/2

　　對此式進行編程，就要比前面的代碼要麻煩些，需要分步對其進行計算，以達到最高的運算效率，減少運算量：

…… for(float x=0;x<0.6;x+=0.1) { r=x+expf(-x); T1=y[i]+0.1*(x-y[i]+1); file://分步進行計算 T2=y[i]+0.1*((x+0.1)-T1+1); y[i+1]=(T1+T2)/2; r=fabs(r-y[i]); str.Format("y[%d]=%f r=%f\r\n",i,y[i],r); i++; msg+=str; } AfxMessageBox(msg);

　　從下表得出的實驗數據可以看出，這種經過改進的歐拉算法所存在的誤差已大為減少，可以直接單獨應用于實際的工程計算。誤差的減少主要是由于先利用了歐拉公式對yI+1的值進行了預估，然后又利用梯形公式對預估值作了校正，從而在預估--校正的過程中減少了誤差。

xI（各分點）	yI （數值解）	y(xi) （真實值）	| y(xi)- yI | (誤差)
0.0	1.000000	1.000000	0.000000
0.1	1.005000	1.004837	0.000163
0.2	1.019025	1.018731	0.000294
0.3	1.041218	1.040818	0.000400
0.4	1.070802	1.070320	0.000482
0.5	1.107076	1.106531	0.000545

　　三、 使用經典龍格-庫塔算法進行高精度求解

　　龍格-庫塔(Runge-Kutta)方法是一種在工程上應用廣泛的高精度單步算法。由于此算法精度高，采取措施對誤差進行抑制，所以其實現原理也較復雜。同前幾種算法一樣，該算法也是構建在數學支持的基礎之上的。對于一階精度的歐拉公式有：

　　yi＋１＝yi＋h*K1
　　K1=f(xi,yi)

　　當用點xi處的斜率近似值K1與右端點xi＋１處的斜率K2的算術平均值作為平均斜率K*的近似值，那么就會得到二階精度的改進歐拉公式：

　　yi＋１＝yi＋h*( K1+ K2)/2
　　K1=f(xi,yi)
　　K2=f(xi+h,yi+h*K1)

　　依次類推，如果在區間[xi,xi+1]內多預估幾個點上的斜率值K1、K2、……Km，并用他們的加權平均數作為平均斜率K*的近似值，顯然能構造出具有很高精度的高階計算公式。經數學推導、求解，可以得出四階龍格－庫塔公式，也就是在工程中應用廣泛的經典龍格－庫塔算法：

　　yi＋１＝yi＋h*( K1+ 2*K2 +2*K3+ K4)/6
　　K1=f(xi,yi)
　　K2=f(xi+h/2,yi+h*K1/2)
　　K3=f(xi+h/2,yi+h*K2/2)
　　K4=f(xi+h,yi+h*K3)

　　下面的具體程序實現同改進的歐拉算法類似，只需作些必要的改動，下面將該算法的關鍵部分代碼清單列出：

…… for(float x=0;x<0.6;x+=0.1) { r=x+expf(-x); K1=x-y[i]+1; file://求K1 K2=(x+(float)(0.1/2))-(y[i]+K1*(float)(0.1/2))+1; file://求K2 K3=(x+(float)(0.1/2))-(y[i]+K2*(float)(0.1/2))+1; file://求K3 K4=(x+0.1)-(y[i]+K3*0.1)+1; file://求K4 y[i+1]=y[i]+(float)(0.1*(K1+2*K2+2*K3+K4)/6); file://求yi+1 r=fabs(r-y[i]); file://計算誤差 str.Format("y[%d]=%f r=%f\r\n",i,y[i],r); i++; msg+=str; } AfxMessageBox(msg); file://報告計算結果及誤差情況

　　從下表記錄的程序運行結果來看，在步長為0.1的情況下所計算出來的常微分方程的數值解是非常精確的，用浮點數進行運算時由近似所引入的誤差幾乎不會對計算結果產生影響，這樣高的精度是非常令人滿意的。無論是從計算速度上還是從計算精度要求上，都能非常好的滿足工程計算的需要。

xI（各分點）	yI （數值解）	y(xi) （真實值）	| y(xi)- yI | (誤差)
0.0	1.000000	1.000000	0.000000
0.1	1.004838	1.004837	0.000001
0.2	1.018731	1.018731	0.000000
0.3	1.040818	1.040818	0.000000
0.4	1.070320	1.070320	0.000000
0.5	1.106531	1.106531	0.000000

　　小結：本文針對工程計算中遇到的常微分方程初值問題的求解，根據實際情況引如了計算機作為輔助計算工具，并根據高等數學的有關知識將歐拉公式、改進的歐拉公式、經典龍格－庫塔方法等融入到計算機程序算法之中，充分利用了計算機的速度優勢，大大減輕了工程技術人員的勞動強度，同時也使計算結果更加可靠、準確。但在實際應用時，針對不同的工程函數選擇合適的求解方法需要有較高的要求，既要考慮到方法的簡易，又要減少計算量，同時又不能讓截斷誤差超出指定范圍。一般來說經典龍格－庫塔算法精確度高又利于計算機編程實現，穩定性也很好，可以考慮作為首選實現算法。

原文轉自：http://www.anti-gravitydesign.com

評論列表（網友評論僅供網友表達個人看法，并不表明本站同意其觀點或證實其描述）