Numerical Integration by Forward Euler Method
There are a lot more ODEs outside of what we have discussed in the previous chapters that cannot be solved analytically than we may have imagined. In fact, if we pick an ODE at random, it will have a zero probability of possessing a closed-form solution. In practice, if we want to extract some useful information about the solution of such an ODE, the only hope is to do so via numerical integration methods on a computer, which approximate the solution values at discrete sampling points.
The simplest numerical integration method is probably the forward Euler method. Take a general first-order ODE as an example:
\begin{equation}
\frac{dy}{dt} = f(t,y) \tag{1}
\end{equation}
We set out to compute the estimated values of \( y \) at an even time step \( \Delta t \), e.g. if \( \Delta t = 0.1 \), then we iterate over \( t = 0, \Delta t, 2\Delta t, \ldots = 0, 0.1, 0.2, \ldots \). The forward Euler method assumes that \( \frac{dy}{dt} \) in (1) remains the same throughout a single time step \( (t, t+\Delta t) \), taking the value at the beginning of that integration step, \( t \). Hence, we can approximate the value of \( y \) at the next time step by a linear extrapolation:
\begin{equation}
y(t + \Delta t) = y(t) + \frac{dy}{dt}|_t \Delta t = y(t) + f(t,y(t))\Delta t \tag{2}
\end{equation}
Then, we repeatedly apply this formula with \( t \to t + \Delta t \) to move forward in time, computing the value of \(y\) at the next time step from the last one until we reach the end point. The forward Euler method is a first-order method, which has a truncation error proportional to the square of the time step size, \( (\Delta t)^2 \).
It can be readily generalized to a system of first-order ODEs, where (2) simply becomes
\begin{equation}
\textbf{y}(t + \Delta t) = \textbf{y}(t) + \frac{d\textbf{y}}{dt}|_t \Delta t = \textbf{y}(t) + \textbf{f}(t,\textbf{y}(t))\Delta t \tag{3}
\end{equation}
by replacing the scalars \( y \) and \( f \) by vector-valued versions \( \textbf{y}, \textbf{t} \). (see exercise below)
Example
Numerically integrate
\begin{equation}
\frac{dy}{dt} = -\alpha y \tag{4}
\end{equation}
where \( \alpha = 0.2 \) and \( y(0) = 1 \), with the forward Euler method at \( \Delta t = 0.1 \) until \( t = 1 \).
Below shows a table listing out the values of \( \frac{dy}{dt} \) and \( y \) computed according to (2), at each time step.
| \(t\) | \( y \) | \( \frac{dy}{dt} \) |
| \( 0.0 \) | \( 1 \) | \( -0.2(1) = -0.2 \) |
| \( 0.1 \) | \( 1 + (-0.2)(0.1) = 0.98 \) | \( -0.2(0.98) = -0.196 \) |
| \( 0.2 \) | \( 0.98 + (-0.196)(0.1) = 0.9604 \) | \( -0.2(0.9604) = -0.1928 \) |
| \( 0.3 \) | \( 0.9604 + (-0.1928)(0.1) = 0.94119 \) | \( -0.2(0.94119) = -0.18824 \) |
| \( 0.4 \) | \( 0.94119 + (-0.18824)(0.1) = 0.92237 \) | \( -0.2(0.92237) = -0.18447 \) |
| \( 0.5 \) | \( 0.92237 + (-0.18447)(0.1) = 0.90392 \) | \( -0.2(0.90392) = -0.18078 \) |
| \( 0.6 \) | \( 0.90392 + (-0.18078)(0.1) = 0.88584 \) | \( -0.2(0.88584) = -0.17717 \) |
| \( 0.7 \) | \( 0.88584 + (-0.17717)(0.1) = 0.86813 \) | \( -0.2(0.86813) = -0.17363 \) |
| \( 0.8 \) | \( 0.86813 + (-0.17363)(0.1) = 0.85076 \) | \( -0.2(0.85076) = -0.17015 \) |
| \( 0.9 \) | \( 0.85076 + (-0.17015)(0.1) = 0.83375 \) | \( -0.2(0.83375) = -0.16675 \) |
| \( 1.0 \) | \(0.83375 + (-0.16675)(0.1) = 0.81707 \) | \( -0.2(0.81707) = -0.16341 \) |
Compared to the analytical solution \( e^{-\alpha t} \), where the true value will be \( e^{-0.2} = 0.81873 \) at \( t = 1.0 \), the numerical output \( 0.81707 \) is quite close, albeit with a bit of underestimation. This is because \( dy/dt = -\alpha y = -e^{-\alpha t} \) is an exponentially decaying function, so by using a fixed value of \( dy/dt \) from the beginning of an integration step, the numerical decrease in \(y\) always overshoots.
Stiffness Problem
Now let’s consider (4) again, but with \( \alpha \) set to \( 25 \). If we keep the same time step size \( \Delta t = 0.1 \), then it is easy to see that the forward Euler scheme will give, at the first few time steps
| \(t\) | \( y \) | \( \frac{dy}{dt} \) |
| \( 0.0 \) | \( 1 \) | \( -25(1) = -25 \) |
| \( 0.1 \) | \( 1 + (-25)(0.1) = -1.5 \) | \( -25(-1.5) = 37.5 \) |
| \( 0.2 \) | \( -1.5 + (37.5)(0.1) = 2.25 \) | \( -25(2.25) = -56.25 \) |
| \( 0.3 \) | \( 2.25 + (-56.25)(0.1) = -3.375 \) | \( -25(-3.375) = 84.375 \) |
and the pattern will keep going on, the sign of \(y\) oscillates while the magnitude grows by a factor of \(1.5\) at each iteration, indefinitely. The solution will hence blow up to infinity (while the actual solution should decay), and the forward Euler scheme is numerically unstable in this set-up. The problem (4) with such a value of \( \alpha \) is called stiff because it now requires a very small time step size, e.g. \( \Delta t = 0.01 \) to remain stable.
In fact, to inquire the criterion for stability, we may let \( y = A \lambda^n \) where \( n \) is the time index and plug this into (2) to get
\begin{align}
A \lambda^{n+1} &= A \lambda^n -\alpha A \lambda^n \Delta t \\
\lambda &= 1 -\alpha \Delta t \tag{5}
\end{align}
We then demand \( |\lambda| < 1 \), and this shows that
\begin{align}
0 &< \alpha \Delta t < 2 \\
0 &< \Delta t < \frac{2}{\alpha} \tag{6}
\end{align}
the time step must be smaller than two times the reciprocal of \( \alpha \). The figure below shows how the numerical integration result varies when different values of \( \Delta t \) are used when \( \alpha = 25\).
Exercise
Numerically integrate
\begin{align}
\left\{\begin{aligned}
\frac{dx}{dt} &= 0.5 x -0.2 y \\
\frac{dy}{dt} &= 0.1 x + 0.3 y
\end{aligned}\right. \tag{7}
\end{align}
where the I.C.s are \( x(0) = 1, y(0) = 0 \), using the forward Euler method with a time step of \( \Delta t = 0.2 \).
Answer
The integration result at the first few time steps is shown in the table below.
| \(t\) | \(x\) | \( \frac{dx}{dt} \) | \( y \) | \( \frac{dy}{dt} \) |
| \( 0.0 \) | \( 1 \) | \( 0.5 (1) -0.2 (0) = 0.5\) | \( 0 \) | \( 0.1(1) + 0.3(0) = 0.1 \) |
| \( 0.2 \) | \( 1 + (0.5)(0.2) = 1.1 \) | \( 0.5 (1.1) -0.2 (0.02) = 0.546\) | \( 0 + (0.1)(0.2) = 0.02 \) | \( 0.1(1.1) + 0.3(0.02) = 0.116 \) |
| \( 0.4 \) | \( 1.1 + (0.546)(0.2) = 1.2092\) | \( 0.5 (1.2092) -0.2 (0.0432) = 0.5960 \) | \( 0.02 + (0.116)(0.2) = 0.0432 \) | \( 0.1(1.2092) + 0.3(0.0432) = 0.1339 \) |
| \( 0.6 \) | \( 1.2092 + (0.5960)(0.2) = 1.3284 \) | \( 0.5 (1.3284) -0.2 (0.0700) = 0.6502 \) | \( 0.0432 + (0.1339)(0.2) = 0.0700 \) | \( 0.1(1.3284) + 0.3(0.0700) = 0.1538 \) |
| \( 0.8 \) | \( 1.3284 + (0.6502)(0.2) = 1.4584 \) | \( 0.5 (1.4584) -0.2 (0.1007) = 0.7091 \) | \( 0.0700 + (0.1538)(0.2) = 0.1007 \) | \( 0.1(1.4584) + 0.3(0.1007) = 0.1761 \) |
| \( 1.0 \) | \( 1.4584 + (0.7091)(0.2) = 1.6002 \) | \( 0.5 (1.6002) -0.2 (0.1360) = 0.7729 \) | \( 0.1007 + (0.1538)(0.2) = 0.1360 \) | \( 0.1(1.6002) + 0.3(0.1360) = 0.2008 \) |
Leave a Reply