import matplotlib
if not hasattr(matplotlib.RcParams, "_get"):
matplotlib.RcParams._get = dict.get
1.2. Initial Value Problem for ODE: single-step methods#
Let’s look back at what we have achieved in the first part of this chapter:
We can use finite differences to numerically approximate derivatives
We can derive these finite difference approximations from a Taylor series expansion, giving us an idea about the error of the approximation.
Next, let us see how we can apply these concepts in order to numerically integrate (that is, solve) ODEs.
The solution process of an ODE is performed by steps. In a single-step approach, the solution of the following step depends only on the current one. In a multiple step method, the solution of the next step is calculated from several steps. A multiple step approach is more accurate, you can think of a similar construct as the higher-order derivatives treated previously.
Consider a first-order ODE with the general form:
One initial condition is needed to find only one solution. Without it, there would be an infinite number of solutions possible. The initial condition is:
Forward (Explicit) Euler#
The simplest method for integrating initial value problems is the forward Euler scheme. It can be derived from the forward finite difference approximation of the derivative. Recall the forward finite difference formula:
where
For simplicity, we can also define \(y_{i} = y(x_{i})\) and write:
Now we can solve this expression for the state at the new time step, \(y_{i+1}\):
Ignoring the error terms, the forward Euler formula becomes:
where we substituted \(y'(x_i) = f(x_i, y_i)\). Note that the derivative of \(y\) (the slope) is evaluated at the current time step \(x_i\). This makes the forward Euler scheme an explicit method. Although simple, it contains the basic characteristics as more advanced and accurate methods.
The solution starts at \(i=0\) given by the initial condition, then \(i\) is increased to 1 where the values are calculated using the previous equations. This loop continues until the points cover the desired domain.
Fig. 1.4 Illustrating the forward Euler method#
Example#
Consider the following equation:
We want to know the solution in the domain \([a=0,b=30]\). The initial condition is \(y(x_0)=1\) and the step size is \(\Delta x = 0.2\).
The discretization of the differential equation transforms the problem into an algebraic one. Following the formula of explicit Euler for the discretization yields:
The following code snippet implements the forward Euler method for this example.
import numpy as np
import matplotlib.pyplot as plt
dx = 0.2
x = np.arange(0, 30 + dx, dx)
y = np.zeros(len(x))
##-----------------------------
##Forward Euler Implementation
##-----------------------------
alpha = 1
y[0] = 1
for i in range(len(x) - 1):
y[i + 1] = y[i] + dx * (-alpha * y[i])
##------------------------------
y_exact = np.exp(-alpha * x)
plt.plot(x, y_exact)
plt.plot(x, y)
plt.legend(["exact solution", "Forward Euler solution"])
Use the forward Euler to approximate the solution to the initial-value problem.
with \(\Delta t=0.5\)
Solution
\(t_0= 0, t_1= 0.5, t_2= 1\)
Use the forward Euler:
for i =0
for i = 1
The approximated solution is 2.625
Backward (Implicit) Euler#
Just as the forward Euler method is the simplest explicit method, the backward Euler method is the simplest implicit method. Both methods look very similar, but the key difference is the point where the slope, given by \(f(x, y)\), is evaluated. Remember that in the explicit Euler scheme, we evaluate the slope at the current value \(x_i\). In contrast, in the backward Euler scheme, the slope is computed at the next step, \(x_{i+1}\).
That is, the slope depends on the unknown value \(y_{i+1}\). This is why the backward Euler scheme is an implicit method.
Derive the backward Euler formula from the backward finite difference approximation for the first derivative (presented in Section II.ii).
Solution
The backward finite difference approximation for the first derivative of a function \(y\) at point \(x_{i+1}\) is:
Since \(x_{i+1} - \Delta x = x_i\), we can write:
Multiplying by \(\Delta x\) yields:
Now, we solve for \(y(x_{i+1})\):
Defining \(y_i = y(x_i)\) and \(y_{i+1} = y(x_{i+1})\), and considering that \(y'(x_{i+1}) = f(x_{i+1}, y_{i+1})\) we have:
Let’s consider the same problem as before:
We want to know the solution in the domain \([a=0,b=30]\). The initial condition is \(y(x_0)=1\) and the step size is \(\Delta x = 0.2\).
The discretization of the differential equation transforms the problem into an algebraic one. Plugging the ODE into the formula for the implicit Euler scheme yields:
We bring all unknowns to the left side:
Finally, we solve for \(y_{i+1}\):
In the following code, Implicit Euler is implemented.
dx = 0.3
x = np.arange(0, 30 + dx, dx)
y = np.zeros(len(x))
##-----------------------------
##Backward Euler Implementation
##-----------------------------
alpha = 1
y[0] = 1
for i in range(len(x) - 1):
y[i + 1] = y[i] / (1 + dx * alpha)
##------------------------------
y_exact = np.exp(-alpha * x)
plt.plot(x, y_exact)
plt.plot(x, y)
plt.legend(["exact solution", "Backward Euler solution"])
You can see that the result is similar to the forward Euler solution, except that the curve for the backward Euler scheme is slightly above the exact solution. This makes sense as the derivative is taken at the end of the integration interval. For exponential decay, the slope decreases throughout the integration interval but we assume that it is constant. If we evaluate it at the end of the integration interval (backward Euler), we underestimate the slope. In contrast, the explicit (forward Euler) solution overestimates the slope.
Attribution
This chapter is written by Jaime Arriaga Garcia, Anna Störiko, Justin Pittman and Robert Lanzafame. Find out more here.