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

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

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

3、=y()(2-1)事實(shí)上，對(duì)于更復(fù)雜的常微分方程組或者高階常微分方程，只需要將X看做向量，(2-1)就成了一個(gè)一階常微分方程組，而高階常微分方程也可以通過降階化成一個(gè)一階常微分方程組。歐拉方法是解常微分方程初值問題最簡單最古老的一種數(shù)值方法，其基本思路就是把(2-1)中的導(dǎo)數(shù)項(xiàng)y'用差商逼近，從而將一個(gè)微分方程轉(zhuǎn)化為一個(gè)代數(shù)方程，以便求解。設(shè)在a,b】中取等距節(jié)點(diǎn)h,因?yàn)樵诠?jié)點(diǎn)xn點(diǎn)上，由(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),則有計(jì)算yn書的歐拉公式y(tǒng)n1=ynhf(Xn,y(Xn)(2-5)2.1.2歐拉法誤差分析從歐拉公式中可以看出，右端的yn都是近似的，所以用它計(jì)算出來的yn書會(huì)有累計(jì)誤差，累計(jì)誤差比較復(fù)雜，為簡化分析，我們考慮局部截?cái)嗾`差，即認(rèn)為Yn是精確的前提下來估計(jì)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ù)值算法精度的定義,如果一個(gè)數(shù)值方法的局部截?cái)嗾`差sn+=O(hp*)則稱這個(gè)算法具

5、有P階精度，所以，歐拉方法具有一階精度或者稱歐拉方法為一階方法。2.2改進(jìn)的歐拉方法2.2.1改進(jìn)的歐拉法介紹用數(shù)值積分離散化問題(1),兩邊做積分有:Xn1y(Xn1)-y(Xn)=xf(X,y(X)dXXn(2-7)對(duì)右端積分使用梯形積分公式可得: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)得到改進(jìn)的歐拉方法計(jì)算公式：h，一yn1=yn-3(f(xn,yn)f(xn1,丫口/(2-9)一般來說，如果求解yn.時(shí)，算法右端不包含yn,稱為顯性計(jì)算公式，如

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

7、*x+1)A(1/2)它的解析解為y=J1+2x,下面我們分別用歐拉方法和改進(jìn)的歐拉方法來求其數(shù)值解。歐拉方法：創(chuàng)建M文件euler1.m,內(nèi)容如下：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文件，內(nèi)容如下：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,圖形如下：改進(jìn)的歐拉方法：創(chuàng)建M文件eulerprove1.m,內(nèi)容如下: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,為了便于比較兩種方法的誤差，將兩者的誤差作到同一個(gè)圖上，繼續(xù)輸入:plot(x,abs(y-sqrt(1+2*x),'y*-',x,abs(y1-sqrt(1+2*x),'g+-');xlabel('x');ylabel('y');title('誤差曲線');legend('歐拉方法','改進(jìn)的歐拉方法')將圖片保

11、存為error1.fig,圖形如下：從該圖形來看，改進(jìn)的歐拉方法與歐拉方法誤差接近，歐拉方法誤差稍微大些，將x的取值擴(kuò)寬，n取值增大時(shí)，可以發(fā)現(xiàn)改進(jìn)的歐拉方法相比歐拉方法有更高的精度2.3.2 (高階微分方程)對(duì)于二階常微分方程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)分別給出過歐拉方法和改進(jìn)的歐拉方法的算例跟誤差比較,這里我們就用精度更高的改進(jìn)歐拉

12、法進(jìn)行數(shù)值求解。改進(jìn)的歐拉方法：先換元，令x'=y,則原方程可以轉(zhuǎn)化為y'7«x'=yx(0)=1y(0)=2(x0,1),現(xiàn)在，二階常微分方程轉(zhuǎn)化為了一個(gè)一階常微分方程組，同2.3.2的方法，建立M文件eulerprove3.m,內(nèi)容如下: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)通過例子看出，改進(jìn)的歐拉法相比于歐拉法，在每一個(gè)節(jié)點(diǎn)處的誤差值更下，下面，我們來討論節(jié)點(diǎn)的多少(步長大小)對(duì)誤差的影響，創(chuàng)建err03r.m文件，內(nèi)容如下：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é)點(diǎn)個(gè)數(shù)值和Y值，Y代表在N個(gè)節(jié)點(diǎn)時(shí)，數(shù)值解與精確解差的絕對(duì)值的最大值的對(duì)數(shù)(10為底)。:plot(N,Y,'d-.'),xlabel('N'),ylabel('log_1_0max|error|)title('誤差曲線'),將圖片保存為error3.fig,圖形如下：口一loMBVEQdrn-所以，為了得到盡可從這個(gè)圖的誤差曲線，隨著插值節(jié)點(diǎn)個(gè)

16、數(shù)增多，誤差也是越來越小的,能準(zhǔn)確的解，使用改進(jìn)的歐拉法時(shí)插值節(jié)點(diǎn)數(shù)應(yīng)當(dāng)越大越好。但是，隨著節(jié)點(diǎn)個(gè)數(shù)的增多,計(jì)算量也會(huì)隨之增大，所以，在實(shí)際應(yīng)用時(shí)，應(yīng)當(dāng)結(jié)合具體要求靈活選取。3龍格-庫塔法3.1 龍格-庫塔法與泰勒展開3.1.1 泰勒(Taylor)展開首先，考慮考慮微分方程y'=f(x,y)(3-1)引入算子符號(hào)°+f。,二x二y(3-2)于是有(3-3)y''=fxfyy'=fxffy二二ff；x；:yy'"-f(ff)-f2f,；x:y；x；y；xfyy(m)fTm'fFxjy設(shè)不帶括號(hào)的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)方法介紹龍格-庫塔方法的基本思路是想辦法計(jì)算f(x,y)在某些點(diǎn)上的函數(shù)值，然后對(duì)這些函數(shù)值做數(shù)值線性組合，構(gòu)造出一個(gè)近似的計(jì)算公式；再把近似的計(jì)算公式和解的泰勒展開式相比較，使得前面的若

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

19、的f及其各階微商都在（xn,y（xn）處取值，由泰勒展開式（3-4）及截?cái)嗾`差的定義可知：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）有四個(gè)未知元，三個(gè)方程，故有無窮組解。若令立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)這個(gè)公式稱為預(yù)估-校正公式。常用的二階公式還有中點(diǎn)公式（取口2=1/2）,這里不再列舉。3.2.2其他常見公式和ode函數(shù)從二階二段R-K公式的推導(dǎo)可以看出，二段方法每一步需要計(jì)算兩次函數(shù)值，而它也只能達(dá)到2階精度，如果我們提高函數(shù)的計(jì)算次數(shù)，就可以得到精度更高的計(jì)算公式，高階的R-K法的推導(dǎo)與二階方法的推導(dǎo)完全相同，只是隨著階數(shù)的增高計(jì)算量會(huì)逐漸增大。下面給出幾個(gè)常用的計(jì)算公式。（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)標(biāo)準(zhǔn)四階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)上述三階公式具有三階精度，四階公式具有四階精度，要驗(yàn)證這一點(diǎn)，我們觀察如下泰勒展開: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泰勒展開，代入截?cái)嗾`差的計(jì)算公式中即可證明。（詳細(xì)證明見參考文獻(xiàn)2P253）由于R-K法是基于泰勒展開的數(shù)值方法，所以，如果解具有較好的光滑性，則能夠得到較好的計(jì)算精度，反之，則可能四階R-K法的數(shù)值精度還不如低階的數(shù)值方法，這個(gè)需要根據(jù)具體情況自行選擇數(shù)值方法。一般而言，R-K法在求解范圍較大而精度要求較高時(shí)是比較好的方法，四階的R-K法也可以用作阿達(dá)姆斯預(yù)估校正法的前幾步計(jì)算，以保持四階精度（下一節(jié)阿達(dá)姆斯預(yù)估校正法中也提到了這一點(diǎn)）。在MATLABK專門提供了求解微分方程的

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

24、是函數(shù)句柄，可以是函數(shù)文件名或內(nèi)聯(lián)函數(shù)名，tspan是求解區(qū)間t0tf或者一系列散點(diǎn)t0,t1,.,tf,y0是初始值向量，T返回列向量的時(shí)間點(diǎn)，Y返回對(duì)應(yīng)T的求解列向量，options是求解參數(shù)設(shè)置，可以用odeset在計(jì)算前設(shè)定誤差，輸出參數(shù)事件等。3.3算例用自帶ODE45求解高階常微分方程組L3x''+x(x2+y2產(chǎn)=03將常微分方程組初值問題y''+y(x2+y2)W)=0轉(zhuǎn)化為一階常微分方程組并用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'),運(yùn)行大概五分鐘，提示：Warning:Explicitsolutioncouldnotbefound.x=emptysym,y=,產(chǎn)生該問題的原因可能在于該非線性方程無法用solve求解，或該方程的解析解不存在。下