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