Backward Euler Method

• The backward Euler method is an implicit numerical scheme that replaces continuous time derivatives with stable, implicit updates for stiff and nonlinear problems.
• It achieves strong stability and robust convergence properties in both deterministic and stochastic differential equations by solving a nonlinear system at each time step.
• Extensions like the split-step, randomized, and filtered variants enable accurate simulation of complex systems with delays, jumps, and non-globally Lipschitz coefficients.

The backward Euler method is a fully implicit, single-step time discretization scheme for solving ordinary and stochastic differential equations, as well as more general evolution equations and partial differential equations. It is especially prominent in the numerical approximation of stiff, non-globally Lipschitz, or highly nonlinear problems, offering strong stability, robustness, and—under suitable conditions—optimal or provable rates of convergence in the mean square, $L^p$, or weak senses. The method and its variants underpin a broad spectrum of modern simulation and analysis tasks in deterministic and stochastic computation, including applications with memory, delay, jumps, multi-valued maps, and degenerate potentials.

1. Core Formulation and Methodology

The backward Euler method discretizes the time derivative of a differential equation, replacing the continuous evolution by a fully implicit update. For a general system

$u'(t) = f(u(t)), \quad u(0) = u_0,$

the scheme with time step $h$ and $u_n$ approximating $u(nh)$ is

$u_{n+1} = u_n + h\,f(u_{n+1}).$

In practice, at each new time step, a nonlinear system (or inclusion, if $f$ is multivalued) must be solved to obtain $u_{n+1}$.

For stochastic differential equations (SDEs),

$dX_t = b(X_t)\,dt + \sigma(X_t)\,dW_t,$

the corresponding backward Euler–Maruyama (BEM) scheme reads

$$\bar{X}_{n+1} = \bar{X}_n + b(\bar{X}_{n+1})\,h + \sigma(\bar{X}_n) \Delta W_n,$$

where $\Delta W_n = W_{(n+1)h} - W_{nh}$ is the Wiener increment.

Several sophisticated extensions exist:

• Split-step backward Euler methods for SDEs with delays (Wang et al., 2011): The drift is solved implicitly (possibly with delayed arguments interpolated), while the diffusion is handled explicitly.
• Randomized backward Euler (Eisenmann et al., 2017): Randomized integration points are introduced to handle time-irregular coefficients.
• Time-filtered backward Euler (Cibik et al., 2019, Demir et al., 2021): A modular second-order filter is applied after the implicit step for higher temporal accuracy and conservation.

2. Convergence, Stability, and Error Analysis

The backward Euler method’s haLLMark is unconditional stability for a wide range of stiff and monotonic problems. The key convergence results include:

• Mean-square or root-mean-square convergence: For SDEs, SDDEs, and NSDDEs with superlinear, non-globally Lipschitz coefficients or delays, the method attains strong order $1/2$ in mean-square sense for stepsize $h$ (Wang et al., 2011, Liu et al., 2022, Cai et al., 14 Feb 2024).
• Strong order 1 convergence: For drift-dominated SDEs with additive fractional Brownian motion and one-sided Lipschitz drift (Zhou et al., 2022), and for random periodic SDEs (Guo et al., 2023), the scheme achieves $O(h)$ mean-square error.
• Order reduction for low regularity or multivalued drift: For multi-valued SDEs or stochastic gradient flows, the scheme may attain only $1/4$ strong order or worse, depending on drift regularity (Eisenmann et al., 2019).
• Weak error and invariant measures: For ergodic SDEs, the BEM approximates invariant measures with order 1 in the weak sense (Chen et al., 2019, Liu et al., 2022, Jin, 2023).

Stability results include:

• Exponential mean-square stability: The method inherits or preserves mean-square exponential stability from the continuous system, without step size restriction in many cases (Wang et al., 2011, Liu et al., 2022, Cai et al., 14 Feb 2024).
• Positivity preservation: When properly designed with implicit handling of the drift, positivity of the solution is maintained, critical for financial and population models (Zhao et al., 2020).
• Long-time uniform bounds: Uniform Dirichlet or $L^2$ energy estimates hold over arbitrary time intervals for parabolic and viscoelastic equations with the backward Euler method (Goswami et al., 2012, Bir et al., 2021).

3. Extensions and Applications in Complex Systems

The backward Euler method, and its elaborations, have been applied and theoretically analyzed in many advanced contexts:

Stochastic Delay Equations (SDDEs, NSDDEs)

The improved split-step versions handle variable delays, with piecewise (linear or constant) delayed value interpolation ensuring convergence and preserving stability for generic, non-uniform delays (Wang et al., 2011, Liu et al., 2022, Cai et al., 14 Feb 2024).

Evolution Equations and PDEs

For time-discretized parabolic PDEs, including nonlinear heat conduction (Botchev et al., 2022), viscoelastic Oldroyd models (Goswami et al., 2012, Bir et al., 2021), and incompressible Navier–Stokes EMAC (Demir et al., 2021), backward Euler is often embedded in operator splitting, Galerkin, DPG, or structure-preserving spatial approximations. Time stepping proceeds by solving large but robust implicit or linearized systems at each step.

Stochastic Dynamics with Superlinear Growth, Jumps, or Degeneracy

For SDEs with superlinear drifts or diffusion, including multiplicative or jump noise, polynomial growth, and multi-valued subdifferential drifts, the method remains stable and convergent (with possible order reduction), provided monotonicity or coercivity conditions are satisfied (Zhao et al., 2020, Eisenmann et al., 2019, Liu et al., 2022).

Chemotaxis and Conservation Laws

In models such as chemotaxis (Keller–Segel systems), backward Euler is coupled with finite element discretizations and specialized positivity-preserving, flux-limited stabilization to ensure that physically relevant solutions (e.g., nonnegative densities) are maintained (Chatzipantelidis et al., 2022).

Invariant Measures and Long-time Behavior

The BEM numerically reproduces unique invariant (ergodic) measures under dissipativity/monotonicity, superlinear growth, or degenerate noise (Chen et al., 2019, Liu et al., 2022, Jin, 2023), with convergence of order $O(h)$ in Wasserstein or weak topology and, for temporal averages, associated central limit theorems for ergodic estimators (Jin, 2023).

4. Implementation Issues and Algorithmic Details

Key aspects for implementation include:

• Nonlinear solver at each step: Implicitness requires solving a nonlinear or even multi-valued problem at each stage; monotonicity and coercivity underpin convergence and well-posedness.
• Interpolation in delay equations: For variable delays, construction of $\tilde{y}_n^*$ must use careful interpolation (piecewise constant, linear) matching the grid and delay function (Wang et al., 2011).
• Stabilized and positivity-preserving postprocessing: In conservation/balance law PDEs, implicit discretization is often coupled with correction filters or limiters to guarantee mass conservation and avoid unphysical oscillations (Chatzipantelidis et al., 2022, Demir et al., 2021).
• Order-raising time filters: Post-processed solutions (linear combinations of three successive time levels) can increase the temporal order to 2 and restore conservation properties with only minor code changes (Cibik et al., 2019, Demir et al., 2021).
• Randomized node selection: For ODEs/PDEs with time-irregular coefficients, randomizing the integration point in each step (within the time interval) enhances error regularization and prevents aliasing by pathological coefficients (Eisenmann et al., 2017).

5. Numerical Experiments and Empirical Findings

Theoretical results are consistently corroborated by numerical experiments:

• Convergence rates measured versus $h$: Empirical mean-square errors consistently validate theoretical order $1/2$ or 1, even in the presence of highly nonlinear growth, non-globally Lipschitz drifts, or noisy coefficients (Wang et al., 2011, Zhao et al., 2020, Zhou et al., 2022, Guo et al., 2023, Cai et al., 14 Feb 2024).
• Comparison with explicit Euler and alternatives: In stiff or degenerate cases, explicit schemes become unstable or inaccurate; backward Euler–type methods preserve stability and often yield better accuracy for moderate to large $h$ (Wang et al., 2011, Zhao et al., 2020, Zhou et al., 2022).
• Long-time behavior and ergodicity: The method robustly approximates invariant measures and yields stable distributions in ergodic or stationary regimes (Chen et al., 2019, Liu et al., 2022, Jin, 2023).
• Preservation of physical invariants: Time-filtered backward Euler methods exactly or nearly conserve energy, momentum, and angular momentum in MHD or Navier–Stokes simulations, with documented reduction in spurious dissipation over long time intervals (Cibik et al., 2019, Demir et al., 2021).
• Stability on coarse grids: For highly oscillatory or time-irregular test problems, the method—especially when randomized—demonstrates strong error control even for coarse discretization, outperforming classical explicit or non-randomized methods (Eisenmann et al., 2017).

6. Limitations and Trade-offs

Notwithstanding its strengths, the backward Euler method presents several methodological compromises:

• Computational cost per step: Solving implicit or nonlinear equations/inclusions at each time step is more expensive than explicit methods; however, this is offset by larger allowable step sizes for stiff or dissipative systems.
• Order reduction in weak or strong sense: In low regularity, non-smooth multivalued drift, or irregular stochasticity scenarios, the convergence order may drop below $1/2$ (e.g., $1/4$ in MSDEs) (Eisenmann et al., 2019).
• Accuracy vs. stability: The method favors stability and long-time fidelity over high-order accuracy; for fine accuracy or smooth problems, higher-order or exponential integrators may be preferable (Botchev et al., 2022).
• Parameter tuning and initialization: For positivity preservation or long-time stability, explicit bounds on step sizes may be needed; filtering and limiters require careful choice of coefficients for conservation and non-negativity.

7. Broader Implications and Advancements

The backward Euler method, in modern variants and under rigorous analysis, provides a foundation for robust, structure-preserving, and stable time discretization across deterministic and stochastic modeling. Its central role is due to:

• Monotonicity and coercivity handling for stiff, dissipative, and non-linear systems.
• Intrinsic ability to manage non-Lipschitz growth where explicit methods fail.
• Versatility in integrating delay, jump, memory, and randomness, as required in SDEs, SDDEs, NSDDEs, SPDEs, and PDEs.
• Suitability for ergodic and long-run statistics, including central limit behavior of statistical estimators (Jin, 2023).

Continued advances in error analysis, coupling with spatial discretization, and modular higher-order postprocessing methods continue to expand the scope, accuracy, and applicability of the backward Euler methodology in computational mathematics and applied sciences.