2.7. Implicit methods for nonlinear ODE#
In the previous section, the explicit and implicit implementation seemed to only differ in the formulation. However, this changes when treating non-linear equations. Let’s see an example to identify the “implicitness”.
An example: mixed bioreactor#
Consider a bioreactor that is continuously stirred and where microorganisms degrade a chemical compound (see Fig. 2.12). A solution containing the compound at concentration \(c_{\mathrm{in}}\) is fed into the reactor with flow rate \(Q\). The same amount of water is removed from the reactor at the outflow. The concentration inside the reactor is designated by \(c\), and the reactor volume by \(V\).
Fig. 2.12 Conceptual model of a mixed bioreactor#
The reaction rate \(r\) with which the bacteria degrade the compound depends non-linearly on the concentration: at low concentrations the increase is linear, but when the concentration increases, the rate levels off and converges towards a constant maximum rate, \(r_{\mathrm{max}}\). This can be described by the so called Monod kinetics:
where the parameter \(K\) indicates the concentration where the rate reaches half its maximum value.
Based on a mass balance, we can set up the following ODE to describe how the concentration in the reactor changes over time:
The initial concentration in the reactor is \(c(t=0) = c_0\). To predict how the concentration in the reactor evolves over time, we need to solve the ODE. Let’s see how this looks like for an implicit and explicit Euler scheme.
The expression for the explicit Euler scheme is:
The algorithm to solve the explicit Forward Euler is exactly the same as before.
import numpy as np
import matplotlib.pyplot as plt
# set up time vector
dt = 0.2
t_end = 100
n_times = int(t_end / dt)
t = np.linspace(0, t_end, num=n_times)
# parameter values
c0 = 80.0
K = 50.0
c_in = 50.0
r_max = 5.0
Q = 0.1
V = 2.0
##-----------------------------
##Explicit Euler Implementation
##-----------------------------
c = [c0]
for i in range(n_times - 1):
c_i = c[-1]
dcdt = Q / V * (c_in - c_i) - r_max * c_i / (c_i + K)
c.append(c_i + dt * dcdt)
##------------------------------
plt.plot(t, c)
plt.legend(["Forward Euler solution"])
plt.xlabel("time")
plt.ylabel("concentration")
Now, using the implicit Euler scheme, the expression becomes:
Even though this equation looks similar to the explicit Euler scheme at a first glance, it requires a different procedure because now the unknown \(c_{i+1}\) appears on both sides of the equation. Since the equation is non-linear, we cannot solve directly for \(c_{i+1}\). Instead, an iterative procedure is needed, for example, the widely used Newton-Raphson method.
The Newton-Raphson Method#
The Newton-Raphson method is a numerical technique to find the roots of a function, that is, the locations where its value equals zero. How does that help us to solve our implicit Euler expression for \(c_{i+1}\)? Note that we can easily rewrite the equation such that the right-hand-side is zero, by moving all terms from the right to the left:
That is, the expression on the left is the function that we want to find the root of – let’s call it \(q(c_{i+1})\), and the root will equal the new concentration \(c_{i+1}\) that we are looking for. So let’s see how the Newton-Raphson method works.
We consider a function \(g(z)\), and we are looking for the value \(z\) where \(g(z) = 0\). The solution is graphically represented by the intersection of \(g(z)\) with the \(z\) axis (see the Figure below). The key idea of the Newton-Raphson method is to linearize the root finding problem. That is, we find the root of a tangent to the function \(g(z)\) at a point \(z_0\). Since the tangent is a simple line, it is easy to find its intersection with the \(y\)-axis. This will not directly bring us to the correct solution for our original problem (finding the root of \(g(z)\)), but if we repeat the approach, it will iteratively bring us closer. This is illustrated by Fig. 2.13.
Fig. 2.13 Illustration of the Newton-Raphson method#
Step by step, the procedure for the Newton-Raphson method is as follows:
We start with an initial guess \(z_0\) and compute the tangent of the function \(g(z)\) at this point. The tangent \(t(z)\) has a slope of \(g^\prime(z)\) and is described by \(t(z) = g(z_0) + (z-z_0) \,g^\prime(z_0)\).
We follow the tangent to the point where it crosses the \(z\)-axis. Setting \(t(z) = 0\), we can easily see that this is at \(z_1 = z_0 - \frac{g(z_0)}{g^\prime(z_0)}\). Note how this point is already much closer to the root of \(g(z)\).
Now, we repeat the procedure. That is, we compute the tangent at point \(z_1\) and find its root to get to a new point \(z_2\). Each new point \(z_{j+1}\) can be obtained from the previous one as:
\[z_{j+1} = z_j - \frac{g(z_j)}{g^\prime(z_j)}\]We stop the process when the value of \(g(z_j)\) becomes very small: \(|g(z_j)<\epsilon|\). Here, \(\epsilon\) is the tolerance e.g., \(\epsilon=10^{-6}\) (see the red square in Fig. 2.13).
Application to the implicit Euler scheme#
We now apply this procedure to the expression we obtained from the implicit Euler scheme, to find the new concentration \(c_{i+1}\). Remember the function \(q(c_{i+1})\) that we defined earlier:
More generally, we can also write this for any ODE of the form \(\frac{dc}{dt} = f(c, t)\) as:
Since \(q(c_{i+1})\) corresponds to the function \(g(z)\), we can iteratively find new values for \(c_{i+1}^{j+1}\) with the following equation:
Note that the index \(i\) refers to the time step of the implicit Euler method, whereas the index \(j\) refers to the iteration of the Newton-Raphson scheme. To compute the value \(c_{i+1}^{j+1}\), we first need to evaluate \(q^\prime\), the derivative of \(q\) with respect to \(c_{i+1}\).
where \(f^\prime = \frac{\partial f(c, t)}{\partial c}\).
For the ODE describing the bioreactor, this is:
The final expression for \(c_{i+1}^{j+1}\) then becomes:
And for the bioreactor:
This expression might look a bit intimidating when written out fully, but remember that when implementing the procedure in code, you can break it down into several steps. For example, define expressions (or functions) for \(q\) and \(q^\prime\) first and use them for defining \(c_{i+1}^{j+1}\). In practice, the implementation of the implicit Euler method with Newton-Raphson iterations in Python will require two nested loops: an outer loop for time stepping, and an inner loop for the Newton-Raphson iterations.
Advantages and disadvantages of the Newton-Raphson method#
The Newton-Raphson method
converges quickly and can provide highly accurate approximations,
can be used for both real and complex roots and
can be extended to multi-dimensional problems (e.g., a system of coupled ODE).
But the method also has its limitations:
It may not converge if the initial guess is far from the actual value or if the function is discontinuous near the value.
The derivative of the function is needed, which may not always be known.
If no closed form expression is available for the derivative, one alternative is to approximate it with finite differences. That is, we use a small perturbation of the state and evaluate the function twice, once for each value. This approach is shown in the introductory example on numerical solutions in Section 2.1 (see under Extra: More efficient numerical solution approach). The drawback of approximating the derivative numerically is the added computational cost.
For multi-dimensional problems, the Newton-Raphson method requires inverting the Jacobian, which is computationally expensive. The Jacobian is the matrix of partial derivatives for all equations with respect to all states.
Attribution
This chapter is written by Jaime Arriaga Garcia, Anna Störiko, Justin Pittman and Robert Lanzafame. Find out more here.