歐拉法解常微分方程

1、實用標準文案大全數(shù)學與計算科學學院驗報告實驗項目名稱 Eular方法求解一階常微分方程數(shù)值解所屬課程名稱偏微分方程數(shù)值解實驗類型驗證性實驗日期 2015-3-26班級學號姓名成績一、實驗概述: 【實驗目的】熟練掌握應用顯性Eular法和隱式Eular法求解一般一階常微分方程的近似數(shù) 值解?！緦嶒炘怼侩m然求解常微分方程有各種各樣的解析方法，但解析方法只能用來求解一些特 殊類型的方程。求解從實際問題當中歸結(jié)出來的微分方程主要靠數(shù)值解法。歐拉方法是一類重要的數(shù)值解法。這類方法回避解y(x)的函數(shù)表達式，而是尋求它在一系 列離散節(jié)點上的近似值，相鄰的兩個節(jié)點的間距稱作步長。假定步長為定數(shù)。歐拉方法是

2、一類離散化方法，這類方法將尋求解y(x)的分析問題轉(zhuǎn)化為計算離 散值值的代數(shù)問題，從而使問題獲得了實質(zhì)性的簡化。然而隨之帶來的困難是，由于數(shù)據(jù)量往往很大，差分方法所歸結(jié)出的可能是個大規(guī)模的代數(shù)方程組。 【實驗環(huán)境】1.硬件環(huán)境穌處理器安裝內(nèi)存(MM);系統(tǒng)超：葦設低插：系統(tǒng)汾丞詢用Intel(R) Core(TM) i3-311OM CPU 240GH? 2.40 GHz4,00 GB (3.85 G3 可用)64位逐f修珠沒有可用于此顯示器的莖或在控墟人2.2.軟件環(huán)境Windows 版本W(wǎng)indows 7 揖®!版版院侑 © 2009 Microsoft Corpora

3、tion,保翟所Wt又機Service Pack 1MATLAB7.0二、實驗內(nèi)容：【實驗過程】(實驗步驟)(一)實驗任務描述某種化學反應過程的方程，利用顯性和隱形Eualar方法求解下列一階線性微分方程組的近似數(shù)值解：d%=-0.04yi +10丫2jdy%=0.04yi+l04yiy2-3Ml07y2dy =3-y2 dt.(0)=1, y2(0)=Q y3(0)=0(二)求解過程Eular方法：一階線性微分方程初值問題7= f (x,y),a <x <b <、y(a) = y0a=xo<x1<.<xn=b (1)Xn = Xo nh, h為步長方程離散

4、化：差分和差商yi -yoyi - yOy (xo)：Xi -xohf(xo, yo) =* J0 h yi =y0 +hf (x°,yo)(2)Vn 1 =yn hf(x0,yo)通過初始值yo,依據(jù)遞推公式(2)逐步算出yi,y2,., yn就為顯性的Eular方 法。隱形Eular方法:（3）、=y0 +hf（x1，yi）=Jn4t=yn +hf Un由，Yn 十）公式(3)即為隱式Eular公式。(三)程序算法1.利用顯式Eular法方求解利用MATLA進行求解，編寫腳本文件如下:文件名：hql.m% 顯性Eular方法f0=1; g0 =0;z0=0delta=0.01;t

5、ime=1;t=0:delta:time;f=zeros(size(t);g=zeros(size(t);z=zeros(size(t);f1=zeros(size(t);g1=zeros(size(t);z1=zeros(size(t);f(1)=f0;g(1)=g0；z(1)=z0;for i=2:length(t)f1(i-1) = -0.04*f(i-1) + 10000*f(i-1)*g(i-1);f(i)=f(i-1)+f1(i-1)*delta;g1(i-1) = 0.04*f(i-1) - 10000*f(i-1)*g(i-1)-3*10A7*g(i-1)A2;g(i)=g(i

6、-1)+g1(i-1)*delta；z1(i-1)=3*10A7*g(i-1)A2;z(i)=z(i-1)+z(i-1)*delta;Fun=f+g+zendfigureplot(t,f, 'o')；xlabel('t');ylabel('y1');title( 't-y1 變化圖')figureplot(t,g, 'o');xlabel('t');ylabel('y2');title( 't-y2 變化圖')figureplot(t,z, 'o');

7、xlabel('t');ylabel('y3');title( 't-y3 變化圖')figureplot(t,Fun);xlabel('t');ylabel('y1+y2+y3');title( 't-y1+y2+y3 變化圖)【實驗結(jié)論】A步長h=0.001時進行數(shù)據(jù)測試。結(jié)果如下:迭代第一次時,l.OOQOe-* 1x101 de0Columns 12 through 22結(jié)果與方程描述內(nèi)容相符。迭代第二次時,Columns I through 111.0000 L 0000E00000.9959結(jié)果

8、與方程描述內(nèi)容基本相符。迭代三次時,Fun -Columns 1 thraugh 11EOOOO k 0000 L 00000.99590.980800八-alColumns 12 through 22結(jié)果與方程描述內(nèi)容基本相符。5-Q5(H)Q也MM也MM-0, D0007, WOO-0.0M03MlNaXFlaNNaNMMTJaNNaffNaJiNaffNaifMM迭代1000次時,Fun =和 3 *Caliums I through i I6 MM D.g 鐘 MMgColumns 12 thraufh 22二工 lfi-55 -InfColiuicas 23 through 33Hm

9、兵 Colujufts 3l through 44模擬結(jié)果已經(jīng)嚴重脫離事實，故當選擇delta為0.001時，該迭代方法不收斂。時間與個變量直接的變化關(guān)系如圖所示:從上述圖形可以明顯看出，在迭代的不斷進行時，各變量與時間的變化越來越 大，且嚴重脫離了方程所描述的現(xiàn)實意義。B.當選擇h=0.00000001時，模擬結(jié)果如下:迭代第一次,Fun =Columns 1 through IE1100Columns 19 through 350000與A中結(jié)果相同。Fun =Colwhns1Columns迭代第二次，1 through 1811019 through 36跌二次迭代結(jié)果明顯優(yōu)于一中。跌三

10、次迭代結(jié)果,Fun -Columns值l.OOOOe-A- 1x101 di1"01 de1x1 1 10JColumns1 through 1811 I19 through 36并未產(chǎn)生誤差。地1000次迭代結(jié)果,Colums 1 thxouEh 10I.OWl.omLQOQQLOWL MML WQL QMQLMML MQQColuw 11 thrcniEh 20lmmLornl omLornLmolmqqlmq。i , 0000LQggCdliwH 21 thrwili JO構(gòu)AuomLom Lmolmqq lmw l qoooColiwis. 31 thmugh 40lowI,

11、l GmLornk«w。lwoolmel, coool oooo1shihdt W1 throuih 5Ql.omIRDWJLQQMLDmIL WMQLMQQL(NMI, MMLQQgL QEQCglvims 5：l tkiuEh <501通1.八fiJOJ dJ |311有巾J泰1 mi Jpl 3A dt a 上de V >ft結(jié)果明顯是收斂的。時間與個變量直接的變化關(guān)系如圖所示:t-y1 +y2+y凌化圖00.10.20.30.40.50,60 70.80.91t乂d從圖中能夠清晰看出，當 h=0.00000001時，模擬結(jié)果與方程所表示的顯示意義相吻合。說明了顯性

12、Eualr方法的收斂性是與步長的選擇是相關(guān)。這就對我們們選擇步長造成了困難， 由于選擇的步長不合適有可能得出錯誤的結(jié)論。【實驗小結(jié)】（收獲體會）1、軟件使用在寫MATLABS言的時候要深刻理解題的意圖，整理好思緒再做題目，在我運算 的過程中，h取值取得越小、越細微，曲線逼近的越好。2、歐拉法的缺點簡單地取切線的端點作為下一的起點進行計算，當步數(shù)增多時，誤差會因積累 而越來越大。因此歐拉格式一般不用于實際計算。3、實驗感想在這次上機實驗中，我掌握了解決常微分方程的基本方法，同時學會使用計算機 軟件對兩種不同方法得到的結(jié)果進行判斷，對我們以后對數(shù)據(jù)進行分析很有祁助。三、指導教師評語及成績：評語評語

13、等級優(yōu)良中及 格/、及格1.實驗報告按時完成，字跡清楚，文字敘述流暢，邏輯性強2.實驗方案設計合理3.實驗過程（實驗步驟詳細，記錄完整，數(shù)據(jù)合理，分析透徹）4實驗結(jié)論正確.成績：指導教師簽名： 批閱日期：附錄：源程序程序1:%&Ti Eular 方法f0=1; g0 =0;z0=0delta=0.00000001;time=0.00001;t=0:delta:time;f=zeros(size(t);g=zeros(size(t);z=zeros(size(t);f1=zeros(size(t);g1=zeros(size(t);z1=zeros(size(t);f(1)=f0;g(1

14、)=g0；z(1)=z0;for i=2:length(t)f1(i-1) = -0.04*f(i-1) + 10000*f(i-1)*g(i-1);f(i)=f(i-1)+f1(i-1)*delta;g1(i-1) = 0.04*f(i-1) - 10000*f(i-1)*g(i-1)-3*10A7*g(i-1)A2;g(i)=g(i-1)+g1(i-1)*delta；z1(i-1)=3*10A7*g(i-1)A2;z(i)=z(i-1)+z(i-1)*delta;Fun=f+g+zendfigureplot(t,f, 'o')；xlabel( 't');ylabel( 'y1');title( 't-y1 變化圖)figureplot(t,g, 'o')；xlabel( 't');ylabel( '