The phase the rule loses

The ringing that belongs to the rule found the trapezoidal rule’s weakness on a pole it cannot resolve: a factor per step that approaches −1, so a fast pole’s error is flipped rather than killed and rings at half the stepping rate. Backward Euler, whose factor goes to zero, has no such ringing. On that evidence backward Euler is the safer rule.

The same property, pointed at a different pole, reverses the verdict. One step, computed twice said the trapezoidal rule preserves the amplitude of a lossless oscillation exactly, at the cost of slightly mis-stating its frequency, and that backward Euler decays a lossless resonance visibly faster than the real one with nothing in the picture to say so. Both statements were qualitative. This essay marches a resonance that should never decay and measures what each rule does to it — in degrees and in nepers, per cycle, as a function of how many steps a cycle the march takes.

The resonance being marched is the one the frequency field measures in closed form. Resonance, and the bandwidth it sets exactly found its half-power width $f_0/Q$ to every digit, and the Q the components allow found the quality factor capped by the parts’ own losses; a backward-Euler march adds a loss no part has, and caps the simulated quality factor the same way. Where the behaviour is written down puts this circuit’s two poles exactly on the imaginary axis — the one place a rule’s factor per step has to have magnitude one to be right.

A resonance with nothing to lose #

The circuit is an inductor of ten millihenries and a capacitor chosen to resonate with it at one kilohertz, driven by a one-volt step with no resistance anywhere in the loop. The capacitor’s voltage is exactly $1 - \cos\omega t$: it rises to two volts, falls to zero, and does that for ever, because nothing in the circuit dissipates. Any decay a march shows is the march’s, and any drift in phase is the march’s.

Each rule advances the circuit by a fixed linear map per step, so the capacitor’s distance from one volt obeys a two-term recurrence, $e_{k+1} = p\,e_k - q\,e_{k-1}$, with $q$ the square of the magnitude of the rule’s factor per step and $p$ twice that magnitude times the cosine of its angle. The figures fit $p$ and $q$ by least squares over the whole marched record, which is exact for a recurrence of that form and needs no zero crossings or fitted envelopes, and read off each rule’s rotation and decay per step. They are then checked against closed forms to a part in a million.

Ten cycles, marched #

At twenty steps a cycle the two marches look nothing alike. The trapezoidal one is a sinusoid of the right height, swinging from zero to two volts every cycle, slightly slow: it loses 2.918 degrees a cycle, and after ten cycles it is 29.2 degrees behind the circuit. Its amplitude per cycle is 1.00000 — the figure checks it to a part in a billion. The backward-Euler march collapses. It keeps 0.3901 of its swing each cycle, and after ten cycles 0.0082 per cent of the ringing is left: the march has settled at one volt, in a circuit with no resistance, as though the resonance had a quality factor of a little over three.

Both are wrong. Only one of them looks wrong. A designer who sees the trapezoidal waveform would accept it; one who sees the backward-Euler waveform might well accept that too, if the circuit were believed to have some loss, and the invented damping would then be read as a property of the components.

What each rule does to a pole on the axis #

A lossless resonance is a pair of poles at $\pm j\omega$. With a step $h$, the trapezoidal rule’s factor per step for such a pole is $(1 + j\omega h/2)/(1 - j\omega h/2)$, whose magnitude is exactly one and whose angle is $2\arctan(\omega h/2)$. Backward Euler’s is $1/(1 - j\omega h)$: magnitude $1/\sqrt{1 + (\omega h)^2}$, angle $\arctan(\omega h)$.

So the trapezoidal rule rotates by slightly less than $\omega h$ each step and never shrinks: over $N$ steps a cycle, $\omega h = 2\pi/N$, it falls behind by

$$2\pi - 2N\arctan\frac{\pi}{N}$$

radians a cycle, which is 2.918° at twenty steps. Backward Euler both rotates by less — $2\pi - N\arctan(2\pi/N)$ a cycle, 11.19° at twenty steps — and shrinks by $(1 + (2\pi/N)^2)^{-N/2}$ a cycle, 0.3901 at twenty. The figure’s measured values agree with all three to a part in a million.

The two rules are the same kind of approximation to $e^{j\omega h}$ — both rational, both of low order — and they spend their error differently. The trapezoidal rule’s factor lies exactly on the unit circle, where the true factor lies, and is wrong only in angle. Backward Euler’s lies inside the circle and is wrong in both. That one geometric fact is the whole difference between a march that keeps a resonance ringing slightly slow and one that kills it.

At ten steps a cycle — coarse, but not unusual for a long simulation of a slow circuit that happens to contain a resonance — the trapezoidal march loses 11.2 degrees a cycle and is 112 degrees behind after ten cycles, a third of a cycle, with its amplitude still exactly right. Backward Euler keeps 0.1894 a cycle, so after ten cycles the ringing is down by a factor of sixteen million.

At fifty steps a cycle the trapezoidal march is 4.7 degrees behind after ten cycles and backward Euler has 2 per cent of its ringing left. At a hundred, 1.2 degrees and 14 per cent. The trapezoidal rule’s error has fallen by a factor of twenty-five for five times the steps, between twenty and a hundred; backward Euler’s loss per cycle, in nepers, has fallen by a factor of 4.8.

Two errors, two powers #

Swept from eight steps a cycle to 128, the two errors separate onto different slopes. The trapezoidal rule’s phase lost a cycle falls as the square of the steps a cycle — the figure fits −1.99 and checks it within a tenth — because $2\arctan(\pi/N)$ differs from $2\pi/N$ by a cubic term, and $N$ of them make a square. Backward Euler’s amplitude lost a cycle, measured in nepers, falls only as the first power — fitted −0.97 — because $\tfrac{N}{2}\ln(1 + (2\pi/N)^2)$ tends to $2\pi^2/N$. Its phase error falls as the square, like the trapezoidal rule’s, and four times larger.

The powers turn into step counts that are not close. To stay within one per cent of a cycle in phase over a hundred cycles, the trapezoidal rule needs 182 steps a cycle. To keep 99 per cent of the amplitude over the same hundred cycles, backward Euler needs 196,404. The first-order loss is not a detail a finer step removes in practice; it is a loss that needs a thousand times the work to push below one per cent.

The same warping a digital filter pre-corrects #

The trapezoidal rule’s slow frequency has a closed form that deserves a name, because it has one elsewhere. A pole at $j\omega$ marched with a step $h$ oscillates at

$$\omega_{\text{num}} = \frac{2}{h}\arctan\frac{\omega h}{2},$$

and that is exactly the frequency warping of the bilinear transform, which is the trapezoidal rule applied to a filter’s differential equation to turn it into a digital one. The corner that moved measured the same compression from the digital side: a corner designed in continuous time and translated by the bilinear transform lands below where it was put, by the arctangent above, and the standard remedy is to pre-warp the design so that the corner lands where it was meant to.

A circuit simulator does not pre-warp. It cannot, because it does not know which frequencies the circuit will ring at before it marches it, and the warping is a function of each resonance’s own $\omega h$. So a simulated resonance at one kilohertz, stepped at twenty steps a cycle, rings at 991.9 hertz — slow by 2.918 degrees a cycle, which is 0.81 per cent — and one at a tenth of a kilohertz, stepped at the same absolute step, rings at 99.99 hertz. The error is not a fixed fraction of the frequency; it is a fraction that grows as the square of the frequency against the step, and it is largest on exactly the resonances a simulation resolves least well.

Backward Euler has a warping of its own, $\arctan(\omega h)/h$, four times as large at small $\omega h$, and a decay on top of it. No digital filter is designed with backward Euler for that reason, and the same reason argues against using it to simulate anything that rings.

A simulation that runs for a second #

The step counts above are per cycle, and a simulation is set by a duration. Put a realistic resonance into a realistic run and the two rules’ requirements come apart by orders of magnitude.

Take a converter’s output filter resonating at ten kilohertz, simulated for one second to watch a slow control loop settle — ten thousand cycles of the resonance. For the ringing’s timing to stay within a tenth of a cycle at the end, the trapezoidal rule’s phase error of about $(2\pi)^3/(12N^2)$ radians a cycle, accumulated over ten thousand cycles, must stay below $0.2\pi$: that needs about 570 steps a cycle, 5.7 million steps for the second. For backward Euler to keep ninety-nine per cent of the ringing’s amplitude over the same ten thousand cycles, its loss of $2\pi^2/N$ nepers a cycle must total under a hundredth: that needs about twenty million steps a cycle, two hundred billion for the second.

Neither simulation is run that way. The trapezoidal one is usually run coarser, and its resonance drifts in phase by a fraction of a cycle that nobody checks because the waveform still looks like a resonance. The backward-Euler one, when a simulator falls back to it, is run at the same coarse step and simply shows no resonance — the filter appears critically damped, and the control loop’s settling is simulated against a plant that does not exist.

Which error a simulation can live with #

A lost degree of phase and a lost per cent of amplitude are not the same kind of error, and the comparison depends on what the simulation is for.

For anything that measures a decay, the trapezoidal rule is the only choice. A quality factor read off a marched ring-down, a settling time on a lightly damped loop, the tail of a resonance excited by a switching edge — each of these is an amplitude against time, and backward Euler adds a decay rate of $2\pi^2/N$ nepers a cycle to whatever the circuit has. At a hundred steps a cycle that is 0.197 nepers a cycle, equivalent to a quality factor of about sixteen in parallel with the circuit’s own; a resonance whose true quality factor is fifty would be read at about twelve.

For anything that measures a time, both rules carry a phase error, and the trapezoidal rule’s is a quarter of backward Euler’s. A frequency read off zero crossings, a delay through a resonant network, the timing of a switching event triggered by a ringing node: all of these accumulate the phase lost per cycle, and after enough cycles they are wrong by a whole cycle. At twenty steps a cycle the trapezoidal rule loses a whole cycle every 123 cycles.

And a long simulation of a lightly damped resonance should be judged by its phase budget. The trapezoidal rule’s phase error per cycle is about $(2\pi)^3/(12N^2)$ radians, so the number of cycles before it has lost a given fraction of a cycle grows as $N^2$; doubling the steps a cycle quadruples the cycles a simulation can run before its timing is wrong. That is the number to set a step by when a resonance is being simulated for many periods, and it is not the number an error control based on local truncation will choose.

Telling the two errors apart on a real run #

The rule of thumb for false ringing — shorten the step and see whether it changes — has a counterpart here, and the powers make it quantitative. Run the same simulation at a step and at half that step, and compare the resonance’s phase at the end of the run. If the march is trapezoidal, its phase lag should shrink by a factor of four, because the lag per cycle goes as the inverse square of the steps a cycle; if what remains after halving is a quarter of what was there, the lag was the rule’s, and a further halving will leave a sixteenth.

The amplitude says the same about backward Euler with a different number. A decay rate invented by the rule halves when the step halves, because it goes as $2\pi^2/N$; a decay rate that belongs to the circuit does not move. So a ring-down whose decay rate falls by half when the step is halved was mostly the rule’s, and one whose decay rate is unchanged was the circuit’s. Two runs, one comparison each, and the measurement says which part of a simulated resonance is a simulation.

Neither test is free, and both are cheaper than being wrong: halving the step doubles the work once, where an unrecognised invented damping can put a control loop’s design against a plant that does not exist.

Where the two essays meet #

The two measurements on this subject point in opposite directions, and they are the same fact about two factors.

On a real pole far faster than the step, the trapezoidal factor tends to −1 and backward Euler’s to zero, and backward Euler is right: the pole should be dead within a step, and the trapezoidal rule rings instead. On an imaginary pole, the trapezoidal factor has magnitude one and backward Euler’s is less, and the trapezoidal rule is right: the pole should never decay, and backward Euler kills it. A rule whose factor sits on the unit circle along the whole imaginary axis cannot also send the far left of the plane to zero, because a rational function that does both would have to change its magnitude discontinuously at infinity.

On the real axis the ordering is the reverse of the imaginary axis’s. The figure is the earlier measurement’s: the trapezoidal factor tends to −1 for a pole far faster than the step and backward Euler’s to zero, so on a pole that should die the trapezoidal rule is the one that is wrong. Set beside the lossless resonance, the two figures are one statement about two rational approximations to an exponential, read along the two axes of the plane where they disagree most.

That is why simulators mix them. The ringing that belongs to the rule took one or two backward-Euler steps after a discontinuity to kill the fast pole’s false ringing and then returned to the trapezoidal rule; the measurement here says what those steps cost a resonance that is ringing at the time. At twenty steps a cycle one backward-Euler step shrinks a resonance by $(1 + (2\pi/20)^2)^{-1/2}$, 0.954 — a loss of four and a half per cent of amplitude in one step, applied every time a discontinuity occurs. A converter model with a switching event every cycle and a resonant tank ringing through them would lose that fraction of its tank amplitude per switching event to the remedy alone.

Still open: the second-order rules that damp, the step error control chooses, and a resonance with its own loss #

A rule between the two. The second-order backward differentiation formula has a factor that goes to zero for fast poles, like backward Euler, and loses amplitude on the imaginary axis only at second order. Measuring its phase and amplitude errors per cycle beside these two would say how much of the trapezoidal rule’s accuracy on a resonance it keeps while losing the false ringing, which is the trade a circuit simulator makes when it offers it.

Error control on a resonance. A variable-step march sets its step by estimating the local error from two rules at each step. On a lossless resonance the local error per step is tiny — third order in the step for the trapezoidal rule — and the accumulated phase error is not local at all. Whether error control lets a long ringing simulation drift a cycle behind while reporting every step within tolerance is the practical form of the phase budget above.

A resonance with loss of its own. With a resistor in the loop the true decay competes with backward Euler’s invented one, and the trapezoidal rule’s factor moves off the unit circle by exactly the right amount at first order. Two numbers without solving for the waveform marched such a circuit for sixty cycles at two hundred steps a cycle and found the phase error growing as the damping fell; putting a number on how a rule’s error in decay rate compares with a circuit’s true decay rate would say at what quality factor backward Euler’s damping dominates the circuit’s.