1、微分方程數值解法【摘要】自然界與工程技術中的很多現象，可以歸結為微分方程定解問題。其中，常微分方程求解是微分方程的重要基礎內容。但是，對于許多的微分方程，往往很難得到甚至不存在精確的解析表達式，這時候，數值解提供了一個很好的解決思路。，針對于此，本文對常微分方程數值解法進行了簡單研究，主要討論了一些常用的數值解法，如歐拉法、改進的歐拉法、RungeKutta方法、Adams預估校正法以及勒讓德譜方法等，通過具體的算例，結合MATLAB求解畫圖，初步給出了一般常微分方程數值解法的求解過程。同時，通過對各種方法的誤差分析，讓大家對各種方法的特點和適用范圍有一個直觀的感受。【關鍵詞】常微分方程數值解

2、法MATLAB誤差分析引言在我國高校，微分方程數值解法作為對數學基礎知識要求較高且應用非常廣泛的一門課程，不僅在數學專業，其他的理工科專業的本科及研究生教育中開設這門課程.近四十年來，微分方程數值解法不論在理論上還是在方法上都獲得了很大的發展.同時，由于微分方程是描述物理、化學和生物現象的數學模型基礎，且它的一些最新應用已經擴展到經濟、金融預測、圖像處理及其他領域在實際應用中，通過相應的微分方程模型解決具體問題，采用數值方法求得方程的近似解，使具體問題迎刃而解。2歐拉法和改進的歐拉法2.1 歐拉法2.1.1 歐拉法介紹首先，我們考慮如下的一階常微分方程初值問題:y'=f(x,y)M%)

3、=y()(2-1)事實上，對于更復雜的常微分方程組或者高階常微分方程，只需要將X看做向量，(2-1)就成了一個一階常微分方程組，而高階常微分方程也可以通過降階化成一個一階常微分方程組。歐拉方法是解常微分方程初值問題最簡單最古老的一種數值方法，其基本思路就是把(2-1)中的導數項y'用差商逼近，從而將一個微分方程轉化為一個代數方程，以便求解。設在a,b】中取等距節點h,因為在節點xn點上，由(2-1)可得：y'(Xn)=f(Xn,y(Xn)(2-2)又由差商的定義可得:(2-3)所以有(2-4)y'(Xn)y(Xn1)-y(Xn)y(Xn1):y(Xn)hf(Xn，y(X

4、n)用y(Xk)的近似值yk(k=n,n+1)代入(2-4),則有計算yn書的歐拉公式yn1=ynhf(Xn,y(Xn)(2-5)2.1.2歐拉法誤差分析從歐拉公式中可以看出，右端的yn都是近似的，所以用它計算出來的yn書會有累計誤差，累計誤差比較復雜，為簡化分析，我們考慮局部截斷誤差，即認為Yn是精確的前提下來估計y(Xn書)-yn書，記為4*,泰勒展開有h23、近人1)/%)加7()O(h)(2-6)h2聯R(2-5),(2-6)即付與韋="2y'')+O(h3)=O(h2),根據數值算法精度的定義,如果一個數值方法的局部截斷誤差sn+=O(hp*)則稱這個算法具

5、有P階精度，所以，歐拉方法具有一階精度或者稱歐拉方法為一階方法。2.2改進的歐拉方法2.2.1改進的歐拉法介紹用數值積分離散化問題(1),兩邊做積分有:Xn1y(Xn1)-y(Xn)=xf(X,y(X)dXXn(2-7)對右端積分使用梯形積分公式可得:I'f(x,y(x)dx之hf(xn，y(xn)+f(Xn+，y(xn書)】Xn2(2-8)同歐拉方法，用y(xk)的近似值yk(k=n,n+1)代入(2-7),聯立(2-8)得到改進的歐拉方法計算公式：h，一yn1=yn-3(f(xn,yn)f(xn1,丫口/(2-9)一般來說，如果求解yn.時，算法右端不包含yn,稱為顯性計算公式，如

6、果包含，則求解時還需要解方程，這種稱為隱式計算公式。顯然公式(2-9)是一個隱式計算公式，事實上,改進的歐拉方法是用歐拉方法先求一個預測值右，再用這個預測值來計算yn卡,即：yn17nhf(xn,Yn)h%1fn萬(”4,外)f(41,%i)(2-10)2.2.2改進的歐拉法誤差分析同歐拉法誤差分析類似，用泰勒展開容易知道改進的歐拉方法具有二階精度，證明略2.3算例2.3.1 (一階常微分方程)求解初值問題y(0)=12xY=Yx0,1Y解析：在MATLA訃求解這個方程y=dsoke('Dy=y-2*x/y','y(0)=1','x')得y=(2

7、*x+1)A(1/2)它的解析解為y=J1+2x,下面我們分別用歐拉方法和改進的歐拉方法來求其數值解。歐拉方法：創建M文件euler1.m,內容如下：functionx,y=euler1(fun,x0,xfinal,y0,n)ifnargin<5,n=50;endh=(xfinal-x0)/n;fori=1:nx(i+1)=x(i)+h;y(i+1)=y(i)+h*feval(fun,x(i),y(i);End再定義函數方程組中的函數f1,創建f1.m文件，內容如下：functionf=f1(x,y)f=y-2*x/y在MATLA沖輸入：x,y=euler1('f1',0

8、,1,1,20)輸出f,x,y的值，將數值解跟精確解畫圖表示，輸入：plot(x,y,'r*-',x,sqrt(1+2*x),'g+-')xlabel('x');ylabel('y');title('y'=y-2x/y');legend('數值解,精確解')得到圖形,保存為euler1.fig,圖形如下：改進的歐拉方法：創建M文件eulerprove1.m,內容如下:functionx,y=eulerprove1(fun,x0,xfinal,y0,n)ifnargin<5,n=50;e

9、ndh=(xfinal-x0)/n;y(1)=y0,x(1)=x0;fori=1:nx(i+1)=x(i)+h;y1=y(i)+h*feval(fun,x(i),y(i)/2;y2=h*feval(fun,x(i+1),y1)/2;y(i+1)=y1+y2End返回f,x1,y1的值,在MATLAB勺command®口隼俞入:x,y1=eulerprove1('f1',0,1,1,20)plot(x,y1,'r*-',x,sqrt(1+2*x),'g+-')xlabel('x');ylabel('y');

10、title('y=y-2x/y'）;legend（'數值解','精確解'），將圖片保存為圖形eulerprove1.fig,為了便于比較兩種方法的誤差，將兩者的誤差作到同一個圖上，繼續輸入:plot(x,abs(y-sqrt(1+2*x),'y*-',x,abs(y1-sqrt(1+2*x),'g+-');xlabel('x');ylabel('y');title('誤差曲線');legend('歐拉方法','改進的歐拉方法')將圖片保

11、存為error1.fig,圖形如下：從該圖形來看，改進的歐拉方法與歐拉方法誤差接近，歐拉方法誤差稍微大些，將x的取值擴寬，n取值增大時，可以發現改進的歐拉方法相比歐拉方法有更高的精度2.3.2 (高階微分方程)對于二階常微分方程JI_x=_xx'(0)=2x(0)=1(xw0,1),求數值解解析：先算出其解析解，在MATLA訃輸入：y=dsolve('D2x=-x','x(0)=1','Dx(0)=2')得到解為:Y=2*sin(t)+cos(t),前面已經分別給出過歐拉方法和改進的歐拉方法的算例跟誤差比較,這里我們就用精度更高的改進歐拉

12、法進行數值求解。改進的歐拉方法：先換元，令x'=y,則原方程可以轉化為y'7«x'=yx(0)=1y(0)=2(x0,1),現在，二階常微分方程轉化為了一個一階常微分方程組，同2.3.2的方法，建立M文件eulerprove3.m,內容如下:functiont,x,y=eulerprove2(t0,tfinal,x0,y0,n)f1=inline('y');f2=inline('-x')ifnargin<5,n=50;endh=(tfinal-t0)/n;t(1)=t0,x(1)=x0;y(1)=y0;fori=1:nt(

13、i+1)=t(i)+h;x1=x(i)+h*feval(f1,y(i);y1=y(i)+h*feval(f2,x(i);x2=x(i)+h*feval(f1,y1);y2=y(i)+h*feval(f2,x1);x(i+1)=(x1+x2)/2;y(i+1)=(y1+y2)/2;end在command®口隼入t,x,y=eulerprove3(0,1,1,2,10),得到t,x,y的值，其中x就是我們要求的數值解，作圖，輸入：plot(t,x,'r*-',t,2*sin(t)+cos(t),'b+-');xlabel('x');ylab

14、el('y');legend('數值解,精確解')，將圖形保存為eulerprove3.fig,圖形如下：上面，我們已經通過例子看出，改進的歐拉法相比于歐拉法，在每一個節點處的誤差值更下，下面，我們來討論節點的多少(步長大小)對誤差的影響，創建err03r.m文件，內容如下：functionN,Y=error3(n0,nfinal)N(1)=n0,m=fix(nfinal-n0)/4);|fori=1:mN(i+1)=N(i)+4;t,x1,y1=eulerprove3(0,1,1,2,N(i);|Y(i)=log10(max(x1-2*sin(t)-cos(t

15、);t,x2,y1=eulerprove3(0,1,1,2,N(i);|Y(i+1)=log10(max(x2-2*sin(t)-cos(t);end輸入N,Y=error3(4,100),返回節點個數值和Y值，Y代表在N個節點時，數值解與精確解差的絕對值的最大值的對數(10為底)。:plot(N,Y,'d-.'),xlabel('N'),ylabel('log_1_0max|error|)title('誤差曲線'),將圖片保存為error3.fig,圖形如下：口一loMBVEQdrn-所以，為了得到盡可從這個圖的誤差曲線，隨著插值節點個

16、數增多，誤差也是越來越小的,能準確的解，使用改進的歐拉法時插值節點數應當越大越好。但是，隨著節點個數的增多,計算量也會隨之增大，所以，在實際應用時，應當結合具體要求靈活選取。3龍格-庫塔法3.1 龍格-庫塔法與泰勒展開3.1.1 泰勒(Taylor)展開首先，考慮考慮微分方程y'=f(x,y)(3-1)引入算子符號°+f。,二x二y(3-2)于是有(3-3)y''=fxfyy'=fxffy二二ff；x；:yy'"-f(ff)-f2f,；x:y；x；y；xfyy(m)fTm'fFxjy設不帶括號的f及其各階微商都在(xn,y(x

17、n)處取值，于是有泰勒展開式:y(xn1)7函)hfh2L2!；:x一3二h；2fff2fjy3!；:x；:y(3-4)h4r4!Fx5O(h5)7%)hfffyfxx2ffxyf2fyyfxfyffy23!h423'fxxx3ffxxy3ffxyy'ffyyyfyfxx3fxfxy5ffyfxy4!_2_23ffxfyy4ffyfyyfxfy“35ffyO(h)3.1.2龍格-庫塔(R-K)方法介紹龍格-庫塔方法的基本思路是想辦法計算f(x,y)在某些點上的函數值，然后對這些函數值做數值線性組合，構造出一個近似的計算公式；再把近似的計算公式和解的泰勒展開式相比較，使得前面的若

18、干項相吻合，從而達到較高的精度，一般的顯示R-K方法的形式如下:ryn由=yn&iki,iT，k1=hf(xn,yn),ik=hf(xn+sh+E0/)，(2«i«r)jT(3-5)當式(3-5)的局部截斷誤差達到O(hp+)時稱公式(35)為p階r段R-K方法。龍格庫塔方法是應用最廣的求解常微分方程初值問題的單步法，下面我們給出幾個推導公式。3.2 龍格-庫塔法公式與ode函數3.2.1 二階二段R-K法公式推導由式子（3-5）,二階二段的R-K方法可以寫成yn1=Yn1k12k2kl=hf（Xn，Yn）,K=hf（Xn+uh,yn+P21K）,（3-6）不帶括號

19、的f及其各階微商都在（xn,y（xn）處取值，由泰勒展開式（3-4）及截斷誤差的定義可知：n1=丫函1)-y(Xn)-”hf(Xn,Y(Xn)-2hf(4二工2。y(Xn)-21hf(Xn,y(Xn)h2=y(Xn)hf萬ffy)-y(Xn)-1hf-0(f二21%hfy)O(h3)1o1Q=(1一1一，2)hf4-,2：2)hfX(-2：21)hffyO(h)(3-7)由于上式是二階二段的R-K法，即必須有8n書=O（h3）,所以E1+切2=1，920f2=1/2,、。2021=1/2,（3-8）式（3-8）有四個未知元，三個方程，故有無窮組解。若令立2=1,由式子（3-9）可解得以=

20、69;2=1/2,221=1,于1. yn+=yn+一左十一卜2,22<k1=hf(Xn,yn)k2=hf(Xn+h,yn+k)這個公式稱為預估-校正公式。常用的二階公式還有中點公式（取口2=1/2）,這里不再列舉。3.2.2其他常見公式和ode函數從二階二段R-K公式的推導可以看出，二段方法每一步需要計算兩次函數值，而它也只能達到2階精度，如果我們提高函數的計算次數，就可以得到精度更高的計算公式，高階的R-K法的推導與二階方法的推導完全相同，只是隨著階數的增高計算量會逐漸增大。下面給出幾個常用的計算公式。（1）庫塔公式（3階3段）1yn書=V。+(ki+4k2+k3)6，ki=hf(x

21、n,y。)k2=hf(Xn+h/2,yn+ki/2)k3=hf(Xnh,Vn-ki2k2)(2)標準四階R-K公式11”一、y5=yn+(ki+2k2+2k3+k4)6ki=hf(Xn,yn)*2=hf(Xn+h/2,yn+ki/2)k3=hf(Xn+h/2,yn+k2/2)k4=hf(Xn+h,yn+k3)上述三階公式具有三階精度，四階公式具有四階精度，要驗證這一點，我們觀察如下泰勒展開:nf(Xh,yh-裨三k.)if(X，y)i=e(ni)!一一nH(hk)f(x,y)二x二y22尸i3尸=fhfXkfy藥(h2fXX2hkfXyk2fyy)(h3fXXX3h2kfXXy3hk2fXyy

22、k3fyyy).(3-9)其中x,y落在平面上連接（x,y）和（x+h,y+k）的線段上。只需要按照式（3.9）將k泰勒展開，代入截斷誤差的計算公式中即可證明。（詳細證明見參考文獻2P253）由于R-K法是基于泰勒展開的數值方法，所以，如果解具有較好的光滑性，則能夠得到較好的計算精度，反之，則可能四階R-K法的數值精度還不如低階的數值方法，這個需要根據具體情況自行選擇數值方法。一般而言，R-K法在求解范圍較大而精度要求較高時是比較好的方法，四階的R-K法也可以用作阿達姆斯預估校正法的前幾步計算，以保持四階精度（下一節阿達姆斯預估校正法中也提到了這一點）。在MATLABK專門提供了求解微分方程的

23、ode函數，如ode45,ode23,ode113,ode15s,ode23s等等。ode求解器分為變步長和固定步長兩種求解模式，其中，ode45是采用R-K法的變步長求解器，和它一樣的還有ode23。目前，ode45是使用最多的求解函數，主要用于求解非剛性常微分方程。(如果求解時長時間沒有響應，則方程是剛性的，可以換用ode23求解)，ode函數的調用方式大都類似，以ode45為例，常用的語法格式有：T,Y=ode45(odefun,tspan,y0);T,Y=ode45(odefun,tspan,y0,options);sol=ode45(odefun,t0tf,y0.),其中odefun

24、是函數句柄，可以是函數文件名或內聯函數名，tspan是求解區間t0tf或者一系列散點t0,t1,.,tf,y0是初始值向量，T返回列向量的時間點，Y返回對應T的求解列向量，options是求解參數設置，可以用odeset在計算前設定誤差，輸出參數事件等。3.3算例用自帶ODE45求解高階常微分方程組L3x''+x(x2+y2產=03將常微分方程組初值問題y''+y(x2+y2)W)=0轉化為一階常微分方程組并用x(0)=0.5,x'(0)=0.75y(0)=0.25,y'(0)=1.0ode45求解x(1),y(1).解析：首先我們嘗試用MATLA球出該方程組的解析解，在command®口輸入：x,y=dsolve('D2x+x*(xA2+yA2)A(-3/2)=0','D2y+y*(xA2+yA2)A(-3/2)=0','x(0)=0.5','y(0)=0.25','Dx(0)=0.75','Dy(0)=1'),運行大概五分鐘，提示：Warning:Explicitsolutioncouldnotbefound.x=emptysym,y=,產生該問題的原因可能在于該非線性方程無法用solve求解，或該方程的解析解不存在。下