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}_{i=2,4,6}2f(x_i)}}$$

这里，重要的是要记住“n”是偶数。

假设我们想要执行以下积分：

$$\mathrm{\int^{\frac{\pi}{2}}_0x\:sin(x)dx}$$

现在可以编写编程过程：

导入模块

# Importing module for array
from numpy import * 

定义需要进行积分的函数

# Defining the function
def f(x):
   return x*sin(x)

输入函数的左右极限(𝑎, 𝑏)以及要将面积划分的梯形数𝑛 (**𝑛必须是偶数**)。

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

计算每个梯形的宽度(ℎ = (𝑏 − 𝑎)/𝑛)

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

创建𝑥数组。由于梯形数为n，所以𝑥的个数为𝑛 + 1，即“𝑥”数组将有𝑛 + 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}(𝑓(𝑎) + 𝑓(𝑏))}$并将结果赋值给变量“I”。

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

运行一个循环，用于计算𝑥 = 1到𝑥 = 𝑛的公式1（因为有𝑛 + 1个节点）。这里需要注意的是，当主循环中的索引j为偶数时，将计算$\mathrm{2𝑓(𝑥_𝑖)}$，否则将计算$\mathrm{4𝑓(𝑥_𝑖)}$项。

# Summing over remaining n-2 trapezoids
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)
# Left and Right limits
a=0
b=pi/2
# Number of trapezoids
n=4
# Width of trapezoid
h=(b-a)/n
# array of x
# For n trapezoids (no. of nodes=n+1)
x=linspace(a,b,n+1)
# First term of I
I=(f(a)+f(b))
# Summing over remaining n-2 trapezoids
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必须是3的倍数。我们将使用两个循环，一个用于第二项，一个用于第三项，因为从上面的公式可以明显看出，第二项从1开始，第三项从3开始。循环计数器将以三递增。

在这里，我们将采用与之前的1/3情况相同的示例：

$$\mathrm{\int^\frac{\pi}{2}_0x\:sin(x)dx}$$

与1/3情况相比，算法的唯一区别在于求解公式2的方式。如果你仔细观察，你会理解到右边第二项的索引从1开始，而第三项的索引从3开始。因此，循环必须分成两部分，如下所示：

# Summing over remaining n-2 trapezoids
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)
# Left and Right limits
a=0
b=pi/2

# Number of trapezoids
n=6

# Width of trapezoid
h=(b-a)/n

# array of x
# For n trapezoids (no. of nodes=n+1)
x=linspace(a,b,n+1)

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

# Summing over remaining n-2 trapezoids
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法则的算法。