Second-order time integration, and the trap underneath it

Nearly every transient magnetic FEA solver integrates in time with backward Euler. It is unconditionally stable, it is trivial to implement, and it is first-order accurate — halving the time step only halves the error. Its stability also comes from numerical damping, which in an eddy-current problem is indistinguishable from physical damping you did not ask for.

The scheme

The transient field equation is σ·∂A/∂t + curl(ν(B) curl A) = J. Write it as a θ-scheme, so backward Euler stays available as θ = 1 and Crank–Nicolson is θ = ½:

M (Aₙ₊₁ − Aₙ)/dt + θ[G(Aₙ₊₁) − Fₙ₊₁] + (1−θ)[G(Aₙ) − Fₙ] = 0

where G(A) = K_secant(A)·A is the internal magnetic force. Keeping θ general rather than hard-coding ½ is worth the small extra effort: at θ = 1 every formula must collapse exactly to the existing backward-Euler code, which turns the entire existing regression suite into a test of the new one. Our own gate was bit-identical solution files at θ = 1, and it caught a summation-order change on the way through.

One implementation note that is not obvious: store the history force Gₙ captured at the end of step n, rather than re-assembling Kₙ·Aₙ at the start of step n+1. With motion in the model the sliding band re-stitches every step, so the matrix at the old rotor position is no longer reconstructible once the rotor has moved.

The trap: a FEM magnetics system is a DAE, not an ODE

Here is the part that costs people weeks. Nodes inside a conductor carry eddy-current mass (σ > 0) and are genuinely differential. Nodes in air carry no mass at all: their equations are algebraic constraints, to be satisfied at the new time level, not averaged across the step. The system is a differential-algebraic system, and a θ-scheme applied uniformly to all of it quietly loses its order.

We measured this the slow way. Weighting θ per element — the obvious first fix — still came out at about first order on a smooth test case. The reason is the conductor surface: those rows carry mass but assemble stiffness from both conductor elements (θ-weighted) and air elements (fully implicit) in the same equation, which is O(dt) right at the boundary, and that poisons the whole solution. The correct treatment is per row: a node with eddy mass is differential and gets the trapezoidal rule; a mass-free row is algebraic and is enforced at tₙ₊₁.

And the naive way to do that — scaling individual rows — destroys the symmetry of the matrix, at which point CG and LDLT explode, as we observed. The form that works divides each differential row by θ, which leaves the matrix K + M/(θ·dt): symmetric, and the same assembly as backward Euler with an effective step of θ·dt. The θ-dependence moves entirely into the right-hand side, on the masked rows only.

What it is worth

Copper cylinder in a 250 Hz AC field, 10k nodes, σ = 5.96e7, 2 ms of simulated time. Measured convergence order p from a step-refinement study:

Source	Backward Euler	Crank–Nicolson
Smooth start, consistent initial condition	p = 1.00	p = 1.99–2.06
Current jump at t = 0	p ≈ 1.1–1.2	p ≈ 1.1–1.2

On the smooth problem, second order is confirmed in A at interior, surface and air probes alike, with errors 100 to 500 times below backward Euler at the same time step. Put the other way round: to match a Crank–Nicolson result, backward Euler needs somewhere between 100 and 500 times as many steps.

The second row is just as important to state. With a current that jumps at t = 0, both schemes drop to about first order. That is the classic order reduction of a suddenly-applied source — a property of the problem, not of the code — and Crank–Nicolson also rings on such a start. The standard remedy is Rannacher startup: take the first couple of steps with backward Euler to damp the transient, then switch. Nabla does this by default, with the number of startup steps exposed as a setting.

One honest limitation: the eddy-current density reported at tₙ₊₁ is still computed from a backward difference, which is midpoint-accurate under Crank–Nicolson and so stays first order at the grid times, in both schemes. The field is second order; that one derived output is not, and it is documented rather than papered over.

Two bugs worth naming

Both were found by validation, not by review, and both are the kind that produce a plausible answer rather than a crash.

• Frozen permeability. The linear path assembles K once with μ frozen. Re-assembling it each step under Crank–Nicolson picked up a small per-step μ drift on BH-labelled regions, so the two schemes diverged for reasons that had nothing to do with time integration.
• A double-counted back-EMF in the trapezoidal circuit companion model, which made rotor-bar currents jump by a factor of 2.1 on the first Crank–Nicolson step. The check that caught it was comparing against a scalar RL circuit integrated with the trapezoidal rule by hand — a two-line reference model.

When to use which

Use Crank–Nicolson when the time accuracy is the answer: eddy-current diffusion, loss in a solid conductor, an AC transient, a switching event you want resolved rather than averaged. Stay on backward Euler when you are marching to a steady state and only the endpoint matters, or when the source is genuinely discontinuous and you would spend the accuracy on startup ringing anyway. In Nabla it is a per-model setting, with backward Euler as the default, so switching costs one dropdown and a re-run.

Crank–Nicolson time integration is part of Nabla Core. Related: second-order elements in space and the validation dossier.