scipy 模块参考书

SciPy Tutorial

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

scipy 主要有15个子模块

from scipy import cluster, constants, fftpack, integrate, interpolate, io, linalg, ndimage, odr, optimize, signal, sparse, spatial, special, stats

常量 Constants(scipy.constants)

主要包含数学常数\(\pi\)(pi)以及一些物理常数


基本函数

np.r_ np.c_ 拼接

>>> np.r_[3,[0]*5,-1:1:10j]
array([ 3.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        , -1.        , -0.77777778, -0.55555556, -0.33333333,
       -0.11111111,  0.11111111,  0.33333333,  0.55555556,  0.77777778,
        1.        ])

将切片对象转换为沿第一轴的连接。(横向拼接)

>>> np.c_[np.array([1,2,3]), np.array([4,5,6])]
array([[1, 4],
       [2, 5],
       [3, 6]])

将切片对象转换为沿第二轴的连接。(竖向拼接)
复数10j作为切片语​​法中的步长。这种非标准用法允许将数字解释为在范围内产生的点数而不是步长

np.mgrid 网格生成

>>> np.mgrid[0:5:4j,0:5:4j]
array([[[ 0.    ,  0.    ,  0.    ,  0.    ],
        [ 1.6667,  1.6667,  1.6667,  1.6667],
        [ 3.3333,  3.3333,  3.3333,  3.3333],
        [ 5.    ,  5.    ,  5.    ,  5.    ]],
       [[ 0.    ,  1.6667,  3.3333,  5.    ],
        [ 0.    ,  1.6667,  3.3333,  5.    ],
        [ 0.    ,  1.6667,  3.3333,  5.    ],
        [ 0.    ,  1.6667,  3.3333,  5.    ]]])

积分 (scipy.integrate)

此模块可用用于各种包括反常积分在内的定积分的求解
\(I(a, b)=\int_{0}^{1} a x^{2}+b d x\)

>>> from scipy.integrate import quad
>>> def integrand(x, a, b):
...     return a*x**2 + b
...
>>> a = 2
>>> b = 1
>>> I = quad(integrand, 0, 1, args=(a,b))
>>> I
(1.6666666666666667, 1.8503717077085944e-14)

常微分方程求解 scipy.integrate.ode
BerkeleyMadonna 推荐使用,常微分方程求解器

线性代数相关

from scipy import linalg

范数

向量的范数
\(\|\mathbf{x}\|=\left\{\begin{array}{cl}{\max \left|x_{i}\right|} & {\text { ord }=\inf } \\ {\min \left|x_{i}\right|} & {\text { ord }=-\inf } \\ {\left(\sum_{i}\left|x_{i}\right|^{\text {ord }}\right)^{1 / \operatorname{ord}}} & { | \text { ord } |<\infty}\end{array}\right.\)
矩阵的范数
\(\|\mathbf{A}\|=\left\{\begin{aligned} \max _{i} \sum_{j}\left|a_{i j}\right| & \text { ord }=\text { inf } \\ \min _{i} \sum_{j}\left|a_{i j}\right| & \text { ord }=-\inf \\ \max _{j} \sum_{i}\left|a_{i j}\right| & \text { ord }=1 \\ \min _{j} \sum_{i}\left|a_{i j}\right| & \text { ord }=-1 \\ \max \sigma_{i} & \text { ord }=2 \\ \min \sigma_{i} & \text { ord }=-2 \\ \sqrt{\operatorname{trace}\left(\mathbf{A}^{H} \mathbf{A}\right)} & \text { ord }=\text { 'fro }^{\prime} \end{aligned}\right.\)
例:

>>> import numpy as np
>>> from scipy import linalg
>>> A=np.array([[1,2],[3,4]])
>>> A
array([[1, 2],
      [3, 4]])
>>> linalg.norm(A)
5.4772255750516612
>>> linalg.norm(A,'fro') # frobenius norm is the default
5.4772255750516612
>>> linalg.norm(A,1) # L1 norm (max column sum)
6
>>> linalg.norm(A,-1)
4
>>> linalg.norm(A,np.inf) # L inf norm (max row sum)
7

特征值与特征向量

>>> import numpy as np
>>> from scipy import linalg
>>> A = np.array([[1, 2], [3, 4]])
>>> la, v = linalg.eig(A)
>>> l1, l2 = la
>>> print(l1, l2)   # eigenvalues
(-0.3722813232690143+0j) (5.372281323269014+0j)
>>> print(v[:, 0])   # first eigenvector
[-0.82456484  0.56576746]
>>> print(v[:, 1])   # second eigenvector
[-0.41597356 -0.90937671]

LU分解

\(A=PLU\)

>>> from scipy.linalg import lu
>>> A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])
>>> p, l, u = lu(A)

追赶法

\(a_{i} x_{i-1}+b_{i} x_{i}+c_{i} x_{i+1}=d_{i}\)
where \(a_{1}=0\) and \(c_{n}=0\)
\(\left[\begin{array}{ccccc}{b_{1}} & {c_{1}} & {} & {} & {0} \\ {a_{2}} & {b_{2}} & {c_{2}} & {} & {} \\ {} & {a_{3}} & {b_{3}} & {\ddots} & {} \\ {} & {} & {\ddots} & {\ddots} & {c_{n-1}} \\ {0} & {} & {} & {a_{n}} & {b_{n}}\end{array}\right]\left[\begin{array}{c}{x_{1}} \\ {x_{2}} \\ {x_{3}} \\ {\vdots} \\ {x_{n}}\end{array}\right]=\left[\begin{array}{c}{d_{1}} \\ {d_{2}} \\ {d_{3}} \\ {\vdots} \\ {d_{n}}\end{array}\right]\)

import numpy as np

## Tri Diagonal Matrix Algorithm(a.k.a Thomas algorithm) solver
def TDMAsolver(a, b, c, d):
    '''
    TDMA solver, a b c d can be NumPy array type or Python list type.
    refer to http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm
    and to http://www.cfd-online.com/Wiki/Tridiagonal_matrix_algorithm_-_TDMA_(Thomas_algorithm)
    '''
    nf = len(d) # number of equations
    ac, bc, cc, dc = map(np.array, (a, b, c, d)) # copy arrays
    for it in range(1, nf):
        mc = ac[it-1]/bc[it-1]
        bc[it] = bc[it] - mc*cc[it-1] 
        dc[it] = dc[it] - mc*dc[it-1]
                
    xc = bc
    xc[-1] = dc[-1]/bc[-1]

    for il in range(nf-2, -1, -1):
        xc[il] = (dc[il]-cc[il]*xc[il+1])/bc[il]

    return xc

例:

A = np.array([[10,2,0,0],[3,10,4,0],[0,1,7,5],[0,0,3,4]],dtype=float)   

a = np.array([3.,1,3]) 
b = np.array([10.,10.,7.,4.])
c = np.array([2.,4.,5.])
d = np.array([3,4,5,6.])

print(TDMAsolver(a, b, c, d))
>> [ 0.14877589  0.75612053 -1.00188324  2.25141243]
#compare against numpy linear algebra library
print(np.linalg.solve(A, d))
>> [ 0.14877589  0.75612053 -1.00188324  2.25141243]

对角线元素

numpy.diag(v,k=0)

如果v是一个二维的矩阵,则返回对角线元素组成的向量
如果v是一个一维的向量,则返回对角线为该向量的矩阵
k表示偏移量,正数为右移,负数为左移