利用Python实现数值积分的方法

2024-04-02 19:04:59 518人浏览薄情痞子

Python 官方文档：入门教程 => 点击学习

摘要

目录1. 栗子2. 矩形计算面积2.1 左侧边长计算面积2.2 右侧边长计算面积2.3 中值边长计算面积3. 梯形计算面积4. 真值比对5. 总结1. 栗子为了加深大家的印象，首先

1. 栗子

为了加深大家的印象，首先我们来看个例子：

图示如下：

2. 矩形计算面积

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

2.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

2.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

2.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

3. 梯形计算面积

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

样例代码如下：

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

4. 真值比对

最后，我们来针对不同的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];

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