春數(shù)學實驗基礎 實驗報告(1)常微分方程.doc

年級、專業(yè) 姓名 學號 名單序號 實驗時間 2013年 3月 日 使用設備、軟件 PC, MATLAB 注： 實驗報告的最后一部分是實驗小結與收獲 實驗一 常微分方程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法只有一階精度，所以實際應用效率比較差，而ode45的效果比較好，很接近真實值。 (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. 求一通過原點的曲線，它在處的切線斜率等于若上限增為1.58，1.60會發(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);結論：隨著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. (溫度過程)夏天把開有空調的室內一支讀數(shù)為20的溫度計放到戶外，10分鐘后讀25.2, 再過10分鐘后讀數(shù)28.32。建立一個較合理的模型來推算戶外溫度。設：t時刻溫度計的讀數(shù)為T，戶外溫度為c，T的增速與室內外溫差（c-T）成正比，由此建立微分方程，其中k為比例系數(shù)%首先，計算解析解 y=dsolve(DT=k*(c-T),T(0)=20,t)y =c - (c - 20)/exp(k*t)%又已知，用非線性最小二乘擬合該函數(shù)，調用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. (廣告效應)某公司生產一種耐用消費品，市場占有率為5%時開始做廣告，一段時間的市場跟蹤調查后，該公司發(fā)現(xiàn)：單位時間內購買人口百分比的相對增長率與當時還沒有買的百分比成正比，且估得此比例系數(shù)為0.5。(1) 建立該問題的數(shù)學模型，并求其數(shù)值解與模擬結果作以比較；設：市場占有率為，市場占有率的增長速度為，則相對增長率為，由此建立微分方程為%首先，計算解析解 s=dsolve(Dx=(0.5*(1-x)*x,x(0)=0.05,t)s=1/(exp(log(19) - t/2) + 1)%再調用ode45計算數(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) 廠家問：要做多少時間廣告，可使市場購買率達到80%？t_min=min(find(x0.8);t(t_min)ans = 8.8543結果：大約8.8543個單位的時間后，可使市場購買率達到80%。6. (腫瘤生長) 腫瘤大小V生長的速率與V的a次方成正比，其中a為形狀參數(shù)，;而其比例系數(shù)K隨時間減小，減小速率又與當時的K值成正比，比例系數(shù)為環(huán)境參數(shù)b。設某腫瘤參數(shù)a=1, b=0.1, K的初始值為2，V的初始值為1。問(1)此腫瘤生長不會超過多大？由已知條件建立微分方程：%先編寫上述函數(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%故腫瘤生長不會超過(2)過多長時間腫瘤大小翻一倍？%從圖像可以看出，大約在（0，1）內，腫瘤大小翻一倍，以此求解t2,y2=ode45(ex6_1fun,0 1,1;2);t_min=min(find(y2(:,1)2);t2(t_min)ans =0.3750答：大約0.4個單位時間后腫瘤大小翻一倍。(3)何時腫瘤生長速率由遞增轉為遞減？dv=y1(:,2).*y1(:,1);Vn,tn=max(dv);t1(tn)ans = 29.5466答：大約30個單位時間后腫瘤生長速率由遞增轉為遞減。(4)若參數(shù)a=2/3呢？%重新編寫微分方程ex6_4fun，并依次計算上述三個問題function f=ex6_4fun(t,y)f(1)=y(2)*y(1).(2/3);f(2)=-0.1*y(2);f=f(:);%畫圖像，腫瘤生長不會超過450.7959%求多長時間腫瘤大小翻一倍，大約為0.4ans = 0.4000%何時腫瘤生長速率由遞增轉為遞減，大約為10ans = 9.5718選做題：1.(生態(tài)系統(tǒng)的振蕩現(xiàn)象)第一次世界大戰(zhàn)中，因為戰(zhàn)爭很少捕魚，按理戰(zhàn)后應能捕到更多的魚才是。可是大戰(zhàn)后，在地中海卻捕不到鯊魚，因而漁民大惑不解。令x1為魚餌的數(shù)量，x2為鯊魚的數(shù)量，t為時間。微分方程為 (5.20)式中a1, a2, b1, b2都是正常數(shù)。第一式魚餌x1的增長速度大體上與x1成正比，即按a1x1比率增加, 而被鯊魚吃掉的部分按b1x1x2的比率減少；第二式中鯊魚的增長速度由于生存競爭的自然死亡或互相咬食按a2x2的比率減少，但又根據(jù)魚餌的量的變化按b2x1x2的比率增加。對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法程序實驗小結與收獲：完成第二份實驗報告后，我