NumPy Ndarray 分割及數值積分

作者：王一哲
日期：2020/5/15

產生 Ndarray 的方法

由 list 轉為 ndarray

首先引入函式庫

import numpy as np

可以手動輸入資料儲存為 list ，也可以用 for 迴圈產生 list，再用 np.array 轉為 ndarray，語法如下

arr1 = np.array([i for i in range(10)])

也可以改用 np.asarray 轉為 ndarray，語法如下

arr1 = np.asarray([i for i in range(10)])

計算結果為

array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

使用 NumPy 內建函式

NumPy 內建的 arange 語法與 Python 預設的 range 很像，語法為

np.arange(首項, 末項, 增量)

輸出值不包含末項，第3個參數增量可以是小數，例如

arr2 = np.arange(0, 1, 0.1)

計算結果為

array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9])

另一個常用的函式為 linspace，語法為

np.linspace(首項, 末項, 分割數量)

將首項與末項之間平均分割為指定的數量，輸出值包含末項，例如

arr3 = np.linspace(0, 1, 11)

計算結果為

array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ])

分割 Ndarray 的方法

假設產生的 ndarray 名稱為 arr

arr = np.linspace(0, 1, 11)

輸出為

array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ])

ndarray 中的元素索引值從0開始編號，可以用索引值取出指定的元素，因此

arr[1:]

輸出不包含索引值為0的元素

array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ])

如果使用

arr[:-1]

輸出不包含最後一個元素

array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9])

也可以取出指定間隔的元素

arr[::2]

輸出為

array([0. , 0.2, 0.4, 0.6, 0.8, 1. ])

可以在第一個冒號前面加上起始的索引值

arr[1::2]

輸出為

array([0.1, 0.3, 0.5, 0.7, 0.9])

使用 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
$$

定積分

使用 SymPy 的作法如下

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))

計算結果為

sym = 5.333333333333

長方形法

將函數下方的面積拆成 $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
$$

長方形法（中點）

使用 ndarray 計算的方法如下

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))

計算結果為

N = 20, result3 = 5.320000000000

若用長方形的左側計算對應的函數值 $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
$$

長方形法（左側）

使用 ndarray 計算的方法如下

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))

計算結果為

N = 20, rectL = 5.360000000000

若用長方形的右側計算對應的函數值 $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
$$

長方形法（右側）

使用 ndarray 計算的方法如下

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))

計算結果為

N = 20, rectR = 5.360000000000

梯形法

將函數下方的面積拆成 $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
$$

梯形法

使用 ndarray 計算的方法如下

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))

計算結果為

N = 20, trape = 5.360000000000

辛普森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 計算的方法如下

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))

計算結果為

N = 20, simpson = 5.333333333333

辛普森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 計算的方法如下

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))

計算結果為

N = 21, simpson2 = 5.333333333333

結語

我在好幾年前寫過兩篇關於數值積分的文章：〈使用 LibreOffice Calc 或 C 語言計算數值積分〉、〈使用 GeoGebra 繪製數值積分圖形〉，當時的作法比較複雜，但是使用 NumPy ndarray 計算數值積分時連迴圈都不需要用到，而且計算速度很快，處理資料時一定要善用這個工具。

---

HackMD 版本連結：https://hackmd.io/@yizhewang/HJn7vTsqU