当前位置:首页  科学计算

《用Python进行科学计算》——SciPy数值计算库

SciPy函数库在NumPy库的基础上增加了众多的数学、科学以及工程计算中常用的库函数。例如线性代数、常微分方程数值求解、信号处理、图像处理、稀疏矩阵等。

最小二乘拟合

假设有一组实验数据(x[i],y[i]),我们知道他们之间的函数关系:y=f(x),通过这些已知信息,需要确定函数中的一些参数项。例如:如果f是一个线形函数f(x)=k*x+b,那么参数k和b就是我们需要确定的值。如果将这些参数用p表示的话,那么我们就要找到一组p值使得如下公式的S函数最小:

其中y为原函数,f为拟合函数。拟合函数可通过给定值找出表达式,不要求多精确,之要求能反映数据变化的趋势。
这种算法被称为最小二乘拟合。
scipy中的子函数库optimize已经提供了实现最小二乘拟合算法的函数leastsq。

#-*- coding: utf-8 -*-

import numpy as np
from scipy.optimize import leastsq
import matplotlib.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)

#在0~-2*np.pi均匀返回100个数字
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("Real_parameter:",[A,k,theta])
print("Fitting_parameters" ,plsq[0])

 

pl.plot(x,y0,label="RealData")
pl.plot(x,y1,label="NoiseData")
pl.plot(x,func(x,plsq[0]),label="FittingData")

pl.legend()
pl.show()

这里我们要拟合的函数是一个正弦波函数,它有参个参数A,k,theta分别对应振幅,频率,相角。假设实验数据是一组包含噪声的数据x,y1,其中y1是在真实数据的基础上加入噪声得到的。
通过leastsq函数对带噪声的实验数据x,y1进行数据拟合,可以找到x和真实数据y0之间的正弦关系的参数:A,k,theta。
 

Real_parameter: [10, 0.34, 0.5235987755982988]
Fitting_parameters [10.72946651  0.34093584  0.58900944]


函数最小值

optimize库提供了几个求函数最小值的算法:fmin,fmin_powell,fmin_cg,fim_bfgs。下面的程序通过求解卷积的逆运算岩石fmin的功能。

对于一个离散的线性时不变系统h,如果他的输入是x,那么其输出也可以用x和h卷积表示,即y=x*h。
现在的问题是如果已知系统的输入x和输出y,如何计算系统的传递函数h;或者如果已知系统的传递函数h和系统的输出y,如何计算系统的输入x。这种运算被称为反卷积运算。
 

#coding:utf-8

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:",float(np.sum(np.convolve(x,h0)-y)**2)/np.sum(y**2))

    #输出h0和h之间的相对误差
    print("error of h:",float(np.sum((h0-h)**2)/np.sum(h**2)))
    print('/n')

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.03524512840813301
error of h: 0.10872448809534953
/n
Optimization terminated successfully.
         Current function value: 0.214475
         Iterations: 32
         Function evaluations: 5999
fmin_powell
----------------------
error of y: 0.04085253762222568
error of h: 0.00025751019243286225
/n
Optimization terminated successfully.
         Current function value: 0.214441
         Iterations: 19
         Function evaluations: 858
         Gradient evaluations: 39
fmin_cg
----------------------
error of y: 0.04084548316927832
error of h: 0.00025921974956644503
/n
Optimization terminated successfully.
         Current function value: 0.214441
         Iterations: 32
         Function evaluations: 1012
         Gradient evaluations: 46
fmin_bfgs
----------------------
error of y: 0.0408454823459306
error of h: 0.00025922110838010043
/n

非线性方程组求解

optimize库中的fsolve函数可以用来对非线性方程组进行进行求解。
func(x)是计算方程组误差的函数,它的参数x是一个矢量,表示方程组的各个未知数的一组可能解,func返回将x带入方程组之后得到的误差;x0为未知数矢量的初始值。
实例:
求解

    5*x1+3=0
    4x0x0-2sin(x1x2)=0
    x1*x2-1.5=0
 

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

result=fsolve(f,[1,1,1])  #输入一组可能解

print(result)
print(f(result))  #返回误差

[-0.70622057 -0.6        -2.5       ]
[0.0, -9.126033262418787e-14, 5.329070518200751e-15]

在对方程组进行求解时,fsolve会自动计算方程组的雅可比矩阵,如果方程组中的未知数很多,而与每个方程有关的未知数较少时,传递一个计算雅可比矩阵的函数将能大幅度提高运算速度。
什么是雅可比矩阵?
雅可比矩阵是一阶偏导数以一定方式排列的矩阵,它给出可微分方程与给定点的最优线性逼近,因此类似于多元函数的导数。
使用雅可比矩阵的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 matplotlib.pylab as pl
import matplotlib.pyplot as plt
from scipy import interpolate

plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签


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)

#计算出B-Spline曲线的参数
tck=interpolate.splrep(x,y)

#将参数传递给splev函数计算出各个取样点的插值结果
y_bspline=interpolate.splev(x_new,tck)


pl.plot(x,y,"o",label="orginallyData")
pl.plot(x_new,f_linear(x_new),label="LineInsert")
pl.plot(x_new,y_bspline,label="B-splineInsert")

pl.legend()
pl.show()


数值积分

数值积分是对定积分的数值求解,例如可以利用数值积分计算某个形状的面积。
实例:
计算半径为1的半圆的面积。
单位半圆曲线可以用下面函数表示。

#coding:utf-8
import numpy as np

def half_circle(x):
    return (1-x**2)**0.5


#使用经典的分小矩形计算面积总和,计算出单位面积。

N=1000
x=np.linspace(-1,1,N)
y=half_circle(x)

#使用分割法计算
dx=2.0/N
S=dx*np.sum(y[:-1]+y[1:])

#进行数值积分
S1=np.trapz(y,x) *2
#此函数计算的是以x,y为顶点坐标的折线与X所夹的面积。同样的分割点数,结果更精确


#使用quad函数积分
from scipy import integrate

pi_half,err=integrate.quad(half_circle,-1,1)
S2=pi_half*2


print("Split_run: %s" %S)
print("Trapz: %s" %S1)
print("Quad: %s" %S2)

 

Split_run: 3.138345831727239
Trapz: 3.141487319046285
Quad: 3.1415926535897967

可以看出quad函数是非常精确的。
多重定积分的求值可以通过多次调用quad函数实现,为了方便,integrate库提供了dblquad函数进行二重积分,tplquad函数进行三重积分。
实例:
计算单位半球体积。
符合x**2+y**2+z**2=1

from scipy import integrate
import numpy as np

def half_circle(x):
    return (1-x**2)**0.5

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)
#进行积分

V=integrate.dblquad(half_sphere,-1,1,
                    lambda x:-half_circle(x),
                    lambda x:half_circle(x))

print(V)

#通过球体体积公式计算
print(np.pi*4/3/2)

 

(2.0943951023931984, 1.0002354500215915e-09)
2.0943951023931953

dblquad函数的调用方式:

dblquad(func2d,a,b,gfun,hfun)

对于func2d(x,y)函数进行二重积分,其中a,b为变量x的积分期间。

解常微分方程组

scipy.integrate库提供了数值积分和常微分方程组求解算法odeint。
实例:
计算洛伦兹吸引子的轨迹,洛伦兹吸引子有三个微分方程定义:

dx/dt=σ(y-x)
dy/dt=x(ρ-z)-y
dz/dt=xy-βz

这三个方程定义了三维空间中各个坐标点上的速度矢量。从某个坐标开始沿着速度矢量进行积分,就可以计算出无质量点在此空间中的运动轨迹。其中 σ, ρ, β 为三个常数,不同的参数可以计算出不同的
运动轨迹: 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/dtt,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],track1[:,1],track1[:,2])

plt.show()

在程序中先定义一个lozren函数,它的任务是计算出某个位置的各个方向的微分值,这个计算直接根据洛伦兹吸引子公式得出。但后调用odeint,对微分方程求解:
参数祥解:

    lorenz:计算某个位移上的各个方向的速度;
    (0.0,1.0,0.0),位移初始值。计算常微分方程所需的各个变量的初始值。
    t:表示时间的数组,odeint对于此数组中的每个时间点进行求解,的出所有时间点的位置。
    args:这些参数直接传递给lorenz函数,因此他们都是常量。

滤波器设计

scipy.signal库提供了西多信号处理反面的函数。如何利用signal库设计滤波器,查看滤波器的频率响应,以及如何使用滤波器对信号进行滤波。
设计一个带通滤波器:

import scipy.signal as signal

b,a=signal.iirdesign([0.2,0.5],[0.1,0.6],2,40)


#参数解析:[0.2,0.5]通带范围
#       [0.1,0.6]阻带小于0.1*f0大于0.6*f0
#       通带的最大增益衰减为2dB,阻带的最小增益衰减为40dB,
# 其中f0为1/2的信号取样频率
# iirdesign返回的两个数组b和a,他们分别是IIR滤波器的分子和分母部分的系数.
# 其中a[0]恒等于1

通过调用freqz计算所得到的滤波器的频率响应:

w,h=signal.freqz(b,a)

frez返回两个数组w和h,其中w是圆频率数组,通过w/pi*f0可以计算出其对应的实际频率.h是w中的对应频率点的响应.其幅值为滤波器的增益,相角为滤波器的相位特性.
计算h的增益特性,并转换为dB度量,由于h中存在幅值几乎为0的值,因此先用clip函数对其进行裁剪,在调用对数函数,避免计算出错。
 

#绘制出滤波器的增益特性图
pl.plot(w/np.pi*4000,power)

pl.show()


测量未知系统的频率特性
将频率扫描到系统中,观察系统的输出,从而计算其频率特性。

# 为了调用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代表输入信号,y代表输出信号:

y[n]=b[0]x[n]+b[1]x[n-1]+...+b[P]x[n-P]
-a[1]y[n-1]-a[2]y[n-2]-...-a[Q]y[n-Q]

获取输出波形的包络,所以先将输出波形数据转换为能量值:

# 为了调用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)


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

pl.show()


用Weave嵌入C语言

Python作为动态语言其功能虽然强大,但是在数值计算方面有一个最大的缺点:速度不够快。在Python级别的循环和计算的速度只有C语言程序的百分之一。因此才有了NumPy, SciPy这样的函数库,将高度优化的C、Fortran的函数库进行包装,以供Python程序调用。如果这些高度优化的函数库无法实现我们的算法,必须从头开始写循环、计算的话,那么用Python来做显然是不合适的。因此SciPy提供了快速调用C++语言程序的方法-- Weave。