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

浙江大學(xué)城市學(xué)院實驗報告課程名稱數(shù)值計算方法實驗項目名稱常微分方程初值問題的數(shù)值解法實驗成績指導(dǎo)老師（簽名）日期2015/12/16一.實驗?zāi)康暮鸵笥肕atlab軟件掌握求微分方程數(shù)值解的歐拉方法和龍格－庫塔方法；通過實例學(xué)習(xí)用微分方程模型解決簡化的實際問題。二.實驗內(nèi)容和原理編程題2-1要求寫出Matlab源程序(m文件)，并有適當?shù)淖⑨屨Z句；分析應(yīng)用題2-2，2-3，2-4，2-5要求將問題的分析過程、Matlab源程序和運行結(jié)果和結(jié)果的解釋、算法的分析寫在實驗報告上。編程編寫用向前歐拉公式和改進歐拉公式求微分方程數(shù)值解的Matlab程序，問題如下：在區(qū)間內(nèi)個等距點處，逼近下列初值問題的解，并對程序的每一句添上注釋語句。Euler法y=euler(a,b,n,y0,f,f1,b1)改進Euler法y=eulerpro(a,b,n,y0,f,f1,b1)分析應(yīng)用題假設(shè)等分區(qū)間數(shù)，用歐拉法和改進歐拉法在區(qū)間內(nèi)求解初值問題并作出解的曲線圖形，同時將方程的解析解也畫在同一張圖上，并作比較，分析這兩種方法的精度。分析應(yīng)用題用以下三種不同的方法求下述微分方程的數(shù)值解，取畫出解的圖形，與精確值比較并進行分析。1）歐拉法；2）改進歐拉法；3）龍格－庫塔方法；分析應(yīng)用題考慮一個涉及到社會上與眾不同的人的繁衍問題模型。假設(shè)在時刻(單位為年)，社會上有人口人，又假設(shè)所有與眾不同的人與別的與眾不同的人結(jié)婚后所生后代也是與眾不同的人。而固定比例為的所有其他的后代也是與眾不同的人。如果對所有人來說出生率假定為常數(shù)，又如果普通的人和與眾不同的人的婚配是任意的，則此問題可以用微分方程表示為：其中變量表示在時刻社會上與眾不同的人的比例，表示在時刻人口中與眾不同的人的數(shù)量。1）假定和，當步長為年時，求從到解的近似值，并作出近似解的曲線圖形。2）精確求出微分方程的解，并將你當時在分題(b)中得到的結(jié)果與此時的精確值進行比較?！綧ATLAB相關(guān)函數(shù)】求微分方程的解析解及其數(shù)值的代入dsolve(‘egn1’，‘egn2’，‘’)subs(expr,{x,y,…},{x1,y1,…})其中‘egn’表示第個方程，‘’表示微分方程中的自變量，默認時自變量為。subs命令中的expr、x、y為符合型表達式，x、y分別用數(shù)值x1、x2代入。>>symsxyz>>subs('x+y+z',{x,y,z},{1,2,3})ans=6>>symsx>>subs('x^2',x,2)ans=4s=dsolve(‘’,‘’,‘’)ans=>>symsx>>subs(s,x,2)ans=-0.3721右端函數(shù)的自動生成f=inline(‘expr’,’var1’,‘var2’,其中’expr’表示函數(shù)的表達式，’var1’,‘var2’表示函數(shù)表達式中的變量，運行該函數(shù)，生成一個新的函數(shù)表達式為f(var1,var2,>>f=inline('x+3*y','x','y')f=Inlinefunction: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ù)值(列向量)，t的步長是程序根據(jù)誤差限自動選定的。若ts=[t0,t1,t2,…,tf]，則輸出在自變量指定值，等步長時用ts=t0:k:tf，輸出在等分點；options用于設(shè)定誤差限(可以缺省，缺省時設(shè)定為相對誤差，絕對誤差)，程序為：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=01.00000.05021.04900.10051.09590.15071.14083.85072.95033.90052.96723.95022.98394.00003.0006plot(,,‘o-’,)%解函數(shù)圖形表示%不用輸出變量，則直接輸出圖形ans=01.00001.00001.73212.00002.23613.00002.64584.00003.0006三.操作方法與實驗步驟（包括實驗數(shù)據(jù)記錄和處理）2-1編程編寫用向前歐拉公式和改進歐拉公式求微分方程數(shù)值解的Matlab程序，問題如下：在區(qū)間內(nèi)個等距點處，逼近下列初值問題的解，并對程序的每一句添上注釋語句。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,y0,f,f1,b1)y=zeros(1,n+1);y(1)=y0;h=(b-a)/n;x=a:h:b;fori=1:n;y(i+1)=y(i)+h*f(x(i),y(i));endplot(x,y)holdon%求微分方程的精確解x1=linspace(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,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)holdon%求微分方程的精確解x1=linspace(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分析應(yīng)用題假設(shè)等分區(qū)間數(shù)，用歐拉法和改進歐拉法在區(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*Columns1through80.00010.00010.00010.00010.00010.00000.00000.0000Columns9through16-0.0000-0.0000-0.0001-0.0001-0.0001-0.0001-0.0002-0.0002Columns17through24-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.0039Columns41through48-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.0867Columns73through80-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改進歐拉法>>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-0.0003-0.0003-0.0004-0.0005-0.0005-0.0006-0.0007-0.0008Columns25through32-0.0009-0.0010-0.0011-0.0013-0.0014-0.0016-0.0018-0.0020Columns33through40-0.0022-0.0025-0.0028-0.0031-0.0034-0.0038-0.0042-0.0047Columns41through48-0.0052-0.0058-0.0064-0.0071-0.0079-0.0087-0.0097-0.0107Columns49through56-0.0119-0.0131-0.0145-0.0161-0.0178-0.0197-0.0218-0.0241Columns57through64-0.0266-0.0294-0.0325-0.0360-0.0398-0.0440-0.0486-0.0537Columns65through72-0.0594-0.0656-0.0726-0.0802-0.0886-0.0980-0.1083-0.1197Columns73through80-0.1323-0.1462-0.1615-0.1785-0.1973-0.2180-0.2409-0.2663Columns81through88-0.2942-0.3251-0.3593-0.3971-0.4388-0.4849-0.5358-0.5921Columns89through96-0.6543-0.7230-0.7989-0.8828-0.9755-1.0780-1.1912-1.3163Columns97through101-1.4545-1.6073-1.7760-1.9626-2.1686

改進歐拉法的精度比向前歐拉法更高。2-3分析應(yīng)用題用以下三種不同的方法求下述微分方程的數(shù)值解，取畫出解的圖形，與精確值比較并進行分析。1）歐拉法；2）改進歐拉法；2-4分析應(yīng)用題考慮一個涉及到社會上與眾不同的人的繁衍問題模型。假設(shè)在時刻(單位為年)，社會上有人口人，又假設(shè)所有與眾不同的人與別的與眾不同的人結(jié)婚后所生后代也是與眾不同的人。而固定比例為的所有其他的后代也是與眾不同的人。如果對所有人來說出生率假定為常數(shù)，又如果普通的人和與眾不同的人的婚配是任意的，則此問題可以用微分方程表示為：其中變量表示在時刻社會上與眾不同的人的比例，表示在時刻人口中與眾不同的人的數(shù)量。1）假定和，當步長為年時，求從到解的近似值，并作出近似解的曲線圖形。2）精確求出微分方程的解，并將你當時在分題(b)中得到的結(jié)果與此時的精確值進行比較。

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=Columns1through80.01000.01200.01400.0159