数值积分是对定积分的数值求解,例如可以利用数值积分计算某个形状的面积。下面让我们来考虑一 下如何计算半径为1的半圆的面积,根据圆的面积公式,其面积应该等于PI/2。单位半圆曲线可以用下 面的函数表示:
import scipy.optimize as opt
import numpy as np
def half_circle(x):
return (1-x**2)**0.5
下面的程序使用经典的分小矩形计算面积总和的方式,计算出单位半圆的面积:
N = 10000
x = np.linspace(-1, 1, N)
dx = 2.0/N
y = half_circle(x)
dx * np.sum(y[:-1] + y[1:]) # 面积的两倍
np.float64(3.1412751679989044)
利用上述方式计算出的圆上一系列点的坐标,还可以用numpy.trapz进行数值积分:
import numpy as np
np.trapz(y, x) * 2 # 面积的两倍
/tmp/ipykernel_1672/3394964639.py:2: DeprecationWarning: `trapz` is deprecated. Use `trapezoid` instead, or one of the numerical integration functions in `scipy.integrate`. np.trapz(y, x) * 2 # 面积的两倍
np.float64(3.1415893269315975)
此函数计算的是以x,y为顶点坐标的折线与X轴所夹的面积。同样的分割点数,trapz函数的结果更加接 近精确值一些。 如果我们调用scipy.integrate库中的quad函数的话,将会得到非常精确的结果:
from scipy import integrate
pi_half, err = integrate.quad(half_circle, -1, 1)
pi_half*2
3.1415926535897967
多重定积分的求值可以通过多次调用quad函数实现,为了调用方便,integrate库提供了dblquad函数 进行二重定积分,tplquad函数进行三重定积分。下面以计算单位半球体积为例说明dblquad函数的用 法。 单位半球上的点(x,y,z)符合如下方程: x 2 + y 2 + z 2 = 1 因此可以如下定义通过(x,y)坐标计算球面上点的z值的函数:
def half_sphere(x, y):
return (1-x**2-y**2)**0.5
X-Y轴平面与此球体的交线为一个单位圆,因此积分区间为此单位圆,可以考虑为X轴坐标从-1到1进 行积分,而Y轴从 -half_circle(x) 到 half_circle(x) 进行积分,于是可以调用dblquad函数:
integrate.dblquad(half_sphere, -1, 1, lambda x:-half_circle(x),lambda x:half_circle(x))
(2.0943951023931984, 1.0002354500215915e-09)
通过球体体积公式计算的半球体积：
np.pi*4/3/2
2.0943951023931953
dblquad函数的调用方式为:
dblquad(func2d, a, b, gfun, hfun)对于func2d(x,y)函数进行二重积分,其中a,b为变量x的积分区间,而gfun(x)到hfun(x)为变量y的积分 区间。 半球体积的积分的示意图如下: 图 3.3 - 半球体积的双重定积分示意图 X轴的积分区间为-1.0到1.0,对于X=x0时,通过对Y轴的积分计算出切面的面积,因此Y轴的积分区间 如图中红色点线所示。