微積分模型ppt

1、模型分析和建立模型分析和建立先看如下三個例子先看如下三個例子1. 衛(wèi)星軌道的長度2. 射擊命中概率3. 人口增長模型1、衛(wèi)星軌道的長度 人造地球衛(wèi)星的軌道可以視為平面上的橢圓，中國第一顆人造衛(wèi)星近地點距離地球表面439km，遠(yuǎn)地點距離地球表面2384km，地球半徑為6371m，試求該衛(wèi)星的軌道長度？分析：2、射擊命中概率 炮彈射擊的目標(biāo)為一正橢圓形區(qū)域，當(dāng)瞄準(zhǔn)目標(biāo)的中心發(fā)射時，在縱多因素的影響下，彈著點與目標(biāo)中心有隨機偏差?？梢院侠淼丶僭O(shè)彈著點圍繞中心呈二維正態(tài)分布，且偏差在x方向和y方向相互獨立。若橢圓區(qū)域在x方向半軸長120m,y方向半軸80m，設(shè)彈著點偏差的均方差在x方向和y方向均為10

2、0m,試求 炮 彈 落 在 橢 圓 形 區(qū) 域 內(nèi) 的 概 率 ？分析：3、人口增長率20世紀(jì)美國人口統(tǒng)計數(shù)據(jù)如表6.5所示，計算表中這些年份的人口增長率。又已知某地區(qū)20世紀(jì)70年代的人口增長率如表6.6所示，且1970年人口為210萬，試估計該地區(qū)1980年的人口？分析:小小結(jié)結(jié)模型求解u 數(shù)值微分相關(guān)知識回顧u 數(shù)值積分相關(guān)知識u Matlab 相關(guān)命令u 程序?qū)崿F(xiàn)u 數(shù)值微分相關(guān)知識回顧hhxfhxfhhxfxfhxfhxfxfhhh2)()(lim)()(lim)()(lim)( 0001. 函數(shù)f(x)以離散點列給出時，而要求我們給出導(dǎo)數(shù)值，2. 函數(shù)f(x)過于復(fù)雜這兩種情況都要

3、求我們用數(shù)值的方法求函數(shù)的導(dǎo)數(shù)值微積分中，關(guān)于導(dǎo)數(shù)的定義如下：自然，而又簡單的方法就是，取極限的近似值，即差商 向前差商公式()( )( ),0,1iiif xhf xfxih 向后差商公式 中心差商公式( )()( ),1,2iiif xf xhfxih111()()()2f xhf xhfxh給定點列20)(,(iiixfx且2110,xxxxh)( ),( ),( 012xfxfxf求兩點差商公式：兩點差商公式：三點差商公式三點差商公式0000012201()3 ()4 ()(2 )21 3 ()4 ()()fxf xf xhf xhhf xf xf xxx1002201()()(2 )

4、21 ()()fxf xf xhhf xf xxx2000012201()()4 ()3 (2 )21 =()4 ()3 ()fxf xf xhf xhhf xf xf xxxu 數(shù)值積分相關(guān)知識回顧)()()(aFbFdxxfba關(guān)于積分，有Newton-Leibniz公式但是，在很多情況下，還是要數(shù)值積分：1、函數(shù)有離散數(shù)據(jù)組成2、F(x)求不出3、F(x)非常復(fù)雜定義數(shù)值積分如下：是離散點上的函數(shù)值的線性組合)()(0iniinxfafI稱為求積系數(shù)，與f(x)無關(guān)，與積分區(qū)間和積分點有關(guān)矩形公式矩形公式左矩形公式左矩形公式( )() ( )baf x dxba f a右矩形公式右矩形公

5、式( )() ( )baf x dxb a f b中矩形公式中矩形公式( )() ()2babaf x dxba f梯形公式梯形公式( )( )( )2babafx dxf af b辛普森辛普森(Simpson)公式公式( ) ( )4 ()( )62babaabf x dxf aff b復(fù)合梯形公式復(fù)合梯形公式11011101( )( ) ( )2()( )22kknbxaxknnkkkkkfx dxfx dxhhfxfxf afxf b復(fù)合辛普森復(fù)合辛普森(Simpson)公式公式111012112( ) ( )4()2()( )6=2nnbkakkkkkkbaf x dxf af xf

6、xf bxxx其中 將將 分為分為 等份，步長等份，步長 ，分點，分點, a bbahnn,0,1,kxakh knMatlab相關(guān)命令 x=linspace(0,pi,50);y=sin(x);R1=trapz(x,y)；% 梯形公式計算R2=quad(sin,0,pi,1.0e06)；% 辛普森公式計算結(jié)果R1=1.999314849324062R2=1.9999934965349640sin( )x dx程序?qū)崿F(xiàn)模型實例講解 森林救火問題 水箱的水流問題 森林救火問題森林救火問題森林失火后，要確定派出消防隊員的數(shù)量。森林失火后，要確定派出消防隊員的數(shù)量。隊員多，森林損失小，救援費用大；隊

7、員多，森林損失小，救援費用大；隊員少，森林損失大，救援費用小。隊員少，森林損失大，救援費用小。所以需要綜合考慮森林損失費和救援費與消防隊員所以需要綜合考慮森林損失費和救援費與消防隊員之間的關(guān)系，以總費用最小來決定派出隊員的數(shù)量。之間的關(guān)系，以總費用最小來決定派出隊員的數(shù)量。問題問題分析分析問題問題記隊員人數(shù)記隊員人數(shù)x, 失火時刻失火時刻t=0, 開始救火時刻開始救火時刻t1, 滅火時刻滅火時刻t2, 時刻時刻t森林燒毀面積森林燒毀面積B(t). 損失費損失費f1(x)是是x的減函數(shù)的減函數(shù), 由燒毀面積由燒毀面積B(t2)決定決定. 救援費救援費f2(x)是是x的增函數(shù)的增函數(shù), 由隊員人數(shù)

8、和救火時間決定由隊員人數(shù)和救火時間決定.找到找到恰當(dāng)?shù)那‘?dāng)?shù)膞，使，使f1(x), f2(x)之和最小之和最小 關(guān)鍵是對關(guān)鍵是對B(t)作出合理的簡化假設(shè)作出合理的簡化假設(shè).問題問題分析分析失火時刻失火時刻t=0, 開始救火時刻開始救火時刻t1, 滅火時刻滅火時刻t2, 畫出時刻畫出時刻 t 森林燒毀面積森林燒毀面積B(t)的大致圖形的大致圖形t1t20tBB(t2)分析分析B(t)比較困難比較困難,轉(zhuǎn)而討論森林燒毀轉(zhuǎn)而討論森林燒毀速度速度dB/dt. （其中（其中dB/dt表示表示單位時間燒毀的面單位時間燒毀的面積）積）模型假設(shè)模型假設(shè) 3）f1(x)與與B(t2)成正比，系數(shù)成正比，系數(shù)c

9、1 (燒毀單位面積損失費）燒毀單位面積損失費） 1）0 t t1, dB/dt 與與 t成正比，系數(shù)成正比，系數(shù) (火勢蔓延速度）火勢蔓延速度） 2）t1 t t2, 降為降為 - x ( 為隊員的平均滅火為隊員的平均滅火速度）速度） 4）每個）每個隊員的單位時間滅火費用隊員的單位時間滅火費用c2, 一次性費用一次性費用c3 (運運送消防隊員和器材等一次性支出送消防隊員和器材等一次性支出)假設(shè)假設(shè)1）的解釋的解釋 rB火勢以失火點為中心，火勢以失火點為中心，均勻向四周呈圓形蔓延，均勻向四周呈圓形蔓延，半徑半徑 r與與 t 成正比成正比面積面積 B與與 t2成正比，成正比， dB/dt與與 t成

10、正比成正比.xbtt12202)()(tdttBtB模型建立模型建立dtdBb0t1tt2x假設(shè)假設(shè)1）,1tbxcttxcxftBcxf31222211)()(),()(目標(biāo)函數(shù)目標(biāo)函數(shù)總費用總費用)()()(21xfxfxC假設(shè)假設(shè)3）4）xttt112假設(shè)假設(shè)2）)(222212212xttbt0dxdCxcxxtcxtctcxC3122121211)(22)(模型建立模型建立目標(biāo)函數(shù)目標(biāo)函數(shù)總費用總費用模型求解模型求解求求 x使使 C(x)最小最小231221122ctctcx結(jié)果解釋結(jié)果解釋 / 是火勢不繼續(xù)蔓延的最少隊員數(shù)是火勢不繼續(xù)蔓延的最少隊員數(shù)dtdBb0t1t2tx其中其中

11、 c1,c2,c3, t1, , 為已知參數(shù)為已知參數(shù)模型模型應(yīng)用應(yīng)用c1,c2,c3已知已知, t1可估計可估計, c2 x c1, t1, x c3 , x 結(jié)果結(jié)果解釋解釋231221122ctctcxc1燒毀單位面積損失費燒毀單位面積損失費, c2每個每個隊員單位時間滅火費隊員單位時間滅火費, c3每個每個隊員一次性費用隊員一次性費用, t1開始救火時刻開始救火時刻, 火火勢蔓延速度勢蔓延速度, 每個每個隊員平均滅火隊員平均滅火速度速度.為什么為什么? ? , 可可設(shè)置一系列數(shù)值設(shè)置一系列數(shù)值由模型決定隊員數(shù)量由模型決定隊員數(shù)量x練習(xí)題 在森林救火模型中，如果考慮消防隊員的滅火速度 與

12、開始救火時的火勢b有關(guān)，試假設(shè)一個合理的函數(shù)關(guān)系，重新求解模型。程序?qū)崿F(xiàn) 美國某州的用水管理機構(gòu)要求各社區(qū)提供以每小時多少加美國某州的用水管理機構(gòu)要求各社區(qū)提供以每小時多少加侖計的用水量以及每天所用的總水量。許多社區(qū)沒有測量流入侖計的用水量以及每天所用的總水量。許多社區(qū)沒有測量流入或流出水塔水量的裝置，只能代之以每小時測量水塔中的水位，或流出水塔水量的裝置，只能代之以每小時測量水塔中的水位，其誤差不超過其誤差不超過5%。需要注意的是，當(dāng)水塔中的水位下降到最低。需要注意的是，當(dāng)水塔中的水位下降到最低水位水位L時，水泵就自動向水塔輸水直到最高水位時，水泵就自動向水塔輸水直到最高水位H，此期間不能，

13、此期間不能測量水泵的供水量，因此，當(dāng)水泵正在輸水時不容易建立水塔測量水泵的供水量，因此，當(dāng)水泵正在輸水時不容易建立水塔中水位和用水量之間的關(guān)系。水泵每天輸水一次或兩次，每次中水位和用水量之間的關(guān)系。水泵每天輸水一次或兩次，每次約約2小時。小時。水箱的水流問題水箱的水流問題 試估計任何時刻（包括水泵正在輸水時間）從水塔試估計任何時刻（包括水泵正在輸水時間）從水塔流出的水流量流出的水流量f(t)，并估計一天的總用水量。已知該水塔，并估計一天的總用水量。已知該水塔是一個高為是一個高為40ft（英尺），直徑為（英尺），直徑為57ft（英尺）的正圓（英尺）的正圓柱，表柱，表6-1給出了某個小鎮(zhèn)一天水塔水

14、位的真實數(shù)據(jù)，給出了某個小鎮(zhèn)一天水塔水位的真實數(shù)據(jù)，水位降至約水位降至約27.00ft水泵開始工作，水位升到水泵開始工作，水位升到35.50ft停止停止工作。（注：工作。（注：1ft（英尺）（英尺）=0.3024m（米）（米）表表6.1 某小鎮(zhèn)某天水塔水位某小鎮(zhèn)某天水塔水位時間時間/s水位水位/0.01ft時間時間/s水位水位/0.01ft031754663633503316311049953326066353054539363167106192994572543087139372947605743012179212892645542927212402850685352842252232795

15、71854276728543275275021269732284269779254水泵開動水泵開動35932水泵開動水泵開動82649水泵開動水泵開動39332水泵開動水泵開動8596834753943535508995333974331834459327033402.2 問題分析問題分析 模型中所指流量可視為單位時間內(nèi)流出水的模型中所指流量可視為單位時間內(nèi)流出水的體積體積。由于水塔是正圓柱形，橫截面積是常數(shù)，所以在水泵由于水塔是正圓柱形，橫截面積是常數(shù)，所以在水泵不工作時段，流量很容易根據(jù)水位相對時間的變化算不工作時段，流量很容易根據(jù)水位相對時間的變化算出。問題的難點在于如何估計水泵供水時段

16、的流量。出。問題的難點在于如何估計水泵供水時段的流量。 水泵供水時段的流量只能靠供水時段前后的流量經(jīng)水泵供水時段的流量只能靠供水時段前后的流量經(jīng)插值或擬合得到。作為用于插值或擬合的原始數(shù)據(jù)，我插值或擬合得到。作為用于插值或擬合的原始數(shù)據(jù)，我們希望水泵不工作時段的流量越準(zhǔn)確越好。這此流量大們希望水泵不工作時段的流量越準(zhǔn)確越好。這此流量大體上可由兩種方法計算，一是直接對表體上可由兩種方法計算，一是直接對表6-1中的水量用中的水量用數(shù)值微分算出各時段的流量，用它們擬合其它時刻或連數(shù)值微分算出各時段的流量，用它們擬合其它時刻或連續(xù)時間的流量；二是先用表中數(shù)據(jù)擬合水位一時間函數(shù)，續(xù)時間的流量；二是先用表

17、中數(shù)據(jù)擬合水位一時間函數(shù)，求導(dǎo)數(shù)即可得到連續(xù)時間的流量。求導(dǎo)數(shù)即可得到連續(xù)時間的流量。 有了任何時刻的流量，就不難計算一天的總用水有了任何時刻的流量，就不難計算一天的總用水量。其實，水泵不工作時段的用水量可以由測量記錄量。其實，水泵不工作時段的用水量可以由測量記錄直接得到，由表直接得到，由表6-1中下降水位乘以水搭的截面積就是中下降水位乘以水搭的截面積就是這一時段的用水量。這個數(shù)值可以用來檢驗數(shù)據(jù)插值這一時段的用水量。這個數(shù)值可以用來檢驗數(shù)據(jù)插值或擬合的結(jié)果?；驍M合的結(jié)果。 在具體給出本問題的解答之前，先介紹一個簡單在具體給出本問題的解答之前，先介紹一個簡單的數(shù)據(jù)插值方法。的數(shù)據(jù)插值方法。2.

18、3 拉格朗日插值拉格朗日插值1、線性插值、線性插值)(xf假設(shè)已知假設(shè)已知在區(qū)間在區(qū)間,ba上的兩個結(jié)點上的兩個結(jié)點10, xx和它們的函數(shù)值和它們的函數(shù)值 10 xxx )()()(10 xfxfxf求一個一次多項式求一個一次多項式xaaxP101)(，使得多項式，使得多項式)(1xP在結(jié)點上滿足條件在結(jié)點上滿足條件 1, 0),()(1ixfxPii 這種插值方法稱為線性插值方法（也稱兩點插值）。這種插值方法稱為線性插值方法（也稱兩點插值）。 可以求出：可以求出： )()()(101001011xfxxxxxfxxxxxP2、拋物插值、拋物插值)(xf已知已知在區(qū)間在區(qū)間,ba上的三個結(jié)點

19、上的三個結(jié)點210,xxx和它們的函數(shù)值和它們的函數(shù)值 210 xxxx )()()()(210 xfxfxfxf求一個次數(shù)不超過求一個次數(shù)不超過2的多項式的多項式22102)(xaxaaxP，使得它在結(jié)點上滿足條件，使得它在結(jié)點上滿足條件2, 1, 0),()(2ixfxPii這種插值方法稱為拋物線插值法，這種插值方法稱為拋物線插值法， 可求出：可求出：)()( )()( )()()( )()( )()(121012002010212xfxxxxxxxxxfxxxxxxxxxP)()( )()( )(2021210 xfxxxxxxxx3、n次拉格朗日插值次拉格朗日插值,ba假設(shè)取區(qū)間假設(shè)取

20、區(qū)間 上的上的1n個結(jié)點個結(jié)點nxxx,10，并且已知函數(shù)，并且已知函數(shù))(xf在這此點的函數(shù)值在這此點的函數(shù)值 nxxxxx210 )()()()()(210nxfxfxfxfxf現(xiàn)在求一個次數(shù)不超過現(xiàn)在求一個次數(shù)不超過n的多項式的多項式nnnxaxaxaaxP2210)(，使得滿足條件，使得滿足條件 nixfxPiin, 1, 0),()( 這種插值方法稱為這種插值方法稱為n次多項式插值（或稱代數(shù)插值），次多項式插值（或稱代數(shù)插值）， 利用拉格朗日插值插值方法可得利用拉格朗日插值插值方法可得)()()()()()()()()(,2, 21, 10, 0nnnnnnnxfxlxfxlxfxl

21、xfxlxP0( )0,1,njknjkjj kxxlxknxx其中上述多項式稱為上述多項式稱為 n次拉格朗日（次拉格朗日（Lagrange）插值多項式，）插值多項式， 函數(shù)函數(shù)nkxlnk, 1, 0),(,稱為拉格朗日插值基函數(shù)。稱為拉格朗日插值基函數(shù)。當(dāng)當(dāng)n=1,2時，時，n次拉格朗日（次拉格朗日（Lagrange）插值多項式即）插值多項式即為線性插值多項式和拋物插值多項式。為線性插值多項式和拋物插值多項式。4321x例例6.1 已知函數(shù)發(fā)已知函數(shù)發(fā)f(x)的函數(shù)表如下：的函數(shù)表如下： 31 . 28 . 15 . 2)(xf求其拉格朗日插值多項式，并求求其拉格朗日插值多項式，并求)5

22、. 2(f的近似值。的近似值。 解解 由于給出了由于給出了4個插值結(jié)點，所以可做出次數(shù)不超個插值結(jié)點，所以可做出次數(shù)不超過過3的拉格朗日插值多項式。的拉格朗日插值多項式。43132361)41 ( )31 ( )21 ()4( )3( )2()(233 , 0 xxxxxxxl6219421)42( )32( )22()4( )3( )2()(233 , 1xxxxxxxl472721)43( )33( )23()4( )3( )2()(233 , 2xxxxxxxl161161)44( )34( )24()4( )3( )2()(233 , 3xxxxxxxl 將上列將上列4式代入式代入n=

23、3 的拉格朗日插值公式，可得所的拉格朗日插值公式，可得所要求的插值多項式為要求的插值多項式為6 . 4933. 29 . 00667. 0)(233xxxxP將將x=2.5代入可得代入可得f(2.5)的近似值為的近似值為1.8496。 拉格朗日插值法適合節(jié)點較少的情況，當(dāng)節(jié)點較多拉格朗日插值法適合節(jié)點較少的情況，當(dāng)節(jié)點較多的大范圍高次插值的逼近效果往往并不理想且當(dāng)插值的大范圍高次插值的逼近效果往往并不理想且當(dāng)插值結(jié)點增加時，計算越來越繁。為了提高精度和減少計結(jié)點增加時，計算越來越繁。為了提高精度和減少計算，還有牛頓插值法下、三次樣條插值等，具體可參算，還有牛頓插值法下、三次樣條插值等，具體可參

24、閱有關(guān)書籍。閱有關(guān)書籍。1一維插值命令一維插值命令interp1的具體使用格式的具體使用格式y(tǒng)y=interp1(x,y,xx,method) 其中其中x,y是插值結(jié)點的數(shù)據(jù)向量，如果是插值結(jié)點的數(shù)據(jù)向量，如果y是矩陣，是矩陣，則對矩陣則對矩陣y的每一列相對的每一列相對x進(jìn)行插值，進(jìn)行插值，xx是待求函數(shù)值是待求函數(shù)值的插值結(jié)點向量，可以缺省。的插值結(jié)點向量，可以缺省。method是可選項，是可選項，說明插值使用的方法。對于一維插值，說明插值使用的方法。對于一維插值，MATLAB提供提供可選的方法有：可選的方法有：nearest，linear，spline，cubic，它們分別表示最近插值、線

25、性插值、三次樣條插值和它們分別表示最近插值、線性插值、三次樣條插值和三次插值。三次插值。2.4 MATLAB軟件實現(xiàn)插值法軟件實現(xiàn)插值法 MATLAB軟件提供了專門做各種插值的命令：軟件提供了專門做各種插值的命令：interp1（一維插值），（一維插值），interp2（二維插值），（二維插值），interp3（三維插值），（三維插值），interpn（n維插值），維插值），spline（樣條插值）等。（樣條插值）等。 2二維插值命令二維插值命令interp2的具體使用格式的具體使用格式zz=interp2(x,y,z,xx,yy,method) 該指令的意思是根據(jù)數(shù)據(jù)向量該指令的意思是根據(jù)數(shù)

26、據(jù)向量x，y，z按按method指定的方法來做插值，然后將指定的方法來做插值，然后將xx，yy處插值處插值函數(shù)的插值結(jié)點向量，如果函數(shù)的插值結(jié)點向量，如果xx，yy在插值范圍之內(nèi)，在插值范圍之內(nèi)，則返回值在則返回值在zz中，否則返回值為空中，否則返回值為空NaN。method是插值方法可選項，具體要求同一是插值方法可選項，具體要求同一維插值的情況。維插值的情況。 該命令還有以下幾種省略格式：該命令還有以下幾種省略格式：zz=interp2(z,xx,yy)zz=interp2(x,y,z,xx,yy)zz=interp2(z,ntimes)3三維插值命令三維插值命令interp3的具體使用格式

27、的具體使用格式vi=interp3(x,y,z,v,xi,yi,zi,method) 它的具體含義跟前面的一、二維插值是相似的，在它的具體含義跟前面的一、二維插值是相似的，在此不作解釋，讀者可在此不作解釋，讀者可在MATLAB工作空間中用工作空間中用help interp3命令獲得。命令獲得。4樣條插值命令樣條插值命令spline的具體使用格式的具體使用格式y(tǒng)y=spline(x,y,xx)它的意思等同于命令它的意思等同于命令yy=interp1(x,y,xx,spline)例例6.2 在用外接電源給電容器充電時，電容器兩端在用外接電源給電容器充電時，電容器兩端的電壓的電壓V將會隨著充電時間將

28、會隨著充電時間t發(fā)生變化，已知在某一次發(fā)生變化，已知在某一次實驗時，通過測量得到下列觀測值，分別用拉格朗日實驗時，通過測量得到下列觀測值，分別用拉格朗日插值法、分段線性插值法、三次樣條插值法畫出插值法、分段線性插值法、三次樣條插值法畫出V隨隨著時間著時間t變化的曲線圖，分別計算當(dāng)時間變化的曲線圖，分別計算當(dāng)時間t=7s時，三時，三種插值法各自算得電容器兩端電壓的近似值。種插值法各自算得電容器兩端電壓的近似值。1295 . 64321st 4 .101 .106 . 90 . 92 . 83 . 72 . 6vV解解:由于由于MATLAB沒有提供現(xiàn)成的拉格朗日插值命令，沒有提供現(xiàn)成的拉格朗日插值

29、命令，我們可以編寫一個函數(shù)我們可以編寫一個函數(shù)lglrcz.m來完成，其他兩種插值來完成，其他兩種插值法可用現(xiàn)成的命令。法可用現(xiàn)成的命令。 用用MATLAB軟件進(jìn)行三種插值計算的程序為軟件進(jìn)行三種插值計算的程序為Three_interpmethods.m。程序程序lagrange.m：function y=lagrange(x0,y0,x) n=length(x0)； m=length(x) ； for i=1:m z=x(i) ； s=0.0； for k=1:n p=1.0； for j=1:n if j=k p=p*(z-x0(j)/(x0(k)-x0(j) ； end end s=p*

30、y0(k)+s； end y(i)=s； end程序為程序為Three_interpmethods.mt=1,2,3,4,6.5,9,12 ；v=6.2,7.3,8.2,9.0,9.6,10.1,10.4 ；t0=0.2:0.1:12.5；lglr=lagrange(t,v,t0) ；laglr=lagrange(t,v,7) ；fdxx=interp1(t,v,t0) ；fendxx=interp1(t,v,7) ；scyt=interp1(t,v,t0,spline) ；sancyt=interp1(t,v,7,spline)plot(t,v,*,t0,lglr,r,t0,fdxx,g,t

31、0,scyt,b)gtext(lglr)gtext(fdxx)gtext(scyt)執(zhí)行結(jié)果是執(zhí)行結(jié)果是laglr= 9.52988980716254fendxx= 9.70000000000000sancyt= 9.67118039327444圖形如圖圖形如圖6.1所示。所示。圖中曲線圖中曲線lglr表示拉格朗日插值曲線，表示拉格朗日插值曲線，scyt表示三次樣表示三次樣條插值曲線，條插值曲線，fdxx表示分段線性插值曲線。表示分段線性插值曲線。 2.5 問題求解問題求解 為了表示方便，我們將為了表示方便，我們將2.1節(jié)問題中所給表節(jié)問題中所給表6-1中的中的數(shù)據(jù)全部化為國際標(biāo)準(zhǔn)單位（表數(shù)據(jù)

32、全部化為國際標(biāo)準(zhǔn)單位（表6-2），時間用小時），時間用小時（h），高度用米（），高度用米（m）：）：表表6-2 一天內(nèi)水塔水位記錄一天內(nèi)水塔水位記錄(單位：米單位：米)時間時間/h水位水位/m時間時間/h水位水位/m09.6812.9510.210.929.4513.889.941.849.3114.989.652.959.1315.909.413.878.9816.839.184.988.8117.938.925.908.6919.048.667.008.5219.968.437.938.3920.848.228.978.2222.02水泵開動水泵開動9.98水泵開動水泵開動22.96水泵開

33、動水泵開動10.93水泵開動水泵開動23.8810.5910.9510.8224.9910.3512.0310.5025.9110.181模型假設(shè)模型假設(shè)（1）流量只取決于水位差，與水位本身無關(guān)，故由物）流量只取決于水位差，與水位本身無關(guān)，故由物理學(xué)中理學(xué)中Torriceli定律：小孔流出的液體的流速正比于定律：小孔流出的液體的流速正比于水面高度的平方根。問題中給出水塔的最低和最高水水面高度的平方根。問題中給出水塔的最低和最高水位分別是位分別是8.1648m)3024. 027(和和10.7352m)3024. 050.35(（設(shè)出口的水位為零），因為（設(shè)出口的水位為零），因為sqrt1467

34、. 1)1648. 87352.10(，約為，約為1，所以可忽略水位對流速的影響。，所以可忽略水位對流速的影響。（2）將流量看作是時間的連續(xù)函數(shù)，為計算簡單，）將流量看作是時間的連續(xù)函數(shù)，為計算簡單，不妨將不妨將流量定義成單位時間流出水的高度流量定義成單位時間流出水的高度，即水位對，即水位對時間變化率的絕對值（水位是下降的），水塔截面積時間變化率的絕對值（水位是下降的），水塔截面積為為3475.2334*)3024. 0*57(2S（m2），得到結(jié)果），得到結(jié)果后乘以后乘以s即可得真實的流量。即可得真實的流量。 2流量估計方法流量估計方法 首先依照表首先依照表6-2所給數(shù)據(jù)，用所給數(shù)據(jù)，用MA

35、TLAB作出時間作出時間水水位散點圖（圖位散點圖（圖6.2）。）。下面來計算水箱下面來計算水箱流量流量與與時間時間的關(guān)系。的關(guān)系。 根據(jù)圖根據(jù)圖6.2，一種簡單的處理方法為，將表，一種簡單的處理方法為，將表6-2中的中的數(shù)據(jù)分為三段，然后對每一段的數(shù)據(jù)做如下處理：數(shù)據(jù)分為三段，然后對每一段的數(shù)據(jù)做如下處理： 設(shè)某段數(shù)據(jù)設(shè)某段數(shù)據(jù)),(,),(),(1100nnyxyxyx，相鄰數(shù)據(jù)中點的平均流速用下面的公式（流速，相鄰數(shù)據(jù)中點的平均流速用下面的公式（流速=（左（左端點的水位右端點的水位）端點的水位右端點的水位）/區(qū)間長度）：區(qū)間長度）：iiiiiixxyyxxv111)2(每段數(shù)據(jù)首尾點的流速

36、用下面的公式計算：每段數(shù)據(jù)首尾點的流速用下面的公式計算：)()43()(022100 xxyyyxv)()43()(221nnnnnnxxyyyxv用以上公式求得時間與流速之間的數(shù)據(jù)如表用以上公式求得時間與流速之間的數(shù)據(jù)如表6-3。時間/h流速/cmh-1時間/h流速/ cmh-1029.8912.4931.520.4621.7413.4229.031.3818.4814.4326.362.39516.2215.4426.093.4116.3016.3724.734.42515.3217.3823.645.4413.0418.4923.426.4515.4519.5025.007.46513.

37、9820.4023.868.4516.3520.8422.178.9719.2922.02水泵開動9.98水泵開動22.96水泵開動10.93水泵開動23.8827.0910.9533.5024.4321.6211.4929.6325.4518.4825.9113.30由表由表6-3作出時間作出時間流速散點圖如圖流速散點圖如圖6.3。1）插值法）插值法 由表由表6-3，對水泵不工作時段，對水泵不工作時段1，2采取插值方法，可以得采取插值方法，可以得到任意時刻的流速，從而可以知道任意時刻的流量。到任意時刻的流速，從而可以知道任意時刻的流量。 我們分別采取拉格朗日插值法，分段線性插值法及我們分別采

38、取拉格朗日插值法，分段線性插值法及三次樣條插值法；對于水泵工作時段三次樣條插值法；對于水泵工作時段1應(yīng)用前后時期應(yīng)用前后時期的流速進(jìn)行插值，由于最后一段水泵不工作時段數(shù)據(jù)的流速進(jìn)行插值，由于最后一段水泵不工作時段數(shù)據(jù)太少，我們將它與水泵工作時段太少，我們將它與水泵工作時段2合并一同進(jìn)行插值合并一同進(jìn)行插值處理（該段簡稱混合時段）。處理（該段簡稱混合時段）。 我們總共需要對四段數(shù)據(jù)（第我們總共需要對四段數(shù)據(jù)（第1，2未供水時段，第未供水時段，第1供水時段，混合時段）進(jìn)行插值處理，下面以第供水時段，混合時段）進(jìn)行插值處理，下面以第1未供未供水時段數(shù)據(jù)為例分別用三種方法算出水時段數(shù)據(jù)為例分別用三種方

39、法算出流量函數(shù)流量函數(shù)和和用水用水量量（用水高度）。（用水高度）。 下面是用下面是用MATLAB實現(xiàn)該過程的程序?qū)崿F(xiàn)該過程的程序(第第1未供水時段未供水時段)。 t=0,0.46,1.38,2.395,3.41,4.425,12.44,6.45,7.465,8.45,8.97；v=29.89,21.74,18.48,16.22,16.30,15.32,13.04,15.45,13.98,16.35,19.29； t0=0:0.1:8.97；lglr=lagrange(t,v,t0); /*注：注：lagrange為一函數(shù)，程序同為一函數(shù)，程序同 lagrange.m*/lglrjf=trapz

40、(t0,lglr)fdxx=interpl(t,v,t0)；fdxxjf=trapz(t0,fdxx)scyt=interpl(t,v,t0,spline)；sancytjf=trapz(t0,scyt)plot(t,v,o,t0,lglr,r,t0,fdxx,g,t0,scyt,b)gtext(lglr)gtext(fdxx)gtext(spline)運行結(jié)果lglrjf=145.6231 ； fdxxjf=147.1430； sancytjf=145.6870 圖圖6.4是對第是對第1未供水段數(shù)據(jù)用三種不同方法得到的未供水段數(shù)據(jù)用三種不同方法得到的插值函數(shù)圖，插值函數(shù)圖， 圖中曲線圖中曲線

41、lglr、fdxx和和spline分別表示用拉格朗日插分別表示用拉格朗日插值法，分段線性插值法及三次樣條插值法得到的曲線。值法，分段線性插值法及三次樣條插值法得到的曲線。 由表由表6-2知，第知，第1未供水時段的總用水高度為未供水時段的總用水高度為146(=968822)，可見上述三種插值方法計算的結(jié)果，可見上述三種插值方法計算的結(jié)果與實際值（與實際值（146）相比都比較接近?？紤]到三次樣條插）相比都比較接近。考慮到三次樣條插值方法具有更加良好的性質(zhì)，建議采取該方法。值方法具有更加良好的性質(zhì)，建議采取該方法。 其他三段的處理方法與第其他三段的處理方法與第1未供水時段的處理方法未供水時段的處理方

42、法類似，這里不再詳細(xì)敘述，只給出數(shù)值結(jié)果和函數(shù)圖類似，這里不再詳細(xì)敘述，只給出數(shù)值結(jié)果和函數(shù)圖像（圖像（圖6.5圖圖6.7），圖中曲線標(biāo)記同圖），圖中曲線標(biāo)記同圖6.4。 圖圖6.5 第一供水段時間第一供水段時間流速示意圖流速示意圖圖圖12.6 第第2未供水段時間未供水段時間流速示意圖流速示意圖圖圖6.7 混合時段時間混合時段時間流速示意圖流速示意圖圖圖6.8是用分段線性及三次樣條插值方法得到的整個過是用分段線性及三次樣條插值方法得到的整個過程的時間程的時間流速函數(shù)示意圖。流速函數(shù)示意圖。表表6-4 各時段及一天的總用水量（用水高度）各時段及一天的總用水量（用水高度）第第1未供水段未供水段第第

43、2未供水段未供水段第第3供水段供水段混合時段混合時段全天全天拉格朗日拉格朗日插值法插值法145.6231258.866454.068992.1337550.6921分段線性分段線性插值法插值法147.1430258.969749.605176.4688532.1866三次樣條三次樣條插值法插值法145.6870258.654753.333481.7699539.4450 表表6-5是對一天中任取的是對一天中任取的4個時刻分別用個時刻分別用3種方法得到種方法得到的水塔水流量近似值。的水塔水流量近似值。時間時間6.8810.8815.8822.8815.9826671234851433.74260

44、09085346325.5662241818047734.7099621055169414.8272413793103432.9976262626262725.4465591397849525.471578947368415.0527858196582033.7089553595325925.5490892055742329.41733175863551注：注：拉格朗日插值法拉格朗日插值法分段線性插值法分段線性插值法三次樣條插值法三次樣條插值法2）擬合法）擬合法（1）擬合水位）擬合水位時間函數(shù)時間函數(shù) 從表從表6-2中的測量記錄看，一天有兩次供水時段和三中的測量記錄看，一天有兩次供水時段和三次

45、未供水時段，分別對第次未供水時段，分別對第1，2未供水時段的測量數(shù)據(jù)未供水時段的測量數(shù)據(jù)直接作多項式擬合，可得到水位函數(shù)（注意，根據(jù)多直接作多項式擬合，可得到水位函數(shù)（注意，根據(jù)多項式擬合的特點，此處擬合多項式的次數(shù)不宜過高，項式擬合的特點，此處擬合多項式的次數(shù)不宜過高，一般以一般以36次為宜）。對第次為宜）。對第3未供水時段來說，數(shù)據(jù)未供水時段來說，數(shù)據(jù)過少不能得到很好的擬合。過少不能得到很好的擬合。 設(shè)設(shè)t，h分別為已輸入的時刻和水位測量記錄（由表分別為已輸入的時刻和水位測量記錄（由表6.2提供，水泵啟動的提供，水泵啟動的4個時刻不輸入），這樣第個時刻不輸入），這樣第1未未供水時段各時刻的

46、水位可由如下供水時段各時刻的水位可由如下MATLAB程序完成：程序完成：t=0,0.92,1.84,2.95,3.87,4.98,5.90,7.00,7.93,8.97,10.95,12.03,12.95,13.88,14.98,15.90,16.83,17.93,19.04,19.96,20.84,23.88, 24.99,25.91h=9.68,9.45,9.31,9.13,8.98,8.81,8.69,8.52,8.39,8.22,10.82,10.50,10.21,9.94,9.65,9.41,9.18,8.92,8.66,8.43,8.22,10.59,10.35,10.18；c1=

47、polyfit(t(1:10),h(1:10),3)；% 用用3次多項式擬合次多項式擬合tp1=0:0.1:8.9；x1=polyval(c1,tp1)；plot(tp1,x1) 圖圖6.9給出的是第給出的是第1未供水時段的時間未供水時段的時間水位擬合水位擬合函數(shù)圖形。函數(shù)圖形。變量變量x1中存放了以中存放了以0.1為步長算出的各個時刻的水位高為步長算出的各個時刻的水位高度。同樣地，第度。同樣地，第2未供水時段時間未供水時段時間水位圖可由如下水位圖可由如下MATLAB程序完成，讀者可自己上機運行查看。程序完成，讀者可自己上機運行查看。c2=polyfit(t(11:21),h(11:21),3

48、)；tp2=10.9:0.1:20.9；x2=-polyval(c2,tp2)；plot(tp2,x2)（2）確定流量）確定流量時間函數(shù)時間函數(shù)對于第對于第1，2未供水時段的流量可直接對水位函數(shù)求導(dǎo)未供水時段的流量可直接對水位函數(shù)求導(dǎo)，程序如下：，程序如下：c1=polyfit(t(1:10),h(1:10),3)；c2=polyfit(t(11:21),h(11:21),3)；a1=polyder(c1)；a2=polyder(c2)； % 對擬合得到的水位函數(shù)求導(dǎo)，得流量對擬合得到的水位函數(shù)求導(dǎo)，得流量tp1=0:0.01:8.97；tp2=10.95:0.01:20.84；x13=-po

49、lyval(a1,tp1)； % Water_total1=100*trapz(tp1,x13); % 計算第計算第1未供水時段的總用水量未供水時段的總用水量x14=-polyval(a1,7.93,8.97)； % 求某些點的流量求某些點的流量 x23=-polyval(a2,tp2)；Water_total2=100*trapz(tp2,x114); %計算第計算第2未供水時段的總用水量未供水時段的總用水量x24=-polyval(a2,10.95,12.03)； x25=-polyval(a2,19.96,20.84)；subplot(1,2,1)plot(tp1,x13*100) %與

50、圖與圖6.10單位保持一致單位保持一致 % 倍數(shù)倍數(shù)100是使單位由米變成厘米是使單位由米變成厘米subplot(1,2,2)plot(tp2,x23*100) %與圖與圖6.10單位保持一致單位保持一致 程序運行得到第程序運行得到第1，2未供水時段的時間未供水時段的時間流量圖如圖流量圖如圖6.10，可以看到與圖，可以看到與圖6.8中用插值給出的曲線比較吻合。中用插值給出的曲線比較吻合。 如果用如果用5次多項式擬合則得圖次多項式擬合則得圖6.11的圖形，顯然較的圖形，顯然較三次擬合的效果好。三次擬合的效果好。 而第而第1供水時段的流量則用前后時期的流量進(jìn)行擬供水時段的流量則用前后時期的流量進(jìn)行

51、擬合得到。為使流量函數(shù)在合得到。為使流量函數(shù)在t=9和和t=11連續(xù)，我們只取連續(xù)，我們只取4個點，用三次多項式擬合得到第個點，用三次多項式擬合得到第1供水時段的時間供水時段的時間流量圖形如圖流量圖形如圖6.12，可以看到與圖，可以看到與圖6.8中的相應(yīng)部分中的相應(yīng)部分比較吻合。比較吻合。圖圖6.12圖圖6.8程序如下：程序如下： 在第在第2供水時段之前取供水時段之前取t=19.96,20.84兩點的流量，兩點的流量，用第用第3未供水時段的未供水時段的3個記錄做差分得到兩個流量數(shù)個記錄做差分得到兩個流量數(shù)據(jù)據(jù)21.62，18.48，然后用這，然后用這4個數(shù)據(jù)做三次多項式擬個數(shù)據(jù)做三次多項式擬合

52、得到第合得到第2供水時段與第供水時段與第3未供水時段的時間未供水時段的時間流量流量圖如圖圖如圖6.13，可以看到與圖，可以看到與圖6.8中的相應(yīng)部分也比較中的相應(yīng)部分也比較吻合吻合。圖圖6.13，圖圖6.8程序如下：程序如下：t3=19.96,20.84,t(22),t(23)；x3=x25*100,21.62,18.48；x_flow3=polyfit(t3,x3,3)；tp3=19.96:0.01:25.91；xx3=polyval(x_flow3,tp3)；tp4=20.84:0.01:23.88；gs2=polyval(x_flow3,tp4)；gsysl2=trapz(tp4,gs2); %該語句計算第該語句計算第2供水時段的總用水量供水時段的總用水量figure;plot(tp3,xx3)（3）一天總用水量的估計）一天總用水量的估計 分別對供水的兩個時段和不供水的兩個時段積分分別對供水的兩個時段和不供水的兩個時段積分（流量對時間）并求和得到一天的總用水量約為（流量對時間）并求和得到一天的總用水量約為526.8935（此數(shù)據(jù)是總用水高度，單位為（此數(shù)據(jù)是總用水高度，單位為cm）。表）。表6-6列出各段用水量，與插值法算得的表列出各段用水量，與插值法算得的表6-4相比，二相比，二者較為吻合。者較為吻合。表表6-6時段時段第第