# SciPy-數值計算庫
SciPy函數庫在NumPy庫的基礎上增加了眾多的數學、科學以及工程計算中常用的庫函數。例如線性代數、常微分方程數值求解、信號處理、圖像處理、稀疏矩陣等等。由于其涉及的領域眾多、本書沒有能力對其一一的進行介紹。作為入門介紹，讓我們看看如何用SciPy進行插值處理、信號濾波以及用C語言加速計算。
## 最小二乘擬合
假設有一組實驗數據(x[i], y[i])，我們知道它們之間的函數關系:y = f(x)，通過這些已知信息，需要確定函數中的一些參數項。例如，如果f是一個線型函數f(x) = k*x+b，那么參數k和b就是我們需要確定的值。如果將這些參數用 **p** 表示的話，那么我們就是要找到一組 **p** 值使得如下公式中的S函數最小：
這種算法被稱之為最小二乘擬合(Least-square fitting)。
scipy中的子函數庫optimize已經提供了實現最小二乘擬合算法的函數leastsq。下面是用leastsq進行數據擬合的一個例子：
```
# -*- coding: utf-8 -*-
import numpy as np
from scipy.optimize import leastsq
import pylab as pl
def func(x, p):
"""
數據擬合所用的函數: A*sin(2*pi*k*x + theta)
"""
A, k, theta = p
return A*np.sin(2*np.pi*k*x+theta)
def residuals(p, y, x):
"""
實驗數據x, y和擬合函數之間的差，p為擬合需要找到的系數
"""
return y - func(x, p)
x = np.linspace(0, -2*np.pi, 100)
A, k, theta = 10, 0.34, np.pi/6 # 真實數據的函數參數
y0 = func(x, [A, k, theta]) # 真實數據
y1 = y0 + 2 * np.random.randn(len(x)) # 加入噪聲之后的實驗數據
p0 = [7, 0.2, 0] # 第一次猜測的函數擬合參數
# 調用leastsq進行數據擬合
# residuals為計算誤差的函數
# p0為擬合參數的初始值
# args為需要擬合的實驗數據
plsq = leastsq(residuals, p0, args=(y1, x))
print u"真實參數:", [A, k, theta]
print u"擬合參數", plsq[0] # 實驗數據擬合后的參數
pl.plot(x, y0, label=u"真實數據")
pl.plot(x, y1, label=u"帶噪聲的實驗數據")
pl.plot(x, func(x, plsq[0]), label=u"擬合數據")
pl.legend()
pl.show()
```
這個例子中我們要擬合的函數是一個正弦波函數，它有三個參數 **A**, **k**, **theta** ，分別對應振幅、頻率、相角。假設我們的實驗數據是一組包含噪聲的數據 x, y1，其中y1是在真實數據y0的基礎上加入噪聲的到了。
通過leastsq函數對帶噪聲的實驗數據x, y1進行數據擬合，可以找到x和真實數據y0之間的正弦關系的三個參數： A, k, theta。下面是程序的輸出：
```
# -*- coding: utf-8 -*-
# 本程序用各種fmin函數求卷積的逆運算
import scipy.optimize as opt
import numpy as np
def test_fmin_convolve(fminfunc, x, h, y, yn, x0):
"""
x (*) h = y, (*)表示卷積
yn為在y的基礎上添加一些干擾噪聲的結果
x0為求解x的初始值
"""
def convolve_func(h):
"""
計算 yn - x (*) h 的power
fmin將通過計算使得此power最小
"""
return np.sum((yn - np.convolve(x, h))**2)
# 調用fmin函數，以x0為初始值
h0 = fminfunc(convolve_func, x0)
print fminfunc.__name__
print "---------------------"
# 輸出 x (*) h0 和 y 之間的相對誤差
print "error of y:", np.sum((np.convolve(x, h0)-y)**2)/np.sum(y**2)
# 輸出 h0 和 h 之間的相對誤差
print "error of h:", np.sum((h0-h)**2)/np.sum(h**2)
print
def test_n(m, n, nscale):
"""
隨機產生x, h, y, yn, x0等數列，調用各種fmin函數求解b
m為x的長度, n為h的長度, nscale為干擾的強度
"""
x = np.random.rand(m)
h = np.random.rand(n)
y = np.convolve(x, h)
yn = y + np.random.rand(len(y)) * nscale
x0 = np.random.rand(n)
test_fmin_convolve(opt.fmin, x, h, y, yn, x0)
test_fmin_convolve(opt.fmin_powell, x, h, y, yn, x0)
test_fmin_convolve(opt.fmin_cg, x, h, y, yn, x0)
test_fmin_convolve(opt.fmin_bfgs, x, h, y, yn, x0)
if __name__ == "__main__":
test_n(200, 20, 0.1)
```
下面是程序的輸出：
```
fmin
ーーーーーーーーーーー
error of y: 0.00568756699607
error of h: 0.354083287918
fmin_powell
ーーーーーーーーーーー
error of y: 0.000116114709857
error of h: 0.000258897894009
fmin_cg
ーーーーーーーーーーー
error of y: 0.000111220299615
error of h: 0.000211404733439
fmin_bfgs
ーーーーーーーーーーー
error of y: 0.000111220251551
error of h: 0.000211405138529
```
## 非線性方程組求解
optimize庫中的fsolve函數可以用來對非線性方程組進行求解。它的基本調用形式如下：
```
fsolve(func, x0)
```
func(x)是計算方程組誤差的函數，它的參數x是一個矢量，表示方程組的各個未知數的一組可能解，func返回將x代入方程組之后得到的誤差；x0為未知數矢量的初始值。如果要對如下方程組進行求解的話：
* f1(u1,u2,u3) = 0
* f2(u1,u2,u3) = 0
* f3(u1,u2,u3) = 0
那么func可以如下定義：
```
def func(x):
u1,u2,u3 = x
return [f1(u1,u2,u3), f2(u1,u2,u3), f3(u1,u2,u3)]
```
下面是一個實際的例子，求解如下方程組的解：
* 5*x1 + 3 = 0
* 4*x0*x0 - 2*sin(x1*x2) = 0
* x1*x2 - 1.5 = 0
程序如下：
```
from scipy.optimize import fsolve
from math import sin,cos
def f(x):
x0 = float(x[0])
x1 = float(x[1])
x2 = float(x[2])
return [
5*x1+3,
4*x0*x0 - 2*sin(x1*x2),
x1*x2 - 1.5
]
result = fsolve(f, [1,1,1])
print result
print f(result)
```
輸出為：
```
[-0.70622057 -0.6 -2.5 ]
[0.0, -9.1260332624187868e-14, 5.3290705182007514e-15]
```
由于fsolve函數在調用函數f時，傳遞的參數為數組，因此如果直接使用數組中的元素計算的話，計算速度將會有所降低，因此這里先用float函數將數組中的元素轉換為Python中的標準浮點數，然后調用標準math庫中的函數進行運算。
在對方程組進行求解時，fsolve會自動計算方程組的雅可比矩陣，如果方程組中的未知數很多，而與每個方程有關的未知數較少時，即雅可比矩陣比較稀疏時，傳遞一個計算雅可比矩陣的函數將能大幅度提高運算速度。筆者在一個模擬計算的程序中需要大量求解近有50個未知數的非線性方程組的解。每個方程平均與6個未知數相關，通過傳遞雅可比矩陣的計算函數使計算速度提高了4倍。
雅可比矩陣
雅可比矩陣是一階偏導數以一定方式排列的矩陣，它給出了可微分方程與給定點的最優線性逼近，因此類似于多元函數的導數。例如前面的函數f1,f2,f3和未知數u1,u2,u3的雅可比矩陣如下：
![\begin{bmatrix}
\dfrac{\partial f1}{\partial u1} & \dfrac{\partial f1}{\partial u2} & \dfrac{\partial f1}{\partial u3} \\[9pt]
\dfrac{\partial f2}{\partial u1} & \dfrac{\partial f2}{\partial u2} & \dfrac{\partial f2}{\partial u3} \\[9pt]
\dfrac{\partial f3}{\partial u1} & \dfrac{\partial f3}{\partial u2} & \dfrac{\partial f3}{\partial u3} \\
\end{bmatrix}](img/0e0589610ac2875dcd7d12db88259d4076bd8fe7.png)
使用雅可比矩陣的fsolve實例如下，計算雅可比矩陣的函數j通過fprime參數傳遞給fsolve，函數j和函數f一樣，有一個未知數的解矢量參數x，函數j計算非線性方程組在矢量x點上的雅可比矩陣。由于這個例子中未知數很少，因此程序計算雅可比矩陣并不能帶來計算速度的提升。
```
# -*- coding: utf-8 -*-
from scipy.optimize import fsolve
from math import sin,cos
def f(x):
x0 = float(x[0])
x1 = float(x[1])
x2 = float(x[2])
return [
5*x1+3,
4*x0*x0 - 2*sin(x1*x2),
x1*x2 - 1.5
]
def j(x):
x0 = float(x[0])
x1 = float(x[1])
x2 = float(x[2])
return [
[0, 5, 0],
[8*x0, -2*x2*cos(x1*x2), -2*x1*cos(x1*x2)],
[0, x2, x1]
]
result = fsolve(f, [1,1,1], fprime=j)
print result
print f(result)
```
## B-Spline樣條曲線
interpolate庫提供了許多對數據進行插值運算的函數。下面是使用直線和B-Spline對正弦波上的點進行插值的例子。
```
# -*- coding: utf-8 -*-
import numpy as np
import pylab as pl
from scipy import interpolate
x = np.linspace(0, 2*np.pi+np.pi/4, 10)
y = np.sin(x)
x_new = np.linspace(0, 2*np.pi+np.pi/4, 100)
f_linear = interpolate.interp1d(x, y)
tck = interpolate.splrep(x, y)
y_bspline = interpolate.splev(x_new, tck)
pl.plot(x, y, "o", label=u"原始數據")
pl.plot(x_new, f_linear(x_new), label=u"線性插值")
pl.plot(x_new, y_bspline, label=u"B-spline插值")
pl.legend()
pl.show()
```
使用interpolate庫對正弦波數據進行線性插值和B-Spline插值
在這段程序中，通過interp1d函數直接得到一個新的線性插值函數。而B-Spline插值運算需要先使用splrep函數計算出B-Spline曲線的參數，然后將參數傳遞給splev函數計算出各個取樣點的插值結果。
## 數值積分
數值積分是對定積分的數值求解，例如可以利用數值積分計算某個形狀的面積。下面讓我們來考慮一下如何計算半徑為1的半圓的面積，根據圓的面積公式，其面積應該等于PI/2。單位半圓曲線可以用下面的函數表示：
```
def half_circle(x):
return (1-x**2)**0.5
```
下面的程序使用經典的分小矩形計算面積總和的方式，計算出單位半圓的面積：
```
>>> N = 10000
>>> x = np.linspace(-1, 1, N)
>>> dx = 2.0/N
>>> y = half_circle(x)
>>> dx * np.sum(y[:-1] + y[1:]) # 面積的兩倍
3.1412751679988937
```
利用上述方式計算出的圓上一系列點的坐標，還可以用numpy.trapz進行數值積分：
```
>>> import numpy as np
>>> np.trapz(y, x) * 2 # 面積的兩倍
3.1415893269316042
```
此函數計算的是以x,y為頂點坐標的折線與X軸所夾的面積。同樣的分割點數，trapz函數的結果更加接近精確值一些。
如果我們調用scipy.integrate庫中的quad函數的話，將會得到非常精確的結果：
```
>>> from scipy import integrate
>>> pi_half, err = integrate.quad(half_circle, -1, 1)
>>> pi_half*2
3.1415926535897984
```
多重定積分的求值可以通過多次調用quad函數實現，為了調用方便，integrate庫提供了dblquad函數進行二重定積分，tplquad函數進行三重定積分。下面以計算單位半球體積為例說明dblquad函數的用法。
單位半球上的點(x,y,z)符合如下方程：
因此可以如下定義通過(x,y)坐標計算球面上點的z值的函數：
```
def half_sphere(x, y):
return (1-x**2-y**2)**0.5
```
X-Y軸平面與此球體的交線為一個單位圓，因此積分區間為此單位圓，可以考慮為X軸坐標從-1到1進行積分，而Y軸從 -half_circle(x) 到 half_circle(x) 進行積分，于是可以調用dblquad函數：
```
>>> integrate.dblquad(half_sphere, -1, 1,
lambda x:-half_circle(x),
lambda x:half_circle(x))
>>> (2.0943951023931988, 2.3252456653390915e-14)
>>> np.pi*4/3/2 # 通過球體體積公式計算的半球體積
2.0943951023931953
```
dblquad函數的調用方式為：
```
dblquad(func2d, a, b, gfun, hfun)
```
對于func2d(x,y)函數進行二重積分，其中a,b為變量x的積分區間，而gfun(x)到hfun(x)為變量y的積分區間。
半球體積的積分的示意圖如下：
半球體積的雙重定積分示意圖
X軸的積分區間為-1.0到1.0，對于X=x0時，通過對Y軸的積分計算出切面的面積，因此Y軸的積分區間如圖中紅色點線所示。
## 解常微分方程組
scipy.integrate庫提供了數值積分和常微分方程組求解算法odeint。下面讓我們來看看如何用odeint計算洛侖茲吸引子的軌跡。洛侖茲吸引子由下面的三個微分方程定義：
洛侖茲吸引子的詳細介紹: [http://bzhang.lamost.org/website/archives/lorenz_attactor](http://bzhang.lamost.org/website/archives/lorenz_attactor)
這三個方程定義了三維空間中各個坐標點上的速度矢量。從某個坐標開始沿著速度矢量進行積分，就可以計算出無質量點在此空間中的運動軌跡。其中 , ,  為三個常數，不同的參數可以計算出不同的運動軌跡： x(t), y(t), z(t)。 當參數為某些值時，軌跡出現餛飩現象：即微小的初值差別也會顯著地影響運動軌跡。下面是洛侖茲吸引子的軌跡計算和繪制程序：
```
# -*- coding: utf-8 -*-
from scipy.integrate import odeint
import numpy as np
def lorenz(w, t, p, r, b):
# 給出位置矢量w，和三個參數p, r, b計算出
# dx/dt, dy/dt, dz/dt的值
x, y, z = w
# 直接與lorenz的計算公式對應
return np.array([p*(y-x), x*(r-z)-y, x*y-b*z])
t = np.arange(0, 30, 0.01) # 創建時間點
# 調用ode對lorenz進行求解, 用兩個不同的初始值
track1 = odeint(lorenz, (0.0, 1.00, 0.0), t, args=(10.0, 28.0, 3.0))
track2 = odeint(lorenz, (0.0, 1.01, 0.0), t, args=(10.0, 28.0, 3.0))
# 繪圖
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
fig = plt.figure()
ax = Axes3D(fig)
ax.plot(track1[:,0], track1[:,1], track1[:,2])
ax.plot(track2[:,0], track2[:,1], track2[:,2])
plt.show()
```
用odeint函數對洛侖茲吸引子微分方程進行數值求解所得到的運動軌跡
我們看到即使初始值只相差0.01，兩條運動軌跡也是完全不同的。
在程序中先定義一個lorenz函數，它的任務是計算出某個位置的各個方向的微分值，這個計算直接根據洛侖茲吸引子的公式得出。然后調用odeint，對微分方程求解，odeint有許多參數，這里用到的四個參數分別為：
1. lorenz， 它是計算某個位移上的各個方向的速度(位移的微分)
2. (0.0, 1.0, 0.0)，位移初始值。計算常微分方程所需的各個變量的初始值
3. t， 表示時間的數組，odeint對于此數組中的每個時間點進行求解，得出所有時間點的位置
4. args， 這些參數直接傳遞給lorenz函數，因此它們都是常量
## 濾波器設計
scipy.signal庫提供了許多信號處理方面的函數。在這一節，讓我們來看看如何利用signal庫設計濾波器，查看濾波器的頻率響應，以及如何使用濾波器對信號進行濾波。
假設如下導入signal庫:
```
>>> import scipy.signal as signal
```
下面的程序設計一個帶通IIR濾波器：
```
>>> b, a = signal.iirdesign([0.2, 0.5], [0.1, 0.6], 2, 40)
```
這個濾波器的通帶為0.2*f0到0.5*f0，阻帶為小于0.1*f0和大于0.6*f0，其中f0為1/2的信號取樣頻率，如果取樣頻率為8kHz的話，那么這個帶通濾波器的通帶為800Hz到2kHz。通帶的最大增益衰減為2dB，阻帶的最小增益衰減為40dB，即通帶的增益浮動在2dB之內，阻帶至少有40dB的衰減。
iirdesgin返回的兩個數組b和a， 它們分別是IIR濾波器的分子和分母部分的系數。其中a[0]恒等于1。
下面通過調用freqz計算所得到的濾波器的頻率響應：
```
>>> w, h = signal.freqz(b, a)
```
freqz返回兩個數組w和h，其中w是圓頻率數組，通過w/pi*f0可以計算出其對應的實際頻率。h是w中的對應頻率點的響應，它是一個復數數組，其幅值為濾波器的增益，相角為濾波器的相位特性。
下面計算h的增益特性，并轉換為dB度量。由于h中存在幅值幾乎為0的值，因此先用clip函數對其裁剪之后，再調用對數函數，避免計算出錯。
```
>>> power = 20*np.log10(np.clip(np.abs(h), 1e-8, 1e100))
```
通過下面的語句可以繪制出濾波器的增益特性圖，這里假設取樣頻率為8kHz：
```
>>> pl.plot(w/np.pi*4000, power)
```
在實際運用中為了測量未知系統的頻率特性，經常將頻率掃描波輸入到系統中，觀察系統的輸出，從而計算其頻率特性。下面讓我們來模擬這一過程。
為了調用chirp函數以產生頻率掃描波形的數據，首先需要產生一個等差數組代表取樣時間，下面的語句產生2秒鐘取樣頻率為8kHz的取樣時間數組：
```
>>> t = np.arange(0, 2, 1/8000.0)
```
然后調用chirp得到2秒鐘的頻率掃描波形的數據：
```
>>> sweep = signal.chirp(t, f0=0, t1 = 2, f1=4000.0)
```
頻率掃描波的開始頻率f0為0Hz，結束頻率f1為4kHz，到達4kHz的時間為2秒，使用數組t作為取樣時間點。
下面通過調用lfilter函數計算sweep波形經過帶通濾波器之后的結果：
```
>>> out = signal.lfilter(b, a, sweep)
```
lfilter內部通過如下算式計算IIR濾波器的輸出：
通過如下算式可以計算輸入為x時的濾波器的輸出，其中數組x代表輸入信號，y代表輸出信號：
![y[n] & = b[0] x[n] + b[1] x[n-1] + \cdots + b[P] x[n-P] \\
& - a[1] y[n-1] - a[2] y[n-2] - \cdots - a[Q] y[n-Q]](img/8c429be55d6571ecd8bf39cf20455defe09eadcd.png)
為了和系統的增益特性圖進行比較，需要獲取輸出波形的包絡，因此下面先將輸出波形數據轉換為能量值：
```
>>> out = 20*np.log10(np.abs(out))
```
為了計算包絡，找到所有能量大于前后兩個取樣點(局部最大點)的下標：
```
>>> index = np.where(np.logical_and(out[1:-1] > out[:-2], out[1:-1] > out[2:]))[0] + 1
```
最后將時間轉換為對應的頻率，繪制所有局部最大點的能量值：
```
>>> pl.plot(t[index]/2.0*4000, out[index] )
```
下圖顯示freqz計算的頻譜和頻率掃描波得到的頻率特性，我們看到其結果是一致的。
帶通IIR濾波器的頻率響應和頻率掃描波計算的結果比較
計算此圖的完整源程序請查看附錄中的 [_帶通濾波器設計_](example_scipy_signal.html) 。
## 用Weave嵌入C語言
Python作為動態語言其功能雖然強大，但是在數值計算方面有一個最大的缺點：速度不夠快。在Python級別的循環和計算的速度只有C語言程序的百分之一。因此才有了NumPy, SciPy這樣的函數庫，將高度優化的C、Fortran的函數庫進行包裝，以供Python程序調用。如果這些高度優化的函數庫無法實現我們的算法，必須從頭開始寫循環、計算的話，那么用Python來做顯然是不合適的。因此SciPy提供了快速調用C++語言程序的方法-- Weave。下面是對NumPy的數組求和的例子：
```
# -*- coding: utf-8 -*-
import scipy.weave as weave
import numpy as np
import time
def my_sum(a):
n=int(len(a))
code="""
int i;
double counter;
counter =0;
for(i=0;i<n;i++){
counter=counter+a(i);
}
return_val=counter;
"""
err=weave.inline(
code,['a','n'],
type_converters=weave.converters.blitz,
compiler="gcc"
)
return err
a = np.arange(0, 10000000, 1.0)
# 先調用一次my_sum，weave會自動對C語言進行編譯，此后直接運行編譯之后的代碼
my_sum(a)
start = time.clock()
for i in xrange(100):
my_sum(a) # 直接運行編譯之后的代碼
print "my_sum:", (time.clock() - start) / 100.0
start = time.clock()
for i in xrange(100):
np.sum( a ) # numpy中的sum，其實現也是C語言級別
print "np.sum:", (time.clock() - start) / 100.0
start = time.clock()
print sum(a) # Python內部函數sum通過數組a的迭代接口訪問其每個元素，因此速度很慢
print "sum:", time.clock() - start
```
此例子在我的電腦上的運行結果為：
```
>>> from sympy import *
```
## 封面上的經典公式
本書的封面上的公式：
叫做歐拉恒等式，其中e是自然指數的底，i是虛數單位，  是圓周率。此公式被譽為數學最奇妙的公式，它將5個基本數學常數用加法、乘法和冪運算聯系起來。下面用SymPy驗證一下這個公式。
載入的符號中，E表示自然指數的底，I表示虛數單位，pi表示圓周率，因此上述的公式可以直接如下計算：
```
>>> E**(I*pi)+1
0
```
歐拉恒等式可以下面的公式進行計算，
為了用SymPy求證上面的公式，我們需要引入變量x。在SymPy中，數學符號是Symbol類的對象，因此必須先創建之后才能使用：
```
>>> x = Symbol('x')
```
expand函數可以將公式展開，我們用它來展開E**(I*pi)試試看：
```
>>> expand( E**(I*x) )
exp(I*x)
```
沒有成功，只是換了一種寫法而已。這里的exp不是math.exp或者numpy.exp，而是sympy.exp，它是一個類，用來表述自然指數函數。
expand函數有關鍵字參數complex，當它為True時，expand將把公式分為實數和虛數兩個部分：
```
>>> expand(exp(I*x), complex=True)
I*exp(-im(x))*sin(re(x)) + cos(re(x))*exp(-im(x))
```
這次得到的結果相當復雜，其中sin, cos, re, im都是sympy定義的類，re表示取實數部分，im表示取虛數部分。顯然這里的運算將符號x當作復數了。為了指定符號x必須是實數，我們需要如下重新定義符號x：
```
>>> x = Symbol("x", real=True)
>>> expand(exp(I*x), complex=True)
I*sin(x) + cos(x)
```
終于得到了我們需要的公式。那么如何證明它呢。我們可以用泰勒多項式展開：
```
>>> tmp = series(exp(I*x), x, 0, 10)
>>> pprint(tmp)
2 3 4 5 6 7 8 9
x I*x x I*x x I*x x I*x
1 + I*x - -- - ---- + -- + ---- - --- - ---- + ----- + ------ + O(x**10)
2 6 24 120 720 5040 40320 362880
```
series是泰勒展開函數，pprint將公式用更好看的格式打印出來。下面分別獲得tmp的實部和虛部，分別和cos(x)和sin(x)的展開公式進行比較：
> ```
> >>> pprint(re(tmp))
> 2 4 6 8
> x x x x
> 1 + re(O(x**10)) - -- + -- - --- + -----
> 2 24 720 40320
>
> ```
>
> ```
> >>> pprint( series( cos(x), x, 0, 10) )
> 2 4 6 8
> x x x x
> 1 - -- + -- - --- + ----- + O(x**10)
> 2 24 720 40320
>
> ```
>
> ```
> >>> pprint(im(tmp))
> 3 5 7 9
> x x x x
> x + im(O(x**10)) - -- + --- - ---- + ------
> 6 120 5040 362880
>
> ```
>
> ```
> >>> pprint(series(sin(x), x, 0, 10))
> 3 5 7 9
> x x x x
> x - -- + --- - ---- + ------ + O(x**10)
> 6 120 5040 362880
>
> ```
## 球體體積
在[_用SciPy數值積分_](scipy_intro.html#sec-spherevolume)一節我們介紹了如何使用數值定積分計算球體的體積，而SymPy的符號積分函數integrate則可以幫助我們進行符號積分。integrate可以進行不定積分：
```
>>> integrate(x*sin(x), x)
-x*cos(x) + sin(x)
```
如果指定x的取值范圍的話，integrate則進行定積分運算：
```
>>> integrate(x*sin(x), (x, 0, 2*pi))
-2*pi
```
為了計算球體體積，首先讓我們來看看如何計算圓形面積，假設圓形的半徑為r，則圓上任意一點的Y坐標函數為：
因此我們可以直接對上述函數在-r到r區間上進行積分得到半圓面積，注意這里我們使用symbols函數一次創建多個符號：
```
>>> x, y, r = symbols('x,y,r')
>>> 2 * integrate(sqrt(r*r-x**2), (x, -r, r))
2*Integral((r**2 - x**2)**(1/2), (x, -r, r))
```
很遺憾，integrate函數沒有計算出結果，而是直接返回了我們輸入的算式。這是因為SymPy不知道r是大于0的，如下重新定義r，就可以得到正確答案了：
```
>>> r = symbols('r', positive=True)
>>> circle_area = 2 * integrate(sqrt(r**2-x**2), (x, -r, r))
>>> circle_area
pi*r**2
```
接下來對此面積公式進行定積分，就可以得到球體的體積，但是隨著X軸坐標的變化，對應的切面的的半徑會發生變化，現在假設X軸的坐標為x，球體的半徑為r，則x處的切面的半徑為可以使用前面的公式y(x)計算出。
球體體積的雙重定積分示意圖
因此我們需要對circle_area中的變量r進行替代：
```
>>> circle_area = circle_area.subs(r, sqrt(r**2-x**2))
>>> circle_area
pi*(r**2 - x**2)
```
用subs進行算式替換
subs函數可以將算式中的符號進行替換，它有3種調用方式：
* expression.subs(x, y) : 將算式中的x替換成y
* expression.subs({x:y,u:v}) : 使用字典進行多次替換
* expression.subs([(x,y),(u,v)]) : 使用列表進行多次替換
請注意多次替換是順序執行的，因此：
```
expression.sub([(x,y),(y,x)])
```
并不能對兩個符號x,y進行交換。
然后對circle_area中的變量x在區間-r到r上進行定積分，得到球體的體積公式：
```
>>> integrate(circle_area, (x, -r, r))
4*pi*r**3/3
```
- 用Python做科學計算
- 軟件包的安裝和介紹
- NumPy-快速處理數據
- SciPy-數值計算庫
- matplotlib-繪制精美的圖表
- Traits-為Python添加類型定義
- TraitsUI-輕松制作用戶界面
- Chaco-交互式圖表
- TVTK-三維可視化數據
- Mayavi-更方便的可視化
- Visual-制作3D演示動畫
- OpenCV-圖像處理和計算機視覺
- Traits使用手冊
- 定義Traits
- Trait事件處理
- 設計自己的Trait編輯器
- Visual使用手冊
- 場景窗口
- 聲音的輸入輸出
- 數字信號系統
- FFT演示程序
- 頻域信號處理
- Ctypes和NumPy
- 自適應濾波器和NLMS模擬
- 單擺和雙擺模擬
- 分形與混沌
- 關于本書的編寫
- 最近更新
- 源程序集
- 三角波的FFT演示
- 在traitsUI中使用的matplotlib控件
- CSV文件數據圖形化工具
- NLMS算法的模擬測試
- 三維標量場觀察器
- 頻譜泄漏和hann窗
- FFT卷積的速度比較
- 二次均衡器設計
- 單擺擺動周期的計算
- 雙擺系統的動畫模擬
- 繪制Mandelbrot集合
- 迭代函數系統的分形
- 繪制L-System的分形圖