日期：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
沒有留言:
張貼留言