27.2.2 First-Order Methods: Euler’s Method

27.2.2.1 Introduction

The first method we will look into is the Euler method that was introduced by Euler in 1768 [1]. This method is one of the methods that are referred to as first-order methods because they use an expansion degree of n_max = 1 for Eq. 27.1. For the first step, we will only restrict ourselves to Eq. 27.1a. Expanding this equation to n_max = 1 yields

$\begin{array}{ll}f\left(x+\mathrm{\Delta }x\right)\hfill & =\frac{1}{0!}\frac{d^{0}f}{d{x}^0}{\left(\mathrm{\Delta }x\right)}^0+\frac{1}{1!}\frac{d{f}^1}{d{x}^1}{\left(\mathrm{\Delta }x\right)}^1+O_2\hfill \\ \hfill & =f\left(x\right)+\frac{df}{dx}\left(x\right)\mathrm{\Delta }x+O_2\hfill \\ \hfill & \approx f\left(x\right)+\frac{df}{dx}\left(x\right)\mathrm{\Delta }x\hfill \end{array}$

As you can see, this approximation states that the value of the function f (x + Δx) is equal to the value of the function f at x plus the slope of the function in x times the distance Δx. If you look at Fig. 27.1, you can see this is identical to putting the tangent to the function in x₀ and extrapolating for a distance Δx. As you can also see, this approximation makes little mistakes if the curve is smooth and the values of Δx are small (see Fig. 27.1a). For more strongly curved functions f and bigger values for Δx, the approximation becomes inexact (see Fig. 27.1b). In this example, the algorithm would predict a lower value of f (x₁) than actually found.

The Euler method will allow us predicting the (unknown) value of a function f (x₁) at the location x₁, starting from an initial (known or assumed) value f (x₀) at x₀ simply by using the gradient or slope of the function at $\frac{df}{dx}\left(x_0\right)$ at x₀. Once this value has been found, we are able to calculate the unknown value of the function f (x₂) at x₂ from the value of f (x₁) and so forth. By repeating this scheme, we are able to derive the entire function, i.e., the values of the function f (x_i) for different discreet values x_i.

27.2.2.2 Using Euler’s Method to Approximate a Known Function

In order to get a better feeling for what the Euler method actually does, we will look at an example which is not yet an ODE. However, it will allow understanding of what exactly Eq. 27.2 does and what the individual parameters influencing the exactness of the solution actually are. For this we will assume an arbitrary function

$f\left(x\right)=\frac{1}{4}{x}^4-\frac{1}{2}{x}^3-\frac{3}{2}{x}^2$

For which we find the first derivative as

$\frac{df}{dx}=x^3-\frac{3}{2}{x}^2-3x$

According to Eq. 27.2 we should be able to reconstruct discreet values of this function using an iterative scheme starting from x₀ = 0 and f (0) = 0. The derivative at x₀ = 0 is $\frac{df}{dx}\left(0\right)=0$. We will assume a step width of Δx = 1. As we will find out in a moment, this value is significantly too big. Applying Eq. 27.2 we find the next value

$f\left(0+\mathrm{\Delta }x\right)=f\left(1\right)=f\left(0\right)+\frac{df}{dx}\left(0\right)\mathrm{\Delta }x=0$

If we look at Eq. 27.3 we can see that the value should be $f\left(1\right)=-\frac{7}{4}=-1.75$. If we continue with the iteration, we find the values diverging further (see Tab. 27.1). If we decrease Δx to a value of 0.4, we start getting significantly better approximation of the function (see Tab. 27.1).

Fig. 27.2 shows the influence Δx has on the exactness of the solution. For Δx = 1 we already found the function to be only poorly approximated, which we can confirm from Fig. 27.2a. Decreasing Δx to Δx = 0.4 already catches on significantly better to the function. Starting from Δx = 0.2 we get a close-to-correct approximation and at Δx = 0.08 there is hardly any error anymore.

This serves to illustrate that the Euler method is very sensitive to values of Δx which are too big. The higher the required accuracy the smaller the value of Δx should be.

Comparison to Higher-Order Methods. Compared to first-order methods, higher-order methods use higher expansion orders, thus we would find a more complicated expression as Eq. 27.2. However, these methods will provide more stable results at larger values of Δx. In general, this is the trade-off: If taking into account higherorder terms of Eq. 27.1a, the solution is more stable and can tolerate larger values for Δx.

Backwards Finite Distance. As we have worked with Eq. 27.1a, the equations we derived are referred to as the forward finite distance approximation of f (x). Obviously, we could just as well have worked with Eq. 27.1b in which case we would have found

$\begin{array}{ll}f\left(x-\mathrm{\Delta }x\right)\hfill & ={\left(-1\right)}^n\frac{1}{0!}\frac{d^{0}f}{d{x}^0}{\left(\mathrm{\Delta }x\right)}^0+\frac{1}{1!}\frac{d{f}^1}{d{x}^1}{\left(\mathrm{\Delta }x\right)}^1+O_2\hfill \\ \hfill & =f\left(x\right)-\frac{df}{dx}\left(x\right)\mathrm{\Delta }x+O_2\hfill \\ \hfill & \approx f\left(x\right)-\frac{df}{dx}\left(x\right)\mathrm{\Delta }x\hfill \end{array}$

instead of Eq. 27.2. Eq. 27.5 is referred to as the backwards finite distance approximation of f (x).

27.2.2.3 Using Euler’s Method to Solve an ODE

We will now use Euler’s method to actually solve an ODE. For this we turn to the following example differential equation that we obtained when deriving the discharging characteristics of a capacitor

$\begin{array}{ll}Q\hfill & =-CU\hfill \\ \hfill & =-CRI\hfill \\ \hfill & =-CR\frac{dQ}{dt}\hfill \\ \frac{dQ}{dt}\hfill & =-\frac{1}{CR}Q\hfill \end{array}$

where C is the capacitance, U the voltage, Q the charge, and R the Ohmic resistance. For the sake of simplicity, we will combine the constants to a fixed value and assume the capacitor was full at the beginning of the experiment, i.e., Q(t = 0) = 1. The equation can thus be rewritten to

$\frac{dQ}{dt}=-\frac{1}{2}Q,Q\left(t=0\right)=1$

Obviously, this equation can be solved immediately by separation of variables and partial integration (see section 8.2.2) yielding

$Q\left(t\right)={\mathrm{e}}^{-\frac{t}{2}}$

However, we will now derive this solution numerically.

Iterative Scheme. For this we will employ the iterative scheme given by Eq. 27.2. We know the first value Q(0) = 1 that was given as initial value. With this value, we can calculate the derivative at t = 0 using Eq. 27.6 to be

$\frac{dQ}{dt}\left(t=0\right)=-\frac{1}{2}$

Employing Eq. 27.2 we can then calculate the next value Q(t + Δt). We will assume a value of Δt = 0.1 in which case we find

$\begin{array}{ll}Q\left(\mathrm{\Delta }t\right)\hfill & =Q\left(0\right)+\frac{dQ}{dt}\left(0\right)0.1\hfill \\ \hfill & =1-\frac{1}{2}0.1=0.95\hfill \end{array}$

We would then proceed to insert the calculated derivative at 2Δt using Eq. 27.6 as

$\begin{array}{ll}\frac{dQ}{dt}\left(\mathrm{\Delta }t\right)\hfill & =-\frac{1}{2}Q\left(\mathrm{\Delta }t\right)\hfill \\ \hfill & =-\frac{1}{2}0.95\hfill \\ \hfill & =-0.475\hfill \end{array}$

whereupon we could calculate the value as

$\begin{array}{ll}Q\left(2\mathrm{\Delta }t\right)\hfill & =Q\left(\mathrm{\Delta }t\right)+\frac{dQ}{dt}\left(\mathrm{\Delta }t\right)0.1\hfill \\ \hfill & =0.95-0.4750.1=0.9025\hfill \end{array}$

We would proceed with this scheme until we have found sufficient discreet values t_i of Q(t) to accept as a solution.

Maple Worksheet. Obviously, at this point it makes sense to turn to Maple and implement Euler’s method. Listing 27.1 shows the implementation. We will briefly walk through the code. As usual, we first restart the Maple server and read Core.txt (see line 1). Next, we define the function DoEuler which implements Euler’s method. It accepts the differential equation df and the numerical range of the solution in the parameters _from (lower bound) and _to (upper bound). It also expects the initial value of the function at the lower bound given as parameter initialvalue. The number of points is given as points. It is important to note that the parameter df requires a function both of the value of the function Q and of the independent variable t, as both may be required to evaluate the derivative. Using our example equation (Eq. 27.6), we actually would not require the independent variable for the right-hand side of the equation. However, it is required in general.

1  restart;read("Core.txt"):
2
3  DoEuler:=proc(df,_from,_to,initialvalue,points) local ret,delta,i:
4    ret:=Matrix(points,2):
5    delta:=(_to-_from)/(points-1):
6    ret[1,1]:=_from:ret[1,2]:=initialvalue:
7    for i from 2 to points do:
8     ret[i,1]:=evalf(ret[i−1,1]+delta):
9     ret[i,2]:=evalf(ret[i−1,2]+delta*df(ret[i−1,1],ret[i−1,2])):
10    end do:
11  return ret;
12 end proc:
13
14 display(quickplot([t->exp(-t/2)],0..4,"x","f(x)"),pointplot(DoEuler((t,Q)->-1/2*Q,0,4,1,41)));

The function first prepares the return value ret that is a points × 2 matrix (see line 4). It then calculates delta and sets the initial values as the first entries in the matrix (see line 6). It then iterates starting from the second line in the matrix (thus the starting value of 2 for i), as we already inserted the initial conditions at the first line. Each iteration calculates another data point.

We now want to use this code for calculating values of t_i (see line 8) and Q_i (see line 9). The latter equation uses the derivative as a function of Q that must be given as the parameter df. As you can see from line 14, we supply Eq. 27.6 for this function. If this equation is evaluated with a given value Q_i at t_i it allows calculating $\frac{dQ}{dt}\left(t_i\right)$. This procedure is repeated for a total of points iterations.

In line 14 we call DoEuler creating a point plot from the resulting matrix. As stated, we supply Eq. 27.6 for df and evaluate the function in the range 0 ≤ t ≤ 4. The initial value for t = 0 is given as 1 and we request a total of 41 points that amounts to a Δt = 0.1 as used in the first step of our manual calculation. We combine this pointplot with the plot of the analytical solution function given by Eq. 27.7 in order to compare the output. The result of this line is shown in Fig. 27.3. The first 10 data points are also shown in Tab. 27.2. As you can see, the numerical solution is pretty close to the analytical solution.

Standard Notation. At this point we will quickly introduce a standard notation for iterative schemes, which makes the formulas a little more compact. In this notation, the values of the function f (x) at point x_i are abbreviated as F⁽ⁱ⁾, whereas the derivative $\frac{df}{dx}\left(x\right)$ at point x_i is abbreviated X⁽ⁱ⁾. The stepwidth Δx is generally abbreviated as h. In this notation Eq. 27.2 would be written as

$F^{\left(i+1\right)}=F^{\left(i\right)}+h{X}^{\left(i\right)}$

Likewise, the equation for backward finite distance Eq. 27.5 would be written as

$F^{\left(i-1\right)}=F^{\left(i\right)}-h{X}^{\left(i\right)}$

27.2.2.4 Summary

As you can see, the Euler method can be used conveniently to solve ODEs. We have discussed a simple example, but more complicated ones would not pose a significant problem for the function DoEuler we implemented. However, as we will see in just a moment, the applicability of the Euler method is quite limited due to the fact that it tends to be numerically not very stable and therefore, often fails to find solutions even to pretty simple ODEs.

27.2.3 First-Order Methods: Backward Euler Method

27.2.3.1 Introduction

We have studied the Euler’s method as a very convenient tool for solving ODEs. However, the simplicity of the method comes at a disadvantage which is bad stability. Stability of a numerical method is strongly dependent on the magnitude of the error. During the iterative solving, the error function may grow substantially as a consequence of parasitic solutions, i.e., solutions which are not sought. We have experienced a similar problem during the study of the LORAN positioning system where we have found that our algorithm aborted for wrong initial values (see Tab. 26.2). This is a consequence of a parasitic solution becoming predominant over the actual solution we are looking for. For first-order methods, we usually try to suppress these solutions by using small stepwidth values and choosing suitable initial values. As a quick example, we will numerically solve the ODE

$\frac{df}{dt}=-af,f\left(0\right)=0$

which has the obvious solution obtained by separation of variables and partial integration (see section 8.2.2)

$f\left(x\right)=-{\mathrm{e}}^{-ax}$

We will apply the function DoEuler from listing 27.1 to solve this equation numerically for different values of a. Fig. 27.4 lists the output of the algorithm for different values of a. As you can see, the solution is relatively exact for a = 1, whereas the results for a = 2 are incorrect. For higher values of a, the algorithm becomes numerically unstable. This is a clear example of the numerical instability of this first-order method. As we will see in a moment, multistep or higher-order methods have smaller error functions, which is why they tend to be less likely to grow due to the parasitic solutions. However, we can apply a small trick to Euler’s method, increasing its stability significantly.

27.2.3.2 Derivation

The starting equation is Eq. 27.1b. We expand this series to n_max = 1 finding

$\begin{array}{ll}f\left(x-\mathrm{\Delta }x\right)\hfill & =f\left(x\right)-\frac{df}{dx}\left(x\right)\mathrm{\Delta }x+O_2\hfill \\ \hfill & \approx f\left(x\right)-\frac{df}{dx}\left(x\right)\mathrm{\Delta }x\hfill \end{array}$

Now if we evaluate Eq. 27.11 at the discreet value x_i we would find

$f\left(x_i-\mathrm{\Delta }x\right)=f\left(x_i\right)-\frac{df}{dx}\left(x_i\right)\left(\mathrm{\Delta }x\right)$

For the value x^(i+ 1) = x_i + ΔxEq. 27.11 would amount to

$f\left(x^{\left(i+1\right)}\right)=f\left(x_i+\mathrm{\Delta }x-\mathrm{\Delta }x\right)=f\left(x_i\right)=f\left(x_i+\mathrm{\Delta }x\right)-\frac{df}{dx}\left(x_i+\mathrm{\Delta }x\right)\mathrm{\Delta }x$

from which we can derive

$f\left(x^{\left(i+1\right)}\right)=f\left(x_i\right)+\frac{df}{dt}\left(x^{\left(i+1\right)}\right)\mathrm{\Delta }x$

where we have derived an implicit equation for the value of the function at the location x^(i+ 1). In standard notation, this equation would be written as

$F^{\left(i+1\right)}=F^{\left(i\right)}+h{X}^{\left(i+1\right)}$

This method is commonly referred to as backward Euler method or implicit Euler method. Even though it has the same order as the regular (often referred to as forward) Euler method, it is vastly superior in terms of stability. The equation is implicit, which means it requires the values of the next step in order to calculate the previous. The forward Euler method proceeds by calculating the next value from the previous. The backward Euler method is implicit because values at x^(i+ 1) appear on both sides of the equation. Thus for each step, one would need to solve an algebraic equation or potentially (if more than one unknown are used) a system of equations. Once this equation is solved, the value f (x^(i+ 1)) is found.

The obvious disadvantage of having to solve an algebraic equation (or a set of equations) comes with one major advantage: stability. The forward Euler method is significantly more stable and thus a very good candidate for ODEs that are troublesome to solve with less stable (but computationally less expensive) algorithms.

27.2.3.3 Applying the Backward Euler Method

Let us return to the example of Eq. 27.10 for which the forward Euler method has failed. We will quickly go through the first steps of solving this ODE using the backward Euler method for a = 10 and Δx = 0.5. We first begin by noting that our ODE is not very complicated and thus, we can write down

$\frac{df}{dt}\left(x^{\left(i+1\right)}\right)=-10f\left(x^{\left(i+1\right)}\right)$

directly from Eq. 27.10. We can therefore rewrite Eq. 27.8 to

$\begin{array}{ll}f\left(x^{\left(i+1\right)}\right)\hfill & =f\left(x_i\right)-10x^{\left(i+1\right)}f\left(x^{\left(i+1\right)}\right)\mathrm{\Delta }x\hfill \\ f\left(x^{\left(i+1\right)}\right)\hfill & =\frac{f\left(x^{\left(i\right)}\right)}{1+10\mathrm{\Delta }x}\hfill \end{array}$

where we have already obtained our iteration scheme as this solution is already solved for f (x_i + Δx). As stated, if the right-hand side of Eq. 27.10 was a bit more complex, we would not be able to directly derive a solution analytically, but would be required to solve the equation for f (x_i + Δx) in each step.

We begin with the first value x₀ = 0 for which the initial condition gives us f (x₀) = 1. From Eq. 27.15 we find for the point f (x_i + Δx) as

$\begin{array}{ll}f\left(0+\mathrm{\Delta }x\right)\hfill & =\frac{f\left(x_0\right)}{1+10\mathrm{\Delta }x}\hfill \\ f\left(0.5\right)\hfill & =\frac{1}{1+100.5}=0.166\hfill \end{array}$

The fact that we are able to write this out so conveniently is by accident. In a moment, we will look at an example of a function where this is not the case, and we need to solve an equation for each iteration.

27.2.3.4 Maple Worksheet

At this point, we will derive a listing for the backward Euler method. The code is shown in listing 27.2. For the most part, this listing is analogous to the forward Euler method (see listing 27.1). The main difference is in line 9 where we inserted a try/catch block. In line 10 we encounter Eq. 27.8 in which f (x^(i+ 1)), stored as the variable Fi, is set equal to the last value f (x⁽ⁱ⁾) plus the evaluated function df at x^(i+ 1). As in listing 27.2, df must be given as a function of the independent variable x and of the function value at x^(i+ 1) (which is f(x^(i+ 1)) or Fi). This will give rise to an algebraic equation that is solved in the next line. In case this solving process is not successful, the algorithm aborts printing out the equation that could not be solved. In line 20 the example function is solved and plotted for a = 10. This is one of the plots shown in Fig. 27.5 where the function was applied to the four equations for which the forward Euler method failed in Fig. 27.4. As you can see, in all of the four cases, the backward Euler method shows stable solutions. This is where the superior stability compared to the forward Euler method becomes obvious. The algorithm is still of first order which means that the solution will not be exceedingly good for large stepwidths, but we will obtain a stable solution in almost all cases.

1 restart;read("Core.txt"):
2
3 DoBackwardsEuler:=proc(df,_from,_to,initialvalue,points) local ret,delta,i,Fi,tosolve:
4  ret:=Matrix(points,2):
5  delta:=(_to-_from)/(points-1):
6  ret[1,1]:=_from:ret[1,2]:=initialvalue:
7  for i from 2 to points do:
8    ret[i,1]:=evalf(ret[i-1,1]+delta):
9    try
10     tosolve:=Fi=ret[i-1,2]+df(ret[i,1],Fi)*delta:
11     ret[i,2]:=rhs(solve({tosolve},{Fi})[1]):
12   catch:
13      print(sprintf("Cannot solve %a",tosolve));
14      return;
15    end try:
16   end do:
17  return ret;
18 end proc:
19
20 display(quickplot([x->exp(-10*x)],0..4,"x","f(x)"),pointplot(DoBackwardsEuler((x,f)->-10*f,0,4,1,9)) );

27.2.3.5 Nonlinear Example

As a short example for a case where we would require to solve individual equations for each iteration, we will consider the following ODE

$\frac{df}{dx}=-x{e}^{-f}$

When applying Eq. 27.8, we would run into equations with f (x^(i+ 1)) on the left-hand side and within the exponent of the right-hand side. Here, we would need to apply numerical methods to solve these equations, e.g., Newton’s method (see section 26.3). Luckily, Maple will take care of finding solutions to these equations as we call solve on them. Therefore, our function DoBackwardsEuler will handle this ODE without any problem. The analytical solution to Eq. 27.16 is the following function

$f\left(x\right)=ln\left(\frac{x^2}{2}+e\right)$

which is best derived using an Algebra tool. If we apply our function DoBackwardsEuler to Eq. 27.16, we obtain the output shown in Fig. 27.6, which is obviously a very good approximation to the analytical result.

27.2.3.6 Conclusion

In this section we have studied the backward Euler method that is a slight modification of the classical forward Euler method. The backward Euler method is a numerically very stable method and can be used to find solutions, even in cases where the forward Euler method fails. The clear disadvantage of the method is the fact that it requires solving an algebraic equation for each iteration, which is computationally more expensive. However, in many applications, especially if less computationally expensive algorithms fail to obtain stable results, the backward Euler method can be used as a very convenient tool.

27.2.4 Other Important First-Order Methods: Forward and Backward Finite Distance

27.2.4.1 Introduction

In this section we will briefly recapitulate two first-order numerical schemes that we already introduced while discussing the Taylor series in section 4.2.4.4. As they are computationally inexpensive and relatively easy to implement, they form the basis of some of the more rudimentary solvers.

27.2.4.2 Forward Finite Distance

The first scheme we will discuss is the forward finite distance scheme given by Eq. 4.7 as

$\frac{df}{dx}\left(x\right)=\frac{f\left(x+\mathrm{\Delta }x\right)-f\left(x\right)}{\mathrm{\Delta }x}$

We originally developed this method with the aim of having access to the first derivative. Obviously, we can also use it to calculate subsequent data points by rewriting as

$f\left(x+\mathrm{\Delta }x\right)=f\left(x\right)+\frac{df}{dx}\left(x\right)\mathrm{\Delta }x$

or in standard notation

$F^{\left(i+1\right)}=F^{\left(i\right)}+h{X}^{\left(i\right)}$

By comparison with Eq. 27.2, you will find that this is effectively the forward Euler scheme. This scheme is a first-order method with an error function of order two.

27.2.4.3 Backward Finite Distance

The second scheme is the backward finite distance scheme given by Eq. 4.8 as

$\frac{df}{dx}\left(x\right)=\frac{f\left(x\right)-f\left(x-\mathrm{\Delta }x\right)}{\mathrm{\Delta }x}$

which we rewrite to

$f\left(x\right)=f\left(x-\mathrm{\Delta }x\right)+\frac{df}{dx}\left(x\right)\mathrm{\Delta }x$

or in standard notation as

$F^{\left(i\right)}=F^{\left(i-1\right)}+h{X}^{\left(i\right)}$

which again yields a first-order method for approximating F⁽ⁱ⁾.

27.2.5 Important Second-Order Methods

27.2.5.1 Introduction

Besides first-order methods, there are a number of important second-order methods that are commonly used. In this section, we will quickly introduce the most common schemes. Although they are computationally somewhat more expensive to implement, the added accuracy and numerical stability of these methods is generally an advantage to consider. In many cases, it makes sense to default to a second-order scheme, especially in cases where the numerical stability may become critical. All of these methods are one-step methods. We will look at some of the more important multistep methods in just a moment.

27.2.5.2 Crank-Nicolson

We will begin with one method which is especially common in computational fluid dynamics (CFD) and is widely used, e.g., in solving the diffusion equation (Eq. 8.4). We have studied the two one-step first-order Euler methods in section 27.2.2 and section 27.2.3. The forward and the backward Euler methods are two special cases of a wider variety of first-order methods that use weighted averages of these methods. In general these schemes are given as

$\begin{array}{ll}{F}^{\left(i+1\right)}& =\left(1-\omega \right)\left(F^{\left(i\right)}+h{X}^{\left(i\right)}\right)+\omega \left(F^{\left(i\right)}+h{X}^{\left(i+1\right)}\right)\\ & =F^{\left(i\right)}+h\left(\left(1-\omega \right)X^{\left(i\right)}+\omega X^{\left(i+1\right)}\right)\end{array}$

where we assume the next value F^(i+ 1) to be a weighted average of the value obtained from the forward and the backward Euler methods. Please note that this procedure is very similar to SOR where we also use a weighted average of the last and the next value. The weighting factor is the relaxation parameter ω. For Eq. 27.20 this parameter is 0 ≤ ω ≤ 1 where we obtain the forward Euler scheme for ω = 0 and the backward Euler scheme for ω = 1.

In general, every value of ω can be used, but by far the most common scheme uses ω = 0.5. This scheme is referred to as the Crank-Nicolson method named after John Crank¹ and Phyllis Nicolson² who introduced the method in 1947 [2]. It is given by

$F^{\left(i+1\right)}=F^{\left(i\right)}+\frac{h}{2}\left(X^{\left(i\right)}+X^{\left(i+1\right)}\right)$

This linear combination of two first-order methods yields a second-order scheme, i.e., it has an error function of order three. As a consequence, the Crank-Nicolson scheme is significantly more stable than the forward Euler method. Being of second-order, it also yields numerically better results than the backward Euler method which, although being stable, suffers from the imprecision inherent to first-order methods.

Maple Worksheet. We will illustrate the Crank-Nicolson method returning to the example we used when studying the forward (see listing 27.1) and the backward Euler method (see listing 27.2). The implementation of the method is shown in listing 27.3. If you look carefully, you will see that the implementation is identical with the implementation of DoBackwardsEuler in listing 27.2 with the exception of a single line, which is line 10. Here you see that we implement Eq. 27.21 by taking the mean value of the values obtained from the forward and the backward Euler methods. As a consequence, we obtain again an implicit scheme that needs to be solved for the value F^(i+ 1).

1 restart;read("Core.txt"):
2
3 DoCrankNicolson:=proc(df,_from,_to,initialvalue,points) local ret,delta,i,Fi,tosolve:
4  ret:=Matrix(points,2):
5  delta:=(_to-_from)/(points-1):
6  ret[1,1]:=_from:ret[1,2]:=initialvalue:
7  for i from 2 to points do:
8    ret[i,1]:=evalf(ret[i-1,1]+delta):
9    try
10     tosolve:=Fi=ret[i-1,2]+(df(ret[i-1,1],ret[i-1,2])+df(ret[i,1],Fi))*delta/2:
11     ret[i,2]:=rhs(solve({tosolve},{Fi})[1]):
12    catch:
13     print(sprintf("Cannot solve %a",tosolve));
14     return;
15    end try:
16   end do:
17  return ret;
18 end proc:
19
20 display(quickplot([x->exp(-2*x)],0..4,"x","f(x)"),pointplot(DoCrankNicolson((x,f)->-2*f,0,4,1,9 )) );

If we use the Crank-Nicolson method for solving our example ODE Eq. 27.10, we obtain the results shown in Fig. 27.7. In all cases, we obtain stable solutions. In direct comparison to the solutions obtained from the backward Euler method (see Fig. 27.5) the Crank-Nicolson method achieves better precision.

The Crank-Nicolson method is in many cases a good balance if good precision, together with numerical stability, are required. Usually, whenever small stepwidths h are acceptable, forward Euler is the method preferred, as it is sufficiently precise. However, if numerical stability is a problem because larger values of h must be used, implicit schemes must be used. In most of these cases, the Crank-Nicolson method is a good compromise between precision and stability and is often preferred over the pure backward Euler method.

27.2.5.3 Central Finite Distance Scheme

We have already derived the central finite distance scheme while studying some of the most important properties of Taylor series in section 4.2. Originally we derived this equation in order to calculate the derivative $\frac{df}{dx}$ at x_i from the values f(x^(i–1)) and f(x^(i+ 1)). For this we derived Eq. 4.6 which reads

$\frac{df}{dx}\left(x\right)=\frac{f\left(x+\mathrm{\Delta }x\right)-f\left(x-\mathrm{\Delta }x\right)}{2\mathrm{\Delta }x}$

which we can rewrite as a function of f (x + ∆x) as

$f\left(x+\mathrm{\Delta }x\right)=f\left(x-\mathrm{\Delta }x\right)+2\mathrm{\Delta }x\frac{df}{dx}\left(x\right)$

or in the standard notation

$F^{\left(i+1\right)}=F^{\left(i-1\right)}+2h{X}^{\left(i\right)}$

Referring to Eq. 4.6 we find that this is a second-order method, i.e., the error function is of order three, which grants the scheme decent numerical stability. However, it comes with the disadvantage of requiring access to the two most recent data points, a problem known from other schemes, e.g., Adams-Bashforth methods. However, in many applications, a central finite distance scheme can be started, e.g., with an Adams-Moulton scheme in order to obtain the first values. It has the advantage of being computationally inexpensive, which makes it a very convenient method, especially for application where a large number of points must be calculated.

27.2.6 Multistep Methods: Adams-Bashforth Methods

27.2.6.1 Introduction

As we have seen, first-order methods have the advantage of being relatively easy to implement, but the disadvantage of requiring small increments ∆x in order to obtain good solution. We have seen that the backward Euler method can be a convenient solution in these cases, as it is (despite the fact that it is first-order) very stable and often yields results where other algorithms fail. However, it is still a first-order method, i.e., it requires relatively small stepwidths in order to obtain results which approximate the solution correctly.

Obviously, by including higher-order terms in the Taylor series expansions we expect to reduce these errors and come up with methods that allow stable solutions with larger stepwidths. These methods are generally referred to as higher-order or multistep methods. The term higher-order refers to the fact that they take more terms of the Taylor series into account. Multistep refers to the fact that, as we will see, these methods require more than one recent points, i.e., multiple steps in order to produce a result.

In this section we will look at a family of such multistep methods, which are referred to as Adams-Bashforth methods. Interestingly, we have already encountered these methods (see section 24.3), i.e., the application for which these methods were originally developed by Adams and Bashforth while studying the shape of drops [3]. Adams and Bashforth used these schemes to derive the drop contour from a rather tricky differential equation (see Eq. 24.10) that was unsolvable analytically in 1833 and still is today.

27.2.6.2 Derivation

In the following we are again considering a function f (x) that we want to solve numerically. For this we expand Eq. 27.1a and Eq. 27.1b to n_max = 2, respectively. We know from Eq. 4.3 that the Taylor series allows us to write the change in f as x is increased by the infinitesimally small value of ∆x as

$f\left(x+\mathrm{\Delta }x\right)=\frac{1}{0!}\frac{d^{0}f}{d{x}^0}{\left(\mathrm{\Delta }x\right)}^0+\frac{1}{1!}\frac{d^{1}f}{d{x}^1}{\left(\mathrm{\Delta }x\right)}^1+\frac{1}{2!}\frac{d^{2}f}{d{x}^2}{\left(\mathrm{\Delta }x\right)}^2+O_3$

$f\left(x-\mathrm{\Delta }x\right)=\frac{1}{0!}\frac{d^{0}f}{d{x}^0}{\left(\mathrm{\Delta }x\right)}^0-\frac{1}{1!}\frac{d^{1}f}{d{x}^1}{\left(\mathrm{\Delta }x\right)}^1+\frac{1}{2!}\frac{d^{2}f}{d{x}^2}{\left(\mathrm{\Delta }x\right)}^2+O_3$

Eq. 27.23 are the equations from which all Adams-Bashforth methods are derived. The trick here is to use different functions f (x) for the expansion in the negative direction. The obvious way would be to use f (x) directly. However, we could also use the derivative $\frac{df}{dx}\left(x\right)$ when performing the first step. That would effectively mean that we make the Taylor-Series expansion of the derivative.

A Simple Second-Order Formula. We start with the most basic method by simply using f (x) for expanding the Taylor series in the negative direction. Expanding Eq. 27.23b to n_max = 2 we find

$\begin{array}{ll}f\left(x-\mathrm{\Delta }x\right)& =f\left(x\right)-\frac{df}{dx}\left(x\right)\left(\mathrm{\Delta }x\right)+\frac{1}{2}\frac{d^{2}f}{d{x}^2}\left(x\right){\left(\mathrm{\Delta }x\right)}^2+O_3\\ \frac{d^{2}f}{d{x}^2}\left(x\right){\left(\mathrm{\Delta }x\right)}^2& =2\left(f\left(x-\mathrm{\Delta }x\right)-f\left(x\right)+\frac{df}{dx}\left(x\right)\left(\mathrm{\Delta }x\right)\right)-O_3\end{array}$

which we reinsert into Eq. 27.23a expanded to n_max = 2 to find

$\begin{array}{ll}f\left(x+\mathrm{\Delta }x\right)& =f\left(x\right)+\frac{df}{dx}\left(x\right)\left(\mathrm{\Delta }x\right)+f\left(x-\mathrm{\Delta }x\right)-f\left(x\right)+\frac{df}{dx}\left(x\right)\left(\mathrm{\Delta }x\right)+O_3\\ & =f\left(x-\mathrm{\Delta }x\right)+2\frac{df}{dx}\left(x\right)\left(\mathrm{\Delta }x\right)+O_3\end{array}$

or in standard notation

$F^{\left(i+1\right)}=F^{\left(i-1\right)}+2h{X}^{\left(i\right)}$

where we now have an equation that allows us to calculate the value of f (x + ∆x) if we know the value of f (x − ∆x) and of the derivative $\frac{df}{dx}\left(x\right)$. This equation is of second order and the error function of third order. This is already one step better than both of the Euler methods. The main disadvantage is that we require information about the two most recent points in order to predict the next value. As we will see in a moment, this is inherent to multistep methods, and it is the main disadvantage of these methods. If you look closely you will see that this is the equation for the central finite distance scheme (see Eq. 27.22).

Second-Order Adams-Bashforth Formulae. Having derived a first second-order formula, we will now turn to the class of formulae referred to as Adams-Bashforth formulae. We will begin by deriving the formulae of second order. For this, we use the first derivative $\frac{df}{dx}\left(x\right)$ instead of f (x) for the Taylor series expansion in the negative step direction. Therefore we expand Eq. 27.23b to n_max = 1 in which case we find

$\begin{array}{ll}\frac{df}{dx}\left(x-\mathrm{\Delta }x\right)& =\frac{df}{dx}\left(x\right)-\frac{d^{2}f}{d{x}^2}\left(x\right)\left(\mathrm{\Delta }x\right)+O_2\\ \frac{d^{2}f}{d{x}^2}\left(x\right)\left(\mathrm{\Delta }x\right)& =\frac{df}{dx}\left(x\right)-\frac{df}{dx}\left(x-\mathrm{\Delta }x\right)+O_2\end{array}$

Now we reinsert this into Eq. 27.23a. We will use the Taylor series of f (x), not of the derivative as we have done for Eq. 27.23b. Please note that this perfectly allowed: We expanded the Taylor series of the derivative. Inserting this equation into a Taylor series of the function itself is obviously perfectly allowed. Again, we expand Eq. 27.23a to n_max = 2 finding

$f\left(x+\mathrm{\Delta }x\right)=f\left(x\right)+\frac{df}{dx}\left(x\right)\left(\mathrm{\Delta }x\right)+\frac{1}{2}\frac{d^{2}f}{d{x}^2}\left(x\right){\left(\mathrm{\Delta }x\right)}^2+O_3$

into which we insert Eq. 27.24 to find

$\begin{array}{ll}f\left(x+\mathrm{\Delta }x\right)& =f\left(x\right)+\frac{df}{dx}\left(x\right)\left(\mathrm{\Delta }x\right)+\frac{\mathrm{\Delta }x}{2}\left(\frac{df}{dx}\left(x\right)-\frac{df}{dx}\left(x-\mathrm{\Delta }x\right)+O_2\right)+O_3\\ & =f\left(x\right)+\frac{\mathrm{\Delta }x}{2}\left(3\frac{df}{dx}\left(x\right)-\frac{df}{dx}\left(x-\mathrm{\Delta }x\right)\right)+O_3\end{array}$

or in standard notation

$F^{\left(i+1\right)}=F^{\left(i\right)}+\frac{h}{2}\left(3X^{\left(i\right)}-X^{\left(i-1\right)}\right)$

Here again, we have derived a formula that is of second-order, which means that the error function is of order 3. As stated, this is an advantage over the Euler method where we found an error of order 2 (see Eq. 27.2). However, it may not be directly obvious why Eq. 27.24 should be helpful at all, as this equation contains the derivative of f (x) that we presumably do not know. In fact this is true, but this is not really a problem. If you look closely, you will see that we require only the value of the derivative at x and x − ∆x, i.e., at positions which we have already calculated. The iterative scheme is set up such that we calculate the next value of f at x + ∆x. This means that we already know the value of f (x) and of f (x − ∆x). As the ODE we want to solve will give us an equation which puts f (x) in correlation with $\frac{df}{dx}\left(x\right)$, we can calculate these values of the derivative if we know the values of the function. We therefore replaced the requirements of having access to the full derivative (as required for the Taylor series expansion) to the necessity of having access to the value of the derivative at the last discreet points. As you can see from Eq. 27.26, we must know these values for x and for x − ∆x. Obviously, once the iterative scheme is running, these values are available. However, what happens at the very beginning? Usually we will be provided with an initial value or a boundary condition, therefore we know the value f (x_min) (if x_min was the lower boundary). Using the ODE, we would be able to calculate the value of the derivative, which we shall call X (x_min). But in order to calculate the next value f (x_min + ∆x), we require not only the value of X (x_min), but also of X (x_min − ∆x). As this is below the lower boundary, we are not likely to obtain physically meaningful values at all. This means that in order to start this Adams-Bashforth scheme, we must be provided with more than just the value f (x_min). We also require the value of f (x_min + ∆x). If this value was given, then we can obtain both X (x_min) and X (x_min + ∆x) and thus be able to initiate the scheme.

The requirement of having to have access to more than one recent value of f (x) is a disadvantage of the Adams-Bashforth methods. This is inherent to all multistep methods. They take into account the derivative of the function at several previous points, thus making better approximations of the next point of the function. This makes for significantly better approximations as we will see in just a moment. The Euler method we have been studying only takes into account the most recent point to make a prediction about the value of the function f (x) for the next point. This, in turn, requires only one initial value, whereas the Adams-Bashforth methods require more than one point.

Local and Global Truncation Error. You may have noticed that we dropped the O₂ ∆x term during the derivation of Eq. 27.26 and only kept the O₃ term in the final equation. There are two types of errors that we usually consider in numerical methods. The local truncation error is the error we make when deriving the next data point from the previous data point, i.e., while advancing the calculation by one increment. The global truncation error is the accumulation of all errors made during the calculation. Obviously, it would have been necessary to combine the two error terms in Eq. 27.26 as they combine to form the local truncation error.

However, the main cause of errors in Eq. 27.26 results from the incomplete Taylor series expansion, which is why this is assumed to be the main cause of truncation errors. The error introduced by O₂ ∆x can be made arbitrarily small by choosing small values for ∆x. By studying the local truncation errors, it is often possible to estimate the number of steps, i.e., the required minimal stepwidth h to keep the result exact to a given degree.

Adams-Bashforth Formulae. We will not derive more than the Adams-Bashforth formula of second order here. It is sufficient to note that these formulae form a series of equations of increasing order. Obviously, the higher the order of the formulae, the smaller the error function. Tab. 27.3 lists the first members of the Adams-Bashforth family. As you can see, these formulae include one additional previous value for each additional degree. The coefficients can be derived numerically [4]. However, in most cases, it is easier to simply look them up.

Tab. 27.3

Family of the Adams-Moulton formulae. The standard notation used is: function f, independent variable x, discreet point x_i, discreet point x^(i± 1) = x_i ± h, function value at discreet point F⁽ⁱ⁾ = f(x_i), derivative of the function at discreet point$X^{\left(i\right)}=\frac{df}{dx}\left(x_i\right)$, stepwidth h

Equation	Formula order	Order of error function
$F^{\left(i+1\right)}=F^{\left(i\right)}+h{X}^{\left(i+1\right)}$	1	2
$F^{\left(i+1\right)}=F^{\left(i\right)}+\frac{h}{2}\left(X^{\left(i+1\right)}-X^{\left(i\right)}\right)$	2	3
$F^{\left(i+1\right)}=F^{\left(i\right)}+\frac{h}{12}\left(5X^{\left(i+1\right)}+8X^{\left(i\right)}-X^{\left(i-1\right)}\right)$	3	4
$F^{\left(i+1\right)}=F^{\left(i\right)}+\frac{h}{24}\left(19X^{\left(i+1\right)}+9X^{\left(i\right)}-5X^{\left(i-1\right)}+X^{\left(i-2\right)}\right)$	4	5
$F^{\left(i+1\right)}=F^{\left(i\right)}+\frac{h}{720}\left(251X^{\left(i+1\right)}+646X^{\left(i\right)}-264X^{\left(i-1\right)}+106X^{\left(i-2\right)}-19X^{\left(i-3\right)}\right)$	5	6
$F^{\left(i+1\right)}=F^{\left(i\right)}+\frac{h}{1440}\left(475X^{\left(i+1\right)}+1427X^{\left(i\right)}-798X^{\left(i-1\right)}+482X^{\left(i-2\right)}-173X^{\left(i-3\right)}+27X^{\left(i-4\right)}\right)$	6	7

The Problem With the First Values. As you can see from Tab. 27.3 the number of recent points required for the Adams-Bashforth formulae increases with the order of the formula. This makes the Adams-Bashforth formulae difficult to apply at the beginning of a calculation where these points are missing. The most commonly employed way around this problem is to use a different (potentially slower or less stable) algorithm to find the first points and then employ an Adams-Bashforth formula to continue the calculation.

27.2.6.3 An Example

Let us proceed by testing Eq. 27.26, i.e., the second-order Adams-Bashforth formula with our previous example of Eq. 27.6. We seek to find a solution to the ODE Eq. 27.6. We have the initial value Q(0) = F⁽⁰⁾ = 1, which was sufficient for the Euler method. For the Adams-Bashforth method, we require one additional initial point. We will assume a stepwidth of ∆t = 1. As we will see, the Adams-Bashforth method yields significantly better results than the Euler method for this comparatively high stepwidth. As stated, we require the additional value Q(1) = F⁽¹⁾ = 0.61 for the Adams-Bashforth method which we obtained using the analytical solution to the ODE.

Having access to F⁽⁰⁾ and F⁽¹⁾ we can calculate X⁽⁰⁾ and X⁽¹⁾ using Eq. 27.6 as

$\begin{array}{l}{X}^{\left(0\right)}=-\frac{1}{2}{F}^{\left(0\right)}=-\frac{1}{2}Q\left(0\right)=-\frac{1}{2}\\ X^{\left(1\right)}=-\frac{1}{2}{F}^{\left(1\right)}=-\frac{1}{2}Q\left(\mathrm{\Delta }t\right)=-\frac{1}{2}0.606531=-0.31\end{array}$

We can now use Eq. 27.26 to find F⁽²⁾ as

$Q\left(2\mathrm{\Delta }t\right)=F^{\left(2\right)}=F^{\left(1\right)}+\frac{\mathrm{\Delta }t}{2}\left(3X^{\left(1\right)}-X^{\left(0\right)}\right)=0.61+\frac{\mathrm{\Delta }t}{2}\left(-0.93+0.5\right)=0.40$

Once this value has been derived, we can proceed calculating Q(3∆t), etc.

Maple Worksheet. At this point we implement this scheme quickly as a Maple function in order to make working with it more convenient. The code is shown in listing 27.4. It is very similar to listing 27.1 in which we illustrated the Euler method. Again we first start the Maple server and include Core.txt (see line 1). The implementation of Eq. 27.26 starts at line 4 where we first prepare the return value ret and calculate delta, which is h of Eq. 27.26. The function expects the derivative of the function we want to solve, i.e., the right-hand side of Eq. 27.6 as parameter df. This function is used to calculate X⁽ⁱ⁾. It also requires the boundaries of the range for which the solution is to be obtained as _from and _to. As stated, Eq. 27.26 requires two initial values F⁽⁰⁾ and F⁽¹⁾ which must be given to the function as parameters F1 and F2. The number of points from which the stepwidth is calculated is given as parameter points.

1 restart;read("Core.txt"):
2
3 #this function implements the second Adams-Bashforth method
4 DoAdamsBashforth2:=proc(df,_from,_to,F0,F1,points) local ret,delta,i:
5   ret:=Matrix(points,2):
6   delta:=(_to-_from)/(points-1):
7   ret[1,1]:=_from:ret[1,2]:=F0:
8   ret[2,1]:=_from+delta:ret[2,2]:=F1:
9   for i from 3 to points do:
10     ret[i,1]:=evalf(ret[i-1,1]+delta):
11     ret[i,2]:=evalf(ret[i-1,2]+delta/2*(3*df(ret[i-1,1],ret[i-1,2])-df(ret[i-2,1],ret[i-2,2]))):
12   end do:
13  return ret;
14 end proc:
15
16 #solve the ODE numerically
17 display(quickplot([t->exp(-t/2)],0..4,"x","f(x)"),pointplot(DoAdamsBashforth2((t,Q)->-1/2*Q  ,0,4,1,0.606531,4)));

Before starting the iteration, the initial values are inserted as the first two lines of entry in the matrix ret (see line 7). In the loop of the iteration, we first calculate the current value for t (see line 10) before applying Eq. 27.26 in line 11. Here we calculate X⁽ⁱ⁾ and X^(i− 1) using the derivative df given as a function. This is analogous to the calculation of the derivative in listing 27.1.

Line 17 calls the function obtaining a plot for a stepwidth of h = ∆t = 1 and plotting it alongside the analytical solution to Eq. 27.6. The values obtained during the iteration are shown in Tab. 27.4 alongside the values obtained from the Euler method and the correct values from the analytical equation. As you can see, even for this relatively large stepwidth, the Adams-Bashforth method does a very good job approximating the solution. Fig. 27.8 shows the two methods side-by-side together with the analytical solution.

Tab 27.5

Family of the Adams-Moulton formulae. The standard notation used is: function f, independent variable x, discreet point x_i, discreet point x^(i± 1) = x_i ± h, function value at discreet point F⁽ⁱ⁾ = f(x_i), derivative of the function at discreet point X⁽ⁱ⁾ = $\frac{df}{dx}\left(x_i\right)$, stepwidth h

Equation	Formula order	Order of error function
$F^{\left(i+1\right)}=F^{\left(i\right)}+h{X}^{\left(i+1\right)}$	1	2
$F^{\left(i+1\right)}=F^{\left(i\right)}+\frac{h}{2}\left(X^{\left(i+1\right)}-X^{\left(i\right)}\right)$	2	3
$F^{\left(i+1\right)}=F^{\left(i\right)}+\frac{h}{12}\left(5X^{\left(i+1\right)}+8X^{\left(i\right)}-X^{\left(i-1\right)}\right)$	3	4
$F^{\left(i+1\right)}=F^{\left(i\right)}+\frac{h}{24}\left(19X^{\left(i+1\right)}+9X^{\left(i\right)}-5X^{\left(i-1\right)}+X^{\left(i-2\right)}\right)$	4	5
$F^{\left(i+1\right)}=F^{\left(i\right)}+\frac{h}{720}\left(251X^{\left(i+1\right)}+646X^{\left(i\right)}-264X^{\left(i-1\right)}+106X^{\left(i-2\right)}-19X^{\left(i-3\right)}\right)$	5	6
$F^{\left(i+1\right)}=F^{\left(i\right)}+\frac{h}{1440}\left(475X^{\left(i+1\right)}+1427X^{\left(i\right)}-798X^{\left(i-1\right)}+482X^{\left(i-2\right)}-173X^{\left(i-3\right)}+27X^{\left(i-4\right)}\right)$	6	7

27.2.6.4 Summary

The obvious advantage of the second-order Adams-Bashforth formula is the fact that it is of second-order compared to the Euler method, which is only of first order. This means that the error function of the Adams-Bashforth formula is of third order (see Eq. 27.26) compared to the Euler function, which has an error term of second order (see Eq. 27.2). Therefore the Adams-Bashforth method will produce results that are closer to the actual solution, even for large stepwidths. The second advantage is the fact that because the error functions are of higher order, the formula is less likely to become numerically unstable as we have seen with the Euler method. The obvious downside of the Adams-Bashforth method is the fact that we need more initial values.

27.2.7 Predictor-Corrector Methods: Adams-Moulton Methods

27.2.7.1 Introduction

As we have seen, multistep methods have the advantage of good conversion and high stability with the downside of requiring multiple recent points, which is a problem especially at the beginning of a calculation where (potentially) only one data point may be available from the initial or boundary conditions. In this case, we would require methods which would to “find those values themselves”. This is the regime of predictor-corrector methods. These methods work in two steps

1. They first predict a value with which a preliminary calculation is carried out, after which

2. the value obtained is corrected in case the prediction was not entirely correct.

27.2.7.2 Derivation

We have already seen that all schemes for numerical solvers are derived from the Taylor series (Eq. 27.1), which is either evaluated for f (x) or for the derivative $\frac{df}{dx}$, as we have done during the derivation of the second-order Adams-Bashforth formula. However, we are left perfectly free to expand these series around any points we want, not just x + ∆x or x − ∆x. As an example, take the expansion of Eq. 27.1 for f (x) around x + ∆x expanded to n_max = 2 and the expansion of $\frac{df}{dx}$ around x + ∆x expanded to n_max = 1 which yields

$f\left(x+\mathrm{\Delta }x\right)=f\left(x\right)+\frac{df}{dx}\left(x\right)\left(\mathrm{\Delta }x\right)+\frac{1}{2}\frac{d^{2}f}{d{x}^2}\left(x\right){\left(\mathrm{\Delta }x\right)}^2+O_3$

$\frac{df}{dx}\left(x+\mathrm{\Delta }x\right)=\frac{df}{dx}\left(x\right)+\frac{d^{2}f}{d{x}^2}\left(x\right)\left(\mathrm{\Delta }x\right)+O_2$

From Eq. 27.28 we find

$\frac{d^{2}f}{d{x}^2}\left(x\right){\left(\mathrm{\Delta }x\right)}^2=\mathrm{\Delta }x\left(\frac{df}{dx}\left(x+\mathrm{\Delta }x\right)-\frac{df}{dx}\left(x\right)-O_2\right)$

If we now subtract Eq. 27.28 from Eq. 27.27 we find

$\begin{array}{ll}f\left(x+\mathrm{\Delta }x\right)& =f\left(x\right)+\frac{df}{dx}\left(x\right)\left(\mathrm{\Delta }x\right)+\frac{\mathrm{\Delta }x}{2}\left(\frac{df}{dx}\left(x+\mathrm{\Delta }x\right)-\frac{df}{dx}\left(x\right)\right)+O_3\\ & =f\left(x\right)+\frac{\mathrm{\Delta }x}{2}\left(\frac{df}{dx}\left(x+\mathrm{\Delta }x\right)+\frac{df}{dx}\left(x\right)\right)+O_3\end{array}$

where we have again dropped O₂. If we convert Eq. 27.29 to the general notation, we find

$\begin{array}{ll}{F}^{\left(i+1\right)}& =F^{\left(i\right)}+\frac{h}{2}\left(X^{\left(i+1\right)}+X^{\left(i\right)}\right)+O_3\\ & \approx F^{\left(i\right)}+\frac{h}{2}\left(X^{\left(i+1\right)}+X^{\left(i\right)}\right)\end{array}$

We note that the equation we obtained is of second-order, which is a good start. However, we also note that we require X^(i+ 1) and F^(i+ 1) for this equation. In all schemes we derived so far, we could use recent points for the derivative, i.e., X⁽ⁱ⁾ or F⁽ⁱ⁾ in order to calculate F^(i+ 1). This seems to be a problem because we do not know F^(i+ 1) and therefore, we also do not know X^(i+ 1). This seems like a circular reference: We need an unknown first value to calculate an unknown second value, which in turns is required to calculate the first value, etc. This is another example of an implicit scheme as we have already encountered for the backward Euler method (see section 27.2.3).

However, if we were able, by any suitable trick, to circumvent this circular reference, we would be left with one obvious advantage: We do not need any previous values. Therefore this scheme could be an ideal starting algorithm at the beginning of a numerical solution process, as we really only require one value pair F⁽⁰⁾ and X⁽⁰⁾ (both of which we can calculate from the initial or boundary values) to start the algorithm.

Untangling the Circular Reference. We begin by first finding a prediction ${\widehat{F}}^{\left(i+1\right)}$ for F^(i+ 1). In this first step, we want to use an equation which does not require values we do not know yet. For this we note from Eq. 27.27 that a somewhat rough first estimation for ${\widehat{F}}^{\left(i+1\right)}$ would be given by

${\stackrel{⌢}{F}}^{\left(i+1\right)}=f\left(x+\mathrm{\Delta }x\right)\approx f\left(x\right)+\frac{df}{dx}\left(x\right)\left(\mathrm{\Delta }x\right)=F^{\left(i\right)}+h{X}^{\left(i\right)}$

where we simply expanded the Taylor series to only n_max = 1. Obviously, we know already that this is not a good approximation, as it is first-order. However, we may start our algorithm using this first guess. Once we calculated the prediction value ${\widehat{F}}_{\left(i+1\right)}$, we can use this value in Eq. 27.30 and calculate the corrected value. The predictor-corrector scheme therefore consists of the following steps

1. Calculate ${\widehat{F}}^{\left(i+1\right)}$ using Eq. 27.31

2. Use this value as value F^(i+ 1) for the right-hand side of Eq. 27.30 to calculate X^(i+ 1) hence, obtaining the corrected value for F^(i+ 1) as the left-hand side of Eq. 27.29

Adams-Moulton Formulae. The method we have just derived is one of a family of methods that are referred to as Adams-Moulton formulae. These methods are predictor-corrector methods, which use information about recent points. They are therefore also multistep methods. The first members of the family of Adams-Moulton formulae are shown in Tab. 27.5. The scheme we just derived is the second-order formula.

Tab. 27.6

Family of the Runge-Kutta formulae. The standard notation used is: function f, independent variable x, discreet point x_i, discreet point x^(i±1) = x_i ± h, function value at discreet point F⁽ⁱ⁾ = f(x_i), derivative of the function at discreet point$X^{\left(i\right)}=\frac{df}{dx}\left(x_i\right)$stepwidth h

Order	Equation	$c_1=\frac{df}{dx}\left(\text{,}\right)$	$c_2=\frac{df}{dx}\left(\text{,}\right)$	$c_3=\frac{df}{dx}\left(\text{,}\right)$	$c_4=\frac{df}{dx}\left(\text{,}\right)$	Order of error function
1	F(i+ 1)= F(i) + h .c1	x(i),F(i)				2
2	F(i+ 1)= F(i) + h .c2	x(i),F(i)	$x^{\left(i\right)}+\frac{h}{2},F\left(i\right)+\frac{c_1}{2}$			3
3	$F^{\left(i+1\right)}=F^{\left(i\right)}+\frac{h}{6}\left(c_1+c_2+c_3\right)$	x(i),F(i)	$x^{\left(i\right)}+\frac{h}{2},F^{\left(i\right)}+\frac{c_1}{2}$	$x^{\left(i\right)}+h,F^{\left(i\right)}-c_1+2c_2$		4
4	$F^{\left(i+1\right)}=F^{\left(i\right)}+\frac{h}{6}\left(c_1+2\cdot c_2+2\cdot c_3+c_4\right)$	x(i),F(i)	$x^{\left(i\right)}+\frac{h}{2},F^{\left(i\right)}+\frac{c_1}{2}$	$x^{\left(i\right)}+\frac{h}{2},F^{\left(i\right)}+\frac{c_2}{2}$	$x^{\left(i\right)}+h,F^{\left(i\right)}+c_3$	5

Adams-Bashforth Predictors. If you look at Eq. 27.31, you can see that this is actually the first-order Adams-Bashforth formula (see Tab. 27.3). This, in fact, is a very common procedure: Combining a Adams- Bashforth formula of order p – 1 as predictor method with an Adams-Moulton formula of order p. Obviously, there are other predictor methods that could also be used, e.g., any first-order method, such as Euler or central difference schemes.

27.2.7.3 Example

In order to illustrate the usage of this scheme we will return to the example that we used during the illustration of the Euler and the Adams-Bashforth methods. Again, we seek to find a solution to the ODE Eq. 27.6, again using a relatively large stepwidth Δt = 1. We have found that this yielded quite bad values when using the Euler method, but decent values when using the second Adams-Bashforth scheme. However, the latter required two points to start the algorithm. Considering that we only evaluated a total of four points (see Tab. 27.4) we already provided the algorithm with half the solution. We expect the implicit scheme we just derived to perform better, as we only require the initial condition which is given.

First Prediction. The initial condition gives us F⁽⁰⁾ = 1 from which we can calculate X⁽⁰⁾ = − 0.5 using Eq. 27.6. For our first prediction we use Eq. 27.31 to find

${\widehat{F}}^{\left(1\right)}=F^{\left(0\right)}+h{X}^{\left(0\right)}=1-0.5=0.5$

We already know from the analytical solution that this is not the correct value. The correct value should be 0.606531 (see Tab. 27.4). However, this was only the prediction.

First Correction. We will now use Eq. 27.30 to correct this first prediction. For this we calculate X⁽¹⁾ = − 0.25 using Eq. 27.6. We then use this value for the right-hand side of Eq. 27.30 to find F⁽¹⁾ as

$F^{\left(1\right)}=F^{\left(0\right)}+\frac{h}{2}\left(X^{\left(1\right)}+X^{\left(0\right)}\right)=1+0.5\left(-0.25-0.5\right)=0.625$

which is significantly closer to the actual solution.

Iterative Scheme. Obviously, we would now use the value of F⁽¹⁾ to reiterate this scheme. At this point, we will once again turn to Maple to create a worksheet that will perform these calculations for us.

27.2.7.4 Maple Worksheet

The listing which implements the second-order Adams-Moulton formula is shown in listing 27.5. As usual, we first restart the Maple server and include Core.txt line 1. The function AdamsMoulton2 implements the scheme we just performed manually. It accepts the differential equation as a function of the independent variable t and the function value Q as the parameter df. It also requires the range given by _from and _to, as well as the initial value F0. Please note that compared to the Adams-Bashforth method we do not require a second initial value here. The number of points required is given as parameter points.

 1 restart;read("Core.txt"):
 2
 3 #this function implements the second Adams-Moulton method
 4 AdamsMoulton2:=proc(df,_from,_to,FO,points) local ret,delta,i,Fpred:
 5  ret:=Matrix(points,2):
 6  delta:=(_to-_from)/(points-1):
 7  ret [1,1] : = _from:ret [1,2]:=F0:
 8  for i from 2 to points do:
 9    ret[i,1]:=evalf(ret[i-1,1]+delta):
10    Fpred:=ret[i-1,2]+delta*df(ret[i-1,1] ,ret[i-1,2]) :
11    ret[i,2]:=evalf(ret[i-1,2]+delta/2*(df(ret[i,1],Fpred)+df(ret[i-1,1],ret[i-1,2]))):
12   end do:
13  return ret;
14 end proc:
15
16 #solve the ODE numerically
17. display(quickplot([t->exp(-t/2)],0..4,"x","f(x)"), pointplot(AdamsMoulton2((t,Q)->-1/2*Q,0,4,1,4))));

The function then proceeds by calculating the stepwidth h as variable delta and inserts the initial value (see line 7). The prediction according to Eq. 27.31 is implemented in line 10. Here, ${\widehat{F}}^{\left(i+1\right)}$ is calculated based on the values t⁽ⁱ⁾ and F⁽ⁱ⁾ both of which we know. The prediction is then corrected in line 11 according to Eq. 27.30.

The listing plots the numerically obtained function in comparison to the analytical solution in line 17. The plot is shown in Fig. 27.9a. As you can see, the solution is very close to the correct solution. Fig. 27.9b shows the results of the Euler, second-order Adams-Bashforth and second-order Adams-Moulton methods in comparison. As you can see, the output of the Adams-Bashforth and the Adams-Moulton methods are very close to the

solution with the Adams-Moulton method yielding slightly better results. In addition this method comes with the benefit of not having to provide two initial values, as we did for the Adams-Bashforth method. Tab. 27.4 shows the output values of the three methods for better comparison.

27.2.7.5 Summary

In this section we have studied predictor-corrector schemes, which are implicit methods. Here, the term implicit refers to the fact that these schemes do not allow us to calculate the next data point directly, as the formulae give the expression of the next point implicitly. This means that we need a first prediction for the value that we correct in a second step. The most important class of corrector formulae are the Adams-Moulton formulae which we have derived and studied. Usually, the predictor methods of choice are Adams-Bashforth formulae of order p – 1, which are then combined with a Adams-Moulton corrector formula of order p. We have studied examples of such a combination.

As we have shown previously, implicit schemes are significantly more stable than explicit schemes, e.g., the Adams-Bashforth methods alone. However, just as with the forward and the backward Euler method, the implicit methods are significantly more expensive computationally. Here, we require two evaluations per step. In contrast, the backward Euler method required finding solutions to an algebraic equation or a system of algebraic equations. The Adams-Moulton schemes are less tricky to use, as we do not have to solve an equation, but merely evaluate each point twice.

27.2.8 Modern One-Step Methods: Runge-Kutta Methods

27.2.8.1 Introduction

One particularly interesting family of schemes for numerical solutions are the Runge-Kutta methods. They originate from the seminal work of Runge, which is considered the starting point of modern one-step methods. The method was significantly improved first by Heun [5] and later by Kutta [6], who derived the most famous of the Runge-Kutta method, i.e., the fourth-order scheme.

Rationale. Being one-step methods, the Runge-Kutta methods only require the most recent point in order to calculate the next. However, compared to the other methods discussed so far, they make use of the fact that the ODE to solve can also be evaluated at points in between the discreet points for which a solution must be obtained. As an example, the forward Euler method (see Eq. 27.8) approximates the next point F^(i+ 1) by taking into account the derivative X⁽ⁱ⁾ at the current point. Obviously, this approximation may become significantly better, if we could take into account the derivative at several points in between x_i and x_i + h. Obviously, this is possible, as the ODE can be evaluated for any intermediary point x_i + with 0 ≤ ≤ h.

The Runge-Kutta formulae are a family of methods which take into account the derivative of f (x) at several intermediary points. The second-order method uses two intermediary points, the third-order method, three, and the fourth-order, four. In general, a Runge-Kutta method of order n will require n separate evaluations, thus the methods become more expensive computationally. However, with increasing order they obviously become more exact and more stable as the degree of the error function increases. For the second-order Runge-Kutta method, you will commonly find the abbreviation “RK2”, for the third-order the abbreviation “RK3”, and “RK4” for the fourth-order equation.

Runge-Kutta Formulae. The derivation of these equations is a bit more complex and beyond the scope of this book. They all follow the same concept of using additional intermediary points when evaluating the gradient of the function f (x). Tab. 27.6 gives the equations of the first Runge-Kutta methods. We already know one of these formulae: If you compare the first-order Runge-Kutta method (see Tab. 27.6) with the forward Euler equation (see Eq. 27.2) you will see that they are identical.

In the following, we will use the fourth-order Runge-Kutta method which is the one you will most likely come across.

27.2.8.2 Example

For this example, we return to the ODE of the capacitor discharge given by Eq. 27.6. We apply the fourth-order Runge-Kutta method to solve this equation. For this, we need to calculate the four coefficients given by Tab. 27.6 using the values F⁽⁰⁾ = 1 at t₀ = 0 and X⁽⁰⁾ = − 0.5, which we can calculate from Eq. 27.6. Again, we will use a stepwidth of h = ∆t = 1. We then proceed to calculate the coefficients as

$\begin{array}{llll}{c}_1\hfill & =\frac{dQ}{dt}\left(t_0,F^{\left(0\right)}\right)\hfill & =\frac{dQ}{dt}\left(0,1\right)\hfill & =-0.5\hfill \\ c_2\hfill & =\frac{dQ}{dt}\left(t_0+\frac{h}{2},F^{\left(0\right)}+\frac{c_1}{2}\right)\hfill & =\frac{dQ}{dt}\left(0.5,1-0.25\right)\hfill & =-0.375\hfill \\ c_3\hfill & =\frac{dQ}{dt}\left(t_0+\frac{h}{2},F^{\left(0\right)}+\frac{c_2}{2}\right)\hfill & =\frac{dQ}{dt}\left(0.5,1-0.189\right)\hfill & =-0.406\hfill \\ c_4\hfill & =\frac{dQ}{dt}\left(t_0+h,F^{\left(0\right)}+c_3\right)\hfill & =\frac{dQ}{dt}\left(0.5,1-0.406\right)\hfill & =-0.297\hfill \end{array}$

Here you can already see the disadvantage of the Runge-Kutta method: We have to calculate four coefficients, evaluating the ODE at four distinct locations. Obviously, this comes with computational overhead. However, the approximation we obtain for the first value are obviously significantly better. We find

$\begin{array}{ll}{F}^{\left(i+1\right)}\hfill & =F^{\left(i\right)}+\frac{h}{6}\left(c_1+2c_2+2c_3+c_4\right)\hfill \\ \hfill & =1+\frac{1}{6}\left(-0.5-20.375-0.297\right)=0.60677\hfill \end{array}$

which is very close to the exact solution 0.606531

27.2.8.3 Maple Worksheet

As for the other schemes, we will now implement this scheme as a Maple worksheet. The listing is shown in listing 27.6. As usual, we first restart the Maple server and include Core.txt. The function that implements the fourth-order Runge-Kutta scheme is called RungeKutta4. It expects the ODE as a function of the independent variable t and of the function value Q as parameter df. The range is again given in the form of the upper and lower boundary _from and _to with a total of points from which the stepwidth h is calculated as variable delta. Being a one-step method, the Runge-Kutta function only requires one initial value F0 as parameter.

 1 restart;read("Core.txt"):
 2
 3 #this function implements the fourth order Runge-Kutta method
 4 RungeKutta4:=proc(df,_from,_to,F0,points) local ret,delta,i,c1,c2,c3,c4:
 5  ret:=Matrix(points,2):
 6  delta:=(_to-_from)/(points-1):
 7  ret [1,1] : = _f rom:ret [1,2] : = F0:
 8  for i from 2 to points do:
 9    ret [i,1]:=evalf(ret[i-1,1]+delta):
10    c1 :=df (ret [i-1,1] ,ret [i-1,2]):
11    c2:=df(ret[i-1,1]+delta/2,ret[i-1,2]+c1/2):
12    c3:=df (ret [i-1, l]+delta/2 ,ret [i-1,2] + c2/2):
13    c4:=df(ret[i-1,1]+delta,ret[i-1,2]+c3):
14    ret[i,2]:=evalf(ret[i-1,2]+delta/6*(c1+2*c2+2*c3+c4)):
15   end do:
16  return ret;
17 end proc:
18
19 #solve the ODE numerically
20 display(quickplot ([t->exp(-t/2)] ,0. .4, "x","f (x)") .pointplot(RungeKutta4((t,Q)->-1/2*Q,0,4,1,5)));

The function then proceeds preparing the return value ret and calculating the stepwidth delta. The initial value is set in line 7. The algorithm then proceeds iteratively in a loop first calculating the parameters c1, c2, c3, and c4 before calculating the next value according to Tab. 27.6.

The function is then called and the result plotted as point plot alongside the analytical solution in line 20. Fig. 27.10a shows this numerical output obtained from RungeKutta4 when applied to Eq. 27.6 alongside the analytical solution. Fig. 27.10b shows a comparison of all the methods we have studied so far applied to the same equation. Tab. 27.4 shows the numerical output in tabular format. As you can see, the results are very close to the analytical solution. Of all the algorithms we have used to far, the Runge-Kutta method gives the closest approximation. However, it comes with the disadvantage of being numerically more expensive. However, in most applications this is not so critical and Runge-Kutta remains one of the best algorithms to try on a numerical problem, being relatively stable and decently quick in converging.

27.2.8.4 Conclusion

The Runge-Kutta methods (foremost the fourth-order method) are among the most important schemes in numerics. The main advantage of the Runge-Kutta methods is the fact that they do not require knowledge of recent points being one-step methods. This is a significant advantage over Adams-Bashforth methods which require more recent points for the higher-order formulae. As we have experienced from working with first-order methods, the predictions using only a single point are usually not very good. However, the Runge-Kutta methods make use of additional intermediary supporting points in between the discreet intervals and thus, exploit the fact that the ODE is not a discreet function and will therefore, yield information about the value of the derivative in supporting points in between the discreet values for which the numerical scheme must derive a solution. This makes these methods significantly more exact.

The fact that they are numerically expensive makes the Runge-Kutta not the method of choice if a numerical problem with a large number of points must be solved. In these scenarios, they are usually combined with a predictor-corrector method. A common combination is the fourth-order Runge-Kutta method which is used only to calculate the first two points of a problem. At this point, the scheme is changed to a third-order Adams-Bashforth scheme predictor method, combined with a fourth-order Adams-Moulton corrector method for completing the calculation. This combined scheme is often used, as it is (in total) a third-order method with a fourth-order error function, which makes for good numerical stability.