Backward Euler Method (BEM)

• Backward Euler Method (BEM) is an implicit time-stepping scheme that replaces the time derivative with a backward difference and is used in both deterministic and stochastic evolution problems.
• The method requires solving a nonlinear algebraic equation at each step to handle stiff or non-globally-Lipschitz dynamics, ensuring properties like positivity and monotonicity.
• Convergence orders and adaptations of BEM vary by model, with variants such as randomized BEM and extensions for multi-valued SDEs and jump-diffusion processes.

Searching arXiv for recent and foundational papers on the backward Euler method and closely related variants/applications. arxiv_search(query="backward Euler method nonlinear evolution equations stochastic differential equations", max_results=10) arxiv_search(query="Backward Euler method semilinear SDE random periodic solutions arXiv", max_results=5) arxiv_search(query="generalized Ait-Sahalia backward Euler Poisson jumps arXiv", max_results=5) The backward Euler method (BEM) is an implicit time-stepping method for evolution problems. For the abstract problem $u'(t)=F(u(t),t)$ on a uniform grid $t_n=nk$, it replaces the time derivative by a backward difference and evaluates the right-hand side at the new time,

$$\frac{u^n-u^{n-1}}{k}=F(u^n,t_n), \qquad u^n=u^{n-1}+k\,F(u^n,t_n),$$

and is therefore an implicit, first-order accurate time-stepping method (Bir et al., 2021). Across the arXiv literature, this template appears in deterministic ODEs, nonlinear evolution equations, semilinear and multi-valued SDEs, jump-diffusions, nonlinear heat conduction, viscoelastic flow, and coupled finite-volume/boundary-element discretizations. The resulting analyses emphasize well-posedness under one-sided Lipschitz or monotonicity hypotheses, preservation of qualitative structure such as positivity or nonnegativity, and convergence rates that depend strongly on the regularity and stochastic structure of the underlying problem (Zhao et al., 2020).

1. Core formulation and representative realizations

The defining feature of BEM is drift-implicit evaluation. In finite-dimensional deterministic form, the step $u^n=u^{n-1}+kF(u^n,t_n)$ requires solution of a nonlinear algebraic equation at every time level (Bir et al., 2021). For nonlinear evolution equations with time-irregular coefficients, the randomized backward Euler method modifies the evaluation time by drawing $\tau_n\sim \mathrm{Uniform}(0,1)$, setting $\xi_n=t_{n-1}+k\tau_n$, and solving

$U^n=U^{n-1}+k\,f(\xi_n,U^n),$

so the method remains implicit in the state variable while randomizing the temporal sampling point (Eisenmann et al., 2017).

In stochastic settings, the same pattern persists. For the generalized Ait-Sahalia-type rate model with Poisson jumps,

$dX_t=f(X_t)\,dt+g(X_t)\,dW_t+h(X_{t^-})\,dN_t,$

the backward Euler discretization is

$$Y_{n+1} = Y_n+h\,f(Y_{n+1}) +g(Y_n)\,\Delta W_n +\varphi(Y_n)\,\Delta N_n, \qquad Y_0=x_0,$$

with the drift handled implicitly and the diffusion and jump increments handled explicitly (Zhao et al., 2020). For semilinear SDEs with additive noise and random periodic solutions, the update takes the form

$$\tilde X_{j+1} = \tilde X_j -h\,\Lambda\,\tilde X_{j+1} +h\,f(t_{j+1},\tilde X_{j+1}) +g(t_j)\,\Delta W_j,$$

again separating an implicit drift from an explicit noise increment (Guo et al., 2023). For additive fractional Brownian motion with Hurst parameter $t_n=nk$0, the backward Euler–Maruyama step is

$t_n=nk$1

which is fully implicit in the drift (Zhou et al., 2022).

Multi-valued stochastic differential equations introduce an inclusion rather than an equation. On an equidistant grid $t_n=nk$2, one seeks $t_n=nk$3 and $t_n=nk$4 such that

$t_n=nk$5

or equivalently $t_n=nk$6, emphasizing the multi-valued drift (Eisenmann et al., 2019).

For semidiscrete PDE systems, BEM produces nonlinear algebraic systems. In nonlinear heat conduction, spatial discretization yields

$t_n=nk$7

and backward Euler gives

$t_n=nk$8

with nonlinearity entering through the stiffness matrix $t_n=nk$9 (Botchev et al., 2022). In the Oldroyd model of order one, backward Euler is combined with finite elements, a divergence-free constraint, a nonlinear convection term, and a quadrature approximation of the memory convolution (Bir et al., 2021). In coupled parabolic-elliptic interface problems, backward-Euler-type time discretizations appear inside a non-symmetric FVM–BEM coupling, with mass-lumping and boundary integral operators entering the fully discrete formulation (Erath et al., 2018).

2. Solvability of the implicit step

Because BEM is implicit, its first analytical question is the existence and uniqueness of each time step. In nonlinear evolution equations with time-irregular coefficients, if $$\frac{u^n-u^{n-1}}{k}=F(u^n,t_n), \qquad u^n=u^{n-1}+k\,F(u^n,t_n),$$0, where $$\frac{u^n-u^{n-1}}{k}=F(u^n,t_n), \qquad u^n=u^{n-1}+k\,F(u^n,t_n),$$1 is the one-sided Lipschitz constant from Assumption 3.1, the implicit equations $$\frac{u^n-u^{n-1}}{k}=F(u^n,t_n), \qquad u^n=u^{n-1}+k\,F(u^n,t_n),$$2 admit a unique adapted solution $$\frac{u^n-u^{n-1}}{k}=F(u^n,t_n), \qquad u^n=u^{n-1}+k\,F(u^n,t_n),$$3 (Eisenmann et al., 2017). This condition is the discrete counterpart of monotonicity-based solvability.

For the generalized Ait-Sahalia-type rate model with Poisson jumps, solvability is tied to positivity. Writing

$$\frac{u^n-u^{n-1}}{k}=F(u^n,t_n), \qquad u^n=u^{n-1}+k\,F(u^n,t_n),$$4

one has

$$\frac{u^n-u^{n-1}}{k}=F(u^n,t_n), \qquad u^n=u^{n-1}+k\,F(u^n,t_n),$$5

provided $$\frac{u^n-u^{n-1}}{k}=F(u^n,t_n), \qquad u^n=u^{n-1}+k\,F(u^n,t_n),$$6. Together with

$$\frac{u^n-u^{n-1}}{k}=F(u^n,t_n), \qquad u^n=u^{n-1}+k\,F(u^n,t_n),$$7

this yields a unique positive solution $$\frac{u^n-u^{n-1}}{k}=F(u^n,t_n), \qquad u^n=u^{n-1}+k\,F(u^n,t_n),$$8 of the algebraic step for every realization of $$\frac{u^n-u^{n-1}}{k}=F(u^n,t_n), \qquad u^n=u^{n-1}+k\,F(u^n,t_n),$$9 (Zhao et al., 2020). The same paper notes that in practice one may solve for $u^n=u^{n-1}+kF(u^n,t_n)$0 by a fixed-point iteration or by a few steps of Newton’s method, and that convergence is rapid under the one-sided-Lipschitz condition on $u^n=u^{n-1}+kF(u^n,t_n)$1 (Zhao et al., 2020).

For multi-valued SDEs, well-posedness is expressed through monotone operator theory. Under Assumptions A–D and the step-size restriction $u^n=u^{n-1}+kF(u^n,t_n)$2, the backward Euler–Maruyama method is well-defined; the analysis uses maximal monotonicity of the drift, coercivity, polynomial growth, and the Savaré–Verdi extra condition (Eisenmann et al., 2019). In nonlinear heat conduction, solvability of the nonlinear system is addressed through Picard iteration: $u^n=u^{n-1}+kF(u^n,t_n)$3 If $u^n=u^{n-1}+kF(u^n,t_n)$4 is Lipschitz and

$u^n=u^{n-1}+kF(u^n,t_n)$5

then the iterates converge to the true solution $u^n=u^{n-1}+kF(u^n,t_n)$6 of the backward Euler system (Botchev et al., 2022).

In the fBm-driven SDE setting, the practical implication of one-sided Lipschitz dissipativity is that BEM remains stable under the mild condition $u^n=u^{n-1}+kF(u^n,t_n)$7, whereas explicit methods can fail for super-linear drifts (Zhou et al., 2022). In the Oldroyd discretization, the nonlinear algebraic problem is solved, for example, by fixed-point or Newton iteration for the convective term together with enforcement of the divergence-free constraint through a pressure solve (Bir et al., 2021). In the FVM–BEM interface setting, well-posedness is derived from ellipticity of the coupled bilinear form under

$u^n=u^{n-1}+kF(u^n,t_n)$8

which yields coercivity in the natural $u^n=u^{n-1}+kF(u^n,t_n)$9 norm (Erath et al., 2018).

3. Structural properties: positivity, monotonicity, boundedness, and stability

A central reason for using BEM is its ability to preserve qualitative features that explicit schemes may destroy. For the generalized Ait-Sahalia-type rate model with Poisson jumps, the method preserves positivity of each $\tau_n\sim \mathrm{Uniform}(0,1)$0; the paper explicitly states that the BEM preserves positivity of the original problem and proves that the algebraic step has a unique positive solution whenever $\tau_n\sim \mathrm{Uniform}(0,1)$1 (Zhao et al., 2020). Numerical tests reinforce the analytical result: the explicit Euler method produces a large fraction of negative samples even for $\tau_n\sim \mathrm{Uniform}(0,1)$2, while BEM remains positive for all tested $\tau_n\sim \mathrm{Uniform}(0,1)$3 (Zhao et al., 2020).

For nonlinear heat conduction, monotonicity and boundedness are obtained under an M-matrix structure. If $\tau_n\sim \mathrm{Uniform}(0,1)$4 is symmetric positive semidefinite with nonpositive off-diagonal entries, $\tau_n\sim \mathrm{Uniform}(0,1)$5, and $\tau_n\sim \mathrm{Uniform}(0,1)$6, then for any $\tau_n\sim \mathrm{Uniform}(0,1)$7 and at every time step all Picard iterates satisfy $\tau_n\sim \mathrm{Uniform}(0,1)$8 elementwise, and

$\tau_n\sim \mathrm{Uniform}(0,1)$9

The inverse positivity of $\xi_n=t_{n-1}+k\tau_n$0 is the mechanism behind this discrete admissibility (Botchev et al., 2022).

In randomized backward Euler for stiff ODEs, the scheme is contrasted with explicit discretization through a classical stability statement: randomized forward Euler oscillates or explodes unless the step is sufficiently small, whereas backward Euler is A-stable (Eisenmann et al., 2017). For the Oldroyd model, the main structural result is long-time control. A weighted Gronwall argument and positivity of the quadrature rule yield a uniform bound in the Dirichlet norm, and under the uniqueness condition

$\xi_n=t_{n-1}+k\tau_n$1

one also obtains $\xi_n=t_{n-1}+k\tau_n$2 uniformly in $\xi_n=t_{n-1}+k\tau_n$3 (Bir et al., 2021).

In the coupled FVM–BEM formulation, structural stability takes the form of discrete energy estimates rather than positivity. The fully discrete solution satisfies

$\xi_n=t_{n-1}+k\tau_n$4

and the proof avoids Grönwall-exponential growth by relying on a telescoping discrete energy identity and coercivity of the coupled bilinear form (Erath et al., 2018). A plausible implication is that BEM is especially attractive when the numerical method is required to preserve positivity, nonnegativity, or dissipative energy structure in addition to approximating the solution.

4. Convergence theory and proved orders

The convergence order of BEM is not universal; it depends on drift structure, noise model, regularity, and the quantity being estimated. The following statements are explicitly proved in the cited works.

Setting	Convergence statement	Source
Generalized Ait-Sahalia-type rate model with Poisson jumps	$\xi_n=t_{n-1}+k\tau_n$5	(Zhao et al., 2020)
Randomized backward Euler, time-irregular coefficients	$\xi_n=t_{n-1}+k\tau_n$6	(Eisenmann et al., 2017)
Multi-valued SDEs	$\xi_n=t_{n-1}+k\tau_n$7	(Eisenmann et al., 2019)
Random periodic solutions of semilinear SDEs	$\xi_n=t_{n-1}+k\tau_n$8	(Guo et al., 2023)
Additive fBm-driven SDEs, $\xi_n=t_{n-1}+k\tau_n$9	$U^n=U^{n-1}+k\,f(\xi_n,U^n),$0	(Zhou et al., 2022)
Oldroyd model, nonsmooth initial data	$U^n=U^{n-1}+k\,f(\xi_n,U^n),$1	(Bir et al., 2021)
FVM–BEM coupling, variant BE	$U^n=U^{n-1}+k\,f(\xi_n,U^n),$2	(Erath et al., 2018)
FVM–BEM coupling, classical BE	$U^n=U^{n-1}+k\,f(\xi_n,U^n),$3	(Erath et al., 2018)

The proof mechanisms are correspondingly diverse. For the jump-diffusion rate model, the order $U^n=U^{n-1}+k\,f(\xi_n,U^n),$4 follows from local truncation estimates $U^n=U^{n-1}+k\,f(\xi_n,U^n),$5, conditional mean estimates $U^n=U^{n-1}+k\,f(\xi_n,U^n),$6, and a discrete Gronwall argument applied to the error recursion (Zhao et al., 2020). For randomized backward Euler, convergence in the root-mean-square norm with rate $U^n=U^{n-1}+k\,f(\xi_n,U^n),$7 is established under only square-integrability of the coefficient function with respect to the temporal parameter, and the stability estimate explicitly contains both the local residual and its conditional expectation (Eisenmann et al., 2017).

The multi-valued SDE analysis is notably different. It combines a deterministic error representation and stability argument for the drift part, based on Nochetto–Savaré–Verdi techniques, with stochastic remainder bounds of Burkholder–Davis–Gundy type; the final result is strong convergence of order at least $U^n=U^{n-1}+k\,f(\xi_n,U^n),$8 (Eisenmann et al., 2019). By contrast, for semilinear SDEs with additive noise and random periodic solutions, improved $U^n=U^{n-1}+k\,f(\xi_n,U^n),$9-bounds and a sharper local truncation estimate,

$dX_t=f(X_t)\,dt+g(X_t)\,dW_t+h(X_{t^-})\,dN_t,$0

lead to first-order mean-square convergence under the simple spectral gap condition $dX_t=f(X_t)\,dt+g(X_t)\,dW_t+h(X_{t^-})\,dN_t,$1 (Guo et al., 2023).

The fBm-driven case also attains order $dX_t=f(X_t)\,dt+g(X_t)\,dW_t+h(X_{t^-})\,dN_t,$2, but for a different reason: additive noise with $dX_t=f(X_t)\,dt+g(X_t)\,dW_t+h(X_{t^-})\,dN_t,$3 permits a proof based on a variation-of-constants representation, Young–Hölder estimates, Malliavin integration by parts, and moment bounds on the exact solution and its Malliavin derivatives (Zhou et al., 2022). That same paper proves an asymptotic error distribution result,

$dX_t=f(X_t)\,dt+g(X_t)\,dW_t+h(X_{t^-})\,dN_t,$4

which shows that the strong-order-$dX_t=f(X_t)\,dt+g(X_t)\,dW_t+h(X_{t^-})\,dN_t,$5 rate is optimal (Zhou et al., 2022).

5. Variants and extensions

Several works study BEM not as a single fixed scheme but as a family of closely related implicit discretizations adapted to irregular data, infinite-dimensional settings, or coupled discretizations. The randomized backward Euler method is designed for Carathéodory-type functions with time-irregular coefficients. It randomizes the evaluation node within each interval, requires no extra quadrature, and has the same nonlinear-solve cost as classical BEM plus one uniform random draw per step (Eisenmann et al., 2017). The same paper extends the method to infinite-dimensional evolution equations on a Gelfand triple $dX_t=f(X_t)\,dt+g(X_t)\,dW_t+h(X_{t^-})\,dN_t,$6, combining randomized backward Euler in time with a Galerkin finite element method in space and deriving an error estimate that reflects temporal Hölder regularity and spatial approximation errors (Eisenmann et al., 2017).

In interface problems, two backward-Euler-type time discretizations are distinguished. The variant backward Euler method uses $dX_t=f(X_t)\,dt+g(X_t)\,dW_t+h(X_{t^-})\,dN_t,$7-weighted temporal averages of the right-hand side, and the paper emphasizes that no further time-regularity on $dX_t=f(X_t)\,dt+g(X_t)\,dW_t+h(X_{t^-})\,dN_t,$8 is needed for the corresponding convergence bound (Erath et al., 2018). The classical backward Euler method replaces the weighted averages by pointwise values $dX_t=f(X_t)\,dt+g(X_t)\,dW_t+h(X_{t^-})\,dN_t,$9, $$Y_{n+1} = Y_n+h\,f(Y_{n+1}) +g(Y_n)\,\Delta W_n +\varphi(Y_n)\,\Delta N_n, \qquad Y_0=x_0,$$0 at time $$Y_{n+1} = Y_n+h\,f(Y_{n+1}) +g(Y_n)\,\Delta W_n +\varphi(Y_n)\,\Delta N_n, \qquad Y_0=x_0,$$1; this is cheaper because no weights are required, but it needs stronger time-regularity such as $$Y_{n+1} = Y_n+h\,f(Y_{n+1}) +g(Y_n)\,\Delta W_n +\varphi(Y_n)\,\Delta N_n, \qquad Y_0=x_0,$$2 and analogous assumptions for the boundary data (Erath et al., 2018).

In nonlinear heat conduction, the backward Euler step is combined with nonlinear iterations, and the implementation may solve the inner linear systems by a Chebyshev-based local-iteration scheme. The total cost is then counted in “mat-vecs,” and the paper compares this strategy to a nonlinear exponential Euler scheme based on restarted Krylov subspace methods (Botchev et al., 2022). In the Oldroyd model, BEM is not merely a time integrator but part of a fully discrete Galerkin formulation involving divergence-free finite element spaces, an LBB-stable pressure space, and a right-rectangle quadrature for the memory term $$Y_{n+1} = Y_n+h\,f(Y_{n+1}) +g(Y_n)\,\Delta W_n +\varphi(Y_n)\,\Delta N_n, \qquad Y_0=x_0,$$3 with $$Y_{n+1} = Y_n+h\,f(Y_{n+1}) +g(Y_n)\,\Delta W_n +\varphi(Y_n)\,\Delta N_n, \qquad Y_0=x_0,$$4 (Bir et al., 2021).

The stochastic literature also exhibits problem-specific adaptations. For the Ait-Sahalia-type model with Poisson jumps, the paper notes that for jump-adapted BEM one can synchronize time-grid points with jump times of $$Y_{n+1} = Y_n+h\,f(Y_{n+1}) +g(Y_n)\,\Delta W_n +\varphi(Y_n)\,\Delta N_n, \qquad Y_0=x_0,$$5, but that uniform grids suffice if $$Y_{n+1} = Y_n+h\,f(Y_{n+1}) +g(Y_n)\,\Delta W_n +\varphi(Y_n)\,\Delta N_n, \qquad Y_0=x_0,$$6 is moderate (Zhao et al., 2020). A plausible implication is that the backward Euler framework is best understood as a structural template—implicit treatment of stiff or monotone drift—onto which problem-specific choices of noise discretization, temporal sampling, and spatial approximation are attached.

6. Numerical behavior, comparisons, and recurrent misconceptions

The numerical evidence reported in the cited works is largely consistent with the theoretical emphasis on stability and structural preservation. In the generalized Ait-Sahalia-type rate model, computing $$Y_{n+1} = Y_n+h\,f(Y_{n+1}) +g(Y_n)\,\Delta W_n +\varphi(Y_n)\,\Delta N_n, \qquad Y_0=x_0,$$7 by $$Y_{n+1} = Y_n+h\,f(Y_{n+1}) +g(Y_n)\,\Delta W_n +\varphi(Y_n)\,\Delta N_n, \qquad Y_0=x_0,$$8 paths and plotting on a log–log scale gives a straight line of slope $$Y_{n+1} = Y_n+h\,f(Y_{n+1}) +g(Y_n)\,\Delta W_n +\varphi(Y_n)\,\Delta N_n, \qquad Y_0=x_0,$$9, confirming the theoretical mean-square rate, while the explicit Euler method produces negative samples whereas BEM remains positive for all tested $$\tilde X_{j+1} = \tilde X_j -h\,\Lambda\,\tilde X_{j+1} +h\,f(t_{j+1},\tilde X_{j+1}) +g(t_j)\,\Delta W_j,$$0 (Zhao et al., 2020). In additive fBm-driven problems with cubic or dissipative polynomial drift, explicit Euler–Maruyama and Crank–Nicolson may blow up at coarse step sizes, whereas BEM remains stable and attains a fitted slope in the vicinity of $$\tilde X_{j+1} = \tilde X_j -h\,\Lambda\,\tilde X_{j+1} +h\,f(t_{j+1},\tilde X_{j+1}) +g(t_j)\,\Delta W_j,$$1–$$\tilde X_{j+1} = \tilde X_j -h\,\Lambda\,\tilde X_{j+1} +h\,f(t_{j+1},\tilde X_{j+1}) +g(t_j)\,\Delta W_j,$$2 in two-dimensional tests (Zhou et al., 2022).

A common misconception is that backward Euler automatically improves convergence order simply because it is implicit. The available results do not support such a blanket statement. For nonlinear evolution equations with time-irregular coefficients, randomized backward Euler converges with rate $$\tilde X_{j+1} = \tilde X_j -h\,\Lambda\,\tilde X_{j+1} +h\,f(t_{j+1},\tilde X_{j+1}) +g(t_j)\,\Delta W_j,$$3 in the root-mean-square norm under only square-integrability in time (Eisenmann et al., 2017). For multi-valued stochastic differential equations, the proved rate is at least $$\tilde X_{j+1} = \tilde X_j -h\,\Lambda\,\tilde X_{j+1} +h\,f(t_{j+1},\tilde X_{j+1}) +g(t_j)\,\Delta W_j,$$4 (Eisenmann et al., 2019). Order one in mean square is obtained in more specialized settings: semilinear SDEs with additive noise and random periodic solutions under $$\tilde X_{j+1} = \tilde X_j -h\,\Lambda\,\tilde X_{j+1} +h\,f(t_{j+1},\tilde X_{j+1}) +g(t_j)\,\Delta W_j,$$5 (Guo et al., 2023), and additive fractional Brownian noise with $$\tilde X_{j+1} = \tilde X_j -h\,\Lambda\,\tilde X_{j+1} +h\,f(t_{j+1},\tilde X_{j+1}) +g(t_j)\,\Delta W_j,$$6 under one-sided Lipschitz and polynomial-growth assumptions (Zhou et al., 2022).

Another misconception is that the classical deterministic grid-point sampling of the right-hand side is always harmless. The randomized backward Euler paper gives a counterexample through the Prothero–Robinson ODE: classical BEM only samples $$\tilde X_{j+1} = \tilde X_j -h\,\Lambda\,\tilde X_{j+1} +h\,f(t_{j+1},\tilde X_{j+1}) +g(t_j)\,\Delta W_j,$$7 at grid points where $$\tilde X_{j+1} = \tilde X_j -h\,\Lambda\,\tilde X_{j+1} +h\,f(t_{j+1},\tilde X_{j+1}) +g(t_j)\,\Delta W_j,$$8, leading to no convergence until $$\tilde X_{j+1} = \tilde X_j -h\,\Lambda\,\tilde X_{j+1} +h\,f(t_{j+1},\tilde X_{j+1}) +g(t_j)\,\Delta W_j,$$9, while randomized backward Euler converges with slope $t_n=nk$00 already for $t_n=nk$01 (Eisenmann et al., 2017). This motivates randomization as a device to protect the method against “fooling” by time-irregular data (Eisenmann et al., 2017).

In PDE applications, computational cost and nonlinear solver behavior are equally important. For nonlinear heat conduction, backward Euler with local iteration may require up to $t_n=nk$02–$t_n=nk$03 nonlinear iterations, with each step invoking $t_n=nk$04–$t_n=nk$05 Chebyshev iterations, and total mat-vec counts ranging from approximately $t_n=nk$06 to $t_n=nk$07 depending on $t_n=nk$08 and $t_n=nk$09 (Botchev et al., 2022). In the Oldroyd computations, P$t_n=nk$10–P$t_n=nk$11 and MINI elements confirm the theoretical optimal and uniform-in-time error estimates, including first-order convergence in $t_n=nk$12 and bounded errors for $t_n=nk$13 up to $t_n=nk$14 (Bir et al., 2021). In the FVM–BEM coupling examples, the observed rate is approximately first order for smooth convection-dominated problems and approximately $t_n=nk$15 on an L-shaped domain with reduced regularity, exactly mirroring the regularity-dependent theory (Erath et al., 2018).

Taken together, these results portray BEM as an implicit discretization paradigm whose principal strengths are robustness under one-sided Lipschitz or monotone drift, compatibility with stiff and non-globally-Lipschitz dynamics, and preservation of structural constraints such as positivity or nonnegativity. The precise accuracy statement, however, is model-dependent rather than universal.