數(shù)學(xué)實(shí)驗(yàn)第二次作業(yè)——常微分方程數(shù)值求解

實(shí)驗(yàn)4常微分方程數(shù)值解實(shí)驗(yàn)?zāi)康模?. 練習(xí)數(shù)值積分的計(jì)算；2. 掌握用MATLAB軟件求微分方程初值問(wèn)題數(shù)值解的方法；3. 通過(guò)實(shí)例學(xué)習(xí)用微分方程模型解決簡(jiǎn)化的實(shí)際問(wèn)題；4. 了解歐拉方法和龍格庫(kù)塔方法的基本思想和計(jì)算公式，及穩(wěn)定性等概念。實(shí)驗(yàn)內(nèi)容：3.小型火箭初始質(zhì)量為1400kg，其中包括1080kg燃料，火箭豎直向上發(fā)射是燃料燃燒率為18kg/s，由此產(chǎn)生32000N的推力，火箭引擎在燃料用盡時(shí)關(guān)閉。設(shè)火箭上升是空氣阻力正比于速度的平方，比例系數(shù)為0.4kg/m，求引擎關(guān)閉瞬間火箭的高度，速度，加速度，及火箭到達(dá)最高點(diǎn)是的高度，速度和加速度，并畫出高度，速度，加速度隨時(shí)間變化的圖形。解答如下：這是一個(gè)典型的牛頓第二定律問(wèn)題，分析火箭受力情況；先規(guī)定向上受力為正數(shù)建立數(shù)學(xué)模型：A燃料未燃盡前，在任意時(shí)刻（t60s）火箭受到向上的-F=32000N，向下的重力G=mg，g=9.8，向下的阻力f=kv2, k=0.4, v表示此時(shí)火箭速度；此時(shí)火箭收到的合力為F1=（F-mg-f）；火箭的初始質(zhì)量為1400kg，燃料燃燒率為-18kg/s;此刻火箭質(zhì)量為m=1400-18*t根據(jù)牛頓第二定律知，加速度a=F1/m=(F-mg-f)/(m-r*t)= (32000-(0.4.*v.2)-9.8.*(1400-18.*t)由此可利用龍格-庫(kù)塔方法來(lái)實(shí)現(xiàn)，程序?qū)崿F(xiàn)如下Function dx=rockett,x %建立名為rocket的方程m=1400;k=0.4;r=-18;g=9.8; %給出題目提供的常數(shù)值dx=x(2);(32000-(k*x(2)2)-g*(m+r*t)/(m+r*t);%以向量的形式建立方程a=(32000-(k*x(2)2)-g*(m+r*t)/(m+r*t); %給出a的表達(dá)式End;ts=0:60; %根據(jù)題目給定燃燒率計(jì)算出燃料燃盡的時(shí)間，確定終點(diǎn)x0=0,0; %輸入x的初始值t,x=ode15s(rocket,ts,x0); %調(diào)用ode15s計(jì)算t,x;h=x(:,1);v=x(:,2);plot(t,x(:,1),grid; %繪出火箭高度與時(shí)間的關(guān)系曲線title(h-t);xlabel(t/s);ylabel(h/m),pause;plot(t,x(:,2),grid ; %繪出火箭速度與時(shí)間的曲線關(guān)系title(v-t);xlabel(t/s);ylabel(v/m/s),pause;a=(32000-(0.4.*v.2)-9.8.*(1400-18.*t)/(1400-18.*t);plot(t,a),grid; %繪出火箭加速度與時(shí)間的曲線關(guān)系title(a-t);xlabel(t/s),ylabel(a/m2/s),pause 火箭高度隨時(shí)間變化的曲線火箭速度隨時(shí)間變化的曲線火箭加速度隨時(shí)間變化的曲線數(shù)據(jù)過(guò)多，故截取部分如下第一列為時(shí)間，第二列為火箭高度，第三列為火箭速度由此可以，在t=60s時(shí)，即火箭燃料燃盡瞬間，引擎關(guān)閉瞬間，火箭將到達(dá)12912m的高度，速度為267,29m，加速度a=0.9m/s2B燃料燃盡之后，與A 類似，分析受力如下火箭受到向上的F=0向下的重力G=mg，g=9.8，向下的阻力f=kv2, k=0.4, v表示此時(shí)火箭速度；此時(shí)火箭收到的合力為F2=（-mg-f）；火箭的初始質(zhì)量為320kg，恒定根據(jù)牛頓第二定律,加速度a=F2/m=-g-0.4v2/320;程序?qū)崿F(xiàn)如下function dx = rocket2( t,x ) %建立以rocket2為名的函數(shù)dx=x(2);-9.8-0.4.*x(2).2/320; %以向量的形式建立方程ts=60:120; %給出初始時(shí)刻，估計(jì)終點(diǎn)時(shí)刻x0=12190,267.26; %給出x初始值t,x=ode15s(rocket2,ts,x0); %調(diào)用ode15s計(jì)算t,xplot(t,x(:,1),grid; %繪出火箭高度隨時(shí)間變化的曲線title(h-t);xlabel(t/s),ylabel(h/m),pause;plot(t,x(:,2),grid; %繪出火箭速度隨時(shí)間的變化曲線title(v-t);xlabel(t/s),ylabel(v/m/s),pause;v=x(:,2);a=-9.8-0.4*v.2/320; %給出加速度的具體表達(dá)式plot(t,a),grid; %繪出火箭加速度隨時(shí)間變化的曲線title(a-t);xlabel(t/s),ylabel(a/m2/s),pause得到的曲線圖形如下火箭高度隨時(shí)間的變化曲線從圖中可以大致看出，最高點(diǎn)在13km左右，火箭速度隨時(shí)間的變化曲線加速度隨時(shí)間變化曲線如下數(shù)據(jù)表格大致如下從圖表中可以看出，在71s左右速度到達(dá)0，即此時(shí)到達(dá)最高處，高度為13117m加速度為-9.8m/m/s2;本題總結(jié)：這道題是典型的物理牛頓力學(xué)的題目，通過(guò)受力的正確分析，可以知道，以h,v為向量建立微分方程即可求解，h的微分是速度v，速度v的微分是加速度a解題過(guò)程中存在的難點(diǎn)是：取值步長(zhǎng)不太容易確定，而且是哪種算法不確定， 先用ode15s速度較快，ode23s速度差不太多，其他兩種速度較慢，等待時(shí)間較長(zhǎng)5一只小船渡過(guò)寬為d 的河流，目標(biāo)是起點(diǎn)A 正對(duì)著的另一岸B 點(diǎn)。已知河水流速為v1 與船在靜水中的速度v2之比為k。（1） 建立描述小船航線的數(shù)學(xué)模型，求其解析解；（2） 設(shè)d=100m，v1=1m/s，v2=2m/s，用數(shù)值解法求渡河所需的時(shí)間、任意時(shí)刻小船的位置及航行曲線，作圖，并與解析解比較。（3） 若流速v1=0，0.5，1.5，2（m/s）,情況又如何建立數(shù)學(xué)模型：在任意時(shí)刻t，小船位于（x，y），此時(shí)速度為v，根據(jù)物理中路程與速度的關(guān)系，知路程的微分為速度v，由此中，小船在x,y方向上的速度分別為：X：dxdt=v1-v2*xx2+y2Y:dydt=-v2*yx2+y2初始條件為dxdtt0=v1 *（1）dydt t0=v2 *（2）現(xiàn)求其解析解,利用微積分知識(shí)*（1）*（2）得dxdy=v1*x2+y2-v2y+xy;*(3)將*（3）右端帶根號(hào)部分分子分母均y得dxdy=v1*(x/y)2+1-v2+xy;令x/y=p，得到dx/dy=dp/dy*y+p;dx/dy=p;分離變量有y(-k)=c(1+p2+p);代入初值可確定當(dāng)t=0時(shí)，y=-d,x=0,p=x/y=0,C=(-d)(-k)y(-k)/d(-k)=1+x2y2 +x/y;根據(jù)隱函數(shù)與顯函數(shù)的關(guān)系可得到解析解：x=-d2(y-d)(-k)-(yd)(1+k)已知d=100m,v1=1m/s,v2=2m/s;先利用龍格庫(kù)塔方法求解渡河時(shí)間，及任意時(shí)刻小船的位置（x，y），及航行曲線，與解析解比較此時(shí)，k=0.5，d=100用MATLB編寫源程序如下：function dx = boat(t,x) %建立名為boat的函數(shù) d=-100;v1=1;v2=2; %給定常數(shù)值s=(x(1).2+x(2).2).0.5; dx=v1-v2.*x(1)./s;-v2.*x(2)./s; %用向量形式建立方程ts=0:70; %大致估算，確定終點(diǎn)值，給定步長(zhǎng)為”1”x0=0,-100; %給出x初始值t,x=ode23s(boat,ts,x0); %調(diào)用ode23s計(jì)算t,x;plot(t,x),grid, %繪出x(t),y(t)的函數(shù)曲線（圖21）gtext(x(t);gtext(y(t);pause;plot(t,x(:,1), %繪出x隨時(shí)間變化的曲線（圖22）grid;xlabel(t/s);ylabel(x/m);pause;plot(t,x(:,2),grid, %繪出y隨時(shí)間變化的曲線（圖23）title(y-t);xlabel(t/s),ylabel(y/m);圖21 圖22圖23得到數(shù)據(jù)如下從表格中數(shù)據(jù)可知，在大約67s時(shí)y=0即船到達(dá)對(duì)岸目的地，為比較，先進(jìn)行解析解的求解設(shè)計(jì)程序如下：function x=f(y)k=0.5;x=-0.5.*(-0.01).k.*y.(k+1)+0.5.*(-0.01).(-k).*y.(-k+1);y=0:-0.1:-100;for i=0:1:1000;x(:,i+1)=xy(-i/10);endplot(x,y);grid;gtext(x);gtext(y);由此可以看出，由數(shù)值解和解析解得到的x-y曲線相差不多，所以可以認(rèn)為解析解正確改變水流速度v1，只要在原有程序基礎(chǔ)上重新復(fù)制給v1=0,0.5,1.5,2,同時(shí)適當(dāng)改變終點(diǎn)值即可現(xiàn)實(shí)現(xiàn)程序如下A:v1=0，d=-100;v1=0;v2=2;s=(x(1).2+x(2).2).0.5;dx=v1-v2.*x(1)./s;-v2.*x(2)./sts=0:60;x0=0,-100;t,x=ode15s(boat21,ts,x0);t,x;plot(x(:,1),x(:,2),grid,title(y-x);pause;plot(t,x(:,1),grid;xlabel(t/s);ylabel(x/m);plot(t,x(:,1),grid;title(x-t);xlabel(t/s),ylabel(x/m),pause;plot(t,x(:,2),grid,title(y-t);xlabel(t/s),ylabel(y/m);圖形如下從此圖形中我們可以看到，船并未偏離x=0的點(diǎn)，我們也可以從直觀想象中的得到，當(dāng)水速v1=0時(shí)，只要出發(fā)時(shí)，船頭對(duì)準(zhǔn)目標(biāo)點(diǎn)，船將一直朝著直線向目的地行進(jìn)從表格中數(shù)據(jù)我們也可以很清楚地看到路程與時(shí)間是成明顯的線性關(guān)系的，這是與我們水速為0的必然結(jié)果，由此也可以驗(yàn)證我們模型基本正確現(xiàn)改變水速B:v1=0.5程序?qū)崿F(xiàn)如下d=-100;v1=0.5;v2=2;s=(x(1).2+x(2).2).0.5;dx=v1-v2.*x(1)./s;-v2.*x(2)./sts=0:60;x0=0,-100;t,x=ode15s(boat22,ts,x0);t,x;plot(t,x),grid,gtext(x(t);gtext(y(t);pause;plot(t,x(:,1),grid;xlabel(t/s);ylabel(x/m);plot(t,x(:,1),grid;title(x-t);xlabel(t/s),ylabel(x/m),pause;plot(t,x(:,2),grid,實(shí)現(xiàn)程序過(guò)程中發(fā)現(xiàn)，終點(diǎn)值設(shè)為60，而在53s之后不再有出現(xiàn)，所以可以認(rèn)為在54s之前就已經(jīng)達(dá)到對(duì)岸現(xiàn)在重新改變水速C:v1=1.5d=-100;v1=1.5;v2=2;s=(x(1).2+x(2).2).0.5;dx=v1-v2.*x(1)./s;-v2.*x(2)./s ts=0:120;x0=0,-100;t,x=ode15s(boat23,ts,x0);t,x;plot(x(:,1),x(:,2),grid,title(y-x);pause;plot(t,x(:,1),grid;xlabel(t/s);ylabel(x/m);plot(t,x(:,1),grid;title(x-t);xlabel(t/s),ylabel(x/m),pause;plot(t,x(:,2),grid,title(y-t);xlabel(t/s),ylabel(y與v1=0.5s類似，在114s之后不再給出數(shù)據(jù)，而我們?cè)O(shè)定的終點(diǎn)值是120，所以可以大致在114s時(shí)到達(dá)B點(diǎn)現(xiàn)改變水速D:v1=2程序?qū)崿F(xiàn)如下d=-100;v1=2;v2=2;s=(x(1).2+x(2).2).0.5;dx=v1-v2.*x(1)./s;-v2.*x(2)./sts=0:300;x0=0,-100;t,x=ode15s(boat24,ts,x0);t,x;plot(x(:,1),x(:,2),grid,title(y-x);pause;plot(t,x(:,1),grid;xlabel(t/s);ylabel(x/m);plot(t,x(:,1),grid;title(x-t);xlabel(t/s),ylabel(x/m),pause;plot(t,x(:,2),grid,title(y-t);xlabel(t/s),ylabel(y/m)曲線如下我們可以從此圖中看出，當(dāng)y=0時(shí)，x=50，也就是說(shuō)，船根本沒本法到達(dá)正對(duì)著的B點(diǎn)，而只能到達(dá)對(duì)岸，我們可以直觀地理解假設(shè)我們的船在開始時(shí)與水速反向，這時(shí)船與水的合速度是0，故船無(wú)法前進(jìn)，而根據(jù)方程知道，在船尚未到達(dá)對(duì)岸之前，只有當(dāng)船的x軸向速度在模式刻超過(guò)水速，才有可能克服因水速而在x軸偏離B點(diǎn)的距離，而重新返回，到達(dá)B點(diǎn)，而在此給定的速度中，我們可以看到，船速的分量不可能超過(guò)水速，因此不可能到達(dá)B本題總結(jié)：這道題跟課本例題緝私船有較大的相似處，大體是模仿例題而做，所以在畫圖上顯得有點(diǎn)繁瑣，如果采用subplot作圖，將會(huì)節(jié)約一些篇幅，使得作業(yè)版面看起來(lái)更整潔一點(diǎn)3.兩種群的競(jìng)爭(zhēng)模型如下：其中x(t)，y(t)分別為甲乙兩種群的數(shù)量；r1，r2為他們的固有增漲率；n1，n2為他們的最大容量。s1的含義是，對(duì)于供養(yǎng)甲的資源而言，單位數(shù)量乙（相對(duì)n2）的消耗為單位數(shù)量甲的s1倍，對(duì)s2可作相應(yīng)的解釋。 該模型無(wú)解析解，試用數(shù)值解法研究一下問(wèn)題： （1） 設(shè)r1=r2=1，n1=n2=100，s1=0.5，s2=2，初值x0=y0=10，計(jì)算x(t)，y(t)，畫出他們的圖形及相圖(x,y)，說(shuō)明時(shí)間t充分大以后x(t)，y(t)的變化趨勢(shì)。 （2） 改變r(jià)1，r2，n1，n2，x0，y0，但s1，s2不變（或保持s11），計(jì)算并分析所得結(jié)果；若s1=1.5，s2=0.7,再分析結(jié)果，由此你能得到什么結(jié)論，請(qǐng)用各參數(shù)生態(tài)學(xué)上的含義做出解釋。 （3） 試驗(yàn)當(dāng)s1=0.8，s2=0.7時(shí)會(huì)有什么結(jié)果；當(dāng)s1=1.5，s2=1.7時(shí)又會(huì)有什么結(jié)果。你能解釋這些結(jié)果嗎？（1）程序如下：函數(shù)：function dx = species1(t,x ); %建立以species1為名的函數(shù)r1=1;r2=1;n1=100;n2=100;s1=0.5;s2=2; %給定常數(shù)值dx=r1*x(1).*(1-x(1)/n1-s1*x(2)/n2);r2*x(2).*(1-s2*x(1)/n1-x(2)/n2);%以向量形式建立方程end代碼：ts=1:0.01:100; %給定終點(diǎn)值，確定步長(zhǎng)x0=10,10; %給定x初始值t,x=ode45(species1,ts,x0); %調(diào)用ode45計(jì)算A=t,x;plot(t,x(:,1),t,x(:,2),grid; %繪出種群x,y隨時(shí)間變化的曲線 gtext(x(t);gtext(y(t);xlabel(t);ylabel(x,y);title(x,y-t圖像);pause;plot(x(:,1),x(:,2),grid; %繪出x y種群數(shù)量隨時(shí)間變化的曲線xlabel(x);ylabel(y);title(y-x圖像);種群x,y隨t的變化圖形如下：讀圖可知，當(dāng)t10后，兩種群數(shù)量基本不變，x=100,y=0滅絕。從生態(tài)學(xué)角度分析，以上反映了優(yōu)勝劣汰的自然規(guī)律。很顯然，種群x更具有競(jìng)爭(zhēng)優(yōu)勢(shì)，隨著自然的長(zhǎng)期演變，x得以繁衍，而y最終被滅絕。種群y與x的數(shù)量關(guān)系如下圖：上圖為相圖(x,y)，由上圖也可以得出，x與y自初始點(diǎn)(10，10)點(diǎn)演變至(100,0)，種群y隨x的增長(zhǎng)，雖然在開始時(shí)也有一定的增長(zhǎng)，但在種群x增長(zhǎng)至一定限度以后，種群y的競(jìng)爭(zhēng)劣勢(shì)開始顯現(xiàn)，隨著x的增長(zhǎng)，y逐漸減少，最終滅絕。（2）當(dāng)s1,s2保持不變：a.其他賦值不變，r1=1,r2=10時(shí)，x,y隨t的變化及y隨x變化的圖像如下：其他賦值不變，r1=10,r2=1時(shí)，x,y隨t的變化及y隨x變化的圖像如下：由上圖我們可以推斷出，r在生態(tài)學(xué)中代表種群的增長(zhǎng)率，而從圖中我們也可以知道，r對(duì)種群的競(jìng)爭(zhēng)能力并沒有起到?jīng)Q定性的作用。雖然在初始過(guò)程中，r對(duì)種群的數(shù)量有顯著的影響，但隨著時(shí)間的推移，達(dá)到穩(wěn)定態(tài)時(shí)，種群的數(shù)量并沒有因?yàn)閞的變化而改變。 b.其他賦值不變，n1=10,n2=100時(shí)，x,y隨t的變化及y隨x變化的圖像如下：其他賦值不變，n1=100,n2=10時(shí)，x,y隨t的變化及y隨x變化的圖像如下：n在生態(tài)學(xué)中的意義為，環(huán)境對(duì)該種群的最大容量，我們可以從圖中知道，n值得大小并不會(huì)改變種群的競(jìng)爭(zhēng)能力。而隨著時(shí)間的推移，競(jìng)爭(zhēng)優(yōu)勢(shì)者的數(shù)量逐漸趨于環(huán)境對(duì)其的最大容量，而競(jìng)爭(zhēng)劣勢(shì)者最終依舊會(huì)滅亡。 值得注意的現(xiàn)象是，當(dāng)初始值大于或等于最大容量時(shí)，無(wú)論競(jìng)爭(zhēng)是否具有優(yōu)勢(shì)，種群數(shù)量都會(huì)在最開始的時(shí)間內(nèi)有所下降。這在生態(tài)學(xué)上也是可以理解的。c.其他賦值不變，x0=5,y0=50時(shí)，x,y隨t的變化及y隨x變化的圖像如下：其他賦值不變，x0=50,y0=5時(shí)，x,y隨t的變化及y隨x變化的圖像如下： 有以上兩圖，我們可以初步判定，初始值的大小只能決定開始階段的競(jìng)爭(zhēng)優(yōu)勢(shì)，但隨著時(shí)間的推移，這種優(yōu)勢(shì)不斷減少，而最終達(dá)到穩(wěn)定狀態(tài)時(shí)與初始值并沒有關(guān)系。綜合以上幾種情況，我們可以推斷，s才是決定種群在生態(tài)中的競(jìng)爭(zhēng)能力的根本因素。當(dāng)s1=1.5,s2=0.7時(shí)a.其他賦值不變，r1=1,r2=10時(shí)，x,y隨t的變化及y隨x變化的圖像如下：其他賦值不變，r1=10,r2=1時(shí)，x,y隨t的變化及y隨x變化的圖像如下：同中情形相同，改變r(jià)值只能改變初始暫態(tài)過(guò)程的斜率，對(duì)最終的穩(wěn)態(tài)沒有作用。b.其他賦值不變，n1=10,n2=100時(shí)，x,y隨t的變化及y隨x變化的圖像如下：其他賦值不變，n1=100,n2=10時(shí)，x,y隨t的變化及y隨x變化的圖像如下：上圖可知，同情形相同，最大容量n不會(huì)對(duì)穩(wěn)定態(tài)的情況產(chǎn)生影響c.其他賦值不變，x0=5,y0=50時(shí)，x,y隨t的變化及y隨x變化的圖像如下：其他賦值不變，x0=50,y0=5時(shí)，x,y隨t的變化及y隨x變化的圖像如下：從上圖可以看出，初始值的差別只會(huì)影響進(jìn)入穩(wěn)定區(qū)域前的變化過(guò)程，與最終的穩(wěn)定狀態(tài)無(wú)關(guān)。與之前的情況是一樣的。(3) 當(dāng)s1=0.8,s2=0.7時(shí)：a.當(dāng)r1=r2=1,n1=n1=100,x0=y0=10時(shí)，x,y隨t的變化及y隨x變化的圖像如下：可以看出這種情況與上文中情況有明顯的差別： 當(dāng)參數(shù)s1和s2均小于1時(shí)，兩種種群的數(shù)量均在初始時(shí)快速上升。兩條曲線上升至某一點(diǎn)后x增長(zhǎng)率下降，直至種群數(shù)量降低至一穩(wěn)定值；而y則持續(xù)上升，但也無(wú)法達(dá)到最大容量；之后兩種種群穩(wěn)定保持在一定數(shù)量，二者平衡共存，生態(tài)系統(tǒng)達(dá)到穩(wěn)定。b.r1=10,r2=1,n1=100,n2=100,x0=y0=10時(shí)，x,y隨t的變化及y隨x變化的圖像如下：r1=1,r2=10,n1=100,n2=100,x0=y0=10時(shí)，x,y隨t的變化及y隨x變化的圖像如下：r1=1,r2=1,n1=100,n2=100,x0=5,y0=50時(shí)，x,y隨t的變化及y隨x變化的圖像如下：r1=1,r2=1,n1=100,n2=100,x0=50,y0=5時(shí)，x,y隨t的變化及y隨x變化的圖像如下：由圖像可知，與之前的情形相同，改變?cè)鲩L(zhǎng)率r和初始值同樣只能改變達(dá)到穩(wěn)定前的初始過(guò)程，并不能影響穩(wěn)定態(tài)的情況。c.r1=1,r2=1,n1=10,n2=100,x0=y0=10時(shí)，x,y隨t的變化及y隨x變化的圖像如下： r1=1,r2=1,n1=100,n2=10,x0=y0=10時(shí)，x,y隨t的變化及y隨x變化的圖像如下：由上圖可知，改變最大容量會(huì)影響到穩(wěn)定態(tài)的情況，由此可得出，當(dāng)s均小于1時(shí)，環(huán)境對(duì)種群的最大容量n和s共同決定穩(wěn)定態(tài)的情況。當(dāng)s1=1.5,s2=1.7時(shí)a.當(dāng)r1=r2=1,n1=n1=100,x0=y0=10時(shí)，x,y隨t的變化及y隨x變化的圖像如下：x最終達(dá)到最大容量，y最終完全滅絕。b.r1=10,r2=1,n1=100,n2=100,x0=y0=10時(shí)，x,y隨t的變化及y隨x變化的圖像如下：r1=1,r2=10,n1=100,n2=100,x0=y0=10時(shí)，x,y隨t的變化及y隨x變化的圖像如下：由以上兩圖可以看出，固有增長(zhǎng)率的差別也會(huì)影響到最終競(jìng)爭(zhēng)的結(jié)果：當(dāng)y0的固有增長(zhǎng)率超過(guò)x0很多時(shí)，最終y就有可能達(dá)到最大容量而x可能會(huì)滅絕。c.r1=1,r2=1,n1=10,n2=100,