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

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