實驗報告七常微分方程初值問題的數(shù)值解法

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

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

3、與精確值比較并進行分析。1)歐拉法；2)改進歐拉法；3)龍格庫塔方法；2-4分析應用題考慮一個涉及到社會上與眾不同的人的繁衍問題模型。假設在時刻t (單位為年)，社會上有人口 x(t)人，又假設所有與眾不同的人與別的與眾不同的人結婚后所生后代也是與眾 不同的人。而固定比例為r的所有其他的后代也是與眾不同的人。如果對所有人來說出生率假定為常數(shù)b,又如果普通的人和與眾不同的人的婚配是任意的，則此問題可以用微分方程表不為：dp(t) rb(1 p(t) dt其中變量p(t) Xi(t)/x(t)表示在時刻t社會上與眾不同的人的比例，X(t)表示在時刻t人口中與眾不同的人的數(shù)量。1)假定p(0) 0.

4、01,b 0.02和r 0.1,當步長為h 1年時，求從t 0到t 50解p(t)的近 似值，并作出近似解的曲線圖形。2)精確求出微分方程的解p(t),并將你當t 50時在分題(b)中得到的結果與此時的精確值進行比較。【MATLAB相關函數(shù)求微分方程的解析解及其數(shù)值的代入dsolve( egn1egn2L x ) subs (expr, x,y, x1,y1,)其中e'gni '表示第i個方程，x '表示微分方程中的自變量，默認時自變量為t。subs命令中的expr、x、y為符合型表達式，x、y分別用數(shù)值x1、x2代入。>> syms x y z>&g

5、t; subs('x+y+z',x,y,z,1,2,3)ans =6>> syms x>> subs('xA2',x,2)ans =4ans =tan(x 1 4 pi)>> syms x>> subs(s,x,2)ans =-0.3721右端函數(shù)f(x, y)的自動生成f= inline( expr, varl' Var2；)其中'expr'表示函數(shù)的表達式，van； v'ar2'表示函數(shù)表達式中的變量，運行該函數(shù)，生 成一個新的函數(shù)表達式為f (varl, var2, )

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

7、對誤差10 3,絕對誤差10 6),程序為：options=odeset( reltol ；rt,'abstol',at),這里rt,at分別為設定的相對誤差和絕對誤差。常用 選項見下表。選項名功能可選值省缺值AbsTol設定絕對誤差正數(shù)1e 6RelTol設定相對誤差正數(shù)1e 3InitialStep設定初始步長正數(shù)自動MaxStep僅7E J K-t界正數(shù)tspan/10MaxOrder設定ode15s的最高階數(shù)1,2,3,4,515Stats顯示計算成本統(tǒng)計on,offoffBDF設定ode15s是否用反向差分on,offoff例：解微分方程y y 2-0 t 4yy(0

8、) 1在命令窗口執(zhí)行odefun = inline ( y 2*t/y： t, y t,y ode45(odefun,0,4,1)；t,y01.00000.05021.04900.10051.09590.15071.14083.85072.95033.90052.96723.95022.98394.00003.0006plop ,y ,o-',)ode45(odefun,0,4,1)t,y ode45(odefun,0 : 4,1)；ans =%解函數(shù)圖形表示%不用輸出變量，則直接輸出圖形t,yans =01.00001.00001.73212.00002.23613.00002.64

9、584.00003.0006.操作方法與實驗步驟(包括實驗數(shù)據(jù)記錄和處理)2-1編程編寫用向前歐拉公式和改進歐拉公式求微分方程數(shù)值解的Matlab程序，問題如下：在區(qū)間a,b內(nèi)(N 1)個等距點處，逼近下列初值問題的解，并對程序的每一句添上注釋語 句。y f(x, y) a x b y yEuler 法 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,y0,f,f1,b1) y=zeros(1,n+1);y(1)=y0;h=(b-a)/n;x=a:h:b;for i=1:n

10、;y(i+1)=y(i)+h*f(x(i),y(i);endplot(x,y)hold on% 求微分方程的精確解x1=linspace(a,b,100);'精確解為 's=dsolve(f1,b1,'x')syms xy1=zeros(1,100);fori=1:100y1(i)=subs(s,x,x1(i);endplot(x1,y1,'r')title(' 紅色代表精確解')改進Euler法y=eulerpro(a,b,n,y0,f,f1,b1)% 求微分方程的數(shù)值解y=zeros(1,n+1);y(1)=y0;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=linspace(a,b,100);'精確解為's=dsolve(f1,b1,'x')syms xy1=zeros(1,100);for i=1:100y1(i)=subs(s,x,x1(i);endplot(x1,y1,'r')title(' 紅色代表精確解')2-2 分析應用題假設等分區(qū)間數(shù)n

12、 100 ，用歐拉法和改進歐拉法在區(qū)間 t 0,10 內(nèi)求解初值問題y (t) y(t) 20y(0) 10并作出解的曲線圖形， 同時將方程的解析解也畫在同一張圖上， 并作比較， 分析這兩種方法 的精度。(1)向前歐拉法>> euler(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 *Columns 1 through 80.00010.00010.00010.

13、00010.00010.00000.00000.0000Columns 9 through 16-0.0000-0.0000-0.0001-0.0002Columns 17 through 24-0.0003-0.0003-0.0004-0.0007Columns 25 through 32-0.0008-0.0009-0.0010-0.0017Columns 33 through 40-0.0019-0.0021-0.0024-0.0039Columns 41 through 48-0.0043-0.0048-0.0053-0.0086Columns 49 through 56-0.0095

14、-0.0105-0.0115-0.0187Columns 57 through 64-0.0001-0.0001-0.0001-0.0002-0.0004-0.0005-0.0005-0.0006-0.0011-0.0012-0.0014-0.0015-0.0026-0.0029-0.0032-0.0035-0.0058-0.0064-0.0071-0.0078-0.0127-0.0140-0.0154-0.0170-0.0206-0.0227-0.0250-0.0275-0.0302-0.0333-0.0366-0.0403Columns 65 through 72-0.0444-0.048

15、8-0.0537-0.0591-0.0651-0.0716-0.0788-0.0867Columns 73 through 80-0.0954-0.1049-0.1154-0.1270-0.1397-0.1537-0.1691-0.1860Columns 81 through 88-0.2046-0.2251-0.2477-0.2724-0.2997-0.3297-0.3627-0.3990Columns 89 through 96-0.4389-0.4828-0.5311-0.5842-0.6427-0.7070-0.7777-0.8555Columns 97 through 101-0.9

16、410-1.0352-1.1387-1.2526-1.3779(2)改進歐拉法>> eulerpro(0,10,100,10,inline('y-20','x',y),'Dy=y-20','y(0)=10') ans =精確解為20 - 10*exp(x)ans =1.0e+005 *Columns 1 through 80.00010.00010.00010.00010.00010.00000.0000-0.0000Columns 9 through 16-0.0000-0.0000-0.0001-0.0001-0.

17、0001-0.0002-0.0002-0.0002Columns 17 through 24-0.0003-0.0003-0.0004-0.0005-0.0005-0.0006-0.0007-0.0008Columns 25 through 32-0.0009-0.0010-0.0011-0.0013-0.0014-0.0016-0.0018-0.0020Columns 33 through 40-0.0022-0.0025-0.0028-0.0031-0.0034-0.0038-0.0042-0.0047Columns 41 through 48-0.0052-0.0058-0.0064-0

18、.0071-0.0079-0.0087-0.0097-0.0107Columns 49 through 56-0.0119-0.0131-0.0145-0.0161-0.0178-0.0197-0.0218-0.0241Columns 57 through 64-0.0266-0.0294-0.0325-0.0360-0.0398-0.0440-0.0486-0.0537Columns 65 through 72-0.0594-0.0656-0.0726-0.0802-0.0886-0.0980-0.1083-0.1197Columns 73 through 80-0.1323-0.1462-

19、0.1615-0.1785-0.1973-0.2180-0.2409-0.2663Columns 81 through 88-0.2942-0.3251-0.3593-0.3971-0.4388-0.4849-0.5358-0.5921Columns 89 through 96-0.6543-0.7230-0.7989-0.8828-0.9755-1.0780 -1.1912-1.3163Columns 97 through 101-1.4545-1.6073-1.7760-1.9626-2.1686改進歐拉法的精度比向前歐拉法更高。2-3分析應用題用以下三種不同的方法求下述微分方程的數(shù)值解，

20、取h 10y y 2x y(0) 1畫出解的圖形，與精確值比較并進行分析。1)歐拉法；2)改進歐拉法； 2-4分析應用題考慮一個涉及到社會上與眾不同的人的繁衍問題模型。假設在時刻t (單位為年)，社會上有人口 x(t)人，又假設所有與眾不同的人與別的與眾不同的人結婚后所生后代也是與眾 不同的人。而固定比例為r的所有其他的后代也是與眾不同的人。如果對所有人來說出生率假定為常數(shù)b,又如果普通的人和與眾不同的人的婚配是任意的，則此問題可以用微分方程表示為：dp(t) rb(1 P(t) dt其中變量p(t) x(t)/x(t)表示在時刻t社會上與眾不同的人的比例，X(t)表示在時刻t人口中與眾不同的

21、人的數(shù)量。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)中得到的結果與此時的精確值進行比較。D >>euler(0,50,50,0.01,inline('0.002-0.002*p','t','p'),'Dp=0.002-0.002*p','p(0)=0.0 1') ans =精確解為1 - 99/(100*exp(x/ 500)ans =Columns 1 through 80.01000.01200.01400.01590.01790.01990.02180.0