SPICE手法統合シミュレータ

Pythonコード

# --- 総合SPICE手法統合シミュレータ / Full SPICE Simulation Python 合体プログラム ---
# Program Name: spice_solver_methods_demo.py
import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import csc_matrix
from scipy.sparse.linalg import splu
from sympy import symbols, diff, lambdify

# -------------------------------
# ① 疎行列LU分解（Sparse LU Decomposition）
# -------------------------------
A_dense = np.array([
    [10.0, 0.0, 2.0],
    [3.0, 9.0, 0.0],
    [0.0, 7.0, 8.0]
])
b = np.array([7.85, -19.3, 71.4])

A_sparse = csc_matrix(A_dense)     # 疎行列に変換 / Convert to sparse format
lu_sparse = splu(A_sparse)         # 疎LU分解 / Sparse LU factorization
x_sparse_lu = lu_sparse.solve(b)   # 解を得る / Solve Ax = b

# -------------------------------
# ② ニュートン・ラフソン法（ヤコビアン自動構築）
# -------------------------------
x_sym = symbols('x')
f_expr = x_sym**3 - x_sym - 2
df_expr = diff(f_expr, x_sym)

f_func = lambdify(x_sym, f_expr, modules='numpy')
df_func = lambdify(x_sym, df_expr, modules='numpy')

x0 = 1.5
tol = 1e-6
max_iter = 10
x_vals = [x0]

for i in range(max_iter):
    x_new = x0 - f_func(x0) / df_func(x0)
    x_vals.append(x_new)
    if abs(x_new - x0) < tol:
        break
    x0 = x_new

# -------------------------------
# ③ 過渡解析：台形法 vs 後退オイラー法 / Trapezoidal vs Backward Euler
# -------------------------------
a = 1.0   # 減衰係数 / Decay constant
dt = 0.1  # ステップ幅 / Time step
T = 1.0   # 総時間 / Total time
N = int(T / dt)
t_vals = np.linspace(0, T, N+1)

# 初期化 / Initialization
x_trap = np.zeros(N+1)
x_euler = np.zeros(N+1)
x_trap[0] = 1.0
x_euler[0] = 1.0

for n in range(N):
    x_trap[n+1] = ((1 - a*dt/2) / (1 + a*dt/2)) * x_trap[n]     # 台形法 / Trapezoidal
    x_euler[n+1] = x_euler[n] / (1 + a*dt)                      # 後退オイラー / Backward Euler

# -------------------------------
# 可視化 / Visualization
# -------------------------------
plt.figure(figsize=(14, 10))

# ① LU分解結果 / Sparse LU Decomposition
plt.subplot(3, 1, 1)
plt.bar(['x1', 'x2', 'x3'], x_sparse_lu)
plt.title("① Sparse LU Decomposition (for large circuit matrices)")
plt.ylabel("Value")
plt.grid(True)

# ② ニュートン法収束過程 / Newton-Raphson Method
plt.subplot(3, 1, 2)
plt.plot(range(len(x_vals)), x_vals, marker='o')
plt.title("② Newton-Raphson Iteration with Auto Jacobian")
plt.xlabel("Iteration")
plt.ylabel("x value")
plt.grid(True)

# ③ 台形法 vs 後退オイラー法 / Trapezoidal vs Backward Euler
plt.subplot(3, 1, 3)
plt.plot(t_vals, x_trap, marker='o', label='Trapezoidal')
plt.plot(t_vals, x_euler, marker='s', label='Backward Euler')
plt.title("③ Transient Analysis: Trapezoidal vs Backward Euler")
plt.xlabel("Time [s]")
plt.ylabel("x(t)")
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

結果