Backward Euler Method University of Massachusetts Amherst
Backward Euler evaluates f f f next time step instead of the
current one, making it implicit : u n + 1 u_{n+1} u n + 1 unconditional
stability : the method is stable for any step size when
Re ( λ ) < 0 \operatorname{Re}(\lambda) < 0 Re ( λ ) < 0
Derivation ¶ Integrate u ′ = f ( t , u ) u' = f(t,u) u ′ = f ( t , u )
u ( t n + 1 ) = u ( t n ) + ∫ t n t n + 1 f ( s , u ( s ) ) d s u(t_{n+1}) = u(t_n) + \int_{t_n}^{t_{n+1}} f(s, u(s))\,ds u ( t n + 1 ) = u ( t n ) + ∫ t n t n + 1 f ( s , u ( s )) d s Forward Euler (Definition 1 left endpoint t n t_n t n right endpoint t n + 1 t_{n+1} t n + 1
∫ t n t n + 1 f ( s , u ( s ) ) d s ≈ h f ( t n + 1 , u ( t n + 1 ) ) \int_{t_n}^{t_{n+1}} f(s, u(s))\,ds \approx h\,f(t_{n+1}, u(t_{n+1})) ∫ t n t n + 1 f ( s , u ( s )) d s ≈ h f ( t n + 1 , u ( t n + 1 )) This gives the backward Euler (or implicit Euler) method:
The Cost of Implicitness ¶ Each step of backward Euler requires solving an equation for u n + 1 u_{n+1} u n + 1
Linear problems. For f ( t , u ) = λ u + g ( t ) f(t, u) = \lambda u + g(t) f ( t , u ) = λ u + g ( t )
u n + 1 = u n + h ( λ u n + 1 + g ( t n + 1 ) ) ⟹ u n + 1 = u n + h g ( t n + 1 ) 1 − h λ u_{n+1} = u_n + h(\lambda u_{n+1} + g(t_{n+1}))
\quad\implies\quad
u_{n+1} = \frac{u_n + hg(t_{n+1})}{1 - h\lambda} u n + 1 = u n + h ( λ u n + 1 + g ( t n + 1 )) ⟹ u n + 1 = 1 − hλ u n + h g ( t n + 1 ) No iteration needed.
Nonlinear problems. For general f f f
G ( u n + 1 ) = u n + 1 − u n − h f ( t n + 1 , u n + 1 ) = 0 G(u_{n+1}) = u_{n+1} - u_n - hf(t_{n+1}, u_{n+1}) = 0 G ( u n + 1 ) = u n + 1 − u n − h f ( t n + 1 , u n + 1 ) = 0 at every time step, using any of the root-finding methods from the
nonlinear equations chapter (e.g.,
Newton’s method or a fixed-point iteration). The previous time step
u n u_n u n h h h
See the Euler’s method notebook for an
interactive implementation.
Error Analysis and Convergence ¶ Following the same analysis as for
forward Euler
For backward Euler applied to u ′ = f ( t , u ) u' = f(t,u) u ′ = f ( t , u ) u ∈ C 2 u \in C^2 u ∈ C 2
τ n = − h 2 u ′ ′ ( t n + 1 ) + O ( h 2 ) = O ( h ) \tau_n = -\frac{h}{2}u''(t_{n+1}) + O(h^2) = O(h) τ n = − 2 h u ′′ ( t n + 1 ) + O ( h 2 ) = O ( h ) Backward Euler has consistency order p = 1 p = 1 p = 1
Proof 1
Taylor expand the exact solution about t n + 1 t_{n+1} t n + 1 t n t_n t n f f f
u ( t n ) = u ( t n + 1 ) − h u ′ ( t n + 1 ) + h 2 2 u ′ ′ ( t n + 1 ) + O ( h 3 ) u(t_n) = u(t_{n+1}) - hu'(t_{n+1}) + \frac{h^2}{2}u''(t_{n+1}) + O(h^3) u ( t n ) = u ( t n + 1 ) − h u ′ ( t n + 1 ) + 2 h 2 u ′′ ( t n + 1 ) + O ( h 3 ) Since u ′ ( t n + 1 ) = f ( t n + 1 , u ( t n + 1 ) ) u'(t_{n+1}) = f(t_{n+1}, u(t_{n+1})) u ′ ( t n + 1 ) = f ( t n + 1 , u ( t n + 1 ))
τ n = u ( t n + 1 ) − u ( t n ) h − f ( t n + 1 , u ( t n + 1 ) ) = − h 2 u ′ ′ ( t n + 1 ) + O ( h 2 ) \tau_n = \frac{u(t_{n+1}) - u(t_n)}{h} - f(t_{n+1}, u(t_{n+1}))
= -\frac{h}{2}u''(t_{n+1}) + O(h^2) τ n = h u ( t n + 1 ) − u ( t n ) − f ( t n + 1 , u ( t n + 1 )) = − 2 h u ′′ ( t n + 1 ) + O ( h 2 ) The sign of the leading term is opposite to forward Euler’s
(+ h 2 u ′ ′ +\frac{h}{2}u'' + 2 h u ′′
Under the same Lipschitz conditions as Theorem 1
∣ E n ∣ ≤ C L ( e L ( t n − t 0 ) − 1 ) h |E_n| \leq \frac{C}{L}\left(e^{L(t_n - t_0)} - 1\right) h ∣ E n ∣ ≤ L C ( e L ( t n − t 0 ) − 1 ) h Backward Euler converges with order 1.
Proof 2
The argument is identical to the forward Euler convergence proof
(Theorem 1 p = 1 p = 1 p = 1
Stability Analysis ¶ Apply backward Euler to the test equation u ′ = λ u u' = \lambda u u ′ = λ u
u n + 1 = u n + h λ u n + 1 ⟹ u n + 1 = 1 1 − h λ u n u_{n+1} = u_n + h\lambda\,u_{n+1}
\quad\implies\quad
u_{n+1} = \frac{1}{1 - h\lambda}\,u_n u n + 1 = u n + hλ u n + 1 ⟹ u n + 1 = 1 − hλ 1 u n The stability function is R ( z ) = 1 1 − z R(z) = \frac{1}{1-z} R ( z ) = 1 − z 1 z = h λ z = h\lambda z = hλ R ( z ) = 1 + z R(z) = 1 + z R ( z ) = 1 + z Definition 6
The stability region of backward Euler is
S = { z ∈ C : ∣ 1 1 − z ∣ ≤ 1 } = { z ∈ C : ∣ 1 − z ∣ ≥ 1 } \mathcal{S} = \left\{z \in \mathbb{C} : \left|\frac{1}{1-z}\right| \leq 1\right\}
= \{z \in \mathbb{C} : |1-z| \geq 1\} S = { z ∈ C : ∣ ∣ 1 − z 1 ∣ ∣ ≤ 1 } = { z ∈ C : ∣1 − z ∣ ≥ 1 } This is everything outside the disk of radius 1 centered at ( 1 , 0 ) (1, 0) ( 1 , 0 ) { Re ( z ) ≤ 0 } \{\operatorname{Re}(z) \leq 0\} { Re ( z ) ≤ 0 }
Stability region of backward Euler: everything outside the disk ∣ 1 − z ∣ ≥ 1 |1 - z| \geq 1 ∣1 − z ∣ ≥ 1 h > 0 h > 0 h > 0 Re ( λ ) < 0 \operatorname{Re}(\lambda) < 0 Re ( λ ) < 0
For Re ( λ ) < 0 \operatorname{Re}(\lambda) < 0 Re ( λ ) < 0 all h > 0 h > 0 h > 0 unconditionally stable .
Compare forward Euler: ∣ 1 + z ∣ ≤ 1 |1+z| \leq 1 ∣1 + z ∣ ≤ 1 h h h ( − 1 , 0 ) (-1,0) ( − 1 , 0 ) ∣ 1 − z ∣ ≥ 1 |1-z| \geq 1 ∣1 − z ∣ ≥ 1 ∣ λ ∣ |\lambda| ∣ λ ∣ h h h
Stiffness ¶ Stiff problems (Definition 7 ODE has
widely separated time scales. The fast components demand small steps for
explicit stability, even after those components have decayed to
negligible levels. This is where backward Euler’s unconditional stability
pays off.
Example 1 (The Prothero--Robinson Equation)
Consider u ′ ( t ) = λ ( u − cos t ) − sin t u'(t) = \lambda(u - \cos t) - \sin t u ′ ( t ) = λ ( u − cos t ) − sin t u ( 0 ) = 1 u(0) = 1 u ( 0 ) = 1
The exact solution is u ( t ) = cos t u(t) = \cos t u ( t ) = cos t λ \lambda λ
With λ = − 1000 \lambda = -1000 λ = − 1000 h < 0.002 h < 0.002 h < 0.002 O ( 1 ) O(1) O ( 1 )
Backward Euler with h = 0.1 h = 0.1 h = 0.1 ∣ R ( z ) ∣ = ∣ 1 / ( 1 − h λ ) ∣ = 1 / 101 ≪ 1 |R(z)| = |1/(1 - h\lambda)| = 1/101 \ll 1 ∣ R ( z ) ∣ = ∣1/ ( 1 − hλ ) ∣ = 1/101 ≪ 1
Example 2 (Heat Equation Semi-Discretization)
Discretize u t = u x x u_t = u_{xx} u t = u xx [ 0 , 1 ] [0,1] [ 0 , 1 ] Δ x \Delta x Δ x
d u j d t = u j + 1 − 2 u j + u j − 1 ( Δ x ) 2 \frac{du_j}{dt} = \frac{u_{j+1} - 2u_j + u_{j-1}}{(\Delta x)^2} d t d u j = ( Δ x ) 2 u j + 1 − 2 u j + u j − 1 This is a system u ′ = A u \mathbf{u}' = A\mathbf{u} u ′ = A u A A A
λ k = − 4 ( Δ x ) 2 sin 2 ( k π 2 ( N + 1 ) ) \lambda_k = -\frac{4}{(\Delta x)^2}\sin^2\!\left(\frac{k\pi}{2(N+1)}\right) λ k = − ( Δ x ) 2 4 sin 2 ( 2 ( N + 1 ) kπ ) The smallest eigenvalue is λ 1 ≈ − π 2 \lambda_1 \approx -\pi^2 λ 1 ≈ − π 2 λ N ≈ − 4 / ( Δ x ) 2 \lambda_N \approx -4/(\Delta x)^2 λ N ≈ − 4/ ( Δ x ) 2 ∣ λ N ∣ / ∣ λ 1 ∣ → ∞ |\lambda_N|/|\lambda_1| \to \infty ∣ λ N ∣/∣ λ 1 ∣ → ∞ Δ x → 0 \Delta x \to 0 Δ x → 0
Forward Euler requires h ≤ ( Δ x ) 2 / 2 h \leq (\Delta x)^2/2 h ≤ ( Δ x ) 2 /2 Δ x = 0.01 \Delta x = 0.01 Δ x = 0.01 h ≤ 5 × 1 0 − 5 h \leq 5 \times 10^{-5} h ≤ 5 × 1 0 − 5 T = 1 T = 1 T = 1
Backward Euler has no stability restriction. Choose h h h