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程式碼。