試驗一常微分方程

1、實用文檔1.分別用Euler法和ode45解下列常微分方程并與解析解比較:y： =x y, y(0) =1,0 :二 x :二 3解析法：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)-1y(k+1)=y(i)+h*feval(odefun,x(i),y(i);endx=x;JOk!by=y;endode45:odefun=inline(x+y,x,y);xspan=0,3

2、;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求得的結果與用解析法求得的結果更接近，故ode45的精度較高,Euler法求得的結果精度較低。y -0.01(y)2 2y =sin(t), y(0) = 0, y (0) =1,0 二 t :二 5實用文檔令yi = y

3、, y2 = y,則原方程等價于方程組:yi =y2 2y = -2yi 0.01y2sin(t)yi（0） =0,y2（0） =1, 0 ct 5 ,不能解析，只能用數值法求解。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)-1k1=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+

4、1)=y(2,i)+h*d2;endt=t;y=y;end ode45:odefun1=inline(0*t1+0*y1+y2);odefun2=inline(-2*y1+0.01*y2A2+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)

5、;dy(1)=y(2);dy(2)=-2*y(1)+0.01*y(2)A2+sin(t);ode45求得的結果精度較高，euler法求得的結果在準確值上下波動。22. 一通過原點的曲線，它在 （x, y）處的切線斜率等于 2x + y ,0 x 1.57.若x上限增為實用文檔1.58 , 1.60會發(fā)生什么？2等價于求解y = 2x + y , y(0) =0，且0 x1.57的初值問題。解析法：y=dsolve(Dy=2*x+yA2,y(0)=0,x)y=(2A(1/3)*airy(3,-2A(1/3)*x)+2A(1/3)*3A(1/2)*airy(1,-2A(1/3)*x)/(3A(1/

6、2)*airy(0, -2A(1/3)*x) + airy(2, -2A(1/3)*x)ode45:odefun=inline(2*x+yA2);subplot(1,3,1);ode45(odefun,0,1.57,0);title(0 x1.57);subplot(1,3,2);ode45(odefun,0,1.58,0);title(0 x1.58);subplot(1,3,3);ode45(odefun,0,1.60,0);title(0 x1.60);.57Kid巾ax1.5之后的斜率增長速度很快，若 x上限增為1.58 ,1.60 ,則相應的y將會出現更大的增長。3.求解剛性方程組：

7、y； = 7000.25yl 999.75y2 05 yl(0) = 1,y2 =999.75yl -1000.25y2 0.5,y1(0) =1,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;ode45:t,y=ode15s(fun,0,50,1,-1);plot(t,y(:,1),o,t,y(:,2),k-,LineWidth,2);legend(y1,y2);4.（溫度過程）夏天把開有空調的室內一支讀數為20c的溫度計

8、放到戶外，10分鐘后讀25.2 C ,再過10分鐘后讀數28.32 C 。建立一個較合理的模型來推算戶外溫度。由題意可知由于隨著時間的增加，溫度增長越來越慢，戶外溫度與溫度計示數之差也越來越小，且溫差為零時溫度的增長率也為零，故可以認為溫度的增長率與溫差成正比，設戶外溫度為m,溫度計的示數為 y,比例系數為k,則可建立模型y =k(m-y), y(0) =20解析法：y=dsolve(Dy=k*(c-y),y(0)=20,t)y = m - (m - 20)/exp(k*t)由 y (10) =25.2,y(20)=28.32建立方程組m-20 mem-20 me= 25.2,消去 k,得：

9、(m-20)= 28.32(m-28.32 ) = (m-25.2) (m-25.2 )解得：m=33 所以戶外溫度約為 33c 。5.（廣告效應）某公司生產一種耐用消費品，市場占有率為5%寸開始做廣告，一段時間的實用文檔市場跟蹤調查后，該公司發(fā)現：單位時間內購買人口百分比的相對增長率與當時還沒有買 的百分比成正比，且估得此比例系數為0.5。(1)建立該問題的數學模型，并求其數值解與模擬結果作以比較；設t時刻該消費品的市場占有率為y,建立方程:y=0.5* (1 -y), y(0) =5%解析解：y=dsolve(Dy=0.5-0.5*y,y(0)=0.05)y = 1 - (19*exp(-

10、t/2)/20數值解：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)廠家問：要做多少時間廣告，可使市場購買率達到 80%4由解析解可列出萬程1 - (19*exp(-t/2)/20=0.8,所以t = 2ln 19解得 t=3.11636.(腫瘤生長)腫瘤大小V生長的速率與 V的a次方成正比，其中a為形狀參數，0芻E1; 而其比例系數K隨時間減小，減小速率又與當時的K

11、值成正比，比例系數為環(huán)境參數 bo設某腫瘤參數a=1, b=0.1, K 的初始值為2, V的初始值為1。問(1)此腫瘤生長不會超過多大？,八，、 k = -bkn k = 0.1k由已知列出方程組 ， a ，代入具體數值，得，v=kvaV = kv實用文檔k(0) =2,v(0) =1, (k,v都是關于時間t的函數)解析法：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

12、,v,LineWidth,2);xlabel(t 軸)；ylabel(v 軸)；因腫瘤不斷長大，故t趨于無窮時，該腫瘤達到最大，此時極限為exp(20)=4.8517*10A8 , 故此腫瘤生長不會超過 4.8517*10A8 。(2)過多長時間腫瘤大小翻一倍？令 exp(20)*exp(-20*exp(-t/10)=2,解得 t=-10*ln(1-1/20*ln2)=0.3527,(3)何時腫瘤生長速率由遞增轉為遞減？由已求得的結果可得V與t的關系為v=2*exp(20-t/10)*exp(-20*exp(-t/10)t1=0:0.1:100;v1=2*exp(20-t1/10).*exp(

13、-20*exp(-t1/10);plot(t1,v1,LineWidth,2);xlabel(t 軸)；ylabel(v 軸)；顯然，最大值處對應的t即為所求：t2=0:0.01:100;v2=2*exp(20-t1/10).*exp(-20*exp(-t1/10);實用文檔m,n=max(v2);t=t2(n)得至ij t = 29.96(4)若參數a=2/3呢？k = -0.1k1、建立方程組2, k(0) =2,v(0) =1v = kv3解析法：k,v=dsolve(Dk=-0.1*k,Dv=k*vA(2/3),k(0)=2,v(0)=1,t)k = 2*exp(-t/10)2*exp

14、(-t/10)2*exp(-t/10)v = -(20*exp(-t/10) - 23)A3/27(37/2 + (3A(1/2)*3*i)/2 - 20*exp(-t/10)F3/27-(20*exp(-t/10) + (3A(1/2)*3*i)/2 - 37/243/27取實解 k = 2*exp(-t/10) v=-(20*exp(-t/10) - 2343/27并畫出v-t圖像：t=0:0.1:100;v=-(20*exp(-t/10) - 23).A3/27;plot(t,v,LineWidth,2);xlabel(t軸)；ylabel(v軸)；顯然，當t趨于無窮時，該腫瘤達到最大，

15、此時極限為 -(-23)A3/27=450.6296,故此腫瘤生長不會超過 450.6296 。2、令-(20*exp(-t/10) - 23)A3/27=2,解得 t=0.3977,3、由已求得的結果可得v與t的關系為實用文檔v =2*exp(-t/10)*(20*exp(-t/10) - 23)A2/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軸)；顯然，最大值處對應的 t即為所求，： t2=0:0.01:100;v2=2*exp

16、(-t2/10).*(20*exp(-t2/10) - 23).人2/9;m,n=max(v2);t=t2(n) t = 9.5900選做題:1.(生態(tài)系統的振蕩現象)第一次世界大戰(zhàn)中，因為戰(zhàn)爭很少捕魚，按理戰(zhàn)后應能捕到更多 的魚才是?？墒谴髴?zhàn)后，在地中海卻捕不到鯊魚，因而漁民大惑不解。令x1為魚餌的數量，x2為鯊魚的數量，t為時間。微分方程為式中a1, a 2, b 1, b 2都是正常數。第一式魚餌x1的增長速度大體上與 玄成正比，即按a比率增加，而被鯊魚吃掉的部分按bx1x2的比率減少；第二式中鯊魚的增長速度由于生存競爭的自然死亡或互相咬食按a2x2的比率減少，但又根據魚餌的量的變化按b2x1x2的比率增加。對a1=3, b 1=2, a 2=2.5, b 2=1, x 1(0)=x 2(0)=1求解。畫出解曲線圖和相軌線圖，可 以觀察到魚餌和鯊魚數量的周期振蕩現象。x1 =3x1