Python科学计算学习

Xin Qiu Nov 07, 2015

之前参加数学建模比赛都是用的 MatLab，然而电脑每次运行 MatLab 都烫的要命，所以决定用自己比较熟悉的 python 来进行科学计算。

NumPy

NumPy提供了两种基本的对象：ndarray（N-dimensional array object）和 ufunc（universal function object）。ndarray(下文统一称之为数组)是存储单一数据类型的多维数组，而ufunc则是能够对数组进行处理的函数。

推荐使用以下方法导入NumPy 函数库：

1	import numpy as np

ndarray对象

创建

通过给array函数传递Python的序列对象创建数组，如果传递的是多层嵌套的序列，将创建多维数组：

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

上面的例子是在 python 交互界面中产生的结果，在 PyCharm 中运行以下代码:

import numpy as np

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

print A
print B

显示的是：

[1 2 3 4]
[[ 1  2  3  4]
 [ 4  5  6  7]
 [ 7  8  9 10]]

数组的shape属性返回的是数组的维度的 Tuple 。如果是以为数组，就省略了维度1，直接显示元素个数，如果是二维数组，返回的 Tuple 第一个元素代表二维数组是几行，第二个元素代表二维数组是几列。多维数组亦可以表示。

>>> A.shape
(4,)
>>> B.shape
(3, 4)
>>>

还可以使用数组的 reshape 方法，创建一个改变了尺寸的新数组。

>>> C = A.reshape((2,2))
>>> C
array([[1, 2],
       [3, 4]])
>>>

数组的 dtype 属性代表数组的元素类型，可以在创建数组的时候，使用 dtype 参数来改变元素类型。float代表浮点数，complex代表复数。

>>> np.array([[1, 2, 3, 4],[4, 5, 6, 7], [7, 8, 9, 10]], dtype=np.complex)
array([[  1.+0.j,   2.+0.j,   3.+0.j,   4.+0.j],
       [  4.+0.j,   5.+0.j,   6.+0.j,   7.+0.j],
       [  7.+0.j,   8.+0.j,   9.+0.j,  10.+0.j]])
>>>

NumPy提供了很多专门用来创建数组的函数。

arange 函数类似于python的range函数，通过指定开始值、终值和步长来创建一维数组，注意数组不包括终值:
1
2
>>> np.arange(0,1,0.1)
array([ 0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9])

linspace函数通过指定开始值、终值和元素个数来创建一维数组，可以通过endpoint关键字指定是否包括终值，缺省设置是包括终值:

 np.linspace(0, 1, 12)
array([ 0.        ,  0.09090909,  0.18181818,  0.27272727,  0.36363636,
        0.45454545,  0.54545455,  0.63636364,  0.72727273,  0.81818182,
        0.90909091,  1.        ])

logspace函数和linspace类似，不过它创建等比数列，下面的例子产生1(10^0)到100(10^2)、有20个元素的等比数列:

	>>> np.logspace(0, 2, 20)
array([   1.        ,    1.27427499,    1.62377674,    2.06913808,
          2.6366509 ,    3.35981829,    4.2813324 ,    5.45559478,
          6.95192796,    8.8586679 ,   11.28837892,   14.38449888,
         18.32980711,   23.35721469,   29.76351442,   37.92690191,
         48.32930239,   61.58482111,   78.47599704,  100.        ])

ufunc运算

ufunc 是 universal function 的缩写，它是一种能对数组的每个元素进行操作的函数。

frompyfunc 是将任意的 Python 函数转为 Numpy ufunc。frompyfunc 的调用格式为 frompyfunc(func, nin, nout)，其中 func 是计算单个元素的函数，nin 是此函数的输入参数的个数，nout 是此函数的返回值的个数。

reduce 方法和Python的reduce函数类似。

矩阵运算

使用matrix类创建的是矩阵对象，它们的加减乘除运算缺省采用矩阵方式计算，用法和matlab十分类似。

dot : 对于两个一维的数组，计算的是这两个数组对应下标元素的乘积和(数学上称之为内积)；对于二维数组，计算的是两个数组的矩阵乘积；对于多维数组，它的通用计算公式如下，即结果数组中的每个元素都是：数组a的最后一维上的所有元素与数组b的倒数第二位上的所有元素的乘积和：
1
dot(a, b)[i,j,k,m] = sum(a[i,j,:] * b[k,:,m])
inner : 和dot乘积一样，对于两个一维数组，计算的是这两个数组对应下标元素的乘积和；对于多维数组，它计算的结果数组中的每个元素都是：数组a和b的最后一维的内积，因此数组a和b的最后一维的长度必须相同：
1
inner(a, b)[i,j,k,m] = sum(a[i,j,:]*b[k,m,:])
outer : 只按照一维数组进行计算，如果传入参数是多维数组，则先将此数组展平为一维数组之后再进行运算。outer乘积计算的列向量和行向量的矩阵乘积。

SciPy

SciPy函数库在NumPy库的基础上增加了众多的数学、科学以及工程计算中常用的库函数。例如线性代数、常微分方程数值求解、信号处理、图像处理、稀疏矩阵等等。由于其涉及的领域众多、本书没有能力对其一一的进行介绍。作为入门介绍，让我们看看如何用SciPy进行插值处理、信号滤波以及用C语言加速计算。

最小二乘拟合

今年的数学建模国赛——太阳影子长度中，就有一问需要用到最小二乘法拟合，可惜当初不是很会。最小二乘法就是假设有一组实验数据(x[i], y[i])，我们知道它们之间的函数关系:y = f(x)，通过这些已知信息，需要确定函数中的一些参数项。例如，如果f是一个线型函数f(x) = k*x+b，那么参数k和b就是我们需要确定的值。如果将这些参数用 p 表示的话，那么我们就是要找到一组 p 值使得如下公式中的S函数最小:

$S(p)=\sum_{i=1}^m[y_{i}-f(x_{i},P)]^{2}$
这种算法被称之为最小二乘拟合(Least-square fitting)。

scipy中的子函数库optimize已经提供了实现最小二乘拟合算法的函数leastsq。下面是用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()

通过leastsq函数对带噪声的实验数据x, y1进行数据拟合，可以找到x和真实数据y0之间的正弦关系的三个参数： A, k, theta。

函数最小值

optimize库提供了几个求函数最小值的算法：fmin, fmin_powell, fmin_cg, fmin_bfgs。就只举一个小例子：

from scipy.optimize import fmin
def myfunc(x):
    return x**2-4*x+8
x0=[2]
xopt=fmin(myfunc,x0)
print xopt

返回结果：

Optimization terminated successfully.
         Current function value: 4.000000
         Iterations: 11
         Function evaluations: 22
[ 2.]

其中 x0 为初始估计值。

SymPy

SymPy是Python的数学符号计算库，用它可以进行数学公式的符号推导。常用的导入方法是：

1	from sympy import *

符号变量

在SymPy中，数学符号是Symbol类的对象.

1	x = Symbol('x')

此时， x 是一个数学符号。

I 在 Sympy 中表示虚数单位

有理数

SymPy 中有三种不同的数值类型：Real，Rational，Integer：

1	r1 = Rational(4,5)

$\frac{4}{5}$

r1就是一个分数。

数值计算

Sumpy 以任意精度的数值库作为后端，预定义了一些数学常量的表达式，如 pi, e, 以及 oo（代表极限）。

为了计算表达式的数值，可以使用 evalf 函数 (或者 N 函数)。

显示pi的50位。n参数决定有效数位。

1	pi.evalf(n=50)

subs函数可以将算式中的符号进行替换:

expression.subs(x, y) : 将算式中的x替换成y
expression.subs({x:y,u:v}) : 使用字典进行多次替换
expression.subs([(x,y),(u,v)]) : 使用列表进行多次替换

1	y = (x + pi)**2

$y=(x+\pi )^{2}$

以上式子中的 x 若要替换成其他值，可以这么写：

1	y.subs(x, 1.5)

代数运算

展开与分解

expand 函数可以将公式展开。

1	expand((x+1)(x+2)(x+3))

$x^{3}+6 x^{2}+11x+6$

expand 有一系列参数能够决定式子的展开形式，trig=True 能展开三角公式：

1	expand(sin(a+b), trig=True)

$sin(a+b)=sin(a)cos(b)+cos(a)sin(b)$

化简

simplify 进行化简操作，特殊的化简还有 trigsimp， powsimp， logcombine 等等。

分式分解与分式展开

1	f1 = 1/((a+1)*(a+2))

$-\frac{1}{a+2} +\frac{1}{a+1}$

微积分

微分

diff 函数用来做微分, $y=(x+\pi )^{2}$ ,对 y 做 x 的微分，

1	diff(y, x)

1	diff(y, x, 2)

这个表示对 y 做 x 的二次微分，第三个参数表示几重微分。

计算多变量式子的多重微分可以这么做：
$\frac{d^{3}f}{dxdy^{2} }$

1	diff(f, x, 1, y, 2)

表示先对 f 做 x 的一重微分，再做 y 的二重微分。

积分

使用 integrate 函数求积分。

1	integrate(f, x)

表示对 f 做 x 的积分。

1	integrate(f, (x, -1, 1))

表示对 f 做 x 的积分，其中 x 属于 -1 到 1。

总和与连乘积

求总和：

$\sum_{n=1}^{10} \frac{1}{n^{2}}$

1 2	n = Symbol("n") Sum(1/n**2, (n, 1, 10))

求连乘积：

$\prod_{n=1}^{10} n$

1	Product(n, (n, 1, 10))

极限

计算极限：

$\lim_{x \rightarrow 0}{\frac{sin(x)}{x} }$

1	limit(sin(x)/x, x, 0)

数列

数列展开的用法：(默认从 x=0 展开)

$1 + x + \frac{x^{2}}{2} + \frac{x^{3}}{6} + \frac{x^{4}}{24} + \frac{x^{5}}{120} + \mathcal{O}\left(x^{6}\right)$

1	series(exp(x), x)

解方程

solve 能够解方程与方程组:

1	solve(x**2 - 1, x)

解方程组：

1	solve([x + y - 1, x - y - 1], [x,y])

参考资料

用Python做科学计算

Numpy and Scipy Documentation