春數(shù)學(xué)實(shí)驗(yàn)基礎(chǔ) 實(shí)驗(yàn)報(bào)告(1)常微分方程.doc

年級(jí)、專業(yè) 姓名 學(xué)號(hào) 名單序號(hào) 實(shí)驗(yàn)時(shí)間 2013年 3月 日 使用設(shè)備、軟件 PC, MATLAB 注： 實(shí)驗(yàn)報(bào)告的最后一部分是實(shí)驗(yàn)小結(jié)與收獲 實(shí)驗(yàn)一 常微分方程1. 分別用Euler法和ode45解下列常微分方程并與解析解比較： (1) 編寫Euler法的matlab函數(shù)，命名為euler.mfunction t,y=euler(odefun,tspan,y0,h)t=tspan(1):h:tspan(2);y(1)=y0;for i=1:length(t)-1y(i+1)=y(i)+h*feval(odefun,t(i),y(i);endt=t;y=y;下面比較三者的差別：% ode45odefun=inline(x+y,x,y);x1,y1=ode45(odefun,0,3,1);plot(x1,y1,ko);pausehold on;% Eulerx2,y2=euler(odefun,0,3,1,0.05);plot(x2,y2,r+);pausehold on;% 解析解y0=dsolve(Dy=t+y,y(0)=1);ezplot(y0,0,3);pausehold off;legend(ode45,euler法,解析解);Euler法只有一階精度，所以實(shí)際應(yīng)用效率比較差，而ode45的效果比較好，很接近真實(shí)值。 (2) 先寫M文件ex1_2fun.mfunction f=ex1_2fun(t,y)f(1)=y(2);f(2)=0.01*y(2).2-2*y(1)+sin(t);f=f(:);% ode45t1,y1=ode45(ex1_2fun,0,5,0;1);plot(t1,y1(:,1),ko);% 解析解s=dsolve(D2y-0.01*(Dy)2+2*y=sin(t),y(0)=0,Dy(0)=1,t)s = empty sym %由此可知該微分方程無解析解2. 求一通過原點(diǎn)的曲線，它在處的切線斜率等于若上限增為1.58，1.60會(huì)發(fā)生什么？ odefun=inline(2*x+y2,x,y);subplot(1,4,1);x1,y1=ode45(odefun,0,1.57,0);plot(x1,y1,r*);title(上限1.57);subplot(1,4,2);x2,y2=ode45(odefun,0,1.58,0);plot(x2,y2,bo);title(上限1.58);subplot(1,4,3);x3,y3=ode45(odefun,0,1.6,0);plot(x3,y3,k);title(上限1.60);subplot(1,4,4);plot(x1,y1,r*);hold on;plot(x2,y2,bo);hold on;plot(x3,y3,k);hold off;legend(上限1.57,上限1.58,上限1.60);結(jié)論：隨著x上界的增加，解趨于無窮大。3. 求解剛性方程組：先寫M函數(shù)ex3fun.mfunction f=ex3fun(t,y)f(1)=-1000.25*y(1)+999.75*y(2)+0.5;f(2)=999.75*y(1)-1000.25*y(2)+0.5;f=f(:);%作圖t,y=ode15s(ex3fun,0,50,1,-1);plot(t,y,*);4. (溫度過程)夏天把開有空調(diào)的室內(nèi)一支讀數(shù)為20的溫度計(jì)放到戶外，10分鐘后讀25.2, 再過10分鐘后讀數(shù)28.32。建立一個(gè)較合理的模型來推算戶外溫度。設(shè)：t時(shí)刻溫度計(jì)的讀數(shù)為T，戶外溫度為c，T的增速與室內(nèi)外溫差（c-T）成正比，由此建立微分方程，其中k為比例系數(shù)%首先，計(jì)算解析解 y=dsolve(DT=k*(c-T),T(0)=20,t)y =c - (c - 20)/exp(k*t)%又已知，用非線性最小二乘擬合該函數(shù)，調(diào)用lsqcurvefit命令：fun=inline(c(1)-(c(1)-20)./exp(c(2)*t),c,t);lsqcurvefit(fun,30 1,10 20,25.2 28.32)ans = 33.0000 0.0511即戶外溫度c=33，比例系數(shù)k=0.05115. (廣告效應(yīng))某公司生產(chǎn)一種耐用消費(fèi)品，市場占有率為5%時(shí)開始做廣告，一段時(shí)間的市場跟蹤調(diào)查后，該公司發(fā)現(xiàn)：單位時(shí)間內(nèi)購買人口百分比的相對(duì)增長率與當(dāng)時(shí)還沒有買的百分比成正比，且估得此比例系數(shù)為0.5。(1) 建立該問題的數(shù)學(xué)模型，并求其數(shù)值解與模擬結(jié)果作以比較；設(shè)：市場占有率為，市場占有率的增長速度為，則相對(duì)增長率為，由此建立微分方程為%首先，計(jì)算解析解 s=dsolve(Dx=(0.5*(1-x)*x,x(0)=0.05,t)s=1/(exp(log(19) - t/2) + 1)%再調(diào)用ode45計(jì)算數(shù)值解，并作圖比較解析解與數(shù)值解的區(qū)別：odefun=inline(0.5*(1-x)*x,t,x);t,x=ode45(odefun,0,20,0.05);plot(t,x,r*);hold on;ezplot(s,0 20);hold off;(2) 廠家問：要做多少時(shí)間廣告，可使市場購買率達(dá)到80%？t_min=min(find(x0.8);t(t_min)ans = 8.8543結(jié)果：大約8.8543個(gè)單位的時(shí)間后，可使市場購買率達(dá)到80%。6. (腫瘤生長) 腫瘤大小V生長的速率與V的a次方成正比，其中a為形狀參數(shù)，;而其比例系數(shù)K隨時(shí)間減小，減小速率又與當(dāng)時(shí)的K值成正比，比例系數(shù)為環(huán)境參數(shù)b。設(shè)某腫瘤參數(shù)a=1, b=0.1, K的初始值為2，V的初始值為1。問(1)此腫瘤生長不會(huì)超過多大？由已知條件建立微分方程：%先編寫上述函數(shù)的ex6_1fun.m文件function f=ex6_1fun(t,y)f(1)=y(2).*y(1);f(2)=-0.1*y(2);f=f(:);%畫出其圖像，并求最大的腫瘤大小Vt1,y1=ode45(ex6_1fun,0,100,1,2);plot(t1,y1(:,1),r*);max(y1(:,1)ans = 4.8567e+008%故腫瘤生長不會(huì)超過(2)過多長時(shí)間腫瘤大小翻一倍？%從圖像可以看出，大約在（0，1）內(nèi)，腫瘤大小翻一倍，以此求解t2,y2=ode45(ex6_1fun,0 1,1;2);t_min=min(find(y2(:,1)2);t2(t_min)ans =0.3750答：大約0.4個(gè)單位時(shí)間后腫瘤大小翻一倍。(3)何時(shí)腫瘤生長速率由遞增轉(zhuǎn)為遞減？dv=y1(:,2).*y1(:,1);Vn,tn=max(dv);t1(tn)ans = 29.5466答：大約30個(gè)單位時(shí)間后腫瘤生長速率由遞增轉(zhuǎn)為遞減。(4)若參數(shù)a=2/3呢？%重新編寫微分方程ex6_4fun，并依次計(jì)算上述三個(gè)問題function f=ex6_4fun(t,y)f(1)=y(2)*y(1).(2/3);f(2)=-0.1*y(2);f=f(:);%畫圖像，腫瘤生長不會(huì)超過450.7959%求多長時(shí)間腫瘤大小翻一倍，大約為0.4ans = 0.4000%何時(shí)腫瘤生長速率由遞增轉(zhuǎn)為遞減，大約為10ans = 9.5718選做題：1.(生態(tài)系統(tǒng)的振蕩現(xiàn)象)第一次世界大戰(zhàn)中，因?yàn)閼?zhàn)爭很少捕魚，按理戰(zhàn)后應(yīng)能捕到更多的魚才是。可是大戰(zhàn)后，在地中海卻捕不到鯊魚，因而漁民大惑不解。令x1為魚餌的數(shù)量，x2為鯊魚的數(shù)量，t為時(shí)間。微分方程為 (5.20)式中a1, a2, b1, b2都是正常數(shù)。第一式魚餌x1的增長速度大體上與x1成正比，即按a1x1比率增加, 而被鯊魚吃掉的部分按b1x1x2的比率減少；第二式中鯊魚的增長速度由于生存競爭的自然死亡或互相咬食按a2x2的比率減少，但又根據(jù)魚餌的量的變化按b2x1x2的比率增加。對(duì)a1=3, b1=2, a2=2.5, b2=1, x1(0)=x2(0)=1求解。畫出解曲線圖和相軌線圖，可以觀察到魚餌和鯊魚數(shù)量的周期振蕩現(xiàn)象。%首次編寫上述微分方程的ex_7fun.m函數(shù)function f=ex_7fun(t,x)a1=3;b1=2;a2=2.5;b2=1;f(1)=x(1)*(a1-b1*x(2);f(2)=-x(2)*(a2-b2*x(1);f=f(:);%畫出解曲線圖和相軌線圖t,x=ode45(ex_7fun,0,4,1;1);subplot(1,2,1);plot(t,x(:,1),r,t,x(:,2),k:);legend(x1,x2);title(解曲線);subplot(1,2,2);plot(x(:,1),x(:,2);title(相軌線);2.編寫四階Runge-Kutta法程序?qū)嶒?yàn)小結(jié)與收獲：完成第二份實(shí)驗(yàn)報(bào)告后，我