# NumPy Ndarray 分割及數值積分
> 作者：王一哲
> 日期：2020/5/15
## 產生 Ndarray 的方法
### 由 list 轉為 ndarray
首先引入函式庫
```python
import numpy as np
```
可以手動輸入資料儲存為 list ，也可以用 for 迴圈產生 list，再用 np.array 轉為 ndarray，語法如下
```python
arr1 = np.array([i for i in range(10)])
```
也可以改用 np.asarray 轉為 ndarray，語法如下
```python
arr1 = np.asarray([i for i in range(10)])
```
計算結果為
```python
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
```
<br />
### 使用 NumPy 內建函式
NumPy 內建的 **arange** 語法與 Python 預設的 range 很像，語法為
```python
np.arange(首項, 末項, 增量)
```
輸出值不包含末項，第3個參數**增量**可以是小數，例如
```python
arr2 = np.arange(0, 1, 0.1)
```
計算結果為
```python
array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9])
```
<br />
另一個常用的函式為 **linspace**，語法為
```python
np.linspace(首項, 末項, 分割數量)
```
將首項與末項之間平均分割為指定的數量，輸出值包含末項，例如
```python
arr3 = np.linspace(0, 1, 11)
```
計算結果為
```python
array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ])
```
<br />
## 分割 Ndarray 的方法
假設產生的 ndarray 名稱為 arr
```python
arr = np.linspace(0, 1, 11)
```
輸出為
```python
array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ])
```
<br />
ndarray 中的元素索引值從0開始編號，可以用索引值取出指定的元素，因此
```python
arr[1:]
```
輸出不包含索引值為0的元素
```python
array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ])
```
<br />
如果使用
```python
arr[:-1]
```
輸出不包含最後一個元素
```python
array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9])
```
<br />
也可以取出指定間隔的元素
```python
arr[::2]
```
輸出為
```python
array([0. , 0.2, 0.4, 0.6, 0.8, 1. ])
```
<br />
可以在第一個冒號前面加上起始的索引值
```python
arr[1::2]
```
輸出為
```python
array([0.1, 0.3, 0.5, 0.7, 0.9])
```
<br />
## 使用 ndarray 計算數值積分
假設要計算以下式子
$$
\int_{-2}^2 \left[ x^2 + \sin(2 \pi x) \right] dx = \left[ \frac{x^3}{3} - \frac{1}{2 \pi} \cos (2 \pi x) \right]_{-2}^{2} \approx 5.333333333333
$$
<img height="100%" width="100%" src="https://lh5.googleusercontent.com/sf9XCQ8tqDEk-isZSw4SRCDlkacbZYiC1YaFmHaj7u_NFhBLrgRSnOmdrN5pnp0rx3gl2y9KUYLgdFD5jZNuHKHdj4Yz1AEgOVie2LpBnMvNMJkH45xJqfP_HX96HrifBdkomhET" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">定積分</div>
<br />
使用 SymPy 的作法如下
```python=
import sympy as sp
import mpmath
x = sp.symbols('x', commutative=True)
x0, x1 = -2, 2
sym = sp.integrate((x**2 + sp.sin(2*sp.pi*x)), (x, x0, x1), manual=True).evalf()
print("sym = {:.12f}".format(sym))
```
<br />
計算結果為
```python
sym = 5.333333333333
```
<br />
### 長方形法
將函數下方的面積拆成 $N$ 個的長方形，若用長方形的中點計算對應的函數值 $f(x)$，則積分值可以近似為
$$
\int_a^b f(x)dx \approx \Delta x \left( f(x_1) + f(x_2) + \dots + f(x_n) \right)
$$
其中
$$
\Delta x = \frac{b-a}{n}
$$
$$
x_i = a + \left( i - \frac{1}{2} \right) \Delta x ~~~~~ i = 1, 2, 3, \dots, n
$$
<img height="100%" width="100%" src="https://lh6.googleusercontent.com/ruki7kpunaMN6IhDi0ndoczLeY6BzOt6dmMYt5TVESwEf9VQf1xhdUa0_7Of8V9jG0LSqgmg39DyGL-K8E9S--ybJOa2-waIDzsKjqqqJxc_OGCAq6SE0wr124FomlaRjwy907wG" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">長方形法（中點）</div>
<br />
使用 ndarray 計算的方法如下
```python=
import numpy as np
x0, x1, N = -2, 2, 20
xs = np.linspace(x0, x1, N+1)
xm = (xs[1:] + xs[:-1]) / 2
def func(x):
return x**2 + np.sin(2*np.pi*x)
rectM = sum(func(xm)*(x1-x0)/N)
print("N = {:d}, result3 = {:.12f}".format(N, rectM))
```
<br />
計算結果為
```python
N = 20, result3 = 5.320000000000
```
<br />
若用長方形的左側計算對應的函數值 $f(x)$，則積分值可以近似為
$$
\int_a^b f(x)dx \approx \Delta x \left( f(x_1) + f(x_2) + \dots + f(x_n) \right)
$$
其中
$$
\Delta x = \frac{b-a}{n}
$$
$$
x_i = a + (i - 1) \Delta x ~~~~~ i = 1, 2, 3, \dots, n
$$
<img height="100%" width="100%" src="https://lh5.googleusercontent.com/Yis5nPN14HseDj_H4O0zIgv0wLYtPQwAcrbfnGXMvDHFhasROmX1G0y8r3Dx7tUa9qnI1ngiXdjpTy9Qs4IKQ4zxa0h5XUf1ja8TddahR4hDu2VcCufxnLQ0yYKsaVSwc7sTIgBM" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">長方形法（左側）</div>
<br />
使用 ndarray 計算的方法如下
```python=
import numpy as np
x0, x1, N, r = -2, 2, 20, 0
xs = np.linspace(x0, x1, N+1)
interval = (x1-x0)/N
xm = xs[:-1] + interval*r
def func(x):
return x**2 + np.sin(2*np.pi*x)
rectL = sum(func(xm)*interval)
print("N = {:d}, rectL = {:.12f}".format(N, rectL))
```
<br />
計算結果為
```python
N = 20, rectL = 5.360000000000
```
<br />
若用長方形的右側計算對應的函數值 $f(x)$，則積分值可以近似為
$$
\int_a^b f(x)dx \approx \Delta x \left( f(x_1) + f(x_2) + \dots + f(x_n) \right)
$$
其中
$$
\Delta x = \frac{b-a}{n}
$$
$$
x_i = a + i \Delta x ~~~~~ i = 1, 2, 3, \dots, n
$$
<img height="100%" width="100%" src="https://lh6.googleusercontent.com/NfVC9A71pqeV9ANrDAxWQSEHiHUyudlvy6dJH0gb3m7ev5KBIjGIG8muexUuKvEYZ_J-R_3dWm6m2HfZUoAV_oobL9t215cmd92g9AFOG2_IHgm41QiauCCr-pIJcwkR6MMN5R_C" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">長方形法（右側）</div>
<br />
使用 ndarray 計算的方法如下
```python=
import numpy as np
x0, x1, N, r = -2, 2, 20, 0
xs = np.linspace(x0, x1, N+1)
interval = (x1-x0)/N
xm = xs[:-1] + interval*1
def func(x):
return x**2 + np.sin(2*np.pi*x)
rectR = sum(func(xm)*interval)
print("N = {:d}, rectR = {:.12f}".format(N, rectR))
```
<br />
計算結果為
```python
N = 20, rectR = 5.360000000000
```
<br />
### 梯形法
將函數下方的面積拆成 $N$ 個的梯形，若用梯形兩側的點計算對應的函數值 $f(x)$ 作為上底及下底，則積分值可以近似為
$$
\int_a^b f(x)dx \approx \frac{\Delta x}{2} \left( f(x_0) + 2f(x_1) + 2f(x_2) + \dots + 2f(x_{n-1}) + f(x_n) \right)
$$
其中
$$
\Delta x = \frac{b-a}{n}
$$
$$
x_i = a + i \Delta x ~~~~~ i = 0, 1, 2, \dots, n
$$
<img height="100%" width="100%" src="https://lh4.googleusercontent.com/tIbe_S7AbFROcxqQQD7aKpD6quhQ7f54O84VpX3z5OM58WP4zSOzr_5hzVgGeMFI1ET3h1j38ttgFQiQfZxvsxrtPDCWrq8hUdjCCPxSNKOvjRhC2tp7hh1OxBo3BVf6Tl42BRm7" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">梯形法</div>
<br />
使用 ndarray 計算的方法如下
```python=
import numpy as np
x0, x1, N = -2, 2, 20
xs = np.linspace(x0, x1, N+1)
interval = (x1-x0)/N
def func(x):
return x**2 + np.sin(2*np.pi*x)
trape = sum(0.5 * (func(xs[1:]) + func(xs[:-1])) * interval)
print("N = {:d}, trape = {:.12f}".format(N, trape))
```
<br />
計算結果為
```python
N = 20, trape = 5.360000000000
```
<br />
### 辛普森1/3法則
將函數下方的面積拆成 $N$ 個，最上方用二次曲線取近似，則積分值可以近似為
$$
\int_a^b f(x)dx \approx \frac{\Delta x}{3} \left( f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + 2f(x_4) + \dots + 4f(x_{n-1}) + f(x_n) \right)
$$
其中
$$
\Delta x = \frac{b-a}{n}
$$
$$
x_i = a + i \Delta x ~~~~~ i = 0, 1, 2, 3, \dots, n
$$
分割數量必須是偶數，首項、末項係數為1，其餘每2項1組，第1項係數為4，第2項係數為2。使用 ndarray 計算的方法如下
```python=
import numpy as np
x0, x1, N = -2, 2, 20
xs = np.linspace(x0, x1, N+1)
interval = (x1-x0)/N
def func(x):
return x**2 + np.sin(2*np.pi*x)
fs = np.zeros(len(xs))
fs[0] = func(xs[0])
fs[-1] = func(xs[-1])
fs[1:-1:2] = 4*func(xs[1:-1:2])
fs[2:-1:2] = 2*func(xs[2:-1:2])
simpson = sum(fs)*interval/3
print("N = {:d}, simpson = {:.12f}".format(N, simpson))
```
<br />
計算結果為
```python
N = 20, simpson = 5.333333333333
```
<br />
### 辛普森3/8法則
將函數下方的面積拆成 $N$ 個，最上方用二次曲線取近似，則積分值可以近似為
$$
\int_a^b f(x)dx \approx \frac{3 \Delta x}{8} \left( f(x_0) + 3f(x_1) + 3f(x_2) + 2f(x_3) + 3f(x_4) + 3f(x_5) + 2f(x_6) + \dots + 3f(x_{n-2}) + 3f(x_{n-1}) + f(x_n) \right)
$$
其中
$$
\Delta x = \frac{b-a}{n}
$$
$$
x_i = a + i \Delta x ~~~~~ n = 0, 1, 2, 3, \dots, n
$$
分割數量必須是3的倍數，首項、末項係數為1，其餘每3項1組，1、2項係數為3，第3項係數為2。使用 ndarray 計算的方法如下
```python=
import numpy as np
x0, x1, N = -2, 2, 21
xs = np.linspace(x0, x1, N+1)
interval = (x1-x0)/N
def func(x):
return x**2 + np.sin(2*np.pi*x)
fs = np.zeros(len(xs))
fs[0] = func(xs[0])
fs[-1] = func(xs[-1])
fs[1:-1:3] = 3*func(xs[1:-1:3])
fs[2:-1:3] = 3*func(xs[2:-1:3])
fs[3:-1:3] = 2*func(xs[3:-1:3])
simpson2 = sum(fs)*interval*3/8
print("N = {:d}, simpson2 = {:.12f}".format(N, simpson2))
```
<br />
計算結果為
```python
N = 21, simpson2 = 5.333333333333
```
<br />
## 結語
我在好幾年前寫過兩篇關於數值積分的文章：〈[使用 LibreOffice Calc 或 C 語言計算數值積分](https://keejko.blogspot.com/2014/08/blog-post.html)〉、〈[使用 GeoGebra 繪製數值積分圖形](https://keejko.blogspot.com/2018/01/blog-post_30.html)〉，當時的作法比較複雜，但是使用 NumPy ndarray 計算數值積分時連迴圈都不需要用到，而且計算速度很快，處理資料時一定要善用這個工具。
<br />
---
###### tags:`Physics`、`Python`# NumPy Ndarray 分割及數值積分
> 作者：王一哲
> 日期：2020/5/15
## 產生 Ndarray 的方法
### 由 list 轉為 ndarray
首先引入函式庫
```python
import numpy as np
```
可以手動輸入資料儲存為 list ，也可以用 for 迴圈產生 list，再用 np.array 轉為 ndarray，語法如下
```python
arr1 = np.array([i for i in range(10)])
```
計算結果為
```python
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
```
<br />
### 使用 NumPy 內建函式
NumPy 內建的 **arange** 語法與 Python 預設的 range 很像，語法為
```python
np.arange(首項, 末項, 增量)
```
輸出值不包含末項，第3個參數**增量**可以是小數，例如
```python
arr2 = np.arange(0, 1, 0.1)
```
計算結果為
```python
array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9])
```
<br />
另一個常用的函式為 **linspace**，語法為
```python
np.linspace(首項, 末項, 分割數量)
```
將首項與末項之間平均分割為指定的數量，輸出值包含末項，例如
```python
arr3 = np.linspace(0, 1, 11)
```
計算結果為
```python
array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ])
```
<br />
## 分割 Ndarray 的方法
假設產生的 ndarray 名稱為 arr
```python
arr = np.linspace(0, 1, 11)
```
輸出為
```python
array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ])
```
<br />
ndarray 中的元素索引值從0開始編號，可以用索引值取出指定的元素，因此
```python
arr[1:]
```
輸出不包含索引值為0的元素
```python
array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ])
```
<br />
如果使用
```python
arr[:-1]
```
輸出不包含最後一個元素
```python
array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9])
```
<br />
也可以取出指定間隔的元素
```python
arr[::2]
```
輸出為
```python
array([0. , 0.2, 0.4, 0.6, 0.8, 1. ])
```
<br />
可以在第一個冒號前面加上起始的索引值
```python
arr[1::2]
```
輸出為
```python
array([0.1, 0.3, 0.5, 0.7, 0.9])
```
<br />
## 使用 ndarray 計算數值積分
假設要計算以下式子
$$
\int_{-2}^2 \left[ x^2 + \sin(2 \pi x) \right] dx = \left[ \frac{x^3}{3} - \frac{1}{2 \pi} \cos (2 \pi x) \right]_{-2}^{2} \approx 5.333333333333
$$
<img height="100%" width="100%" src="https://lh5.googleusercontent.com/sf9XCQ8tqDEk-isZSw4SRCDlkacbZYiC1YaFmHaj7u_NFhBLrgRSnOmdrN5pnp0rx3gl2y9KUYLgdFD5jZNuHKHdj4Yz1AEgOVie2LpBnMvNMJkH45xJqfP_HX96HrifBdkomhET" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">定積分</div>
<br />
使用 SymPy 的作法如下
```python=
import sympy as sp
import mpmath
x = sp.symbols('x', commutative=True)
x0, x1 = -2, 2
sym = sp.integrate((x**2 + sp.sin(2*sp.pi*x)), (x, x0, x1), manual=True).evalf()
print("sym = {:.12f}".format(sym))
```
<br />
計算結果為
```python
sym = 5.333333333333
```
<br />
### 長方形法
將函數下方的面積拆成 $N$ 個的長方形，若用長方形的中點計算對應的函數值 $f(x)$，則積分值可以近似為
$$
\int_a^b f(x)dx \approx \Delta x \left( f(x_1) + f(x_2) + \dots + f(x_n) \right)
$$
其中
$$
\Delta x = \frac{b-a}{n}
$$
$$
x_i = a + \left( i - \frac{1}{2} \right) \Delta x ~~~~~ i = 1, 2, 3, \dots, n
$$
<img height="100%" width="100%" src="https://lh6.googleusercontent.com/ruki7kpunaMN6IhDi0ndoczLeY6BzOt6dmMYt5TVESwEf9VQf1xhdUa0_7Of8V9jG0LSqgmg39DyGL-K8E9S--ybJOa2-waIDzsKjqqqJxc_OGCAq6SE0wr124FomlaRjwy907wG" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">長方形法（中點）</div>
<br />
使用 ndarray 計算的方法如下
```python=
import numpy as np
x0, x1, N = -2, 2, 20
xs = np.linspace(x0, x1, N+1)
xm = (xs[1:] + xs[:-1]) / 2
def func(x):
return x**2 + np.sin(2*np.pi*x)
rectM = func(xm).sum()*(x1-x0)/N
print("N = {:d}, result3 = {:.12f}".format(N, rectM))
```
<br />
計算結果為
```python
N = 20, result3 = 5.320000000000
```
<br />
若用長方形的左側計算對應的函數值 $f(x)$，則積分值可以近似為
$$
\int_a^b f(x)dx \approx \Delta x \left( f(x_1) + f(x_2) + \dots + f(x_n) \right)
$$
其中
$$
\Delta x = \frac{b-a}{n}
$$
$$
x_i = a + (i - 1) \Delta x ~~~~~ i = 1, 2, 3, \dots, n
$$
<img height="100%" width="100%" src="https://lh5.googleusercontent.com/Yis5nPN14HseDj_H4O0zIgv0wLYtPQwAcrbfnGXMvDHFhasROmX1G0y8r3Dx7tUa9qnI1ngiXdjpTy9Qs4IKQ4zxa0h5XUf1ja8TddahR4hDu2VcCufxnLQ0yYKsaVSwc7sTIgBM" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">長方形法（左側）</div>
<br />
使用 ndarray 計算的方法如下
```python=
import numpy as np
x0, x1, N, r = -2, 2, 20, 0
xs = np.linspace(x0, x1, N+1)
interval = (x1-x0)/N
xm = xs[:-1] + interval*r
def func(x):
return x**2 + np.sin(2*np.pi*x)
rectL = func(xm).sum()*interval
print("N = {:d}, rectL = {:.12f}".format(N, rectL))
```
<br />
計算結果為
```python
N = 20, rectL = 5.360000000000
```
<br />
若用長方形的右側計算對應的函數值 $f(x)$，則積分值可以近似為
$$
\int_a^b f(x)dx \approx \Delta x \left( f(x_1) + f(x_2) + \dots + f(x_n) \right)
$$
其中
$$
\Delta x = \frac{b-a}{n}
$$
$$
x_i = a + i \Delta x ~~~~~ i = 1, 2, 3, \dots, n
$$
<img height="100%" width="100%" src="https://lh6.googleusercontent.com/NfVC9A71pqeV9ANrDAxWQSEHiHUyudlvy6dJH0gb3m7ev5KBIjGIG8muexUuKvEYZ_J-R_3dWm6m2HfZUoAV_oobL9t215cmd92g9AFOG2_IHgm41QiauCCr-pIJcwkR6MMN5R_C" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">長方形法（右側）</div>
<br />
使用 ndarray 計算的方法如下
```python=
import numpy as np
x0, x1, N, r = -2, 2, 20, 0
xs = np.linspace(x0, x1, N+1)
interval = (x1-x0)/N
xm = xs[:-1] + interval*1
def func(x):
return x**2 + np.sin(2*np.pi*x)
rectR = func(xm).sum()*interval
print("N = {:d}, rectR = {:.12f}".format(N, rectR))
```
<br />
計算結果為
```python
N = 20, rectR = 5.360000000000
```
<br />
### 梯形法
將函數下方的面積拆成 $N$ 個的梯形，若用梯形兩側的點計算對應的函數值 $f(x)$ 作為上底及下底，則積分值可以近似為
$$
\int_a^b f(x)dx \approx \frac{\Delta x}{2} \left( f(x_0) + 2f(x_1) + 2f(x_2) + \dots + 2f(x_{n-1}) + f(x_n) \right)
$$
其中
$$
\Delta x = \frac{b-a}{n}
$$
$$
x_i = a + i \Delta x ~~~~~ i = 0, 1, 2, \dots, n
$$
<img height="100%" width="100%" src="https://lh4.googleusercontent.com/tIbe_S7AbFROcxqQQD7aKpD6quhQ7f54O84VpX3z5OM58WP4zSOzr_5hzVgGeMFI1ET3h1j38ttgFQiQfZxvsxrtPDCWrq8hUdjCCPxSNKOvjRhC2tp7hh1OxBo3BVf6Tl42BRm7" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">梯形法</div>
<br />
使用 ndarray 計算的方法如下
```python=
import numpy as np
x0, x1, N = -2, 2, 20
xs = np.linspace(x0, x1, N+1)
interval = (x1-x0)/N
def func(x):
return x**2 + np.sin(2*np.pi*x)
trape = 0.5 * (func(xs[1:]) + func(xs[:-1])).sum() * interval
print("N = {:d}, trape = {:.12f}".format(N, trape))
```
<br />
計算結果為
```python
N = 20, trape = 5.360000000000
```
<br />
### 辛普森1/3法則
將函數下方的面積拆成 $N$ 個，最上方用二次曲線取近似，則積分值可以近似為
$$
\int_a^b f(x)dx \approx \frac{\Delta x}{3} \left( f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + 2f(x_4) + \dots + 4f(x_{n-1}) + f(x_n) \right)
$$
其中
$$
\Delta x = \frac{b-a}{n}
$$
$$
x_i = a + i \Delta x ~~~~~ i = 0, 1, 2, 3, \dots, n
$$
分割數量必須是偶數，首項、末項係數為1，其餘每2項1組，第1項係數為4，第2項係數為2。使用 ndarray 計算的方法如下
```python=
import numpy as np
x0, x1, N = -2, 2, 20
xs = np.linspace(x0, x1, N+1)
interval = (x1-x0)/N
def func(x):
return x**2 + np.sin(2*np.pi*x)
fs = np.zeros(len(xs))
fs[0] = func(xs[0])
fs[-1] = func(xs[-1])
fs[1:-1:2] = 4*func(xs[1:-1:2])
fs[2:-1:2] = 2*func(xs[2:-1:2])
simpson = fs.sum()*interval/3
print("N = {:d}, simpson = {:.12f}".format(N, simpson))
```
<br />
計算結果為
```python
N = 20, simpson = 5.333333333333
```
<br />
### 辛普森3/8法則
將函數下方的面積拆成 $N$ 個，最上方用二次曲線取近似，則積分值可以近似為
$$
\int_a^b f(x)dx \approx \frac{3 \Delta x}{8} \left( f(x_0) + 3f(x_1) + 3f(x_2) + 2f(x_3) + 3f(x_4) + 3f(x_5) + 2f(x_6) + \dots + 3f(x_{n-2}) + 3f(x_{n-1}) + f(x_n) \right)
$$
其中
$$
\Delta x = \frac{b-a}{n}
$$
$$
x_i = a + i \Delta x ~~~~~ n = 0, 1, 2, 3, \dots, n
$$
分割數量必須是3的倍數，首項、末項係數為1，其餘每3項1組，1、2項係數為3，第3項係數為2。使用 ndarray 計算的方法如下
```python=
import numpy as np
x0, x1, N = -2, 2, 21
xs = np.linspace(x0, x1, N+1)
interval = (x1-x0)/N
def func(x):
return x**2 + np.sin(2*np.pi*x)
fs = np.zeros(len(xs))
fs[0] = func(xs[0])
fs[-1] = func(xs[-1])
fs[1:-1:3] = 3*func(xs[1:-1:3])
fs[2:-1:3] = 3*func(xs[2:-1:3])
fs[3:-1:3] = 2*func(xs[3:-1:3])
simpson2 = fs.sum()*interval*3/8
print("N = {:d}, simpson2 = {:.12f}".format(N, simpson2))
```
<br />
計算結果為
```python
N = 21, simpson2 = 5.333333333333
```
<br />
## 結語
我在好幾年前寫過兩篇關於數值積分的文章：〈[使用 LibreOffice Calc 或 C 語言計算數值積分](https://keejko.blogspot.com/2014/08/blog-post.html)〉、〈[使用 GeoGebra 繪製數值積分圖形](https://keejko.blogspot.com/2018/01/blog-post_30.html)〉，當時的作法比較複雜，但是使用 NumPy ndarray 計算數值積分時連迴圈都不需要用到，而且計算速度很快，處理資料時一定要善用這個工具。
<br />
---
###### tags:`Physics`、`Python`
Sign in with Wallet
Connect another wallet