1、浙江大學城市學院實驗報告課程名稱數值計算方法實驗項目名稱 常微分方程初值問題的數值解法實驗成績指導老師(簽名)日期 2015/12/16一 實驗目的和要求1. 用Matlab軟件掌握求微分方程數值解的歐拉方法和龍格-庫塔方法；2. 通過實例學習用微分方程模型解決簡化的實際問題。二實驗內容和原理編程題2-1要求寫岀Matlab源程序(m文件)，并有適當的注釋語句；分析應用題2-2 , 2-3 , 2-4, 2-5要求將問題的分析過程、Matlab源程序和運行結果和結果的解釋、算法的分析寫在實驗報告上。2-1編程 編寫用向前歐拉公式和改進歐拉公式求微分方程數值解的Matlab程序，問題如下：在區間

2、b,b 1內(N 1)個等距點處，逼近下列初值問題的解，并對程序的每一句添上注釋語句。Euler 法 y=euler(a,b,n,y0,f,f1,b1) 改進 Euler 法 y=eulerpro(a,b,n,y0,f,f1,b1)2-2分析應用題假設等分區間數 n =100，用歐拉法和改進歐拉法在區間L 0,10內求解初值問題y二 y(t)-207(0) -10并作岀解的曲線圖形，同時將方程的解析解也畫在同一張圖上，并作比較，分析這兩種方法的精度。2-3分析應用題用以下三種不同的方法求下述微分方程的數值解，取h=10畫岀解的圖形，與精確值比較并進行分析。1) 歐拉法；2) 改進歐拉法；3)

3、龍格庫塔方法；2-4分析應用題考慮一個涉及到社會上與眾不同的人的繁衍問題模型。假設在時刻t(單位為年)，社會上有人口 X(t)人，又假設所有與眾不同的人與別的與眾不同的人結婚后所生后代也是與眾不同的人。而固 定比例為r的所有其他的后代也是與眾不同的人。如果對所有人來說岀生率假定為常數b，又如果普通的人和與眾不同的人的婚配是任意的，則此問題可以用微分方程表示為： 其中變量p(t) =Xj(t)；x(t)表示在時刻t社會上與眾不同的人的比例，x(t)表示在時刻t人口中與眾不同的人的數量。1) 假定p(0) =0.01,0.02和r =0.1，當步長為h =1年時，求從t = 0到t =50解p(t

4、)的近似值， 并作岀近似解的曲線圖形。2) 精確求岀微分方程的解p(t)，并將你當t =50時在分題(b)中得到的結果與此時的精確值進行比 較。MATLAB相關函數】求微分方程的解析解及其數值的代入dsolve( 6gnT, egn2' - X'subs(expr,x,y, ,x1,y1, -)其中egni '表示第i個方程，X '表示微分方程中的自變量，默認時自變量為t osubs命令中的expr、x、y為符合型表達式，x、y分別用數值x1、x2代入。>>symsxyz>>subs('x+y+z',x,y,z,1,2,3

5、)ans=6>>symsx>>subs('xA2',x,2)ans=4注 s=dsolve( Dy =1 亠 y 2' y(0) =1', x'ans=>>symsx>>subs(s,x,2)ans=-0.3721右端函數f (x, y)的自動生成f=inline( expr','var1', var2',)其中expr'表示函數的表達式，varl ', var2 '表示函數表達式中的變量，運行該函數，生成一個新的函數表達式為f(var1,var2,)

6、。>>f=i nli ne('x+3*y','x','y')仁In li nefun cti on:f(x,y)=x+3*y>>f(2,3)ans=114，5階龍格-庫塔方法求解微分方程數值解t,x=ode45(f,ts,x0,opti ons)其中f是由待解方程寫成的m文件名；x0為函數的初值；t,x分別為輸岀的自變量和函數值(列向量)，t的步長是程序根據誤差限自動選定的。若ts=t0,t1,t2,tf，則輸岀在自變量指定值，等步長時用ts=t0:k:tf ，輸岀在等分點；options用于設定誤差限(可以缺省，缺省時設

7、定為相對誤差 10"，絕對誤差 10-6)，程序為：options=odeset( ' reltol ' ,rt, ' abstol ' ,at)，這里rt,at 分別為設定的相對誤差和絕對誤差。常用選項見下表。選項名功能可選值省缺值AbsTol設定絕對誤差正數RelTol設定相對誤差正數Ini tialStep設定初始步長正數:自動MaxStep設定步長上界正數MaxOrder設定ode15s的最高階數1,2,3,4,55Stats顯示計算成本統計on ,offoffBDF設定ode15s是否用反向差分on ,offoff例：解微分方程在命令窗口執行

8、注 odefun = inline ( y -2* t. y ', t', y'z t, y 丨-ode45(odefun,0,4,1)；ans=01.00000.05021.04900.10051.09590.15071.14083.85072.95033.90052.96723.95022.98394.00003.0006 plot( t, y, 0-'%解函數圖形表示.ode45(odefu n,0,4,1) %不用輸岀變量，則直接輸岀圖形注 t, y J-ode45(odefun,0 :4,1) ；It, y ans=01.00001.00001.732

9、12.00002.23613.00002.64584.00003.0006三操作方法與實驗步驟(包括實驗數據記錄和處理)2-1編程編寫用向前歐拉公式和改進歐拉公式求微分方程數值解的Matlab程序，問題如下：在區間la,b 內(N1)個等距點處，逼近下列初值問題的解，并對程序的每一句添上注釋語句Euler 法 y=euler(a,b,n,y0,f,f1,b1)改進 Euler 法 y=eulerpro(a,b,n,y0,f,f1,b1)Euler 法y=euler(a,b, n,y 0,f,f1,b1)y=zeros(1, n+1);y(1)=y0;h=(b-a)/n;x=a:h:b;fori

10、=1: n;y(i+1)=y(i)+h*f(x(i),y(i);endplot(x,y)hold on%求微分方程的精確解x1=li nspace(a,b,100);'精確解為s=dsolve(f1,b1,'x')symsxy1=zeros(1,100);fori=1:100y1(i)=subs(s,x,x1(i);endplot(x1,y1,'r')title('紅色代表精確解')改進Euler法y=eulerpro(a,b ,n,y 0,f,f1,b1)%求微分方程的數值解y=zeros(1, n+1);y(i)=yo；h=(b-a)

11、/n;x=a:h:b;fori=1: n;T1=f(x(i),y(i);T2=f(x(i+1),y(i)+h*T1);y(i+1)=y(i)+(h/2)* (T1+T2);endPlot(x,y)hold on%求微分方程的精確解x1=li nspace(a,b,100);'精確解為s=dsolve(f1,b1,'x')symsxy1=zeros(1,100);fori=1:100y1(i)=subs(s,x,x1(i);endplot(x1,y1,'r')title('紅色代表精確解')2-2分析應用題假設等分區間數n =100，用歐拉

12、法和改進歐拉法在區間t 0,10內求解初值問題y(0) -10并作岀解的曲線圖形，同時將方程的解析解也畫在同一張圖上，并作比較，分析這兩種方法的精度。(1)向前歐拉法>>euler(0,10,100,10,i nlin e('y-20','x','y'),'Dy=y-20','y(0)=10')ans=精確解為s=20-10*exp(x)ans=1.0e+005*Colum ns1through80.00010.00010.00010.00010.00010.00000.00000.0000Colum

13、ns9through16-0.0000-0.0000-0.0001-0.0001-0.0001-0.0001-0.0002 -0.0002 Columns17through24 -0.0003-0.0003-0.0004-0.0004-0.0005-0.0005-0.0006 -0.0007Columns25through32 -0.0008-0.0009-0.0010-0.0011-0.0012-0.0014-0.0015 -0.0017Columns33through40 -0.0019-0.0021-0.0024-0.0026-0.0029-0.0032-0.0035 -0.0039Co

14、lumns41through48 -0.0043-0.0048-0.0053-0.0058-0.0064-0.0071-0.0078 -0.0086Columns49through56 -0.0095-0.0105-0.0115-0.0127-0.0140-0.0154-0.0170 -0.0187Columns57through64 -0.0206-0.0227-0.0250-0.0275-0.0302-0.0333-0.0366 -0.0403Columns65through72 -0.0444-0.0488-0.0537-0.0591-0.0651-0.0716-0.0788 -0.08

15、67Columns73through80 -0.0954-0.1049-0.1154-0.1270-0.1397-0.1537-0.1691 -0.1860Columns81through88 -0.2046-0.2251-0.2477-0.2724-0.2997-0.3297-0.3627 -0.3990Columns89through96 -0.4389-0.4828-0.5311-0.5842-0.6427-0.7070-0.7777 -0.8555Columns97through101 -0.9410-1.0352-1.1387-1.2526-1.3779 (2)改進歐拉法 >&

16、gt;eulerpro(0,10,100,10,inline('y-20','x','y'),'Dy=y-20','y(0)=10') ans=精確解為s=20-10*exp(x)ans=1.0e+005*Columns1through80.00010.00010.00010.00010.00010.00000.0000-0.0000Columns9through16-0.0000-0.0000-0.0001-0.0001-0.0001-0.0002-0.0002-0.0002Columns17through24-

17、0.0003-0.0003-0.0004-0.0005-0.0005-0.0006-0.0007-0.0008Colum ns25through32-0.0009-0.0010-0.0011-0.0013-0.0014-0.0016-0.0018-0.0020Colum ns33through40-0.0022-0.0025-0.0028-0.0031-0.0034-0.0038-0.0042-0.0047Colum ns41through48-0.0052-0.0058-0.0064-0.0071-0.0079-0.0087-0.0097-0.0107Colum ns49through56-

18、0.0119-0.0131-0.0145-0.0161-0.0178-0.0197-0.0218-0.0241Colum ns57through64-0.0266-0.0294-0.0325-0.0360-0.0398-0.0440-0.0486-0.0537Colum ns65through72-0.0594-0.0656-0.0726-0.0802-0.0886-0.0980-0.1083-0.1197Colum ns73through80-0.1323-0.1462-0.1615-0.1785-0.1973-0.2180-0.2409-0.2663Colum ns81through88-

19、0.2942-0.3251-0.3593-0.3971-0.4388-0.4849-0.5358-0.5921Colum ns89through96-0.6543-0.7230-0.7989-0.8828-0.9755-1.0780-1.1912-1.3163Colum ns97through101-1.4545-1.6073-1.7760-1.9626-2.1686改進歐拉法的精度比向前歐拉法更高。2-3分析應用題用以下三種不同的方法求下述微分方程的數值解，取h =10畫岀解的圖形，與精確值比較并進行分析。1) 歐拉法；2) 改進歐拉法；2-4分析應用題考慮一個涉及到社會上與眾不同的人的繁衍

20、問題模型。假設在時刻t(單位為年)，社會上有人口 X(t)人，又假設所有與眾不同的人與別的與眾不同的人結婚后所生后代也是與眾不同的人。而固 定比例為r的所有其他的后代也是與眾不同的人。如果對所有人來說岀生率假定為常數b，又如果普通的人和與眾不同的人的婚配是任意的，則此問題可以用微分方程表示為：其中變量p(t) =Xj(t”x(t)表示在時刻t社會上與眾不同的人的比例，x(t)表示在時刻t人口中與眾不同的人的數量。1) 假定p(0) =0.01,b =0.02和r =0.1，當步長為h =1年時，求從t =0到t =50解p(t)的近似值, 并作岀近似解的曲線圖形。2) 精確求岀微分方程的解p(t)，并將你當t =50時在分題(b)中得到的結果與此時的精確值進行比 較。1)>>euler(0,50,50,0.01,i nlin e('0.002-0.002*p','t','p