使用 Python 为数值积分建模辛普森规则
辛普森规则是一种非常强大的数值积分工具。为了在需要更少除法的情况下最大限度地提高准确性，辛普森规则使用加权因子计算积分。梯形规则仅考虑两个点，𝑥𝑖 和 𝑥𝑖+1，来估计梯形的面积，但辛普森规则在每一轮中使用两个以上的点(或多个条带)(参见下图)。通过权衡各种标准来更新各个点的 𝑓(𝑥) 值，以减少误差
辛普森 1/3 规则
此方法同时计算两个条带的面积，因此每次都会考虑三个不同的"x"值。为了覆盖整个域，条带数 n 必须完全为偶数。
对于前两个条带，辛普森 1/3 规则公式如下 −
$$\mathrm{A=\frac{1}{3}h(f(x_0)+4f(x_1)+f(x_2))}$$
因此，条带配对的总数将是−
$$\mathrm{I=\frac{1}{3}h((x_0)+4f(x_1)+f(x_2)))+\frac{1}{3}h(f(x_2)+4f(x_3)+f(x_4))+...}$$
$$\mathrm{+\frac{1}{3}h(f(x_{n-4})+4f(x_{n-3})+f(x_{n-2}))}$$
$$\mathrm{+\frac{1}{3} h(f(x_{n-2})+4f(x_{n-1})+f(x_n))}$$
为了便于编程，简化形式可以写成如下形式−
$$\mathrm{I=\frac{1}{3}h{(f(x_0)+f(x_n))+4(f(x_1)+f(x_3)+...+f(x_{n-3})+f(x_{n-1}))+2(f(x_2)+f(x_4)+...+f(x_{n-4})+f(x_{n-2}))}}$$
最后，将这些项以求和形式表示出来−
$$\mathrm{I=\frac{1}{3}h{(f(a)+f(b))+\sum^{n-1}_{i=1,3,5}4f(x_i)+\sum^{n-2}_{2,4,6}2f(x_i)}}$$
这里，重要的是要记住"n"是偶数。
假设我们要执行以下积分 −
$$\mathrm{\int^{\frac{\pi}{2}}_0x\:sin(x)dx}$$
现在编程过程可以写成 −
导入模块 −
# 导入数组模块 from numpy import *
定义必须执行积分的函数 −
# 定义函数 def f(x): return x*sin(x)
输入函数的左、右极限 (𝑎, 𝑏) 以及要划分面积的梯形数量，即 𝑛 (𝒏 应为偶数)。
# 定义函数域(即左、右极限) a=0 b=pi/2 # 定义梯形数量 n=4
计算每个梯形的宽度 (ℎ = (𝑏 − 𝑎)/𝑛)
# 计算梯形的宽度 h=(b-a)/n
创建一个 𝑥 数组。由于梯形的数量为 n，因此 𝑥 的数量将为 𝑛 + 1， 即"𝑥"数组将具有 𝑛 + 1 个元素。这是因为 NumPy 数组的索引从 0 开始，而不是从 1 开始。
# 创建 x 数组 #(注意:对于 n 个梯形(x 的数量/节点 = n+1)) x=linspace(a,b,n+1)
首先评估 $\mathrm{\frac{1}{2}(𝑓(𝑎) + 𝑓(𝑏))}$ 并将结果分配给变量"I"。
# I 的第一项 I=0.5*(f(a)+f(b))
对从 𝑥 = 1 到 𝑥 = 𝑛 的方程 1 的评估运行循环(因为有 𝑛 + 1 个节点)。这里要注意的重要一点是，当主循环中的索引 j 为偶数时，将评估 $\mathrm{2𝑓(𝑥_𝑖)}$，否则将评估 $\mathrm{4𝑓(𝑥_𝑖)}$ 项。
# 对剩余的 n-2 个梯形求和
for j in range(1,n):
if j%2==0:
I=I+2*f(x[j])
else:
I=I+4*f(x[j])
最后，在循环外将"I"乘以 h/3，并将其分配给"I"并打印结果
I=(h/3)*I
print(f'I = {round(I,6)}')
完整的 Python 程序将是 −
from numpy import *
def f(x):
return x*sin(x)
# 左和右极限
a=0
b=pi/2
# 梯形数量
n=4
# 梯形宽度
h=(b-a)/n
# x 数组
# 对于 n 个梯形(节点数=n+1)
x=linspace(a,b,n+1)
# 的第一项I
I=(f(a)+f(b))
# 对剩余的 n-2 个梯形求和
for j in range(1,n):
if j%2==0:
I=I+2*f(x[j])
else:
I=I+4*f(x[j])
I=(h/3)*I
print(f'I = {round(I,6)}')
输出
对于 n=4，程序返回 I = 0.999591，非常接近准确结果。
辛普森 3/8 规则
它与辛普森 1/3 规则的唯一区别在于，它考虑了三条带，因此每次都会包含四个点。前三个条带的公式如下 −
$$\mathrm{A=\frac{3}{8}h(f(x_0)+3f(x_1)+3f(x_2)+f(x_3))}$$
因此，所有条带配对的总和为。
$$\mathrm{I=\frac{3}{8}[(f(x_0)+3f(x_1)+3f(x_2)+f(x_3))+(f(x_3)+3f(x_4)+3f(x_5)+f(x_6))+(f(x_{n-3})+3f(x_{n-2}+3f(x_n-1)+f(x_n))]}$$
将术语以摘要形式放在最后，以便进行编程更简单。
$$\mathrm{I=\frac{3}{8}[(f(a)+f(b))+\sum^{n-2}_{i=1,4,7}3[f(x_i)+f(x_{i+1})]+\sum^{n-3}_{i=3,6,9}2f(x_i)]}$$
条带的数量 n 必须是三的倍数。我们将使用两个循环，一个用于第二项，一个用于第三项，因为从上面的公式可以明显看出，第二项以 1 开头，第三项以 3 开头。循环计数器将增加三。
在这里，我们将处理与前面 1/3 案例中相同的问题−
$$\mathrm{\int^\frac{\pi}{2}_0x\:sin(x)dx}$$
该算法与 1/3 情况的唯一区别在于方程 2 的求解方式。如果仔细观察，您会明白 RHS 上第二项的索引从 1 开始，而第三项的索引从 3 开始。因此，循环必须分为两部分，如下所示 −
# 对剩余的 n-2 个梯形求和 j=1 while j<n: I=I+3*(f(x[j])+f(x[j+1])) j=j+3 j=3 while j<n: I=I+2*f(x[j]) j=j+3 I=(3*h/8)*I
完整的 Python 程序如下 −
from numpy import *
def f(x):
return x*sin(x)
# 左和右极限
a=0
b=pi/2
# 梯形数量
n=6
# 梯形宽度
h=(b-a)/n
# x 数组
# 对于 n 个梯形(节点数=n+1)
x=linspace(a,b,n+1)
# I 的第一项
I=(f(a)+f(b))
# 对剩余的 n-2 个梯形求和
j=1
while j<n:
I=I+3*(f(x[j])+f(x[j+1]))
j=j+3
j=3
while j<n:
I=I+2*f(x[j])
j=j+3
I=(3*h/8)*I
print(f'I = {round(I,6)}')
输出
对于"n=6"，程序返回 −
I = 0.999819
结论
在本教程中，我们学习了如何在 Python 中对辛普森规则进行建模以执行数值积分。1/3 和 3/8 规则的算法都经过了阐述和编程。
