The rule between the two

One step, computed twice set a step response from the poles beside the same circuit marched forward in time by the trapezoidal rule, and found the difference between them to be the rule’s own error, falling fourfold every time the step was halved. The next two essays found the rule’s two characteristic failures. The ringing that belongs to the rule marched a stiff circuit — a one-microsecond pole driving a millisecond one — at twenty microseconds a step, and the fast node’s error was multiplied by −0.818 every step: an oscillation at half the stepping rate that belonged to the rule, taking twenty-three steps to fall below a per cent. The phase the rule loses marched an undamped inductor and capacitor and found the rule keeping the amplitude exactly and falling behind in phase, 29.2° after ten cycles at twenty steps a cycle, while backward Euler, the rule that cured the stiff ringing, destroyed the amplitude instead.

That essay named the rule a circuit simulator offers as the compromise: the second-order backward differentiation formula, which damps a fast pole the way backward Euler does and loses amplitude on a resonance only at second order. Whether it keeps enough of the trapezoidal rule’s accuracy on a resonance to be worth its price is a measurement, and this page makes it.

The two failures, for reference #

The two rules fail in complementary ways because of where their factors per step sit. For a pole at $\lambda$ and a step $h$, each rule multiplies the solution by a factor $z(\lambda h)$ every step. The trapezoidal rule’s is $(1 + \lambda h/2)/(1 - \lambda h/2)$, which has magnitude exactly one on the imaginary axis — an undamped ringing keeps its amplitude — and tends to −1 as $\lambda h \to -\infty$, so a very fast pole’s error is kept and flipped. Backward Euler’s is $1/(1 - \lambda h)$, which is less than one everywhere off the origin: a fast pole’s error vanishes, and so does an undamped ringing.

A rule that looks back two steps #

The second-order backward differentiation formula, BDF2, approximates the derivative at the new point from three points rather than two:

$$\frac{3x_{n+1} - 4x_n + x_{n-1}}{2h} = f(x_{n+1}).$$

Like backward Euler it evaluates the circuit only at the new point, which is what gives it backward Euler’s strength against stiffness; unlike it, the difference is accurate to second order, which is what the trapezoidal rule has. It needs two previous values, so its first step is taken by backward Euler, the usual start.

Applied to $x' = \lambda x$ it gives $(3 - 2\lambda h)z^2 - 4z + 1 = 0$, a quadratic with two roots. One, the principal root, tends to one as the step shrinks and carries the solution. The other, the parasitic root, is about a third and decays within a few steps. On the imaginary axis, $\lambda h = j\theta$, the principal root is

$$z = \frac{2 + \sqrt{1 + 2j\theta}}{3 - 2j\theta},$$

whose magnitude falls short of one only at fourth order in $\theta$ and whose angle lags $\theta$ at third order, with four times the trapezoidal rule’s coefficient.

The parasitic root deserves a sentence, because it is what makes BDF2 a two-step rule and not a one-step rule with a different factor. Every step the solution is a combination of the principal root, raised to the step number, and the parasitic one. The parasitic root starts at a third for a slow pole and shrinks further for a fast one, so its contribution is divided by at least three every step and is below a part in a million within thirteen. It is excited by the start — by the backward-Euler step, and by any discontinuity — and it is the reason a two-step rule needs a one-step rule to start it. After that it is invisible, and every error below is the principal root’s.

Ten cycles with three rules #

At twenty steps a cycle BDF2 is plainly between the other two. It keeps 0.9616 of the amplitude a cycle where backward Euler keeps 0.3901 and the trapezoidal rule all of it, so after ten cycles 67.6 per cent of the ringing remains, against 0.0082 per cent for backward Euler. And it falls behind in phase by 10.5° a cycle, which is three and a half times the trapezoidal rule’s 2.9° and about the same as backward Euler’s 11.2°.

Measured from the march — by fitting the recurrence the solution satisfies once the parasitic root has died away — the amplitude and phase a cycle agree with the principal root’s magnitude and angle to the digits drawn. The march is doing what its factor says, and the factor is what the rule is.

How each error falls with the step #

Twenty steps a cycle is coarse. The question a simulator’s user cares about is how fast each rule’s errors fall as the step is refined, since that decides how many steps a long simulation needs.

The three rules’ errors fall at three different rates. Backward Euler’s amplitude loss a cycle falls only as the first power of the steps a cycle, which is why it needs 196,404 steps a cycle to keep 99 per cent of a ringing over a hundred cycles. BDF2’s amplitude loss falls as the cube — measured at an exponent of −2.91 over the range drawn, tending to −3 — and so needs only 157. Its phase loss falls as the square, as the trapezoidal rule’s does, with a constant that tends to exactly four times the trapezoidal rule’s; so where the trapezoidal rule stays within one per cent of a cycle over a hundred cycles at 182 steps a cycle, BDF2 needs 363, about twice as many.

The exponents follow from the principal root’s expansion, and the constants can be read off it too. The trapezoidal rule turns by $2\arctan(\theta/2) \approx \theta - \theta^3/12$ a step, so it loses $\theta^3/12$ of phase a step and $N\theta^3/12 = (2\pi)^3/(12N^2)$ a cycle. BDF2 turns by $\theta - \theta^3/3$ a step: the same power with four times the coefficient, which is where the factor of four comes from. Its magnitude is $1 - \theta^4/4$ a step, so it loses $N\theta^4/4 = (2\pi)^4/(4N^3)$ nepers a cycle — one power of $N$ better than its own phase error, and two better than backward Euler’s $2\pi^2/N$. The measured exponents, −2.91 and −1.95, are those powers seen over a range that does not yet reach the asymptote; at 128 steps a cycle the amplitude loss is 1.85 × 10⁻⁴ nepers, which is $(2\pi)^4/(4 \cdot 128^3)$ to three figures.

That is the price, and it is modest. A factor of two in steps buys a rule that is immune to the stiff ringing the next section measures, and the price is paid in phase, which is the error that matters for how long a waveform stays in step with the truth, rather than in amplitude, where BDF2’s error is already smaller than the phase error at any step count where either is small.

The stiff circuit, marched #

The stiff circuit is where the trapezoidal rule’s factor approaching −1 does damage. A step of twenty microseconds is twenty of the fast pole’s time constants, which is a perfectly reasonable step for the slow pole’s millisecond but is wildly too long for the fast one, and the trapezoidal rule then preserves the fast pole’s starting error almost intact and flips it every step.

BDF2 takes three steps to bring the fast node within a per cent, against the trapezoidal rule’s twenty-three and backward Euler’s two. Its error shrinks by 0.1525 a step in magnitude — exactly $1/\sqrt{3 + 2h/\tau}$, the magnitude of its roots for a real pole this fast, which have become a complex pair — and it rotates as it shrinks, so the error changes sign every few steps rather than every one, but it is gone before the rotation matters. The first step, taken by backward Euler, does most of the work.

The first step matters more here than anywhere else. BDF2 needs two previous values, so the step at the discontinuity is taken by backward Euler, and backward Euler’s factor at twenty time constants is 0.0476: the first step alone takes the error from a volt to 48 millivolts. BDF2 then takes over and shrinks what is left by 0.1525 a step, rotating. A BDF2 started instead by a trapezoidal step would begin from the trapezoidal rule’s −0.818 and have that much further to go; a simulator that starts its backward rule with a backward step has chosen the start that suits it.

Why a backward rule cannot ring #

The difference between the rules on a stiff pole has a plain reason. At the instant a step arrives, the fast node’s derivative is enormous — a volt over a microsecond — and a step of twenty microseconds is far too long to follow it. The trapezoidal rule advances by the average of the derivative at the start of the step and at the end, and the start’s derivative is the enormous one. Averaged in, it carries the node a long way past its final value, and the next step’s start is then the overshoot’s large derivative the other way. The error is handed from step to step, flipped each time, and shrinks only by what the average fails to cancel.

A backward rule uses only the derivative at the end of the step — BDF2 with a two-step correction for accuracy — and at the end of a step twenty time constants long, the fast node’s derivative is whatever pulls it towards its final value from where it has landed. There is nothing from the start of the step to carry it past. The rule lands short of the final value, or on it, and never beyond it by the whole amount. That is the entire mechanism, and it is why every backward difference rule, of any order, damps a real pole that is fast against the step: the higher orders give up some of their stability for poles that ring, but none of them carries a fast decaying error from one step into the next.

The factors for a fast pole #

The three curves say, in one picture, what kind of rule each is. The trapezoidal rule’s factor dips to zero at a step of two time constants — the only step at which it damps a fast pole completely — and climbs back towards one in magnitude on either side, so that a step a thousand times the fast pole’s time constant leaves its error at 0.996 of its size a step. Backward Euler’s falls as one over the step. BDF2’s falls as one over the square root of the step: more slowly than backward Euler’s, but to zero, which is the property the numerical-analysis literature calls L-stability and the trapezoidal rule does not have.

So BDF2 is not a halfway house between the other two in every respect. Against a stiff pole it behaves like backward Euler, a little less briskly. Against a resonance it behaves like the trapezoidal rule in its order of accuracy, a little less accurately in phase and much more accurately in amplitude than backward Euler. That combination — damping where the trapezoidal rule rings, accuracy where backward Euler destroys — is why simulators offer it, and the numbers here are the size of the compromise.

What a simulator’s user should take from the numbers #

A circuit with fast poles that the step does not resolve and slow ringing that it does is the ordinary case: a switching converter’s parasitic ringing beside its control loop, a digital edge beside an analogue filter. The trapezoidal rule rings on the first and is excellent on the second. Backward Euler is clean on the first and ruinous on the second. BDF2 is clean on the first and costs about twice the steps of the trapezoidal rule on the second.

The previous essays found the trapezoidal rule’s ringing cured by a few backward-Euler steps at each discontinuity, which is the other common remedy, and a cheaper one if the discontinuities are few. BDF2 is the remedy that needs no knowledge of where the discontinuities are, paid for everywhere in phase. A simulator that switches between them — trapezoidal between discontinuities, backward steps at them — is trying to have both, and each of the numbers above says what it gives up in each regime.

In numbers for a concrete job: a 100 kHz tank simulated for a millisecond is a hundred cycles. To stay within one per cent of a cycle in phase, the trapezoidal rule needs 182 steps a cycle, 18,200 steps in all; BDF2 needs 363, 36,300; backward Euler, which would hold phase only if the amplitude survived, would need close to twenty million to keep the amplitude. If the same simulation also contains a switching edge whose parasitic pole is a nanosecond, every one of those steps is tens of the fast pole’s time constants, and the trapezoidal rule’s 18,200 steps come with a ringing at half the stepping rate after every edge, while BDF2’s 36,300 come clean. Which is cheaper depends on how many edges there are and whether the ringing would be mistaken for the circuit’s own.

How the numbers were obtained #

Every march is the circuit’s nodal equations integrated in time, each reactive element replaced by the companion model its rule implies: for BDF2, the new state is found from $(G + 3C/2h)\,x_{n+1} = (C/2h)(4x_n - x_{n-1}) + b_{n+1}$, after one backward-Euler step. Each rule’s amplitude and phase a cycle are measured from the march itself by fitting, by least squares, the two-term recurrence a single pair of complex roots obeys, over the whole record after the first forty steps, and compared with the closed-form factors. The stiff node’s error factor is found the same way over its first eight steps. The step counts for a hundred cycles are bisected on the closed forms.

What it leaves out #

Higher orders. BDF3 to BDF6 continue the family with more history and higher order, and give up stability progressively for poles near the imaginary axis, which is where a lightly damped resonance lives; how their phase and amplitude errors compare on a resonance is the same measurement with more roots to follow.

Variable steps. A simulator chooses its step from an estimate of the local error, and a variable-step BDF2 is not the fixed-step rule measured here; its behaviour on a long resonance depends on how the step controller reads an error that is mostly phase.

And a resonance with loss of its own. Every resonance here is undamped, so every loss of amplitude is the rule’s. With a real resistor in the circuit the rule’s amplitude error is competing with the circuit’s own decay, and at what quality factor BDF2’s added damping becomes comparable to the circuit’s is the number that decides whether it is noticed.

Still open: the rule’s damping against the circuit’s, the step a controller chooses, and a higher order #

A resonance with its own loss. BDF2’s amplitude loss a cycle at N steps a cycle is known in closed form, and a circuit of quality factor Q loses π/Q nepers a cycle of its own. Setting the two equal gives the quality factor above which the rule’s damping dominates the circuit’s at a given step count, which is the practical limit on simulating a high-Q resonator with a backward rule.

The step a controller chooses. A variable-step BDF2 estimates its local error from the difference between its prediction and its correction. On a lossless resonance that estimate is dominated by the phase error; whether a controller holding a stated tolerance keeps a hundred cycles within one per cent, or drifts a cycle behind while reporting every step within tolerance, is the question the trapezoidal rule’s version of this essay also left.

BDF3. One more step of history raises the order to three and costs some of the stability region. Its amplitude and phase errors a cycle, beside these three, would say whether the next order is worth its history on the resonances circuits actually have.