數(shù)學建模作業(yè)、微分方程實驗、北京工業(yè)大學

1、2微分方程實驗1、微分方程穩(wěn)定性分析繪出下列自治系統(tǒng)相應的軌線，并標出隨 t增加的運動方向，確定平衡點, 并按穩(wěn)定的、漸近穩(wěn)定的、或不穩(wěn)定的進行分類：dx dydtx,y;dx(2)dtdydtdxx,dt2y；崇dx y，dt2x; dydtx+1,2y.解：(1)由 f (x) =x=0, f (y) =y=0;可得平衡點為(0,0),1 0.一系數(shù)矩陣A,求得特征值入1=1,入2=1;0 1p=-(入1+入2)=-2<0 , q=入1入2=1>0;對照穩(wěn)定性的情況表，可知平衡點(0, 0) 是不穩(wěn)定的。圖形如下：(2)如上題可求得平衡點為(0,0),特征值入1=-1,入2=2

2、;p=-(入1+入2)=-1<0 , q=入1入2=-2<0;對照穩(wěn)定性的情況表，可知平衡點(0, 0) 是不穩(wěn)定的。其圖形如下:(3)如上題可求得平衡點為(0,0),特征值入1=0 + 1.4142i,入2=0 - 1.4142i; p=-(入1+入2)= 0, q=X 1入2=1.4142>0;對照穩(wěn)定性的情況表，可知平衡點(0, 0)是不穩(wěn)定的。其圖形如下：(4)如上題可求得平衡點為(1,0),特征值入1=-1,入2=-2;p=-(入1+入2)= 3>0, q=入1入2=2>0;對照穩(wěn)定性的情況表，可知平衡點(1,0) 是穩(wěn)定的。其圖形如下：2、種群增長模型

3、一個片子上的一群病菌趨向于繁殖成一個圓菌落 .設病菌的數(shù)目為N,單位成員的增長率為U ,則由Malthus生長律有dN ri N ,但是，處于周界表面的 dt那些病菌由于寒冷而受到損傷，它們死亡的數(shù)量與N2成比例，其比例系數(shù)為r2, 求N滿足的微分方程.不用求解，圖示其解族.方程是否有平衡解，如果有，是否 為穩(wěn)定的？解：由題意很容易列出N滿足的微分方程：dN r1N 冰1f(N)dt令f(N)=0,可求得方程白兩個平衡點 N1=0,N2=r22/r12211進而求得d2N1252(r1r2N 2) (r1Nr2N 2)dt2令空220可求得N=T2 /4ri則 N=N1別有dN0口,2dt.2

4、2cdN cdNcdN cdN0, ；0,20 ；0,2-dtdtdtdt0。N=N2 N=r22/4r；可以把第一象限劃為三部分，且從下到上三部分中分則可以畫出N (t)的圖形，即微分方程的解族，如下圖所示:由圖形也可以看出，對于方程的兩個平衡點，其中N1=0是不穩(wěn)定的；N2二始/r；是穩(wěn)定的。3、有限資源競爭模型1926年Volterra提出了兩個物種為共同的、有限的食物來源而競爭的模型dx b 1(幾X h2X2)x dtdx2 M2(hlXi h2X2)X2dt假設b1殳，稱.為物種i對食物不足的敏感度， 12i(1)證明當x1(t0)>0時，物種2最終要滅亡；(2)用圖形分析方

5、法來說明物種 2最終要滅亡.解：(1)由上述方程組 f (x1) =b11('X1 h2x2Hxi=0,f (x2) =b22(1X1 h2X2)x2=0,可得方程白平衡點為 R (0,0), P1 (d,0),P2 (0,bM . 2 h203bi 0對平衡點Po (0,0),1X1h2 1X2h2 2X1h1 2X2b2hl 2x1h1 2X2則p=- (b1+b2) <0,所以該平衡點不穩(wěn)定對平衡點Pi(*，0)，系數(shù)矩陣Abih2 1x2h2 2x1bih2 blhihl 2x2b2 hi 2x1hi2x2則 p= n b2 理，q= n(b2 ”)，ii由題意. 邁，x

6、i(t0)>0 ,可以得出p>0,q>0,因此該平衡點是穩(wěn)定的 i 2即t 時，(xi(t),x2(t)(且,0),說明物種2最終要滅亡。ihi對平衡點P2 (0, -b2-),同理可以得到q<0,在該平衡點不穩(wěn)定。9 2最終要滅亡。00生，xi(t0)>0 ,可將第一象限分為2. ''在中間區(qū)域，Xi,X2都小于0,在取右邊，在,xi(t0)>0i 2(2)對于線性方程組bib2的條件下，物彳i(Xi h2X2)2(hlxl h2X2)在平面上匹配兩條直線l i和l 2,由題意立 i三個區(qū)域。在最左邊區(qū)域， Xi,X2都大于0;區(qū)域，Xi,

7、X2分別是大于0和小于0.,由各區(qū)域中x'i,X2的取值可得到如下圖形: 公t+ bi由圖也可以看出，隨著時間的增加，物種 1最終能達到穩(wěn)定值 ,ihi物種2最=£ , &為(3(0),重0.57050.00000.00000.5000-0.0000-0.00000.53520.00000.0000終要滅亡。4、蝴蝶效應與混沌解考慮Lorenz模型Xi (t)Xi(t)X2(t)X3(t)'_ _ 一一 一一 X2(t)X2(t)X3(t)X3(t)Xl(t)X2(t)X2(t)X3(t)其中=10, p =28, B =8/3,且初值為,xi (0) =X2

8、 (0) =0, x3 (0) 一個小常數(shù)，假設e =10-1°,且0&t&100。(1)用函數(shù)ode45求解，并畫出x2x1,x2x3,x3x1的平面圖;(2)適當?shù)卣{整參數(shù)6，p , B值，和初始值x1 (0), x2 (0) =0, 復一的工作，看有什么現(xiàn)象發(fā)生。解：(1)編寫Lorenz函數(shù)，function xdot=lorenz1(t,x,b,a,c) xdot=-b*x(1)+x(2)*x(3);-a*x(2)+a*x(3);-x(1)*x(2)+c*x(2)-x(3);對各參數(shù)賦值并用ode45函數(shù)求解，可得數(shù)值解: Columns 1 through

9、 900.12500.25000.37500.60570.640900.00000.0000-0.00000.00000.000000.00000.0000-0.00000.0000-0.00000.0000-0.0000-0.00000.00000.00000.00000.00000.0000Columns 10 through 180.67610.71140.74660.78180.83080.87970.97761.01050.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.000

10、00.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000Columns 19 through 271.04341.07631.10921.14211.17501.20791.27971.31860.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00010.00010.00000.00000.00000.00000.00000.00010.00020.0002Columns 28 through 361.357

11、51.39641.42461.45281.48101.50921.56561.59380.00000.00000.00000.00000.00000.00000.00000.00000.00020.00030.00040.00050.00080.00110.00210.00290.00040.00060.00080.00120.00160.00230.00450.0063Columns 5590 through 559899.575199.592199.609099.626099.646299.666499.706999.733816.945716.526116.201015.985415.8

12、97816.034817.193318.7894-3.3551-3.7119-4.1098-4.5568-5.1636-5.8601-7.5438-8.8677-5.3476-5.9274-6.5941-7.3519-8.3766-9.5370-12.1944 -14.0725-0.00000.92860.00000.00000.00001.24090.00000.00000.00011.53740.00000.00150.003299.686716.4476-6.6527-10.8209Columns 5599 through 560799.760899.787899.894099.9099

13、21.218624.470936.714937.0528-10.3044-11.7459-13.4302 -12.7919-15.7415-16.8309-10.1013-8.100599.814899.830628.262930.5354-12.9932-13.5361-16.9274-16.364899.846599.862332.653334.4645-13.8800-13.9854-15.3302-13.870799.878235.8466-13.8351-12.0805Columns 5608 through 5613 99.9257 99.941699.956236.8957 36

14、.323935.5248-11.9606 -10.9899 -10.0237 -6.2187-4.5596-3.281899.970899.9854 100.000034.543533.448132.2971-9.0332-8.0563-7.1223-2.2590-1.4835-0.9275x2x1,x2x3,x3x1的平面圖分別如下:x2-x10IIIIIUI-20-15-10-5051015202x3304 口302010TO11111-30-20-100102030(2)令參數(shù)，p , B值各減1,初始值X1 (0), X2 (0)不變，X3 (0) =10-8分別得到得到x2x1,x2

15、x3,x3x1的平面圖如下：x2x150可以看出，參數(shù)6，p , B值各減1,初始值xi (0), X2 (0)不變，X3 (0) 數(shù)值變?yōu)閄3 (0) =108,參數(shù)和初始值很小的改變，就會導致最后圖形發(fā)生很大 的變化。5、用微分方程考察共振現(xiàn)象設物體沿X軸運動(如圖所示)其平衡位置取為原點0,物體的質量為1,在 時間t物體的位置為x(t)其所受的恢復為(如彈性力等)與物體所在位置的坐標成 正比，即k2x,其中常數(shù)k稱為恢復系數(shù)，運動過程所受的阻力(由于介質及摩擦 等)設與速度成正比，即2hdx, h>0,稱為阻尼系數(shù)。dt(1)根據(jù)Newton第二定律，建立相應的微分方程.不妨設初始

16、位置為1, 初始速度為0,取k=2, h=0(當h = 0稱為簡諧振動的方程)和h=0.1,用Matlab軟件 得到相應的數(shù)值解，并在t-x平面上畫出x (t)的圖形。(2)如果物體還受到附加外力的干擾，且外力是一個依據(jù)時間t的函數(shù)f(t)(設f (t)=B sinwt),建立相應的微分方程(該方程稱為強迫振動方程).在上述參 數(shù)不變的情況下，取振幅 B=1,分別取 w=1, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2,2.4, 2.6, 2.8, 3.0,用Matlab軟件得到相應的數(shù)值解，并在t-x平面上畫出x(t)的圖形。(3)分別對上述圖形講行分析，并解釋為什么會出現(xiàn)這些

17、現(xiàn)象。k2xdx2h 出d 2x解：(1)根據(jù)Newton弟二止律，F(xiàn)=ma,可得微分方程：m了 mg dt由題意知邊界值：x(0)=1,x'(0)=0,令 yi=x, y2=一 , dt y2可將二階方程轉化為：dy2 mg k2yl 2hy2dtyi(0) i,y2(0) 0其中，m=1; g=9.8; k=2當h=0時，由matlab解得數(shù)值解：Columns 1 through 900.00000.00000.00000.00000.00010.00010.00020.00021.00001.00001.00001.0000Columns 10 through 181.0000

18、1.00001.00001.00001.00000.00040.00060.00540.01080.00090.00110.00220.00320.00431.00001.00001.00011.0003Columns 19 through 271.00001.00001.00001.00001.00010.01620.02160.21550.29560.02710.05410.08120.10830.13531.00081.00141.13261.2461Columns 100 through 1081.00211.00851.01911.03391.05288.23218.35938.98

19、529.09398.46158.56368.66588.76798.87663.50193.21611.52141.3034Columns 109 through 1172.95052.66412.36882.07691.78349.20259.31429.912410.00009.42609.53779.64949.73719.82471.13931.03451.64121.86341.00011.03831.14671.27731.4438h=0時，x (t)圖形如下所示:3.53x 2.521.5h=o0123456t18910當h=0.1時，由matlab解得數(shù)值解:Columns 1

20、 through 900.00000.00020.00020.00000.00000.00000.00010.00011.00001.00001.00001.0000Columns 10 through 181.00001.00001.00001.00001.00000.00040.00060.00540.01080.00090.00110.00220.00320.00431.00001.00001.00011.0003Columns 19 through 271.00001.00001.00001.00001.00010.01620.02160.21560.29580.02710.05410

21、.08120.10830.13531.00081.00141.13081.24171.00211.00851.01901.03361.0523Columns 109 through 1178.56318.66908.77498.88098.98689.12239.25799.39349.52892.58572.45632.32942.21062.10491.99521.92131.88801.8962Columns 118 through 1219.64679.76459.882210.00001.93532.00162.09092.1979h=0.1 時，x(t)圖形如下所示:h=0 13.

22、5x 2.510(2)如果還受到外力干擾，則微分方程變?yōu)?d2x_, 2dxm2mgk x2h一Bsinwtdt2出將其化為一階方程組:dy1dtdy2dty1(0)y2,2,mg k y1 2hy2 Bsinwt1,(0) 0其中,m=1; g=9.8; k=2; B=1; h=0.1;用Matlab解得數(shù)值解:w=1時Columns 1 through 900.00000.00000.00000.00000.00010.00010.00020.00021.00001.00001.00001.00001.00001.00001.00001.00001.0000Columns 118 thro

23、ugh 1219.82079.88049.94021.89731.92101.9503w=1.2 時Columns 1 through 900.00000.000010.00001.98460.00000.00000.00010.00010.00020.00021.00001.00001.00001.00001.00001.00001.00001.00001.0000Columns 118 through 1219.82819.88549.94271.69991.75681.8199w=1.4 時Columns 1 through 900.00000.000010.00001.88860.00

24、000.00000.00010.00010.00020.00021.00001.00001.00001.00001.00001.00001.00001.00001.0000Columns 118 through 1259.62369.73099.83829.94549.95919.97279.986410.00002.24012.31712.40842.51032.52382.53732.55102.5647w=1.4 時Columns 1 through 900.00000.00000.00000.00000.00010.00010.00020.00021.00001.00001.00001

25、.00001.00001.00001.00001.00001.0000Columns 118 through 1219.78239.85499.927410.00002.11502.06762.02851.9979w=1.6 時Columns 1 through 900.00000.00000.00000.00000.00010.00020.00021.00001.00001.00001.00001.00001.00001.00001.00000.00011.0000Columns 118 through 1219.69339.79559.89780.80620.75030.756210.00

26、000.8231w=1.8 時Columns 1 through 9 00.00000.00020.00021.00001.00001.00001.00000.00000.00000.00000.00011.00001.00001.00001.00000.00011.0000Columns 118 through 1219.77529.85019.925110.00000.86291.08501.33831.6174w=2.0 時Columns 1 through 90.00011.000000.00000.00000.00000.00000.00010.00020.00021.00001.0

27、0001.00001.00001.00001.00001.00001.000010.00003.1243Columns 118 through 1219.75389.83599.91792.32782.60232.8705w=3.0 時Columns 1 through 900.00000.00000.00020.00021.00001.00001.00000.00000.00001.00001.00000.00011.00000.00011.00001.00001.000010.00002.4915如上圖所示,w分別為1.0; 1.2; .; 3.0時x (t)的圖形。Columns 118

28、 through 1219.74959.83309.91652.24142.33032.4145當h=0時，即無阻尼振動的情況下，得到t (x)圖形如下所示:w=l,2w=1 Aw=2 0051005100w=2 46 4 2 05w二2一66 4 2 0I0w 二6 4 2 005106 4 2 0(3)由上述圖形可以看出，不考慮外加作用力時，當沒有摩擦等阻力的時候， 簡諧運動一直進行下去，當有摩擦等阻力的時候，其振幅會逐漸的減小，最后趨 近于零，但是周期不會發(fā)生變化。在外加一個作用力之后，整個運動情況就發(fā)生了很大的變化。首先剛開始的 振幅和周期都變化很大，后來隨著時間的推移，他們也逐漸開始

29、趨于穩(wěn)定。但是， 穩(wěn)定之后的振幅和周期都各不相同， 與外力的頻率有直接的關系，這就是受迫振 動，不僅跟系統(tǒng)本身的性質有關，還和外界干擾的情況有很大關系?；謴拖禂?shù)k與w的值越接近時，隨著時間增加，物體的振動振幅會一直增大，當w=k=2時， 物體的振幅增大的速度最快，可以預見，在t趨于無窮大時，振幅也趨于無窮大， 這就是共振現(xiàn)象。也就是當外力與系統(tǒng)處于共振時，會引起振幅無限增大的振動， 這在機械和建筑中一般是必須嚴格避免的。6、人口模型與曲線擬合問題表2.1列出的是美國1790 1980年的人口統(tǒng)計表.(1)試用 Malthus 人口模型,按三段時間(1790-1850,1850-1910,191

30、0-1970)分別確定其增長率r。并將數(shù)據(jù)和不同r值的三段曲線畫在同一張圖上.(2)利用Logistic模型，重新確定固有增長率r和最大容量Nm,作圖，再 利用該模型得到的結果預測1990年的美國人口數(shù)。解：(1)分段研究，我們先求出增長率，編寫命令如下：t = 10 : 10 : 70;x = 3.9 5.3 7.2 9.6 12.9 17.1 23.2;p = poly攸t, log(x), 1); p得到結果p = 0.0296即17901850年時間段的增長率是0.0296同理可以得到另外兩段的增長率，1850-1910年 為0.022619101970 年為 0.0129每次畫出圖形

31、，使用以下命令plot(t, x, 'bx');hold on在以上每一次求時順便畫出圖形，得到最后的總圖如下2502001501005010203040506070r1 0,0296 2二口必 6 門=0.口 1290利用Logistics模型，將以上所有數(shù)據(jù)進行擬合，編寫程序如下t = 0 : 10 : 190;x = 3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5;f = inline('p(1)./(1+p

32、(2)*exp(-p(3).*t)','p','t');p = lsqcurvefit(f, 300, 50, 0.02, t, x);得到結果人口最大容量Nm是360.3560 (百萬)，增長率r是0.0234如下圖所示：輸入t_pre = 200;x_pre = p(1) ./ (1 + p(2) * exp(-p(3)*t_pre);可得最后的結果，即1990年的人口總數(shù)x_pre =241.7704 萬人。如下圖：一家環(huán)保餐廳用微生物將剩余的食物變成肥料餐廳每天將剩余的食物制成槳狀物并與蔬菜下腳及少量紙片混合成原料， 加入真菌菌種后放入容器內。真

33、菌 消化這此混合原料，變成肥料，由于原料充足，肥料需求旺盛，餐廳希望增加肥 料產量。由于無力購置新設備，餐廳希望用增加真菌活力的辦法來加速肥料生產 試通過分析以前肥料生產的記錄（如表2.2所示），建立反映肥料生成機理的數(shù)學 模型，提出改善肥料生產的建議。表*?以前肥料的生產記錄泥漿（磅）綠葉菜4磅）紙片（磅）喂入日期生成堆肥的Fl期31iJ1)D(J. 07.131990,08.1011279U1.9E)U,O7.17iggu.cs is7121019&0.07.241990.020203S201990.07.271990.08.22792S01990.08.101990.09.121

34、055201 )30.08.131990.09.1812115D1)yO.U8,21J1990.09.241103201990.08.221990.10.08表2；（續(xù)）以前肥料的生產記錄泥漿闕綠葉菜畸紙片（磅）喂入日期生成堆肥的日期824491951.04.301991.03.1857卬761991.05.0i21991.00.20771991.05.07199152(1199 1,05,1111991.06.28解答：1.問題條件假設1）每次堆肥的質量相等2）所給的幾次堆肥混合物的比例僅由當天的實際情況決定3）所有分離肥倉工作條件相同4）每磅蔬菜可提供的氧氣量相等5）細菌消耗的溶解氧氣完全

35、由蔬菜葉提供6）每天提供的廢物混合物中的化學成分大致相同7）廢物混合物在喂給細菌前混合均勻并保持良好的通氣環(huán)境。2.問題分析堆肥是利用微生物的分解作用將有機廢物轉化成無害穩(wěn)定形式的生物化學過 程，要提高堆肥產量方法之一就是通過增強細菌的生長繁殖能力以提高分解率。細菌群體的增長一般要經(jīng)歷延滯期，加速生長期、對數(shù)生長期、減數(shù)生長期、 穩(wěn)定期、加速死亡期和對數(shù)死亡期。另外，對數(shù)生長期培養(yǎng)基中所有養(yǎng)分都過剩， 細菌可以充分繁殖，其倍增速率恒定，取決于底物濃度、溫度、水活度、供氧量。 對于當前該餐廳來說底物濃度由每天的剩余食物、蔬菜葉和碎紙屑決定。碎紙屑是吸收水分的調理劑，微菌呼吸所需要的溶解氧由蔬菜也

36、提供， 水活度可以通過 測定相對濕度來決定，其關系式是：相對濕度：B=P/P0*100%水活度：aw=P/P0 其中P為該溶液蒸汽壓，P0為純水蒸汽壓，在這個堆肥系統(tǒng)中，可供微菌消耗的營養(yǎng)底物和溶解氧都是有限的，他們的消耗會對微菌生長率產生重要影響，一種微菌消耗營養(yǎng)底物的速率和底物濃度 之間的關系式為：ds/dt=KmSX/ （Ks+S）,式中ds/dt表示為底物的有效消耗率，x表示為微菌濃 度，Km表示為最大有效系數(shù)，在高濃度營養(yǎng)底物中最大的物料消耗率（物質質 量/微菌一天的質量）；Ks表示為半速系數(shù)（質量/體積），S表示為有限底物的濃 度（質量/體積）微菌生長過程是一個生物化學過程，其生長

37、率和溫度之間滿足公式：K=AeA（一Ea/RT）,式中K表示為反應速度常數(shù)，T表示為反應的絕對溫度，R表示氣體常 數(shù)，Ea表示反應活化能，A表示頻率分子。我們知道，各種微菌都有一個最適生長溫度，如果溫度控制在最適值是微菌生長速率最高。微菌生長對水活度也有一定的要求， 與微菌最高生長速率相對應 的有一最適水活度。優(yōu)化堆肥意味著盡可能少的時間內生產出高質量的肥料，參數(shù)的優(yōu)化依賴于所應用的系統(tǒng)。3模型的建立假設每天投入的廢料比是隨意的， 僅僅取決于當天的情況，首先從已知數(shù)據(jù)中得 到經(jīng)驗最佳比例，由于假設各次堆肥后肥料質量相同，因此堆肥時間較短就對應 了較好的廢料配比，將12組數(shù)據(jù)按其堆肥日期及完全堆

38、肥時間分成三組，每組 中的較優(yōu)比例見下表：分組每組中的最短天數(shù)比例14No.4: 26203:82:058No.5: 3379:28:0912No.9: 4982:44:9上述模型顯然過于粗糙 模型二：營養(yǎng)底物和氧氣是細菌生長的兩種底物，物耗公式為：dSi/dt=KmiXSi/（Ksi+Si）（1）i=1時代表營養(yǎng)底物中可降解得有機物。i=2時代表氧氣。1）在高濃度底物中，物料的轉化過程很迅速，進一步增加底物濃度就不再引起底物轉化率的提高，可假設S1>>Ks1,S2>>Ks2,則（1）式簡化為：dSi/dt*（1/X）=一Kmi,(2)這是關于底物濃度的零級反應，消去X

39、得到：dS1=Km1/Km2*dS2(3)積分后得：S1(t)=Km1/Km2*S2(t)+S1(0) Km1/Km2*S2(0)(4)2)在低濃度底物中，可假設 S1<<Ks1,S2<<Ks2,則(1)式可簡化為：dSi/dt*(1/X)= Kmi/Ksi*Si(5)消去 X 得：dS1/S1 *1/X=Km1*Ks2/(Km2*Ks1)* dS2/S2(6)積分后得：S1(t)=S1(0)S2(t)/s2(0)AKm1*Ks2/(Km2*Ks1)(7) 3）當Si=Ksi, i=1,2時有：dS1=Km1/Km2*dS2下面的步驟就和1相同了4結果分析：由上面的計算可知，當S>>Ks時，簡化方程為（4）且Km1=AeA（Ea/RT）,其 中，R=8.31J/(mol*K) ， Ea=51.3KJ/molT=25°