MNA Theta Method Capacitor Discretization

This post is on the technical side and has a good bit of math in it. In fact, it will not make much sense to you unless you happen to have done some background study on circuit simulation in the modified nodal analysis (MNA) tradition. But for those that are on that path, it may be helpful (it would have been for me) and so I’m sharing for that reason.

To provide a little background, as I began studying analog modeling using MNA techniques, I of course came across Andy Simper’s well known presentation from ADC20 From Circuit to Code: Under the Hood of Analog Modelling.

In his discussion of capacitor discretization formulas, he mentioned the idea that the Forward Euler (FE), Backward Euler (BE) and trapezoidal rule (TR) numerical integration methods can be combined using “mixing”. In his formulation, this is controlled by a mixing variable $m$. At $m=\frac{1}{2}$, you get TR, $m=1$ is BE, and as $m \to 0$, you approach the FE formula (if we actually set $m$ to zero we would get a divide-by-zero in the update equation).

If you happen to be doing something very similar, you of course can use Andy’s equations directly. But as I have been developing a more general simulation code base, and since I was interested in experimenting with this “mixing” formulation, I had a need to understand where it came from so that I could put it in the context of the more general MNA formulation I was working in.

I had bought the book Circuit Simulation by Farid N. Najm, and was working through it to develop my understanding of MNA. It has good coverage of all three basic numerical integration methods (FE, BE, and TR), but makes no mention of this “mixing” technique. So I searched high and low for a reference and did eventually find one in the 2012 paper “Parallel Spectral Numerical Methods” which was written by Gong Chen, Brandon Cloutier, Ning Li, Benson K. Muite and Paul Rigge (henceforth “Chen et al”), where it is referred to as the Theta Method. In their formulation, we get BE at $\theta=0$ and approximate FE as $\theta\to1$ i.e. it is effectively a flipped version of Andy’s formulation. So let us note here that $\theta$ (or Andy’s $m$) is arbitrary in definition and we are free to select either.

I later picked up a copy of the book Fundamentals of Computer-Aided Circuit Simulation by William J. McCalla. This book is frequently referenced by the Najm book, and covers some of the same material, but in places feels more like a survey of the methods in use at the time (1987), unlike Najm which generally feels like more of rigorous introduction to the topic (the Najm book is also much more recent and was first published in 2010).

The McCalla book actually does mention the Theta Method, but only very briefly, without explaining it in any depth or deriving the element companion model equations. In the formulation that is discussed, the controlling variable is $\mu$ and it is defined similarly to Chen at al’s definition of $\theta$ i.e. BE is given with $\mu=0$. McCalla cites W. Liniger and R.A. Willoughby’s 1970 paper “Efficient Integration Methods for Stiff Systems of Ordinary Differential Equations”. I am left to assume what is now known as the Theta Method may therefore have originated with Liniger and Willoughby, but I must admit I haven’t read the paper.

What follows here is a derivation of the MNA stamping equations for both linear and non-linear capacitors using the Theta Method. To this end I first review the patterns established by Najm in the derivation of the BE and TR equations before applying the same approach to the Theta Method. After that I also derive PLTE estimates (again following Najm’s lead) and finally re-derive the minimal (and efficient) update equations that Andy shows in his presentation.

Please note that, much like how the definition of $\theta$ can vary arbitrarily from one reference to the next, companion model conventions do also. With capacitor companion models this comes up around the polarity of the current source. Indeed, even within the Najm book, the convention used for the BE and TR derivations differs in a way that I found confusing. To minimize this kind of confusion, I will start by defining a single convention that will be applied for all examples, and when presenting Najm’s derivations I make adjustments as needed to match this convention. As such my equations do not match Najm’s exactly in some places. In general I have made an effort to point out where this is the case.

As is probably obvious, I have borrowed very heavily from Najm throughout, and so I think it’s appropriate that I direct you to buy his book if you are interested in learning more about any of these techniques. I must say it is certainly one of the more expensive books I’ve bought, and I do have an occasional quibble prompted by my own confusion, but that is just as likely because my linear algebra background is pretty thin as it is because of how the material is presented (not to mention that I’m surely just not as smart as Najm). Overall the book is of high quality and delivers the material in a thoughtful, organized way, so I’d absolutely recommend it.

Lastly let me just say I do not claim to be an expert on any of this, and indeed I am still working to deepen my understanding of many aspects of MNA. Without following Najm’s lead I doubt I’d have been able to derive anything on my own. Just sharing as I learn is all. If you notice I’ve made a mistake somewhere, please let me know.

So without any further rambling, here is the companion model we will use for all derivations below:

And therefore in all cases we will use the following element stamp:

Likewise will express the linearization at time $t_{n+1}$ in the following form:

Where implicitly $I_{eq} \equiv I_{eq,n+1}$ and (for non-linear capacitors) $G_{eq} \equiv G_{eq,n+1}$.

Linear Capacitor

In the case of a linear capacitor, we will use the usual capacitance-based model, where we are given:

Backward Euler

Given the Backward Euler discretization formula:

And given the initial value problem (IVP) specified by $x'(t) = f(x, t)$ and $x(t_0)=x_0$, we have:

Applying this to the capacitor state:

Leads to the following direct discretization:

And the companion model equations:

Trapezoidal Rule

Given the Trapezoidal Rule discretization formula:

And given the IVP specified by $x'(t) = f(x, t)$ and $x(t_0)=x_0$, we have:

Applying this to capacitor state:

Leads to the following direct discretization, where we’ve also made use of $i_n \approx Cu'(t_n)$:

And the companion model equations:

Theta Method

Given the Theta Method (θM) discretization formula:

Note that we have flipped the definition of $\theta$ in comparison with Chen et al by making the substitution $\theta\mapsto(1-\theta)$. This gives us BE at $\theta=1$ which is consistent with Andy Simper’s definition of $m$.

Given the IVP specified by $x'(t) = f(x, t)$ and $x(t_0)=x_0$, we have:

Applying this to the capacitor state:

Leads to the following direct discretization, where we’ve also made use of $i_n \approx Cu'(t_n)$:

And the companion model equations:

Note that with $\theta=\frac{1}{2}$ we recover the TR equations, and with $\theta=1$ we recover the BE equations.

Non-Linear Capacitor

In the case of a non-linear capacitor, we will use the charge-based model (as presented by Najm) rather than the capacitance-based model. With the charge based model, we are given:

Backward Euler

Recall the BE finite difference equation used previously:

Applied to the capacitor state:

Leads to the following direct discretization, where we’ve also made use of $i(t_n) \approx i_n$ and $q(t_n) \approx q(u_n)$:

Recall that the first order Taylor expansion is written like this:

Linearizing $i$ in the Newton loop at $t_{n+1}$ (we are solving for unknown $u \equiv u_{n+1}$):

Taking the derivative of $g_C(u)$ and evaluating at $u_{n+1}^{(k)}$:

Leads to the companion model equations:

The above form is derived identically for the TR and θM integration methods. Only the derivation of $g_C(u)$ and $g_C'(u)$ differs.

Trapezoidal Rule

Recall the TR finite difference equation used previously:

Applied to the capacitor state:

Leads to the following direct discretization, where we’ve also made use of $i(t_n) \approx i_n$, $q(t_n) \approx q(u_n)$, and $i_n \approx q'(t_n)$:

Taking the derivative of $g_C(u)$ and evaluating at $u_{n+1}^{(k)}$:

Where, having derived the above, we can reuse the companion model equations from the BE method.

Theta Method

Recall the θM finite difference equation used previously:

Applied to the capacitor state:

Leads to the following direct discretization, where we’ve also made use of $i(t_n) \approx i_n$, $q(t_n) \approx q(u_n)$, and $i_n \approx q'(t_n)$:

Taking the derivative of $g_C(u)$ and evaluating at $u_{n+1}^{(k)}$:

Where, having derived the above, we can reuse the companion model equations from the BE method.

PLTE Computation For Theta Method

Given local truncation error for θM as:

Where, once again, we have flipped Chen et al’s definition of $\theta$ via substitution $\theta\mapsto(1-\theta)$.

We can derive the order $p$ and error constants $C_0 \ldots C_{p+1}$ as follows:

Equidistant Data (Fixed Time Step)

Given:

For equidistant data, we can estimate PLTE as follows:

Applying this to the Theta Method, we get:

Note that for $\theta = \frac{1}{2}$ we simply have the TR case, and that when $\theta = 1$ we recover the PLTE estimate for the BE case. However, note also that as $\theta \to \frac{1}{2}$, our PLTE estimate becomes vanishingly small before “switching“ to the TR form. This occurs simply because we generally ignore the 2nd order error term when using a 1st order integration method such as BE (or θM when not in the vicinity of $\theta=\frac{1}{2}$. However, it would likely only cause issues for us if we are implementing a variable $\theta$ algorithm, which I currently am not doing!

Non-Equidistant Data (Variable Time Step)

With $p=1$, we are given:

For $p=2$, we are given:

To apply for the Theta Method, we need only plug in the correct error constants.

For $\theta\neq\frac{1}{2}$ we use the $p=1$ form:

Note that with $\theta=1$ we recover the BE result.

And for $\theta=\frac{1}{2}$ we use the $p=2$ form:

Note that this is simply the TR result.

Deriving Andy’s Update Equations

In his presentation, Andy wrote:

I have adjusted his notation to match ours, except that I’ve used his $m$ instead of our $\theta$ to avoid confusion with our own equations. I have also made the $I_{eq}$ index explicit since we will be expressing $I_{eq,n+1}$ in terms of $I_{eq,n}$. Finally, notice that Andy’s companion model current source is defined with opposite polarity from ours, and as such the $I_{eq}$ expression would have to be negated to match ours exactly (I have not gone this far).

Andy goes on to derive the following:

This is nice for a few reasons. First, it is a little more efficient to compute, both because $m$ appears only once in this formulation, but also because the results from evaluating the companion model equations are reused. And thus, it also only requires us to retain a single value from the previous time step: $I_{eq}$.

However, I will note that this more minimal formulation will not work if we are doing a separate DC Analysis step, since we generally treat capacitors as an open circuit in the DC Analysis. For this reason, we need separate voltage and current history terms if we are doing DC analysis prior to proceeding with Transient Analysis.

Now let’s derive Andy’s minimal $I_{eq}$ update equation, starting with our companion model equations and KCL. Note that I’ve once again made the $I_{eq}$ index explicit here.

We proceed with the following algebraic manipulations:

Where, finally, we note that this is exactly the result we’d have gotten by simply applying a sign flip to each of the $I_{eq}$ variables in Andy’s equation.

References

Andrew Simper. (2020). “From Circuit to Code: Under the Hood of Analog Modelling”.
https://www.youtube.com/watch?v=eGcqomH6aAc

Farid N. Najm. (2010). Circuit Simulation.

William J. McCalla. (1987). Fundamentals of Computer-Aided Circuit Simulation.

Gong Chen, Brandon Cloutier, Ning Li, Benson K. Muite and Paul Rigge. (2012).
“Parallel Spectral Numerical Methods”.

Stefan Jahn, Michael Margraf, Vincent Habchi, Raimund Jacob. (2003).
QUCS Technical Papers and in particular the subsection “Energy Storage Components”
https://qucs.sourceforge.net/tech/node26.html