Python 語言在數值計算中的應用

羅堂耘,陳卓鑫

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

0 引言

數值計算是強大的問題求解方法,它能處理大規模方程組,尤其是非線性系統的問題。 在工程中,用解析方法對其求解幾乎是不可能的,此時計算數值解就顯得尤為重要[1]。

在數值計算領域,MATLAB 是最常用的工具,而Python 語言多用于人工智能、機器學習等方面。Python 入門簡單,有易學、開源等特點,更重要的是MATLAB 中大部分常用的功能都可以在Python 中找到對應的擴展庫,故而Python 在數值計算領域中也頗有建樹[2]。

本文利用Python 語言,對一道現實的數值計算問題作解答研究。

1 問題提出

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

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

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

2 Python 擴展庫

2.1 Numpy 庫

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

2.2 Scipy 庫

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

2.3 Matplotlib 庫

Matplotlib 是Python 中最著名的繪圖庫,可以繪制各種精美的圖表、快速生成柱狀圖、散點圖、餅狀圖、折線圖、等值線圖和三角網格等[5],配合Numpy模塊使用,可以實現科學計算結果的可視化顯示。

3 模型的建立與求解

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

3.1 數值微分

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

利用多點數值微分公式即可算出流速,具體如下。

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

n1=9

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

if i<=1:進行計算;

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])))

第二段與第三段按類似的操作進行即可。

這樣便得到了24 個點,每個點在x 軸上的值對應時間,y 軸上的值對應流速。

3.2 三次樣條插值

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

圖1 流速-時間圖像

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

(1) 輸入節點個數n,節點x0,x1,x2…xn,對應函數值y0,y1,y2…yn。

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

hi= xi- xi-1

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

μi= 1 - λi

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

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

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

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

tck=interpolate.splrep(x,f)

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

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

圖2 三次樣條插值結果

3.3 數值積分

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

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

3.3.1 整體積分

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

復合梯形公式是常用的數值積分方法,為:

利用如下的Python 代碼進行計算。

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)

最終得到測量時間內水的總流量為49 155.641 6 ft3,進而得到一天時間內水的總流量為44 586.333 0 ft3。

3.3.2 分段積分

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

計算得到:測量時間內水的總流量為 48 121.033 9 ft3,進而一天內水的總流量為43 647.898 3 ft3。

4 結語

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