V19 — Backward Euler vs Crank-Nicolson order study

← Back to the validation dossier

Tier A+S (analytic + solver-internal) - priority P1 Reference frozen on 2026-08-15. Nabla 0.1.0, solver licence mode enforcing, run in 126.8 s.

1. Problem

Nabla integrates every transient with the one-parameter theta scheme, and offers two settings of it: Backward Euler (theta = 1) and Crank-Nicolson (theta = 1/2, a Core Pro capability). The first is documented as first-order accurate and the second as second-order. This chapter measures both claims, and measures what the second one is worth.

Two problems, twenty-four solves.

rl_*	a 200-turn coil on a linear (mu_r = 1000) closed core, fed from a 100 V step through 20 ohm, planar, no eddy currents
slab_*	1-D magnetic diffusion into a 10 mm half-slab of copper (sigma = 5.8e7), driven by a step current sheet across a 5 mm air gap
Levels	16, 32, 64, 128, 256 steps over the same window - 16x in dt
Schemes	Backward Euler and Crank-Nicolson, each at every level
Extra	one 2048-step run of the slab, and one run per problem with Crank-Nicolson forced through its Rannacher startup on every step
Element order	P1 throughout - this is a study of time, and the element order is held fixed

2. The trap this case is built around

A temporal order study is easy to run and easy to get wrong, and the way it goes wrong is always the same: the reference carries a spatial error, the spatial error does not shrink when dt does, the convergence curve flattens onto it, and the fit reports the ruler rather than the scheme. A fit that reads 1.7 where it should read 2 looks like a mild implementation defect and is usually not one.

So the two problems are chosen for two different relationships to that trap.

The circuit has no spatial error at all in the quantity it is judged on. For a coil of any geometry on a linear core, the semi-discrete field-circuit system is K A = C i with lambda = C^T A, so

lambda = (C^T K^-1 C) i = L i        identically

L is a constant of the mesh, not an approximation to one, and the loop equation is exactly L di/dt + R i = V. The whole spatial discretization sits inside one number, which the case measures with a single static solve at 1 A on the same mesh: L = 0.423693 H, tau = L/R = 0.0211846 s. The exponential (V/R)(1 - e^(-t/tau)) is therefore the exact solution of the discretized problem the transient is integrating - refine the mesh and the run and its reference move together - and every deviation from it is time-integration error and nothing else.

The slab is the opposite, deliberately. It is a real field problem, its spatial error is real, and this chapter measures it: at 2048 steps, where the temporal error is negligible, the flux density inside the slab sits 0.00124568 of the surface value away from the exact series. That is the floor, and §5 shows what it does to an order fit that does not know it is there.

3. Reference

Derived in full in the case's reference notes. Four things, and only the first two are about the schemes.

The orders. Expanding the theta scheme about t_(n+1/2) gives a local truncation error O(dt^2) for theta != 1/2 and O(dt^3) at theta = 1/2, hence global orders 1 and 2 (Butcher §3.2). Those integers are the references; the tolerance is the +/-0.15 band the validation plan §3.4 fixes for a fitted order.

Crank-Nicolson is A-stable but not L-stable, so the discontinuous initial condition both problems have - a step response from a field-free state - excites a (-1)^n mode that decays only algebraically (Hairer & Wanner §IV.3). The remedy is Rannacher's: take the leading steps fully implicitly. One such step of size dt costs O(dt^2) globally, which is second order and not below it. Nabla does this by default (cnStartupSteps = 2), and the case measures the scheme as shipped, startup included.

The drive on the slab is Ampere's law, exactly. The far edge of the model is natural, so there is no field beyond the current sheet; a loop from there into the air gap encloses the sheet and nothing else, whatever the slab is doing. So B_surface = mu0 J w for all t > 0, independent of the response. Nabla reads 0.0628318 T against 0.0628319 T, -0.000 %.

The response is Carslaw & Jaeger. With mu constant, B = -dA/dx obeys the diffusion equation with alpha = 1/(mu sigma), so this is the slab quenched at both faces (§3.4 eq. 6) - the same series V08 problem 5 already uses, shared as one function rather than transcribed twice.

4. Two instrument checks, before any slope

The circuit really is that ODE. The same elimination that gives the exponential gives the exact discrete solution of the theta scheme, i_n = [i_(n-1)(L/dt - (1-theta)R) + V]/(L/dt + theta R). Nabla's Backward-Euler run reproduces it to 7.45626e-13 of the steady-state current, worst step of all five levels - and Crank-Nicolson reproduces its own recurrence, Rannacher startup included, to 7.61702e-13. Both are round-off.

This row is not decoration. A coupling that carried an inductance the case did not know about - an end-winding term, a stray branch, a factor on the flux linkage - would shift the effective time constant and appear in the order fit only as a small, plausible, wrong number. It is checked first because after it there is nothing left in the circuit problem except the time integrator.

theta = 1 inside the Crank-Nicolson path is Backward Euler. The validation plan asks for this invariant by name. It is reachable through the shipped API by asking for more Rannacher startup steps than the run has, since theta_eff = 1 while simStep < max(1, cnStartupSteps). Both problems return exactly zero difference: 0 on the circuit and 0 on the field. Bit-identical, which matters because the field path assembles a second, frozen matrix for Crank-Nicolson (K_cn = K_stiff + M/(theta dt)) and is supposed to leave it unused at theta_eff = 1.

5. Result

Every judged row passes. The four orders:

Problem	Scheme	Reference	Fitted	Theory
circuit	Backward Euler	exact exponential	0.959764	1
circuit	Crank-Nicolson	exact exponential	1.9637	2
slab	Backward Euler	2048-step run	0.924403	1
slab	Crank-Nicolson	2048-step run	2.09515	2

One feature of the figures needs saying, because it looks like a mistake and is not: at the coarsest level the Crank-Nicolson error equals the Backward-Euler error, on both problems. The error is a maximum over the shared time stations, and at 16 steps the coarsest station is the single Rannacher startup step, which both schemes take fully implicitly. From the second level on, that station is a Crank-Nicolson step and the curves separate. The startup is the reason Crank-Nicolson converges here at all - without it a trapezoidal rule meets a discontinuous initial condition with a (-1)^n mode it cannot damp - and the price it charges is exactly this one point.

And the trap, measured. The slab's orders above are fitted against the 2048-step run, because there the spatial error is identical in every run and cancels exactly in the difference. Fit the same five solutions against the exact series instead and Backward Euler is unmoved - 0.943266 against 0.924403 - while Crank-Nicolson reads 1.67998, and that number is published without a verdict because it is not about the time integrator at all.

The five-point fit understates how sharply it happens, because it averages the levels that were still converging with the ones that have stopped. The last interval alone is the honest statement, and it costs no extra run:

Finest pair of levels (128 -> 256 steps)	vs the 2048-step run	vs the exact series
Backward Euler	0.997557	1.0408
Crank-Nicolson	1.97001	0.0992163

Same two solutions in both columns. Halving the time step still buys Crank-Nicolson the full second order - and buys it nothing at all when the error is read against the exact PDE solution, because by then it is below the 0.00124568 spatial error of the mesh. Backward Euler, whose error at 256 steps is still an order of magnitude above that floor, does not notice which reference it is measured against.

This is the same warning V09 carries about a mesh convergence study: a convergence study that shows the wrong order may only have found its own ruler. Here it is not a warning, it is a row.

6. What second order is worth

An order is an asymptotic statement. The number a reader can use is the step size it buys at a stated accuracy, so the case solves the fitted err = C dt^p at 1 % error and reports it in units of each problem's own time constant.

Problem	Backward Euler	Crank-Nicolson	Ratio
circuit (dt/tau)	0.0564884	0.176667	3.12749x
slab (dt/t1)	0.0066173	0.0312465	4.72194x

Three things a user should take from that table.

1. The advantage is a ratio, not a factor. It grows as the accuracy demanded tightens - at 10 % error the two schemes are much closer, and at 0.1 % much further apart. Quoting "Crank-Nicolson is twice as good" is meaningless without an accuracy attached.
2. The two problems disagree, and are entitled to: a circuit is one mode and a diffusion equation is a stiff spectrum, so the constant in front of dt^p differs and the crossover moves with it.
3. Both schemes need a small step in absolute terms. 1 % on the diffusion problem costs 0.0312465 of a diffusion time constant even at second order. A second-order scheme is not permission to take a coarse step; it is permission to take a less fine one.

At the finest step of the study the accuracy delivered is 0.00355529 (BE) against 0.000130951 (CN) on the circuit, and 0.00857612 against 0.000320472 on the slab - a factor of 27 in both cases, from the same mesh, the same matrix and the same number of solves.

7. Limitations

• Linear, both problems. The theta scheme's order is a statement about the time discretization; a saturating core adds a Newton loop whose tolerance becomes part of the answer. Nothing here says what Crank-Nicolson does when the nonlinear residual is only driven to 1e-6.
• No motion. A sliding band rebuilds its constraints each step; whether the second order survives that is V20's question, not this one.
• P1 only. The element order is held fixed on purpose - mixing a spatial and a temporal refinement in one fit is how a study loses the ability to say which one it measured.
• The slab's spatial floor is this mesh's, not Nabla's. A finer mesh or P2 would push it down and let the exact series see further; the point of §5 is that a study has to know where its floor is, not that this one is low.

8. All quantities

Quantity	Metric	Nabla	Reference	Error	Tolerance	Verdict
rl_recurrence_be_max_abs	abs	7.45626e-13 -	0 -	7.456e-13	1e-07	PASS
rl_recurrence_cn_max_abs	report	7.61702e-13 -	- -	-	-	REPORT
rl_order_be	slope	0.959764 -	1 -	-0.0402	+/-0.15	PASS
rl_order_cn	slope	1.9637 -	2 -	-0.0363	+/-0.15	PASS
rl_error_be_finest	abs	0.00355529 -	0 -	0.003555	0.01	PASS
rl_error_cn_finest	abs	0.000130951 -	0 -	0.000131	0.001	PASS
rl_theta1_vs_be_max_abs	abs	0 -	0 -	0	1e-12	PASS
rl_L_H	report	0.423693 H	- H	-	-	REPORT
rl_tau_s	report	0.0211846 s	- s	-	-	REPORT
rl_dt_over_tau_1pct_be	report	0.0564884 -	- -	-	-	REPORT
rl_dt_over_tau_1pct_cn	report	0.176667 -	- -	-	-	REPORT
rl_dt_ratio_cn_over_be	report	3.12749 -	- -	-	-	REPORT
slab_surface_B_T	rel	0.0628318 T	0.0628319 T	-0.000 %	1.000 %	PASS
slab_fine_vs_series_max_abs	abs	0.00124568 -	0 -	0.001246	0.01	PASS
slab_order_be	slope	0.924403 -	1 -	-0.0756	+/-0.15	PASS
slab_order_cn	slope	2.09515 -	2 -	0.0951	+/-0.15	PASS
slab_order_be_vs_series	slope	0.943266 -	1 -	-0.0567	+/-0.15	PASS
slab_order_cn_vs_series	report	1.67998 -	- -	-	-	REPORT
slab_last_order_cn	report	1.97001 -	- -	-	-	REPORT
slab_last_order_cn_vs_series	report	0.0992163 -	- -	-	-	REPORT
slab_last_order_be	report	0.997557 -	- -	-	-	REPORT
slab_last_order_be_vs_series	report	1.0408 -	- -	-	-	REPORT
slab_error_be_finest	abs	0.00857612 -	0 -	0.008576	0.02	PASS
slab_error_cn_finest	abs	0.000320472 -	0 -	0.0003205	0.001	PASS
slab_theta1_vs_be_max_abs	abs	0 -	0 -	0	1e-12	PASS
slab_dt_over_t1_1pct_be	report	0.0066173 -	- -	-	-	REPORT
slab_dt_over_t1_1pct_cn	report	0.0312465 -	- -	-	-	REPORT
slab_dt_ratio_cn_over_be	report	4.72194 -	- -	-	-	REPORT
rl_dt_over_tau	report	(series)	-	-	-	REPORT
rl_dt_s	report	(series)	-	-	-	REPORT
rl_err_be	report	(series)	-	-	-	REPORT
rl_err_cn	report	(series)	-	-	-	REPORT
slab_dt_over_t1	report	(series)	-	-	-	REPORT
slab_dt_s	report	(series)	-	-	-	REPORT
slab_err_be	report	(series)	-	-	-	REPORT
slab_err_be_vs_series	report	(series)	-	-	-	REPORT
slab_err_cn	report	(series)	-	-	-	REPORT
slab_err_cn_vs_series	report	(series)	-	-	-	REPORT

• rl_recurrence_cn_max_abs: the same comparison for Crank-Nicolson, including the Rannacher startup steps; published without a verdict because the startup rule the reference assumes is read off the solver source rather than from a specification
• rl_L_H: the finite-element inductance this mesh has, which is half of the ODE both the run and its reference obey
• rl_tau_s: L/R, the circuit's only time scale and the unit every step size below is quoted in
• rl_dt_over_tau_1pct_be: the practical number: the largest step, in circuit time constants, at which Backward Euler holds 1 % error
• rl_dt_over_tau_1pct_cn: the same for Crank-Nicolson
• rl_dt_ratio_cn_over_be: how many times larger a step Crank-Nicolson buys at 1 % error on this problem
• slab_order_cn_vs_series: published without a verdict BY DECLARATION: this fit is expected to read well below 2, and the reason is the case's own spatial error rather than anything the time integrator did
• slab_last_order_cn: the order the FINEST pair of levels alone shows, against the 2048-step run - the reference in which the spatial error cancels
• slab_last_order_cn_vs_series: the same finest pair, measured against the exact series instead - this is the number that shows the floor rather than describing it
• slab_last_order_be: the same finest-pair order for Backward Euler, against the 2048-step run
• slab_last_order_be_vs_series: the same, against the exact series
• slab_dt_over_t1_1pct_be: the largest step, in diffusion time constants, at which Backward Euler holds 1 % error on the field problem
• slab_dt_over_t1_1pct_cn: the same for Crank-Nicolson
• slab_dt_ratio_cn_over_be: how many times larger a step Crank-Nicolson buys at 1 % error on the field problem
• rl_dt_over_tau: extracted but not frozen in expected.json - reported only
• rl_dt_s: extracted but not frozen in expected.json - reported only
• rl_err_be: extracted but not frozen in expected.json - reported only
• rl_err_cn: extracted but not frozen in expected.json - reported only
• slab_dt_over_t1: extracted but not frozen in expected.json - reported only
• slab_dt_s: extracted but not frozen in expected.json - reported only
• slab_err_be: extracted but not frozen in expected.json - reported only
• slab_err_be_vs_series: extracted but not frozen in expected.json - reported only
• slab_err_cn: extracted but not frozen in expected.json - reported only
• slab_err_cn_vs_series: extracted but not frozen in expected.json - reported only

Mesh, per sub-model:

Sub-model	Element order	Nodes	Elements
rl_static	P1	1832	3554
rl_be_L0	P1	1832	3554
rl_be_L1	P1	1832	3554
rl_be_L2	P1	1832	3554
rl_be_L3	P1	1832	3554
rl_be_L4	P1	1832	3554
rl_cn_L0	P1	1832	3554
rl_cn_L1	P1	1832	3554
rl_cn_L2	P1	1832	3554
rl_cn_L3	P1	1832	3554
rl_cn_L4	P1	1832	3554
rl_theta1	P1	1832	3554
slab_fine	P1	1263	2317
slab_be_L0	P1	1263	2317
slab_be_L1	P1	1263	2317
slab_be_L2	P1	1263	2317
slab_be_L3	P1	1263	2317
slab_be_L4	P1	1263	2317
slab_cn_L0	P1	1263	2317
slab_cn_L1	P1	1263	2317
slab_cn_L2	P1	1263	2317
slab_cn_L3	P1	1263	2317
slab_cn_L4	P1	1263	2317
slab_theta1	P1	1263	2317

9. Frozen-reference changelog

No entry: the reference values and tolerances are as first frozen.

10. Reproduce

cd validation
python -m harness.run_case V19

The twenty-four models are published beside this chapter in models/.