Transient Analysis & Time-Stepping Architecture

This document explains the mathematical architecture behind PyCircuit’s transient simulation engine, specifically detailing the integrators, nonlinear solvers, and adaptive time-stepping methods.

Integrators (Differentiators)

PyCircuit solves Differential Algebraic Equations (DAEs) of the form:

\[f(v(t)) + \frac{d}{dt} q(v(t)) + u(t) = 0\]

To solve this numerically, the continuous time derivative must be replaced by a discrete-time approximation. This is the job of the Integrator (often referred to as a differentiator in MNA context).

PyCircuit supports two primary implicit integration methods:

  • Trapezoidal Rule: A 2nd-order method that uses a linear interpolation of the derivative. It is highly accurate and preserves energy in non-dissipative circuits (like LC tanks) but can suffer from numerical “ringing” (oscillations) during sharp discontinuous transitions.

  • Gear-2 (BDF2): A 2nd-order Backward Differentiation Formula. It is strictly A-stable and heavily damps numerical ringing, making it the industry standard for stiff circuits with sharp edges, though it slightly dampens physical oscillations in lossless circuits.

Nonlinear Solvers

At every discrete time step \(t_n\), the discretized DAE becomes a system of nonlinear algebraic equations:

\[F(x_n) = 0\]

To find the circuit state \(x_n\) (voltages and branch currents), PyCircuit uses the Newton-Raphson (NR) method.

  • The solver computes the Jacobian matrix \(J = \partial F / \partial x\).

  • It iteratively updates the state: \(\Delta x = -J^{-1} F(x)\).

  • The NR loop terminates when both the absolute and relative changes in \(x\) fall below user-defined tolerances (vabstol, iabstol, reltol).

Adaptive Time-Stepping Methods

Choosing the right step size \(h = t_n - t_{n-1}\) is critical. Small steps are needed during fast transitions to maintain accuracy, while large steps are needed during inactive periods for simulation speed.

Option B: Standard Predictive Time-Stepping (Legacy)

This is the traditional “trial and error” approach used by SPICE:

  1. The solver predicts a step size \(h\) based on the Local Truncation Error (LTE) of the previous time steps.

  2. It solves the nonlinear NR system for the current state \(x_n\) using \(h\).

  3. After convergence, it calculates the actual LTE for the current step.

  4. If the LTE exceeds the error tolerance (TRTOL), the step is rejected (backed up). \(h\) is reduced, and the solver tries again.

Drawbacks: For highly stiff circuits, the predicted \(h\) is often overly optimistic, leading to frequent backups and wasted matrix inversions.

Option A: Coupled Schur Complement Time-Stepping (New)

This method (based on G. Peter Fang’s “A New Time-Stepping Method for Circuit Simulation”) eliminates rejected time steps entirely by treating the timestep \(h\) as an independent variable to be solved simultaneously with the circuit state \(x\).

The solver couples the LTE error equation \(E(x_n, h) = 0\) directly into the DAE system, forming an augmented \((N+1) \times (N+1)\) system:

  1. Approximate Newton Update: To avoid solving a massive \((N+1)\) matrix, PyCircuit exploits the analytical relationship between LTE and \(h\) (where LTE \(\propto h^3\) for Trapezoidal).

  2. Analytical Gradient (\(E_h\)): The sensitivity of the error to the timestep is calculated analytically as \(E_h = \frac{p \cdot (E + TRTOL)}{h}\), which provides a perfectly smooth, mathematically stable Newton update for \(h\).

  3. Golden Window (\(\gamma\) bounds): Following Section 3.3 of Fang’s paper, the solver defines an acceptable error window (\(0.7\tau \le \epsilon \le 3.0\tau\)). If the current error falls inside this window, the solver accepts the step without trying to over-optimize \(h\), saving thousands of redundant iterations.

This approach guarantees that the solver scales \(h\) aggressively upward when the circuit is quiet, and shrinks \(h\) stably without getting trapped by numerical noise during severe stiffness.