import json
import warnings
import numpy as np
import plotly.graph_objects as go
from clerq import PropagationSolver, plot # noqa: F401 (registers the clerq Plotly template)
from clerq.exceptions import PropagationInstabilityErrorNumerical Stability and Step Sizes
PropagationSolver marches the fields through the medium in steps of \(\Delta z\) (z_steps recorded steps, each split into z_steps_inner finite-difference steps), and records them on a time grid of t_steps steps. Nothing in the physics tells you how fine these need to be. Too coarse a \(z\)-step gives a wrong answer that may grow without bound, and too coarse a time grid gives a wrong answer that looks perfectly reasonable.
solve() checks for both. This page shows what the checks catch, how to read them and how to choose the step sizes.
The test problem is a \(2\pi\) hyperbolic-secant pulse in a resonant two-level medium. This is self-induced transparency: the pulse passes through unchanged apart from a delay, so any change in its peak is numerical error.
problem = {
"system": {
"num_states": 2,
"fields": [
{
"coupled_levels": [[0, 1]],
"rabi_freq_t_func": "sech",
"rabi_freq_t_args": {
"ampl": 0.8384014365421667,
"centre": 0.0,
"width": 0.3796628587572578,
},
}
],
},
"t_min": -2.0,
"t_max": 14.0,
"t_steps": 160,
"z_min": 0.0,
"z_max": 1.0,
"z_steps": 10,
"interaction_strengths": [10.0],
}
def make(**changes):
return PropagationSolver(**{**problem, **changes})
def peak_field(solver):
"""Peak |Omega| of the field at each recorded z, as a cyclic frequency."""
return np.abs(solver.rabi_freqs_zt[0]).max(axis=1)Too coarse a z-step
With z_steps_inner=2 the solver runs without complaint, but the pulse should leave the medium with the peak it went in with.
coarse = make(z_steps_inner=2)
coarse.solve(progress=False)
print(f"peak in {peak_field(coarse)[0]:.3f}, out {peak_field(coarse)[-1]:.3f}")peak in 0.838, out 3.926/var/folders/41/hkxt57911pzcfzrv2p66w1z80000gn/T/ipykernel_90155/2736311170.py:2: UserWarning: The z-step may be too coarse: estimated accumulated relative error is 466.0% (tolerance 10.0%). Try z_steps_inner=14 (currently 2) or z_steps_inner='auto'.
coarse.solve(progress=False)
/var/folders/41/hkxt57911pzcfzrv2p66w1z80000gn/T/ipykernel_90155/2736311170.py:2: UserWarning: The time grid may be too coarse: the fields or coherences are interpolated from every other time sample with a relative error of 102.6% (tolerance 10.0%). Try t_steps=1491 (currently 160), and repeat if the warning remains; suggest_t_steps() gives an estimate from the problem's frequency scales.
coarse.solve(progress=False)
/var/folders/41/hkxt57911pzcfzrv2p66w1z80000gn/T/ipykernel_90155/2736311170.py:2: UserWarning: A field is still changing near the edge of the time window (by 41.0% of its peak), so a pulse or its ringing may be cut off. Try widening t_min and t_max.
coarse.solve(progress=False)The pulse that goes in with a peak of 0.84 comes out with 3.9, nearly five times too large, and before the checks nothing said so. The first warning is the one to act on: the \(z\)-step is too coarse. The other two (the time grid and the window edge) are consequences. The wrong field has structure that the time grid cannot resolve and spreads to the edge of the window, so fixing the \(z\)-step usually clears them. Here is the peak along the medium as the \(z\)-step is refined:
fig = go.Figure()
with warnings.catch_warnings():
warnings.simplefilter("ignore")
for n in [1, 2, 4, 16, 64]:
s = make(z_steps_inner=n)
s.solve(progress=False)
fig.add_scatter(
x=s.zlist, y=peak_field(s), mode="lines", name=f"z_steps_inner={n}"
)
fig.update_layout(
xaxis_title="z",
yaxis_title="Peak Rabi frequency",
yaxis_type="log",
width=700,
height=400,
)
fig.show(renderer="notebook_connected")The step error estimate
While solving, clerq estimates the error of the \(z\)-stepping from how much the derivative \(\partial\Omega/\partial z\) changes between steps: the difference between the Euler and Adams-Bashforth predictions. This costs nothing extra. It is summed over the steps and kept relative to the peak field, as solver.step_error_estimate. solve() warns when it exceeds step_error_warn (0.1 by default).
Here it is against the true error, measured from a very finely stepped solution:
reference = make(z_steps_inner=256)
reference.solve(progress=False, check_stability=False)
exact = peak_field(reference)[-1]
print("z_steps_inner estimate true error")
with warnings.catch_warnings():
warnings.simplefilter("ignore")
for n in [1, 2, 4, 8, 16, 32, 64, 128]:
s = make(z_steps_inner=n)
s.solve(progress=False)
true_err = abs(peak_field(s)[-1] - exact) / exact
print(f"{n:13d} {s.step_error_estimate:8.3g} {true_err:10.3g}")z_steps_inner estimate true error
1 2.42 6.19
2 4.66 3.72
4 5.52 0.881
8 3.08 0.168
16 0.484 0.0268
32 0.0837 0.00964
64 0.0207 0.000783
128 0.0052 9.33e-05The estimate is a heuristic, not a stability limit like a Courant number, so read it as follows:
- It is conservative for the default Adams-Bashforth scheme (
step="ab"). In the table, an estimate of 0.08 goes with a true error of 1%, and 0.5 with 3%. Above about 3 the result is wrong by order one. - It falls with
z_steps_innerat the rate of the scheme: as the square for Adams-Bashforth and linearly for Euler. The suggestion in the warning uses this, and it can be a little off the first time. Repeat if the warning remains. - For
step="euler"it is close to the true error rather than conservative (see below), and Euler is much less accurate at the same step size. Prefer the default. - Whether \(\Delta z\) is small enough depends on the density, the pulse and the decay rates together, which is why there is no single safe default.
with warnings.catch_warnings():
warnings.simplefilter("ignore")
s = make(z_steps_inner=128)
s.solve(step="euler", progress=False)
true_err = abs(peak_field(s)[-1] - exact) / exact
print(f"Euler, 128 inner steps: estimate {s.step_error_estimate:.3g}, true error {true_err:.3g}")Euler, 128 inner steps: estimate 0.192, true error 0.142Choosing the step automatically
Pass z_steps_inner="auto" and solve() starts from 2 inner steps and repeats the solve with more until the estimate is under step_error_warn, up to max_z_steps_inner (512).
auto = make(z_steps_inner="auto")
auto.solve(progress=False)
print(f"z_steps_inner chosen: {auto.z_steps_inner}")
print(f"estimate {auto.step_error_estimate:.3g}; peak out {peak_field(auto)[-1]:.3f}")z_steps_inner chosen: 41
estimate 0.0505; peak out 0.836The repeated solves cost time, so for a calculation you will run many times, find the value once and set it. With a savefile the chosen value is stored with the result, and a later run loads both without solving again.
Runs that blow up
Some problems are unstable rather than just inaccurate. A dense medium with a coarse step can make the field grow without limit, or make the master-equation integrator fail. solve() stops as soon as the field is not finite, exceeds blowup_factor (1000 by default) times its initial peak, or the integrator fails, and raises PropagationInstabilityError:
dense = make(interaction_strengths=[3000.0], z_steps_inner=2)
try:
dense.solve(progress=False)
except PropagationInstabilityError as e:
print(e)
print("stopped at z =", e.z, "; suggested z_steps_inner:", e.suggested_z_steps_inner)Propagation became unstable at z=0.1: the master-equation integrator failed (Excess work done on this call. Try to increasing the nsteps parameter in the Options class). The z-step is probably too coarse for this medium; try z_steps_inner=4 (currently 2), z_steps_inner='auto', or more t_steps. Pass check_stability=False to solve() to suppress this check.
stopped at z = 0.1 ; suggested z_steps_inner: 4/Users/tpo/repos/clerq/.claude/worktrees/cyclic-output-parameters-9652ec/.venv/lib/python3.14/site-packages/scipy/integrate/_ode.py:418: UserWarning: _zvode: Excess work done on this call. (Perhaps wrong MF.)
self._y, self.t = mth(self.f, self.jac or (lambda: None),The fields up to the point of failure are still in the solver. The suggestion in the error is a starting point, not a guarantee. Here it is too small. Four inner steps still fail at this density, because the medium is far too optically thick for the pulse to be followed at all with a step of this kind. Try the suggestion, then lower the density or shorten the medium, or split it into more z_steps.
A medium with real gain, such as an inverted population, can legitimately amplify the field by more than 1000. Raise blowup_factor for those. check_stability=False turns off every check on this page.
Too coarse a time grid
Time is not marched like \(z\) is. The master equation is integrated with its own adaptive steps, and the time grid only fixes where the fields and coherences are recorded and interpolated. So a coarse grid cannot blow up. It makes the answer wrong in a way that looks reasonable, which is the harder case.
After solving, solve() takes the fields and the summed coherences at every \(z\), drops every other time sample, interpolates them back and compares. This is solver.t_step_error_estimate, relative to the peak, and it warns above t_step_error_warn.
print("t_steps estimate true error")
ref_t = make(t_steps=1280, z_steps_inner=128)
ref_t.solve(progress=False, check_stability=False)
exact_t = peak_field(ref_t)[-1]
with warnings.catch_warnings():
warnings.simplefilter("ignore")
for n in [20, 40, 80, 160, 320]:
s = make(t_steps=n, z_steps_inner=32)
s.solve(progress=False)
true_err = abs(peak_field(s)[-1] - exact_t) / exact_t
print(f"{n:7d} {s.t_step_error_estimate:8.3g} {true_err:10.3g}")t_steps estimate true error
20 1.47 0.258
40 0.929 0.0735
80 0.212 0.00326
160 0.0134 0.00321
320 0.00371 0.00318The true error stops falling at about 0.3% because that is the remaining \(z\)-step error at z_steps_inner=32. The estimate of the time grid keeps falling.
You can estimate a grid before solving. suggest_t_steps() finds the largest frequency in the problem (the bandwidth of the input fields, their Rabi frequencies, the detunings and the Doppler width) and returns enough steps for 10 samples per period of it:
make().suggest_t_steps()135This is an estimate from the input fields. The medium can add structure that is not there at the start, such as sidebands and ringing, so check it against the warning after solving.
The time window
The window t_min to t_max has to contain the whole pulse at every \(z\), including the delay the medium adds. solve() warns when a field is still changing in the first or last 5% of the window (window_edge_warn, window_edge_frac). In this problem the pulse is delayed by about 9 time units over the medium, so a window ending at 9 cuts it off:
clipped = make(t_max=9.0, z_steps_inner=32)
clipped.solve(progress=False)
print(f"field still changing at the window edge by {clipped.window_edge_estimate:.0%} of its peak")field still changing at the window edge by 58% of its peak/var/folders/41/hkxt57911pzcfzrv2p66w1z80000gn/T/ipykernel_90155/2846482427.py:2: UserWarning: A field is still changing near the edge of the time window (by 58.0% of its peak), so a pulse or its ringing may be cut off. Try widening t_min and t_max.
clipped.solve(progress=False)Summary
| To do | Use |
|---|---|
| Find out if a result is under-resolved | Read the warnings; check step_error_estimate, t_step_error_estimate and window_edge_estimate after solve() |
| Let clerq choose the \(z\)-step | z_steps_inner="auto" |
| Get a starting time grid | solver.suggest_t_steps() |
| Handle a run that blew up | Catch PropagationInstabilityError; raise z_steps_inner (or z_steps) or lower the density |
| Allow genuine gain | solve(blowup_factor=...) |
| Silence the checks | solve(check_stability=False) |
The warning thresholds are arguments of solve(): step_error_warn, t_step_error_warn and window_edge_warn. See the API reference.