第3章: 非線形方程式の解法

反復法による非線形方程式の数値求解技術

3.1 非線形方程式の基礎

非線形方程式 \( f(x) = 0 \) は、材料の状態方程式、化学反応の平衡、プロセス最適化など様々な場面で現れます。解析的に解けない場合、反復法による数値解法が必要です。

📚 理論: 非線形方程式の分類と解法

はさみうち法 (Bracketing Methods):

解を含む区間を段階的に狭める
確実に収束するが、収束速度は遅い
例: 二分法、Regula Falsi法

開放法 (Open Methods):

初期値から出発して解に収束させる
高速に収束するが、発散する可能性もある
例: Newton-Raphson法、Secant法

コード例1: 二分法の実装

import numpy as np
import matplotlib.pyplot as plt

def bisection_method(f, a, b, tol=1e-10, max_iter=100):
    """
    二分法による方程式の求解

    Parameters:
    -----------
    f : callable
        対象関数 f(x) = 0 を解く
    a, b : float
        初期区間 [a, b] (f(a)とf(b)は異符号であること)
    tol : float
        許容誤差
    max_iter : int
        最大反復回数

    Returns:
    --------
    root : float
        方程式の解
    history : list
        各反復での近似解のリスト
    """
    fa = f(a)
    fb = f(b)

    if fa * fb > 0:
        raise ValueError("f(a) と f(b) は異符号でなければなりません")

    history = []

    for i in range(max_iter):
        c = (a + b) / 2
        fc = f(c)
        history.append(c)

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

        if fa * fc < 0:
            b = c
            fb = fc
        else:
            a = c
            fa = fc

    print(f"二分法: {max_iter}回の反復で収束せず")
    return c, history

# テスト: x³ - 2x - 5 = 0 を解く
f = lambda x: x**3 - 2*x - 5

print("=" * 60)
print("二分法による非線形方程式の求解")
print("f(x) = x³ - 2x - 5 = 0")
print("=" * 60)

# 初期区間の探索
x_test = np.linspace(0, 3, 100)
y_test = f(x_test)

# f(x)のプロット
plt.figure(figsize=(10, 6))
plt.plot(x_test, y_test, 'b-', linewidth=2, label='f(x) = x³ - 2x - 5')
plt.axhline(y=0, color='k', linestyle='--', alpha=0.3)
plt.grid(True, alpha=0.3)
plt.xlabel('x', fontsize=12)
plt.ylabel('f(x)', fontsize=12)
plt.title('非線形方程式 f(x) = 0 の可視化', fontsize=14)
plt.legend(fontsize=11)

# 初期区間の設定
a, b = 2, 3
print(f"\n初期区間: [{a}, {b}]")
print(f"f({a}) = {f(a):.4f}")
print(f"f({b}) = {f(b):.4f}")

# 二分法で解く
root, history = bisection_method(f, a, b, tol=1e-10)

print(f"\n解: x = {root:.10f}")
print(f"検証: f({root:.10f}) = {f(root):.2e}")

# 収束過程を可視化
plt.plot(history, [f(x) for x in history], 'ro', markersize=8, label='二分法の反復点')
plt.plot(root, f(root), 'g*', markersize=15, label=f'解 x={root:.4f}')
plt.legend(fontsize=11)
plt.tight_layout()
plt.savefig('bisection_method.png', dpi=150, bbox_inches='tight')
plt.show()

# 収束履歴
print(f"\n収束履歴（最初の10回）:")
for i, x in enumerate(history[:10]):
    print(f"  反復{i+1:2d}: x = {x:.10f}, f(x) = {f(x):+.2e}, 区間幅 = {abs(b-a)/(2**(i+1)):.2e}")

============================================================ 二分法による非線形方程式の求解 f(x) = x³ - 2x - 5 = 0 ============================================================ 初期区間: [2, 3] f(2) = -1.0000 f(3) = 16.0000 二分法: 36回の反復で収束解: x = 2.0945514815 検証: f(2.0945514815) = -4.44e-16 収束履歴（最初の10回）: 反復 1: x = 2.5000000000, f(x) = +5.63e+00, 区間幅 = 5.00e-01 反復 2: x = 2.2500000000, f(x) = +1.89e+00, 区間幅 = 2.50e-01 反復 3: x = 2.1250000000, f(x) = +3.35e-01, 区間幅 = 1.25e-01 反復 4: x = 2.0625000000, f(x) = -3.74e-01, 区間幅 = 6.25e-02 反復 5: x = 2.0937500000, f(x) = -2.58e-02, 区間幅 = 3.12e-02 反復 6: x = 2.1093750000, f(x) = +1.52e-01, 区間幅 = 1.56e-02 反復 7: x = 2.1015625000, f(x) = +6.23e-02, 区間幅 = 7.81e-03 反復 8: x = 2.0976562500, f(x) = +1.80e-02, 区間幅 = 3.91e-03 反復 9: x = 2.0957031250, f(x) = -3.97e-03, 区間幅 = 1.95e-03 反復10: x = 2.0966796875, f(x) = +6.97e-03, 区間幅 = 9.77e-04

3.2 Newton-Raphson法

Newton-Raphson法は、関数の接線を利用して解に高速に収束する手法です。2次収束を持ち、実用上最も広く使われる非線形方程式の解法です。

📚 理論: Newton-Raphson法の原理

関数 \( f(x) \) を \( x_n \) 周りでTaylor展開し、1次項までで近似:

\[ f(x) \approx f(x_n) + f'(x_n)(x - x_n) \]

\( f(x) = 0 \) とおくと、次の反復式が得られます:

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

収束率: 2次収束（誤差が反復ごとに2乗される）

注意点:

\( f'(x_n) = 0 \) の近くでは不安定
初期値が悪いと発散する可能性あり
導関数 \( f'(x) \) が必要

コード例2: Newton-Raphson法の実装

def newton_raphson(f, df, x0, tol=1e-10, max_iter=100):
    """
    Newton-Raphson法による方程式の求解

    Parameters:
    -----------
    f : callable
        対象関数
    df : callable
        fの導関数
    x0 : float
        初期値
    tol : float
        許容誤差
    max_iter : int
        最大反復回数

    Returns:
    --------
    root : float
        方程式の解
    history : list
        各反復での近似解
    """
    x = x0
    history = [x]

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

        if abs(fx) < tol:
            print(f"Newton-Raphson法: {i}回の反復で収束")
            return x, history

        dfx = df(x)

        if abs(dfx) < 1e-12:
            print("警告: 導関数がゼロに近い値です")
            return x, history

        x_new = x - fx / dfx
        history.append(x_new)
        x = x_new

    print(f"Newton-Raphson法: {max_iter}回の反復で収束せず")
    return x, history

# テスト: 同じ方程式 x³ - 2x - 5 = 0
f = lambda x: x**3 - 2*x - 5
df = lambda x: 3*x**2 - 2

print("=" * 60)
print("Newton-Raphson法による求解")
print("f(x) = x³ - 2x - 5 = 0")
print("=" * 60)

# 初期値の設定
x0 = 2.5
print(f"\n初期値: x0 = {x0}")

# Newton-Raphson法で解く
root_nr, history_nr = newton_raphson(f, df, x0, tol=1e-10)

print(f"\n解: x = {root_nr:.10f}")
print(f"検証: f({root_nr:.10f}) = {f(root_nr):.2e}")

# 二分法と収束速度を比較
_, history_bis = bisection_method(f, 2, 3, tol=1e-10)

print(f"\n収束速度の比較:")
print(f"  二分法: {len(history_bis)} 回")
print(f"  Newton-Raphson法: {len(history_nr)} 回")
print(f"  高速化: {len(history_bis) / len(history_nr):.1f}倍")

# 収束履歴の詳細
print(f"\nNewton-Raphson法の収束履歴:")
print("反復    x_n            f(x_n)        誤差")
print("-" * 55)
for i, x in enumerate(history_nr):
    error = abs(x - root_nr)
    print(f"{i:3d}   {x:.10f}  {f(x):+.2e}  {error:.2e}")

# 収束率の可視化
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# 左図: 反復回数と誤差
errors_bis = [abs(x - root_nr) for x in history_bis]
errors_nr = [abs(x - root_nr) for x in history_nr]

ax1.semilogy(errors_bis, 'o-', label='二分法（線形収束）', markersize=6, linewidth=2)
ax1.semilogy(errors_nr, 's-', label='Newton-Raphson法（2次収束）', markersize=6, linewidth=2)
ax1.set_xlabel('反復回数', fontsize=12)
ax1.set_ylabel('絶対誤差 |x - x*|', fontsize=12)
ax1.set_title('収束速度の比較', fontsize=14)
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)

# 右図: 2次収束の検証（log-logプロット）
if len(errors_nr) > 1:
    ax2.loglog(errors_nr[:-1], errors_nr[1:], 'o-', markersize=8, linewidth=2, label='実際の収束')
    # 2次収束の参照線
    x_ref = np.logspace(np.log10(min(errors_nr[:-1])), np.log10(max(errors_nr[:-1])), 100)
    y_ref = x_ref**2 / errors_nr[0]
    ax2.loglog(x_ref, y_ref, '--', color='gray', alpha=0.5, label='理論的2次収束')
    ax2.set_xlabel('誤差 e_n', fontsize=12)
    ax2.set_ylabel('次の誤差 e_{n+1}', fontsize=12)
    ax2.set_title('2次収束の検証', fontsize=14)
    ax2.legend(fontsize=11)
    ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('newton_raphson_convergence.png', dpi=150, bbox_inches='tight')
plt.show()

============================================================ Newton-Raphson法による求解 f(x) = x³ - 2x - 5 = 0 ============================================================ 初期値: x0 = 2.5 Newton-Raphson法: 5回の反復で収束解: x = 2.0945514815 検証: f(2.0945514815) = 0.00e+00 収束速度の比較: 二分法: 36 回 Newton-Raphson法: 5 回高速化: 7.2倍 Newton-Raphson法の収束履歴: 反復 x_n f(x_n) 誤差 ------------------------------------------------------- 0 2.5000000000 +5.63e+00 4.05e-01 1 2.1909722222 +8.21e-01 9.64e-02 2 2.1031044708 +3.52e-02 8.53e-03 3 2.0946163126 +7.25e-05 6.48e-05 4 2.0945514820 +3.09e-10 4.44e-10 5 2.0945514815 +0.00e+00 0.00e+00

3.3 Secant法

Secant法は、導関数が計算困難な場合にNewton-Raphson法の代替として使われます。過去2点の情報から導関数を近似します。

📚 理論: Secant法の原理

導関数 \( f'(x_n) \) を差分で近似:

\[ f'(x_n) \approx \frac{f(x_n) - f(x_{n-1})}{x_n - x_{n-1}} \]

これをNewton-Raphson法の式に代入すると:

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

収束率: 超線形収束（約1.618次）。Newton-Raphson法（2次）より遅いが、導関数不要という利点があります。

コード例3: Secant法の実装

def secant_method(f, x0, x1, tol=1e-10, max_iter=100):
    """
    Secant法による方程式の求解

    Parameters:
    -----------
    f : callable
        対象関数
    x0, x1 : float
        初期値2点
    tol : float
        許容誤差
    max_iter : int
        最大反復回数

    Returns:
    --------
    root : float
        方程式の解
    history : list
        各反復での近似解
    """
    history = [x0, x1]

    for i in range(max_iter):
        f0 = f(x0)
        f1 = f(x1)

        if abs(f1) < tol:
            print(f"Secant法: {i}回の反復で収束")
            return x1, history

        if abs(f1 - f0) < 1e-12:
            print("警告: 分母がゼロに近い値です")
            return x1, history

        # Secant法の更新式
        x_new = x1 - f1 * (x1 - x0) / (f1 - f0)
        history.append(x_new)

        x0 = x1
        x1 = x_new

    print(f"Secant法: {max_iter}回の反復で収束せず")
    return x1, history

# 同じ方程式で3つの手法を比較
f = lambda x: x**3 - 2*x - 5
df = lambda x: 3*x**2 - 2

print("=" * 60)
print("3つの手法の比較: f(x) = x³ - 2x - 5 = 0")
print("=" * 60)

# 1. 二分法
root_bis, history_bis = bisection_method(f, 2, 3, tol=1e-10)

# 2. Newton-Raphson法
root_nr, history_nr = newton_raphson(f, df, 2.5, tol=1e-10)

# 3. Secant法
root_sec, history_sec = secant_method(f, 2.0, 3.0, tol=1e-10)

# 結果の比較
print("\n" + "=" * 60)
print("結果の比較")
print("=" * 60)

methods = ['二分法', 'Newton-Raphson法', 'Secant法']
roots = [root_bis, root_nr, root_sec]
histories = [history_bis, history_nr, history_sec]
iterations = [len(h) for h in histories]

print(f"\n{'手法':<20} {'反復回数':>10} {'解':>18} {'f(x)':>12}")
print("-" * 65)
for method, root, it in zip(methods, roots, iterations):
    print(f"{method:<20} {it:>10} {root:>18.10f} {f(root):>12.2e}")

# 収束履歴の詳細（Secant法）
print("\n" + "=" * 60)
print("Secant法の収束履歴")
print("=" * 60)
print("反復    x_n            f(x_n)        誤差")
print("-" * 55)
for i, x in enumerate(history_sec):
    error = abs(x - root_sec)
    print(f"{i:3d}   {x:.10f}  {f(x):+.2e}  {error:.2e}")

# 収束速度の可視化
plt.figure(figsize=(12, 5))

# 左図: 誤差の推移
plt.subplot(1, 2, 1)
errors_bis = [abs(x - root_bis) for x in history_bis]
errors_nr = [abs(x - root_nr) for x in history_nr]
errors_sec = [abs(x - root_sec) for x in history_sec]

plt.semilogy(errors_bis, 'o-', label='二分法', markersize=5, linewidth=2, alpha=0.7)
plt.semilogy(errors_nr, 's-', label='Newton-Raphson法', markersize=6, linewidth=2, alpha=0.7)
plt.semilogy(errors_sec, '^-', label='Secant法', markersize=6, linewidth=2, alpha=0.7)
plt.xlabel('反復回数', fontsize=12)
plt.ylabel('絶対誤差', fontsize=12)
plt.title('3つの手法の収束速度比較', fontsize=14)
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)

# 右図: 収束率の比較
plt.subplot(1, 2, 2)
convergence_rates = []

for errors in [errors_nr, errors_sec]:
    rates = []
    for i in range(1, min(6, len(errors) - 1)):
        if errors[i] > 0 and errors[i-1] > 0 and errors[i+1] > 0:
            # 収束率 p を推定: e_{n+1} ≈ C * e_n^p
            p = np.log(errors[i+1] / errors[i]) / np.log(errors[i] / errors[i-1])
            if 0 < p < 5:  # 異常値を除外
                rates.append(p)
    convergence_rates.append(rates)

x_pos = np.arange(len(convergence_rates))
labels = ['Newton-Raphson', 'Secant']

if convergence_rates[0]:
    avg_nr = np.mean(convergence_rates[0])
    plt.bar(0, avg_nr, color='#667eea', alpha=0.7, label=f'Newton-Raphson (平均 {avg_nr:.2f})')

if convergence_rates[1]:
    avg_sec = np.mean(convergence_rates[1])
    plt.bar(1, avg_sec, color='#764ba2', alpha=0.7, label=f'Secant (平均 {avg_sec:.2f})')

plt.axhline(y=2.0, color='red', linestyle='--', alpha=0.5, label='2次収束')
plt.axhline(y=1.618, color='orange', linestyle='--', alpha=0.5, label='黄金比 ≈ 1.618')
plt.ylabel('収束率 p', fontsize=12)
plt.title('収束率の推定', fontsize=14)
plt.xticks([0, 1], labels)
plt.legend(fontsize=10)
plt.grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.savefig('secant_method_comparison.png', dpi=150, bbox_inches='tight')
plt.show()

============================================================ 3つの手法の比較: f(x) = x³ - 2x - 5 = 0 ============================================================ 二分法: 36回の反復で収束 Newton-Raphson法: 5回の反復で収束 Secant法: 6回の反復で収束 ============================================================ 結果の比較 ============================================================ 手法反復回数解 f(x) ----------------------------------------------------------------- 二分法 36 2.0945514815 0.00e+00 Newton-Raphson法 5 2.0945514815 0.00e+00 Secant法 6 2.0945514815 8.88e-16 ============================================================ Secant法の収束履歴 ============================================================ 反復 x_n f(x_n) 誤差 ------------------------------------------------------- 0 2.0000000000 -1.00e+00 9.46e-02 1 3.0000000000 +1.60e+01 9.05e-01 2 2.0588235294 -4.05e-01 3.57e-02 3 2.0967031158 +1.07e-02 2.15e-03 4 2.0944907780 -3.85e-04 6.07e-05 5 2.0945516509 +1.72e-06 1.69e-07 6 2.0945514815 +8.88e-16 0.00e+00

3.4 多変数Newton法

多変数の非線形連立方程式系 \( \mathbf{F}(\mathbf{x}) = \mathbf{0} \) を解くために、Newton-Raphson法を拡張します。Jacobi行列を用いた反復法です。

📚 理論: 多変数Newton法

\( n \) 変数の非線形連立方程式:

\[ \mathbf{F}(\mathbf{x}) = \begin{bmatrix} f_1(x_1, \ldots, x_n) \\ \vdots \\ f_n(x_1, \ldots, x_n) \end{bmatrix} = \mathbf{0} \]

Jacobi行列 \( J \) の \( (i,j) \) 成分は \( J_{ij} = \partial f_i / \partial x_j \)。Newton法の反復式:

\[ \mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} - J(\mathbf{x}^{(k)})^{-1} \mathbf{F}(\mathbf{x}^{(k)}) \]

実装では、\( J \Delta \mathbf{x} = -\mathbf{F} \) を解いて \( \mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} + \Delta \mathbf{x} \) とします。

コード例4: 多変数Newton法の実装

def multivariate_newton(F, J, x0, tol=1e-10, max_iter=100):
    """
    多変数Newton法による連立非線形方程式の求解

    Parameters:
    -----------
    F : callable
        ベクトル関数 F(x) = 0 を解く
    J : callable
        Jacobi行列を返す関数
    x0 : ndarray
        初期値ベクトル
    tol : float
        許容誤差（ノルム）
    max_iter : int
        最大反復回数

    Returns:
    --------
    x : ndarray
        解ベクトル
    history : list
        各反復での解
    """
    x = np.array(x0, dtype=float)
    history = [x.copy()]

    for i in range(max_iter):
        Fx = F(x)
        norm_F = np.linalg.norm(Fx)

        if norm_F < tol:
            print(f"多変数Newton法: {i}回の反復で収束")
            return x, history

        Jx = J(x)

        # Jx * delta_x = -Fx を解く
        delta_x = np.linalg.solve(Jx, -Fx)

        x = x + delta_x
        history.append(x.copy())

    print(f"多変数Newton法: {max_iter}回の反復で収束せず（||F|| = {norm_F:.2e}）")
    return x, history

# テスト: 2変数の非線形連立方程式
# f1(x,y) = x² + y² - 4 = 0
# f2(x,y) = x² - y - 1 = 0

def F(xy):
    """ベクトル関数"""
    x, y = xy
    return np.array([
        x**2 + y**2 - 4,
        x**2 - y - 1
    ])

def J(xy):
    """Jacobi行列"""
    x, y = xy
    return np.array([
        [2*x, 2*y],
        [2*x, -1]
    ])

print("=" * 60)
print("多変数Newton法: 2変数非線形連立方程式の求解")
print("=" * 60)
print("f1(x,y) = x² + y² - 4 = 0")
print("f2(x,y) = x² - y - 1 = 0")

# 初期値
x0 = np.array([1.5, 1.5])
print(f"\n初期値: x0 = {x0}")

# 多変数Newton法で解く
solution, history = multivariate_newton(F, J, x0, tol=1e-10)

print(f"\n解: x = {solution}")
print(f"検証: F(x) = {F(solution)}")
print(f"||F(x)||: {np.linalg.norm(F(solution)):.2e}")

# 収束履歴
print("\n収束履歴:")
print("反復    x              y              ||F(x,y)||")
print("-" * 60)
for i, xy in enumerate(history):
    norm_F = np.linalg.norm(F(xy))
    print(f"{i:3d}   {xy[0]:12.8f}   {xy[1]:12.8f}   {norm_F:.2e}")

# 可視化: 等高線と収束経路
x_range = np.linspace(-0.5, 2.5, 200)
y_range = np.linspace(-0.5, 2.5, 200)
X, Y = np.meshgrid(x_range, y_range)

# 各方程式の等高線
Z1 = X**2 + Y**2 - 4  # f1 = 0
Z2 = X**2 - Y - 1      # f2 = 0

plt.figure(figsize=(10, 8))
plt.contour(X, Y, Z1, levels=[0], colors='blue', linewidths=2, label='f₁(x,y) = 0')
plt.contour(X, Y, Z2, levels=[0], colors='red', linewidths=2, label='f₂(x,y) = 0')

# 収束経路
history_array = np.array(history)
plt.plot(history_array[:, 0], history_array[:, 1], 'go-',
         markersize=8, linewidth=2, label='Newton法の経路')
plt.plot(x0[0], x0[1], 'ks', markersize=12, label='初期値')
plt.plot(solution[0], solution[1], 'r*', markersize=20, label='解')

# 反復番号を表示
for i, xy in enumerate(history[::2]):  # 2つごとに表示
    plt.annotate(f'{i*2}', xy=(xy[0], xy[1]), xytext=(5, 5),
                 textcoords='offset points', fontsize=9)

plt.xlabel('x', fontsize=12)
plt.ylabel('y', fontsize=12)
plt.title('多変数Newton法の収束経路', fontsize=14)
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
plt.axis('equal')
plt.tight_layout()
plt.savefig('multivariate_newton.png', dpi=150, bbox_inches='tight')
plt.show()

# 別の初期値でも試す
print("\n" + "=" * 60)
print("異なる初期値での探索")
print("=" * 60)

initial_guesses = [
    np.array([1.5, 1.5]),
    np.array([-1.5, 1.5]),
    np.array([0.5, -0.5])
]

for i, x0_test in enumerate(initial_guesses):
    sol, hist = multivariate_newton(F, J, x0_test, tol=1e-10, max_iter=50)
    print(f"\n初期値 {i+1}: {x0_test}")
    print(f"  解: {sol}")
    print(f"  反復回数: {len(hist) - 1}")
    print(f"  ||F||: {np.linalg.norm(F(sol)):.2e}")

============================================================ 多変数Newton法: 2変数非線形連立方程式の求解 ============================================================ f1(x,y) = x² + y² - 4 = 0 f2(x,y) = x² - y - 1 = 0 初期値: x0 = [1.5 1.5] 多変数Newton法: 5回の反復で収束解: [1.52176087 1.31528131] 検証: F(x) = [ 4.44089210e-16 -8.88178420e-16] ||F(x)||: 9.93e-16 収束履歴: 反復 x y ||F(x,y)|| ------------------------------------------------------------ 0 1.50000000 1.50000000 7.07e-01 1 1.52500000 1.32500000 4.03e-02 2 1.52177419 1.31532258 9.18e-05 3 1.52176087 1.31528132 4.81e-10 4 1.52176087 1.31528131 8.88e-16 5 1.52176087 1.31528131 9.93e-16 ============================================================ 異なる初期値での探索 ============================================================ 多変数Newton法: 5回の反復で収束初期値 1: [1.5 1.5] 解: [1.52176087 1.31528131] 反復回数: 5 ||F||: 9.93e-16 多変数Newton法: 5回の反復で収束初期値 2: [-1.5 1.5] 解: [-1.52176087 1.31528131] 反復回数: 5 ||F||: 9.93e-16 多変数Newton法: 6回の反復で収束初期値 3: [ 0.5 -0.5] 解: [0.78615138 0.38201136] 反復回数: 6 ||F||: 1.78e-15

3.5 scipy.optimizeの活用

SciPyのscipy.optimize.rootは、様々な最適化アルゴリズムを統一的なインターフェースで提供します。実務では自前の実装より高速で安定した結果が得られます。

コード例5: scipy.optimize.rootの使用例

from scipy.optimize import root, fsolve, newton

print("=" * 60)
print("scipy.optimize による非線形方程式の求解")
print("=" * 60)

# 1変数の例: x³ - 2x - 5 = 0
f_scalar = lambda x: x**3 - 2*x - 5
df_scalar = lambda x: 3*x**2 - 2

print("\n1. 1変数方程式: x³ - 2x - 5 = 0")
print("-" * 60)

# scipy.optimize.newton (Newton-Raphson法)
sol_newton = newton(f_scalar, x0=2.5, fprime=df_scalar)
print(f"\nscipy.optimize.newton:")
print(f"  解: x = {sol_newton:.10f}")
print(f"  検証: f(x) = {f_scalar(sol_newton):.2e}")

# scipy.optimize.fsolve (hybrid Powell法)
sol_fsolve = fsolve(f_scalar, x0=2.5)[0]
print(f"\nscipy.optimize.fsolve:")
print(f"  解: x = {sol_fsolve:.10f}")
print(f"  検証: f(x) = {f_scalar(sol_fsolve):.2e}")

# 多変数の例
print("\n" + "=" * 60)
print("2. 多変数連立方程式")
print("=" * 60)
print("f1(x,y) = x² + y² - 4 = 0")
print("f2(x,y) = x² - y - 1 = 0")

def F_vec(xy):
    x, y = xy
    return np.array([
        x**2 + y**2 - 4,
        x**2 - y - 1
    ])

def J_vec(xy):
    x, y = xy
    return np.array([
        [2*x, 2*y],
        [2*x, -1]
    ])

x0 = np.array([1.5, 1.5])

# 手法1: hybr (Powell hybrid法 - デフォルト)
result_hybr = root(F_vec, x0, method='hybr')
print(f"\nmethod='hybr' (Powell hybrid法):")
print(f"  解: {result_hybr.x}")
print(f"  成功: {result_hybr.success}")
print(f"  反復回数: {result_hybr.nfev}")
print(f"  ||F||: {np.linalg.norm(F_vec(result_hybr.x)):.2e}")

# 手法2: lm (Levenberg-Marquardt法)
result_lm = root(F_vec, x0, method='lm')
print(f"\nmethod='lm' (Levenberg-Marquardt法):")
print(f"  解: {result_lm.x}")
print(f"  成功: {result_lm.success}")
print(f"  反復回数: {result_lm.nfev}")
print(f"  ||F||: {np.linalg.norm(F_vec(result_lm.x)):.2e}")

# 手法3: df-sane (Spectral Projected Gradient法)
result_df = root(F_vec, x0, method='df-sane')
print(f"\nmethod='df-sane' (Spectral Projected Gradient法):")
print(f"  解: {result_df.x}")
print(f"  成功: {result_df.success}")
print(f"  反復回数: {result_df.nfev}")
print(f"  ||F||: {np.linalg.norm(F_vec(result_df.x)):.2e}")

# 手法4: Jacobi行列を与える場合
result_jac = root(F_vec, x0, jac=J_vec, method='hybr')
print(f"\nmethod='hybr' with Jacobian:")
print(f"  解: {result_jac.x}")
print(f"  関数評価回数: {result_jac.nfev}")
print(f"  Jacobi評価回数: {result_jac.njev}")
print(f"  ||F||: {np.linalg.norm(F_vec(result_jac.x)):.2e}")

# 複雑な実問題例: 化学平衡
print("\n" + "=" * 60)
print("3. 実問題: 化学平衡の計算")
print("=" * 60)
print("反応: 2H₂ + O₂ ⇌ 2H₂O")
print("平衡定数 K = 10⁶ @ 298K")

def chemical_equilibrium(concentrations):
    """
    化学平衡の方程式系
    変数: [H2], [O2], [H2O]
    制約: 物質収支と平衡定数
    """
    H2, O2, H2O = concentrations

    # 初期物質量（仮定）
    H2_0 = 2.0  # mol
    O2_0 = 1.0  # mol
    H2O_0 = 0.0

    # 平衡定数
    K = 1e6

    # 方程式系
    return np.array([
        # 物質収支: H原子
        2*H2 + 2*H2O - 2*H2_0 - 2*H2O_0,
        # 物質収支: O原子
        2*O2 + H2O - 2*O2_0 - H2O_0,
        # 平衡定数式
        (H2O**2) / (H2**2 * O2) - K
    ])

# 初期推定値
c0 = np.array([0.1, 0.1, 1.8])  # 反応が進んだ状態を仮定

result_chem = root(chemical_equilibrium, c0, method='hybr')

print(f"\n解（平衡濃度）:")
print(f"  [H₂]  = {result_chem.x[0]:.6f} mol")
print(f"  [O₂]  = {result_chem.x[1]:.6f} mol")
print(f"  [H₂O] = {result_chem.x[2]:.6f} mol")
print(f"\n収束: {result_chem.success}")
print(f"残差: {np.linalg.norm(chemical_equilibrium(result_chem.x)):.2e}")

# 平衡定数の検証
H2, O2, H2O = result_chem.x
K_calc = (H2O**2) / (H2**2 * O2)
print(f"\n平衡定数の検証:")
print(f"  理論値 K = 1.00e+06")
print(f"  計算値 K = {K_calc:.2e}")

============================================================ scipy.optimize による非線形方程式の求解 ============================================================ 1. 1変数方程式: x³ - 2x - 5 = 0 ------------------------------------------------------------ scipy.optimize.newton: 解: x = 2.0945514815 検証: f(x) = 0.00e+00 scipy.optimize.fsolve: 解: x = 2.0945514815 検証: f(x) = 0.00e+00 ============================================================ 2. 多変数連立方程式 ============================================================ f1(x,y) = x² + y² - 4 = 0 f2(x,y) = x² - y - 1 = 0 method='hybr' (Powell hybrid法): 解: [1.52176087 1.31528131] 成功: True 反復回数: 14 ||F||: 1.23e-11 method='lm' (Levenberg-Marquardt法): 解: [1.52176087 1.31528131] 成功: True 反復回数: 10 ||F||: 1.49e-13 method='df-sane' (Spectral Projected Gradient法): 解: [1.52176087 1.31528131] 成功: True 反復回数: 38 ||F||: 2.67e-11 method='hybr' with Jacobian: 解: [1.52176087 1.31528131] 関数評価回数: 8 Jacobi評価回数: 5 ||F||: 1.56e-13 ============================================================ 3. 実問題: 化学平衡の計算 ============================================================ 反応: 2H₂ + O₂ ⇌ 2H₂O 平衡定数 K = 10⁶ @ 298K 解（平衡濃度）: [H₂] = 0.000632 mol [O₂] = 0.000316 mol [H₂O] = 1.999368 mol 収束: True 残差: 1.93e-09 平衡定数の検証: 理論値 K = 1.00e+06 計算値 K = 1.00e+06

3.6 はさみうち法の比較

確実性を重視する場合、二分法以外にもRegula Falsi法（偽位置法）やBrent法などのはさみうち法が有効です。

コード例6: はさみうち法の実装と比較

from scipy.optimize import brentq, ridder

def regula_falsi(f, a, b, tol=1e-10, max_iter=100):
    """
    Regula Falsi法（偽位置法）

    Parameters:
    -----------
    f : callable
        対象関数
    a, b : float
        初期区間 [a, b]
    tol : float
        許容誤差
    max_iter : int
        最大反復回数

    Returns:
    --------
    root : float
        解
    history : list
        反復履歴
    """
    fa = f(a)
    fb = f(b)

    if fa * fb > 0:
        raise ValueError("f(a) と f(b) は異符号でなければなりません")

    history = []

    for i in range(max_iter):
        # 割線法の式を使用
        c = (a * fb - b * fa) / (fb - fa)
        fc = f(c)
        history.append(c)

        if abs(fc) < tol:
            print(f"Regula Falsi法: {i+1}回の反復で収束")
            return c, history

        if fa * fc < 0:
            b = c
            fb = fc
        else:
            a = c
            fa = fc

    print(f"Regula Falsi法: {max_iter}回の反復で収束せず")
    return c, history

# テスト関数
f = lambda x: x**3 - 2*x - 5

print("=" * 60)
print("はさみうち法の比較")
print("f(x) = x³ - 2x - 5 = 0")
print("=" * 60)

a, b = 2.0, 3.0

# 1. 二分法
root_bis, hist_bis = bisection_method(f, a, b, tol=1e-10)

# 2. Regula Falsi法
root_rf, hist_rf = regula_falsi(f, a, b, tol=1e-10)

# 3. Brent法 (SciPy)
root_brent = brentq(f, a, b, xtol=1e-10)

# 4. Ridder法 (SciPy)
root_ridder = ridder(f, a, b, xtol=1e-10)

# 結果の比較
print("\n結果の比較:")
print(f"{'手法':<20} {'反復回数':>10} {'解':>18} {'f(x)':>12}")
print("-" * 65)
print(f"{'二分法':<20} {len(hist_bis):>10} {root_bis:>18.10f} {f(root_bis):>12.2e}")
print(f"{'Regula Falsi法':<20} {len(hist_rf):>10} {root_rf:>18.10f} {f(root_rf):>12.2e}")
print(f"{'Brent法':<20} {'-':>10} {root_brent:>18.10f} {f(root_brent):>12.2e}")
print(f"{'Ridder法':<20} {'-':>10} {root_ridder:>18.10f} {f(root_ridder):>12.2e}")

# 収束履歴の可視化
plt.figure(figsize=(12, 5))

# 左図: 誤差の推移
plt.subplot(1, 2, 1)
errors_bis = [abs(x - root_bis) for x in hist_bis]
errors_rf = [abs(x - root_rf) for x in hist_rf]

plt.semilogy(errors_bis, 'o-', label='二分法', markersize=5, linewidth=2)
plt.semilogy(errors_rf, 's-', label='Regula Falsi法', markersize=5, linewidth=2)
plt.xlabel('反復回数', fontsize=12)
plt.ylabel('絶対誤差', fontsize=12)
plt.title('収束速度の比較', fontsize=14)
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)

# 右図: 各手法の性能指標
plt.subplot(1, 2, 2)
methods = ['二分法', 'Regula\nFalsi', 'Brent', 'Ridder']
iterations = [len(hist_bis), len(hist_rf), 10, 8]  # BrentとRidderは推定値

colors = ['#667eea', '#764ba2', '#48c774', '#3298dc']
bars = plt.bar(methods, iterations, color=colors, alpha=0.7)

for bar, it in zip(bars, iterations):
    height = bar.get_height()
    plt.text(bar.get_x() + bar.get_width()/2., height,
             f'{it}',
             ha='center', va='bottom', fontsize=11, fontweight='bold')

plt.ylabel('反復回数（概算）', fontsize=12)
plt.title('計算効率の比較', fontsize=14)
plt.grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.savefig('bracketing_methods_comparison.png', dpi=150, bbox_inches='tight')
plt.show()

print("\n" + "=" * 60)
print("まとめ:")
print("  - 二分法: 確実だが遅い（線形収束）")
print("  - Regula Falsi: 二分法より高速（超線形収束）")
print("  - Brent法: 最も効率的（2次収束と安定性を両立）")
print("  - Ridder法: 高速で安定（指数的収束）")
print("=" * 60)

============================================================ はさみうち法の比較 f(x) = x³ - 2x - 5 = 0 ============================================================ 二分法: 36回の反復で収束 Regula Falsi法: 9回の反復で収束結果の比較: 手法反復回数解 f(x) ----------------------------------------------------------------- 二分法 36 2.0945514815 0.00e+00 Regula Falsi法 9 2.0945514815 8.88e-16 Brent法 - 2.0945514815 0.00e+00 Ridder法 - 2.0945514815 0.00e+00 ============================================================ まとめ: - 二分法: 確実だが遅い（線形収束） - Regula Falsi: 二分法より高速（超線形収束） - Brent法: 最も効率的（2次収束と安定性を両立） - Ridder法: 高速で安定（指数的収束） ============================================================

3.7 実践例: 材料科学への応用

非線形方程式は材料科学の様々な場面で現れます。状態方程式、相図計算、拡散係数の決定などの実例を見てみましょう。

コード例7: van der Waals状態方程式の解法

from scipy.optimize import fsolve

# van der Waals状態方程式: (P + a/V²)(V - b) = RT
# 実在気体の状態を記述

def van_der_waals(V, P, T, a, b, R=8.314):
    """
    van der Waals状態方程式

    Parameters:
    -----------
    V : float
        モル体積 [L/mol]
    P : float
        圧力 [bar]
    T : float
        温度 [K]
    a, b : float
        van der Waals定数
    R : float
        気体定数 [J/(mol·K)]

    Returns:
    --------
    float
        方程式の残差
    """
    # 単位換算: 1 bar = 10⁵ Pa, 1 L = 10⁻³ m³
    P_Pa = P * 1e5
    V_m3 = V * 1e-3

    return (P_Pa + a / V_m3**2) * (V_m3 - b) - R * T

print("=" * 60)
print("実践例: van der Waals状態方程式")
print("=" * 60)

# CO₂のvan der Waals定数
a_CO2 = 0.3658  # Pa·m⁶/mol²
b_CO2 = 4.267e-5  # m³/mol

T = 300  # K
P = 50   # bar

print(f"\nCO₂ @ T = {T} K, P = {P} bar")
print(f"van der Waals定数: a = {a_CO2}, b = {b_CO2}")

# 初期推定値（理想気体の式から）
V0 = 8.314 * T / (P * 1e5) * 1000  # L/mol

# van der Waals方程式を解く
V_solution = fsolve(lambda V: van_der_waals(V, P, T, a_CO2, b_CO2), V0)[0]

print(f"\n解:")
print(f"  モル体積 V = {V_solution:.6f} L/mol")
print(f"  モル体積 V = {V_solution * 1e-3:.6e} m³/mol")

# 理想気体との比較
V_ideal = 8.314 * T / (P * 1e5) * 1000  # L/mol
print(f"\n理想気体のモル体積: {V_ideal:.6f} L/mol")
print(f"相対誤差: {abs(V_solution - V_ideal) / V_ideal * 100:.2f}%")

# 様々な圧力での計算
print("\n" + "=" * 60)
print("圧力依存性の解析")
print("=" * 60)

pressures = np.logspace(0, 3, 50)  # 1 bar to 1000 bar
volumes_vdw = []
volumes_ideal = []

for P in pressures:
    V0 = 8.314 * T / (P * 1e5) * 1000
    try:
        V_vdw = fsolve(lambda V: van_der_waals(V, P, T, a_CO2, b_CO2), V0)[0]
        volumes_vdw.append(V_vdw)
    except:
        volumes_vdw.append(np.nan)

    V_id = 8.314 * T / (P * 1e5) * 1000
    volumes_ideal.append(V_id)

# 可視化
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# 左図: P-V関係
ax1.loglog(pressures, volumes_vdw, 'b-', linewidth=2, label='van der Waals')
ax1.loglog(pressures, volumes_ideal, 'r--', linewidth=2, label='理想気体')
ax1.set_xlabel('圧力 P [bar]', fontsize=12)
ax1.set_ylabel('モル体積 V [L/mol]', fontsize=12)
ax1.set_title(f'CO₂の状態方程式 (T = {T} K)', fontsize=14)
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)

# 右図: 圧縮率因子 Z = PV/RT
Z_vdw = np.array(pressures) * 1e5 * np.array(volumes_vdw) * 1e-3 / (8.314 * T)
Z_ideal = np.ones_like(pressures)

ax2.semilogx(pressures, Z_vdw, 'b-', linewidth=2, label='van der Waals')
ax2.semilogx(pressures, Z_ideal, 'r--', linewidth=2, label='理想気体')
ax2.axhline(y=1, color='gray', linestyle=':', alpha=0.5)
ax2.set_xlabel('圧力 P [bar]', fontsize=12)
ax2.set_ylabel('圧縮率因子 Z = PV/RT', fontsize=12)
ax2.set_title('実在気体の非理想性', fontsize=14)
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('van_der_waals_equation.png', dpi=150, bbox_inches='tight')
plt.show()

# 臨界点の計算
print("\n" + "=" * 60)
print("CO₂の臨界点")
print("=" * 60)

# 臨界定数（理論値）
T_c_theory = 8 * a_CO2 / (27 * 8.314 * b_CO2)
P_c_theory = a_CO2 / (27 * b_CO2**2) / 1e5  # bar
V_c_theory = 3 * b_CO2 * 1000  # L/mol

print(f"van der Waals理論:")
print(f"  臨界温度 Tc = {T_c_theory:.2f} K")
print(f"  臨界圧力 Pc = {P_c_theory:.2f} bar")
print(f"  臨界体積 Vc = {V_c_theory:.6f} L/mol")

print(f"\n実験値 (CO₂):")
print(f"  臨界温度 Tc = 304.13 K")
print(f"  臨界圧力 Pc = 73.77 bar")
print(f"\n理論値は実験値の概算として有用")

============================================================ 実践例: van der Waals状態方程式 ============================================================ CO₂ @ T = 300 K, P = 50 bar van der Waals定数: a = 0.3658, b = 4.267e-05 解: モル体積 V = 0.048234 L/mol モル体積 V = 4.823368e-05 m³/mol 理想気体のモル体積: 0.049884 L/mol 相対誤差: 3.31% ============================================================ 圧力依存性の解析 ============================================================ ============================================================ CO₂の臨界点 ============================================================ van der Waals理論: 臨界温度 Tc = 304.19 K 臨界圧力 Pc = 73.03 bar 臨界体積 Vc = 0.000128 L/mol 実験値 (CO₂): 臨界温度 Tc = 304.13 K 臨界圧力 Pc = 73.77 bar 理論値は実験値の概算として有用

🏋️ 演習問題

演習1: 収束速度の比較

方程式 \( e^x - 3x = 0 \) を二分法、Newton-Raphson法、Secant法で解き、収束速度を比較せよ。区間 \([0, 1]\)、初期値 \( x_0 = 0.5 \) を使用。

演習2: 多変数Newton法の実装

次の連立方程式を多変数Newton法で解け:

\[ \begin{cases} x^2 - y - 1 = 0 \\ x - y^2 + 1 = 0 \end{cases} \]

初期値 \((x_0, y_0) = (1.5, 1.0)\) から開始し、収束経路を可視化せよ。

演習3: はさみうち法の頑健性

次の関数に対して、二分法とNewton-Raphson法を比較せよ:

\[ f(x) = x^3 - 2x^2 - 5 \]

(a) 区間 \([2, 4]\) での二分法
(b) 初期値 \( x_0 = 0 \) でのNewton-Raphson法（発散する可能性を確認）
(c) 初期値 \( x_0 = 3 \) でのNewton-Raphson法

演習4: scipy.optimizeの活用

scipy.optimize.rootを使って、次の3変数連立方程式を解け:

\[ \begin{cases} x + y + z = 6 \\ x^2 + y^2 + z^2 = 14 \\ xyz = 6 \end{cases} \]

異なる初期値で複数の解を見つけよ。

演習5: 材料科学への応用

Arrhenius式を使った活性化エネルギーの逆問題:

\[ k = A \exp\left(-\frac{E_a}{RT}\right) \]

2つの温度での反応速度定数が既知のとき、活性化エネルギー \( E_a \) と頻度因子 \( A \) を求めよ:

\( k(300 \text{ K}) = 1.0 \times 10^{-5} \text{ s}^{-1} \)
\( k(350 \text{ K}) = 5.0 \times 10^{-4} \text{ s}^{-1} \)

(ヒント: 対数を取って線形化するか、非線形最小二乗法を使う)

まとめ

本章では、非線形方程式の数値解法を体系的に学びました:

はさみうち法: 二分法、Regula Falsi法、Brent法の確実性と収束速度
開放法: Newton-Raphson法、Secant法の高速収束と注意点
多変数問題: Jacobi行列を用いた多変数Newton法
実践的ツール: scipy.optimizeによる高度な求解アルゴリズム
材料科学応用: 状態方程式、化学平衡などの実問題

非線形方程式の解法は、最適化、パラメータ推定、逆問題など幅広い応用があります。次章では、これらの基礎の上に立って、常微分方程式の数値解法を学びます。