數(shù)學(xué)建模流水問題的實(shí)驗(yàn)報(bào)告.doc

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