1、清華大學數學實驗報告實驗4：常微分方程數值解化學工程系 分2 安振華 2012011837【實驗目的】1. 練習數值微分的計算；2. 掌握用MATLAB 軟件求微分方程初值問題數值解的方法；3. 通過實例學習用微分方程模型解決簡化的實際問題；4. 了解歐拉方法和龍格-庫塔方法的基本思想和計算公式及穩定性等概念。【實驗內容】【實驗四：習題5】射性廢物的處理：有一段時間，美國原子能委員會（現為核管理委員會）處理濃縮放射性廢物時，把它們裝入密封性很好的圓桶中然后扔到水深300ft的大海中。這種做法是否會造成放射性污染，自然引起生物學家及社會各界的關注。原子能委員會一再保證，圓桶非常堅固，絕不會破漏，

2、這種做法是絕對安全的。然而一些工程師們卻對此表示懷疑，認為圓桶在和海底相撞時有可能發生破裂。于是雙方展開了一場筆墨官司。究竟誰的意見正確呢？原子能委員會使用的是55gal的圓桶，裝滿放射性廢物時的圓桶重量為527.436 lbf,在海水中受到的浮力為470.327 lbf。此外，下沉時圓桶還要受到海水的阻力，阻力與下沉速度成正比，工程師們做了大量實驗，測得其比例系數為0.08lbf s/ft。同時，大量破壞性試驗發現當圓桶速度超過40ft/s時，就會因與海底沖撞而發生破裂。（1）建立解決上述問題的微分方程模型；（2）用數值和解析兩種方法求解微分方程，并回答誰贏了這場官司。【分析】為了便于分析，

3、首先，假設圓桶下沉過程無障礙且做直線運動，圓桶質量為m，裝滿放射性廢物時的圓桶重量為M,下沉速度為v，極限速度為u，加速度為a, 重力加速度為g，海水深度為H，圓桶下沉時距海面距離為h.，比例系數為k。GFf受力分析：圓桶下落過程受重力G、浮力F、阻力f三個力的作用，受力情況如圖所示圓桶下落過程分析：由于圓桶剛開始時所受重力和浮力均不變，而所受阻力與速度成正比關系，速度增加時，阻力也隨之增加，因此圓桶在此過程做加速度減小的加速直線運動。當阻力f增加到G-F-f=0，即加速度a=0時，阻力和速度同時達到最大，此后速度不變，即圓桶做勻速直線運動。碰撞過程分析：通過以上分析我們可以猜想圓桶與海底碰撞

4、時的速度有兩種可能。第一，圓桶在與海底發生碰撞前就已達到最大速度，也就是說，此時圓桶是以最大速度與海底發生碰撞的；第二，圓桶有可能未達到最大速度就已海底發生碰撞。【解答】由已知條件換算得：M=527.436 lbf= 239.241kg，m=470.327lbf=213.337kg，g=9.800m/s2，k=0.08 lbf.s/ft=0.119 kg.s/m，u=40ft/s=12.192m/s，H=300ft=91.440m 通過以上方程，我們可以計算出圓桶到達海底時的速度v，然后與極限速度u作比較，當v>u時，圓桶與海底發生碰撞后破裂，當v<u時，圓桶與海底發生碰撞后不會破

5、裂。構建函數：function dv=fun1(t,v) %構造“時間-速度”函數M=239.241;m=213.337;g=9.8;k=0.119; %根據題意選取模型參數dv=(M-m-k*v)*g/M; %表示微分方程（1）MATLAB程序:clf ;ts=linspace(0,1500,15000); %考察圓桶下落25分鐘內速度變化情況v0=0; %假設初速度為0t,v=ode45(fun1,ts,v0); %解微分方程（1）t,v; %輸出結果plot(t,v,'b'); %作出“時間-速度”圖象grid on; %加網格線xlabel('t-時間'

6、); %橫坐標表示時間ylabel('v-速度'); %縱坐標表示速度title('圓桶下落速度隨時間變化情況');運行結果：結果分析：圖像顯示結果與前面分析大致相同，即圓桶在此過程做加速度減小的加速運動，直到當a=0時，v達到最大。即 此時 這就是前面對圓桶與海底碰撞過程分析猜想的第一種可能情況，即圓桶有可能以最大速度與海底發生碰撞，此時圓桶的速度遠遠大于極限速度，圓桶與海底發生碰撞后會發生破裂。但由于該計算過程并未考慮海底深度，也就是說圓桶有可能并未達到最大速度就與海底發生碰撞，那么在這種情況下圓通是否能夠達到極限速度還無法判斷，因此還必須分析速度與圓桶下沉

7、距離的關系，求出圓桶下落到海底時的速度，然后與極限速度作比較，進而得出結論。clf;ts=linspace(0,300,3000); %考察圓桶下落5分鐘內速度變化情況v0=0; %假設初速度為0t,v=ode45(fun1,ts,v0); %解微分方程（1）t,v; %輸出結果h=cumtrapz(t,v); %利用梯形積分公式求圓桶距海面距離H=91.44;u=12.192; %根據題意選取模型參數plot(h,v,'r'); %作出“速度-距離”圖象hold on;plot(0,100,12.192,12.192,'-g'); %設定水平線-表示極限速度p

8、lot(91.44,91.44,0,15,'-b'); %設定鉛垂線-表示海底grid on; axis(0,100,0,15); %設定圖像范圍text(91.44,6,'海底'); %標注“海底”的位置text(35,12.192,'極限速度'); %標注“極限速度”的位置xlabel('h-圓桶到海面的距離'); ylabel('v-速度'); title('速度隨距離變化情況');MATLAB程序:運行結果：結果分析：由圖像可知，到達海底時，圓通速度已經超過了極限速度，因此，圓桶與海底碰撞后

9、會發生破裂，里面的放射性廢物會泄露。 從而解得這與數值法求解微分方程的結果相一致。另外,由（4）式，令h(t)=91.440m，可解出t=13.2697s，代入（3）式解得，v=13.6348m/s>12.192m/s，即圓桶到達海底時的速度大于極限速度，因此圓桶與海底碰撞后會發生破裂。由上面的分析可知，工程師們的結論是正確的。【結論】通過以上兩種方法得出結論：圓桶到達海底時的速度大于極限速度，因此圓桶與海底碰撞后會發生破裂，所以工程師將贏得這場官司。 【實驗四：習題6】一只小船渡過寬為d 的河流，目標是起點A 正對著的另一岸B 點。已知河水流速v1 與船在靜水中的速度v2之比為k。（1

10、） 建立描述小船航線的數學模型，求其解析解；（2） 設d=100m，v1=1m/s，v2=2m/s，用數值解法求渡河所需的時間、任意時刻小船的位置及航行曲線，作圖，并與解析解比較。（3） 若流速v1=0，0.5，1.5，2（m/s），結果將如何？【分析】小船運動速度矢量分解圖如下：vV2V1根據實際情況，分析可得到較為符合實際的數學模型：船并不知道水速，否則只要調整合適角度，即可沿直線通過即可。現小船不知道水速，應該通過船的控制系統對速度進行xy兩個方向的分解，可列出常微分方程如下：dxdt=v1-v2xx2+y2dydt=-v2yx2+y2其初值條件為x(0)=0,y(0)=-d.1、常微分

11、方程解析解兩式相除得dxdy=-v1x2+y2v2y+xy=-kxy2+12+xy令u=xy 則有dxdy=ydudy+u代入ydudy=-ku2+12分離變量積分可得y-k=c(u2+12+u)代入初值，可得c=-d-k整理得x=-d2(y-d)1-k-(y-d)1+k2、 MATLAB數值解【解答】用matlab求上述方程組的數值解。記x(1)=x, x(2)=y, x=(x(1), x(2)T。編寫M文件如下：function dx=guohe(t,x)v1=1;v2=2;s=(x(1)2+x(2)2)0.5;dx=v1-x(1)*v2/s;-x(2)*v2/s;endMATLAB程序：

12、x0=0,-100;ts=0:0.1:70; %經過試驗設定終值為70t,x=ode15s(guohe,ts,x0); %求解方程組plot(t,x),grid; %畫圖同時做出x、y方向的位移title('圖7.分位移-時間')xlabel('t/s');ylabel('x,y/m');gtext('x(t)');gtext('y(t)');pause;plot(x(:,1),x(:,2),grid; %畫出航行曲線title('航行曲線');xlabel('x');ylabel(

13、'y');t,x(:,1),x(:,2) %輸出數據表運行結果： 由圖7可以看出，在前30s中，y(t)有較好的的線性，因此，在這個時間段內，小船沿y軸方向的為勻速。從圖8可以看出，船合速度的方向為軌跡的切線方向。在剛開始階段，可以認為傳順水而行，小船在剛開始的速度方向的確變化不大，與上圖相吻合。小船在達到x最大值時開始轉向，著這一時刻附近的合速度最小。部分數據：t/sx/my/m00-100.00000.10000.0999-99.80000.20000.1996-99.60000.30000.2991-99.40000.40000.3984-99.20000.50000.4

14、975-99.00000.60000.5964-98.80000.70000.6951-98.60000.80000.7936-98.40000.90000.8919-98.20001.00000.9900-98.00001.10001.0879-97.80001.20001.1856-97.60001.30001.2831-97.40001.40001.3804-97.20001.50001.4775-97.00011.60001.5744-96.80011.70001.6710-96.60011.80001.7675-96.400165.10001.5882-0.098765.20001.

15、4886-0.086665.30001.3889-0.075465.40001.2892-0.064965.50001.1894-0.055265.60001.0896-0.046365.70000.9898-0.038265.80000.8899-0.030865.90000.7900-0.024366.00000.6901-0.018566.10000.5902-0.013566.20000.4902-0.009366.30000.3902-0.005966.40000.2903-0.003366.50000.1903-0.001466.60000.0903-0.000366.70000.

16、0000066.80000.00000解析解與數值解的比較MATLAB程序：y0=0:-1:-100;x0=50.*(y0./(-100).0.5-(y0./(-100).1.5);plot(x(:,1),x(:,2),'r',x0,y0,'b'),grid;legend('數值解','解析解')title('航行曲線比較');xlabel('x');ylabel('y');運行結果：由圖9可以看出，數值解得到的曲線與解析解的到的曲線幾乎完全重合，說明數值解得到的結果還是很準確的。改

17、變流速（1）v1=0m/s時,更改程序中的數值，可得如下分位移-時間圖像及航行曲線：由圖像可知，v1=0m/s時，水的流速為0，即小船的合速度等于v2，因此小船沿y軸方向勻速前進。此時小船的位置的解析解為x=0.5*y-0.5*y=0，所以小船過河軌跡為沿y軸方向的一條直線。小船沿y方向的速度dydt=-v2yx2+y2=-v2y|y|=v2，故小船勻速過河。過河所用時間為50s。（2）v1=0.5s時，由數據表可知，t=53.3s時，y=-0.0131,可近似認為已到達岸邊。將V1=0.5m/s與v1=1m/s時相比，河水流速減小后，船在x方向上的分速度減小了，因此小船在x方向上行駛的最大值

18、減少。這與途中x最大值減小了一半左右相符。與v1=0時比較，小船過河花的時間并沒有增加很多，但路程長了很多，這說明小船的平均合速度大于v2。這是因為實際速度為河水與船水速度的矢量疊加結果。（3）v1=1.5m/s時，由數據表可知，t=114.3s時，x=0.0007m,y=0.000m，可認為小船到達對岸。從圖中可以看出，小船在30s，40s時間范圍內達到x方向的最大值，然后逆流駛向B點。所以v1越小，船回到B的時間就越短，這一點與常識相符。將解析式x=d2(y-d)1-k-(y-d)1+k求導并令其為零，可以求出x在y=-100*2k1-k1+k處，取得最大值，所以k在0,1區間內越大，x取得最大值的點就越靠近x軸。（4）v1=2m/s時利用程序可以求出，小船最終在下游50m處靠岸，不能到達正對岸。也就是說，無論多長時間小船都無法到達正對岸。由解析解也可以證明這一點。dxdt=v1-