Backward Euler Method as an Implicit Scheme
The forward Euler method in the previous tutorial is an explicit numerical scheme that uses the value of the derivative available at the current time step. On the contrary, we also have the backward Euler method as an implicit numerical scheme, which involves the value of the derivative at the next time step. Consider a general first-order ODE again
\begin{equation}
\frac{dy}{dt} = f(t,y) \tag{1}
\end{equation}
To integrate over a single time step, this time we assume that \( \frac{dy}{dt} \) holds the value that should have been obtained at the end of that integration step \( t + \Delta t \). The formula is almost the same as (2) of the last tutorial:
\begin{equation}
y(t + \Delta t) = y(t) + \frac{dy}{dt}|_{t+\Delta t} \Delta t = y(t) + f(t+\Delta t,y(t+\Delta t))\Delta t \tag{2}
\end{equation}
The only difference is that \( \frac{dy}{dt}|_t \) and \( f(t,y(t))\) have been replaced by \( \frac{dy}{dt}|_{t+\Delta t} \) and \( f(t+\Delta t,y(t+\Delta t)) \). Hence, it bears the name “backward” as it borrows future information for current inference. Then, as the unknown \( y(t+\Delta t) \) appears on both sides of (2), we will need to solve this algebraic equation to retrieve it. The backward Euler method is also a first-order method and has a truncation error proportional to the square of the time step size, \( (\Delta t)^2 \).
For a system of first-order ODEs, (2) now becomes
\begin{equation}
\textbf{y}(t + \Delta t) = \textbf{y}(t) + \frac{d\textbf{y}}{dt}|_{t+\Delta t} \Delta t = \textbf{y}(t) + \textbf{f}(t+\Delta t,\textbf{y}(t+\Delta t))\Delta t \tag{3}
\end{equation}
and we will have to solve for a system of algebraic equations. (see exercise)
Example
Numerically integrate the prototype equation
\begin{equation}
\frac{dy}{dt} = -\alpha y \tag{4}
\end{equation}
again with \( \alpha = 0.2 \) and \( y(0) = 1 \) as in the last tutorial, but with the backward Euler method now.
(2) in this case becomes
\begin{align}
y(t + \Delta t) &= y(t) -\alpha y(t + \Delta t) \Delta t \\
[y(t + \Delta t)] (1+\alpha \Delta t) &= y(t) \\
y(t + \Delta t) &= \frac{y(t)}{1+\alpha \Delta t} \tag{5}
\end{align}
If \( \Delta t = 0.1 \), then
\begin{align}
y(t + \Delta t) = \frac{y(t)}{1+(0.2)(0.1)} = \frac{y(t)}{1.02} \tag{6}
\end{align}
Hence at the \(n\)-th time step, the output will be
\begin{equation}
y(n \Delta t) = \frac{y(t)}{1.02^n} \tag{7}
\end{equation}
When \( t = 1.0 \), \( n = 10 \), we have \( y = 0.8203 \). The exponential decay simulated by the backward Euler method is slower than the true solution.
Stability
The backward Euler method is unconditionally stable. To see this, again we let \( y = A \lambda^n \). Substituting this into (5), we have
\begin{align}
A \lambda^{n+1} &= \frac{A \lambda^n}{1+\alpha \Delta t} \\
\lambda &= \frac{1}{1+\alpha \Delta t} \tag{8}
\end{align}
So \( 0 < \lambda < 1 \). No matter what time step size we choose, the solution will not blow up. The figure below shows the same plot as in the last tutorial under \( \alpha = 25\) but with the backward Euler scheme.
One may wonder, since the backward Euler method, or any implicit scheme in general, is more numerically stable and resistant to stiff problems, then why do we still sometimes choose to employ an explicit scheme? The reason is that inverting an algebraic equation or system like (3) in a backward fashion is computationally more difficult than direct computations, and frequently relies on other numerical techniques for that.
Exercise
Numerically integrate the exercise problem in the last tutorial using the backward Euler method with the same settings.
Answer
(3) for this scenario gives
\begin{align}
x(t + \Delta t) = x(t) + (0.5 x(t + \Delta t) -0.2 y(t + \Delta t) )\Delta t \\
(1-0.5\Delta t)x(t + \Delta t) + 0.2\Delta t y(t + \Delta t) = x(t)
\end{align}
and
\begin{align}
y(t + \Delta t) = y(t) + (0.1 x(t + \Delta t) +0.3 y(t + \Delta t) )\Delta t \\
-0.1\Delta t x(t + \Delta t) + (1-0.3\Delta t)y(t + \Delta t) = y(t)
\end{align}
With \( \Delta t = 0.2 \), the system of equations read
\begin{align}
\begin{aligned}
0.9x(t + \Delta t) + 0.04 y(t + \Delta t) &= x(t) \\
-0.02 x(t + \Delta t) + 0.94 y(t + \Delta t) &= x(t) \\
\end{aligned}
\end{align}
or in matrix form
\begin{align}
\begin{bmatrix}
0.9 & 0.04 \\
-0.02 & 0.94
\end{bmatrix}
\begin{bmatrix}
x(t + \Delta t) \\
y(t + \Delta t)
\end{bmatrix} &=
\begin{bmatrix}
x(t) \\
y(t)
\end{bmatrix} \\
\begin{bmatrix}
x(t + \Delta t) \\
y(t + \Delta t)
\end{bmatrix}
&=
\begin{bmatrix}
0.9 & 0.04 \\
-0.02 & 0.94
\end{bmatrix}^{-1}
\begin{bmatrix}
x(t) \\
y(t)
\end{bmatrix}
\end{align}
Hence
\begin{equation}
\begin{bmatrix}
x(n\Delta t) \\
y(n\Delta t)
\end{bmatrix}
= Q^n
\begin{bmatrix}
1 \\
0
\end{bmatrix}
\end{equation}
where
\begin{equation}
Q = \begin{bmatrix}
0.9 & 0.04 \\
-0.02 & 0.94
\end{bmatrix}^{-1}
\end{equation}
At \( t = 1.0 \), \( n = 5 \), it should give \( x = 1.67090, y = 0.16437 \) approximately.
Leave a Reply