SciPy 數值計算庫（一）

假設有一組實驗資料(x[i], y[i])，我們知道它們之間的函式關係:y = f(x)，通過這些已知資訊，需要確定函式中的一些引數項。例如，如果f是乙個線型函式f(x) = k*x+b，那麼引數k和b就是我們需要確定的值。如果將這些引數用 p 表示的話，那麼我們就是要找到一組 p 值使得如下公式中的s函式最小： s(

p)=∑

i=1m

[yi−

f(xi

,p)]

2 s(p

)=∑i

=1m[

yi−f

(xi,

p)]2

這種演算法被稱之為最小二乘擬合(least-square fitting)。

scipy中的子函式庫optimize已經提供了實現最小二乘擬合演算法的函式leastsq。下面是用leastsq進行資料擬合的乙個例子,**如下

# -*- coding: "utf-8 -*-
#coding=utf-8
import numpy as np
from scipy.optimize import leastsq
import pylab as pl
from pylab import mpl 
mpl.rcparams['font.sans-serif'] = ['simhei'] 
deffunc
(x, p):
a, k, theta = p
return a*np.sin(2*np.pi*k*x+theta)
defresiduals
(p, y, x):
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()

將程式拷貝到test.py檔案中，執行此程式：

發現執行結果，圖中的中文全部出現方框，無法正常顯示漢字，在程式開頭加上

from pylab import mpl 
mpl.rcparams['font.sans-serif'] = ['simhei'] 
mpl.rcparams['axes.unicode_minus'] = false

即可解決圖表中出現的漢字亂碼問題，再次執行結果如下所示：

我們看到擬合引數雖然和真實引數完全不同，但是由於正弦函式具有週期性，實際上擬合引數得到的函式和真實引數對應的函式是一致的。

optimize庫中的fsolve函式可以用來對非線性方程組進行求解。它的基本呼叫形式如下：

fsolve
(func, x0)

func(x)是計算方程組誤差的函式，它的引數x是乙個向量，表示方程組的各個未知數的一組可能解，func返回將x代入方程組之後得到的誤差；x0為未知數向量的初始值。如果要對如下方程組進行求解的話：

那麼func可以如下定義：

def
func
(x):
u1,u2,u3 = x
return [f1(u1,u2,u3), f2(u1,u2,u3), f3(u1,u2,u3)]

下面舉例，求解如下方程的解：

程式如下：

# -*- coding: utf-8 -*- 
from scipy.optimize import fsolve
from math import sin,cos
deff
(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)

將此程式拷貝到fslove.py檔案中，執行結果如下：

由於fsolve函式在呼叫函式f時，傳遞的引數為陣列，因此如果直接使用陣列中的元素計算的話，計算速度將會有所降低，因此這裡先用float函式將陣列中的元素轉換為python中的標準浮點數，然後呼叫標準math庫中的函式進行運算。

interpolate庫提供了許多對資料進行插值運算的函式。下面是使用直線和b-spline對正弦波上的點進行插值的例子。

程式如下：

# -*- coding: utf-8 -*-
import numpy as np
import pylab as pl
from scipy import interpolate
from pylab import mpl #解決圖中漢字亂碼問題
mpl.rcparams['font.sans-serif'] = ['simhei'] 
mpl.rcparams['axes.unicode_minus'] = false
#使圖表中負號正常輸出
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函式計算出各個取樣點的插值結果。將此程式儲存在inter.py檔案中，執行結果如下：

SciPy 數值計算庫（一）

Scipy 高階科學計算

Scipy空間計算凸包（convexHull

scipy庫中的odeint函式

SciPy 數值計算庫（一）

Scipy 高階科學計算

Scipy空間 計算凸包（convexHull

scipy庫中的odeint函式

相關推薦

Scipy空間計算凸包（convexHull