數(shù)學建模流水問題的實驗報告.doc

數(shù)學建模實驗報告-流水問題一 問題描述 三個橫截面積為常數(shù)A，高分別為H1，H2，H3的水池內(nèi)都盛滿了水，都由池底一橫截面積為B的小孔放水。設水從小孔流出的速度為v（i）=sqrt（2*g*h（i），求水流空所需的時間。二 前提假設1. 假設在一段極微小的時間間隔dt內(nèi)，三個浴缸的高度變化速率以及三個排水口的排水速率是一個不變化的定值。 2. 排水速率僅與水池高度有關。3. 排水口的高度為水池最低處，即不會出現(xiàn)因水位低于排水口而無法排完水的現(xiàn)象。三 問題分析 將此問題抽象成數(shù)學問題：求初值分別為H1，H2，H3的函數(shù)h1，h2，h3隨時間變化的函數(shù)，以及他們變?yōu)?所需的時間（設水池1，2，3的流出速度為v1，v2，v3）。根據(jù)任何一個水池內(nèi)水量在一個時間微元內(nèi)的減少量等于流出量減去流入量可以得到如下關系：浴缸 1：。水池2：。兩邊取極限后得-dh2*A= ds2*B- dh1*A。注意到|dh1*A|= |ds1*B|,除以dt可得(-dh2/dt)*A= (v2-v1)*B，化簡并帶入v1的函數(shù)表達式可得，再代入h1的表達式可以得到如下的常微分方程。這是一個非線性常微分方程，難以得到解析解（并非不可求，可用待定系數(shù)法等求解，但是在此處求出解析解并不是建模的重點，因為即使h2存在解析解，h3也不一定存在，而對于求出排水時間圖，求出近似解更為重要），在這里我們采用計算方法中的一些數(shù)值計算方法求出幾組h2和t2，v2的近似解。浴缸 3：，由此遞推公式可得如下微分方程：，可以用Euler法求出近似解 。四 問題求解 對于h1，書上已給出解法，即用一個微分方程求解，在此不做累述。在代碼中，我們用一個數(shù)組arrayt記錄時間，一個數(shù)組h1記錄水位高度，然后對應畫圖。 本部分Matlab代碼如下： A=2;B=1;H=10;g=10;h=H;temp=0;t=0;tf=0i=1;h1=;h2=;h3=;temp=;tempt=;arrayt=;arrayt2=;arrayt3=;sqrtH=sqrt(H);h1=h1,H;arrayt=arrayt,t;while (h0)t=t+0.01;arrayt=arrayt,t;h=sqrt(H)-(B/2/A)*(sqrt(2*g )*t);h1=h1,h*h;end 對于h2，由于我們得到了一個遞推公式，并得到了一個常微分方程，本來可以直接求解，但由于此方程屬于不可解的一類，故僅求出近似解。由于第二步的數(shù)值會對后續(xù)的浴缸產(chǎn)生影響，我們選用了精度較高的Runge-Kutta法（龍格-庫塔法，以下簡稱RK法，具體方法參見數(shù)值計算方法第六章），其原理是改進了的歐拉法，將x的區(qū)間分成很多小間隔來分段應用歐拉法，matlab中給出了實現(xiàn)變步長Runge-Kutta法的函數(shù)：t,y=ODE45(odefun,tspan,y0)其中t，y用于輸出廣義時間與相空間向量，odefun為要求的微分方程，tspan為廣義時間區(qū)間，y0為初值向量。由于Runge-Kutta法是一種近似解法，僅在局部區(qū)間收斂，所以會有虛解及負解出現(xiàn)，而求解h3時會需要用到h2的解，所以我們要對得到的解進行篩選，只留下非負實數(shù)解。 本部分Matlab代碼如下： funt.M：function y=funt(t,x) A=2;B=1;g=10;H=10; if(x0)arrayt2=arrayt2,tempt(i);h2=h2,temp(i);i=i+1;end對于h3，浴缸3： -Dh3*A= Ds3*B- Dh2*A，由此遞推公式可得如下微分方程：由于只有三步，且得到的h2的點是離散的，Runge-Kutta法使用不便。此處我們采用簡單一些的歐拉法(Euler法，具體方法見數(shù)值計算方法第六章)進行數(shù)值求解。采用此方法可以得到h3的數(shù)值解，即h3關于時間的離散函數(shù)。 h3的解不需要篩選。本部分Matlab代碼如下：Euler法的實現(xiàn)： j=1;h3(1)=H/2;arrayt3(1)=0;temp3=0;while(ji-1) temp3=h3(j)+(h2(j)-h2(j+1)*(-1*B/A*(sqrt(2*g*h3(j)-sqrt(2*g*h2(j); h3=h3,temp3; arrayt3=arrayt3,arrayt2(j+1); j=j+1;end 至于求解水流完的時間，可采用遍歷h1，h2，h3數(shù)組的方法，直到找到0點所在，然后取出對應的t即可。在計算h3時要注意對模型進行修正，由于h3所有的分量均由h2參與迭代運算產(chǎn)生，故當h2為0之后，h3也停止了迭代。故在我們繪制的圖中h3沒有到0就停止了?？梢愿倪M為，當t時刻h2為0之后，h3采用h3(t)為初值，帶入h1的表達式進行計算來計算零點。本部分Matlab代碼如下：t1=0;t2=0;t3=0; i=length(arrayt);t1=arrayt(i); i=length(arrayt2);t2=arrayt2(i); i=length(arrayt3);t3=arrayt3(i)h3(i)t_extra=A/B*sqrt(2*h3(i)/g);t3=t3+t_extra; t1t2t3五 數(shù)值模擬結果取H0=10，討論以下初始情況下三個浴缸水位的關系：分別是三個浴缸初始量相等的情況，一多兩少的情況和一少二多的情況.H1,H2,H3,t1,t2,t3分別表示三浴缸的初始水量以及全部流出時間排水