工業数学第4編：数値解析

工学の問題を解析的に解けない場合、あるいは複雑なシステムの数値シミュレーションが必要な場合に**数値解析(Numerical Analysis)**が欠かせません。トランジスタの非線形特性分析、複雑な回路網の過渡応答、有限要素解析など、すべて数値解析に基づいています。本記事では核心的な数値手法を原理からPython実装まで完全解説します。

1. 誤差解析

1.1 誤差の種類

数値計算において「完全な解」は存在しません。誤差の種類を理解することが、信頼性の高い数値計算の出発点です。

丸め誤差(Round-off Error): コンピュータが有限桁数で実数を表現するために生じます。

IEEE 754 倍精度: 約15〜16桁の有効数字

$1/3 = 0.333333\ldots \to 0.3333333333333333 \text{ (有限表現)}$

打ち切り誤差(Truncation Error): 無限級数を有限項で近似するときに生じます。

$e^x = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \cdots \approx 1 + x + \frac{x^2}{2!} \text{ (高次項を打ち切り)}$

絶対誤差(Absolute Error): $E_a = |x_{\text{true}} - x_{\text{approx}}|$

相対誤差(Relative Error): $E_r = \frac{|x_{\text{true}} - x_{\text{approx}}|}{|x_{\text{true}}|}$

百分率誤差: $\epsilon_r = E_r \times 100\%$

有効数字(Significant Figures): $\epsilon_r < 0.5\%$ のとき、有効数字が $d$ 桁あることが保証されます: $d = 2 - \log_{10}(E_r \times 100)$

1.2 数値安定性(Numerical Stability)

誤差が計算過程で増幅されてはいけません。

良条件(Well-conditioned): 入力誤差が出力に大きく増幅されない

悪条件(Ill-conditioned): 入力の小さな変化が出力を大きく変化させる

条件数(Condition Number): $\kappa(A) = \|A\| \cdot \|A^{-1}\|$

線形システム $Ax = b$ において $\kappa(A)$ が大きいと、数値解は不安定になります。

import numpy as np

# ヒルベルト行列（悪名高い悪条件行列）
n = 6
H = np.array([[1/(i+j+1) for j in range(n)] for i in range(n)])
kappa = np.linalg.cond(H)
print(f"ヒルベルト行列 ({n}x{n}) 条件数: {kappa:.2e}")

# 簡単な行列
A = np.array([[2, 1], [1, 3]], dtype=float)
print(f"簡単な行列の条件数: {np.linalg.cond(A):.2f}")

# 誤差伝播の実演
b_true = np.array([3, 4], dtype=float)
x_true = np.linalg.solve(A, b_true)

# 少しの誤差を追加
b_perturbed = b_true + np.array([0.01, 0.0])
x_perturbed = np.linalg.solve(A, b_perturbed)

print(f"\n入力誤差: {np.linalg.norm(b_perturbed - b_true):.4f}")
print(f"出力誤差: {np.linalg.norm(x_perturbed - x_true):.4f}")

1.3 ビッグO記法と収束次数

収束次数(Order of Convergence):

$E_{n+1} \leq C \cdot E_n^p$

$p = 1$ : 1次（線形）収束
$p = 2$ : 2次（2乗）収束（はるかに速い）
$p = 1.618$ : 超線形収束（割線法）

ビッグO記法:

$f(h) = O(h^p) \Leftrightarrow |f(h)| \leq C h^p \text{ as } h \to 0$

例: 前進差分の誤差は $O(h)$ 、中心差分の誤差は $O(h^2)$

2. 非線形方程式の求解

2.1 二分法(Bisection Method)

原理: ボルツァーノの定理 — $f(a) \cdot f(b) < 0$ なら $[a, b]$ の間に根が存在します。

アルゴリズム:

$f(a) \cdot f(b) < 0$ を確認
中点 $c = (a + b)/2$ を計算
$f(c) = 0$ なら終了; $f(a) \cdot f(c) < 0$ なら $b = c$ ; それ以外は $a = c$
$|b - a| < \epsilon$ なら終了

収束: 毎回の反復で区間が半分になるため1次収束。 $n$ 回反復後の誤差: $E_n \leq (b-a)/2^n$

必要反復回数: $n \geq \log_2\!\left(\frac{b-a}{\epsilon}\right)$

def bisection(f, a, b, tol=1e-10, max_iter=100):
    """二分法で f(x) = 0 の根を求める"""
    if f(a) * f(b) > 0:
        raise ValueError("f(a) と f(b) は異符号でなければなりません")

    history = []
    for i in range(max_iter):
        c = (a + b) / 2
        history.append({'iter': i+1, 'a': a, 'b': b, 'c': c,
                        'f(c)': f(c), 'error': abs(b-a)/2})

        if abs(f(c)) < tol or (b - a) / 2 < tol:
            print(f"収束: {i+1} 回反復")
            return c, history

        if f(a) * f(c) < 0:
            b = c
        else:
            a = c

    return (a + b) / 2, history

# 例: x^3 - 2x - 5 = 0
def equation(x):
    return x**3 - 2*x - 5

root, history = bisection(equation, 1, 3)
print(f"根: {root:.10f}")
print(f"検証: f({root:.6f}) = {equation(root):.2e}")

print("\n反復 |   近似値    |   誤差")
for h in history[:8]:
    print(f"  {h['iter']:2d} | {h['c']:.8f} | {h['error']:.2e}")

2.2 ニュートン-ラフソン法(Newton-Raphson Method)

原理: $x_n$ でのテイラー1次近似

$f(x) \approx f(x_n) + f'(x_n)(x - x_n) = 0$

$x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$

幾何的解釈: $x_n$ での接線がx軸と交わる点が $x_{n+1}$

収束: 理想的には2次（2乗）収束 — 正確な桁数が毎回の反復で2倍になります！

$E_{n+1} \approx \frac{f''(x^*)}{2f'(x^*)}E_n^2$

限界と注意点:

$f'(x_n) = 0$ のとき発散（接線がx軸と平行）
収束は保証されない（初期点の選択が重要）
重根では1次収束に低下

def newton_raphson(f, df, x0, tol=1e-10, max_iter=50):
    """ニュートン-ラフソン法"""
    x = x0
    history = [{'iter': 0, 'x': x, 'f(x)': f(x), 'error': float('inf')}]

    for i in range(max_iter):
        fx = f(x)
        dfx = df(x)

        if abs(dfx) < 1e-15:
            print("警告: 導関数が0に近い、発散の危険")
            break

        x_new = x - fx / dfx
        error = abs(x_new - x)
        history.append({'iter': i+1, 'x': x_new,
                        'f(x)': f(x_new), 'error': error})

        if error < tol:
            print(f"収束: {i+1} 回反復")
            return x_new, history

        x = x_new

    return x, history

# f(x) = x^3 - 2x - 5,  f'(x) = 3x^2 - 2
def f(x):
    return x**3 - 2*x - 5

def df(x):
    return 3*x**2 - 2

root, history = newton_raphson(f, df, x0=2.0)
print(f"ニュートン-ラフソン根: {root:.10f}")

print("\n反復 |    近似値      |   誤差")
for h in history:
    err_str = f"{h['error']:.2e}" if h['error'] != float('inf') else "   -"
    print(f"  {h['iter']:2d} | {h['x']:.10f} | {err_str}")

収束速度比較:

import numpy as np
import matplotlib.pyplot as plt

true_root = 2.0945514815423265  # x^3 - 2x - 5 = 0 の真の根

_, hist_bisect = bisection(f, 1, 3, tol=1e-12, max_iter=50)
_, hist_newton = newton_raphson(f, df, x0=2.0, tol=1e-12, max_iter=20)

errors_b = [abs(h['c'] - true_root) for h in hist_bisect]
errors_n = [abs(h['x'] - true_root) for h in hist_newton if h['error'] != float('inf')]

plt.figure(figsize=(10, 5))
plt.semilogy(range(1, len(errors_b)+1), errors_b, 'b-o', label='二分法 (1次収束)')
plt.semilogy(range(1, len(errors_n)+1), errors_n, 'r-s', label='ニュートン-ラフソン (2次収束)')
plt.xlabel('反復回数')
plt.ylabel('絶対誤差 (対数スケール)')
plt.title('二分法 vs ニュートン-ラフソン 収束速度')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('convergence_comparison.png', dpi=150)
plt.show()

2.3 割線法(Secant Method)

$f'(x)$ が未知のとき、有限差分で近似します:

$x_{n+1} = x_n - f(x_n)\cdot\frac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})}$

収束次数: $p = (1+\sqrt{5})/2 \approx 1.618$ （黄金比、超線形）

def secant_method(f, x0, x1, tol=1e-10, max_iter=50):
    """割線法"""
    for i in range(max_iter):
        f0, f1 = f(x0), f(x1)
        if abs(f1 - f0) < 1e-15:
            break
        x2 = x1 - f1 * (x1 - x0) / (f1 - f0)
        if abs(x2 - x1) < tol:
            print(f"収束: {i+1} 回反復")
            return x2
        x0, x1 = x1, x2
    return x1

root_s = secant_method(f, 1.0, 3.0)
print(f"割線法の根: {root_s:.10f}")

2.4 scipy.optimize の活用

from scipy.optimize import fsolve, brentq
import numpy as np

def equation(x):
    return x**3 - 2*x - 5

# ブレント法（二分法 + 割線法 + 逆2次補間のハイブリッド）
root = brentq(equation, 1, 3, xtol=1e-12)
print(f"brentq 根: {root:.12f}")

# 連立非線形方程式
def system(vars):
    x, y = vars
    eq1 = x**2 + y**2 - 4   # 円
    eq2 = x*y - 1             # 双曲線
    return [eq1, eq2]

solutions = []
for x0, y0 in [(1, 1), (-1, 1), (1, -1), (-1, -1)]:
    sol = fsolve(system, [x0, y0], full_output=True)
    if sol[2] == 1:
        solutions.append(sol[0])

print("\n連立方程式の解:")
for sol in solutions:
    print(f"  (x, y) = ({sol[0]:.6f}, {sol[1]:.6f})")

# ダイオード特性方程式（非線形）
def diode_equation(V, V_s=5, R=1000, I_s=1e-12, n=1, V_T=0.02585):
    I = I_s * (np.exp(V / (n * V_T)) - 1)
    return V_s - V - R * I

V_diode = brentq(diode_equation, 0, 1)
I_diode = 1e-12 * (np.exp(V_diode / 0.02585) - 1)
print(f"\nダイオード動作点: V = {V_diode:.4f} V, I = {I_diode*1000:.4f} mA")

3. 数値微分と数値積分

3.1 有限差分(Finite Difference)

テイラー級数から導関数の数値近似を導出します。

前進差分(Forward Difference):

$f'(x) \approx \frac{f(x+h) - f(x)}{h} + O(h)$

1次精度（誤差が $h$ に比例）

後退差分(Backward Difference):

$f'(x) \approx \frac{f(x) - f(x-h)}{h} + O(h)$

中心差分(Central Difference):

$f'(x) \approx \frac{f(x+h) - f(x-h)}{2h} + O(h^2)$

2次精度 — 同じ $h$ で前進差分よりはるかに正確です。

2階導関数:

$f''(x) \approx \frac{f(x+h) - 2f(x) + f(x-h)}{h^2} + O(h^2)$

これは熱方程式・波動方程式の有限差分解法に使われます。

import numpy as np
import matplotlib.pyplot as plt

def numerical_derivative(f, x, h=1e-5, method='central'):
    """数値微分"""
    if method == 'forward':
        return (f(x + h) - f(x)) / h
    elif method == 'backward':
        return (f(x) - f(x - h)) / h
    elif method == 'central':
        return (f(x + h) - f(x - h)) / (2 * h)

def f(x):
    return np.sin(x)

def df_true(x):
    return np.cos(x)

x0 = np.pi / 4
true_deriv = df_true(x0)

h_values = np.logspace(-16, 0, 100)
errors_fwd = []
errors_cen = []

for h in h_values:
    e_fwd = abs(numerical_derivative(f, x0, h, 'forward') - true_deriv)
    e_cen = abs(numerical_derivative(f, x0, h, 'central') - true_deriv)
    errors_fwd.append(e_fwd if e_fwd > 0 else 1e-17)
    errors_cen.append(e_cen if e_cen > 0 else 1e-17)

plt.figure(figsize=(10, 6))
plt.loglog(h_values, errors_fwd, 'b-', linewidth=2, label='前進差分 O(h)')
plt.loglog(h_values, errors_cen, 'r-', linewidth=2, label='中心差分 O(h^2)')
plt.xlabel('h (ステップサイズ)')
plt.ylabel('絶対誤差')
plt.title('数値微分の誤差解析')
plt.legend()
plt.grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.savefig('numerical_derivative_error.png', dpi=150)
plt.show()

3.2 数値積分(Numerical Integration)

台形則(Trapezoidal Rule):

$\int_a^b f(x)\,dx \approx h\left[\frac{f(a)}{2} + f(x_1) + \cdots + f(x_{n-1}) + \frac{f(b)}{2}\right]$

ここで $h = (b-a)/n$ 。誤差: $O(h^2)$

シンプソンの1/3則(Simpson's 1/3 Rule):

各2区間を2次多項式で近似:

$\int_a^b f(x)\,dx \approx \frac{h}{3}\left[f(a) + 4f(x_1) + 2f(x_2) + 4f(x_3) + \cdots + f(b)\right]$

$n$ は偶数である必要があります。誤差: $O(h^4)$ — 台形則よりはるかに正確！

ガウス求積法(Gaussian Quadrature):

最適な積分点と重みを選択して精度を最大化します。

$n$ 点ガウス-ルジャンドル求積法は最高 $(2n-1)$ 次の多項式を正確に積分します。

$\int_{-1}^{1}f(x)\,dx \approx \sum_{i=1}^n w_i f(x_i)$

import numpy as np
from scipy import integrate

def trapezoidal(f, a, b, n):
    """台形則"""
    x = np.linspace(a, b, n+1)
    h = (b - a) / n
    y = f(x)
    return h * (y[0]/2 + np.sum(y[1:-1]) + y[-1]/2)

def simpsons(f, a, b, n):
    """シンプソンの1/3則（n は偶数）"""
    if n % 2 != 0:
        n += 1
    x = np.linspace(a, b, n+1)
    h = (b - a) / n
    y = f(x)
    return h/3 * (y[0] + 4*np.sum(y[1::2]) + 2*np.sum(y[2:-1:2]) + y[-1])

# 例: sin(x) の 0 から pi までの積分 = 2
f_test = np.sin
a, b = 0, np.pi
true_val = 2.0

print("数値積分の精度比較")
print(f"正確な値: {true_val}")
print("-" * 60)
print(f"{'手法':<20} {'n':>6} {'近似値':>15} {'相対誤差':>12}")
print("-" * 60)

for n in [4, 8, 16, 32]:
    trap = trapezoidal(f_test, a, b, n)
    simp = simpsons(f_test, a, b, n)
    print(f"{'台形則':<20} {n:>6} {trap:>15.10f} {abs(trap-true_val)/true_val:>12.2e}")
    print(f"{'シンプソン':<20} {n:>6} {simp:>15.10f} {abs(simp-true_val)/true_val:>12.2e}")

# scipy.integrate.quad — 適応型求積法
result, error = integrate.quad(f_test, a, b)
print(f"{'scipy.quad':<20} {'適応型':>6} {result:>15.10f} {'<'+str(error):>12}")

# 難しい積分: 0 から inf の exp(-x^2) dx = sqrt(pi)/2
def gaussian(x):
    return np.exp(-x**2)

result_gauss, _ = integrate.quad(gaussian, 0, np.inf)
print(f"\nガウス積分: {result_gauss:.10f}")
print(f"正確な値 (sqrt(pi)/2): {np.sqrt(np.pi)/2:.10f}")

4. ODE 数値解法

4.1 オイラー法(Euler's Method)

アイデア: 現在の傾きで次の点を予測します。

$y_{n+1} = y_n + hf(x_n, y_n)$

局所打ち切り誤差: $O(h^2)$ 、全体誤差: $O(h)$ （1次法）

import numpy as np

def euler_method(f, x0, y0, h, x_end):
    """オイラー法"""
    x_vals = [x0]
    y_vals = [y0]
    x = x0
    y = y0

    while x < x_end - 1e-10:
        h_actual = min(h, x_end - x)
        y = y + h_actual * f(x, y)
        x = x + h_actual
        x_vals.append(x)
        y_vals.append(y)

    return np.array(x_vals), np.array(y_vals)

4.2 改良オイラー法 / ホイン法(Heun's Method)

アイデア: 予測(predict)してから修正(correct)します。

$k_1 = f(x_n, y_n)$ $\tilde{y}_{n+1} = y_n + hk_1 \quad \text{(予測)}$ $k_2 = f(x_n + h, \tilde{y}_{n+1})$ $y_{n+1} = y_n + \frac{h}{2}(k_1 + k_2) \quad \text{(修正)}$

全体誤差: $O(h^2)$ — オイラー法よりはるかに正確です。

4.3 4次ルンゲ-クッタ法(RK4)

工学で最も広く使用される ODE 数値解法です。

$k_1 = hf(x_n, y_n)$ $k_2 = hf\!\left(x_n + \frac{h}{2},\, y_n + \frac{k_1}{2}\right)$ $k_3 = hf\!\left(x_n + \frac{h}{2},\, y_n + \frac{k_2}{2}\right)$ $k_4 = hf(x_n + h,\, y_n + k_3)$ $y_{n+1} = y_n + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4)$

直感: 区間内の4つの傾きの加重平均を取ります。

$k_1$ : 始点での傾き
$k_2$ , $k_3$ : 中点での傾き
$k_4$ : 終点での傾き

全体誤差: $O(h^4)$ — 誤差が $h^4$ に比例します

比較: $h$ を半分にすると

オイラー法: 誤差 1/2 倍
RK4: 誤差 1/16 倍

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

def rk4_method(f, x0, y0, h, x_end):
    """4次ルンゲ-クッタ法"""
    x_vals = [x0]
    y_vals = [y0]
    x = x0
    y = y0

    while x < x_end - 1e-10:
        h_actual = min(h, x_end - x)
        k1 = h_actual * f(x, y)
        k2 = h_actual * f(x + h_actual/2, y + k1/2)
        k3 = h_actual * f(x + h_actual/2, y + k2/2)
        k4 = h_actual * f(x + h_actual, y + k3)
        y = y + (k1 + 2*k2 + 2*k3 + k4) / 6
        x = x + h_actual
        x_vals.append(x)
        y_vals.append(y)

    return np.array(x_vals), np.array(y_vals)

# RLC 回路シミュレーション
R, L, C = 50, 0.2, 100e-6
Vs = 100

def rlc_ode(t, state):
    q, dq = state
    d2q = (Vs - R * dq - q/C) / L
    return [dq, d2q]

y0 = [0.0, 0.0]
t_end = 0.02

# scipy solve_ivp（基準解）
t_ref = np.linspace(0, t_end, 5000)
sol_ref = solve_ivp(rlc_ode, [0, t_end], y0,
                    t_eval=t_ref, method='RK45',
                    rtol=1e-10, atol=1e-12)
v_ref = sol_ref.y[0] / C

def rlc_2nd_order(t, y):
    q, dq = y[0], y[1]
    return [dq, (Vs - R*dq - q/C) / L]

fig, axes = plt.subplots(2, 1, figsize=(12, 10))
axes[0].plot(sol_ref.t * 1000, v_ref, 'k-', linewidth=2.5, label='基準 (RK45 高精度)')

for h in [0.001, 0.0005, 0.0001]:
    n_steps = int(t_end / h)
    t_rk4 = np.linspace(0, t_end, n_steps + 1)
    y = np.array(y0, dtype=float)
    y_hist = [y.copy()]
    for i in range(n_steps):
        t_i = t_rk4[i]
        k1 = h * np.array(rlc_2nd_order(t_i, y))
        k2 = h * np.array(rlc_2nd_order(t_i + h/2, y + k1/2))
        k3 = h * np.array(rlc_2nd_order(t_i + h/2, y + k2/2))
        k4 = h * np.array(rlc_2nd_order(t_i + h, y + k3))
        y = y + (k1 + 2*k2 + 2*k3 + k4) / 6
        y_hist.append(y.copy())
    y_arr = np.array(y_hist)
    v_rk4 = y_arr[:, 0] / C
    axes[0].plot(t_rk4 * 1000, v_rk4, '--',
                 label=f'RK4 (h={h*1000:.1f}ms)', alpha=0.8)

axes[0].set_xlabel('時間 (ms)')
axes[0].set_ylabel('キャパシタ電圧 (V)')
axes[0].set_title('RLC 直列回路 - RK4 数値解法')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

i_ref = sol_ref.y[1]
axes[1].plot(sol_ref.t * 1000, i_ref * 1000, 'k-', linewidth=2, label='電流 i(t)')
axes[1].set_xlabel('時間 (ms)')
axes[1].set_ylabel('電流 (mA)')
axes[1].set_title('RLC 回路電流')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('rlc_numerical.png', dpi=150)
plt.show()

4.4 硬い ODE と適応ステップサイズ

硬い ODE(Stiff ODE): 時定数が非常に異なるモードが共存するシステムです。

例: 高速な RC 過渡応答と低速な温度変化が同時に存在するシステム

from scipy.integrate import solve_ivp
import numpy as np

# 硬い方程式の例: Robertson 化学反応
def robertson(t, y):
    dy0 = -0.04*y[0] + 1e4*y[1]*y[2]
    dy1 = 0.04*y[0] - 1e4*y[1]*y[2] - 3e7*y[1]**2
    dy2 = 3e7*y[1]**2
    return [dy0, dy1, dy2]

y0 = [1, 0, 0]
t_span = (0, 1e11)
t_eval = np.logspace(-6, 11, 1000)

# 非硬い手法 (RK45): 遅すぎる
# 硬い手法 (Radau, BDF): 効率的
sol = solve_ivp(robertson, t_span, y0,
                method='Radau',
                t_eval=t_eval,
                rtol=1e-6, atol=1e-9)

print(f"関数評価回数: {sol.nfev}")
print(f"ヤコビアン評価回数: {sol.njev}")

5. 線形システムの数値解法

5.1 ガウス消去法(Gaussian Elimination)

前進消去: 行演算で上三角行列を作成

後退代入: 上三角システムを解く

演算回数: $O(n^3)$

部分ピボット選択(Partial Pivoting): 各列で絶対値最大の要素をピボットとして使用し、数値安定性を向上させます。

import numpy as np

def gaussian_elimination_pivoting(A, b):
    """部分ピボット選択付きガウス消去法"""
    n = len(b)
    A = np.array(A, dtype=float)
    b = np.array(b, dtype=float)

    for i in range(n):
        max_row = i + np.argmax(abs(A[i:, i]))
        A[[i, max_row]] = A[[max_row, i]]
        b[[i, max_row]] = b[[max_row, i]]

        if abs(A[i, i]) < 1e-15:
            raise ValueError("特異行列！")

        for j in range(i+1, n):
            factor = A[j, i] / A[i, i]
            A[j, i:] -= factor * A[i, i:]
            b[j] -= factor * b[i]

    x = np.zeros(n)
    for i in range(n-1, -1, -1):
        x[i] = (b[i] - np.dot(A[i, i+1:], x[i+1:])) / A[i, i]

    return x

# 例: 3ループ回路網 (KVL)
A = np.array([
    [10, -3, -1],
    [-3, 12, -2],
    [-1, -2,  8]
], dtype=float)
b = np.array([100, 50, -30], dtype=float)

x_gauss = gaussian_elimination_pivoting(A, b)
x_numpy = np.linalg.solve(A, b)

print("回路ループ電流（ガウス消去法）:")
for i, v in enumerate(x_gauss):
    print(f"  I{i+1} = {v:.6f} A")
print(f"\nnumpy との最大誤差: {np.max(abs(x_gauss - x_numpy)):.2e}")

5.2 LU 分解(LU Decomposition)

$A = LU$

$L$ : 下三角行列、 $U$ : 上三角行列

利点: 同じ $A$ に対して複数の $b$ を効率的に解くことができます。

$Ax = b \to Ly = b$ （前進代入）、 $Ux = y$ （後退代入）

from scipy.linalg import lu, lu_factor, lu_solve
import numpy as np

A = np.array([[2, 1, -1], [-3, -1, 2], [-2, 1, 2]], dtype=float)
b = np.array([8, -11, -3], dtype=float)

P, L, U = lu(A)
print("L 行列:")
print(np.round(L, 4))
print("U 行列:")
print(np.round(U, 4))

lu_and_piv = lu_factor(A)
x = lu_solve(lu_and_piv, b)
print(f"\n解: {x}")

# 複数の右辺ベクトルに対して効率的に解く
for b_new in [b, b+1, b*2]:
    x_new = lu_solve(lu_and_piv, b_new)  # L, U を再利用

5.3 反復法(Iterative Methods)

大規模な疎行列では直接法より効率的です。

ヤコビ法(Jacobi Method): $x_i^{(k+1)} = \frac{1}{a_{ii}}\left(b_i - \sum_{j \neq i}a_{ij}x_j^{(k)}\right)$

すべての変数を同時に更新します。

ガウス-ザイデル法(Gauss-Seidel Method): $x_i^{(k+1)} = \frac{1}{a_{ii}}\left(b_i - \sum_{j < i}a_{ij}x_j^{(k+1)} - \sum_{j > i}a_{ij}x_j^{(k)}\right)$

更新した値をすぐに使用します → ヤコビ法より速い収束。

収束条件: 対角優勢行列

$|a_{ii}| > \sum_{j \neq i}|a_{ij}| \quad \forall i$

import numpy as np

def gauss_seidel(A, b, x0=None, tol=1e-10, max_iter=1000):
    """ガウス-ザイデル反復法"""
    n = len(b)
    x = np.zeros(n) if x0 is None else x0.copy()

    for iteration in range(max_iter):
        x_old = x.copy()
        for i in range(n):
            sum1 = np.dot(A[i, :i], x[:i])
            sum2 = np.dot(A[i, i+1:], x_old[i+1:])
            x[i] = (b[i] - sum1 - sum2) / A[i, i]

        residual = np.linalg.norm(x - x_old)
        if residual < tol:
            print(f"収束: {iteration+1} 回反復、残差: {residual:.2e}")
            return x

    print(f"最大反復 {max_iter} で終了、残差: {residual:.2e}")
    return x

# 対角優勢行列の例（ノード電圧解析）
A = np.array([
    [ 4, -1, -1,  0],
    [-1,  4,  0, -1],
    [-1,  0,  4, -1],
    [ 0, -1, -1,  4]
], dtype=float)
b = np.array([5, 10, 15, 20], dtype=float)

x_gs = gauss_seidel(A, b)
x_exact = np.linalg.solve(A, b)
print(f"ガウス-ザイデル: {x_gs}")
print(f"正確な解: {x_exact}")
print(f"最大誤差: {np.max(abs(x_gs - x_exact)):.2e}")

6. 補間と近似

6.1 ラグランジュ補間(Lagrange Interpolation)

$n+1$ 個の点 $(x_0, y_0), \ldots, (x_n, y_n)$ を通る $n$ 次多項式:

$P_n(x) = \sum_{k=0}^n y_k L_k(x)$

基底多項式（ラグランジュ基底）: $L_k(x) = \prod_{j=0, j\neq k}^n \frac{x - x_j}{x_k - x_j}$

ルンゲ現象(Runge's Phenomenon): 等間隔点での高次補間時に両端で振動が発生します。

解決策: チェビシェフ点(Chebyshev nodes)を使用します

$x_k = \cos\!\left(\frac{(2k+1)\pi}{2(n+1)}\right), \quad k = 0, 1, \ldots, n$

6.2 スプライン補間(Spline Interpolation)

各区間で低次多項式を使用します。3次スプラインが最も一般的です。

隣接区間で関数値・1階・2階微分が連続 → 滑らかな曲線。

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import lagrange, CubicSpline, interp1d

# 測定データ（温度 vs 抵抗）
x_data = np.array([0, 20, 40, 60, 80, 100], dtype=float)   # 摂氏
y_data = np.array([100, 108, 116, 125, 134, 143], dtype=float)  # オーム

x_fine = np.linspace(0, 100, 300)

poly = lagrange(x_data, y_data)
y_lagrange = poly(x_fine)

cs = CubicSpline(x_data, y_data)
y_spline = cs(x_fine)

f_linear = interp1d(x_data, y_data)
y_linear = f_linear(x_fine)

plt.figure(figsize=(10, 6))
plt.scatter(x_data, y_data, s=100, zorder=5, label='測定データ', color='black')
plt.plot(x_fine, y_linear, 'g--', linewidth=1.5, label='線形補間')
plt.plot(x_fine, y_lagrange, 'r-', linewidth=2, label='ラグランジュ補間')
plt.plot(x_fine, y_spline, 'b-', linewidth=2, label='3次スプライン')
plt.xlabel('温度 (C)')
plt.ylabel('抵抗 (オーム)')
plt.title('温度-抵抗 補間比較')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('interpolation.png', dpi=150)
plt.show()

print(f"50度での抵抗予測:")
print(f"  線形: {f_linear(50):.2f} オーム")
print(f"  スプライン: {cs(50):.2f} オーム")

6.3 最小二乗法(Least Squares Method)

$m$ 個のデータ点を $n$ 次多項式（ $m > n+1$ ）で最もよく近似します:

$\min \sum_{i=1}^m [y_i - p(x_i)]^2$

正規方程式(Normal Equations): $A^T A \mathbf{c} = A^T \mathbf{y}$

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

np.random.seed(42)
t_data = np.linspace(0, 5, 30)
y_true = 3 * np.exp(-0.7 * t_data)
y_noisy = y_true + 0.2 * np.random.randn(len(t_data))

t_fine = np.linspace(0, 5, 200)

plt.figure(figsize=(12, 8))
plt.scatter(t_data, y_noisy, s=50, color='black', zorder=5, label='測定データ')
plt.plot(t_fine, 3*np.exp(-0.7*t_fine), 'k-', linewidth=2, label='真の信号')

for deg in [1, 2, 3, 5]:
    coeffs = np.polyfit(t_data, y_noisy, deg)
    y_fit = np.polyval(coeffs, t_fine)
    plt.plot(t_fine, y_fit, '--', linewidth=1.5, label=f'{deg}次多項式')

plt.ylim([-1, 4])
plt.xlabel('時間 (s)')
plt.ylabel('振幅')
plt.title('最小二乗多項式フィッティング')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('polynomial_fit.png', dpi=150)
plt.show()

# 非線形最小二乗: 指数関数フィッティング
def exp_model(t, A, k):
    return A * np.exp(-k * t)

popt, pcov = curve_fit(exp_model, t_data, y_noisy, p0=[3, 0.5])
A_fit, k_fit = popt
A_std, k_std = np.sqrt(np.diag(pcov))

print(f"指数フィッティング結果:")
print(f"  A = {A_fit:.4f} +/- {A_std:.4f}  (真値: 3)")
print(f"  k = {k_fit:.4f} +/- {k_std:.4f}  (真値: 0.7)")

y_pred = exp_model(t_data, *popt)
ss_res = np.sum((y_noisy - y_pred)**2)
ss_tot = np.sum((y_noisy - np.mean(y_noisy))**2)
r_squared = 1 - ss_res / ss_tot
print(f"  R^2 = {r_squared:.6f}")

7. 固有値問題

7.1 べき乗法(Power Method)

行列 $A$ の最大絶対固有値 $\lambda_1$ と対応する固有ベクトル $\mathbf{v}_1$ を求めます。

アルゴリズム:

初期ベクトル $\mathbf{q}_0$ を選択
$\mathbf{z}_{k+1} = A\mathbf{q}_k$
$\mathbf{q}_{k+1} = \mathbf{z}_{k+1}/\|\mathbf{z}_{k+1}\|$
$\lambda_{k+1} = \mathbf{q}_{k+1}^T A \mathbf{q}_{k+1}$ （レイリー商）
収束するまで繰り返す

収束速度: $|\lambda_2/\lambda_1|^k$

import numpy as np

def power_method(A, tol=1e-10, max_iter=1000):
    """べき乗法（最大固有値）"""
    n = A.shape[0]
    q = np.ones(n) / np.sqrt(n)

    for k in range(max_iter):
        z = A @ q
        eigenvalue = np.dot(q, z)
        q_new = z / np.linalg.norm(z)

        if np.linalg.norm(q_new - q) < tol or np.linalg.norm(q_new + q) < tol:
            print(f"収束: {k+1} 回反復")
            return eigenvalue, q_new

        q = q_new

    return eigenvalue, q

# 振動系の固有周波数解析
k_spring = 1000  # N/m
m_mass = 1       # kg

K = np.array([[2*k_spring, -k_spring], [-k_spring, k_spring]], dtype=float)
M = np.array([[m_mass, 0], [0, m_mass]], dtype=float)

M_inv_sqrt = np.diag(1/np.sqrt(np.diag(M)))
A_sym = M_inv_sqrt @ K @ M_inv_sqrt

lam_max, v_max = power_method(A_sym)
omega_max = np.sqrt(lam_max)

eigenvalues_all = np.linalg.eigvalsh(A_sym)
omega_all = np.sqrt(eigenvalues_all)

print("振動系の固有角周波数:")
for i, omega in enumerate(omega_all):
    print(f"  モード {i+1}: {omega/(2*np.pi):.3f} Hz")

print(f"\nべき乗法の最大固有周波数: {omega_max/(2*np.pi):.3f} Hz")

7.2 QR アルゴリズム

すべての固有値を同時に求めます。現代数値線形代数の核心アルゴリズムです。

アイデア: 繰り返しの QR 分解により行列が上三角形に収束すると、対角要素が固有値になります。

import numpy as np

def qr_algorithm(A, max_iter=1000, tol=1e-10):
    """基本 QR アルゴリズム"""
    A_k = A.copy().astype(float)

    for k in range(max_iter):
        Q, R = np.linalg.qr(A_k)
        A_k_new = R @ Q

        off_diag_norm = np.max(np.abs(np.tril(A_k_new, -1)))
        if off_diag_norm < tol:
            print(f"QR 収束: {k+1} 回反復")
            break

        A_k = A_k_new

    return np.diag(A_k)

A = np.array([
    [4, 1, 0],
    [1, 3, 1],
    [0, 1, 2]
], dtype=float)

eig_qr = qr_algorithm(A)
eig_numpy = np.linalg.eigvalsh(A)

print("QR アルゴリズム固有値:", np.sort(eig_qr))
print("numpy 固有値:         ", np.sort(eig_numpy))

7.3 構造共振解析

import numpy as np
import matplotlib.pyplot as plt
from scipy import linalg

def vibration_analysis(K, M, n_modes=None):
    """
    一般化固有値問題: K*v = omega^2 * M * v
    K: 剛性行列、M: 質量行列
    """
    eigenvalues, eigenvectors = linalg.eigh(K, M)
    omega = np.sqrt(np.abs(eigenvalues))
    freq = omega / (2 * np.pi)

    if n_modes is None:
        n_modes = len(eigenvalues)

    print("固有周波数解析結果:")
    print("-" * 40)
    for i in range(n_modes):
        print(f"モード {i+1}: {freq[i]:.3f} Hz ({omega[i]:.2f} rad/s)")

    return omega[:n_modes], eigenvectors[:, :n_modes]

# 3自由度振動系
k1, k2, k3 = 1000, 1500, 800   # N/m
m1, m2, m3 = 1.0, 1.5, 0.8    # kg

K = np.array([
    [k1+k2,   -k2,    0],
    [-k2,   k2+k3,  -k3],
    [0,       -k3,   k3]
], dtype=float)

M = np.diag([m1, m2, m3])

omegas, modes = vibration_analysis(K, M)

positions = np.arange(1, 4)
fig, axes = plt.subplots(1, 3, figsize=(14, 5))

for i, ax in enumerate(axes):
    mode_shape = modes[:, i] / np.max(abs(modes[:, i]))
    ax.plot([0] + list(positions), [0] + list(mode_shape),
            'bo-', linewidth=2, markersize=10)
    ax.axhline(y=0, color='gray', linestyle='--', alpha=0.5)
    ax.set_title(f'モード {i+1}\n{omegas[i]/(2*np.pi):.2f} Hz')
    ax.set_xlabel('自由度')
    ax.set_ylabel('正規化モード形状')
    ax.grid(True, alpha=0.3)
    ax.set_xticks([0, 1, 2, 3])
    ax.set_xticklabels(['固定端', '質量1', '質量2', '質量3'])

plt.suptitle('3自由度振動系のモード形状')
plt.tight_layout()
plt.savefig('vibration_modes.png', dpi=150)
plt.show()

8. AI/ML における数値解析：自動微分

8.1 数値微分 vs 自動微分

数値微分は有限差分で導関数を近似します。一方、自動微分（autograd）は計算グラフを通じて連鎖律を機械的に適用し、正確な導関数を計算します。

import numpy as np

def numerical_grad(f, x, h=1e-5):
    """中心差分によるグラジエント（ベクトル入力）"""
    grad = np.zeros_like(x, dtype=float)
    for i in range(len(x)):
        x_plus  = x.copy().astype(float); x_plus[i]  += h
        x_minus = x.copy().astype(float); x_minus[i] -= h
        grad[i] = (f(x_plus) - f(x_minus)) / (2 * h)
    return grad

# f(x, y) = x^2 + 3*x*y + y^2,  真のグラジエント: [2x+3y, 3x+2y]
def f_quad(x):
    return x[0]**2 + 3*x[0]*x[1] + x[1]**2

x0 = np.array([1.0, 2.0])
grad_num  = numerical_grad(f_quad, x0)
grad_true = np.array([2*x0[0] + 3*x0[1], 3*x0[0] + 2*x0[1]])

print("(1, 2) でのグラジエント比較:")
print(f"  数値微分: {grad_num}")
print(f"  解析値:   {grad_true}")
print(f"  最大誤差: {np.max(abs(grad_num - grad_true)):.2e}")

主な違い:

特性	数値微分	自動微分 (AD)
精度	$O(h^2)$ 、浮動小数点で制限	機械精度
コスト（グラジエント）	$n$ 回の順伝播 ( $n$ パラメータ)	逆伝播1回 (逆モード)
ステップサイズ調整	必要	不要
ML での用途	有限差分による勾配チェック	PyTorch / JAX での誤差逆伝播

8.2 勾配降下法のバリエーション

数値最適化は ML モデルの学習の核心です。各オプティマイザの背後にある数学を理解すると、適切なものを選べます。

import numpy as np
import matplotlib.pyplot as plt

# Rosenbrock 関数: f(x,y) = (1-x)^2 + 100*(y-x^2)^2
def rosenbrock(w):
    x, y = w[0], w[1]
    return (1 - x)**2 + 100*(y - x**2)**2

def rosenbrock_grad(w):
    x, y = w[0], w[1]
    gx = -2*(1 - x) - 400*x*(y - x**2)
    gy = 200*(y - x**2)
    return np.array([gx, gy])

w0 = np.array([-1.5, 1.5])
lr = 0.001
n_iter = 3000

# 標準 SGD
w_sgd = w0.copy()
history_sgd = [w_sgd.copy()]
for _ in range(n_iter):
    w_sgd -= lr * rosenbrock_grad(w_sgd)
    history_sgd.append(w_sgd.copy())

# モーメンタム SGD
w_mom = w0.copy()
v = np.zeros(2)
beta = 0.9
history_mom = [w_mom.copy()]
for _ in range(n_iter):
    v = beta * v - lr * rosenbrock_grad(w_mom)
    w_mom += v
    history_mom.append(w_mom.copy())

# Adam
w_adam = w0.copy()
m_adam, v_adam = np.zeros(2), np.zeros(2)
beta1, beta2, eps = 0.9, 0.999, 1e-8
history_adam = [w_adam.copy()]
for t in range(1, n_iter + 1):
    g = rosenbrock_grad(w_adam)
    m_adam = beta1 * m_adam + (1 - beta1) * g
    v_adam = beta2 * v_adam + (1 - beta2) * g**2
    m_hat = m_adam / (1 - beta1**t)
    v_hat = v_adam / (1 - beta2**t)
    w_adam -= lr * m_hat / (np.sqrt(v_hat) + eps)
    history_adam.append(w_adam.copy())

loss_sgd  = [rosenbrock(w) for w in history_sgd]
loss_mom  = [rosenbrock(w) for w in history_mom]
loss_adam = [rosenbrock(w) for w in history_adam]

plt.figure(figsize=(10, 5))
plt.semilogy(loss_sgd,  label='SGD')
plt.semilogy(loss_mom,  label='モーメンタム')
plt.semilogy(loss_adam, label='Adam')
plt.xlabel('反復回数')
plt.ylabel('損失 (対数スケール)')
plt.title('Rosenbrock 関数でのオプティマイザ収束比較')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('optimizer_comparison.png', dpi=150)
plt.show()

print(f"最終損失 - SGD: {loss_sgd[-1]:.4f}, モーメンタム: {loss_mom[-1]:.4f}, Adam: {loss_adam[-1]:.6f}")

まとめ

本ガイドで扱った内容:

誤差解析: 丸め/打ち切り誤差、絶対/相対誤差、数値安定性、条件数
非線形方程式: 二分法（1次）、ニュートン-ラフソン（2次）、割線法（1.618次）、scipy.optimize
数値微分: 前進/後退/中心差分、誤差解析、最適ステップサイズ
数値積分: 台形則、シンプソン則、ガウス求積法、scipy.integrate
ODE 数値解法: オイラー法、ホイン法、RK4 完全実装、硬い ODE
線形システム: ガウス消去法、LU 分解、反復法、numpy.linalg
補間: ラグランジュ、スプライン、ルンゲ現象、最小二乗、非線形カーブフィッティング
固有値問題: べき乗法、QR アルゴリズム、振動解析
AI/ML への応用: 数値微分 vs 自動微分、勾配降下法のバリエーション

クイズ

Q1. ニュートン-ラフソン法が収束しないのはどのような場合ですか？

答え: ニュートン-ラフソン法はいくつかの状況で失敗します:

導関数がゼロ: ある反復で $f'(x_n) = 0$ になると、更新式が未定義（ゼロ除算）となります。
初期点の選択が悪い: $x_0$ が根から遠いと、反復が振動したり発散したりすることがあります。
重根: $f'(x^*) = 0$ となる重根では、収束が2次から1次に低下します。
平坦な領域 / 極値の近く: 局所的な極小や極大の近くでは、接線がx軸から遠い位置で交差します。

解説: 2次収束の証明には $f'(x^*) \neq 0$ と根に十分近い初期点が必要です。実用上は二分法でまず根を挟んでから、ニュートン-ラフソン法に切り替えて高速な最終収束を得る方法が安全です。

Q2. 複数の右辺ベクトルを解く場合、なぜ LU 分解がガウス消去法より優れているのですか？

答え: LU 分解は高コストの因数分解ステップと安価な解法ステップを分離します。

解説:

ガウス消去法は1システムあたり $O(n^3)$ かかります。同じ行列 $A$ に対して $k$ 個の異なる右辺ベクトル $b_1, \ldots, b_k$ を解く場合、総コストは $O(kn^3)$ です。
LU 分解は $A = LU$ の因数分解を一度だけ $O(n^3)$ で行い、各新しい右辺には2つの三角解法（ $Ly = b$ 、 $Ux = y$ ）のみ必要で、それぞれ $O(n^2)$ です。
総コスト: $O(n^3) + k \cdot O(n^2)$ — $k$ が大きいとき（例：複数荷重ケースの有限要素解析）に劇的に安くなります。

Q3. ノイズのあるデータに対して、ラグランジュ補間とスプライン補間の違いは何ですか？

答え: ノイズのあるデータにはスプライン補間が圧倒的に優れています。

解説:

ラグランジュ補間は $n+1$ 個のデータ点すべてを通る $n$ 次多項式を1つ求めます。ノイズのあるデータでは、多項式がノイズで汚染された測定値を正確に通るため、ルンゲ現象の振動が生じ、高次になるほど悪化します。
3次スプラインは各部分区間で次数3の多項式を使用し、ノットでの滑らかさ条件（1階・2階微分の連続性）を強制します。これにより振動を避け、個々のデータ点にノイズがあっても物理的に合理的な滑らかな曲線を生成します。
ノイズを補間したくない場合は最小二乗フィッティング（補間ではなく近似）が最善です: 各点を正確に通るのではなく残差を最小化する滑らかなモデルを当てはめます。

Q4. なぜルンゲ-クッタ4次法は4次精度を持つのですか？

答え: RK4 は4つの関数評価を行い、テイラー展開の $h^4$ 項まで一致する重み付き組み合わせを使用します。

解説: オイラー法はテイラー級数を1次項で打ち切り、大域誤差は $O(h)$ です。RK4 は区間の始点・2つの中点・終点の4つの傾きを計算し:

$y_{n+1} = y_n + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4)$

この組み合わせは $y(x_{n+1})$ のテイラー展開と $h^4$ 項まで一致します。局所誤差は $O(h^5)$ 、大域誤差は $O(h^4)$ です。ステップサイズを半分にすると大域誤差は $2^4 = 16$ 分の1になります。

Q5. 自動微分（autograd）と数値微分はどう違うのですか？

答え: 自動微分は計算グラフを通じた連鎖律の機械的適用により（浮動小数点精度まで）正確な導関数を計算します。一方、数値微分は有限差分で近似するため、打ち切り誤差と丸め誤差のトレードオフが生じます。

解説:

数値微分: $f'(x) \approx [f(x+h) - f(x)] / h$ 。最適な $h \approx \sqrt{\epsilon_{\text{machine}}} \approx 10^{-8}$ で約8桁の正確さが得られますが、精度は本質的に限られています。
逆モード自動微分（誤差逆伝播）: 関数を素子演算に分解し、各ステップで連鎖律を正確に適用します。全パラメータに対するグラジエントを、順伝播とほぼ同コストの1回の逆伝播で計算できます — 数百万のパラメータを持つニューラルネットワークに不可欠です。
実用上の主な違い: AD はステップサイズ調整不要で機械精度のグラジエントを提供するため、PyTorch、JAX、TensorFlow の基盤となっています。

参考資料

Burden, R. & Faires, J. "Numerical Analysis", 10th Edition, Cengage
Heath, M. "Scientific Computing: An Introductory Survey", 2nd Edition
Press, W. et al. "Numerical Recipes: The Art of Scientific Computing", 3rd Edition
Kreyszig, E. "Advanced Engineering Mathematics", 10th Edition, Wiley
NumPy 線形代数ドキュメント
SciPy 最適化ドキュメント
SciPy 積分ドキュメント

工業数学 第4編：数値解析