微分方程數(shù)值解法

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

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

3、=y()(2-1)事實上，對于更復雜的常微分方程組或者高階常微分方程，只需要將X看做向量，(2-1)就成了一個一階常微分方程組，而高階常微分方程也可以通過降階化成一個一階常微分方程組。歐拉方法是解常微分方程初值問題最簡單最古老的一種數(shù)值方法，其基本思路就是把(2-1)中的導數(shù)項y'用差商逼近，從而將一個微分方程轉化為一個代數(shù)方程，以便求解。設在a,b】中取等距節(jié)點h,因為在節(jié)點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書的歐拉公式y(tǒng)n1=ynhf(Xn,y(Xn)(2-5)2.1.2歐拉法誤差分析從歐拉公式中可以看出，右端的yn都是近似的，所以用它計算出來的yn書會有累計誤差，累計誤差比較復雜，為簡化分析，我們考慮局部截斷誤差，即認為Yn是精確的前提下來估計y(Xn書)-yn書，記為4*,泰勒展開有h23、近人1)/%)加7()O(h)(2-6)h2聯(lián)R(2-5),(2-6)即付與韋="2y'')+O(h3)=O(h2),根據(jù)數(shù)值算法精度的定義,如果一個數(shù)值方法的局部截斷誤差sn+=O(hp*)則稱這個算法具

5、有P階精度，所以，歐拉方法具有一階精度或者稱歐拉方法為一階方法。2.2改進的歐拉方法2.2.1改進的歐拉法介紹用數(shù)值積分離散化問題(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),聯(lián)立(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,下面我們分別用歐拉方法和改進的歐拉方法來求其數(shù)值解。歐拉方法：創(chuàng)建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再定義函數(shù)方程組中的函數(shù)f1,創(chuàng)建f1.m文件，內容如下：functionf=f1(x,y)f=y-2*x/y在MATLA沖輸入：x,y=euler1('f1',0

8、,1,1,20)輸出f,x,y的值，將數(shù)值解跟精確解畫圖表示，輸入：plot(x,y,'r*-',x,sqrt(1+2*x),'g+-')xlabel('x');ylabel('y');title('y'=y-2x/y');legend('數(shù)值解,精確解')得到圖形,保存為euler1.fig,圖形如下：改進的歐拉方法：創(chuàng)建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（'數(shù)值解','精確解'），將圖片保存為圖形eulerprove1.fig,為了便于比較兩種方法的誤差，將兩者的誤差作到同一個圖上，繼續(xù)輸入: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取值增大時，可以發(fā)現(xiàn)改進的歐拉方法相比歐拉方法有更高的精度2.3.2 (高階微分方程)對于二階常微分方程JI_x=_xx'(0)=2x(0)=1(xw0,1),求數(shù)值解解析：先算出其解析解，在MATLA訃輸入：y=dsolve('D2x=-x','x(0)=1','Dx(0)=2')得到解為:Y=2*sin(t)+cos(t),前面已經(jīng)分別給出過歐拉方法和改進的歐拉方法的算例跟誤差比較,這里我們就用精度更高的改進歐拉

12、法進行數(shù)值求解。改進的歐拉方法：先換元，令x'=y,則原方程可以轉化為y'7«x'=yx(0)=1y(0)=2(x0,1),現(xiàn)在，二階常微分方程轉化為了一個一階常微分方程組，同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就是我們要求的數(shù)值解，作圖，輸入：plot(t,x,'r*-',t,2*sin(t)+cos(t),'b+-');xlabel('x');ylab

14、el('y');legend('數(shù)值解,精確解')，將圖形保存為eulerprove3.fig,圖形如下：上面，我們已經(jīng)通過例子看出，改進的歐拉法相比于歐拉法，在每一個節(jié)點處的誤差值更下，下面，我們來討論節(jié)點的多少(步長大小)對誤差的影響，創(chuàng)建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),返回節(jié)點個數(shù)值和Y值，Y代表在N個節(jié)點時，數(shù)值解與精確解差的絕對值的最大值的對數(shù)(10為底)。:plot(N,Y,'d-.'),xlabel('N'),ylabel('log_1_0max|error|)title('誤差曲線'),將圖片保存為error3.fig,圖形如下：口一loMBVEQdrn-所以，為了得到盡可從這個圖的誤差曲線，隨著插值節(jié)點個

16、數(shù)增多，誤差也是越來越小的,能準確的解，使用改進的歐拉法時插值節(jié)點數(shù)應當越大越好。但是，隨著節(jié)點個數(shù)的增多,計算量也會隨之增大，所以，在實際應用時，應當結合具體要求靈活選取。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)在某些點上的函數(shù)值，然后對這些函數(shù)值做數(shù)值線性組合，構造出一個近似的計算公式；再把近似的計算公式和解的泰勒展開式相比較，使得前面的若

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函數(shù)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函數(shù)從二階二段R-K公式的推導可以看出，二段方法每一步需要計算兩次函數(shù)值，而它也只能達到2階精度，如果我們提高函數(shù)的計算次數(shù)，就可以得到精度更高的計算公式，高階的R-K法的推導與二階方法的推導完全相同，只是隨著階數(shù)的增高計算量會逐漸增大。下面給出幾個常用的計算公式。（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法是基于泰勒展開的數(shù)值方法，所以，如果解具有較好的光滑性，則能夠得到較好的計算精度，反之，則可能四階R-K法的數(shù)值精度還不如低階的數(shù)值方法，這個需要根據(jù)具體情況自行選擇數(shù)值方法。一般而言，R-K法在求解范圍較大而精度要求較高時是比較好的方法，四階的R-K法也可以用作阿達姆斯預估校正法的前幾步計算，以保持四階精度（下一節(jié)阿達姆斯預估校正法中也提到了這一點）。在MATLABK專門提供了求解微分方程的

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

24、是函數(shù)句柄，可以是函數(shù)文件名或內聯(lián)函數(shù)名，tspan是求解區(qū)間t0tf或者一系列散點t0,t1,.,tf,y0是初始值向量，T返回列向量的時間點，Y返回對應T的求解列向量，options是求解參數(shù)設置，可以用odeset在計算前設定誤差，輸出參數(shù)事件等。3.3算例用自帶ODE45求解高階常微分方程組L3x''+x(x2+y2產(chǎn)=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=,產(chǎn)生該問題的原因可能在于該非線性方程無法用solve求解，或該方程的解析解不存在。下