- 数值积分
- 导入模块
- 数值积分
- 积分法则
- 辛普森求积公式 Simpson's rule
- 高斯求积公式 Gaussian quadrature
- 使用SciPy进行数值积分
- 高斯积分
- 无穷积分
- 发散函数积分
- 列表积分
- 多重积分
- 维数灾难
- 符号积分和任意精度积分
- 多重积分
- 曲线积分
- 积分变换
相对微分而言，积分的难度要大得多。虽然有很多可以用解析方法来计算和的积分，大部分情况下，我们需要使用数值的方法。
连续函数以及有限积分域的积分在单一维度上可以有效计算，但是对于带奇点或无限积分域的可积函数，即使是一维数值积分也很困难。二维积分和多重积分可以通过重复一维积分来进行计算，但是计算量会随着维度上升急剧增长。高维积分需要使用蒙特卡罗采样算法等技术。
除了进行数值计算外，积分还有很多其他同途：
- 积分方程（含有未知函数的积分进行运算的方程）经常在科学和工程中使用。积分方程一般很难求解，通常可以将它们离散化，然后转换为线性方程组。
- 积分变换，可用于不同域之间的函数和方程变换，例如傅里叶变换。
本部分我们将主要使用SciPy的integrate模块进行数值积分。针对高精度浮点运算的需要，我们也会搭配使用SymPy和mpmath进行任意精度积分。
import matplotlib.pyplot as plt import matplotlib as mpl import numpy as np from scipy import integrate import sympy import mpmath %reload_ext version_information %version_information numpy, matplotlib, scipy, sympy, mpmath数值积分
我们将计算形式为 I ( f ) = ∫ a b f ( x ) d x I(f) = int_a^b f(x) dx I(f)=∫abf(x)dx的定积分，积分的上下限分别为a和b，区间可以是有限的、半无穷的和无穷的。
积分 I ( f ) I(f) I(f)可以理解为积分函数 f ( x ) f(x) f(x)的曲线和x轴之间的面积。
x = np.linspace(0, 5, 100) plt.plot(x, np.sin(x)) plt.fill_between(x[np.sin(x) > 0], 0, np.sin(x[np.sin(x) > 0]), alpha = 0.6) plt.fill_between(x[np.sin(x) < 0], 0, np.sin(x[np.sin(x) < 0]), alpha = 0.6)
一种计算上述形式积分
I
(
f
)
I(f)
I(f)的策略是，将积分写成积分函数值的离散和：
I ( f ) = ∑ i = 1 n w i f ( x i ) + r n I(f) = sum_{i=1}^{n} w_i f(x_i)+r_n I(f)=i=1∑nwif(xi)+rn
其中 w i w_i wi是函数 f ( x ) f(x) f(x)在点 x i x_i xi处的权重， r n r_n rn是近似误差。
I ( f ) I(f) I(f)这个求和公式被称为n点求积法则，其中n的选择、点的位置、权重因子都会对计算的准确性和复杂性产生影响。
积分法则积分法则可以通过在区间[a, b]上对函数 f ( x ) f(x) f(x)进行插值推导而来。如果 x i x_i xi在区间[a, b]上是均匀间隔的，那么可以使用多项式插值，得到的公式被称为牛顿-科特斯求积公式。
使用中间点的零阶多项式逼近
f
(
x
)
f(x)
f(x)，可以得到：
∫
a
b
f
(
x
)
d
x
≈
f
(
a
+
b
2
)
∫
a
b
d
x
=
(
b
−
a
)
a
+
b
2
int_a^b f(x) dx approx f(frac{a+b}{2}) int_a^b dx = (b-a)frac{a+b}{2}
∫abf(x)dx≈f(2a+b)∫abdx=(b−a)2a+b
这被称为中点公式。
使用一阶多项式（线性）逼近
f
(
x
)
f(x)
f(x)，可以得到：
∫
a
b
f
(
x
)
d
x
≈
b
−
a
2
(
f
(
a
)
+
f
(
b
)
)
int_a^b f(x) dx approx frac{b-a}{2}(f(a)+f(b))
∫abf(x)dx≈2b−a(f(a)+f(b))
这被称为梯形公式公式。
使用二阶插值多项式，将会得到辛普森公式：
∫
a
b
f
(
x
)
d
x
≈
b
−
a
6
(
f
(
a
)
+
4
f
(
a
+
b
2
)
+
f
(
b
)
)
int_a^b f(x) dx approx frac{b-a}{6}(f(a)+4f(frac{a+b}{2})+f(b))
∫abf(x)dx≈6b−a(f(a)+4f(2a+b)+f(b))
我们这里使用SymPy符号化推导该公式。
a, b, X = sympy.symbols("a, b, x")
f = sympy.Function("f")
x = a, (a+b)/2, b # simpson's rule
w = [sympy.symbols("w_%d" % i) for i in range(len(x))]
q_rule = sum([w[i] * f(x[i]) for i in range(len(x))])
q_rule
为了得到适合权重因子
w
i
w_i
wi的值，我们使用多项式基函数
{
ϕ
n
(
x
)
=
x
n
}
n
=
0
2
left{ phi_n(x)=x^n right}_{n=0}^2
{ϕn(x)=xn}n=02对
f
(
x
)
f(x)
f(x)进行插值。
phi = [sympy.Lambda(X, X**n) for n in range(len(x))] phi
将求和公式中的
f
(
x
)
f(x)
f(x)替换为基函数
ϕ
n
(
x
)
phi_n(x)
ϕn(x)，对权重因子进行解析求解：
∑ i = 1 2 w i ϕ n ( x i ) = ∫ a b ϕ n ( x ) d x sum_{i=1}^{2} w_i phi_n(x_i) = int_a^b phi_n(x) dx i=1∑2wiϕn(xi)=∫abϕn(x)dx
eqs = [q_rule.subs(f, phi[n]) - sympy.integrate(phi[n](X), (X, a, b)) for n in range(len(phi))] eqs
通过求解该线性方程组可以得到权重因子的解析表达式。
w_sol = sympy.solve(eqs, w) w_sol
得到辛普森求积公式的解析表达式。
q_rule.subs(w_sol).simplify()高斯求积公式 Gaussian quadrature
牛顿-科特斯求积公式的采样点在积分区间上是均匀分布的，对于非均匀分布采样，可以使用高斯求积公式。
我们可以把函数 f ( x ) f(x) f(x)写作$ f(x) = W(x)g(x) ， 其 中 ，其中 ，其中g(x)$是近似多项式， W ( x ) W(x) W(x)是已知的权重函数，这样我们就有
∫ − 1 1 f ( x ) d x = ∫ − 1 1 W ( x ) g ( x ) d x ≈ ∑ i = 1 n w i ′ g ( x i ) int _{-1}^{1}f(x)dx=int _{-1}^{1}W(x)g(x)dxapprox sum _{i=1}^{n}w_{i}'g(x_{i}) ∫−11f(x)dx=∫−11W(x)g(x)dx≈i=1∑nwi′g(xi)
常见的权重函数有 W ( x ) = ( 1 − x 2 ) − 1 / 2 W(x)=(1-x^{2})^{-1/2} W(x)=(1−x2)−1/2（高斯切比雪夫）、 W ( x ) = e − x 2 W(x)=e^{{-x^{2}}} W(x)=e−x2（高斯埃米特）等。
权重函数 W ( x ) = 1 W(x)=1 W(x)=1时，关联多项式为勒让得多项式 P n ( x ) P_{n}(x) Pn(x)，这种方法通常称为高斯勒让德求积。
使用SciPy进行数值积分SciPy的integrate模块中的数值求积函数可以分为两类：一类将被积函数作为Python函数传入，另一类将被积函数在给定点的样本值以数组的形式传入。第一类函数使用高斯求积法（quad、quadrature、fixed_quad），第二类函数使用牛顿-科斯特求积法（trapezoid、simpson、romb）。
高斯积分quadrature函数是一个使用Python实现的自适应高斯求积程序。quadrature函数会重复调用fixed_quad函数，并不断增加多项式的次数，直到满足所需的精度。quad函数是对Fortran库QUADPACK的封装，有更好的性能。一般情况下优先使用quad函数。
考虑定积分
∫
−
1
1
e
−
x
2
d
x
int _{-1}^{1}e^{-x^2}dx
∫−11e−x2dx
def f(x):
return np.exp(-x**2)
val, err = integrate.quad(f, -1, 1)
val, err
quad函数返回一个元组，包含积分的数值结果和绝对误差估计。可以使用参数epsabs和epsrel来设置绝对误差和相对误差的容忍度。
quad函数的关键词参数args可以将函数的参数值传递给被积函数。
考虑含参数定积分
∫
−
1
1
a
e
−
(
x
−
b
)
2
/
c
2
d
x
int _{-1}^{1}ae^{-(x-b)^2/c^2}dx
∫−11ae−(x−b)2/c2dx
def f(x, a, b, c):
return a * np.exp(-((x-b)/c)**2)
val, err = integrate.quad(f, -1, 1, args=(1, 2, 3))
val, err
当被积分函数的变量不是第一个参数是，可以使用lambda函数来调整参数的顺序。
例如我们使用scipy.special模块的jv函数计算0阶贝塞尔函数的积分，jv函数的第一个参数是贝塞尔函数的阶数，第二个参数是变量x。
from scipy.special import jv val, err = integrate.quad(lambda x: jv(0, x), 0, 5) val, err无穷积分
quad函数支持无穷积分。浮点数中的无穷表达式可以使用np.inf获得。
考虑使用quad函数计算 ∫ − ∞ ∞ e − x 2 d x int _{-infty}^{infty}e^{-x^2}dx ∫−∞∞e−x2dx
f = lambda x: np.exp(-x**2) val, err = integrate.quad(f, -np.inf, np.inf) val, err发散函数积分
通过一些额外信息，quad函数也可以处理带可积奇点的函数。
考虑积分 ∫ − 1 1 1 ∣ x ∣ d x int _{-1}^{1} frac{1}{sqrt{|x|}}dx ∫−11∣x∣ 1dx，被积函数在 x = 0 x=0 x=0处是发散的。
f = lambda x: 1/np.sqrt(abs(x))
fig, ax = plt.subplots(figsize=(8, 3))
a, b = -1, 1
x = np.linspace(a, b, 10000)
ax.plot(x, f(x), lw=2)
ax.fill_between(x, f(x), color='green', alpha=0.5)
ax.set_xlabel("$x$", fontsize=18)
ax.set_ylabel("$f(x)$", fontsize=18)
ax.set_ylim(0, 25)
integrate.quad(f, a, b) # 直接使用quad函数求积分，可能会失败
integrate.quad(f, a, b, points=[0])
quad函数的points参数可以设置需要绕过的点，从而正确地计算积分。
quad函数只适合于被积函数可以被Python函数表示的积分。如果被积函数来自实验或者观察数据，其可能只可以在某些预先确定的点求值。这种情况下，可以使用牛顿-科特斯求积法，如之前介绍的中间点公式、梯形公式或辛普森公式。
在SciPy的integrate模块中，trapezoid和simpson函数分别实现了复化梯形公式和辛普森公式。这些函数的第一个参数是数组y，第二个参数是采样点数组x或者采样点间隔dx。
考虑积分
∫
0
2
x
d
x
int _{0}^{2} sqrt{x} dx
∫02x
dx，积分区间为[0, 2]，采样点数目25。
f = lambda x: np.sqrt(x) a, b = 0, 2 x = np.linspace(a, b, 25) y = f(x) fig, ax = plt.subplots(figsize=(8, 3)) ax.plot(x, y, 'bo') xx = np.linspace(a, b, 500) ax.plot(xx, f(xx), 'b-') ax.fill_between(xx, f(xx), color='green', alpha=0.5) ax.set_xlabel(r"$x$", fontsize=18) ax.set_ylabel(r"$f(x)$", fontsize=18)
通过解析积分可以得到积分的精确值。
val_exact = 2.0/3.0 * (b-a)**(3.0/2.0) val_exact
要计算这个积分，将数组y和x传递给trapezoid函数或simpson函数。
val_trapz = integrate.trapz(y, x) val_trapz - val_exact
val_simps = integrate.simps(y, x) val_simps - val_exact
trapezoid函数和simpson函数无法提供对误差的估计。提供准确度的方法是增加样品点的数量或者使用更高阶的方法。
integrate模块中的romb函数实现了Romberg方法。这种方法使用均匀间隔的采样点（数目为
2
n
+
1
2^n+1
2n+1），同时使用Richardson外推算法来加速梯形法的收敛。
x = np.linspace(a, b, 1 + 2**6) y = f(x) val_romb = integrate.romb(y, dx=(x[1]-x[0])) val_romb - val_exact
一般情况下，推荐使用simpson函数。
多重积分，例如二重积分
∫
a
b
∫
c
d
f
(
x
,
y
)
d
x
d
y
int _{a}^{b}int _{c}^{d}f(x, y)dxdy
∫ab∫cdf(x,y)dxdy和三重积分
∫
a
b
∫
c
d
∫
e
f
f
(
x
,
y
,
z
)
d
x
d
y
d
z
int _{a}^{b}int _{c}^{d}int _{e}^{f}f(x, y, z)dxdydz
∫ab∫cd∫eff(x,y,z)dxdydz，可以使用SciPy的integrate模块的dblquad和tplquad函数。n个变量的积分也可以使用nquad函数来计算。这些函数都是对单变量求积函数quad的封装，沿着被积函数每个维度重复调用quad函数。
考虑二重积分
∫
0
1
∫
0
1
e
−
x
2
−
y
2
d
x
d
y
int _{0}^{1}int _{0}^{1}e^{-x^2-y^2}dxdy
∫01∫01e−x2−y2dxdy
def f(x, y):
return np.exp(-x**2-y**2)
fig, ax = plt.subplots(figsize=(6, 5))
x = y = np.linspace(-1.25, 1.25, 75)
X, Y = np.meshgrid(x, y)
c = ax.contour(X, Y, f(X, Y), 15, cmap=mpl.cm.RdBu, vmin=-1, vmax=1)
bound_rect = plt.Rectangle((0, 0), 1, 1,
facecolor="grey")
ax.add_patch(bound_rect)
ax.axis('tight')
ax.set_xlabel('$x$', fontsize=18)
ax.set_ylabel('$y$', fontsize=18)
本例中，x和y的积分限都是固定的。dblquad函数要求y变量的积分限用变量为x的函数表示，因此我们需要定义两个函数
g
(
x
)
g(x)
g(x)和
h
(
x
)
h(x)
h(x)，二者返回常量。
a, b = 0, 1 g = lambda x: 0 h = lambda x: 1 val, err = integrate.dblquad(f, a, b, g, h) # 参数g和h为函数 val, err
对于形如
∫
0
1
∫
0
1
∫
0
1
e
−
x
2
−
y
2
−
z
2
d
x
d
y
d
z
int _{0}^{1}int _{0}^{1}int _{0}^{1}e^{-x^2-y^2-z^2}dxdydz
∫01∫01∫01e−x2−y2−z2dxdydz的三重积分，可以使用nquad函数计算。除了继续使用
g
(
x
)
g(x)
g(x)和
h
(
x
)
h(x)
h(x)表示y轴的积分限之外，我们需要额外提供两个函数
q
(
x
,
y
)
q(x, y)
q(x,y)和
r
(
x
,
y
)
r(x, y)
r(x,y)表示z轴的积分限。注意这里使用了x和y两个坐标。
def f(x, y, z):
return np.exp(-x**2-y**2-z**2)
val, err = integrate.tplquad(f, 0, 1, lambda x : 0, lambda x : 1, lambda x, y : 0, lambda x, y : 1)
val, err
对于任意维度数目的积分，可以使用nquad函数。它的第二个参数是用于指定积分限的列表，列表里面包含每个积分变量的积分限元组，或者包含能够返回积分限的函数。
integrate.nquad(f, [(0, 1), (0, 1), (0, 1)]) # 元素相同的列表可以使用 [(0, 1)] * 3 生成维数灾难
在数值积分中，多重积分的计算复杂度和维数成指数关系。
def f(*args):
return np.exp(-np.sum(np.array(args)**2))
for dim in range(1, 5):
print(dim)
%time integrate.nquad(f, [(0,1)] * dim)
随着维数增加，高维积分的直接求积方法变得不切实际。如果对计算精度要求不是非常高，可以使用蒙特卡罗积分。蒙特卡罗积分在维度上的扩展性非常好，对于高维积分是一种有效的工具。
在符号计算的章节，我们已经演示了使用SymPy的sympy.integrate函数来计算符号函数的有限积分和无限积分。
例如，计算积分 ∫ − 1 1 2 1 − x 2 d x int_{-1}^{1} 2 sqrt{1-x^2} dx ∫−1121−x2 dx
x = sympy.symbols("x")
f = 2 * sympy.sqrt(1-x**2)
a, b = -1, 1
sympy.plot(f, (x, -2, 2))
val_sym = sympy.integrate(f, (x, a, b)) val_sym
存在精确解析解的问题是一种特例。数值方法的精度受算法和浮点精度制约。mpmath库为任意精度计算提供了数值方法。为了按照给定的精度计算积分，可以使用mpmath.quad函数。为了设置精度，我们把变量mpmath.mp.dps设置为所需精度的小数位数。
mpmath.mp.dps = 75 # 75位小数的精度
被积分的Python函数可以使用mpmath库中的数学函数。可以使用sympy.lambdify从一个SymPy表达式中创建这样的函数。
f_mpmath = sympy.lambdify(x, f, 'mpmath') val = mpmath.quad(f_mpmath, (a, b)) val # 返回的对象类型是高精度浮点数
sympy.sympify(val)
我们可以将这个结果与解析结果进行对比。
sympy.N(val_sym, mpmath.mp.dps+1) - val
SciPy的integrate模块中的quad函数无法达到这样的精度，因为受浮点数精度的限制。
mpmath库中的quad函数也可以用于计算多重积分。计算这样的积分，只需要把拥有多个变量参数的被积函数传给quad函数，并未每个积分变量传入积分限的元组。
计算下面的双重积分 ∫ 0 1 ∫ 0 1 cos ( x ) cos ( y ) e − x 2 − y 2 d x int_{0}^{1} int_{0}^{1} cos(x) cos(y) e^{-x^2-y^2} dx ∫01∫01cos(x)cos(y)e−x2−y2dx
def f2(x, y):
return np.cos(x)*np.cos(y)*np.exp(-x**2-y**2)
integrate.dblquad(f2, 0, 1, lambda x : 0, lambda x : 1)
x, y, z = sympy.symbols("x, y, z")
f2 = sympy.cos(x)*sympy.cos(y)*sympy.exp(-x**2-y**2)
f2_mpmath = sympy.lambdify((x, y), f2, 'mpmath')
mpmath.mp.dps = 30 # 指定计算精度
res = mpmath.quad(f2_mpmath, (0, 1), (0, 1))
res
由于高精度浮点运算是在CPU上模拟实现的，缺点是非常慢。
SymPy可以使用line_integral函数来计算形如 ∫ C f ( x , y ) d s int_{C} f(x, y) ds ∫Cf(x,y)ds的曲线积分，其中C是x-y平面上的曲线。该函数的第一个参数是SymPy表示的被积函数，第二个参数是一个sympy.Curve实例，第三个参数是积分变量的列表。
例如，创建Curve实例来表示单位圆的路径。
t, x, y = sympy.symbols("t, x, y")
C = sympy.Curve([sympy.cos(t), sympy.sin(t)], (t, 0, 2 * sympy.pi))
积分路径指定之后，我们对被积函数 f ( x , y ) = 1 f(x, y)=1 f(x,y)=1进行积分，积分结果就是圆的周长。
sympy.line_integrate(1, C, [x, y])积分变换
积分变换是将一个函数作为输入，然后输出另有一个函数的过程。这里我们介绍使用SymPy支持的两种积分变换，拉普拉斯变换和傅里叶变换。这两种变换有很多应用，例如可以使用拉普拉斯变换将微分方程转换为代数方程，或者使用傅里叶变换将时域问题转换为频域问题。
一般来说，函数 f ( t ) f(t) f(t)的积分变换可以写为：
T f ( u ) = ∫ t 1 t 2 K ( t , u ) f ( t ) d t T_f(u) = int_{t_1}^{t_2} K(t, u) f(t) dt Tf(u)=∫t1t2K(t,u)f(t)dt
其中 T f ( u ) T_f(u) Tf(u)是变换后的函数， K ( t , u ) K(t, u) K(t,u)是变换的核函数。
积分变换的逆变换是：
f ( t ) = ∫ u 1 u 2 K − 1 ( t , u ) T f ( u ) d u f(t) = int_{u_1}^{u_2} K^{-1}(t, u) T_f(u) du f(t)=∫u1u2K−1(t,u)Tf(u)du
其中 K − 1 ( t , u ) K^{-1}(t, u) K−1(t,u)是核函数的逆变换。
拉普拉斯变换
F
(
s
)
=
∫
0
∞
f
(
t
)
e
−
s
t
d
t
F(s)=int _{0}^{infty }f(t)e^{-st},mathrm {d} t
F(s)=∫0∞f(t)e−stdt
及其逆变换
f
(
t
)
=
L
−
1
{
F
}
(
t
)
=
1
2
π
i
lim
T
→
∞
∫
γ
−
i
T
γ
+
i
T
e
s
t
F
(
s
)
d
s
f(t)={mathcal {L}}^{-1}{F}(t)={frac {1}{2pi i}}lim _{Tto infty }int _{gamma -iT}^{gamma +iT}e^{st}F(s),mathrm {d} s
f(t)=L−1{F}(t)=2πi1T→∞lim∫γ−iTγ+iTestF(s)ds
以及傅里叶变换
f
^
(
ξ
)
=
∫
−
∞
∞
f
(
x
)
e
−
2
π
i
x
ξ
d
x
hat{f}(xi) = int_{-infty}^infty f(x) e^{- 2pi i x xi},dx
f^(ξ)=∫−∞∞f(x) e−2πixξdx
及其逆变换
f
(
x
)
=
∫
−
∞
∞
f
^
(
ξ
)
e
2
π
i
ξ
x
d
ξ
f(x) = int_{-infty}^infty hat f(xi) e^{2 pi i xi x},dxi
f(x)=∫−∞∞f^(ξ) e2πiξxdξ
在SymPy中，拉普拉斯变换和傅里叶变换对应的函数为sympy.laplace_transform和sympy.fourier_transform，其逆变换为sympy.inverse_laplace_transform和sympy.inverse_fourier_transform。
示例 计算函数
f
(
t
)
=
sin
(
a
t
)
f(t) = sin(at)
f(t)=sin(at)的拉普拉斯变换
s = sympy.symbols("s")
a, t = sympy.symbols("a, t", positive=True)
f = sympy.sin(a*t)
sympy.laplace_transform(f, t, s)
默认情况下，laplace_transform函数返回一个元组，包含转换结果和变换的收敛条件。如果只需要转换结果，可以使用参数noconds=True去除返回中的条件。
F = sympy.laplace_transform(f, t, s, noconds=True) F
逆变换将得到原始函数：
sympy.inverse_laplace_transform(F, s, t, noconds=True)
我们将在微分方程章节中拉普拉斯变换的应用。
示例 计算函数
f
(
t
)
=
e
−
a
t
2
f(t) = e^{-at^2}
f(t)=e−at2的傅里叶变换
w = sympy.symbols("omega")
a, t = sympy.symbols("a, t", positive=True)
f = sympy.exp(-a*t**2)
F = sympy.fourier_transform(f, t, w)
F
sympy.inverse_fourier_transform(F, w, t)
我们在将在信号处理章节中演示傅里叶变换在降噪中的应用。
参考文献来自桑鸿乾老师的课件：科学计算和人工智能
