【翻譯搬運】SciPy-Python科學演算法庫

01-24

SciPy庫提供了大量有用的函數和類，用來解決各種專業領域的問題。本文翻譯自Jupyter nbviewer中的第三講。首先，介紹了一些特殊函數，如貝塞爾函數，這對物理學問題的計算提供了方便；之後是各種數值積分問題，常微分方程求解問題以及傅里葉變換，這些也可以描述並求解一些諸如復擺運動、阻尼震蕩等複雜的物理過程；同時，該庫還可以高效處理線性代數問題，如矩陣的運算、特徵值和特徵向量以及稀疏矩陣的存儲和運算；最優化問題，即尋找函數極值和零點的問題，和插值問題，即用多項式擬合曲線的問題，在SciPy庫中也可以找到相應的函數；最後介紹了統計上的應用，包括各種分布的密度函數、分布函數及其圖像繪製，以及一些統計檢驗問題。

作者：J.R. Johansson (郵箱：jrjohansson@gmail.com)

譯者：一路向上

最新版本的用法介紹見網站http://github.com/jrjohansson/scientific-python-lectures.其他相關介紹見http://jrjohansson.github.io.

簡介

SciPy框架建立於低一級的NumPy框架的多維數組之上，並且提供了大量的高級的科學演算法。一些SciPy的應用如下：

特殊函數 (scipy.special)
積分 (scipy.integrate)
優化 (scipy.optimize)
插值 (scipy.interpolate)
傅里葉變換 (scipy.fftpack)
信號處理 (scipy.signal)
線型代數 (scipy.linalg)
稀疏特徵值問題 (scipy.sparse)
統計 (scipy.stats)
多維圖像處理 (scipy.ndimage)
文件輸入輸出 (http://scipy.io)

這些模塊提供了大量的函數和類，可以用來解決各自領域的問題。在本節中我們將看到如何使用這些子函數包。我們首先導入scipy程序包：

如果我們只需要用scipy框架中的一部分，我們也可以選擇性的導入。例如，只導入線性代數函數包la，我們可以：

特殊函數

大量的數學函數對許多物理問題的計算是非常重要的。SciPy提供了一系列非常廣泛的特殊函數。詳見http://docs.scipy.org/doc/scipy/reference/special.html#module-scipy.special.

為了說明指定的特殊函數的用法，我們先看一下貝塞爾函數的使用細節：

積分

數值積分：求積公式

對f(x)從a到b的積分叫做數值積分，也叫簡單積分。SciPy提供了一系列計算不同數值積分的函數，包括quad,dblquad和tplquad，分別包含二重和三重積分。

quad有一系列的可供選擇的參數，可以用來調節函數的各種行為（輸入help(quad)獲取更多信息）。

其基本用途如下：

如果我們需要添加更多對於被積函數的參數，可以使用args關鍵字參數：

對於簡單函數而言，對於被積函數，我們可以用λ函數（無名稱的函數）來代替清晰定義的函數：

如上所示，我們可以用"Inf"或者"-Inf"作為積分上下限。高維積分用法相同：

注意，我們需要用λ函數對於y積分的極限，因為它們可以看做是x的函數。

常微分方程（ODE）

SciPy提供了兩種不同的方法來解決常微分方程：函數odeint的API，和函數類ode面向對象的API。通常odeint比較容易上手，但是ode函數類能夠更好的控制函數。

這裡我們使用odeint函數，如需了解更多ode函數類的信息，請輸入help(ode)。它和odeint很像，但卻是面向對象的函數。

使用odeint之前，首先從scipy.integrate中調用它：

常微分方程系通常寫作其一般形式：

y" = f(y, t)

其中

y = [y1(t), y2(t), ..., yn(t)]

f的微分是yi(t)。為了解決常微分方程，我們需要知道函數f和初始條件y(0).

高階微分方程可以通過引進新的變數作為中間變數。

當我們定義了Python函數f和數組y_0(f和y(0)都是數學函數），我們調用odeint函數：

y_t = odeint(f, y_0, t)

t是解決ODE問題需要的時間坐標數組，y_t是對於給定點在時間t的一行數組，每一列代表在給定時間t所對應的一個解y_i(t)。我們下面將會看到如何設置f和y_0.

例：復擺

我們考慮一個物理問題：復擺。定義詳見http://en.wikipedia.org/wiki/Double_pendulum.

為了讓Python代碼看起來更簡潔，我們引進新的變數，並規定：

在matplotlib函數的應用中，會介紹如何繪製復擺運動的動圖。

ps:這裡的結果不是報錯，是因為最後一行代碼是每0.1s更新一次狀態，為了後面的函數能夠正常運行，我把它停掉了。

例：阻尼諧波振蕩器

常微分方程問題對計算物理非常重要，下面我們來看另一個例子：阻尼諧波振蕩器。關於概念的解釋詳見 http://en.wikipedia.org/wiki/Damping.

阻尼諧波振蕩器的方程為：

x是振蕩器的位置， $omega_{0}$ 是頻率， $zeta$ 是阻尼比。為了寫出標準形式的二階常微分方程，我們引入 $p=frac{dx}{dt}$ ：

在這個例子中，我們將為常微分方程等號右邊的函數添加額外的參數，而不是像前面的例子那樣使用全局變數。作為等號右邊函數的額外參數的結果，我們需要將一個關鍵字參數args傳遞給odeint函數：

傅里葉變換

傅里葉變換是計算物理中的一種通用工具，它在不同文章中都會反覆出現。SciPy提供能夠從NetLib中接入經典FFTPACK庫的函數，它是由FORTRAN語言編寫的一個行之有效的FFT庫。SciPy API有一些額外的便利功能，但總的來說，API和原來的FORTRAN庫密切相關。

為了調用fftpack，請輸入：

為演示如何用SciPy做一個快速傅里葉變換，讓我們來看看用FFT如何解決之前討論的阻尼震蕩問題：

由於信號是實數，所以譜線圖是對稱的。我們因此只需要畫出對應正頻率的部分。為了提取w和F的部分，我們可以運用NumPy庫：

和預期一樣，我們看到譜線圖在1處達到最高點，這正是在阻尼震蕩這個例子中所採用的頻率。

線性代數

線性代數部分包含了許多矩陣函數，包括線性方程的解，特徵值的解，矩陣函數（例如矩陣指數函數），許多不同的分解（SVD，LU，cholesky）等等。詳見：http://docs.scipy.org/doc/scipy/reference/linalg.html.

下面我們看看如何使用這些函數：

1、線性方程組

線性方程組的矩陣形式：

A是矩陣，x,b是向量。可以求解如下：

我們也可以：

這裡A,B,X都是矩陣：

2、特徵值和特徵向量

矩陣A的特徵值問題是：

這裡 $v_{n}$ 是第n個特徵向量， $lambda_{n}$ 是第n個特徵值：

用eigvals計算矩陣的特徵值，用eig計算特徵值和特徵向量：

第n個特徵值（儲存在evals[n]中）對應的特徵向量是evecs的第n列，也就是evecs[:,n]. 為了驗證它，我們嘗試把矩陣和特徵向量相乘，並與特徵向量和特徵值的乘積做比較：

還有許多其他的特殊的本徵解，如埃爾米特矩陣（用eigh實現）

3、矩陣運算

4、稀疏矩陣

稀疏矩陣在處理巨型系統的數值模擬時非常有用，如果描述該問題的矩陣或向量大部分的元素為0。Scipy對於稀疏矩陣有很多處理方法，包括基礎的線性代數處理（包括解方程，計算特徵值等等）

許多方法都能有效存儲稀疏矩陣，一些常用的方法包括坐標形式（COO）,列表的列表形式（LIL），壓縮稀疏列（CSC）和壓縮稀疏行（CSR）。每種方法都有優勢和不足。大多數的計算演算法（解方程，矩陣和矩陣相乘等等）都能用CSR或者CSC形式處理，但是它們並不直觀，也不太容易進行初始化。所以通常來說，稀疏矩陣採用COO或者LIL進行初始化（我們可以在稀疏矩陣中有效添加元素），然後再轉換為CSC或者CSR並進行計算。

更多關於稀疏矩陣的信息，詳見： http://en.wikipedia.org/wiki/Sparse_matrix.

當我們創建了一個稀疏矩陣，我們要選擇其存儲形式，如：

創建稀疏矩陣更有效的方法是，建立一個空矩陣，並用所在矩陣的位置填充（避免創建大的稠密矩陣）：

最優化

最優化問題（尋找函數的最大值或者最小值）是數學的一大領域，複雜函數的最優化問題，或者多變數的最優化問題可能會非常複雜。下面我們看一些很簡單的例子，更多詳細的對於使用SciPy處理最優化問題的介紹，請見：http://scipy-lectures.github.com/advanced/mathematical_optimization/index.html.

首先調用optimize：