實驗報告七 常微分方程初值問題的數(shù)值解法[教育相關(guān)]

1、浙江大學(xué)城市學(xué)院實驗報告課程名稱 數(shù)值計算方法 實驗項目名稱 常微分方程初值問題的數(shù)值解法 實驗成績 指導(dǎo)老師（簽名 ） 日期 2015/12/16 一. 實驗?zāi)康暮鸵? 用Matlab軟件掌握求微分方程數(shù)值解的歐拉方法和龍格庫塔方法；2 通過實例學(xué)習(xí)用微分方程模型解決簡化的實際問題。二. 實驗內(nèi)容和原理編程題2-1要求寫出Matlab源程序(m文件)，并有適當(dāng)?shù)淖⑨屨Z句；分析應(yīng)用題2-2，2-3，2-4，2-5要求將問題的分析過程、Matlab源程序和運行結(jié)果和結(jié)果的解釋、算法的分析寫在實驗報告上。2-1 編程編寫用向前歐拉公式和改進(jìn)歐拉公式求微分方程數(shù)值解的Matlab程序，問題如下：在

2、區(qū)間內(nèi)個等距點處，逼近下列初值問題的解，并對程序的每一句添上注釋語句。 Euler法 y=euler(a,b,n,y0,f,f1,b1)改進(jìn)Euler法 y=eulerpro(a,b,n,y0,f,f1,b1)2-2 分析應(yīng)用題假設(shè)等分區(qū)間數(shù)，用歐拉法和改進(jìn)歐拉法在區(qū)間內(nèi)求解初值問題并作出解的曲線圖形，同時將方程的解析解也畫在同一張圖上，并作比較，分析這兩種方法的精度。2-3 分析應(yīng)用題用以下三種不同的方法求下述微分方程的數(shù)值解，取 畫出解的圖形，與精確值比較并進(jìn)行分析。 1）歐拉法； 2）改進(jìn)歐拉法；3）龍格庫塔方法； 2-4 分析應(yīng)用題考慮一個涉及到社會上與眾不同的人的繁衍問題模型。假設(shè)在

3、時刻(單位為年)，社會上有人口人，又假設(shè)所有與眾不同的人與別的與眾不同的人結(jié)婚后所生后代也是與眾不同的人。而固定比例為的所有其他的后代也是與眾不同的人。如果對所有人來說出生率假定為常數(shù)，又如果普通的人和與眾不同的人的婚配是任意的，則此問題可以用微分方程表示為： 其中變量表示在時刻社會上與眾不同的人的比例，表示在時刻人口中與眾不同的人的數(shù)量。 1）假定和，當(dāng)步長為年時，求從到解的近似值，并作出近似解的曲線圖形。2）精確求出微分方程的解，并將你當(dāng)時在分題(b)中得到的結(jié)果與此時的精確值進(jìn)行比較。【MATLAB相關(guān)函數(shù)】n 求微分方程的解析解及其數(shù)值的代入 dsolve(egn1， egn2， )

4、subs (expr, x,y, x1,y1, )其中egn表示第個方程，表示微分方程中的自變量，默認(rèn)時自變量為。subs命令中的expr、x、y為符合型表達(dá)式，x、y分別用數(shù)值x1、x2代入。 syms x y z subs(x+y+z,x,y,z,1,2,3)ans = 6 syms x subs(x2,x,2)ans = 4 s=dsolve(, , )ans = syms x subs(s,x,2)ans = -0.3721n 右端函數(shù)的自動生成 f= inline(expr, var1, var2,) 其中expr表示函數(shù)的表達(dá)式，var1, var2 表示函數(shù)表達(dá)式中的變量，運行該

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

6、差)，程序為：options=odeset(reltol,rt,abstol,at)，這里rt,at分別為設(shè)定的相對誤差和絕對誤差。常用選項見下表。選項名功能可選值省缺值A(chǔ)bsTol設(shè)定絕對誤差正數(shù)RelTol設(shè)定相對誤差正數(shù)InitialStep設(shè)定初始步長正數(shù)自動MaxStep設(shè)定步長上界正數(shù)MaxOrder設(shè)定ode15s的最高階數(shù)1,2,3,4,55Stats顯示計算成本統(tǒng)計on,offoffBDF設(shè)定ode15s是否用反向差分on,offoff例：解微分方程 在命令窗口執(zhí)行 = (, , )； ； ans = 0 1.0000 0.0502 1.0490 0.1005 1.0959

7、0.1507 1.1408 3.8507 2.9503 3.9005 2.9672 3.9502 2.9839 4.0000 3.0006 plot(,o-,) %解函數(shù)圖形表示 %不用輸出變量，則直接輸出圖形 ； ans = 0 1.0000 1.0000 1.7321 2.0000 2.2361 3.0000 2.6458 4.0000 3.0006三. 操作方法與實驗步驟（包括實驗數(shù)據(jù)記錄和處理）2-1 編程編寫用向前歐拉公式和改進(jìn)歐拉公式求微分方程數(shù)值解的Matlab程序，問題如下：在區(qū)間內(nèi)個等距點處，逼近下列初值問題的解，并對程序的每一句添上注釋語句。 Euler法 y=euler(

8、a,b,n,y0,f,f1,b1)改進(jìn)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; y(i+1)=y(i)+h*f(x(i),y(i); end plot(x,y) hold on % 求微分方程的精確解 x1=linspace(a,b,100); 精確解為 s=dsolve(f1,b1,x) syms xy1=zeros(1,100); for i=1:100 y1(i)=subs(s,x,x

9、1(i); end plot(x1,y1,r) title(紅色代表精確解)改進(jìn)Euler法 y=eulerpro(a,b,n,y0,f,f1,b1) % 求微分方程的數(shù)值解 y=zeros(1,n+1); y(1)=y0; h=(b-a)/n; x=a:h:b; for i=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); end plot(x,y) hold on% 求微分方程的精確解 x1=linspace(a,b,100);精確解為 s=dsolve(f1,b1,x) syms x y1=ze

10、ros(1,100); for i=1:100 y1(i)=subs(s,x,x1(i); end plot(x1,y1,r)title(紅色代表精確解)2-2分析應(yīng)用題假設(shè)等分區(qū)間數(shù)，用歐拉法和改進(jìn)歐拉法在區(qū)間內(nèi)求解初值問題并作出解的曲線圖形，同時將方程的解析解也畫在同一張圖上，并作比較，分析這兩種方法的精度。(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 8 0.0001 0.0001 0.

11、0001 0.0001 0.0001 0.0000 0.0000 0.0000 Columns 9 through 16 -0.0000 -0.0000 -0.0001 -0.0001 -0.0001 -0.0001 -0.0002 -0.0002 Columns 17 through 24 -0.0003 -0.0003 -0.0004 -0.0004 -0.0005 -0.0005 -0.0006 -0.0007 Columns 25 through 32 -0.0008 -0.0009 -0.0010 -0.0011 -0.0012 -0.0014 -0.0015 -0.0017 Col

12、umns 33 through 40 -0.0019 -0.0021 -0.0024 -0.0026 -0.0029 -0.0032 -0.0035 -0.0039 Columns 41 through 48 -0.0043 -0.0048 -0.0053 -0.0058 -0.0064 -0.0071 -0.0078 -0.0086 Columns 49 through 56 -0.0095 -0.0105 -0.0115 -0.0127 -0.0140 -0.0154 -0.0170 -0.0187 Columns 57 through 64 -0.0206 -0.0227 -0.0250

13、 -0.0275 -0.0302 -0.0333 -0.0366 -0.0403 Columns 65 through 72 -0.0444 -0.0488 -0.0537 -0.0591 -0.0651 -0.0716 -0.0788 -0.0867 Columns 73 through 80 -0.0954 -0.1049 -0.1154 -0.1270 -0.1397 -0.1537 -0.1691 -0.1860 Columns 81 through 88 -0.2046 -0.2251 -0.2477 -0.2724 -0.2997 -0.3297 -0.3627 -0.3990 C

14、olumns 89 through 96 -0.4389 -0.4828 -0.5311 -0.5842 -0.6427 -0.7070 -0.7777 -0.8555 Columns 97 through 101 -0.9410 -1.0352 -1.1387 -1.2526 -1.3779（2） 改進(jìn)歐拉法 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 * Columns 1 through 8 0.0001 0.0001 0.0001

15、0.0001 0.0001 0.0000 0.0000 -0.0000 Columns 9 through 16 -0.0000 -0.0000 -0.0001 -0.0001 -0.0001 -0.0002 -0.0002 -0.0002 Columns 17 through 24 -0.0003 -0.0003 -0.0004 -0.0005 -0.0005 -0.0006 -0.0007 -0.0008 Columns 25 through 32 -0.0009 -0.0010 -0.0011 -0.0013 -0.0014 -0.0016 -0.0018 -0.0020 Columns

16、 33 through 40 -0.0022 -0.0025 -0.0028 -0.0031 -0.0034 -0.0038 -0.0042 -0.0047 Columns 41 through 48 -0.0052 -0.0058 -0.0064 -0.0071 -0.0079 -0.0087 -0.0097 -0.0107 Columns 49 through 56 -0.0119 -0.0131 -0.0145 -0.0161 -0.0178 -0.0197 -0.0218 -0.0241 Columns 57 through 64 -0.0266 -0.0294 -0.0325 -0.

17、0360 -0.0398 -0.0440 -0.0486 -0.0537 Columns 65 through 72 -0.0594 -0.0656 -0.0726 -0.0802 -0.0886 -0.0980 -0.1083 -0.1197 Columns 73 through 80 -0.1323 -0.1462 -0.1615 -0.1785 -0.1973 -0.2180 -0.2409 -0.2663 Columns 81 through 88 -0.2942 -0.3251 -0.3593 -0.3971 -0.4388 -0.4849 -0.5358 -0.5921 Colum

18、ns 89 through 96 -0.6543 -0.7230 -0.7989 -0.8828 -0.9755 -1.0780 -1.1912 -1.3163 Columns 97 through 101 -1.4545 -1.6073 -1.7760 -1.9626 -2.1686改進(jìn)歐拉法的精度比向前歐拉法更高。2-3分析應(yīng)用題用以下三種不同的方法求下述微分方程的數(shù)值解，取 畫出解的圖形，與精確值比較并進(jìn)行分析。 1）歐拉法； 2）改進(jìn)歐拉法；2-4分析應(yīng)用題考慮一個涉及到社會上與眾不同的人的繁衍問題模型。假設(shè)在時刻(單位為年)，社會上有人口人，又假設(shè)所有與眾不同的人與別的與眾不同的人結(jié)

19、婚后所生后代也是與眾不同的人。而固定比例為的所有其他的后代也是與眾不同的人。如果對所有人來說出生率假定為常數(shù)，又如果普通的人和與眾不同的人的婚配是任意的，則此問題可以用微分方程表示為： 其中變量表示在時刻社會上與眾不同的人的比例，表示在時刻人口中與眾不同的人的數(shù)量。 1）假定和，當(dāng)步長為年時，求從到解的近似值，并作出近似解的曲線圖形。2）精確求出微分方程的解，并將你當(dāng)時在分題(b)中得到的結(jié)果與此時的精確值進(jìn)行比較。1） euler(0,50,50,0.01,inline(0.002-0.002*p,t,p),Dp=0.002-0.002*p,p(0)=0.01) ans = 精確解為 s = 1 - 99/(100*exp(x/500) ans = Columns 1