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

✅ 数値シミュレーション手法一覧表（陰的法・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整合初期値生成として不可欠

---

# 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()