實驗一常微分方程解析

1、年級、專業(yè) 2012級數(shù)學與應用數(shù)學姓名 張旭 學號 12160011063 名單序號16實驗一常微分方程ode45的精度較高,實驗時間2014年3月19日使用設備、軟件PC, MATLAB注：實驗報告的最后一部分是實驗小結與收獲1.分別用Euler法和ode45解下列常微分方程并與解析解比較:(1) 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(od efun ,xspa n,y 0,h) x=xspa n(1):h:xspa n( 2

2、);y(1)=y0;for i=1:le ngth(x)-1y(k+1)=y(i)+h*feval(odefun ,x(i),y(i);endx=x;y=y;endode45: odefu n=in li ne(x+y,x,y);xspa n=0,3;y0=1;h=0.1;x1,y1=euler(od efun,xspan,y0,h); x2,y2=od e45(od efun,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

3、,od e45,dsolve);ode45求得的結果與用解析法求得的結果更接近，故Euler法求得的結果精度較低。3 / 92014春數(shù)學實驗實驗一常微分方程注:實驗報告的最后一部分是實驗小結與收獲 y 0.01(y)2 2y =sin(t), y(0) = 0,y (0) = 1,0 : t 5令丫勺=y, y2 = y：則原方程等價于方程組：丿2y2 二 一2%0.01y2 sin(t)%(0) =0,y2(0) =1，0 :：: t 5，不能解析，只能用數(shù)值法求解。Euler: function t,y=eul er2(od efun 1,odefun2,tspan,y0,h) t=ts

4、pa n(1):h:tspa n( 2);y(1,1)=y0(1);y(2,1)=y0(2);for i=1:le ngth(t)-1k1=od efun 1(t(i),y(1,i),y(2,i);k2=od efun2(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;endt=t;y=y;endode45 :odefu n1=i nli ne(0*t1+0*y1+y2);odefu n2=i nlin e(-2*y1+0.01*y2A2+si n( t1); t1,y1=euler2(od efun 1,od efun

5、2,0,5,0,1,0.1); t2,y2=od e45(eu,0,5,0,1);plot(t1,y1(:,1),o,t2,y2(:,1),Li neWidth,2); xlabel(t 軸);ylabel(y 軸);legend(euler,od e45);ode45 中 eu: function dy=eu(t,y)實驗時間2014年3月19日使用設備、軟件PC, MATLAB注：實驗報告的最后一部分是實驗小結與收獲dy=zeros(2,1);dy(1)=y(2);dy(2)=-2*y(1)+0.01*y(2)A2+s in (t);ode45求得的結果精度較高，euler法求得的結果在準

6、確值上下波動。2. 一通過原點的曲線，它在(x,y)處的切線斜率等于 2xy2,： x ： 1.57.若x上限增為1.58，1.60 會發(fā)生什么？等價于求解y丄2x y2，y(0) =0,且0 ： x ： 1.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/2)*airy(0,-2A (1/3)*x) + airy(2, -2 A(1/3)*x)ode45 :odefu n=inlin e(2*x+yA2);subpl

7、ot(1,3,1);od e45(od efun,0,1.57,0);title(0x1.57);subplot(1,3,2);od e45(od efun,0,1.58,0);title(0x1.58);subplot(1,3,3);od e45(od efun,0,1.60,0);title(0x1.60);Dix-tl .57口 a1.5之后的斜率增長速度很快，若x上限增為1.58 , 1.60，則相應的y將會出現(xiàn)更大的增長。年級、專業(yè) 2012級數(shù)學與應用數(shù)學姓名 張旭 學號 12160011063 名單序號16實驗時間2014年3月19 日使用設備、軟件PC, MATLAB實驗一常微

8、分方程注:實驗報告的最后一部分是實驗小結與收獲3. 求解剛性方程組:0 ： x ： 50.1000.25 y 999.75y20.5, %(0)=1,y2 =999.75力-1000.25 y20.5,y1(0 1,fun ctio n dy=fu n( 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=od e15s(fun,0,50,1,-1); plot(t,y(:,1),o,t,y(:,2),k-,Li neWidth,2); le

9、ge nd(y1,y2);5 / 92014春數(shù)學實驗實驗一常微分方程4. (溫度過程)夏天把開有空調的室內(nèi)一支讀數(shù)為20 C的溫度計放到戶外，10分鐘后讀25.2 C,再過10分鐘后讀數(shù)28.32 C。建立一個較合理的模型來推算戶外溫度。由題意可知由于隨著時間的增加，溫度增長越來越慢，戶外溫度與溫度計示數(shù)之差也越來越小，且溫差為零時溫度的增長率也為零，故可以認為溫度的增長率與溫差成正比，設戶外溫度為m，溫度計的示數(shù)為 y，比例系數(shù)為k，則可建立模型y =k(m -y), y(0) =20解析法:y=dsolve(Dy=k*(c-y),y(0)=20,t)注:實驗報告的最后一部分是實驗小結與收

10、獲y = m - (m - 20)/exp(k*t)由 y (10)=25.2,y(20)=28.32建立方程組m20“ cm 10k25-2em 20m -20k28-32e消去k，得:(m-20)(m-28.32)= ( m-25.2)( m-25.2 )解得：m=33所以戶外溫度約為 33C5. (廣告效應)某公司生產(chǎn)一種耐用消費品，市場占有率為5%時開始做廣告，一段時間的市場跟蹤調查后， 該公司發(fā)現(xiàn)：單位時間內(nèi)購買人口百分比的相對增長率與當時還沒有 買的百分比成正比，且估得此比例系數(shù)為0.5。(1)建立該問題的數(shù)學模型，并求其數(shù)值解與模擬結果作以比較;設t時刻該消費品的市場占有率為y,

11、建立方程：/-0.5* (1 y), y(0) =5%解析解：y=dsolve(Dy=0.5-0.5*y,y(0)=0.05)y = 1 - (19*exp(-t/2)/20數(shù)值解：odefu n=inlin e(0.5-0.5*y,t,y);t1,y1=od e45(od efun,0,10,0.05); t2=0:0.1:10;y2=1-(19*exp(-t2/2)/20;plot(t1,y1,o,t2,y2,k);lege nd(od e,dsolve);年級、專業(yè) 2012級數(shù)學與應用數(shù)學姓名 張旭 學號 12160011063 名單序號16實驗一常微分方程實驗時間2014年3月19日

12、使用設備、軟件PC, MATLAB注：實驗報告的最后一部分是實驗小結與收獲(2) 廠家問：要做多少時間廣告，可使市場購買率達到 80% ？4 由解析解可列出方程1 - (19*exp(-t/2)/20=0.8 ，所以t - -2ln -19解得 t=3.11636. (腫瘤生長)腫瘤大小V生長的速率與 V的a次方成正比，其中a為形狀參數(shù)，0蟲乞1; 而其比例系數(shù)K隨時間減小，減小速率又與當時的 K值成正比，比例系數(shù)為環(huán)境參數(shù) b。 設某腫瘤參數(shù)a=1, b=0.1, K的初始值為2，V的初始值為1。問(1)此腫瘤生長不會超過多大？k =k = -0.1k由已知列出方程組丿，代入具體數(shù)值，得丿，

13、v=kvav = kvk(0) =2,v(0) -1，(k,v都是關于時間t的函數(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,L in eWidth,2); xlabel(t 軸);ylabel(v 軸);因腫瘤不斷長大，故t趨于無窮時，該腫瘤達到最大，此時極限為exp(20)=4.8517*10A8，故此腫瘤生長不會超過4.8517*10人89 / 92

14、014春數(shù)學實驗實驗一常微分方程實驗時間2014年3月19日使用設備、軟件PC, MATLAB注：實驗報告的最后一部分是實驗小結與收獲(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(-20*exp(-t1/10); plot(t1,v1,Li neWidth,2);

15、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)得到 t = 29.96若參數(shù)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(-t/10)2*exp(-t/10)v = -(20*exp(-t/10) - 23)A

16、3/27年級、專業(yè) 2012級數(shù)學與應用數(shù)學姓名 張旭 學號 12160011063 名單序號16實驗時間2014年3月19 日使用設備、軟件PC, MATLAB實驗一常微分方程注:實驗報告的最后一部分是實驗小結與收獲(37/2 + (3人(1/2)*3許)/2 - 20*exp(-t/10)A3/27-(20*exp(-t/10) + (3A(1/2)*3*i)/2 - 37/2)人3/27取實解 k = 2*exp(-t/10)v=-(20*exp(-t/10) - 23)人3/27并畫出v-t圖像：t=0:0.1:100;v=-(20*exp(-t/10) - 23).A3/27;plo

17、t(t,v,L in eWidth,2);xlabel(t 軸);ylabel(v 軸);顯然，當t趨于無窮時，該腫瘤達到最大，此時極限為-(-23)人3/27=450.6296，故此腫瘤 生長不會超過450.62962、令-(20*exp(-t/10) - 23)A3/27=2 ，解得 t=0.3977,3、由已求得的結果可得v 與t的關系為v =2*exp(-t/10)*(20*exp(-t/10) - 23)人2/9t1=0:0.1:100;v1= 2*exp(-t1/10).*(20*exp(-t1/10) - 23).人2/9; plot(t1,v1,Li neWidth,2);xl

18、abel(t 軸);ylabel(v 軸);顯然，最大值處對應的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)中，因為戰(zhàn)爭很少捕魚，按理戰(zhàn)后應能捕到更11 / 92014春數(shù)學實驗實驗一常微分方程注:實驗報告的最后一部分是實驗小結與收獲多的魚才是。可是大戰(zhàn)后，在地中海卻捕不到鯊魚，因而漁民大惑不解。令X1為魚餌的數(shù)量，X2為鯊魚的數(shù)量，t為時間。微分方程為dt式中a1, a 2, b 1, b 2都是正常數(shù)。第一式魚餌X1的增長速度大體上與 X1成正比，即按a1X1比率增加，而被鯊魚吃掉的部分按 b1X1X2的比率減少；第二式中鯊魚的增長速度由 于生存競爭的自然死亡或互相咬食按a2X2的比率減少，但又根據(jù)魚餌的量的變化按b2X1 X2 的比率增加。對 a 1=3, b 1=2, a 2=2.5, b2=1, x 1（0）=x