Scipy在Numpy的基础上则加了众多的数学计算,科学计算以及工程计算中常用的模块,例如线性代数,常微分方程的数值求解,信号处理,图像处理,系数矩阵等。在本章中,将通过实例介绍Scipy中常用的的一些模块。为了方便读者理解,在示例程序中使用matplotlib,TVTK以及Mayavi等扩展绘制二维以及三维图表。
Scipy的special模块是一个非常完整的函数库,其中包含了基本数学函数,特殊数学函数以及Numpy中出现的所有函数。
伽马(gamma)函数γ时概率统计学中经常出现的一个函数,它计算公式如下:
显然这样计算起来特别的麻烦,幸运的是。scipy.special中有内置的gamma模块。
>>> import scipy.special as S
>>> S.gamma(4)
6.0
>>> S.gamma(0.5)
1.7724538509055159
>>> S.gamma(1+1j)
(0.4980156681183554-0.15494982830181081j) >>> S.gamma(1000) inf
Γ(z)函数时结成函数在实数和复数系上的扩展,他的增长速度特别的块,1000的阶乘就超过了双精度浮点数的表示范围,因此结果就是无穷大。为了计算更大的范围可以使用gammaln()计算ln(|Γ(x)|)的只,它使用特殊的算法,能够直接计算Γ函数的对数值,因此可以表示更大的范围。
Scipy的constants模块中包含了众多的物理常数:
>>> from scipy import constants as C
>>> C.c # 真空中的光速
299792458.0
>>> C.h #普朗克常数 6.62607004e-34
special模块中的某些函数并不是数学意义上的特殊函数,例如log1p(x)计算log(1+x)的值。这是由于浮点数的精度有限,无法很精确的表示非常接近1 的常数。例如无法用浮点数表示1+1e-20的值,因此
>>> import math
>>> math.log(a)
0.0
,然而当使用log1p的时候,则可以很精确的计算。实际上当x非常小的时候,log1p(x)约等于x,这可以通过对log(1+x)进行泰勒级数展开进行证明。在后续的章节我们会学习如何使用符号计算库SymPy进行泰勒级数展开。
这些特殊函数与NumPy中的一般函数一样,都是ufunc函数,支持数组的广播运算。例如elipj(u,m)计算雅克比矩阵,他有两个参数u和m,返回四个值sn = sinΦ,cn = cosΦ,dn = (1-msin2Θ)-1/2 由于ellipj()支持广播运算。因此下面的程序在调用他时传递的两个参数的形状分别为(200,1)和(1,5),于是得到的四个结果数组的形状都是(200,5),图3-1显示了这些曲线。
Scipy的optimize模块提供了许多数值优化算法,本节对其中的非线性方程组求解,数据拟合,函数值等进行简单介绍
非线性方程组求解
fsolve()可以对非线性方程组进行求解,他的基本调用格式为fsolve(func,x0).其中func是计算方程组误差的函数,他的参数x是一个数组,其值为方程组的一组可能的解func返回将x带入方程组之后得到的每个方程的误差x0为未知数的一组初始值,假设要对下面的方程组进行求解:
f1(u1,u2,u3) = 0,f2(u1,u2,u3) = 0,f3(u1,u2,u3) = 0
那么finc函数可以定义为:
from math import sin,cos
from scipy import optimize def f(x): x0,x1,x2 = x.tolist() return [5*x1+3,4*x0*x0-2*sin(x1*x2),x1*x2-1.5] result = optimize.fsolve(f,[1,1,1]) print(result) print(f(result))
[-0.70622057 -0.6 -2.5 ]
[0.0, -9.126033262418787e-14, 5.329070518200751e-15]
f()时计算方程组的误差的函数,x参数是一组可能的解,fsolve()在调用f()时,传递给f()的参数是一个数组,先调用tolist,将其转化为标准浮点数列表,然后调用math函数中的模块进行计算。因为在进行单个数值的运算的时候,标准浮点类型比NumPy的浮点类型要快许多,所以把数值转化成标准浮点数类型,能缩短一些计算时间。调用fsolve()时,传递计算误差的函数f()以及未知数的初始值。
在对方程组进行求解的时候,fsolve()会自动计算方程组在某点对各个未知数变量的偏导数,这些偏导数组成一个二维数组,数学上称之为雅克比矩阵。如果方程组的未知数很多,而与每个方程有关联的未知数较少,即雅克比矩阵比较稀疏的时候,将计算雅克比矩阵的函数作为参数传递给fsolve(),这样可以大幅度提高运算速度,逼着在一个模拟计算的程序中需要求解有50个未知数的非线性方程组每个方程平均与6个未知数相关,通过传递计算雅克比矩阵的函数时fsolve()的计算速度提高了4 倍。
最小二乘拟合
假设有一组实验数据(xi,yj),我们实现知道他们应该满足某种函数关系yi = f(xi)。通过这些已知信息,需要确定函数f(x)的一些参数。例如函数f()时线性函数f(x) = kx+b那么参数k和b就是需要确定的值。如果用p表示函数中需要确定的参数,则目标是找到一组p使得函数S的值最小。
如果用p表示函数中需要确定的参数,则目标时找到一组p使得函数S的值最小:
这种方法叫做最小二乘拟合。在optimize模块中,可以受用leastsq()对数据进行最小二乘拟合运算。leastsq()的用法非常的简单。只需要将计算误差的函数和待确定参数的初始值传递给它即可。
import numpy as np
from scipy import optimize
x = np.array([8.19 , 2.72 , 6.39 , 8.71 , 4.7 , 2.66 , 3.78]) y = np.array([7.01 , 2.78 , 6.47 , 6.71 , 4.1 , 4.23 , 4.05]) def residuals(p): k,b = p return y - (k*x + b) r = optimize.leastsq(residuals,[1,0]) k,b = r[0] print("k = ",k,"b = ",b)
接下来看一个对正弦波进行拟合的例子。
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize #import pylab as pl def func(x,p): """ 数据拟合用的函数 """ 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]) # 真实数据 广播运算 np.random.seed(0) y1 = y0 + 2 * np.random.randn(len(x)) p0 = [7 , 0.4 , 0] # 第一次猜测的函数拟合参数 plsq = optimize.leastsq(residuals , p0 , args = (y1,x)) print("真实参数",[A,k,theta]) print("拟合参数",plsq[0])#实验数据 拟合之后的参数 plt.plot(x , y1 , "o" , label = u"Noisy data") plt.plot(x , y0 , label = u'Real data' ) plt.plot(x,func(x,plsq[0]),label = u'Fitting data') plt.legend(loc = 'best') plt.show()
真实参数 [10, 0.34, 0.5235987755982988]
拟合参数 [ 10.25218748 0.3423992 0.50817425]
和上一个不同的是,这里的leatsq有一个args参数。因为residuals函数 有三个参数,所以这里调用的时候也需要将这三个参数传入,当然也可以改变一下采用上面没有args的写法。