1. 程式人生 > >Scipy——數值計算庫

Scipy——數值計算庫

SciPy函式庫在NumPy庫的基礎上增加了眾多的數學、科學以及工程計算中常用的庫函式。例如線性代數、常微分方程數值求解、訊號處理、影象處理、稀疏矩陣等等。

1.最小二乘擬合

假設有一組實驗資料(x[i], y[i]),我們知道它們之間的函式關係:y = f(x),通過這些已知資訊,需要確定函式中的一些引數項。例如,如果f是一個線型函式f(x) = k*x+b,那麼引數k和b就是我們需要確定的值。如果將這些引數用 p 表示的話,那麼我們就是要找到一組 p 值使得如下公式中的S函式最小:
這裡寫圖片描述
這種演算法被稱之為最小二乘擬合(Least-square fitting)。

2.函式最小值

optimize庫提供了幾個求函式最小值的演算法:fmin, fmin_powell, fmin_cg, fmin_bfg。
可以使用fmin計算反捲積。

3.非線性方程組

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

在對方程組進行求解時,fsolve會自動計算方程組的雅可比矩陣,如果方程組中的未知數很多,而與每個方程有關的未知數較少時,即雅可比矩陣比較稀疏時,傳遞一個計算雅可比矩陣的函式將能大幅度提升運算速度。
雅可比矩陣:是一階偏導數以一定方式排列的矩陣,它給出了可微分方程與給定點的最優線性逼近,因此類似於多元函式的導數。例如函式f1,f2,f3和未知數u1,u2,u3的雅可比矩陣如下:)
這裡寫圖片描述

4.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

這裡寫圖片描述
在這段程式中,通過interp1d函式直接得到一個新的線性插值函式。而B-Spline插值運算需要先使用splrep函式計算出B-Spline曲線的引數,然後將引數傳遞給splev函式計算出各個取樣點的插值結果。

5.數值積分

數值積分是對定積分的數值求解,例如可以利用數值積分計算某個形狀的面積。
多重定積分的求值可以通過多次呼叫quad函式實現,為了呼叫方便,integrate庫提供了dblquad函式進行二重定積分,tplquad函式進行三重積分。

6.解常微分方程組

scipy.integrate庫提供了數值積分和常微分方程組求解演算法odeint。用odeint
計算洛侖茲吸引子的軌跡。洛侖茲吸引子由下面的三個微分方程定義:
這裡寫圖片描述
這三個方程定義了三維空間中各個座標點上的速度向量。從某個座標開始沿著速度向量進行積分,就可以計算出無質量點在此空間中的運動軌跡。其中 σ, ρ, β 為三個常數,不同的引數可以計算出不同的運動軌跡: 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()

這裡寫圖片描述
我們看到即使初始值只相差0.01,兩條運動軌跡也是完全不同的。
在程式中先定義一個lorenz函式,它的任務是計算出某個位置的各個方向的微分值,這個計算直接根據洛侖茲吸引子的公式得出。然後呼叫odeint,對微分方程求解,odeint有許多引數,這裡用到的四個引數分別為:
① lorenz, 它是計算某個位移上的各個方向的速度(位移的微分)
② (0.0, 1.0, 0.0),位移初始值。計算常微分方程所需的各個變數的初始值
③ t, 表示時間的陣列,odeint對於此陣列中的每個時間點進行求解,得出所有時間點的位置
④args, 這些引數直接傳遞給lorenz函式,因此它們都是常量。

7.濾波器設計

scipy.signal庫提供了許多訊號處理方面的函式。我們可以利用signal庫設計濾波器,檢視濾波器的頻率響應,以及使用濾波器對訊號進行濾波。
####8.用Weave嵌入C語言
Python作為動態語言其功能雖然強大,但是在數值計算方面有一個最大的缺點:速度不夠快。在Python級別的迴圈和計算的速度只有C語言程式的百分之一。因此才有了NumPy, SciPy這樣的函式庫,將高度優化的C、Fortran的函式庫進行包裝,以供Python程式呼叫。如果這些高度優化的函式庫無法實現我們的演算法,必須從頭開始寫迴圈、計算的話,那麼用Python來做顯然是不合適的。因此SciPy提供了快速呼叫C++語言程式的方法-- Weave。
用Weave編譯的C語言程式比numpy自帶的sum函式還要快。而Python的內部函式sum使用陣列的迭代器介面進行運算,因此速度是Python語言級別的,只有Weave版本的1/300。