计算物理（六）Python 代数运算与高等数学

科学运算包括数值运算和符号运算，数值运算可以使用 Numpy 库和 Scipy 库。本篇文章主要介绍如何使用 Python 进行基础的数学运算，例如方程求根、数值积分、常微分方程、偏微分方程等。

数值微分

数值微分使用 numdifftools 包中的 Derivative 函数

import numdifftools as nd

Derivative 可返回导数函数，而非离散的导数值

nd.Derivative(fun, step=None, method='central', order=2, n=1)

• step 为差分步长，默认为 infer
• method 为差分方法，默认为中心差分
• order 为误差阶数，默认为 $\omicron(x^2)$
• n 为求导阶数，默认为一阶导数

数值积分

数值积分使用 Scipy 中的 integrate 模块，调用方法为

from scipy import integrate

一重积分

基本语法

一重积分使用 integrate.quad 函数

integrate.quad(func, a, b, args=(tuple), epsabs, epsrel)

• func 可以是单参数函数 $f(x)$ 或多参数函数 $f(x,t_1,t_2,\cdots)$，如果是后者默认只对第一个变量积分

• a 和 b 分别是积分上限和积分下限

• args 是可选参数，若被积函数为多参数函数 $f(x,t_1,t_2,\cdots)$ ，则须按顺序传入除被积变量以外的参数的值

  /* 例：将 $2x+3$在 $[0,1]$上积分*/

  f = lambda x,t1,t2: t1*x+t2
  integrate.quad(f,0,1,args=(2,3))

• epsabs 为绝对误差限，默认值为 1.49e-8

• epsrel 为相对误差限，默认值为 1.49e-8

正常情况下，函数的返回值是一个元组，第 0 项为积分值，第 1 项为误差。

对任一变量积分

若要对多参数函数 $f(x,y,z,\cdots)$ 的非第一个变量 (x) 积分，可以在 func 处使用 lambda 临时定义一个函数。

例：考虑对 $f(x,y)=x(y+1)$ 的两个变量分别积分，容易算得

$$\int_0^1 f(x,1) \, \mathrm{d}x = 1 \qquad \int_0^1 f(1,y) \, \mathrm{d}y = 1.5$$

对于前者，很容易可以想到通过如下代码实现

# The first integral
f = lambda x,y: x*(y+1)
integrate.quad(f,0,1,args=(1,))

对于后者，可以通过

# The second integral, method one
f = lambda x,y: x*(y+1)
integrate.quad(lambda y,x:f(x,y),0,1,args=(1,))

或者定义临时的匿名函数时直接给定 $x=1$ ，避免使用 args

# The second integral, method two
f = lambda x,y: x*(y+1)
integrate.quad(lambda y:f(1,y),0,1)

显然第二种方法的代码可读性更高。

二重积分

二重积分使用 integrate.dblquad 函数

integrate.dblquad(func, a, b, gfun, hfun, args=(), epsabs=1.49e-8, epsrel=1.49e-8)

将计算函数 $f(y,x)$ 的如下积分

$$\int_a^b \mathrm{d} x \int_{g(x)}^{h(x)} \mathrm{d}y \, f(y,x)$$

三重积分

三重积分使用 integrate.tplquad 函数

integrate.tplquad(func, a, b, gfun, hfun, qfun, rfun, args=(), epsabs=1.49e-08, epsrel=1.49e-08)

将计算函数 $f(z,y,x)$ 的如下积分

$$\int_a^b \mathrm{d} x \int_{g(x)}^{h(x)} \mathrm{d}y \int_{q(x,y)}^{r(x,y)} \mathrm{d}z \, f(z,y,x)$$

方程的数值解

方程数值求根使用 Scipy 中的 optimize 模块

from scipy.optimize import root_scalar

标量方程

标量方程使用 root_scalar 求根，默认对第一个变量求根

root_scalar(f, args=(), method=None, bracket=None, fprime=None, fprime2=None, x0=None, x1=None)

• 二分法 ‘bisect’ ：braket = [a,b] (A sequence of 2 floats)
• 牛顿法 ‘newton’ ：fprime, x0
• 哈雷法 ‘halley’ ：fprime, fprime2, x0
• 弦割法 ‘secant’ ：x0, x1

线性向量方程

线性方程求解使用 Scipy 中的 linalg 线性代数模块

一般线性方程

使用 scipy.linalg.solve 求解矩阵方程A @ x == b

scipy.linalg.solve(A, b, assume_a='gen', transposed=False)

说明：assume_a 可指定矩阵类型

matrix type	assume_a
generic matrix	‘gen’
symmetric	‘sym’
hermitian	‘her’
positive definite	‘pos’

若 transposed = True 则求解矩阵方程 A.T @ x == b 。

带状矩阵方程

使用 scipy.linalg.solve_banded 求解带状矩阵方程，使用 scipy.linalg.solveh_banded 求解 Hermite 正定矩阵方程

scipy.linalg.solve_banded(lu, ab, b)

scipy.linalg.solveh_banded(ab, b)

lu 是一个元组，表示主对角线下 (lower) 和主对角线上 (upper) ，分别有几行非零次对角线。ab 是带状矩阵的非零矩阵元，遵循 ab[u + i - j, j] == a[i,j] ，是一个 $l+u+1$ 行 $m$ 列的矩阵。例如对于矩阵方程：

$$A = \begin{bmatrix} 5 & 2 & -1 & 0 & 0 \\ 1 & 4 & 2 & -1 & 0 \\ 0 & 1 & 3 & 2 & -1 \\ 0 & 0 & 1 & 2 & 2 \\ 0 & 0 & 0 & 1 & 1 \\ \end{bmatrix} \quad \boldsymbol{b} = \begin{bmatrix} 0 \\ 1 \\ 2 \\ 3 \\ 3 \end{bmatrix} \quad A\boldsymbol{x}=\boldsymbol{b}$$

ab 矩阵应当设置为

$$[ab] = \begin{bmatrix} 0 & 0 & -1 & -1 & -1 \\ 0 & 2 & 2 & 2 & 2 \\ 5 & 4 & 3 & 2 & 1 \\ 1 & 1 & 1 & 1 & 0 \end{bmatrix}$$

非线性向量方程

非线性向量方程组的求解使用 scipy.optimize.root 函数，语法为

scipy.optimize.root(fun, x0, args=(), method='hybr', jac=None)

求解实数非线性方程使用默认方法 'hybr' 即可，求解复数方程可使用 ’krylov’ 方法。

jac 参数为待求解方程 $\boldsymbol{f}(\boldsymbol{x})=\bm{0}$ 的 Jacobi 矩阵，若传入则显著提升计算速度。

预搜索

预搜索可通过以下代码实现

def bisection(f,a,b,xtol,maxiter):
    iter = 0
    while np.abs(b-a) > xtol and iter <= maxiter:
        iter = iter + 1
        mid = (a+b)/2
        if f(mid)*f(a)<0:
            b = mid
        else:
            a = mid
    root = mid
    dev = np.abs(f(a)-f(b))
    return root,dev
                
def prebisect(func,x0,step,maxdev,xtol=2e-12,maxiter=100):
    while func(x0+step)*func(x0)>0:
        x0 = x0 + step
    [root,dev] = bisection(func,x0,x0+step,xtol,maxiter)
    # If deviation > delta even though |a-b| < epislon, it indicates a leap here
    # When finding a leap, skip it, start new searching from root+step
    if dev > maxdev:
        root = prebisect(func,x0+step,step,maxdev,xtol,maxiter)
    return root

常微分方程的初值问题

简介

一阶非线性常微分方程的普遍形式为

$$\left\{ \begin{matrix} \dfrac{\mathrm{d}\boldsymbol{y}}{\mathrm{d}t} = f(t,\boldsymbol{y}) \\ \boldsymbol{y}(t_0) = \boldsymbol{y}_0 \end{matrix} \right. \implies \boldsymbol{y}=\boldsymbol{y}(t)$$

对于高阶常微分方程，总可以通过降阶化为一阶常微分方程组，例如牛顿第三定律

$$m \frac{\mathrm{d}^2 x}{\mathrm{d} t^2}=F(x, t) \iff \begin{cases} \frac{\mathrm{d} x}{\mathrm{d} t} = \frac{p}{m} \\ \frac{\mathrm{d}p}{\mathrm{d} t} = F(x,t) \end{cases} \iff \dfrac{\mathrm{d}\boldsymbol{y}}{\mathrm{d}t} = f(t,\boldsymbol{y})$$

其中 $\boldsymbol{y}(t)=[x(t),p(t)], \ f(t,\boldsymbol{y}) = [p(t)/m, F(x,t)]$ 。因此，关键在于求解一阶非线性常微分方程。

求解器

常微分方程的初值问题使用 scipy.integrate.solve_ivp ，其中 ivp 代表 initial value problem

scipy.integrate.solve_ivp(fun, t_span, y0, method='RK45', t_eval=None, dense_output=False)

各参数的意义如下

• fun 为导函数 $f(t,\boldsymbol{y})$，t_span 为 A sequence of 2 floats，y0 为初始值
• method 默认为 4/5 阶 Runge-Kutta 法
• t_eval 为时间节点，必须已排序切若不指定则由求解器自动分割 t_span 并求解
• dense_output 选择是否对结果进行插值

返回的结果为一个对象，属性包括

• t ：时间节点
• y ：函数值
• sol ：Odesolution 类