1、常微分方程初值問(wèn)題的歐拉方法及其改進(jìn)的歐拉方法的 Matlab 實(shí)現(xiàn)紀(jì)秀浩遼寧工程技術(shù)大學(xué)理學(xué)院,遼寧阜新(123000E-mail :摘 要 :歐拉 (Euler方法及改進(jìn)的歐拉方法是解決常微分方程初值問(wèn)題常用的數(shù)值解法, 但 Matlab 的工具箱中沒(méi)有 Euler 方法的功能函數(shù)。本文在簡(jiǎn)要介紹 Euler 方法及其改進(jìn)的 Euler 方法的基礎(chǔ)上, 通過(guò)編寫(xiě) Matlab 程序?qū)崿F(xiàn)兩種數(shù)值解法, 并通過(guò)作圖形式對(duì)比其精度, 加深對(duì)兩種方法的認(rèn)識(shí)。關(guān)鍵詞:歐拉方法;改進(jìn)的歐拉方法; matlab 實(shí)現(xiàn)1.引言常微分方程是解決工程實(shí)例的常用的工具 1,建立微分方程只是解決問(wèn)題的第一步,通

2、常需要求出方程的解來(lái)說(shuō)明實(shí)際現(xiàn)象, 并加以檢驗(yàn)。 如果能得到解析形式的解固然是便于分 析和應(yīng)用的, 但是我們知道, 雖然求解常微分方程有各種各樣的解析方法, 但解析方法只能 用來(lái)求解一些特殊類(lèi)型的方程,實(shí)際問(wèn)題終歸結(jié)出來(lái)的微分方程主要靠數(shù)值解法 2。數(shù)值解 法就是一個(gè)十分重要的手段,而歐拉方法又是數(shù)值解法最基礎(chǔ)最常用的方法。2. 歐拉方法、改進(jìn)的歐拉方法及 Matlab 實(shí)現(xiàn)下面主要討論一階常微分方程的初值問(wèn)題,其一般形式為:' 00(, ( y f x y y x y = (1 我們知道,只要函數(shù) (, f x y 適當(dāng)光滑 譬如關(guān)于 y 滿(mǎn)足利普希茨 (Lipschitz條件(,

3、(, f x y f x y L y y 理論上就可以保證初值問(wèn)題 (1的解 ( y y x =存在并且唯一 3。所 謂 數(shù) 值 解 法 , 就 是 求 問(wèn) 題 (1 在 某 些 離 散 點(diǎn) 01N a x x x b =<<=" 的 近 似 值012, , , N y y y y " 的方法。 012, , , N y y y y " 就稱(chēng)為問(wèn)題 (1的數(shù)值解。 1n n n h x x +=成為 n x 到 n x 的步長(zhǎng),我們?yōu)榱朔奖闳槌Ａ?h 。 2.1. 歐拉方法2.1.1 歐拉方法將微分方程離散化,用向前差商 1( ( n n y x y

4、 x h+代替微分 '( n y x ,代入 (1中的微分方程,可得1( (, (1, 2,3, n n n n y x y x f x y x n h+="化簡(jiǎn)可得1( ( (, ( (1, 2,3, n n n n y x y x f x y x hn +=+="中國(guó) 科技論文在線(xiàn)如果用 n y 近似 ( n y x 代入上式便可得到 1( n y x +的近似值 1n y +,計(jì)算式為:1(, (1, 2,3, n n n n y y f x y h n +=+=" (2 這樣問(wèn)題 (1的近似解可通過(guò)求解下面的差分初值問(wèn)題:10(, (1,2,3,

5、(n n n n y y f x y h n y y a +=+=" (3得到,按 (3式由初值 0y 可逐次求出 12, , y y " 。Eule 方法就是用差分方程初值問(wèn)題 (3的解來(lái)近似微分方程初值問(wèn)題 (1的解。即由公式 (3算出 ( n y x 的近似值。這組公式求問(wèn)題 (1的數(shù)值解就是著名的歐拉 (Euler公式。 2.1.2歐拉方法的誤差估計(jì)對(duì)于 Euler 公式 (3我們看到,當(dāng) 1, 2, n =" 時(shí)公式右端的 n y 都是近似的,所以用它計(jì)算 的 1n y +會(huì)有累積誤差,分析累積誤差比較復(fù)雜,這里先討論比較簡(jiǎn)單的所謂局部截?cái)嗾`差。假定用

6、(3式時(shí)右端的 n y 沒(méi)有誤差,即 ( nn y y x =那么由此算出1( (, ( n n n n y y x f x y x h +=+ (4局部截?cái)嗾`差指的是, 按 (4式計(jì)算由 n x 到 1n x +這一步的計(jì)算值 1n y +與精確值 1( n y x +之差11( n n y x y +。為了估計(jì)它,由 Taylor 展開(kāi)得到的精確值 1( n y x +是2'''31( ( ( ( ( 2n n n n h y x y x hy x y x O h +=+ (5(4、 (5兩式相減(注意到 '(, y f x y =得2''3

7、211( ( ( ( 2n n n h y x y y x O h O h +=+ (6即局部截?cái)嗾`差是 2h 階的,而數(shù)值算法的精度定義為:若一種算法的局部截?cái)嗾`差為1( p O h +,則稱(chēng)該算法具有 p 階精度。顯然 p 越大,方法的精度越高。式 (6說(shuō)明, Euler方法是一階方法,因此它的精度不高。2.2 改進(jìn)的歐拉方法2.2.1 改進(jìn)的歐拉方法用數(shù)值積分方法離散化問(wèn)題 (1,兩端積分可得11( ( (, ( (0,1, 2, n nx n n x y x y x f x y x dx n +="對(duì)右端積分使用梯形公式可得 ,111(, ( (, ( (, (2n nx n

8、 n n n x hf x y x dx f x y x f x y x + 再用 1, n n y y +代替 1(, ( n n y x y x +,則得計(jì)算公式111(, (, 2n n n n n n hy y f x y f x y +=+ (7很明顯可以注意到 (7式為隱式形式。改進(jìn)的歐拉方法是先用歐拉公式求 1( n y x +的一個(gè)近似值 1n y +,稱(chēng)為預(yù)測(cè)值,然后用梯 形公式進(jìn)行矯正并求得近似值 1n y +。即1111(, (, (, 2n n n n n n n n n n y y f x y h hy y f x y f x y +=+=+(8 2.2.2 改進(jìn)的歐

9、拉方法的誤差估計(jì)由 (7式可知,11232434( ( '( '(2'( "( "'(23!( '( ''( '''( 22'''( ( 12n n n n n n n n n n n n hy x y x y x y x h h hy x y x y x h hy x y x hy x y x O h h y x O h +=+=+所以改進(jìn)的歐拉方法是二階的,精度較歐拉方法要高,實(shí)用性更加廣泛。3.實(shí)驗(yàn)及結(jié)果分析3.1 歐拉方法的 matlab 實(shí)現(xiàn)及實(shí)例歐拉方法的 M

10、atlab 實(shí)現(xiàn)程序如下: function x,y=euler(fun,x0,xfinal,y0,n; if nargin<5,n=50; endh=(xfinal-x0/n; x(1=x0;y(1=y0; for i=1:nx(i+1=x(i+h;y(i+1=y(i+h*feval(fun,x(i,y(i; end例 4 求解初值問(wèn)題'2(01(01x y y x y y =<<=(9解 編寫(xiě)函數(shù)文件 doty.m 如下: function f=doty(x,y; f=y-2*x/y;在 Matlab 命令窗口輸入: x,y=euler('doty'

11、;,0,1,1,10便可得到結(jié)果:x=0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 y=1 1.1000 1.1918 1.2774 1.3582 1.435 11.5090 1.5803 1.6498 1.7178 1.7848 初值問(wèn)題 (6 的精確解為:y =, 此時(shí)對(duì)應(yīng) n x 的精確解 ( n y x 可編程:y1=sqrt(1+2.*x得到 y1=1 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492 1.6125 1.6733 1.732

12、1下面通過(guò)作圖比較歐拉數(shù)值解與精確解之間的誤差, 編程:plot(x,abs(y-y1,'*', 可得到 圖 2-1: 圖 2-1 歐拉數(shù)值解與精確解的誤差圖從圖中可以看出歐拉方法精度很差。3.1 改進(jìn)歐拉方法的 matlab 實(shí)現(xiàn)為了編程方便,常把 (5式改寫(xiě)為:11(, (, (, 21( 2p n n n q nn n n p n p q y y f x y h h y y f x y f x y y y y +=+=+=+則 Matlab 實(shí)現(xiàn)程序?yàn)?function x,y=eulerpro(fun,x0,xfinal,y0,n; if nargin<5,n=5

13、0; endh=(xfinal-x0/n; x(1=x0;y(1=y0; for i=1:nx(i+1=x(i+h;y1=y(i+h*feval(fun,x(i,y(i; y2=y(i+h*feval(fun,x(i+1,y1;中國(guó) 科技論文在線(xiàn)y(i+1=(y1+y2/2; end仍然選用上面的實(shí)例 (6,函數(shù) doty.m 文件相同,在 Matlab 命令窗口中輸入: x,y=eulerpro('doty',0,1,1,10便可得到x=0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.000

14、0 y=1 1.0959 1.1841 1.2662 1.3434 1.4164 1.4860 1.5525 1.6165 1.6782 1.7379 y1=1 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492 1.6125 1.6733 1.7321 y 為改進(jìn)的歐拉數(shù)值解, y1為精確解。 從上面的數(shù)據(jù)可看出, 改進(jìn)的歐拉方法有較好的精確 度。作圖 plot(x,abs(y-y1,'*'可得到誤差如圖 2-2: x 10-3圖 2-2 改進(jìn)的歐拉方法數(shù)值解與精確解的誤差圖4.總結(jié)本文通過(guò)編程實(shí)現(xiàn)了常微分方程初值問(wèn)題數(shù)值解法中的

15、歐拉方法及其改進(jìn)的算法, 并比 較了其數(shù)值解與精確解之間的誤差。 可以看出歐拉方法得到的數(shù)值解精確度較差, 而改進(jìn)的 歐拉方法得到的結(jié)果則相對(duì)較好。中國(guó) 科技論文在線(xiàn)中國(guó)科技論文在線(xiàn) 參考文獻(xiàn) 1 2 3 4 王高熊，周之銘，朱思銘，等數(shù)值計(jì)算原理M 北京：清華大學(xué)出版社，2000 李慶揚(yáng)，王能超，易大義 數(shù)值分析M 北京：清華大學(xué)出版社，施普林格出版社，2001 李慶揚(yáng)，關(guān)治，白峰杉數(shù)值計(jì)算原理M 北京：清華大學(xué)出版社，2000 李慶揚(yáng)，王能超，易大義數(shù)值分析M 武漢：華中科技大學(xué)出版社，2001 Matlab Programing For Eulers Methord and Its improved Methord Ji Xiuhao College of Science, Liaoning Technical University, Fuxin, Liaoning (123000 Abstract This paper makes matlab programs to solve the initial value problem of ordinary differential equations using Euler's methord and Euler's improved methord.It shows the accurations