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

1、課程名稱數(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程序，問題如下：在區(qū)間a,b內(nèi)(N 1)

2、個等距點處，逼近下列初值問題的解，并對程序的每一句添上注釋語 句。y f (x, y) a x b y(a) y°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ù) n 100，用歐拉法和改進(jìn)歐拉法在區(qū)間t 0,10內(nèi)求解初值問題y(t) y(t) 20y(0)10并作出解的曲線圖形，同時將方程的解析解也畫在同一張圖上，并作比較，分析這兩種方法的精度。2-3分析應(yīng)用題用以下三種不同的方法求下述微分方程的數(shù)值解，取h 10y y 2xy(0) 1畫出解的圖形，與精

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

4、)0.01,b 0.02和r 0.1，當(dāng)步長為h 1年時，求從t 0到t 50解p(t)的近似值，并作出近似解的曲線圖形。2) 精確求出微分方程的解p(t),并將你當(dāng)t 50時在分題(b)中得到的結(jié)果與此時的精確值 進(jìn)行比較?！綧ATLAB相關(guān)函數(shù)】求微分方程的解析解及其數(shù)值的代入dsolve( subs (expr, x,y, egn1', egn2', L x'), x1,y1,)其中egni '表示第i個方程， x '表示微分方程中的自變量，默認(rèn)時自變量為t。subs命令中的expr、x、y為符合型表達(dá)式，x、y分別用數(shù)值x1、x2代入。>&

5、gt; syms x y z>> subs('x+y+z',x,y,z,1,2,3)ans =6>> syms x>> subs('xA2',x,2) ans =s=dsolve( Dy 1 y 2' , y(0)1' , x')ans =tan(x 1 4 pi)4>> syms x>> subs(s,x,2) ans =右端函數(shù)f (x, y)的自動生成f= in li ne('expr' , ' var1 ' ,' var2 '

6、;,)其中expr '表示函數(shù)的表達(dá)式，varl ' ,var2 '表示函數(shù)表達(dá)式中的變量，運行該函數(shù)，生成一個新的函數(shù)表達(dá)式為 f (varl, var2, )。>> f=i nlin e('x+3*y','x','y')Inline fun cti on: f(x,y) = x+3*y>> f(2,3) ans =114，5階龍格-庫塔方法求解微分方程數(shù)值解t,x=ode45(f,ts,x0,optio ns)其中f是由待解方程寫成的m文件名；x0為函數(shù)的初值；t,x分別為輸出的自變量和函數(shù)值(

7、列向量)，t的步長是程序根據(jù)誤差限自動選定的。若ts=t0,t1,t2,tf，則輸出在自變量指定值，等步長時用ts=tO:k:tf ,輸出在等分點；options用于設(shè)定誤差限(可以缺省，缺省時設(shè)定為相對誤差10 3，絕對誤差10 6)，程序為：options=odeset( reltol ' ,rt, ' abstol ' ,at)，這里 rt,at分別為設(shè)定的相對誤差和絕對誤差。常用選項見下表。選項名功能可選值省缺值A(chǔ)bsTol設(shè)定絕對誤差正數(shù)1e 6RelTol設(shè)定相對誤差正數(shù)1e 3In itialStep設(shè)疋初始步長正數(shù)自動MaxStep2 rc 宀 rH x

8、 r. m 設(shè)定步長上界正數(shù)tspan/10MaxOrder設(shè)定ode15s的最高階數(shù)1,2,3,4,55Stats顯示計算成本統(tǒng)計on,offoffBDF設(shè)定ode15s是否用反向差分on,offoff例：解微分方程y y 20 t 4yy(0) 1在命令窗口執(zhí)行odefun = inline ( ' y 2* t y ' ,' t' ,' y ');t,y ode45(odefun,0,4,1) ；t,y ans =0plot( t, y, o- ' ,)% 解函數(shù)圖形表示ode45(odefun,0,4,1) % 不用輸出變量，則直

9、接輸出圖形 t,y ode45(odefun,0 : 4,1) ； t,yans =0操作方法與實驗步驟(包括實驗數(shù)據(jù)記錄和處理)2-1 編程編寫用向前歐拉公式和改進(jìn)歐拉公式求微分方程數(shù)值解的 Matlab 程序，問題如下： 在區(qū)間 a,b 內(nèi) (N 1)個等距點處，逼近下列初值問題的解，并對程序的每一句添上注釋語 句。y f(x, y) a x b y(a) y°Euler 法 y=euler(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=ze

10、ros(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 x y1=zeros(1,100);fori=1:100y1(i)=subs(s,x,x1(i);endplot(x1,y1,'r')title(' 紅色代表精確解 ')改進(jìn) Euler 法y=eulerpro(a,

11、b,n,y0,f,f1,b1) % 求微分方程的數(shù)值解 y=zeros(1,n+1);y(1)=y0;h=(b-a)/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 x y1=zeros(1,100);for i=1:100y1(i)=subs(s,x,x1(i);endplo

12、t(x1,y1,'r')title(' 紅色代表精確解 ')2-2 分析應(yīng)用題0,10 內(nèi)求解初值問題并作比較， 分析這兩種方法假設(shè)等分區(qū)間數(shù) n 100 ，用歐拉法和改進(jìn)歐拉法在區(qū)間 ty (t) y(t) 20 y(0) 10并作出解的曲線圖形， 同時將方程的解析解也畫在同一張圖上， 的精度。(1) 向前歐拉法>> euler(0,10,100,10,inline('y-20','x','y'),'Dy=y-20','y(0)=10') ans =精確解為s =20

13、- 10*exp(x)ans =Columns 1 through 8+005 *Columns 9 through 16Columns 17 through 24Columns 25 through 32Columns 33 through 40Columns 41 through 48Columns 49 through 56Columns 57 through 64Columns 65 through 72Columns 73 through 80Columns 81 through 88Columns 89 through 96Columns 97 through 101(2)改進(jìn)歐拉

14、法>> eulerpro(0,10,100,10,inline('y-20','x','y'),'Dy=y-20','y(0)=10') ans =精確解為s = 20 - 10*exp(x)ans = +005 *Columns 1 through 8Columns 9 through 16Columns 17 through 24Columns 25 through 32Columns 33 through 40Columns 41 through 48Columns 49 through 56Co

15、lumns 57 through 64Columns 65 through 72Columns 73 through 80Columns 81 through 88Columns 89 through 96Columns 97 through 101改進(jìn)歐拉法的精度比向前歐拉法更高。2-3分析應(yīng)用題用以下三種不同的方法求下述微分方程的數(shù)值解，取h 10y y 2x 0 x 1 y(0) 1畫出解的圖形，與精確值比較并進(jìn)行分析。1) 歐拉法;2) 改進(jìn)歐拉法；2-4分析應(yīng)用題考慮一個涉及到社會上與眾不同的人的繁衍問題模型。假設(shè)在時刻t(單位為年)，社會上有人口 x(t)人，又假設(shè)所有與眾不同的人

16、與別的與眾不同的人結(jié)婚后所生后代也是與眾 不同的人。而固定比例為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.01,b 0.02和r 0.1，當(dāng)步長為h 1年時，求從t 0到t 50解p(t)的近 似值，并作出近似解的曲線圖形。2) 精確求出微分方程的解p(t)，并將你當(dāng)t 50時在分題(b)中得到的結(jié)果與此時的精確值 進(jìn)行比較。1)>>euler(0,50,50,i nlin e(”,'t','p'),'Dp=','p(0)=1')ans =精確解為1 - 99/(100*exp(x/500) ans =Columns 1 through 8Columns 9 t