SciPy - 科學計算庫（上）

一、實驗說明

SciPy 庫建立在 Numpy 庫之上，提供了大量科學算法，主要包括這些主題：

• 特殊函數 (scipy.special)
• 積分 (scipy.integrate)
• 最優化 (scipy.optimize)
• 插值 (scipy.interpolate)
• 傅立葉變換 (scipy.fftpack)
• 信號處理 (scipy.signal)
• 線性代數 (scipy.linalg)
• 稀疏特征值 (scipy.sparse)
• 統計 (scipy.stats)
• 多維圖像處理 (scipy.ndimage)
• 文件 IO (scipy.io)

在本實驗中我們將了解其中一些包的使用方法。

（ps：因本節只講工具的用法，對這些科學主題不展開討論，所以根據自己所學的知識挑選食用就好了，強迫症不要糾結哈～）

1. 環境登錄

無需密碼自動登錄，系統用戶名shiyanlou

2. 環境介紹

本課程實驗環境使用Spyder。首先打開terminal，然后輸入以下命令：

spyder -w scientific-python-lectures 

關於Spyder的使用可參考文檔：https://pythonhosted.org/spyder/

本實驗基本在控制台下進行，可關閉其余窗口，只保留控制台。如需要調出窗口，可以通過 view->windows and toolbar 調出。比如希望在py文件中編寫代碼，可以 view->windows and toolbar->Editor 調出編輯器窗口。

二、實驗內容

讓我們先導入必要的庫

from numpy import * from scipy import * 

特定函數

在計算科學問題時，常常會用到很多特定的函數，SciPy 提供了一個非常廣泛的特定函數集合。函數列表可參考：http://docs.scipy.org/doc/scipy/reference/special.html#module-scipy.special

為了演示特定函數的一般用法我們拿貝塞爾函數舉例：

# # The scipy.special module includes a large number of Bessel-functions # Here we will use the functions jn and yn, which are the Bessel functions # of the first and second kind and real-valued order. We also include the # function jn_zeros and yn_zeros that gives the zeroes of the functions jn # and yn. # %matplotlib qt from scipy.special import jn, yn, jn_zeros, yn_zeros import matplotlib.pyplot as plt n = 0 # order x = 0.0 # Bessel function of first kind print "J_%d(%f) = %f" % (n, x, jn(n, x)) x = 1.0 # Bessel function of second kind print "Y_%d(%f) = %f" % (n, x, yn(n, x)) => J_0(0.000000) = 1.000000 Y_0(1.000000) = 0.088257 x = linspace(0, 10, 100) fig, ax = plt.subplots() for n in range(4): ax.plot(x, jn(n, x), label=r"$J_%d(x)$" % n) ax.legend(); fig 

# zeros of Bessel functions n = 0 # order m = 4 # number of roots to compute jn_zeros(n, m) => array([ 2.40482556, 5.52007811, 8.65372791, 11.79153444]) 

積分

數值積分: 求積

 被稱作 數值求積，Scipy提供了一些列不同類型的求積函數，像是 quad, dblquad 還有 tplquad 分別對應單積分，雙重積分，三重積分。

from scipy.integrate import quad, dblquad, tplquad 

quad 函數有許多參數選項來調整該函數的行為（詳情見help(quad)）。

一般用法如下:

# define a simple function for the integrand def f(x): return x x_lower = 0 # the lower limit of x x_upper = 1 # the upper limit of x val, abserr = quad(f, x_lower, x_upper) print "integral value =", val, ", absolute error =", abserr => integral value = 0.5 , absolute error = 5.55111512313e-15 

如果我們需要傳遞額外的參數，可以使用 args 關鍵字：

def integrand(x, n): """ Bessel function of first kind and order n. """ return jn(n, x) x_lower = 0 # the lower limit of x x_upper = 10 # the upper limit of x val, abserr = quad(integrand, x_lower, x_upper, args=(3,)) print val, abserr => 0.736675137081 9.38925687719e-13 

對於簡單的函數我們可以直接使用匿名函數：

val, abserr = quad(lambda x: exp(-x ** 2), -Inf, Inf) print "numerical =", val, abserr analytical = sqrt(pi) print "analytical =", analytical => numerical = 1.77245385091 1.42026367809e-08 analytical = 1.77245385091 

如例子所示，'Inf' 與 '-Inf' 可以表示數值極限。

高階積分用法類似:

def integrand(x, y): return exp(-x**2-y**2) x_lower = 0 x_upper = 10 y_lower = 0 y_upper = 10 val, abserr = dblquad(integrand, x_lower, x_upper, lambda x : y_lower, lambda x: y_upper) print val, abserr => 0.785398163397 1.63822994214e-13 

注意到我們為y積分的邊界傳參的方式，這樣寫是因為y可能是關於x的函數。

常微分方程 (ODEs)

SciPy 提供了兩種方式來求解常微分方程：基於函數 odeint 的API與基於 ode 類的面相對象的API。通常 odeint 更好上手一些，而 ode 類更靈活一些。

這里我們將使用 odeint 函數，首先讓我們載入它：

from scipy.integrate import odeint, ode 

常微分方程組的標准形式如下:

當

為了求解常微分方程我們需要知道方程  與初始條件 注意到高階常微分方程常常寫成引入新的變量作為中間導數的形式。 一旦我們定義了函數 f 與數組y_0 我們可以使用 odeint 函數：

y_t = odeint(f, y_0, t) 

我們將會在下面的例子中看到 Python 代碼是如何實現 f 與 y_0 。

示例: 雙擺

讓我們思考一個物理學上的例子：雙擺

關於雙擺，參考：http://en.wikipedia.org/wiki/Double_pendulum

Image(url='http://upload.wikimedia.org/wikipedia/commons/c/c9/Double-compound-pendulum-dimensioned.svg') 

維基上已給出雙擺的運動方程:

為了使 Python 代碼更容易實現，讓我們介紹新的變量名與向量表示法： 

g = 9.82 L = 0.5 m = 0.1 def dx(x, t): """ The right-hand side of the pendulum ODE """ x1, x2, x3, x4 = x[0], x[1], x[2], x[3] dx1 = 6.0/(m*L**2) * (2 * x3 - 3 * cos(x1-x2) * x4)/(16 - 9 * cos(x1-x2)**2) dx2 = 6.0/(m*L**2) * (8 * x4 - 3 * cos(x1-x2) * x3)/(16 - 9 * cos(x1-x2)**2) dx3 = -0.5 * m * L**2 * ( dx1 * dx2 * sin(x1-x2) + 3 * (g/L) * sin(x1)) dx4 = -0.5 * m * L**2 * (-dx1 * dx2 * sin(x1-x2) + (g/L) * sin(x2)) return [dx1, dx2, dx3, dx4] # choose an initial state x0 = [pi/4, pi/2, 0, 0] # time coodinate to solve the ODE for: from 0 to 10 seconds t = linspace(0, 10, 250) # solve the ODE problem x = odeint(dx, x0, t) # plot the angles as a function of time fig, axes = plt.subplots(1,2, figsize=(12,4)) axes[0].plot(t, x[:, 0], 'r', label="theta1") axes[0].plot(t, x[:, 1], 'b', label="theta2") x1 = + L * sin(x[:, 0]) y1 = - L * cos(x[:, 0]) x2 = x1 + L * sin(x[:, 1]) y2 = y1 - L * cos(x[:, 1]) axes[1].plot(x1, y1, 'r', label="pendulum1") axes[1].plot(x2, y2, 'b', label="pendulum2") axes[1].set_ylim([-1, 0]) axes[1].set_xlim([1, -1]); fig 

我們將在第四節課看到如何做出更好的演示動畫。

from IPython.display import clear_output
import time fig, ax = plt.subplots(figsize=(4,4)) for t_idx, tt in enumerate(t[:200]): x1 = + L * sin(x[t_idx, 0]) y1 = - L * cos(x[t_idx, 0]) x2 = x1 + L * sin(x[t_idx, 1]) y2 = y1 - L * cos(x[t_idx, 1]) ax.cla() ax.plot([0, x1], [0, y1], 'r.-') ax.plot([x1, x2], [y1, y2], 'b.-') ax.set_ylim([-1.5, 0.5]) ax.set_xlim([1, -1]) display(fig) clear_output() time.sleep(0.1) fig 

示例：阻尼諧震子

常微分方程問題在計算物理學中非常重要，所以我們接下來要看另一個例子：阻尼諧震子。wiki地址：http://en.wikipedia.org/wiki/Damping

阻尼震子的運動公式：

其中  是震子的位置,  是頻率,  是阻尼系數. 為了寫二階標准行事的 ODE 我們引入變量：

:

在這個例子的實現中，我們會加上額外的參數到 RHS 方程中：

def dy(y, t, zeta, w0): """ The right-hand side of the damped oscillator ODE """ x, p = y[0], y[1] dx = p dp = -2 * zeta * w0 * p - w0**2 * x return [dx, dp] # initial state: y0 = [1.0, 0.0] # time coodinate to solve the ODE for t = linspace(0, 10, 1000) w0 = 2*pi*1.0 # solve the ODE problem for three different values of the damping ratio y1 = odeint(dy, y0, t, args=(0.0, w0)) # undamped y2 = odeint(dy, y0, t, args=(0.2, w0)) # under damped y3 = odeint(dy, y0, t, args=(1.0, w0)) # critial damping y4 = odeint(dy, y0, t, args=(5.0, w0)) # over damped fig, ax = plt.subplots() ax.plot(t, y1[:,0], 'k', label="undamped", linewidth=0.25) ax.plot(t, y2[:,0], 'r', label="under damped") ax.plot(t, y3[:,0], 'b', label=r"critical damping") ax.plot(t, y4[:,0], 'g', label="over damped") ax.legend(); fig 

傅立葉變換

傅立葉變換是計算物理學所用到的通用工具之一。Scipy 提供了使用 NetLib FFTPACK 庫的接口，它是用FORTRAN寫的。Scipy 還另外提供了很多便捷的函數。不過大致上接口都與 NetLib 的接口差不多。

讓我們加載它：

from scipy.fftpack import * 

下面演示快速傅立葉變換，例子使用上節阻尼諧震子的例子：

N = len(t)
dt = t[1]-t[0] # calculate the fast fourier transform # y2 is the solution to the under-damped oscillator from the previous section F = fft(y2[:,0]) # calculate the frequencies for the components in F w = fftfreq(N, dt) fig, ax = plt.subplots(figsize=(9,3)) ax.plot(w, abs(F)); fig 

既然信號是實數，同時頻譜是對稱的。那么我們只需要畫出正頻率所對應部分的圖：

indices = where(w > 0) # select only indices for elements that corresponds to positive frequencies
w_pos = w[indices]
F_pos = F[indices]


fig, ax = subplots(figsize=(9,3))
ax.plot(w_pos, abs(F_pos))
ax.set_xlim(0, 5);

fig

正如預期的那樣，我們可以看到頻譜的峰值在1處。1就是我們在上節例子中所選的頻率。