星期六, 4月 27, 2024
最新:

数值分析大巴

专注 Numerical Analysis 数值分析学习, 科研与应用 (NumAnal.com)

编程语言 软件工具 

Python 科学计算: NumPy

numanal 0 个评论 ,

科学计算中包含大量对数组 (包括向量, 矩阵和高维数组等) 的生成, 访问及操作等运算过程, 简洁高效的数组实现和运算库变得极其重要. Numpy 是 Python 语言中专门为数组运算提供支持的库, 它提供了强大且易于使用的数组对象, 因为其底层由 C 语言实现, 克服了原生 Python 语言执行速度慢的缺点, 现在已经处于 Python 科学计算生态中最为核心的位置.

预备知识

Python 的基础语法:

  • 列表的索引从零开始, a[0] 表示引用数组中的第一个元素.
  • 列表切片, a[start : end : step]; 切片的结果不包含结束索引, 即不包含最后的一位.
  • 列表的索引, a[-1] 引用列表的最后一个元素.
  • for 循环和 range 的用法; break 和 continue 的区别.

高等数学中的基本知识, 如向量和矩阵的概念, 基本运算. 在科学计算中, 一个向量可以代表多种含义: 常微分方程的近似解, 有限差分和有限元方法中微分方程的近似解等等; 同样地, 一个矩阵也可以代表多种意义: 一个线性变换, 一个投影算子, 或者有限差分和有限元方法中的刚度矩阵. 也就是说, 向量, 矩阵或者张量在数值分析中需要大量使用, numpy 可以很好的满足这种需要.

入门

数组是 Numpy 中的核心数据结构. 导入 numpy 库:

import numpy as np

然后就可以通过 np.xxx 来使用 numpy 中的数组, 函数和其它对象. 例如, numpy 内置了如下常数

  • 数圆周率 np.pi
  • 自然对数 np.e
  • 无穷大浮点数 np.inf
  • None 的别名 np.newaxis. 可用于将一个一维数组变成一个只有一列的二维数组:
x = np.arange(3)
# array([0, 1, 2])
x[:, newaxis]
# array([[0],
# [1],
# [2]])

数组生成

生成整数向量和实数向量:

a = np.array([1, 2, 3, 4]) # or np.array([1, 2, 3, 4], dtype=int)
b = np.array([1., 2., 3., 4.]) # or np.array([1., 2., 3., 4.], dtype=float)
print(a, a.dtype, b, b.dtype)
# [1 2 3 4] int64 [1. 2. 3. 4.] float64

生成数组时如果没有指定元素类型, 那么 numpy 会默认生成 64 位整数或者实数 (int64 或 float64). 另外, 还可以生成整数或实数高维数组:

# 矩阵
a = np.array([[1, 2], [3, 4]])
# a = np.array([[1, 2], [3, 4]], dtype=int) # same
b = np.array([[1., 2.], [3., 4.]]) 
# b = np.array([[1., 2.], [3., 4.]], dtype=float) # same
print(a, a.dtype, b, b.dtype)

c = np.array([ [[0,  1,  2],  # a 3D integer array, integer
    [ 10, 12, 13]],
    [[100, 101, 102],
    [110, 112, 113]] ])
d = np.arange(24, dtype=float).reshape(2, 3, 4) # a 3d float array
# array([[[ 0,  1,  2,  3],
#         [ 4,  5,  6,  7],
#         [ 8,  9, 10, 11]],

#        [[12, 13, 14, 15],
#         [16, 17, 18, 19],
#         [20, 21, 22, 23]]])

上面创建的数组对象包含属性和方法, 例如常用的属性有 dtype, shape, ndim, size 等, 常用的方法有 T, reshape, ravel 等. 例如:

# methods of array
M = np.arange(12).reshape(3,4)
# [[ 0  1  2  3]
#  [ 4  5  6  7]
#  [ 8  9 10 11]]
print("dtype=", M.dtype, ", shape=", M.shape, ", ndim=", M.ndim, ", size=", M.size)
# dtype= int64 , shape= (3, 4) , ndim= 2 , size= 12

M.ravel()  # returns the array, flattened
# array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11])

N = M.reshape(6, 2)  # returns the array with a modified shape
# array([[ 0,  1],
#        [ 2,  3],
#        [ 4,  5],
#        [ 6,  7],
#        [ 8,  9],
#        [10, 11]])
N.T  # returns the array, transposed
# array([[ 0,  2,  4,  6,  8, 10],
#        [ 1,  3,  5,  7,  9, 11]])
N.T.shape
# (2, 6)
N.shape
# (6, 2)

常用生成函数. Numpy 内置了众多快速生成简单矩阵的函数, 例如:

  • np.zeros: 零矩阵
  • np.ones: 全一矩阵
  • np.random: 随机矩阵
  • np.arange: 功能相当于 python 中的 range
  • np.linspace:
np.zeros((3, 4))
# array([[0., 0., 0., 0.],
#        [0., 0., 0., 0.],
#        [0., 0., 0., 0.]])
np.ones((2, 3, 4), dtype=np.int32)
# array([[[1, 1, 1, 1],
#         [1, 1, 1, 1],
#         [1, 1, 1, 1]],

#        [[1, 1, 1, 1],
#         [1, 1, 1, 1],
#         [1, 1, 1, 1]]], dtype=int32)
np.arange(10, 30, 5)     # start:end:step
# array([10, 15, 20, 25])
np.arange(0, 2, 0.3)     # generate float array
# array([0. , 0.3, 0.6, 0.9, 1.2, 1.5, 1.8])
np.linspace(0, 2, 9)     # 9 numbers from 0 to 2
# array([0.  , 0.25, 0.5 , 0.75, 1.  , 1.25, 1.5 , 1.75, 2.  ])

数组索引

数组中的元素或部分通过索引或切片的方式来访问. Numpy 中的数组索引与 C/C++ 相同, 从 0 开始.

向量: 创建向量 a = np.array([1, 2, 3, 4]) # [1 2 3 4]. 经常用到的索引语法如下:

语句作用输出
a引用数组 a[1 2 3 4]
a[0]访问第一个元素1
a[-1]访问最后一个元素4
a[:]等同于 a[1 2 3 4]
a[0:2]子向量 (注意不包括 a[2])[1 2]
a[0:-1]子向量 (注意不包括 a[-1])[1 2 3]

矩阵: 创建矩阵数组

M = M = np.arange(12).reshape(3,4) 
# [[ 0  1  2  3]
#  [ 4  5  6  7]
#  [ 8  9 10 11]]

经常用到的关于矩阵切片的语法如下:

M[0,0]
# 0, 第一行第一列的元素
M[-1,-1]
# 11, 最后一行最后一列的元素
M[0,:]
# array([0, 1, 2, 3]), 第一行向量
M[:, 1]
# array([1, 5, 9]), 第二列向量
M[:, :]
# 矩阵本身, array([[ 0,  1,  2,  3],
#               [ 4,  5,  6,  7],
#               [ 8,  9, 10, 11]])
M[:, 0:4:2]
# 特定的列 (start:end:step)) 构成的数组, 
# array([[ 0,  2],
#        [ 4,  6],
#        [ 8, 10]])
M[0:4:2, 0:4:2]
# 子矩阵
# array([[ 0,  2],
#        [ 8, 10]])
M[[0, 2], [0, 2]]
# 注意, 这样得到的并不是子矩阵! 而是向量 [M[0,0], M[2,2]]
# 不能使用两个向量来方程子矩阵
# array([ 0, 10])
M[np.ix_([0, 2], [0, 2])]
# 子矩阵
# array([[ 0,  2],
#        [ 8, 10]])

注. 有限元中经常遇到的操作就是将局部刚度矩阵添加到全局刚度矩阵中, 这也意味着需要将局部基函数对应的编号对应到全局编号, 从而正确的整合刚度矩阵. 这个过程可以使用矩阵索引方法 M[np.ix_([0, 2], [0, 2])] 来实现. 另外, 在处理 Dirichlet 边界条件时也会用到这个语法.

数组运算

相同维度的数组之间做初等运算 (+, -, *, /, +=, -=, *=, /=) 等同于各个元素做相应的运算. 需要注意的是, 如果两个二维数组 (矩阵) 做乘法 (*) 运算, 并不是线性代数中的矩阵乘法. 矩阵乘法运算符号为 @ 或者利用函数 np.dot().

初等函数. Numpy 内置了初等函数, 如 np.sin, np.cos, np.tan, np.exp, np.log 等, 当使用这些函数作用在数组上时, 就是分别对数组中每个元素进行作用, 即广播.

# 初等函数
v = np.arange(4)*np.pi/4
# [0, pi/4, pi/2, 3pi/4]
np.sin(v)
# array([0.        , 0.70710678, 1.        , 0.70710678])
np.tan(v)
array([ 0.00000000e+00,  1.00000000e+00,  1.63312394e+16, -1.00000000e+00])

向量与向量. 向量之间可以进行点积 (dot product)、内积 (inner product)、外积 (outer product)、叉乘 (cross) 等运算. 假设 $$ a=\begin{bmatrix} a_0\\ a_1\\ \vdots \\ a_n \end{bmatrix}, b=\begin{bmatrix} b_0\\ b_1\\ \vdots \\ b_n \end{bmatrix}, c=\begin{bmatrix} c_0\\ c_1\\ \vdots \\ c_m \end{bmatrix} $$ 那么

  • a*b = \([a_0b_0, a_1b_1, \dots, a_nb_n]\) (逐元素相乘).
  • np.dot(a,b) = np.vdot(a,b) = np.inner(a,b) = \(a\cdot b=\sum_{i=1}^n a_ib_i\).
  • np.outer(a,c) = $$ \begin{bmatrix} a_0c_0 & a_0c_1 & \cdots & a_0c_m \\ a_1c_0 & a_1c_1 & \cdots & a_1c_m \\ \vdots & \vdots & \ddots & \vdots \\ a_nc_0 & a_nc_1 & \cdots & a_nc_m \end{bmatrix}.$$
  • np.cross(x,y) = $$x\times y = \begin{bmatrix} x_2y_3 – x_3y_2 \\ x_3y_1 – x_1y_3 \\ x_1y_2 – x_2y_1 \end{bmatrix},\quad \text{其中 } x = \begin{bmatrix} x_1\\ x_2\\ x_3 \end{bmatrix}, y = \begin{bmatrix} y_1\\ y_2\\ y_3 \end{bmatrix}.$$

矩阵与向量. 矩阵 \(A\in \mathbb{R}^{n\times n}, x \in \mathbb{R}^n\), 则

  • A @ x = np.matmul(A,x) = np.dot(A,x) 就是线性代数的矩阵乘法.
  • 向量与矩阵: x @ A = np.dot(x,A).
  • 矩阵与矩阵: A @ A = np.dot(A,A) = np.matmul(A,A).

Kronecker 积. \(A\in \mathbb{R}^{m\times n}, B\in \mathbb{R}^{p\times q}\), 那么 np.kron(A,B) = $$ A\otimes B = \begin{bmatrix} a_{11}B & \cdots & a_{1n}B \\ \vdots & \ddots & \vdots \\ a_{m1}B & \cdots & a_{mn}B \end{bmatrix}\in \mathbb{R}^{pm\times nq}.$$

进阶

更多的数组操作

  • 数组求和 np.sum.
A = np.arange(16).reshape(4,4)
# [[ 0  1  2  3]
#  [ 4  5  6  7]
#  [ 8  9 10 11]
#  [12 13 14 15]]
np.sum(A)
# 120
np.sum(A,0)
# array([24, 28, 32, 36])
np.sum(A,1)
# array([ 6, 22, 38, 54])
  • 数组重复 np.repeat, np.tile.
b = np.array([[1, 2], [3, 4]])
np.tile(b, (1, 3))  # repeat array (1,3)
# array([[1, 2, 1, 2, 1, 2],
#        [3, 4, 3, 4, 3, 4]])
np.tile(b, (2, 1)) # repeat array (2,1)
# array([[1, 2],
#        [3, 4],
#        [1, 2],
#        [3, 4]])
  • 数组取整 np.floor, np.ceil, np.round, np.trunc.
  • 数组极值 np.maximum, np.mininum, np.max.
  • 数组逆序 np.flip.
  • 数组增删 np.delete, np.insert, np.append.

线性代数运算

线性代数中更多的运算, 如向量和矩阵范数, 矩阵分解, 线性方程组求解, 特征值和特征向量等等都由线性代数子库 np.linalg (linear algebra) 中的函数提供的.

  • 范数 np.linalg.norm(a)
  • 条件数 np.linalg.cond(a)
  • 行列式 np.linalg.det(a)
  • 矩阵的秩 np.linalg.matrix_rank(a)
  • 矩阵的迹 np.linalg.trace(a)
  • Cholesky 分解 linalg.cholesky(a)
  • QR 分解 linalg.qr(a)
  • SVD 分解 linalg.svd(a)
  • 特征值和特征向量 linalg.eig(a)
  • 线性方程组 linalg.solve(A, b)
  • 矩阵的逆 linalg.inv(A)
  • 矩阵的最小二乘解 linalg.lstsq(A, b)

总结

Numpy 提供了完备的矩阵运算对象和函数, 内容几乎囊括了所有的数值代数中的常见算法. 更多的内容, 如矩阵的 LU 分解, 稀疏矩阵的建立和应用, 广义特征值问题等等更多复杂的问题和算法可以在 Python 的科学计算库 SciPy 中找到.

参考文献

  1. NumPy 官网
  2. NumPy 快速入门 (NumPy quickstart)
  3. Numpy 和 Matlab 语法对照表 (NumPy for Matlab users).
  4. C. Harris, etc. (2020). Array programming with NumPy. Nature Review, 585.
  5. Sympy 使用简介.

分享到:

你也可能喜欢

DistMesh – 一个简单的 MATLAB 网格生成工具

numanal 0

C语言的GCC编译过程和静态库的建立

numanal 0

为什么使用 Python 做科学计算

numanal 0
5 1 投票
文章评分
订阅评论
提醒
guest

0 评论
内联反馈
查看所有评论