Python数值微积分，摆脱被高数支配的恐惧

Python科学计算：数组💯数据生成

差分和累加

微积分是现代科学最基础的数学工具，但其应用对象往往是连续函数，而其在非连续函数的类比，便是差分与累加。在【numpy】中，可通过【diff】和【cumsum】来完成这两项任务。

以$y=\sin 2x$为例，其导数为$\frac{\text{d}y}{\text{d}x}=2\cos x$，积分则为$\int y\text{d}x=-\frac{1}{2}\cos 2x+C$，$C$是某个常数。这三个函数的曲线分别为

绘图函数如下

import matplotlib.pyplot as plt
import numpy as np
dx = 0.1
x = np.arange(100)*dx
y = np.sin(2*x)
plt.plot(x, y, label="y=sin(2x)")
plt.plot(x[1:], np.diff(y)/dx, label="diff(y)/dx")
plt.plot(x, np.cumsum(y)*dx, label="cumsum(y)*dx")plt.grid()
plt.legend()
plt.show()

其中，diff用于求差分，其输入参数除了待差分数组之外，还有n和axis，比较常用，n为差分的阶数，默认为1；axis用于高维数组中，表示计算的方向，默认-1表示最后一个轴。

cumsum用于累加，对于输入数组$y$，其返回数组为$S$，则$S_n={\sum }_{i=0}^n{y}_i$。

无论diff还是cumsum，均只针对输入数组进行操作，而不会考虑微积分中至关重要的$\text{d}x$，所以绘图时对这一部分进行了补全。

此外，由于差分的实质是后一个减去前一个，所以元素个数必然会减少，所以在绘图时，令$x$从1开始。这是一个在编程时很容易出错的地方，故而numpy还提供了另一个函数【ediff1d】，这是一个只做一阶差分计算的函数，但提供了to_end和to_begin参数，分别用于在diff计算结果的后面或前面补充数值。

积分

积分一开始被引入教材，是以梯形求和为示例的：将函数$y=f(x)$无限分割，然后对相邻两点取平均，再乘以$\text{d}x$之后进行求和，即${\lim }_{{\delta }_x\to 0}\sum \frac{y_i+y_{i+1}}{2}{\delta }_x$。

【trapz】可实现上述过程，但要求$y$是一个给定的数组，且${\delta }_x$为1。很显然，这个过程只能称之为梯形求和，毕竟积分的要求是${\delta }_x\to 0$，$1$和$0$有着本质的区别。

为此，【scipy.intergrate】作为顾名思义的积分模块，提供了真真正正的积分。为了行文简洁，后文将此模块简称为【si】模块。

【quad】是【si】中最常用的积分函数，以函数$x^2$和$\sin x$为例，其使用流程如下

import numpy as np
from scipy.integrate import quadfunc = lambda x: x**2
quad(func, 0, 4)        # (21.33, 2.37-13)
quad(np.sin, 0, np.pi)  # (2.0, 2.22e-14)

其中，quad共输入了三个参数，分别是待积分函数、积分下界与积分上界，其返回值有二，分别为积分结果和计算误差。

这两个测试函数的解析形式如下，可见计算结果基本温和。

$$\begin{gathered} {\int }_0^4{x}^2\text{d}x=\frac{1}{3}{x}^3{|}_0^4=\frac{64}{3}\approx 21.3 \\ {\int }_0^{\pi }\sin x\text{d}x=-\cos x{|}_0^{\pi }=2 \end{gathered}$$

除了三个必须输入的参数之外，下列参数也较为常用

• args为func函数中，除待求积分参数之外的其他参数，默认为空
• epsabs, epsrel 分别为绝对和相对误差，默认为$1.49\times 10^{-8}$
• limit 自适应算法中子区间的个数，默认50
• points 断点位置，默认为None
• weight, wvar 定义域区间内的权重类型和权重，默认为None
• wopts, maxp1 切比雪夫矩及其上限，默认为None和50
• full_output=0, limlst=50, complex_func=False

其中，weight和wvar参数的具体取值如下。

weight	wvar	函数
“cos”	$w$	$\cos wx$
“sin”	$w$	$\sin wx$
“alg”	$\alpha ,\beta$	$g(x)$
“alg-loga”	$\alpha ,\beta$	$g(x)\log (x-a)$
“alg-logb”	$\alpha ,\beta$	$g(x)\log (b-x)$
“alg-log”	$\alpha ,\beta$	$g(x)\log (x-a)\log (b-x)$
“cauchy”	$c$	$\frac{1}{x-c}$

其中，$g(x)=(x-a{)}^{\alpha }\ast (b-x{)}^{\beta }$

设func为$f(x)=x$，若weight参数为cos，而wvar取值为$w$，则实际计算的积分表达式为

$${\int }_a^b\cos wf(x)\text{d}x$$

示例如下

func = lambda x : x
quad(func, 0, np.pi)    # (4.935, 5.478e-14)
quad(func, 0, np.pi, weight='cos', wvar=1)  # (-2.00, 1.926e-13)

多重积分

在【si】中，除了quad之外，还提供了二重、三重以及N重积分的API，分别是【dblquad, tplquad, nquad】，三者所需参数如下

MIN = 1.49e-08
dblquad(func, a, b, gfun, hfun, args=(), epsabs=MIN, epsrel=MIN)
tplquad(func, a, b, gfun, hfun, qfun, rfun, args=(), epsabs=MIN, epsrel=MIN)
nquad(func, ranges, args=None, opts=None, full_output=False)

dblquad

以二重积分为例，其对应的问题可表述为下式

$${\int }_a^b{\int }_{y_g(x)}^{y_h(x)}f(y,x)\text{d}x\text{d}y$$

在函数dblquad中，func对应$f(y,x)$，a,b对那个上式的$a,b$，gfun, hfun对应上式的$y_g(x),y_h(x)$。

接下来求解下面的积分

$$\begin{array}{rlr}& {\int }_1^2{\int }_{x^2}^{x^3}xy\text{d}y\text{d}x& \\ =& {\int }_1^2\frac{1}{2}(x{y}^2){|}_{x^2}^{x^3}\text{d}x=& {\int }_1^2\frac{1}{2}(x^7-x^5)\text{d}x\\ =& \frac{1}{2}(\frac{1}{8}{x}^8-\frac{1}{6}{x}^6){|}_1^2=& \frac{1}{2}(\frac{2^8}{8}-\frac{2^6}{6})+\frac{1}{48}\\ =& \frac{513}{48}& \end{array}$$

Python代码如下

from scipy.integrate import dblquad
func = lambda x,y : x*y
gf = lambda x: x**2
hf = lambda x: x**3
dblquad(func, 1, 2, gf, hf)
# (10.6875, 5.284867210146833e-13)

计算结果与$\frac{513}{48}$一致。

与二重积分相比，三重积分tplquad只是多了一组qfun和rfun，相当于z处于qfun(x,y)和rfun(x,y)之间。

【nquad】貌似不支持回调函数，其参数ranges是元组的列表，每个元组代表对应未知量的取值范围。若将其映射为三重积分函数，则ranges可表示为$((a_1,b_1),(a_2,b_2),\cdots ,(a_n,b_n))$

下面仍以函数func为例，用nquad得出结果

from scipy.integrate import nquad
nquad(func, [[1,2], [3, 4]])
#(0.39276170758930756, 4.91851540406507e-15)