---
title: Backward Euler Method (BEM)
url: https://www.emergentmind.com/topics/backward-euler-method-bem
type: topic
---

# Backward Euler Method (BEM)

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 [2106.16052]. 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 [2002.09203].

## 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 [2106.16052]. 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 [1709.01018].

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 [2002.09203]. 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 [2306.06689]. For additive fractional Brownian motion with Hurst parameter \(H>1/2\), the backward Euler–Maruyama step is
\[
Y_{k+1}=Y_k+h\,b(Y_{k+1})+\Delta B_k^H,
\]
which is fully implicit in the drift [2205.13659].

Multi-valued stochastic differential equations introduce an inclusion rather than an equation. On an equidistant grid \(t_n=nk\), one seeks \(X_n\in D(f)\) and \(\eta_n\in f(X_n)\) such that
\[
X_n+k\,\eta_n
=
X_{n-1}
+k\,b(X_n)
+g(X_{n-1})\,\Delta W_n,
\]
or equivalently \(X_n+k\,f(X_n)\ni X_{n-1}+k\,b(X_n)+g(X_{n-1})\Delta W_n\), emphasizing the multi-valued drift [1906.11538].

For semidiscrete PDE systems, BEM produces nonlinear algebraic systems. In nonlinear heat conduction, spatial discretization yields
\[
y'(t)=-A(y(t))\,y(t)+g(t),
\]
and backward Euler gives
\[
(I+\Delta t\,A(y^{n+1}))\,y^{n+1}=y^n+\Delta t\,g^{n+1},
\]
with nonlinearity entering through the stiffness matrix \(A(y^{n+1})\) [2211.04227]. 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 [2106.16052]. 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 [1805.05142].

## 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 \(\nu k<1\), where \(\nu\) is the one-sided Lipschitz constant from Assumption 3.1, the implicit equations \(U^n=U^{n-1}+k\,f(\xi_n,U^n)\) admit a unique adapted solution \(\{U^n\}\in L^2(\Omega;\mathbb R^d)\) [1709.01018]. 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
\[
\Phi(y)=y-h\,f(y),
\]
one has
\[
\Phi'(y)=1-h\,f'(y)\ge 1-h\,a_1>0,
\]
provided \(h\,a_1<1\). Together with
\[
\lim_{y\to 0^+}\Phi(y)=+\infty,
\qquad
\lim_{y\to +\infty}\Phi(y)=+\infty,
\]
this yields a unique positive solution \(Y_{n+1}\) of the algebraic step for every realization of \((\Delta W_n,\Delta N_n)\) [2002.09203]. The same paper notes that in practice one may solve for \(Y_{n+1}\) 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 \(f\) [2002.09203].

For multi-valued SDEs, well-posedness is expressed through monotone operator theory. Under Assumptions A–D and the step-size restriction \(8L_bk\le 1\), 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 [1906.11538]. In nonlinear heat conduction, solvability of the nonlinear system is addressed through Picard iteration:
\[
(I+\Delta t\,A(y^{(m)}))\,y^{(m+1)}=y^n+\Delta t\,g^{n+1}.
\]
If \(A(\cdot)\) is Lipschitz and
\[
0<\Delta t<\frac{1}{L(\|y^n\|+\Delta t\|g^{n+1}\|)},
\]
then the iterates converge to the true solution \(y^{n+1}\) of the backward Euler system [2211.04227].

In the fBm-driven SDE setting, the practical implication of one-sided Lipschitz dissipativity is that BEM remains stable under the mild condition \(\kappa h<1\), whereas explicit methods can fail for super-linear drifts [2205.13659]. 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 [2106.16052]. In the FVM–BEM interface setting, well-posedness is derived from ellipticity of the coupled bilinear form under
\[
\lambda_{\min}(A)-\tfrac14 C_K>0,
\qquad C_K<1,
\]
which yields coercivity in the natural \(H^1(\Omega)\times H^{-1/2}(\Gamma)\) norm [1805.05142].

## 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 \(Y_n\); 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 \(h\,a_1<1\) [2002.09203]. Numerical tests reinforce the analytical result: the explicit Euler method produces a large fraction of negative samples even for \(h=2^{-4}\), while BEM remains positive for all tested \(h\) [2002.09203].

For nonlinear heat conduction, monotonicity and boundedness are obtained under an M-matrix structure. If \(A(y)\) is symmetric positive semidefinite with nonpositive off-diagonal entries, \(g(t)\ge 0\), and \(y^0\ge 0\), then for any \(\Delta t>0\) and at every time step all Picard iterates satisfy \(y^{(m)}\ge 0\) elementwise, and
\[
\|y^{(m+1)}\|\le \|y^n\|+\Delta t\,\|g^{n+1}\|.
\]
The inverse positivity of \(I+\Delta t\,A(y^{(m)})\) is the mechanism behind this discrete admissibility [2211.04227].

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 [1709.01018]. 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
\[
\mu-2N\nu^{-1}\|f\|_\infty>0,
\qquad \nu=\mu+\gamma/\delta,
\]
one also obtains \(\|\nabla U_h^n\|\le C\) uniformly in \(n\) [2106.16052].

In the coupled FVM–BEM formulation, structural stability takes the form of discrete energy estimates rather than positivity. The fully discrete solution satisfies
\[
\|u_{h,\tau}\|_{H_T}+\|\phi_{h,\tau}\|_{B_T}
\lesssim
\|f\|_{L_T^2L_\Omega^2}
+\|q\|_{L^2(\Omega)}
+\|g_2\|_{L_T^2L_\Gamma^2}
+\|g_1\|_{L_T^2(H^{-1/2})'},
\]
and the proof avoids Grönwall-exponential growth by relying on a telescoping discrete energy identity and coercivity of the coupled bilinear form [1805.05142]. 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 | \(\sup_{0\le n\le N}\|Y_n-X_{t_n}\|_{L^2}\le C\,h^{1/2}\) | [2002.09203] |
| Randomized backward Euler, time-irregular coefficients | \(\max_{0\le n\le N}\|U^n-u(t_n)\|_{L^2(\Omega)}\le C\,k^{1/2}\) | [1709.01018] |
| Multi-valued SDEs | \(\max_n\|X(t_n)-X_n\|_{L^2(\Omega;\mathbb R^d)}\le C\,k^{1/4}\) | [1906.11538] |
| Random periodic solutions of semilinear SDEs | \(\sup_j E\|X_{t_j}-\tilde X_j\|^2\le C\,h^2\) | [2306.06689] |
| Additive fBm-driven SDEs, \(H>1/2\) | \(\sup_{t\in[0,T]}E|Y_t-X_t|^2\le C\,h^2\) | [2205.13659] |
| Oldroyd model, nonsmooth initial data | \(\|e^n\|_{L^2}\le K\,t_n^{-1/2}[h^2+k(1+\log(1/k))^{1/2}]\) | [2106.16052] |
| FVM–BEM coupling, variant BE | \(\|u-u_{h,\tau}\|_{H_T}=O(h^s+\tau^r)\) | [1805.05142] |
| FVM–BEM coupling, classical BE | \(\big[\sum_n\tau^n\|(u-u_h^n,\phi-\phi_h^n)\|_H^2\big]^{1/2}=O(h+\tau)\) | [1805.05142] |

The proof mechanisms are correspondingly diverse. For the jump-diffusion rate model, the order \(1/2\) follows from local truncation estimates \(\|\mathcal M_{n+1}\|_{L^2}\le C\,h\), conditional mean estimates \(\|\mathbb E[\mathcal M_{n+1}\mid\mathcal F_{t_n}]\|_{L^2}\le C\,h^{3/2}\), and a discrete Gronwall argument applied to the error recursion [2002.09203]. For randomized backward Euler, convergence in the root-mean-square norm with rate \(0.5\) 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 [1709.01018].

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 \(1/4\) [1906.11538]. By contrast, for semilinear SDEs with additive noise and random periodic solutions, improved \(L^p\)-bounds and a sharper local truncation estimate,
\[
\|\mathcal R_j\|_{L^2}=O(h^{3/2}),
\qquad
\|\mathbb E[\mathcal R_j\mid\mathcal F_{t_{j-1}}]\|_{L^2}=O(h^2),
\]
lead to first-order mean-square convergence under the simple spectral gap condition \(c_f<\lambda_1\) [2306.06689].

The fBm-driven case also attains order \(1\), but for a different reason: additive noise with \(H>1/2\) 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 [2205.13659]. That same paper proves an asymptotic error distribution result,
\[
\frac{Z_T}{h}\longrightarrow U_T \quad \text{in }L^p,\ 1\le p<2,
\]
which shows that the strong-order-\(1\) rate is optimal [2205.13659].

## 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 [1709.01018]. The same paper extends the method to infinite-dimensional evolution equations on a Gelfand triple \(V\subset H\cong H^*\subset V^*\), 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 [1709.01018].

In interface problems, two backward-Euler-type time discretizations are distinguished. The variant backward Euler method uses \(\omega\)-weighted temporal averages of the right-hand side, and the paper emphasizes that no further time-regularity on \(u\) is needed for the corresponding convergence bound [1805.05142]. The classical backward Euler method replaces the weighted averages by pointwise values \(f^n\), \(g_i^n\) at time \(t^n\); this is cheaper because no weights are required, but it needs stronger time-regularity such as \(f\in H^1(0,T;L^2(\Omega))\) and analogous assumptions for the boundary data [1805.05142].

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 [2211.04227]. 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 \(q_r^n(w)=k\sum_{j=1}^n\beta(t_n-t_j)w^j\) with \(\beta(t)=\gamma e^{-\delta t}\) [2106.16052].

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 \(N_t\), but that uniform grids suffice if \(\lambda\) is moderate [2002.09203]. 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 \(\mathbb E|Y_N-X_T|^2\) by \(10^4\) paths and plotting on a log–log scale gives a straight line of slope \(0.5\), confirming the theoretical mean-square rate, while the explicit Euler method produces negative samples whereas BEM remains positive for all tested \(h\) [2002.09203]. 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 \(1.0\)–\(1.2\) in two-dimensional tests [2205.13659].

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 \(0.5\) in the root-mean-square norm under only square-integrability in time [1709.01018]. For multi-valued stochastic differential equations, the proved rate is at least \(1/4\) [1906.11538]. Order one in mean square is obtained in more specialized settings: semilinear SDEs with additive noise and random periodic solutions under \(c_f<\lambda_1\) [2306.06689], and additive fractional Brownian noise with \(H>1/2\) under one-sided Lipschitz and polynomial-growth assumptions [2205.13659].

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 \(g\) at grid points where \(g=0\), leading to no convergence until \(k<p\), while randomized backward Euler converges with slope \(\approx 0.5\) already for \(k>p\) [1709.01018]. This motivates randomization as a device to protect the method against “fooling” by time-irregular data [1709.01018].

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 \(10^3\)–\(10^4\) nonlinear iterations, with each step invoking \(3\)–\(10\) Chebyshev iterations, and total mat-vec counts ranging from approximately \(5\,000\) to \(50\,000\) depending on \(\Delta t\) and \(N\) [2211.04227]. In the Oldroyd computations, P\(_2\)–P\(_0\) and MINI elements confirm the theoretical optimal and uniform-in-time error estimates, including first-order convergence in \(k\) and bounded errors for \(T\) up to \(O(50)\) [2106.16052]. In the FVM–BEM coupling examples, the observed rate is approximately first order for smooth convection-dominated problems and approximately \(h^{2/3}\) on an L-shaped domain with reduced regularity, exactly mirroring the regularity-dependent theory [1805.05142].

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.

Source: https://www.emergentmind.com/topics/backward-euler-method-bem