3.5.1 球的体积

数值积分是对定积分的数值求解，例如可以利用数值积分计算某个形状的面积。让我们先考虑一下如何计算半径为1的半圆的面积。根据圆的面积公式，其面积应该等于π/2。单位半圆的曲线方程为，可以通过下面的half_circle( )进行计算：

    def half_circle(x):
        return (1-x**
2)**
0.5

最简单的数值积分算法就是将要积分的面积分为许多小矩形，然后计算这些矩形的面积之和。下面使用这种方法，将X轴上-1到1的区间分为10000等份，然后计算面积和：

    N = 10000
    x = np.linspace(-1, 1, N)
    dx = x[1] - x[0]
    y = half_circle(x)
    2 *
 dx *
 np.sum(y) # 面积的两倍
    3.1415893269307373

也可以用NumPy的trapz( )计算半圆上的各点所构成的多边形的面积，trapz( )计算的是以(x,y)为顶点坐标的折线与X轴所夹的面积：

    np.trapz(y, x) *
 2 # 面积的两倍
    3.1415893269315975

如果使用integrate.quad( )进行数值积分，就能得到非常精确的结果：

    from scipy import integrate
    pi_half, err = integrate.quad(half_circle, -1, 1)
    pi_half *
 2
    3.141592653589797

计算多重定积分可以通过多次调用quad( )实现，为了调用方便，integrate模块提供了dblquad( )进行二重定积分，提供了tplquad( )进行三重定积分。下面以计算单位半球体积为例，说明dblquad( )的用法。

单位半球面上的点(x,y,z)满足方程x2 +y2 +z2 =1，因此下面的half_sphere( )可以通过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( )计算半球体积：

    volume, error = integrate.dblquad(half_sphere, -1, 1, 
            lambda x:-half_circle(x), 
            lambda x:half_circle(x))
    
    print volume, error, np.pi*
4/3/2
    2.09439510239 2.32524566534e-14 2.09439510239

dblquad( )的调用参数为：dblquad(func2d, a, b, gfun, hfun)。其中，func2d是需要进行二重积分的函数，它有两个参数，假设分别为x和y。a和b参数指定被积分函数的第一个变量x的积分区间，而gfun和hfun参数指定第二个变量y的积分区间。gfun和hfun是函数，它们通过变量x计算出变量y的积分区间，这样可以对X-Y平面上的任何区间对func2d进行积分。

图3-24是半球体积的积分的示意图。从此示意图可以看出，X轴的积分区间为-1.0到1.0，对于X轴上的某点x0，通过对Y轴的积分可以计算出图中深色的垂直切面的面积，因此Y轴的积分区间如图中的点线所示。

图3-24 半球体二重积分示意图