Python中梯形法则数值积分建模

---

积分（定积分）的目的是计算函数曲线在两个极限值a和b之间区域的面积。下图将进一步阐明这个概念。

求积法，也常称为数值积分，是一种用于评估给定函数曲线下面积的方法。过程非常简单，即首先我们将限定区域分成几个区域或条带。然后，借助简单矩形的数学公式评估面积。然后将所有条带的面积相加，以获得函数曲线在a点和b点之间区域的总面积。

在以下情况下，数值积分方法变得非常重要：

• 如果函数的积分不存在直接/显式公式，或者

• 必须对从经验数据获得的某些曲线进行积分。

数值积分方法的精度取决于恰当覆盖曲线下方区域的条带/分割数量。如果能够用最少的条带精确地评估面积，则数值积分方法被称为最佳方法。

梯形法则

梯形法则是由工程师们几十年来用于数值执行积分任务的最基本技术。在这种方法中，给定曲线下的区域被分成具有相同宽度h的垂直梯形，如下图所示。还可以观察到梯形的顶部刚好接触曲线。

第一个梯形的面积可以计算如下：

$$\mathrm{A_{1} \: = \: \frac{h}{2} \lgroup f(x_{0}) \: + \: f(x_{1}) \rgroup}$$

以类似的方式，可以评估其他矩形的连续面积。现在，积分的基本原理指出，必须将每个梯形的面积从$\mathrm{x_{0} \: = \: a}$加到$\mathrm{x_{n} \: = \: b}$以获得积分I。

数学上，可以写成：

$$\mathrm{I \: = \: \frac{h}{2} \lgroup f(x_{0}) \: + \: f(x_{1}) \rgroup \: + \: \frac{h}{2} \lgroup f(x_{1}) \: + \: f(x_{2}) \rgroup \: + \: \: \dotso \: + \frac{h}{2} \lgroup f(x_{n-2}) \: + \: f(x_{n-1}) \rgroup \: + \frac{h}{2} \lgroup f(x_{n-1}) \: + \: f(x_{n}) \rgroup}$$

现在需要注意的是，从$\mathrm{f(x_{1})}$到$\mathrm{f(x_{n-1})}$，函数出现了两次，因此它们都将具有系数h，除了第一项$\mathrm{(f(x_{0}) \: = \: f(a))}$和最后一项$\mathrm{(f(x_{1}) \: = \: f(b))}$。因此，它们可以写成：

$$\mathrm{I \: = \: h(\frac{1}{2} \lgroup f(a) \: + \: f(b) \rgroup \: + \: f(x_{1}) \: + \: f(x_{2}) \: + \: \: \dotso \: + \: f(x_{n-2}) \: + \: f(x_{n-1}) \rgroup}$$

在求和符号表示法中（这很容易在任何编程语言中实现和建模），上述方程也可以写成：

$$\mathrm{I \: = \: h(\frac{1}{2} \lgroup f(a) \: + \: f(b) \rgroup \: + \: \sum_{i=1}^{i=n-1} f(x_{i}) )}$$

Python中梯形法则的实现

为了理解Python中梯形法则的实现，让我们考虑以下问题：

$$\mathrm{\int_{0}^{\frac{\pi}{2}} x \cos(x) \:dx}$$

现在，步骤可以写成：

• 导入模块

# Importing module for array 
from numpy import *

• 定义需要进行积分的函数

# Defining function
def f(x):
   return x*cos(x)

• 输入函数的左右极限 (a, b) 以及要将区域划分的梯形数 n。

# Defining function domain (i.e. Left and Right limits)
a=0
b=pi/2

# Defining number of trapezoids
n=5

• 计算每个梯形的宽度 (h = (b − a)/n)

# Evaluating the width of trapezoids
h=(b-a)/n

• 创建x数组。由于梯形数为n，因此x的数量将为n+1，即x数组将包含n+1个元素。这是因为NumPy数组的索引从0开始，而不是从1开始。

# Creating array of x
# (Note: For n trapezoids (no. of x's/nodes = n+1))
x=linspace(a,b,n+1)

• 首先计算$\mathrm{\frac{1}{2} \lgroup f(a) \: + \: f(b) \rgroup}$并将结果赋值给变量I。

# First term of I
I=0.5*(f(a)+f(b))

• 运行循环以计算x = 1到x = n的$\mathrm{\sum f(x_{i})}$（因为有n + 1个节点）。然后通过$\mathrm{I \: + \: \sum f(x_{i})}$更新I的值。最后，在循环外部，将I乘以h并将其赋值给I。

# Summing over remaining n-2 trapezoids
for j in range(1,n): 
   I=I+f(x[j]) 

I=h*I

示例

完整的程序如下：

# Importing module for array 
from numpy import *

# Defining function
def f(x):
   return x*cos(x)

# Defining function domain (i.e. Left and Right limits)
a=0
b=pi/2

# Defining number of trapezoids
n=5

# Evaluating the width of trapezoids
h=(b-a)/n

# Creating array of x
# (Note: For n trapezoids (no. of x's/nodes = n+1))
x=linspace(a,b,n+1)

# First term of I
I=0.5*(f(a)+f(b))

# Summing over remaining n-2 trapezoids
for j in range(1,n): 
   I=I+f(x[j]) 

I=h*I
print(f'I = {round(I,6)}')

输出

上述程序的输出将是：

I = 0.54959

如果要进一步提高结果的精度，则必须增加梯形的数量。

如果我们针对n=10运行上述代码，则程序输出将是：

I = 0.570585

此结果精确到小数点后5位。

如果要应用上述程序，则必须在上述代码中更改函数、其极限 (a,b) 和梯形数 (n)。

结论

在本教程中，我们使用Python演示了用于评估数值积分的梯形法。我们详细解释了方法和Python代码。