SciPy-数值计算库

SciPy简介

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

Scipy函数库的应用及实例

  • 最小二乘拟合

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

《SciPy-数值计算库》

这种算法被称之为最小二乘拟合(Least-square fitting)。

  • scipy中的子函数库optimize已经提供了实现最小二乘拟合算法的函数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。下面是程序的输出:

>>> 真实参数: [10, 0.34000000000000002, 0.52359877559829882]
>>> 拟合参数 [-9.84152775  0.33829767 -2.68899335]

《SciPy-数值计算库》 调用leastsq函数对噪声正弦波数据进行曲线拟合

  • 函数最小值

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

对于一个离散的线性时不变系统如果它的输入是x,那么其输出y可以用x和h的卷积表示:y = x * h
现在的问题是如果已知系统的输入x和输出y,如何计算系统的传递函数h;或者如果已知系统的传递函数h和系统的输出y,如何计算系统的输入x。这种运算被称为反卷积运算,是十分困难的,特别是在实际的运用中,测量系统的输出总是存在误差的。

下面用fmin计算反卷积,这种方法只能用在很小规模的数列之上。

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

下面是一个实际的例子,求解如下方程组的解:

《SciPy-数值计算库》

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会自动计算方程组的雅可比矩阵,如果方程组中的未知数很多,而与每个方程有关的未知数较少时,即雅可比矩阵比较稀疏时,传递一个计算雅可比矩阵的函数将能大幅度提高运算速度。
  • 雅可比矩阵

雅可比矩阵是一阶偏导数以一定方式排列的矩阵,它给出了可微分方程与给定点的最优线性逼近,因此类似于多元函数的导数。例如前面的函数f1,f2,f3和未知数u1,u2,u3的雅可比矩阵如下:

《SciPy-数值计算库》

使用雅可比矩阵的fsolve实例如下

# -*- 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)

输出为:

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

计算雅可比矩阵的函数j通过fprime参数传递给fsolve,函数j和函数f一样,有一个未知数的解矢量参数x,函数j计算非线性方程组在矢量x点上的雅可比矩阵。由于这个例子中未知数很少,因此程序计算雅可比矩阵并不能带来计算速度的提升。

  • 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()

《SciPy-数值计算库》 使用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)符合如下方程:

《SciPy-数值计算库》

因此可以如下定义通过(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的积分区间。
半球体积的积分的示意图如下:

《SciPy-数值计算库》 半球体积的双重定积分示意图

X轴的积分区间为-1.0到1.0,对于X=x0时,通过对Y轴的积分计算出切面的面积,因此Y轴的积分区间如图中红色点线所示。

  • 解常微分方程组

scipy.integrate库提供了数值积分和常微分方程组求解算法odeint。
下面让我们来看看如何用odeint计算洛仑兹吸引子的轨迹。洛仑兹吸引子由下面的三个微分方程定义:

《SciPy-数值计算库》

[洛仑兹吸引子的详细介绍]这三个方程定义了三维空间中各个坐标点上的速度矢量。从某个坐标开始沿着速度矢量进行积分,就可以计算出无质量点在此空间中的运动轨迹。其中sigma, rho, beta 为三个常数,不同的参数可以计算出不同的运动轨迹: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()

《SciPy-数值计算库》 用odeint函数对洛仑兹吸引子微分方程进行数值求解所得到的运动轨迹

我们看到即使初始值只相差0.01,两条运动轨迹也是完全不同的。
在程序中先定义一个lorenz函数,它的任务是计算出某个位置的各个方向的微分值,这个计算直接根据洛仑兹吸引子的公式得出。然后调用odeint,对微分方程求解,odeint有许多参数,这里用到的四个参数分别为:

  • lorenz, 它是计算某个位移上的各个方向的速度(位移的微分)
  • (0.0, 1.0, 0.0),位移初始值。计算常微分方程所需的各个变量的初始值
  • t, 表示时间的数组,
  • odeint对于此数组中的每个时间点进行求解,得出所有时间点的位置
    原文作者:羽恒
    原文地址: https://www.jianshu.com/p/be1252b0059e
    本文转自网络文章,转载此文章仅为分享知识,如有侵权,请联系博主进行删除。
点赞