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[x^2+\sin (2\pi x)]dx={[\frac{x^3}{3}-\frac{1}{2\pi }\cos (2\pi x)]}_{-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 \mathrm{\Delta }x(f(x_1)+f(x_2)+\cdots +f(x_n))$$

其中

$$\mathrm{\Delta }x=\frac{b-a}{n}$$

$$x_i=a+(i-\frac{1}{2})\mathrm{\Delta }xi=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 = sum(func(xm)*(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 \mathrm{\Delta }x(f(x_1)+f(x_2)+\cdots +f(x_n))$$

其中

$$\mathrm{\Delta }x=\frac{b-a}{n}$$

$$x_i=a+(i-1)\mathrm{\Delta }xi=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 = sum(func(xm)*interval)

print("N = {:d}, rectL = {:.12f}".format(N, rectL))

計算結果為

N = 20, rectL = 5.360000000000

若用長方形的右側計算對應的函數值

$f(x)$，則積分值可以近似為

$${\int }_a^{b}f(x)dx\approx \mathrm{\Delta }x(f(x_1)+f(x_2)+\cdots +f(x_n))$$

其中

$$\mathrm{\Delta }x=\frac{b-a}{n}$$

$$x_i=a+i\mathrm{\Delta }xi=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 = sum(func(xm)*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{\mathrm{\Delta }x}{2}(f(x_0)+2f(x_1)+2f(x_2)+\cdots +2f(x_{n-1})+f(x_n))$$

其中

$$\mathrm{\Delta }x=\frac{b-a}{n}$$

$$x_i=a+i\mathrm{\Delta }xi=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 = sum(0.5 * (func(xs[1:]) + func(xs[:-1])) * 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{\mathrm{\Delta }x}{3}(f(x_0)+4f(x_1)+2f(x_2)+4f(x_3)+2f(x_4)+\cdots +4f(x_{n-1})+f(x_n))$$

其中

$$\mathrm{\Delta }x=\frac{b-a}{n}$$

$$x_i=a+i\mathrm{\Delta }xi=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 = sum(fs)*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\mathrm{\Delta }x}{8}(f(x_0)+3f(x_1)+3f(x_2)+2f(x_3)+3f(x_4)+3f(x_5)+2f(x_6)+\cdots +3f(x_{n-2})+3f(x_{n-1})+f(x_n))$$

其中

$$\mathrm{\Delta }x=\frac{b-a}{n}$$

$$x_i=a+i\mathrm{\Delta }xn=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 = sum(fs)*interval*3/8

print("N = {:d}, simpson2 = {:.12f}".format(N, simpson2))

計算結果為

N = 21, simpson2 = 5.333333333333

結語

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

---

tags:Physics、Python# 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)])

計算結果為

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[x^2+\sin (2\pi x)]dx={[\frac{x^3}{3}-\frac{1}{2\pi }\cos (2\pi x)]}_{-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 \mathrm{\Delta }x(f(x_1)+f(x_2)+\cdots +f(x_n))$$

其中

$$\mathrm{\Delta }x=\frac{b-a}{n}$$

$$x_i=a+(i-\frac{1}{2})\mathrm{\Delta }xi=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 \mathrm{\Delta }x(f(x_1)+f(x_2)+\cdots +f(x_n))$$

其中

$$\mathrm{\Delta }x=\frac{b-a}{n}$$

$$x_i=a+(i-1)\mathrm{\Delta }xi=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 \mathrm{\Delta }x(f(x_1)+f(x_2)+\cdots +f(x_n))$$

其中

$$\mathrm{\Delta }x=\frac{b-a}{n}$$

$$x_i=a+i\mathrm{\Delta }xi=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{\mathrm{\Delta }x}{2}(f(x_0)+2f(x_1)+2f(x_2)+\cdots +2f(x_{n-1})+f(x_n))$$

其中

$$\mathrm{\Delta }x=\frac{b-a}{n}$$

$$x_i=a+i\mathrm{\Delta }xi=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{\mathrm{\Delta }x}{3}(f(x_0)+4f(x_1)+2f(x_2)+4f(x_3)+2f(x_4)+\cdots +4f(x_{n-1})+f(x_n))$$

其中

$$\mathrm{\Delta }x=\frac{b-a}{n}$$

$$x_i=a+i\mathrm{\Delta }xi=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\mathrm{\Delta }x}{8}(f(x_0)+3f(x_1)+3f(x_2)+2f(x_3)+3f(x_4)+3f(x_5)+2f(x_6)+\cdots +3f(x_{n-2})+3f(x_{n-1})+f(x_n))$$

其中

$$\mathrm{\Delta }x=\frac{b-a}{n}$$

$$x_i=a+i\mathrm{\Delta }xn=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 計算數值積分時連迴圈都不需要用到，而且計算速度很快，處理資料時一定要善用這個工具。

---

tags:Physics、Python