1、精選優(yōu)質(zhì)文檔-傾情為你奉上常微分方程的幾種數(shù)值解法數(shù)學(xué)與應(yīng)用數(shù)學(xué) 肖振華 指導(dǎo)教師 張秀艷【摘要】自然界與工程技術(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ì)各種方法的誤差分析，讓

2、大家對(duì)各種方法的特點(diǎn)和適用范圍有一個(gè)直觀的感受。【關(guān)鍵詞】 常微分方程 數(shù)值解法 MATLAB 誤差分析【Abstract】 Many phenomena in nature and engineering can be attributed to the definite solution of the problem for differential equations. Among them, the ordinary differential equation solving is an important foundation for the content of the diffe

3、rential equations. However, many of the differential equations are often difficult to obtain accurate analytical expression .At this time, the numerical solution provides a good idea. For the Numerical Solution of Ordinary Differential Equations in this article, we focuses on some commonly used nume

4、rical solution, such as the Euler method, improved Euler method, Runge-Kutta method, Adams predictor corrector method as well as newer spectral methods. Through specific examples, combined with MATLAB solving and drawing, we initially know the solution process of general numerical solution of ordina

5、ry differential equations . At the same time, according to the error analysis of various methods , everyone has an intuitive feel of the characteristics and scope of the various methods.【Keywords】Ordinary Differential Equations Numerical Solution MATLAB error analysis目錄1 前言1.1 常微分方程概述方程是一個(gè)在數(shù)學(xué)中非常熟悉的名

6、詞，在初等數(shù)學(xué)里，我們將我們要研究的問題作為一個(gè)或幾個(gè)未知量，通過觀察事物的規(guī)律，得出這些未知量與已知量之間的等式關(guān)系，這樣就得到了一個(gè)簡單的方程或方程組當(dāng)然，這只是一個(gè)很淺顯粗略的定義。在數(shù)學(xué)上，物質(zhì)的運(yùn)動(dòng)和變化規(guī)律是通過函數(shù)關(guān)系來表示的，在一些復(fù)雜的現(xiàn)象中，比如火箭的運(yùn)動(dòng)軌跡，我們要求的未知量就變成了滿足特定條件的一個(gè)或一些未知函數(shù)（不再是簡單的數(shù)值）。有時(shí)候，我們需要用到導(dǎo)數(shù)很微分的知識(shí)，即這些未知函數(shù)的導(dǎo)數(shù)與自變量滿足某種關(guān)系，這種方程，我們稱之為微分方程，即可以這樣定義，凡是含有參數(shù)，未知函數(shù)和未知函數(shù)導(dǎo)數(shù)（或微分）的方程稱之為微分方程。未知函數(shù)是一元函數(shù)的微分方程稱為常微分方程（O

7、DE），未知函數(shù)是多元函數(shù)的微分方程稱為偏微分方程（PDE），我們要研究的，就是常微分方程。1.2 常微分方程解與數(shù)值解法微分方程的解，就是找出一個(gè)代入方程使之成為恒等式的函數(shù)。.若微分方程的解中含有任意的個(gè)數(shù)與方程的階數(shù)相同，且任意常數(shù)之間不能合并，則稱此解為該方程的通解(或一般解).當(dāng)通解中的各任意常數(shù)都取特定值時(shí)所得到的解，稱為方程的特解。在實(shí)際問題中，這些函數(shù)往往還滿足一些特定條件，這個(gè)稱之為定解問題，如常微分方程的初值問題和邊值問題，滿足條件的解就是特解。在理想的狀況下，我們解這個(gè)方程如果能夠得出通解的解析表達(dá)式 ，那么代入數(shù)值很容易得到問題所要求的特解，但事實(shí)上，很多常微分方程得不

8、到通解的解析表達(dá)式，或者表達(dá)式過于復(fù)雜。而且，常微分方程的特解是否存在，有幾個(gè)解，牽涉到解的存在唯一性定理。在實(shí)際應(yīng)用中，我們對(duì)解的精度要求往往并沒有那么嚴(yán)格。在這種情況下，數(shù)值解應(yīng)運(yùn)而生，通俗的講，數(shù)值解就是所求函數(shù)在特定點(diǎn)的函數(shù)值，通過函數(shù)在特定點(diǎn)的取值，我們可以反過來觀察這個(gè)函數(shù)的某些性質(zhì)，所以，求常微分方程的數(shù)值解是很有意義的，事實(shí)上，解常微分方程很大一部分就是解數(shù)值解。 下面我們要討論幾種常見的常微分方程數(shù)值解法，目的是展示常微分方程數(shù)值解法的一般過程。以下章節(jié)，由于本人水平有限，只對(duì)數(shù)值解法的思路和過程進(jìn)行討論，而對(duì)于解的存在唯一性，方法精度的具體證明，不做過多討論。一般而言，構(gòu)造

9、數(shù)值解法都是基于兩個(gè)想法：一種是基于數(shù)值積分的構(gòu)造方法，一種是基于泰勒展開的構(gòu)造法。下面我們要討論的是幾種常見的數(shù)值解法，如歐拉法、改進(jìn)的歐拉法、龍格庫塔法、阿達(dá)姆斯預(yù)估校正法以及基于勒讓德多項(xiàng)式的譜方法。這些方法既有基于數(shù)值積分的構(gòu)造方法，如歐拉法，也有基于泰勒展開的構(gòu)造法，如龍格庫塔法。下面的章節(jié)，我們會(huì)一一介紹這些方法的思路，并通過一些簡單的算例，讓大家了解這些解法的一般過程。為了讓大家看到各種數(shù)值解法的精確性和有效性，我們可能會(huì)盡量選取一些通解存在的算例來求數(shù)值解，這樣做并不意味著我們實(shí)際中求數(shù)值解時(shí)要先求出精確解，而只是出于對(duì)方法對(duì)比的需要，望不致引起混淆。2 歐拉法和改進(jìn)的歐拉法2

10、.1 歐拉法 2.1.1 歐拉法介紹首先，我們考慮如下的一階常微分方程初值問題 （2-1）事實(shí)上，對(duì)于更復(fù)雜的常微分方程組或者高階常微分方程，只需要將看做向量，（2-1）就成了一個(gè)一階常微分方程組，而高階常微分方程也可以通過降階化成一個(gè)一階常微分方程組。歐拉方法是解常微分方程初值問題最簡單最古老的一種數(shù)值方法，其基本思路就是把（2-1）中的導(dǎo)數(shù)項(xiàng)用差商逼近，從而將一個(gè)微分方程轉(zhuǎn)化為一個(gè)代數(shù)方程，以便求解。設(shè)在中取等距節(jié)點(diǎn)，因?yàn)樵诠?jié)點(diǎn)點(diǎn)上，由（2-1）可得： ， （2-2） 又由差商的定義可得： （2-3） 所以有 （2-4）用的近似值代入（2-4），則有計(jì)算的歐拉公式 （2-5）2.1.2歐拉

11、法誤差分析從歐拉公式中可以看出，右端的都是近似的，所以用它計(jì)算出來的會(huì)有累計(jì)誤差，累計(jì)誤差比較復(fù)雜，為簡化分析，我們考慮局部截?cái)嗾`差，即認(rèn)為是精確的前提下來估計(jì)，記為，泰勒展開有 （2-6）聯(lián)立（2-5），（2-6）即得=+=，根據(jù)數(shù)值算法精度的定義，如果一個(gè)數(shù)值方法的局部截?cái)嗾`差=則稱這個(gè)算法具有P階精度，所以，歐拉方法具有一階精度或者稱歐拉方法為一階方法。2.2 改進(jìn)的歐拉方法2.2.1 改進(jìn)的歐拉法介紹用數(shù)值積分離散化問題（1），兩邊做積分有： （2-7）對(duì)右端積分使用梯形積分公式可得： （2-8）同歐拉方法，用的近似值代入（2-7），聯(lián)立（2-8）得到改進(jìn)的歐拉方法計(jì)算公式： （2-9

12、）一般來說，如果求解時(shí)，算法右端不包含，稱為顯性計(jì)算公式，如果包含，則求解時(shí)還需要解方程，這種稱為隱式計(jì)算公式。顯然公式（2-9）是一個(gè)隱式計(jì)算公式，事實(shí)上，改進(jìn)的歐拉方法是用歐拉方法先求一個(gè)預(yù)測(cè)值，再用這個(gè)預(yù)測(cè)值來計(jì)算，即： （2-10）2.2.2改進(jìn)的歐拉法誤差分析同歐拉法誤差分析類似，用泰勒展開容易知道改進(jìn)的歐拉方法具有二階精度，證明略。2.3 算例2.3.1（一階常微分方程）求解初值問題 解析：在MATLAB中求解這個(gè)方程y=dsolve('Dy=y-2*x/y','y(0)=1','x') 得y =(2*x+1)(1/2) 它的解析解為

13、，下面我們分別用歐拉方法和改進(jìn)的歐拉方法來求其數(shù)值解。歐拉方法：創(chuàng)建M文件euler1.m,內(nèi)容如下：function x,y=euler1(fun,x0,xfinal,y0,n)if nargin<5,n=50;endh=(xfinal-x0)/n;x(1)=x0;y(1)=y0;for i=1:n x(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)容如下：function f=f1(x,y)f=y-2*x/y在MATLAB中輸入：x,y=euler1('f1',0,

14、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)容如下：function x,y=eulerprove1(fun,x0,xfinal,y0,n)if nargin<

15、5,n=50;endh=(xfinal-x0)/n;y(1)=y0,x(1)=x0;for i=1:n x(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在MATLAB的command窗口輸入: x,y1=eulerprove1('f1',0,1,1,20)返回f,x1,y1的值，作圖，輸入：plot(x,y1,'r*-',x,sqrt(1+2*x),'g+-')； xlabel('x');ylab

16、el('y');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('歐拉方法',

17、9;改進(jìn)的歐拉方法')將圖片保存為error1.fig，圖形如下：從該圖形來看，改進(jìn)的歐拉方法與歐拉方法誤差接近，歐拉方法誤差稍微大些，將x的取值擴(kuò)寬，n取值增大時(shí)，可以發(fā)現(xiàn)改進(jìn)的歐拉方法相比歐拉方法有更高的精度。2.3.2（一階常微分方程組）初值問題： ，先求出精確解，再用數(shù)值方法求解。解析：在MATLAB中輸入：x,y=dsolve('Dx=x+4*y-exp(t)','Dy=x+y+2*exp(t)','x(0)=4','y(0)=5/4','t')得：x =4*exp(3*t)+2*exp(-t)-

18、2*exp(t)，y =2*exp(3*t)-exp(-t)+1/4*exp(t)，下面繼續(xù)用歐拉方法和改進(jìn)的歐拉方法求數(shù)值解。歐拉方法：創(chuàng)建M文件euler2.m，內(nèi)容如下:function t,x,y=euler2(t0,tfinal,x0,y0,n)f1=inline('x+4*y-exp(t)');f2=inline('x+y+2*exp(t)')if nargin<5,n=50;endh=(tfinal-t0)/n;t(1)=t0,x(1)=x0;y(1)=y0;for i=1:n t(i+1)=t(i)+h; x(i+1)=x(i)+h*fev

19、al(f1,t(i),x(i),y(i); y(i+1)=y(i)+h*feval(f2,t(i),x(i),y(i);End在command窗口輸入t,x,y=euler2(0,1,4,5/4,100)，得到t，x,y的值，作圖，輸入命令如下： plot(t,x,'r*-',t,y,'g*-',t,4*exp(3*t)+2*exp(-t)-2*exp(t),'b+-',t,2*exp(3*t)-exp(-t)+1/4*exp(t),'y+-')；legend('x數(shù)值解','y數(shù)值解','

20、x精確解','y精確解')； xlabel('t'),ylabel('x(y)')，將圖形保存為euler2.fig圖形如下：改進(jìn)的歐拉方法：創(chuàng)建M文件eulerprove2.m文件，內(nèi)容如下：function t,x,y=eulerprove2(t0,tfinal,x0,y0,n)f1=inline('x+4*y-exp(t)');f2=inline('x+y+2*exp(t)')if nargin<5,n=50;endh=(tfinal-t0)/n;t(1)=t0,x(1)=x0;y(1)=y0;

21、for i=1:n t(i+1)=t(i)+h; x1=x(i)+h*feval(f1,t(i),x(i),y(i); y1=y(i)+h*feval(f2,t(i),x(i),y(i); x2=x(i)+h*feval(f1,t(i+1),x1,y1); y2=y(i)+h*feval(f2,t(i+1),x1,y1); x(i+1)=(x1+x2)/2; y(i+1)=(y1+y2)/2;End在command窗口輸入t,x1,y1=euler2(0,1,4,5/4,100)，得到t，x,y的值，作圖，輸入命令如下： plot(t,x1,'r*-',t,y1,'g*

22、-',t,4*exp(3*t)+2*exp(-t)-2*exp(t),'b+-',t,2*exp(3*t)-exp(-t)+1/4*exp(t),'y+-')；legend('x數(shù)值解','y數(shù)值解','x精確解','y精確解')； xlabel('t'),ylabel('x(y)')，將圖形保存為eulerprove2.fig，圖形如下： 我們?cè)賮肀容^改進(jìn)的歐拉方法和歐拉方法的誤差。輸入： plot(t,abs(x-4*exp(3*t)+2*exp(-t)-

23、2*exp(t),'r*-',t,abs(x1-4*exp(3*t)+2*exp(-t)-2*exp(t),'g*-',t,abs(y-4*exp(3*t)+2*exp(-t)-2*exp(t),'b+-',t,abs(y1-2*exp(3*t)-exp(-t)+1/4*exp(t),'y+-')，legend('x歐拉方法','x改進(jìn)的歐拉方法','y歐拉方法','y改進(jìn)的歐拉方法')；xlabel('t');ylabel('x(y)'

24、;)；xlabel('t');ylabel('x(y)')；title('誤差曲線')，將圖形保存為error2.fig，圖形如下：從這個(gè)例子我們可以很直觀的發(fā)現(xiàn)改進(jìn)的歐拉方法相比歐拉方法具有更高的精度。注意，這里誤差曲線的x軸變量為節(jié)點(diǎn)參數(shù)t，即對(duì)應(yīng)的是固定步長下每一個(gè)節(jié)點(diǎn)處的誤差值。2.3.3（高階微分方程）對(duì)于二階常微分方程，求數(shù)值解解析：先算出其解析解，在MATLAB中輸入：y=dsolve('D2x=-x','x(0)=1','Dx(0)=2')得到解為：Y=2*sin(t)+cos(t)

25、，前面已經(jīng)分別給出過歐拉方法和改進(jìn)的歐拉方法的算例跟誤差比較，這里我們就用精度更高的改進(jìn)歐拉法進(jìn)行數(shù)值求解。改進(jìn)的歐拉方法：先換元，令，則原方程可以轉(zhuǎn)化為，現(xiàn)在，二階常微分方程轉(zhuǎn)化為了一個(gè)一階常微分方程組，同2.3.2的方法，建立M文件eulerprove3.m,內(nèi)容如下：function t,x,y=eulerprove2(t0,tfinal,x0,y0,n)f1=inline('y');f2=inline('-x')if nargin<5,n=50;endh=(tfinal-t0)/n;t(1)=t0,x(1)=x0;y(1)=y0;for i=1:n

26、 t(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');

27、ylabel('y）；legend('數(shù)值解','精確解')，將圖形保存為eulerprove3.fig,圖形如下：上面，我們已經(jīng)通過例子看出，改進(jìn)的歐拉法相比于歐拉法，在每一個(gè)節(jié)點(diǎn)處的誤差值更下，下面，我們來討論節(jié)點(diǎn)的多少（步長大小）對(duì)誤差的影響，創(chuàng)建erro3r.m文件，內(nèi)容如下：function N,Y=error3(n0,nfinal)N(1)=n0,m=fix(nfinal-n0)/4);for i=1:m N(i+1)=N(i)+4; t,x1,y1=eulerprove3(0,1,1,2,N(i); Y(i)=log10(max(x1-2*

28、sin(t)-cos(t); 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,圖形如下:從這個(gè)圖的誤差曲線，隨著插值節(jié)點(diǎn)個(gè)數(shù)增多，誤差也是

29、越來越小的，所以，為了得到盡可能準(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é)合具體要求靈活選取。注：這三個(gè)例子，一個(gè)是一階常微分方程，一個(gè)是二元常微分方程組，第三個(gè)是二階常微分方程，都是初值問題求數(shù)值解。通過這三個(gè)例子可以看出：用歐拉方法和改進(jìn)的歐拉方法求數(shù)值解的一般思路是，通過換元，將高階常微分方程轉(zhuǎn)換為一階常微分方程組，再用2.3.2的方法求解這個(gè)方程組。3龍格-庫塔法3.1 龍格-庫塔法與泰勒展開3.1.1泰勒（Taylor）展開首先，考慮考慮微分方程 （3-1）引入算子符號(hào)， （3-2）于是有 （3-3）設(shè)

30、不帶括號(hào)的及其各階微商都在處取值，于是有泰勒展開式： （3-4）3.1.2 龍格-庫塔（R-K）方法介紹龍格-庫塔方法的基本思路是想辦法計(jì)算在某些點(diǎn)上的函數(shù)值，然后對(duì)這些函數(shù)值做數(shù)值線性組合，構(gòu)造出一個(gè)近似的計(jì)算公式；再把近似的計(jì)算公式和解的泰勒展開式相比較，使得前面的若干項(xiàng)相吻合，從而達(dá)到較高的精度，一般的顯示R-K方法的形式如下： （3-5）當(dāng)式（3-5）的局部截?cái)嗾`差達(dá)到。龍格庫塔方法是應(yīng)用最廣的求解常微分方程初值問題的單步法，下面我們給出幾個(gè)推導(dǎo)公式。3.2 龍格-庫塔法公式與ode函數(shù)3.2.1 二階二段R-K法公式推導(dǎo)由式子（3-5），二階二段的R-K方法可以寫成 （3-6）不帶括

31、號(hào)的及其各階微商都在處取值，由泰勒展開式（3-4）及截?cái)嗾`差的定義可知： （3-7）由于上式是二階二段的R-K法，即必須有，所以 （3-8）式（3-8）有四個(gè)未知元，三個(gè)方程，故有無窮組解。若令，由式子（3-9）可解得于是 這個(gè)公式稱為預(yù)估-校正公式。常用的二階公式還有中點(diǎn)公式（?。?，這里不再列舉。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ì)算公式。（

32、1） 庫塔公式（3階3段） （2）標(biāo)準(zhǔn)四階R-K公式 上述三階公式具有三階精度，四階公式具有四階精度，要驗(yàn)證這一點(diǎn)，我們觀察如下泰勒展開： (3-9)其中落在平面上連接的線段上。只需要按照式（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ì)算，以保持四階精

33、度（下一節(jié)阿達(dá)姆斯預(yù)估校正法中也提到了這一點(diǎn)）。在MATLAB中，專門提供了求解微分方程的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,tsp

34、an,y0,options) ；sol = ode45(odefun,t0tf,y0.) ，其中odefun是函數(shù)句柄，可以是函數(shù)文件名或內(nèi)聯(lián)函數(shù)名，tspan是求解區(qū)間 t0 tf 或者一系列散點(diǎn)t0,t1,.,tf，y0 是初始值向量，T 返回列向量的時(shí)間點(diǎn)，Y 返回對(duì)應(yīng)T的求解列向量 ，options 是求解參數(shù)設(shè)置,可以用odeset在計(jì)算前設(shè)定誤差,輸出參數(shù),事件等 。3.3算例3.3.1（三階三段R-K法MATLAB編程解常微分方程組初值問題）計(jì)算初值問題： 其中解析：先求其解析解，在MATLAB輸入：x,y,z=dsolve('Dx=t*x+Dy-t2','

35、;Dy=t*y+3*t','Dz=t*z-Dy+6*t3','x(0)=1','y(0)=2','z(0)=3','t')，解得：x=5/2*exp(1/2*t2)*t2+t-1/2*exp(1/2*t2)*pi(1/2)*2(1/2)*erf(1/2*2(1/2)*t)+exp(1/2*t2);y =-3+5*exp(1/2*t2)；Z=-5/2*exp(1/2*t2)*t2-6*t2-12+15*exp(1/2*t2).下面我們用三階三段的龍格-庫塔公式來求數(shù)值解，建立名為rk33.m的M文件，內(nèi)容如下

36、：function t,x,y,z=rk33(t0,tfinal,x0,y0,z0,n)f1=inline('t*x+t*y+3*t-t2','t','x','y');f2=inline('t*y+3*t','t','y');f3=inline('t*z-t*y-3*t+6*t3','t','z','y');if nargin<6,n=50;endt(1)=t0,x(1)=1,y(1)=2,z(1)=3;h=(tf

37、inal-t0)/n;for i=1:n t(i+1)=t(i)+h; x1=h*feval(f1,t(i),x(i),y(i); x2=h*feval(f1,t(i)+h/2,x(i)+x1/2,y(i)+h/2); x3=h*feval(f1,t(i)+h,x(i)-x1+2*x2,y(i)+h); x(i+1)=x(i)+(x1+4*x2+x3)/6; y1=h*feval(f2,t(i),y(i); y2=h*feval(f2,t(i)+h/2,y(i)+y1/2); y3=h*feval(f2,t(i)+h,y(i)-y1+2*y2); y(i+1)=y(i)+(y1+4*y2+y3

38、)/6; z1=h*feval(f3,t(i),z(i),y(i); z2=h*feval(f3,t(i)+h/2,z(i)+z1/2,y(i)+h/2); z3=h*feval(f3,t(i)+h,z(i)-z1+2*z2,y(i)+h); z(i+1)=z(i)+(z1+4*z2+z3)/6;End在command窗口輸入t,x,y,z=rk33(0,0.75,1,2,3,20)，返回t,x,y,z的值，我們分別將x,y,z的解析解跟數(shù)值解作在一張圖上，則得到三幅圖形，分別保存為rk33x.fig, rk33y.fig, rk33z.fig.作圖命令如下：（以y曲線為例，另兩幅類似） pl

39、ot(t,-3+5*exp(1/2*t.2),'r*-',t,y,'g+-')；legend('精確解','數(shù)值解')；xlabel('t');ylabel('y') ;title('函數(shù)y曲線圖')。圖形分別如下：下面我們來分析R-K法的誤差跟數(shù)字節(jié)點(diǎn)個(gè)數(shù)（步長）之間的關(guān)系，任取上述三條曲線中的一條進(jìn)行分析，以y曲線為例，建立errorrk.m文件，內(nèi)容如下：function N,Y=errorrk(n0,nfinal)N(1)=n0,m=fix(nfinal-n0)/2);for

40、 i=1:m N(i+1)=N(i)+2; t,x,y,z=rk33(0,0.75,1,2,3,N(i); Y(i)=log10(max(abs(y+3-5*exp(1/2*t.2); t,x1,y1,z1=rk33(0,0.75,1,2,3,N(i+1); Y(i+1)=log10(max(abs(y1+3-5*exp(1/2*t.2);End在窗口輸入：N,Y=errorrk(5,40);plot(N,Y,'+-.');xlabel('N');ylabel('log_1_0max|error|');title('y的誤差曲線'

41、)，將圖片保存為rk33error.fig，圖形如下：從誤差圖形可以看出，R-K法隨著數(shù)值節(jié)點(diǎn)的增多（步長減?。?，精度是不斷提高的。 3.3.2（用自帶ODE45求解高階常微分方程組）將常微分方程組初值問題轉(zhuǎn)化為一階常微分方程組并用ode45求解解析：首先我們嘗試用MATLAB求出該方程組的解析解，在command窗口輸入：x,y=dsolve('D2x+x*(x2+y2)(-3/2)=0','D2y+y*(x2+y2)(-3/2)=0','x(0)=0.5','y(0)=0.25','Dx(0)=0.75',&#

42、39;Dy(0)=1')，運(yùn)行大概五分鐘，提示：Warning: Explicit solution could not be found.x = empty sym ，y =，產(chǎn)生該問題的原因可能在于該非線性方程無法用solve求解，或該方程的解析解不存在。下面，我們用數(shù)值方法求解該方程組，調(diào)用的是函數(shù)庫中的ode45函數(shù)。首先，我們可以做變量替換，令，則該方程組可以化為一個(gè)一階常微分方程組初值問題： 現(xiàn)在，可以直接調(diào)用ode45函數(shù)進(jìn)行數(shù)值計(jì)算。先建立一個(gè)名為fun45.m的M文件來定義方程組，內(nèi)容如下：function y=fun45(t,x)y(1)=-x(3)*(x(3)2+

43、x(4)2)(-3/2);y(2)=-x(4)*(x(3)2+x(4)2)(-3/2);y(3)=x(1);y(4)=x(2);y=y(:)在command窗口輸入： x0=0.75,1,0.5,0.25;t,x=ode45('fun45',0,1,x0),輸出x的值，x是一個(gè)四列的矩陣，四列分別代表a,b,x,y在區(qū)間0,1上的值，可以讀出 ：x(1)=0.5499, y(1)=0.7407，或者作圖：plot(t,x)；legend('a','b','x','y')，xlabel('t')將圖片

44、保存為fun45.fig.圖形如下：通過Data Cursor一樣可以讀出x(1)=0.5499, y(1)=0.7407。注：這兩個(gè)例子一個(gè)是自己編的三階三段的R-K法的一般固定步長程序，一個(gè)是系統(tǒng)自帶的ode45函數(shù)求解，事實(shí)上MATLAB提供了各種R-K法的算法函數(shù)，不僅使用方便，而且步長能夠根據(jù)精度需要自動(dòng)調(diào)整，這一節(jié)，我們重點(diǎn)需要了解ode函數(shù)的使用方法。 4 阿達(dá)姆斯預(yù)估校正法4.1 阿達(dá)姆斯（Adams）公式 在第二節(jié)講述改進(jìn)的歐拉方法時(shí)，我們先通過積分得到了表達(dá)式（2-7），移項(xiàng)后如下: （4-1）在改進(jìn)的歐拉方法中，我們是使用梯形公式來計(jì)算上式中的積分項(xiàng)，事實(shí)上，采用不同的數(shù)

45、值方法計(jì)算它就能導(dǎo)出不同的計(jì)算方法，普遍采用的方法是用的插值多項(xiàng)式來代替，用來代替，這種做法也是線性多步法的一般思路，計(jì)算公式如下： （4-2）取等距節(jié)點(diǎn)，設(shè)已經(jīng)求得這些點(diǎn)上的數(shù)值解，從而有，其中。 令，構(gòu)造r次牛頓向后插值多項(xiàng)式如下： （4-3）整理可以得到遞推式如下（具體推導(dǎo)過程參見參考文獻(xiàn)2）： (4-4)其中。下表給出了的值 012313/2-1/223/12-16/125/1255/24-59/2437/24-9/24當(dāng)時(shí)，(4-4)就是歐拉公式。當(dāng)時(shí)，(4-4)可以寫為： (4-5)式(4-5)就稱為阿達(dá)姆斯-巴士福德公式。用牛頓插值多項(xiàng)式的余項(xiàng)估算可證明該公式是一個(gè)4階方法（證明

46、略，詳見參考文獻(xiàn)2）若取個(gè)節(jié)點(diǎn)為，對(duì)應(yīng)的函數(shù)為，則類似的可推得： （4-6）其中下表給出了的值. 012311/21/25/128/12-1/129/2419/24-5/241/24當(dāng)時(shí)，式（4-6）為 （4-7）稱為梯形公式（即改進(jìn)的歐拉方法）當(dāng)時(shí)，式（4-6）為 （4-8）用牛頓插值多項(xiàng)式的余項(xiàng)可以上述公式是階方法（證明略，詳見參考文獻(xiàn)2），即改進(jìn)的歐拉方法是二階方法（這與第二節(jié)證明相符），（4-8）是一個(gè)4階方法，這個(gè)式子稱為阿達(dá)姆斯-莫爾德公式。 從上述公式的形式可以看出，要求，歐拉方法與改進(jìn)的歐拉方法只需要知道即可，這種公式稱為單步法，第三節(jié)的R-K方法也是一個(gè)單步法，而用阿達(dá)姆斯公

47、式計(jì)算，需要知道，這種方法稱為多步法。4.2 預(yù)估校正方法 所謂預(yù)估校正方法，在數(shù)值算法中，主要就是用于隱式計(jì)算公式的計(jì)算的，可以看出式（4-8）就是一個(gè)隱式計(jì)算公式?；仡櫟诙?jié)的改進(jìn)歐拉方法。為了能夠使改進(jìn)的歐拉方法在迭代時(shí)不需要每一步都解方程，我們用歐拉方法先求一個(gè)預(yù)測(cè)值，再用這個(gè)預(yù)測(cè)值來計(jì)算，即： 這其實(shí)就是一個(gè)預(yù)估校正算法。同樣的，我們可以用公式（4-5）來預(yù)估，再用公式（4-8）來校正，則可以得到一個(gè)一個(gè)預(yù)估校正方法如下： （4-9）其中，（4-9）也是一個(gè)4階算法。同4階R-K算法相比較，阿達(dá)姆斯預(yù)估校正算法精度相同，但是迭代公式更為簡單，計(jì)算量減少較多，但是，阿達(dá)姆斯預(yù)估算法無法

48、自動(dòng)進(jìn)行，其開始的幾個(gè)初值需要通過其他的方法計(jì)算得出（如泰勒展開），一般來說，多步法的前面幾步通常是通過單步法來計(jì)算的，為了保證阿達(dá)姆斯預(yù)估算法的精度，我們可以用龍格-庫塔法來計(jì)算前面幾步。4.3算例4.3.1 （常微分方程初值問題）用Adams預(yù)估校正算法在上求解初值問題： 將結(jié)果與ode45及改進(jìn)的歐拉法求解結(jié)果想比較。解析：先求出解析解，輸入：dsolve('Dy=x*y2','y(0)=1','x')，解得：y =-2/(x2-2)。下面用阿達(dá)姆斯預(yù)估校正法求數(shù)值解。創(chuàng)建M文件adams1.m，內(nèi)容如下：function x,y=adam

49、s1(fun,x0,xfinal,y0,n)h=(xfinal-x0)/n,x(1)=x0,y(1)=y0;for i=1:3 x(i+1)=x(i)+h; y1=h*feval(fun,x(i),y(i); y2=h*feval(fun,x(i)+h/2,y(i)+y1/2); y3=h*feval(fun,x(i)+h/2,y(i)+y2/2); y4=h*feval(fun,x(i)+h,y(i)+y3); y(i+1)=y(i)+1/6*(y1+2*y2+2*y3+y4);endfor i=4:n x(i+1)=x(i)+h; y5=y(i)+h/24*(55*feval(fun,x(

50、i),y(i)-59*feval(fun,x(i-1),y(i-1)+37*feval(fun,x(i-2),y(i-2)-9*feval(fun,x(i-3),y(i-3); y(i+1)=y(i)+h/24*(9*feval(fun,x(i+1),y5)+19*feval(fun,x(i),y(i)-5*feval(fun,x(i-1),y(i-1)+feval(fun,x(i-2),y(i-2);End（前幾步使用的是4階RK標(biāo)準(zhǔn)公示，下同）定義函數(shù)adams1f如下：function y=f1(x,y);y=x*y2.在command窗口輸入：x,y=adams1('adams

51、1f',0,1,1,40)，返回x,y和函數(shù)的值，作圖，命令如下： plot(x,y,'r*-',x,-2./(x.2-2),'g+-');legend('數(shù)值解','精確解');title('y=x*y2')，得到圖形如下，保存為adams1a.fig調(diào)用ode45函數(shù)，輸入x1,y1=ode45(''adams1f',0,1,1),得到用RK法的數(shù)值解（40個(gè)節(jié)點(diǎn)），現(xiàn)在來比較二者之間的誤差，x2,y2=eulerprove1('adams1f',0,1,1,40

52、)輸入： subplot(311);plot(x,abs(y+2./(x.2-2);title('阿達(dá)姆斯法誤差曲線');xlabel('x');ylabel('|error|')；subplot(312);plot(x1,abs(y1+2./(x1.2-2);title('R-K（ode45）誤差曲線');xlabel('x');ylabel('|error|');subplot(313);plot(x2,abs(y2+2./(x2.2-2);title('改進(jìn)的歐拉法誤差曲線')

53、;xlabel('x');ylabel('|error|')得到如下誤差曲線圖（都是40個(gè)數(shù)值節(jié)點(diǎn)），保存為adams1b.fig.從誤差曲線可以直觀的看出，阿達(dá)姆斯預(yù)估校正法與4階RK法精度相當(dāng)，相比改進(jìn)的歐拉法具有明顯的精度提升，事實(shí)上，阿達(dá)姆斯法與ode45都是四階方法，改進(jìn)的歐拉法是二階方法，在相同數(shù)目的取值節(jié)點(diǎn)情況下，前兩者精度更高，這個(gè)圖形也證明了這一觀點(diǎn)。4.3.2（常微分方程組初值問題）初值問題： ，先求出精確解，再用數(shù)值方法求解。解析：在MATLAB中輸入：y,z=dsolve('Dy=y+4*z-exp(x)','Dz

54、=y+z+2*exp(x)','y(0)=2','z(0)=3','x')得：y =-3/4*exp(-x)+19/4*exp(3*x)-2*exp(x)；z =3/8*exp(-x)+19/8*exp(3*x)+1/4*exp(x)。下面我們用阿達(dá)姆斯預(yù)估校正法來求其數(shù)值解，建立文件adams2.m，內(nèi)容如下：function x,y,z=adams2(x0,xfinal,y0,z0,n)f1=inline('y+4*z-exp(x)','x','z','y');f2=in

55、line('y+z+2*exp(x)','x','y','z');x(1)=x0,y(1)=y0,z(1)=z0;h=(xfinal-x0)/n;for i=1:3 x(i+1)=x(i)+h; y1=h*feval(f1,x(i),z(i),y(i); y2=h*feval(f1,x(i)+h/2,z(i)+h/2,y(i)+y1/2); y3=h*feval(f1,x(i)+h/2,z(i)+h/2,y(i)+y2/2); y4=h*feval(f1,x(i)+h,z(i)+h,y(i)+y3); y(i+1)=y(i)+1/

56、6*(y1+2*y2+2*y3+y4); z1=h*feval(f2,x(i),y(i),z(i); z2=h*feval(f2,x(i)+h/2,y(i)+h/2,z(i)+z1/2); z3=h*feval(f2,x(i)+h/2,y(i)+h/2,z(i)+z2/2); z4=h*feval(f2,x(i)+h,y(i)+h,z(i)+z3); z(i+1)=z(i)+1/6*(z1+2*z2+2*z3+z4);endfor i=4:n x(i+1)=x(i)+h; y5=y(i)+h/24*(55*feval(f1,x(i),z(i),y(i)-59*feval(f1,x(i-1),z(i-1),y(i-1)+37*feval