實(shí)驗(yàn)一 常微分方程

1、年級(jí)、專業(yè) 2012級(jí)數(shù)學(xué)與應(yīng)用數(shù)學(xué) 姓名 張旭 學(xué)號(hào) 12160011063 名單序號(hào) 16 實(shí)驗(yàn)時(shí)間 2014年 3月 19 日 使用設(shè)備、軟件 PC, MATLAB 注： 實(shí)驗(yàn)報(bào)告的最后一部分是實(shí)驗(yàn)小結(jié)與收獲 實(shí)驗(yàn)一 常微分方程1. 分別用Euler法和ode45解下列常微分方程并與解析解比較： (1) 解析法：y=dsolve(Dy=x+y,y(0)=1,x)y = 2*exp(x) - x - 1Euler:function x,y=euler(odefun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;for i=1:length(x)-1 y(

2、k+1)=y(i)+h*feval(odefun,x(i),y(i); end x=x ; y=y ;endode45:odefun=inline(x+y,x,y);xspan=0,3;y0=1;h=0.1;x1,y1=euler(odefun,xspan,y0,h);x2,y2=ode45(odefun,xspan,y0);x3=0:0.1:3;y3=2*exp(x3)-x3-1;plot(x1,y1,k,x2,y2,ko,x3,y3,k*);xlabel(x軸);ylabel(y軸);legend(euler,ode45,dsolve);ode45求得的結(jié)果與用解析法求得的結(jié)果更接近，故o

3、de45的精度較高，Euler法求得的結(jié)果精度較低。 (2) 令 則原方程等價(jià)于方程組：， ，不能解析，只能用數(shù)值法求解。Euler：function t,y=euler2(odefun1,odefun2,tspan,y0,h)t=tspan(1):h:tspan(2);y(1,1)=y0(1);y(2,1)=y0(2);for i=1:length(t)-1 k1=odefun1(t(i),y(1,i),y(2,i); k2=odefun2(t(i),y(1,i),y(2,i); y(1,i+1)=y(1,i)+h*d1; y(2,i+1)=y(2,i)+h*d2; end t=t; y=y

4、;endode45： odefun1=inline(0*t1+0*y1+y2); odefun2=inline(-2*y1+0.01*y22+sin(t1); t1,y1=euler2(odefun1,odefun2,0,5,0,1,0.1); t2,y2=ode45(eu,0,5,0,1); plot(t1,y1(:,1),o,t2,y2(:,1),LineWidth,2); xlabel(t軸); ylabel(y軸); legend(euler,ode45);ode45中 eu：function dy=eu(t,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=-2*y(

5、1)+0.01*y(2)2+sin(t);ode45求得的結(jié)果精度較高，euler法求得的結(jié)果在準(zhǔn)確值上下波動(dòng)。2. 一通過(guò)原點(diǎn)的曲線，它在處的切線斜率等于若上限增為1.58，1.60會(huì)發(fā)生什么？ 等價(jià)于求解 ，且的初值問(wèn)題。解析法： y=dsolve(Dy=2*x+y2,y(0)=0,x)y =(2(1/3)*airy(3,-2(1/3)*x)+2(1/3)*3(1/2)*airy(1,-2(1/3)*x)/(3(1/2)*airy(0, -2(1/3)*x) + airy(2, -2(1/3)*x)ode45： odefun=inline(2*x+y2); subplot(1,3,1);o

6、de45(odefun,0,1.57,0);title(0x1.57); subplot(1,3,2);ode45(odefun,0,1.58,0);title(0x1.58); subplot(1,3,3);ode45(odefun,0,1.60,0);title(0x1.5之后的斜率增長(zhǎng)速度很快，若上限增為1.58，1.60，則相應(yīng)的y將會(huì)出現(xiàn)更大的增長(zhǎng)。3. 求解剛性方程組：function dy=fun(t,y)dy=zeros(2,1);dy(1)=-1000.25*y(1)+999.75*y(2)+0.5;dy(2)=999.75*y(1)-1000.25*y(2)+0.5;ode

7、45： t,y=ode15s(fun,0,50,1,-1); plot(t,y(:,1),o,t,y(:,2),k-,LineWidth,2); legend(y1,y2);4. (溫度過(guò)程)夏天把開(kāi)有空調(diào)的室內(nèi)一支讀數(shù)為20 的溫度計(jì)放到戶外，10分鐘后讀25.2 , 再過(guò)10分鐘后讀數(shù)28.32 。建立一個(gè)較合理的模型來(lái)推算戶外溫度。由題意可知由于隨著時(shí)間的增加，溫度增長(zhǎng)越來(lái)越慢，戶外溫度與溫度計(jì)示數(shù)之差也越來(lái)越小，且溫差為零時(shí)溫度的增長(zhǎng)率也為零，故可以認(rèn)為溫度的增長(zhǎng)率與溫差成正比，設(shè)戶外溫度為m，溫度計(jì)的示數(shù)為y，比例系數(shù)為k，則可建立模型 解析法：y=dsolve(Dy=k*(c-y)

8、,y(0)=20,t)y = m - (m - 20)/exp(k*t)由y（10）=25.2,y(20)=28.32建立方程組，消去k，得： （m-20）（m-28.32）=（m-25.2）（m-25.2）解得：m=33 所以戶外溫度約為33 。5. (廣告效應(yīng))某公司生產(chǎn)一種耐用消費(fèi)品，市場(chǎng)占有率為5%時(shí)開(kāi)始做廣告，一段時(shí)間的市場(chǎng)跟蹤調(diào)查后，該公司發(fā)現(xiàn)：?jiǎn)挝粫r(shí)間內(nèi)購(gòu)買(mǎi)人口百分比的相對(duì)增長(zhǎng)率與當(dāng)時(shí)還沒(méi)有買(mǎi)的百分比成正比，且估得此比例系數(shù)為0.5。(1) 建立該問(wèn)題的數(shù)學(xué)模型，并求其數(shù)值解與模擬結(jié)果作以比較；設(shè)t時(shí)刻該消費(fèi)品的市場(chǎng)占有率為y，建立方程：解析解：y=dsolve(Dy=0.5-0

9、.5*y,y(0)=0.05)y = 1 - (19*exp(-t/2)/20數(shù)值解： odefun=inline(0.5-0.5*y,t,y); t1,y1=ode45(odefun,0,10,0.05); t2=0:0.1:10; y2=1-(19*exp(-t2/2)/20; plot(t1,y1,o,t2,y2,k); legend(ode,dsolve);(2) 廠家問(wèn)：要做多少時(shí)間廣告，可使市場(chǎng)購(gòu)買(mǎi)率達(dá)到80%？由解析解可列出方程 1 - (19*exp(-t/2)/20=0.8，所以解得t=3.11636. (腫瘤生長(zhǎng)) 腫瘤大小V生長(zhǎng)的速率與V的a次方成正比，其中a為形狀參數(shù)，

10、0a1;而其比例系數(shù)K隨時(shí)間減小，減小速率又與當(dāng)時(shí)的K值成正比，比例系數(shù)為環(huán)境參數(shù)b。設(shè)某腫瘤參數(shù)a=1, b=0.1, K的初始值為2，V的初始值為1。問(wèn)(1) 此腫瘤生長(zhǎng)不會(huì)超過(guò)多大？由已知列出方程組 ，代入具體數(shù)值，得 ，解析法：k,v=dsolve(Dk=-0.1*k,Dv=k*v,k(0)=2,v(0)=1,t)k =2*exp(-t/10)v =exp(20)*exp(-20*exp(-t/10)t=0:0.1:100; v=exp(20)*exp(-20*exp(-t/10); plot(t,v,LineWidth,2); xlabel(t軸); ylabel(v軸);因腫瘤不斷

11、長(zhǎng)大，故t趨于無(wú)窮時(shí)，該腫瘤達(dá)到最大，此時(shí)極限為exp(20)=4.8517*108，故此腫瘤生長(zhǎng)不會(huì)超過(guò)4.8517*108 。(2) 過(guò)多長(zhǎng)時(shí)間腫瘤大小翻一倍？令exp(20)*exp(-20*exp(-t/10)=2，解得t=-10*ln(1-1/20*ln2)=0.3527,(3) 何時(shí)腫瘤生長(zhǎng)速率由遞增轉(zhuǎn)為遞減？由已求得的結(jié)果可得與的關(guān)系為=2*exp(20-t/10)*exp(-20*exp(-t/10)，t1=0:0.1:100; v1=2*exp(20-t1/10).*exp(-20*exp(-t1/10); plot(t1,v1,LineWidth,2); xlabel(t軸

12、); ylabel(v軸);顯然，最大值處對(duì)應(yīng)的t即為所求： t2=0:0.01:100; v2=2*exp(20-t1/10).*exp(-20*exp(-t1/10); m,n=max(v2); t=t2(n)得到t = 29.96 (4) 若參數(shù)a=2/3呢？1、建立方程組， 解析法：k,v=dsolve(Dk=-0.1*k,Dv=k*v(2/3),k(0)=2,v(0)=1,t)k = 2*exp(-t/10) 2*exp(-t/10) 2*exp(-t/10)v = -(20*exp(-t/10) - 23)3/27 (37/2 + (3(1/2)*3*i)/2 - 20*exp(-

13、t/10)3/27 -(20*exp(-t/10) + (3(1/2)*3*i)/2 - 37/2)3/27取實(shí)解k = 2*exp(-t/10) v=-(20*exp(-t/10) - 23)3/27并畫(huà)出v-t圖像： t=0:0.1:100; v=-(20*exp(-t/10) - 23).3/27; plot(t,v,LineWidth,2); xlabel(t軸); ylabel(v軸);顯然，當(dāng)t趨于無(wú)窮時(shí)，該腫瘤達(dá)到最大，此時(shí)極限為-(- 23)3/27=450.6296，故此腫瘤生長(zhǎng)不會(huì)超過(guò)450.6296 。2、令-(20*exp(-t/10) - 23)3/27=2，解得t=

14、0.3977,3、由已求得的結(jié)果可得與的關(guān)系為=2*exp(-t/10)*(20*exp(-t/10) - 23)2/9，t1=0:0.1:100; v1=2*exp(-t1/10).*(20*exp(-t1/10) - 23).2/9; plot(t1,v1,LineWidth,2); xlabel(t軸); ylabel(v軸);顯然，最大值處對(duì)應(yīng)的t即為所求，： t2=0:0.01:100; v2=2*exp(-t2/10).*(20*exp(-t2/10) - 23).2/9; m,n=max(v2); t=t2(n)t = 9.5900選做題：1.(生態(tài)系統(tǒng)的振蕩現(xiàn)象)第一次世界大戰(zhàn)

15、中，因?yàn)閼?zhàn)爭(zhēng)很少捕魚(yú)，按理戰(zhàn)后應(yīng)能捕到更多的魚(yú)才是?？墒谴髴?zhàn)后，在地中海卻捕不到鯊魚(yú)，因而漁民大惑不解。令x1為魚(yú)餌的數(shù)量，x2為鯊魚(yú)的數(shù)量，t為時(shí)間。微分方程為 式中a1, a2, b1, b2都是正常數(shù)。第一式魚(yú)餌x1的增長(zhǎng)速度大體上與x1成正比，即按a1x1比率增加, 而被鯊魚(yú)吃掉的部分按b1x1x2的比率減少；第二式中鯊魚(yú)的增長(zhǎng)速度由于生存競(jìng)爭(zhēng)的自然死亡或互相咬食按a2x2的比率減少，但又根據(jù)魚(yú)餌的量的變化按b2x1x2的比率增加。對(duì)a1=3, b1=2, a2=2.5, b2=1, x1(0)=x2(0)=1求解。畫(huà)出解曲線圖和相軌線圖，可以觀察到魚(yú)餌和鯊魚(yú)數(shù)量的周期振蕩現(xiàn)象。代入具體數(shù)值后，原方程組