前言 #
在科學計算、機器學習、數據分析的世界裡，我們常常會遇到「積分」這個老朋友。但現實很殘酷：大部分的積分根本沒有解析解（也就是沒有漂亮的公式可以直接算）！
這時候就需要用「數值積分」來拯救我們了。數值積分的核心概念很簡單：既然算不出精確值，那就用「逼近」的方式來算出一個夠準確的近似值！
這篇文章拍拍君會帶你：
- 了解什麼是數值積分
- 認識常見的數值積分演算法
- 用 Python 的 NumPy 和 SciPy 實際操作
準備好了嗎？Let’s go！
為什麼需要數值積分？ #
想像一下，你要計算這個積分：
$$ \int_{-1}^1 \sqrt{1-x^2} , dx $$
這個積分其實有解析解（答案是 \(\pi/2\) )，因為這是半圓的面積），但如果遇到更複雜的函數，像是：
$$ \int_0^1 e^{-x^2} , dx $$
這種積分就沒有簡單的公式可以算了！這時候數值積分就派上用場。
數值積分的基本想法：把積分想像成「曲線下方的面積」，然後用一堆小矩形、梯形或其他形狀去逼近這個面積。
常見的數值積分演算法 #
1. 矩形法（Riemann Sum） #
最直覺的方法：把積分區間切成很多小段，每一段用矩形來近似。
$$ \int_a^b f(x) , dx \approx \sum_{i=0}^{n-1} f(x_i) \cdot \Delta x $$
其中 \(\Delta x = \frac{b-a}{n}\) 是每個小區間的寬度。
優點：超級簡單好理解 缺點：精度不高，需要切很多段才準確
2. 梯形法（Trapezoidal Rule） #
比矩形法更聰明一點：用梯形來逼近曲線。
$$ \int_a^b f(x) , dx \approx \frac{\Delta x}{2} \left[ f(x_0) + 2f(x_1) + 2f(x_2) + \cdots + 2f(x_{n-1}) + f(x_n) \right] $$
優點：比矩形法準確 缺點：還是需要蠻多分割點
3. 辛普森法（Simpson’s Rule） #
更進階的方法：用二次曲線（拋物線）來逼近每一小段的函數。
$$ \int_a^b f(x) , dx \approx \frac{\Delta x}{3} \left[ f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + \cdots + f(x_n) \right] $$
優點：精度高，分割點不用太多 缺點：需要偶數個區間（奇數個點）
4. 高斯積分（Gaussian Quadrature） #
最強大的方法之一：不是均勻分割區間，而是根據數學上的最佳點來取樣。
優點：超高精度，特別適合多項式函數 缺點：實作比較複雜（好在 SciPy 都幫我們包好了！）
Python 實作：NumPy 版本 #
先來看看如何用 NumPy 手刻一些基本的數值積分方法。
矩形法實作 #
import numpy as np
def rectangle_rule(f, a, b, n=1000):
"""
矩形法數值積分
f: 要積分的函數
a, b: 積分區間 [a, b]
n: 分割數量
"""
x = np.linspace(a, b, n, endpoint=False)
dx = (b - a) / n
return np.sum(f(x) * dx)
# 測試：計算 ∫₋₁¹ √(1-x²) dx（半圓面積，應該是 π/2 ≈ 1.5708）
f = lambda x: np.sqrt(1 - x**2)
result = rectangle_rule(f, -1, 1, n=10000)
print(f"矩形法結果: {result:.6f}")
print(f"真實值 π/2: {np.pi/2:.6f}")
梯形法實作 #
def trapezoidal_rule(f, a, b, n=1000):
"""
梯形法數值積分
"""
x = np.linspace(a, b, n+1)
y = f(x)
dx = (b - a) / n
return dx * (np.sum(y) - 0.5 * (y[0] + y[-1]))
result = trapezoidal_rule(f, -1, 1, n=1000)
print(f"梯形法結果: {result:.6f}")
print(f"真實值 π/2: {np.pi/2:.6f}")
辛普森法實作 #
def simpson_rule(f, a, b, n=1000):
"""
辛普森法數值積分（n 必須是偶數）
"""
if n % 2 == 1:
n += 1 # 確保是偶數
x = np.linspace(a, b, n+1)
y = f(x)
dx = (b - a) / n
return dx / 3 * (y[0] + y[-1] + 4*np.sum(y[1:-1:2]) + 2*np.sum(y[2:-1:2]))
result = simpson_rule(f, -1, 1, n=1000)
print(f"辛普森法結果: {result:.6f}")
print(f"真實值 π/2: {np.pi/2:.6f}")
NumPy 的 trapz 方法 #
NumPy 本身就有內建梯形法！
x = np.linspace(-1, 1, 1000)
y = np.sqrt(1 - x**2)
result = np.trapz(y, x)
print(f"np.trapz 結果: {result:.6f}")
Python 實作：SciPy 版本（推薦！） #
SciPy 的 scipy.integrate 模組提供了更強大、更精確的積分工具。
1. quad：最常用的數值積分 #
from scipy.integrate import quad
# 定義函數
f = lambda x: np.sqrt(1 - x**2)
# 使用 quad 進行積分
result, error = quad(f, -1, 1)
print(f"SciPy quad 結果: {result:.10f}")
print(f"估計誤差: {error:.2e}")
print(f"真實值 π/2: {np.pi/2:.10f}")
優點：
- 自動選擇最佳演算法（通常是高斯積分）
- 自動調整精度
- 會回傳誤差估計值
2. 更複雜的例子：無解析解的積分 #
# 計算 ∫₀¹ e^(-x²) dx（高斯誤差函數相關）
g = lambda x: np.exp(-x**2)
result, error = quad(g, 0, 1)
print(f"結果: {result:.10f}")
print(f"誤差: {error:.2e}")
3. 無窮區間積分 #
# 計算 ∫₀^∞ e^(-x) dx = 1
h = lambda x: np.exp(-x)
result, error = quad(h, 0, np.inf)
print(f"結果: {result:.10f}")
print(f"理論值: 1.0")
4. 多重積分：dblquad 和 tplquad #
from scipy.integrate import dblquad
# 計算二重積分 ∫₀¹ ∫₀¹ (x² + y²) dy dx
f = lambda y, x: x**2 + y**2
result, error = dblquad(f, 0, 1, 0, 1)
print(f"二重積分結果: {result:.6f}")
注意：dblquad 的參數順序是 (f, x_lower, x_upper, y_lower, y_upper)，而函數 f 的參數順序是 f(y, x)（y 在前！）
5. 向量化積分：適合處理多筆資料 #
from scipy.integrate import cumulative_trapezoid
# 累積積分（例如計算速度 → 位移）
t = np.linspace(0, 10, 100)
v = np.sin(t) # 速度函數
s = cumulative_trapezoid(v, t, initial=0) # 位移（積分後）
import matplotlib.pyplot as plt
plt.plot(t, v, label='速度')
plt.plot(t, s, label='位移（積分）')
plt.legend()
plt.show()
各方法精度比較 #
來實際測試一下各種方法的精度：
from scipy.integrate import quad
import numpy as np
# 測試函數：半圓面積
f = lambda x: np.sqrt(1 - x**2)
exact = np.pi / 2
# 不同方法的結果
methods = {
"矩形法 (n=100)": rectangle_rule(f, -1, 1, 100),
"矩形法 (n=1000)": rectangle_rule(f, -1, 1, 1000),
"梯形法 (n=100)": trapezoidal_rule(f, -1, 1, 100),
"梯形法 (n=1000)": trapezoidal_rule(f, -1, 1, 1000),
"辛普森法 (n=100)": simpson_rule(f, -1, 1, 100),
"辛普森法 (n=1000)": simpson_rule(f, -1, 1, 1000),
"SciPy quad": quad(f, -1, 1)[0],
}
print(f"真實值: {exact:.10f}\n")
for name, result in methods.items():
error = abs(result - exact)
print(f"{name:20s}: {result:.10f} (誤差: {error:.2e})")
預期結果：
- 矩形法誤差最大
- 辛普森法比梯形法好很多
- SciPy quad 精度最高（通常到機器精度等級！）
小技巧與注意事項 #
1. 處理奇異點（Singularities） #
如果函數在積分區間內有奇異點（例如 \(\frac{1}{x}\) 在 \(x=0\)），可以用 quad 的 points 參數：
from scipy.integrate import quad
# ∫₋₁¹ 1/x² dx（在 x=0 有奇異點）
f = lambda x: 1 / x**2 if x != 0 else np.inf
result, error = quad(f, -1, 1, points=[0]) # 告訴 quad 在 x=0 要特別處理
2. 選擇適當的分割數 #
- 矩形法、梯形法：至少 1000 個點以上
- 辛普森法：通常 100-1000 個點就很準
- SciPy quad：不用管！它會自動調整
3. 處理振盪函數 #
# 高頻振盪函數
f = lambda x: np.sin(100*x)
result, error = quad(f, 0, 1, limit=100) # 增加 limit 讓它多試幾次
print(f"結果: {result:.6f}")
結語 #
數值積分是科學計算的基石之一，從物理模擬到機器學習都會用到。這篇文章拍拍君帶你從最基本的矩形法，一路到強大的 SciPy 積分工具。
實務建議：
- 學習用：手刻矩形法、梯形法、辛普森法，理解原理
- 工作用：直接用
scipy.integrate.quad，省時省力又準確
數值積分就像是數學世界的「近似藝術」，雖然得不到完美答案，但只要方法對了，就能無限逼近真相！
下次遇到積不出來的函數，別慌張，打開 Python 跑個 quad 就對了！