spiceの数値計算をPythonコードで


https://qiita.com/sukimaengineer/items/a0694120001060ca0055

### ✅ 数値シミュレーション手法一覧表（陰的法・DAE中心）

| 分類         | 手法名                    | 主な特徴                   | 用途・備考                     |
| ---------- | ---------------------- | ---------------------- | ------------------------- |
| **数値積分法**  | 後退オイラー法                | 一次精度・陰的・安定性が高い         | 硬い問題に適用しやすく、SPICEの基本手法    |
|            | 台形公式（Trapezoidal Rule） | 二次精度・平均化を利用・比較的高精度     | 高精度解析向けだが振動的になることがある      |
|            | Gear法（BDF法）            | 多段・陰的・L-安定性・非常に硬い問題に有効 | 電気回路（SPICE）で頻繁に使われる       |
| **DAEソルバ** | IDA（SUNDIALS）          | Index-1対応・陰的法・ヤコビ評価に柔軟 | 汎用DAEソルバ、工学的モデリングツールで広く利用 |
|            | DASSL                  | Fortran由来・古典的DAE解法     | 数値計算書で頻出                  |
| **反復法**    | ニュートン・ラフソン法            | 非線形方程式用・局所二次収束         | ヤコビ行列が必要、初期値依存が大          |
|            | 修正ニュートン法               | ヤコビ行列を固定・反復回数節約        | 大規模系で計算軽減に有効              |
| **ステップ制御** | 予測子-修正子方式              | 誤差推定 → ステップ幅調整         | Gear法やSPICEで標準的に採用        |
| **行列解法**   | LU分解                   | 線形系の直接解法・安定            | 一度分解すれば多回代入に向いている         |
|            | ガウス消去法                 | 線形方程式の基本解法             | 小規模システムに有効                |
| **応用系**    | MNA法                   | 電気回路→Index-1 DAEへ導出    | SPICEの基礎理論                |
|            | DCOP解析                 | 静的状態から初期条件決定           | DAE整合初期値生成として不可欠          |

---

```Python
# Program Name: first_order_compare_error.py
# Creation Date: 20250609
# Overview: Compare 3 integration methods for RC circuit including error analysis
# Usage: Run to see method comparison, absolute error plots, and dt-dependence

!pip install numpy matplotlib

import numpy as np
import matplotlib.pyplot as plt

# --- Simulation parameters / パラメータ ---
R = 1.0
C = 1.0
tau = R * C
T_end = 5.0
V_in = 1.0
V0 = 0.0

# --- 固定ステップでの誤差評価 ---
dt = 0.1
time = np.arange(0, T_end + dt, dt)
N = len(time)

# 初期化 / Initialize
V_be = np.zeros(N)
V_tp = np.zeros(N)
V_gear = np.zeros(N)
V_ana = V_in * (1 - np.exp(-time / tau))
V_be[0] = V_tp[0] = V_gear[0] = V0
V_gear[1] = (V0 + dt / tau * V_in) / (1 + dt / tau)

# Backward Euler
for n in range(1, N):
    V_be[n] = (V_be[n-1] + dt / tau * V_in) / (1 + dt / tau)

# Trapezoidal
for n in range(1, N):
    a = (1 - dt / (2 * tau))
    b = dt / tau
    V_tp[n] = (a * V_tp[n-1] + b * V_in) / (1 + dt / (2 * tau))

# Gear BDF2
for n in range(2, N):
    alpha0, alpha1, alpha2 = 3/2, -2, 1/2
    beta = dt * V_in / tau
    V_gear[n] = (-alpha1 * V_gear[n-1] - alpha2 * V_gear[n-2] + beta) / alpha0

# --- 絶対誤差のプロット / Plot absolute errors ---
plt.plot(time, np.abs(V_be - V_ana), label='Error: Backward Euler')
plt.plot(time, np.abs(V_tp - V_ana), label='Error: Trapezoidal')
plt.plot(time, np.abs(V_gear - V_ana), label='Error: Gear BDF2')
plt.title('Absolute Error vs Time')
plt.xlabel('Time [s]')
plt.ylabel('Absolute Error [V]')
plt.legend()
plt.grid(True)
plt.show()

# --- ステップ幅と誤差の関係 / dt vs error plot ---
dt_list = [0.4, 0.2, 0.1, 0.05, 0.01]
err_be = []
err_tp = []
err_gear = []

for dt in dt_list:
    time = np.arange(0, T_end + dt, dt)
    N = len(time)
    V_ana = V_in * (1 - np.exp(-time / tau))
    V_be = np.zeros(N)
    V_tp = np.zeros(N)
    V_gear = np.zeros(N)
    V_be[0] = V_tp[0] = V_gear[0] = V0
    if N > 1:
        V_gear[1] = (V0 + dt / tau * V_in) / (1 + dt / tau)

    for n in range(1, N):
        V_be[n] = (V_be[n-1] + dt / tau * V_in) / (1 + dt / tau)
        a = (1 - dt / (2 * tau))
        b = dt / tau
        V_tp[n] = (a * V_tp[n-1] + b * V_in) / (1 + dt / (2 * tau))

    for n in range(2, N):
        alpha0, alpha1, alpha2 = 3/2, -2, 1/2
        beta = dt * V_in / tau
        V_gear[n] = (-alpha1 * V_gear[n-1] - alpha2 * V_gear[n-2] + beta) / alpha0

    err_be.append(np.abs(V_be - V_ana).max())
    err_tp.append(np.abs(V_tp - V_ana).max())
    err_gear.append(np.abs(V_gear - V_ana).max())

# --- プロット（対数軸） / Error vs dt ---
plt.loglog(dt_list, err_be, 'o-', label='Backward Euler')
plt.loglog(dt_list, err_tp, 's--', label='Trapezoidal')
plt.loglog(dt_list, err_gear, 'x-.', label='Gear BDF2')
plt.title('Max Error vs Time Step')
plt.xlabel('Time Step Δt [s]')
plt.ylabel('Max Absolute Error [V]')
plt.grid(True, which='both')
plt.legend()
plt.show()

```