羅堂耘,陳卓鑫

(成都理工大學(xué),四川 成都 610059)

0 引言

數(shù)值計(jì)算是強(qiáng)大的問題求解方法,它能處理大規(guī)模方程組,尤其是非線性系統(tǒng)的問題。 在工程中,用解析方法對其求解幾乎是不可能的,此時(shí)計(jì)算數(shù)值解就顯得尤為重要[1]。

在數(shù)值計(jì)算領(lǐng)域,MATLAB 是最常用的工具,而Python 語言多用于人工智能、機(jī)器學(xué)習(xí)等方面。Python 入門簡單,有易學(xué)、開源等特點(diǎn),更重要的是MATLAB 中大部分常用的功能都可以在Python 中找到對應(yīng)的擴(kuò)展庫,故而Python 在數(shù)值計(jì)算領(lǐng)域中也頗有建樹[2]。

本文利用Python 語言,對一道現(xiàn)實(shí)的數(shù)值計(jì)算問題作解答研究。

1 問題提出

某地有一個(gè)高40 ft、底面直徑57 ft 的圓柱形水塔,每天向小鎮(zhèn)上的居民供應(yīng)生活用水。 每隔一段時(shí)間測量一次水塔中的水位,測得數(shù)據(jù)如表1 所示(兩次水泵啟動向水塔中加水的過程沒有水位記錄)。 試估計(jì)一天的總用水量。

表1 某地某天水塔水位測量記錄

求解該問題,首先需要建立水位對于時(shí)間的函數(shù),隨后根據(jù)離散點(diǎn)的公式得到各時(shí)間點(diǎn)的流速,再利用3 次樣條插值得到流速與時(shí)間的關(guān)系,對流速函數(shù)作數(shù)值積分,進(jìn)而估計(jì)出問題給定時(shí)間內(nèi)的總用水量,最后得到一天的總用水量。

2 Python 擴(kuò)展庫

2.1 Numpy 庫

Numpy 是Python 中的科學(xué)計(jì)算庫,是Python 可以高效處理數(shù)組并進(jìn)行矢量運(yùn)算的保障。 Numpy 庫使得運(yùn)算效率比用戶手動作數(shù)組運(yùn)算要高許多[3]。

2.2 Scipy 庫

Scipy 在NumPy 的基礎(chǔ)上增加了許多科學(xué)計(jì)算、工程計(jì)算的模塊,例如線性代數(shù)、微分方程和曲線擬合等[4],利用這些功能可以更高效地作數(shù)值分析。

2.3 Matplotlib 庫

Matplotlib 是Python 中最著名的繪圖庫,可以繪制各種精美的圖表、快速生成柱狀圖、散點(diǎn)圖、餅狀圖、折線圖、等值線圖和三角網(wǎng)格等[5],配合Numpy模塊使用,可以實(shí)現(xiàn)科學(xué)計(jì)算結(jié)果的可視化顯示。

3 模型的建立與求解

為了簡化模型的復(fù)雜程度,但凡是本文未提及的細(xì)節(jié),一切都視作理想情況。 如假定水泵抽水的速度是均勻的,水流速度是連續(xù)變化的。

3.1 數(shù)值微分

除去水泵啟動的4 個(gè)時(shí)刻點(diǎn),還有24 個(gè)時(shí)刻點(diǎn)可供用于數(shù)值微分。 將水泵工作看作是每兩段時(shí)間的分隔標(biāo)識,即將第0 s 到第32 284 s 這段時(shí)間記為第一段,第39 435 s 到第75 021 s 這段時(shí)間記為第二段,第85 968 s 到第93 270 s 這段時(shí)間記為第三段。

利用多點(diǎn)數(shù)值微分公式即可算出流速,具體如下。

這不需要手動計(jì)算,以第一段為例,可以利用如下的Python 代碼快速得到結(jié)果。

n1=9

for i in range(0,n1+1):

if i<=1:進(jìn)行計(jì)算;

f[i] = abs((-3×v[i] + 4×v[i+1]-v[i+2])/(2×(t[i+1]-t[i])))

elif i<=n1-3:

f[i] = abs((-v[i+2]+8×v[i+1]-8×v[i-1]+v[i-2])/(12×(t[i+1]-t[i])))

elif i>=n1-2:

f[i] = abs((3×v[i]-4×v[i-1]+v[i-2])/(2×(t[i]-t[i-1])))

第二段與第三段按類似的操作進(jìn)行即可。

這樣便得到了24 個(gè)點(diǎn),每個(gè)點(diǎn)在x 軸上的值對應(yīng)時(shí)間,y 軸上的值對應(yīng)流速。

3.2 三次樣條插值

樣條(Spline)一詞來源于繪圖技術(shù),在繪圖中常用一根窄的、柔軟的條子畫出一條經(jīng)過一組已知點(diǎn)的平滑曲線,這根條子就是樣條。 在所有的樣條插值方法中,最常用的是三次樣條插值法,將一個(gè)函數(shù)分為若干小區(qū)間,在每一個(gè)小區(qū)間上都有一個(gè)三次多項(xiàng)式。

圖1 流速-時(shí)間圖像

前面已經(jīng)得到了24 個(gè)點(diǎn)的水流速f(t)的值,故只需根據(jù)這24 個(gè)值作三次樣條插值即可。 三次樣條插值的基本思想如下[6]。

(1) 輸入節(jié)點(diǎn)個(gè)數(shù)n,節(jié)點(diǎn)x0,x1,x2…xn,對應(yīng)函數(shù)值y0,y1,y2…yn。

(2) 計(jì)算步長hi及差商f[xi-1,xi],i = 0,1,2…n。

hi= xi- xi-1

(3) 計(jì)算參數(shù)λi,μi,di,i=0,1,2…n。

μi= 1 - λi

di= 6f[xi-1,xi,xi+1]

(4) 得到各參數(shù)后,可將待解式寫作如下的矩陣形式。

(5) 解上式得到:M0,M1,M2…Mn-1。

通過手動計(jì)算或編程來實(shí)現(xiàn)三次樣條插值十分復(fù)雜煩瑣,可以借助Python 中的Scipy 庫為接下來的數(shù)值計(jì)算研究節(jié)省時(shí)間。 利用如下的Python 代碼即可求解。

tck=interpolate.splrep(x,f)

xx=np.linspace(min(x),max(x),360)

yy=interpolate.splev(xx,tck,der=0)

圖2 三次樣條插值結(jié)果

3.3 數(shù)值積分

正如人們所熟知的牛頓-萊布尼茨公式,若F(x)是f(x)的原函數(shù),則有積分式I=∫baf(x)dx=F(b)-F(a)。

但實(shí)際上這是有困難的,如本問題:計(jì)算流速f(x) 的原函數(shù)F(x) 很困難,甚至是無法實(shí)現(xiàn)的。 數(shù)值積分方法由此被提出。

3.3.1 整體積分

在得到f(x)之后,再對其整體進(jìn)行一次數(shù)值積分——以0 為積分下限,96 270 為積分上限,180 為步長。 這樣便得到了測量時(shí)間內(nèi)的總水流量,進(jìn)而也就可以得到一天(86 400 s)內(nèi)水的總流量。

復(fù)合梯形公式是常用的數(shù)值積分方法,為:

利用如下的Python 代碼進(jìn)行計(jì)算。

tck=interpolate.splrep(x,f)

t2=np.linspace(start = 0, stop = 96 270, num =510)

nn=len(t2)

f2=interpolate.splev(t2,tck,der=0)dt=180

s=(f2[0]+f2[nn-1]+2×sum(f2[np. arange(1,nn-2)]))×dt/(2×3 600)

最終得到測量時(shí)間內(nèi)水的總流量為49 155.641 6 ft3,進(jìn)而得到一天時(shí)間內(nèi)水的總流量為44 586.333 0 ft3。

3.3.2 分段積分

鑒于水塔水泵一共工作了兩次,故可使用數(shù)值積分分別計(jì)算3 段時(shí)間內(nèi)的流量,再將各段流量相加得到總流量。 換言之,就是將原本的一次整體數(shù)值積分改成了3 次分段數(shù)值積分,并將各段結(jié)果相加。 此操作細(xì)節(jié)與整體積分類似,不再贅述。

計(jì)算得到:測量時(shí)間內(nèi)水的總流量為 48 121.033 9 ft3,進(jìn)而一天內(nèi)水的總流量為43 647.898 3 ft3。

4 結(jié)語

本文討論了使用Python 進(jìn)行數(shù)值計(jì)算的優(yōu)勢和應(yīng)用。 Python 語言語法簡單,其可讀性和可維護(hù)性使其成為一種理想的工具,且配備有極其強(qiáng)大的第三方庫,可以滿足各種數(shù)值計(jì)算需求。 本文以一道現(xiàn)實(shí)生活問題為背景,介紹了諸多數(shù)值計(jì)算領(lǐng)域中的經(jīng)典方法,使用Python 語言將這些方法依次實(shí)現(xiàn),并串聯(lián)起來,希望能為讀者在相關(guān)問題上的研究提供一定的參考及幫助。