Python科学计算之Scipy

随笔 2017-01-18

起步

SciPy 库建立在 Numpy 库之上,提供了大量科学算法,例如线性代数、常微分方程数值求解、信号处理、图像处理、稀疏矩阵等等,主要包括这些主题:

  • 特殊函数 (scipy.special)
  • 积分 (scipy.integrate)
  • 最优化 (scipy.optimize)
  • 插值 (scipy.interpolate)
  • 傅立叶变换 (scipy.fftpack)
  • 信号处理 (scipy.signal)
  • 线性代数 (scipy.linalg)
  • 稀疏特征值 (scipy.sparse)
  • 统计 (scipy.stats)
  • 多维图像处理 (scipy.ndimage)
  • 文件 IO (scipy.io)

特殊函数

函数列表可参考: http://docs.scipy.org/doc/scipy/reference/special.html#module-scipy.special

积分

数值求积,Scipy提供了一些列不同类型的求积函数,像是 quad, dblquad 还有 tplquad 分别对应单积分,双重积分,三重积分。

wm.png

假设 f(x) = x^2 要计算[1, 2]区间的面积,不引用scipy我们是这样计算的:

# coding: utf-8
import numpy as np
def f(x):
    return x * x

N = 10000

# 原生方法
total = 0
for i in range(1, 10001):
    total += f(2 + i / 10000.0) * (1 / 10000.0)
print total
# output: 6.333583335

# 使用numpy
x = np.linspace(2, 3, N)
dx = 1.0 / N
y = f(x)

print dx * np.sum(y)
# output: 6.33335000167

如果使用scipy:

from scipy.integrate import quad
x_lower = 2 # the lower limit of x
x_upper = 3 # the upper limit of x

val, abserr = quad(f, x_lower, x_upper)
print "integral value =", val, ", absolute error =", abserr
# output: integral value = 6.33333333333 , absolute error = 7.03141248929e-14

微分

SciPy 提供了两种方式来求解常微分方程:基于函数 odeint 的API与基于 ode 类的面相对象的API。通常 odeint 更好上手一些,而 ode 类更灵活一些。

这里我们将使用 odeint 函数,首先让我们载入它:

from scipy.integrate import odeint, ode

常微分方程组的标准形式如下:

wm (1).png

wm (2).png

为了求解常微分方程我们需要知道方程 f 与初始条件y(0) 注意到高阶常微分方程常常写成引入新的变量作为中间导数的形式。 一旦我们定义了函数 f 与数组 y_0 我们可以使用 odeint 函数:

y_t = odeint(f, y_0, t)

傅里叶变换

傅立叶变换是计算物理学所用到的通用工具之一。Scipy 提供了使用 NetLib FFTPACK 库的接口,它是用FORTRAN写的。Scipy 还另外提供了很多便捷的函数。不过大致上接口都与 NetLib 的接口差不多。

fft.gif

Numpy也有一个FFT实现(numpy.fft)。然而,通常scipy的应该优先使用,因为它使用了更有效率的底层实现。

参考:https://docs.scipy.org/doc/scipy/reference/tutorial/fftpack.html

线性方程

线性代数模块包含了大量矩阵相关的函数,包括线性方程求解,特征值求解,矩阵函数,分解函数(SVD, LU, cholesky)等等。 详情见:http://docs.scipy.org/doc/scipy/reference/linalg.html

线性方程组

线性方程组的矩阵形式:

wm (3).png

A是矩阵,x,b 是向量:

# coding: utf-8
import numpy as np
from scipy.linalg import solve

A = np.array([[1,2,3], [4,5,6], [7,8,9]])
b = np.array([1,2,3])

x = solve(A, b)
print x
# output: [-0.33333333,  0.66666667,  0.        ])

# check
print np.dot(A, x)
# output : [ 1.  2.  3.]

最优化

最优化 (找到函数的最大值或最小值) 问题是数学中比较大的话题, 复杂的函数与变量的增加会使问题变得更加困难。这里我们只看一些简单的例子。

想知道更多,详情见: http://scipy-lectures.github.com/advanced/mathematical_optimization/index.html

更多


本文由 hongweipeng 创作,采用 知识共享署名 3.0,可自由转载、引用,但需署名作者且注明文章出处。

赏个馒头吧