用Python实现数值积分

1. 引言

本文主要用于对比使用Python来实现数学中积分的几种计算方式，并和真值进行对比，加深大家对积分运算实现方式的理解。

闲话少说，我们直接开始吧。 ：）

2. 举个栗子

为了加深大家的印象，首先我们来看个例子：
对下述三次函数
$$f(x)=x^3-4x^2+4x+2$$
进行积分运算：
$$I={\int }_a^{b}f(x)dx={\int }_a^b(x^3-4x^2+4x+2)dx$$
图示如下：

如果 $a=1/2$ , $b=5/2$,则上述积分的值为：

$$I=(\frac{x^4}{4}-\frac{4x^3}{3}+2x^2+2x){|}_a^b=\frac{61}{12}\approx 5.0833$$

3. 矩形计算面积

我们知道，在数学中，积分运算表示上述曲线和x轴围成的封闭区域的面积，为此，我们在数值预算中，来近似计算上述区域的面积，最直观的想法就是拆成一个个小的矩形来计算对应的面积。

3.1 左侧边长计算面积

为了计算每个小矩形的面积，设计到边长高的选择，这里我么以左侧函数取值作为对应矩形的高来计算相应的小矩形的面积，图示如下：

对应的代码如下：

import numpy as np
x = np.linspace(0, 3, 1001)
f = lambda x: x**3 - 4*x**2 + 4*x + 2
a = 0.5
b = 2.5
Ax = np.linspace(a, b, 101)
Ay = f(Ax)
def defInt_left(f, a, b, N):
    # left-hand point
    result = 0; FX = []; Xn = []
    dx = abs(b - a)/N
    while a < b:
        result += f(a)*dx
        FX += [f(a)]
        Xn += [a]
        a += dx
    return result, FX, Xn, dx
N = 4
I_left, FX, Xn, dx = defInt_left(f, a, b, N)
print(I_left)

上述代码中，我们将横坐标拆分为4小份，也就是拆分成4个小矩形，然后使用函数左侧的点坐位小矩形的高，上述代码的运行结果如下：

5.25

3.2 右侧边长计算面积

这里和上述原理类似，只不过每个小矩形的高采用右侧边长函数取值来近似计算，图例如下：

样例代码如下：

def defInt_right(f, a, b, N):
    # right-hand point
    result = 0; FX = []; Xn = []
    dx = abs(b - a)/N
    while a < b:
        result += f(a + dx)*dx
        FX += [f(a + dx)]
        Xn += [a]
        a += dx
    return result, FX, Xn, dx

N = 4
I_right, FX, Xn, dx = defInt_right(f, a, b, N)
print(I_right)

运行结果如下：

5.0

3.3 中值边长计算面积

看了上述两种近似计算方式，有同学就说有取左侧点算出来面积大的，有取右侧点算出来面积小的，那干脆折中一下，我们来以中值坐位矩形的高来计算对应的面积。图例如下：

代码实现如下：

def defInt_middle(f, a, b, N):
    # middle point
    result = 0; FX = []; Xn = []
    dx = abs(b - a)/N
    while a < b:
        result += f(a + dx/2)*dx
        FX += [f(a + dx/2)]
        Xn += [a]
        a += dx
    return result, FX, Xn, dx

N = 4
I_mid, FX, Xn, dx = defInt_middle(f, a, b, N)
print(I_mid)

运行结果如下：

5.0625

4. 梯形计算面积

读到这里的同学可能会思考，既然可以将封闭区域划分成一个个的小矩形，那当然也可以将其划分成梯形来近似计算相应的面积，图例如下：

样例代码如下：

def defInt_trapezoid(f, a, b, N):
    # trapezoidal rule
    result = 0; FXa, FXb = [], []; Xn = []
    dx = abs(b - a)/N
    while a < b:
        result += (f(a) + f(a + dx))*dx/2
        FXa += [f(a)]; FXb += [f(a + dx)]
        Xn += [a]
        a += dx
    return result, FXa, FXb, Xn, dx

N = 4
I_trap, FXa, FXb, Xn, dx = defInt_trapezoid(f, a, b, N)
print(I_trap)

运行结果如下：

5.125

5. 真值比对

最后，我们来针对不同的N来讲封闭区域划分成对应的小份，分别针对性的计算上述四种方式的积分值，样例代码如下：

Nx = range(1, 11)
I1, I2, I3, I4 = [], [], [], []
for Ni in Nx:
    i1, *_ = defInt_left(f, a, b, Ni); I1 += [i1];
    i2, *_ = defInt_right(f, a, b, Ni); I2 += [i2];
    i3, *_ = defInt_middle(f, a, b, Ni); I3 += [i3];
    i4, *_ = defInt_trapezoid(f, a, b, Ni); I4 += [i4];

最后将其与真值进行对比，如下：

可以看出，随着划分区域的增多，采用梯形计算面积方式最逼近真值。

6. 总结

本文重点介绍了使用不同面积划分方法来近似计算积分取值的原理和相应的代码实现,其中采用梯形计算面积的方式随着划分子区域数目的增加最接近真值，推荐大家使用。

您学废了吗？

关注公众号《AI算法之道》，获取更多AI算法资讯。