實(shí)驗(yàn)4 常微分方程數(shù)值解

1、實(shí)驗(yàn) 4 常微分方程數(shù)值解實(shí)驗(yàn)?zāi)康?：1、 用 MATLAB軟件掌握求微分方程數(shù)值解的方法，并對(duì)結(jié)果作初步分析；2、 通過(guò)實(shí)例學(xué)習(xí)用微分方程模型解決簡(jiǎn)化的實(shí)際問(wèn)題?！绢}目 1】放射性廢物的處理：有一段時(shí)間，美國(guó)原子能委員會(huì)（現(xiàn)為核管理委員會(huì)）處理濃縮放射性廢物時(shí)， 把它們裝入密封性能很好的圓桶中， 然后扔到水深 300ft 的大海中。這種做法是否會(huì)造成放射性污染，自然引起生態(tài)學(xué)家及社會(huì)各界的關(guān)注。 原子能委員會(huì)一再保證圓桶非常堅(jiān)固， 絕不會(huì)破漏， 這種做法是絕對(duì)安全的。然而一些工程師們卻對(duì)此表示懷疑， 認(rèn)為圓桶在海底相撞時(shí)有可能發(fā)生破裂。于是雙方展開了一場(chǎng)筆墨官司。究竟誰(shuí)的意見正確呢？原子能委

2、員會(huì)使用的是 55gal 的圓桶，裝滿放射性廢物時(shí)的圓桶重量為 527.436lbf ，在海水中受到的浮力為 470.327lbf 。此外，下沉?xí)r圓桶還要受到水的阻力， 阻力與下沉速度成正比， 工程師們做了大量實(shí)驗(yàn)， 測(cè)得其比例系數(shù)為 0.08lbf s/ft 。同時(shí)，大量破壞性實(shí)驗(yàn)發(fā)現(xiàn)當(dāng)圓桶速度超過(guò) 40ft/s 時(shí)，就會(huì)因?yàn)榕c海底沖撞而發(fā)生破裂。（ 1） 建立解決上述問(wèn)題的微分方程模型。（ 2） 用數(shù)值和解析兩種方法求解微分方程，并回答誰(shuí)贏了這場(chǎng)官司?！締?wèn)題分析】 本題應(yīng)該求出圓桶下沉到海底時(shí)的速度，比較其與40ft/s的大小，來(lái)決定誰(shuí)最后贏了這場(chǎng)官司。已知水深 h=300ft ，裝滿廢物

3、的圓桶質(zhì)量為 m=527.436lbf ，在海水中受到的浮力為 B=470.327lbf ，圓桶受到的海水阻力 D與下降速度 v 成正比，比例系數(shù)K=0.08lbfs/ft。設(shè) t 時(shí)刻圓桶下落深度為 y，圓桶速度為 v，加速度為 a，所受重力為 G=mg，所受浮力為 f ，所受阻力為 K*v, 則圓桶所受合力為：F=G-f-Kv依據(jù)牛頓第二定律有：由此可以利用常微分方程求解，其中v(0)=0,y(0)=0。【程序運(yùn)算】一、數(shù)值法記 y(1)=y,y(2)=v, y=(y(1),y(2) T. 編寫如下 M文件：functiondy=fangshe(t,y)%建立名為 fangshe 的函數(shù)

4、M文件m=527.436*0.4536;G=527.436*0.4536*9.8;f=470.327*0.4536*9.8;Kv=0.08*0.4536*9.8*y(2)/0.3048;dy=y(2);(G-f-Kv)/m%以向量性格是表示方程編寫如下運(yùn)行程序： ts=0:2000; y0=0,0; %給定初值 opt=odeset(reltol,1e-3,abstol,1e-6); t,y=ode45(fangshe,ts,y0,opt); plot(t,y(:,2)，r) %畫出 v(t)-t圖象%給定精確度 gtext(t); gtext(v(t);運(yùn)行結(jié)果如下： plot(t,y(:,

5、1)， axis(0 25 0 350),grid;%畫出 y(t)-t圖像 hold on, ymax=300*0.3048; depth=linspace(ymax,ymax,2001); plot(t,depth,b)， axis(0 25 0 350),grid;%畫出最深高度圖像gtext(t),gtext(y(t),gtext(ymax)運(yùn)行結(jié)果如下：由圖像知， t 大概等于 13s時(shí)，圓桶到達(dá)海底，結(jié)合以上兩圖可以得知，圓桶落入海底時(shí)，會(huì)因?yàn)樗俣冗^(guò)大而破裂。二、解析法將模型分離變量得：對(duì)兩邊分別不定積分得：代入已知數(shù)據(jù)得：當(dāng) v=40ft/s 時(shí)，可利用上式算出 y值（即 v的一

6、階常微分方程）：齊次解的形式為，其中 q=-4.87676A =-217.59 ;特解為 v=217.59 ；綜上，故有：V(t)=13.8m/s=44.0ft/s【結(jié)果分析】結(jié)合題目的約束條件作圖分析，判斷下沉距離 S=300ft 時(shí)速度是否大于 40ft/s 或判斷速度達(dá)到 40ft/s 時(shí)下沉的距離是否大于 300ft （我選擇用速度是否大于 40ft/s 判斷），從而判斷圓桶與海底沖撞時(shí)是否會(huì)破裂。從依據(jù)數(shù)值解圖的中可見， 圓桶還沒(méi)沉到海底時(shí)， 其速度已達(dá)到規(guī)定的圓桶不破裂的最大速度。由解析解算出的數(shù)據(jù)也說(shuō)明，當(dāng)下落高度達(dá)到 300ft 時(shí)，速度已達(dá) 44.0ft/s 超過(guò)規(guī)定速度，故

7、圓桶會(huì)破裂。工程師們贏了官司。從速度 - 時(shí)間曲線可以看到，開始時(shí)圓桶的速度隨著時(shí)間逐漸增大，加速度逐漸減小，可見圓桶做的是加速度逐漸減小的變加速運(yùn)動(dòng)。時(shí)間充分大后（ t1500s ）后，圓桶的加速度減至0，速度趨于一個(gè)定值，約為275m/s?！绢}目 2】一只小船渡過(guò)寬為 d 的河流，目標(biāo)是起點(diǎn) A 正對(duì)著的另一岸 B 點(diǎn)。已知河水流速為 v1 與船在靜水中的速度之比為k。（ 1） 建立描述小船航線的數(shù)學(xué)模型，求其解析解；（ 2） 設(shè) d=100m，v1=1m/s， v2=2m/s，用數(shù)值解法求渡河所需的時(shí)間、任意時(shí)刻小船的位置及航行曲線，作圖，并與解析解比較。（ 3） 若流速 v1=0，0.

8、5 ，1.5 ， 2（ m/s），結(jié)果將如何？【問(wèn)題分析】假設(shè)航行者并不知道水流速度。設(shè)時(shí)間為t ，小船坐標(biāo)為（ x,y ）。據(jù)題目所給公式，可以列如下方程：初值為：（ x,y ）=(0,d)【程序運(yùn)算】（ 1） 運(yùn)用解析法將上兩式相除得：令 x/y=p y=0:-0.01:-100; x=xiaochuan(y,0.5); plot(x,y);%畫出圖像 gtext(x/m) gtext(y/m)(2) 運(yùn)用解析法記 x(1)=x,x(2)=y, x=(x(1),x(2) T。functiondx=xiaochuan2(t,x,v1 ,v2)s=(x(1)2+x(2)2)0.5;dx=v1-

9、x(1)/s*v2;-x(2)/s*v2;%建立函數(shù)文件 ts=0:0.1:100; x0=0,-100; %給定初始值 opt=odeset(reltol,1e-3,abstol,1e-9);%給定精確度 t,x=ode15s(xiaochuan2,ts,x0,opt,1,2); plot(t,x)%畫出函數(shù)圖像 gtext(t/s)gtext(x,/m)gtext( y/m )t/sx(t)/my(t)/m00-1000.10.099885-99.80.20.199576-99.60.30.29907-99.40.40.398364-99.20.50.497459-990.60.59635

10、2-98.80.70.695043-98.60.80.793532-98.40.90.891817-98.210.989899-981.11.087753-97.800166.60.23916-0.0022666.70.139166-0.0007666.80.039167-6.00E-05從上圖可以看出，當(dāng) t=66.8s時(shí)， y 接近 0，故可認(rèn)為小船到岸。再運(yùn)行如下程序： plot(x(:,1),x(:,2); xlabel(x/m) ylabel(y/m)（ 3）改變流速 v1 ，得到以下幾種情況： v1=0 ，此時(shí)為靜止的湖面。顯然船將沿直線直接從A 到B，渡河所需時(shí)間為100/2=5

11、0s v1=0.5m/s時(shí)，結(jié)果如下：跟 v1=1時(shí)的圖形基本相似，僅在趨勢(shì)上有微小不同。在t=53.3 s時(shí)，y=-0.0099m可以認(rèn)為船到達(dá)對(duì)岸。 v1=1.5m/s時(shí)，運(yùn)算結(jié)果如下：當(dāng) t=114.0s時(shí)， y=0.000m，認(rèn)為船到達(dá)對(duì)岸。與 v1=0.5， 1m/s 時(shí)比較，小船幾乎都是在 30s， 40s時(shí)間范圍內(nèi)到達(dá) x 方向的最大值，然后逆流駛向 B 點(diǎn)。 v1 越小，小船在逆水回到 B 的時(shí)間就越短。v1=2m/s 時(shí)，運(yùn)算結(jié)果如下：可以發(fā)現(xiàn)曲線的變化趨勢(shì)明顯不同于前幾個(gè)，且最終小船無(wú)法到達(dá) B 點(diǎn)，而到達(dá)了距離 B點(diǎn)50米的位置。因?yàn)?v1=v2 ，故 k=1, 將 k

12、代入得一拋物線，將d=100,y=0代入解得 x=50 ，這與實(shí)驗(yàn)結(jié)果相符?！窘Y(jié)果分析】通過(guò)改變 v1 的值我們可以得到以下結(jié)論：1、在 0v1v2的情況下，隨著 v1的增加，小船到達(dá) B 點(diǎn)的時(shí)間增加，與此同時(shí)，在整個(gè)過(guò)河的過(guò)程中，小船的最大水平距離也將增加。2、在 v1 v2時(shí)，小船始終無(wú)法到達(dá) B 點(diǎn)，且 v1的取值越大，最終小船到達(dá)對(duì)岸時(shí)距離 B 點(diǎn)的水平距離越遠(yuǎn)?！绢}目 3】?jī)煞N群相互競(jìng)爭(zhēng)模型如下：其中 x(t),y(t)分別為甲乙兩種群的數(shù)量， r 1,r 2 為它們的固有增長(zhǎng)率， n1,n 2 為它們的最大容量， s1 的含義是，對(duì)于供養(yǎng)甲的資源而言，單位數(shù)量乙（相對(duì)于n2）的消

13、耗為單位數(shù)量甲（相對(duì)于n1 ）消耗的 s1 倍，對(duì) s2 可作相應(yīng)的解釋 .該模型無(wú)解析解，使用數(shù)值解法研究以下問(wèn)題：(1) 設(shè) r 1=r 2=1, n 1=n2=100, s 1=0.5, s 2=2, 初值 x0 =y0=10，計(jì)算 x(t), y(t)畫出它們的圖形及相圖（ x, y）, 說(shuō)明時(shí)間 t 充分大以后 x(t),y(t) 的變化趨勢(shì) ( 人們今天看到的已經(jīng)是自然界長(zhǎng)期演變的結(jié)局) 。(2) 改變 r 1, r 2 , n 1, n 2,x 0, y 0,但 s1, s 2 保持不變 ( 或保持 s11),計(jì)算并分析所得結(jié)果；若s1=1.5(1), s2=0.7(1),再分析

14、結(jié)果。由此你能得到什么結(jié)論，請(qǐng)用生態(tài)學(xué)上的含義作解釋。(3) 試驗(yàn)當(dāng) s1=0.8(1),s2=0.7(1),s2=1.7(1)時(shí)又會(huì)有什么結(jié)果。能解釋這些結(jié)果嗎？【問(wèn)題分析】此題是很好的生態(tài)模擬題， 用matlab 來(lái)說(shuō)明不同的物種特性在競(jìng)爭(zhēng)中的優(yōu)勢(shì)與劣勢(shì)。(1) 、【運(yùn)算程序】functiondx =compete( t,x )r1=1;r2=1;n1=100;n2=100;s1=0.5;s2=2;dx=r1.*x(1).*(1-x(1)./n1-s1.*x(2)./n2);r2.*x(2).*(1-s2.*x(1)./n1-x(2)./n2)ts=0:0.1:30;x0=10;10;t,

15、x=ode45(compete,ts,x0);t,xplot(x(:,1),x(:,2),grid,gtext(x)gtext(y)pauseplot(t,x(:,2),grid,gtext(t)gtext(y)pauseplot(t,x(:,1),grid,gtext(t)gtext(x)運(yùn)行結(jié)果如下：txy010100.100010.880510.71200.200011.823511.44540.300012.830912.19620.400013.904412.95940.500015.045313.72940.600016.254414.49990.700017.532115.265

16、10.800018.878316.01730.900020.292416.74701.000021.772917.44541.100023.318218.10431.200024.926018.71581.300026.593319.27311.400028.316919.76971.500030.092720.19981.600031.916320.55811.700033.782720.84011.800035.686221.04181.900037.621321.15962.000039.583921.19142.100041.567621.13952.200043.565721.006

17、92.300045.571820.79712.400047.579720.51412.500049.583620.16212.600051.577919.74592.700053.557419.27072.800055.517118.74182.900057.452418.16543.000059.358917.54763.100061.232416.89543.200063.069316.21583.300064.865915.51633.400066.619114.80513.500068.326114.09023.600069.984513.37603.700071.592212.666

18、13.800073.147011.96453.900074.647311.27504.000076.091710.60094.100077.47949.94514.200078.81029.31034.300080.08398.69854.400081.30128.11154.500082.46297.55074.600083.57007.01684.700084.62296.51004.800085.62236.03024.900086.56955.57745.000087.46585.15125.100088.31264.75145.200089.11144.37725.300089.86

19、394.02785.400090.57173.70245.500091.23673.39985.600091.86083.11875.700092.44602.85795.800092.99392.61645.900093.50642.39336.000093.98512.18766.100094.43201.99826.200094.84871.82426.300095.23681.66436.400095.59801.51766.500095.93391.38306.600096.24601.25946.700096.53581.14636.800096.80471.04286.90009

20、7.05400.94837.000097.28500.86217.100097.49880.78367.200097.69660.71217.300097.87960.64697.400098.04870.58757.500098.20500.53327.600098.34930.48387.700098.48250.43897.800098.60530.39817.900098.71850.36108.000098.82280.32748.100098.91890.29698.200099.00750.26928.300099.08900.24408.400099.16400.22118.5

21、00099.23310.20038.600099.29660.18148.700099.35490.16448.800099.40850.14898.900099.45780.13499.000099.50300.12229.100099.54460.11079.200099.58280.10029.300099.61780.09079.400099.65000.08219.500099.67950.07449.600099.70650.06739.700099.73130.06109.800099.75410.05529.900099.77490.050010.000099.79410.04

22、5210.100099.81160.040910.200099.82760.037010.300099.84230.033510.400099.85580.030310.500099.86810.027510.600099.87940.024910.700099.88970.022510.800099.89920.020410.900099.90790.018411.000099.91580.016711.100099.92310.015111.200099.92970.013711.300099.93580.012411.400099.94130.011211.500099.94640.01

23、0111.600099.95100.009211.700099.95530.008311.800099.95920.007511.900099.96270.006812.000099.96590.006112.100099.96890.005612.200099.97160.005012.300099.97410.004612.400099.97630.004112.500099.97840.003712.600099.98030.003412.700099.98200.003112.800099.98360.002812.900099.98500.002513.000099.98630.00

24、2313.100099.98750.002013.200099.98860.001913.300099.98960.001713.400099.99050.001513.500099.99140.001413.600099.99210.001213.700099.99280.001113.800099.99350.001013.900099.99400.000914.000099.99460.000814.100099.99500.000814.200099.99550.000714.300099.99590.000614.400099.99620.000614.500099.99660.00

25、0514.600099.99690.000514.700099.99720.000414.800099.99740.000414.900099.99760.000315.000099.99780.000315.100099.99800.000315.200099.99820.000315.300099.99840.000215.400099.99850.000215.500099.99860.000215.600099.99880.000215.700099.99890.000215.800099.99900.000115.900099.99910.000116.000099.99920.00

26、0116.100099.99920.000116.200099.99930.000116.300099.99940.000116.400099.99940.000116.500099.99950.000116.600099.99950.000116.700099.99960.000116.800099.99960.000116.900099.99960.000017.000099.99970.000017.100099.99970.000017.200099.99970.000017.300099.99970.000017.400099.99980.000017.500099.99980.00

27、0017.600099.99980.000017.700099.99980.000017.800099.99980.000017.900099.99990.000018.000099.99990.000018.100099.99990.000018.200099.99990.000018.300099.99990.000018.400099.99990.000018.500099.99990.000018.600099.99990.000018.700099.99990.000018.800099.99990.000018.900099.99990.000019.000099.99990.00

28、0019.1000 100.00000.000019.2000 100.00000.000019.3000 100.00000.000019.4000 100.00000.000019.5000 100.00000.000019.6000 100.00000.000019.7000 100.00000.000019.8000 100.00000.000019.9000 100.00000.000020.0000 100.00000.000020.1000 100.00000.000020.2000 100.00000.000020.3000 100.00000.000020.4000 100.

29、00000.000020.5000 100.00000.000020.6000 100.00000.000020.7000 100.00000.000020.8000 100.00000.000020.9000 100.00000.000021.0000 100.00000.000021.1000 100.00000.000021.2000 100.00000.000021.3000 100.00000.000021.4000 100.00000.000021.5000 100.00000.000021.6000 100.00000.000021.7000 100.00000.000021.8

30、000 100.00000.000021.9000 100.00000.000022.0000 100.00000.000022.1000 100.00000.000022.2000 100.00000.000022.3000 100.00000.000022.4000 100.00000.000022.5000 100.00000.000022.6000 100.00000.000022.7000 100.00000.000022.8000 100.00000.000022.9000 100.00000.000023.0000 100.00000.000023.1000 100.00000.

31、000023.2000 100.00000.000023.3000 100.00000.000023.4000 100.00000.000023.5000 100.00000.000023.6000 100.00000.000023.7000 100.00000.000023.8000 100.00000.000023.9000 100.00000.000024.0000 100.00000.000024.1000 100.00000.000024.2000 100.00000.000024.3000 100.00000.000024.4000 100.00000.000024.5000 10

32、0.00000.000024.6000 100.00000.000024.7000 100.00000.000024.8000 100.00000.000024.9000 100.00000.000025.0000 100.00000.000025.1000 100.00000.000025.2000 100.00000.000025.3000 100.00000.000025.4000 100.00000.000025.5000 100.0000-0.000025.6000 100.0000-0.000025.7000 100.0000-0.000025.8000 100.0000-0.00

33、0025.9000 100.0000-0.000026.0000 100.0000-0.000026.1000 100.0000-0.000026.2000 100.0000-0.000026.3000 100.0000-0.000026.4000 100.0000-0.000026.5000 100.0000-0.000026.6000 100.0000-0.000026.7000 100.0000-0.000026.8000 100.00000.000026.9000 100.00000.000027.0000 100.00000.000027.1000 100.00000.000027.

34、2000 100.00000.000027.3000 100.00000.000027.4000 100.00000.000027.5000 100.00000.000027.6000 100.00000.000027.7000 100.00000.000027.8000 100.00000.000027.9000 100.00000.000028.0000 100.00000.000028.1000 100.00000.000028.2000 100.00000.000028.3000 100.00000.000028.4000 100.00000.000028.5000 100.00000

35、.000028.6000 100.00000.000028.7000 100.00000.000028.8000 100.00000.000028.9000 100.00000.000029.0000 100.00000.000029.1000100.00000.000029.2000100.00000.000029.3000100.00000.000029.4000100.00000.000029.5000100.00000.000029.6000100.00000.000029.7000100.00000.000029.8000100.00000.000029.9000100.00000.

36、000030.0000100.00000.0000最初甲群數(shù)量隨時(shí)間的推移而增長(zhǎng)。但甲種群種數(shù)量隨時(shí)間推移達(dá)到飽和，保持在 100 左右，不再隨時(shí)間變化。最初乙群數(shù)量隨時(shí)間的推移而增長(zhǎng)。對(duì)于乙種群，隨著時(shí)間的推移，在第2個(gè)單位時(shí)間時(shí)，種群數(shù)量達(dá)到最大值，之后種群數(shù)量便開始下降，直至滅絕?！窘Y(jié)果分析 】在甲和乙的種群競(jìng)爭(zhēng)中， 由于供養(yǎng)甲所需的資源少。 供養(yǎng)甲所需的資源少。 最初環(huán)境適宜食物充足， 兩種群數(shù)量都在都增長(zhǎng)， 但經(jīng)過(guò)一定時(shí)間的資源消耗， 甲乙的環(huán)境、食物等競(jìng)爭(zhēng)變得激烈， 在種群競(jìng)爭(zhēng)中， 占有優(yōu)勢(shì)的甲種群數(shù)量增加至達(dá)飽和，而乙種群數(shù)量大幅下降直至滅絕。（ 2） 改變 r 1 , r2, n1 , n2 ,x0 , y0 , s1 , s2 保持不變 ),計(jì)算并分析所得結(jié)果；A. 只改變固有增長(zhǎng)率r1=0.5; r2=0.8時(shí)，源程序如下：functiondx =compete( t,x )r1=0.5;r2=0.8;n1=100;n2=100;s1=0.5;s2=2;dx=r1.*x(1).*(1-x(1)./n1-s1.*x(2)./n2);r2.*x(2).*(1-s2.*x(1)./n1-x(2)./n2)ts=0:0.1:30;x0=10;10;t,x=ode45(c