Euler方法與改進(jìn)的Euler方法的應(yīng)用

1、CENTRAL SOUTH UNIVERSITY數(shù)值分析實驗報告Euler方法與改進(jìn)的Euler方法的應(yīng)用一、問題背景在工程和科學(xué)技術(shù)的實際問題中，常需求解微分方程，但常微分方程中往往只有少數(shù)較簡單和典型的常微分方程（例如線性常系數(shù)常微分方程等）可求出其解析解,對于變系數(shù)常微分方程的解析求解就比較困難，而一般的非線性常微分方程的求解困難就更不用說了。大多數(shù)情況下，常微分方程只能用近似方法求解。這種近似解法可分為兩大類：一類是近似解析法，如級數(shù)解法、逐次逼近法等;另一類是數(shù)值解法，它給出方程在一些離散點上的近似值。二、數(shù)學(xué)模型在具體求解微分方程時，需具備某種定解條件，微分方程和定解條件合在一起組

2、成定解問題。定解條件有兩種：一種是給出積分曲線在初始點的狀態(tài)，稱為初始條件，相應(yīng)的定解問題稱為初值問題。另一類是給出積分曲線首尾兩端的狀態(tài)，稱為邊界條件，相應(yīng)的定解問題稱為邊值問題。在本文中主要討論的是給定初值條件的簡單Euler方法和改進(jìn)的Euler方法來求解常微分方程。三、算法及流程Euler方法是最簡單的一種顯式單步法。對于方程考慮用差商代替導(dǎo)數(shù)進(jìn)行計算，取離散化點列，則得到方程的近似式即得到簡單Euler方法。具體計算時由出發(fā)，根據(jù)初值，逐步遞推二得到系列離散數(shù)值。簡單Euler方法計算量小，然而精度卻不高，因而我們可以構(gòu)造梯形公式其中。這是一個二階方法，比Euler方法精度高。但是上

3、述公式右邊有，因而是隱式差分方程，可以用迭代方法計算。初值可以由Euler公式提供，一般而言迭代一兩次即可。在計算中的迭代公式為容易看出這實際上是一種預(yù)估-校正方法。改進(jìn)的Euler方法又稱為Henu方法。MATLAB實現(xiàn)過程：（1）簡單Euler方法函數(shù)functiont,x=Euler(fun,t0,tt,x0,N) % MyEuler用前向差分的歐拉方法解微分方程 % fun表示f（t,x） % t0,tt表示自變量的初值和終值 % x0表示函數(shù)在x0處的值，其可以為向量形式 % N表示自變量在t0,tt 上取得點數(shù) h=(tt-t0)/N; % 步長h t=t0+0:N'*h;

4、 % 時間t x(1,:) = x0' % 賦初值 for k = 1:N f=feval(fun,t(k),x(k,:); f=f' x(k+1,:) =x(k,:)+h*f; end 將文件以文件名Euler.m保存。（2）改進(jìn)的Euler方法函數(shù)function t,x=GjEuler(fun,t0,tt,x0,N) %改進(jìn)的Euler方法解微分方程 h=(tt-t0)/N;%計算所取的兩離散點之間的距離 t=t0+0:N'*h;%表示出離散的自變量x x(1,:)=x0' for i=1:N f1=h*feval(fun,t(i),x(i,:); f1=

5、f1' f2=h*feval(fun,t(i+1),x(i,:)+f1); f2=f2' x(i+1,:)=x(i,:)+1/2*(f1+f2);end將文件以文件名GjEuler.m保存。四、計算結(jié)果與分析計算常微分方程的初值問題編寫函數(shù)文件，打開Editor編輯器，輸入以下語句并以文件名GjEuler.m保存。function f=GjEuler(t,x) f=1/t*(x2+x);再編寫主函數(shù)文件，打開Editor編輯器，輸入下列語句并以文件名GjEuler_main.m保存文件。function GjEuler_main() % 比較改進(jìn)Euler法、簡單Euler法及

6、文芬方程符號解t,x=Euler('GjEuler_fun',1,3,-2,15); % Euler方法¨ tgj,xgj=GjEuler('GjEuler_fun',1,3,-2,15); % 改進(jìn)的Euler方法的解sh=dsolve('Dx=1/t*(x2+x)','x(1)=-2','t'); % 符號計算for k=1:16 % 時間離散化 st(k)=t(k); % 取數(shù)值方法的離散的時間 sx(k)=subs(sh,st(k); end plot(t,x,'*',tgj,xg

7、j,'+',st,sx) t,x,xgj,sx' % t 時間% x 簡單Euler方法計算結(jié)果 % xgj 改進(jìn)Euler方法的計算結(jié)果% sx 符號計算結(jié)果 t,sx'-x,sx'-xgj % 兩種方法的計算誤差 %plot(t,sx'-x,'+',t,sx'-xgj,'*') 運行程序，在MATLAB命令窗口輸入：>>GjEuler_main運行結(jié)果如下所示：ans = 1, -2, -2, -2 17/15, -26/15, -6854/3825, -34/19 19/15, -38/2

8、37/5，-6694439885567781/4503599627370496,-43-14/923/15,-3-33468975972395-46/31669915380311/4503599627370496,-10/79/5,-5972253212408169/4503599627370496,-6244734289756621/4503599627370496,-18/1329/15,-14531/15,-285-297-62/4711/5,7/3,-5527499547085519/4503599627370496, -5738095816108559/450359962737049

9、6, -14/1137/15,-5455688972551323/4503599627370496,-5654358103897161/4503599627370496,-74/5913/5,41/15,-5338702448959659/4503599627370496,-275-82/6743/15,-529-272-86/713, -13521-6/5其中時間和改進(jìn)的Euler方法兩列的結(jié)果序列就是我們要的結(jié)果。為了分析問題，下圖給出了分析解和Euler方法的結(jié)果進(jìn)行比較。其中“*”線代表的是Euler法，“+”線代表的是改進(jìn)Euler法，“”線代表的是分析解?？梢钥闯?，改進(jìn)的Euler

10、方法比原方法要精確許多，在本例中的區(qū)間內(nèi)，幾乎改進(jìn)后的方法和分析方法重合，而Euler方法則有比較大的誤差。運行結(jié)果中還有下列數(shù)據(jù)：ans = 1, 0, 017/15, -16/285, 176/7267519/15, -6016/87975, 7/5, -2800435813076915/40532396646334464, 6797504630227/253327479039590423/15,-578292948083527/8725724278030336,5/3,-984409/5, -3425501531362731/58546795155816448, 29/15,-2651331/15, 11/5, -36899857630325880151/382805968326492167/3, -2247899765246235/49539595901075456, 68659194007205/4953959590107545637/15,-11380713/5, -3833349331259339/94575592174780416, 56376328559315/47287796087390208 3, -39224285523514680241/5629499534213120上面結(jié)果給出的是簡單E