第2章: 積分の基礎と数値積分

Fundamentals of Integration and Numerical Integration

基礎数理道場 > 微積分とベクトル解析入門 > 第2章

2.1 定積分と不定積分

📐 定義: 定積分
関数 f(x) の区間 [a, b] における定積分は、リーマン和の極限として定義されます： $$\int_a^b f(x) dx = \lim_{n \to \infty} \sum_{i=1}^n f(x_i^*) \Delta x$$ 幾何学的には、曲線 y = f(x) と x 軸で囲まれた面積を表します。

積分は「累積量」を計算する操作です。材料科学では熱量計算、プロセス工学では反応物の総量、 機械学習では確率分布の正規化定数などに使われます。

💻 コード例1: リーマン和による定積分の近似

import numpy as np import matplotlib.pyplot as plt def f(x): """被積分関数: f(x) = x^2""" return x**2 # リーマン和による積分の近似 def riemann_sum(f, a, b, n, method='midpoint'): """ method: 'left' (左端), 'right' (右端), 'midpoint' (中点) """ dx = (b - a) / n x = np.linspace(a, b, n+1) if method == 'left': x_sample = x[:-1] elif method == 'right': x_sample = x[1:] else: # midpoint x_sample = (x[:-1] + x[1:]) / 2 return np.sum(f(x_sample) * dx) # ∫₀¹ x² dx = 1/3 を計算 a, b = 0, 1 exact_value = 1/3 n_values = [4, 10, 50, 100, 500] print("リーマン和による積分近似:") print(f"解析解: ∫₀¹ x² dx = {exact_value:.10f}\n") for n in n_values: approx = riemann_sum(f, a, b, n, 'midpoint') error = abs(approx - exact_value) print(f"n={n:3d}: 近似値 = {approx:.10f}, 誤差 = {error:.2e}") # 可視化: n=10 の場合の中点リーマン和 n = 10 x = np.linspace(a, b, 1000) dx = (b - a) / n x_rect = np.linspace(a, b, n+1) x_mid = (x_rect[:-1] + x_rect[1:]) / 2 fig, ax = plt.subplots(figsize=(10, 6)) ax.plot(x, f(x), 'b-', linewidth=2, label='f(x) = x²') ax.fill_between(x, 0, f(x), alpha=0.2) # 長方形を描画 for i in range(n): height = f(x_mid[i]) rect = plt.Rectangle((x_rect[i], 0), dx, height, edgecolor='red', facecolor='none', linewidth=1.5) ax.add_patch(rect) ax.plot([x_mid[i]], [height], 'ro', markersize=6) ax.set_xlabel('x', fontsize=12) ax.set_ylabel('f(x)', fontsize=12) ax.set_title(f'中点リーマン和 (n={n})', fontsize=14) ax.legend() ax.grid(True, alpha=0.3) plt.tight_layout() plt.show()
📐 定理: 微積分の基本定理
f(x) が連続関数で F(x) が f(x) の原始関数（不定積分）のとき： $$\int_a^b f(x) dx = F(b) - F(a)$$ この定理により、微分と積分が互いに逆操作であることがわかります。

💻 コード例2: SymPyによる記号積分

import sympy as sp x = sp.Symbol('x') # 不定積分（原始関数） functions = [ x**2, sp.exp(x), sp.sin(x), 1/x, x * sp.exp(x) ] print("不定積分の例:") for func in functions: integral = sp.integrate(func, x) print(f"∫ {func} dx = {integral} + C") print("\n定積分の例:") # ∫₀¹ x² dx result = sp.integrate(x**2, (x, 0, 1)) print(f"∫₀¹ x² dx = {result}") # ∫₀^π sin(x) dx result = sp.integrate(sp.sin(x), (x, 0, sp.pi)) print(f"∫₀^π sin(x) dx = {result}")

2.2 台形則による数値積分

台形則は、曲線を直線で近似して面積を計算する最も基本的な数値積分法です。

💻 コード例3: 台形則の実装

def trapezoidal_rule(f, a, b, n): """ 台形則による数値積分 精度: O(h²) where h = (b-a)/n """ x = np.linspace(a, b, n+1) y = f(x) h = (b - a) / n # 台形則の公式: h * [f(a)/2 + f(x₁) + ... + f(xₙ₋₁) + f(b)/2] integral = h * (y[0]/2 + np.sum(y[1:-1]) + y[-1]/2) return integral # テスト: ∫₀^π sin(x) dx = 2 f = np.sin a, b = 0, np.pi exact = 2.0 print("台形則による積分:") print(f"解析解: ∫₀^π sin(x) dx = {exact}\n") for n in [10, 50, 100, 500]: approx = trapezoidal_rule(f, a, b, n) error = abs(approx - exact) print(f"n={n:3d}: 近似値 = {approx:.10f}, 誤差 = {error:.2e}") # 収束性の確認 n_values = np.array([10, 20, 50, 100, 200, 500, 1000]) errors = [] for n in n_values: approx = trapezoidal_rule(f, a, b, n) errors.append(abs(approx - exact)) plt.figure(figsize=(8, 6)) plt.loglog(n_values, errors, 'o-', linewidth=2, markersize=8, label='台形則の誤差') plt.loglog(n_values, 1/n_values**2, '--', label='O(1/n²)', alpha=0.5) plt.xlabel('分割数 n', fontsize=12) plt.ylabel('絶対誤差', fontsize=12) plt.title('台形則の収束性', fontsize=14) plt.legend() plt.grid(True, alpha=0.3) plt.show()

2.3 Simpson則による高精度積分

Simpson則は、曲線を2次関数（放物線）で近似する方法で、台形則よりも高精度です。

💻 コード例4: Simpson則の実装

def simpson_rule(f, a, b, n): """ Simpson則による数値積分 n は偶数である必要がある 精度: O(h⁴) where h = (b-a)/n """ if n % 2 != 0: raise ValueError("n must be even for Simpson's rule") x = np.linspace(a, b, n+1) y = f(x) h = (b - a) / n # Simpson則の公式 integral = h/3 * (y[0] + 4*np.sum(y[1:-1:2]) + 2*np.sum(y[2:-1:2]) + y[-1]) return integral # 台形則とSimpson則の比較 print("台形則 vs Simpson則:") print(f"解析解: ∫₀^π sin(x) dx = {exact}\n") n = 100 trap_result = trapezoidal_rule(f, a, b, n) simp_result = simpson_rule(f, a, b, n) print(f"n={n}:") print(f"台形則: {trap_result:.10f}, 誤差 = {abs(trap_result - exact):.2e}") print(f"Simpson則: {simp_result:.10f}, 誤差 = {abs(simp_result - exact):.2e}") # 精度比較の可視化 n_values = np.array([10, 20, 50, 100, 200, 500]) errors_trap = [] errors_simp = [] for n in n_values: errors_trap.append(abs(trapezoidal_rule(f, a, b, n) - exact)) errors_simp.append(abs(simpson_rule(f, a, b, n if n%2==0 else n+1) - exact)) plt.figure(figsize=(10, 6)) plt.loglog(n_values, errors_trap, 'o-', linewidth=2, markersize=8, label='台形則') plt.loglog(n_values, errors_simp, 's-', linewidth=2, markersize=8, label='Simpson則') plt.loglog(n_values, 1/n_values**2, '--', label='O(1/n²)', alpha=0.5) plt.loglog(n_values, 1/n_values**4, '--', label='O(1/n⁴)', alpha=0.5) plt.xlabel('分割数 n', fontsize=12) plt.ylabel('絶対誤差', fontsize=12) plt.title('数値積分法の精度比較', fontsize=14) plt.legend() plt.grid(True, alpha=0.3) plt.show()

2.4 SciPyによる高度な数値積分

SciPyライブラリは、適応的刻み幅調整や高次の求積法を実装した高性能な積分関数を提供します。

💻 コード例5: SciPyの数値積分関数

from scipy import integrate # 様々な積分法の比較 def test_function(x): """振動の激しい関数""" return np.sin(x) * np.exp(-x/10) a, b = 0, 20 exact_value = integrate.quad(test_function, a, b, limit=1000)[0] # 高精度参照値 print("SciPyの数値積分:") print(f"参照値 (quad): {exact_value:.10f}\n") # quad: 適応的Gauss-Kronrod求積法 result, error = integrate.quad(test_function, a, b) print(f"quad (適応的): {result:.10f}, 推定誤差 = {error:.2e}") # fixed_quad: 固定次数Gauss-Legendre求積法 result, _ = integrate.fixed_quad(test_function, a, b, n=50) print(f"fixed_quad (n=50): {result:.10f}, 誤差 = {abs(result-exact_value):.2e}") # romberg: Romberg積分（Richardson外挿の拡張） result = integrate.romberg(test_function, a, b, show=False) print(f"romberg: {result:.10f}, 誤差 = {abs(result-exact_value):.2e}") # 可視化 x = np.linspace(a, b, 1000) y = test_function(x) plt.figure(figsize=(10, 6)) plt.plot(x, y, linewidth=2) plt.fill_between(x, 0, y, alpha=0.3) plt.xlabel('x', fontsize=12) plt.ylabel('f(x)', fontsize=12) plt.title('被積分関数: f(x) = sin(x)·exp(-x/10)', fontsize=14) plt.grid(True, alpha=0.3) plt.axhline(y=0, color='k', linewidth=0.5) plt.show()

2.5 材料科学への応用: 熱容量の計算

🔬 応用例: 温度 T₁ から T₂ まで物質を加熱するのに必要な熱量 Q は、 熱容量 C_p(T) を温度で積分して求められます： $$Q = \int_{T_1}^{T_2} C_p(T) dT$$

💻 コード例6: 熱量計算

# 銅の熱容量の温度依存性（実験式） def copper_heat_capacity(T): """ 銅の定圧熱容量 Cp [J/(mol·K)] T: 温度 [K] """ # Shomate equation (簡略版) A, B, C = 22.64, 6.28e-3, 0 return A + B*T + C*T**2 # 300K から 600K まで加熱するのに必要な熱量 T1, T2 = 300, 600 # 数値積分 Q_numerical, error = integrate.quad(copper_heat_capacity, T1, T2) # 解析積分（検証用） def Q_analytical(T1, T2): A, B = 22.64, 6.28e-3 return A*(T2-T1) + B/2*(T2**2-T1**2) Q_exact = Q_analytical(T1, T2) print(f"銅 1 mol を {T1}K → {T2}K に加熱する熱量:") print(f"数値積分: Q = {Q_numerical:.2f} J/mol (推定誤差: {error:.2e})") print(f"解析解: Q = {Q_exact:.2f} J/mol") print(f"差: {abs(Q_numerical - Q_exact):.2e} J/mol") # Cp(T) の可視化 T = np.linspace(200, 800, 100) Cp = copper_heat_capacity(T) fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5)) # 左図: 熱容量の温度依存性 ax1.plot(T, Cp, linewidth=2) ax1.axvspan(T1, T2, alpha=0.2, color='orange', label='積分区間') ax1.set_xlabel('温度 T (K)', fontsize=12) ax1.set_ylabel('熱容量 Cp (J/(mol·K))', fontsize=12) ax1.set_title('銅の熱容量', fontsize=14) ax1.legend() ax1.grid(True, alpha=0.3) # 右図: 熱量の累積 T_range = np.linspace(300, 800, 100) Q_cumulative = [integrate.quad(copper_heat_capacity, 300, T)[0] for T in T_range] ax2.plot(T_range, Q_cumulative, linewidth=2) ax2.axhline(y=Q_numerical, color='red', linestyle='--', label=f'Q(300→600K) = {Q_numerical:.0f} J/mol') ax2.set_xlabel('温度 T (K)', fontsize=12) ax2.set_ylabel('累積熱量 (J/mol)', fontsize=12) ax2.set_title('300Kからの累積加熱量', fontsize=14) ax2.legend() ax2.grid(True, alpha=0.3) plt.tight_layout() plt.show()

2.6 広義積分と特異点の扱い

📝 注意: 被積分関数が無限大に発散する点（特異点）や、 積分区間が無限大の場合（広義積分）は、特別な注意が必要です。

💻 コード例7: 広義積分の計算

# 例1: 積分区間が無限大 def integrand1(x): return np.exp(-x**2) # ∫₀^∞ exp(-x²) dx = √π/2 result, error = integrate.quad(integrand1, 0, np.inf) exact = np.sqrt(np.pi) / 2 print("広義積分の例:") print(f"∫₀^∞ exp(-x²) dx = {result:.10f} (解析解: {exact:.10f})") print(f"推定誤差: {error:.2e}, 実誤差: {abs(result-exact):.2e}\n") # 例2: 特異点がある場合 def integrand2(x): """x=0 で特異点を持つ""" return 1/np.sqrt(x) # ∫₀¹ 1/√x dx = 2 result, error = integrate.quad(integrand2, 0, 1) exact = 2.0 print(f"∫₀¹ 1/√x dx = {result:.10f} (解析解: {exact:.10f})") print(f"推定誤差: {error:.2e}, 実誤差: {abs(result-exact):.2e}") # 特異点を明示的に指定 result2, error2 = integrate.quad(integrand2, 0, 1, points=[0]) print(f"特異点指定: {result2:.10f}, 誤差: {abs(result2-exact):.2e}")

2.7 練習問題

✏️ 演習1: ∫₀² (x³ - 2x² + x) dx を (1) 解析的に、(2) 台形則 (n=100) で、 (3) Simpson則 (n=100) で求めよ。
✏️ 演習2: 反応速度が r(t) = 0.5 exp(-0.1t) [mol/L/s] で与えられる。 0 ≤ t ≤ 20 秒における総反応量を数値積分で求めよ。

まとめ

← 第1章: 微分の基礎 第3章: 多変数関数の微積分 →