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法則的演算法。