#### **7.1.2 Brute Force Optimization**

**Numerical Algorithm** We realize that *xi* corresponds to a maximum point if *yi*−1 *< yi > yi*+1. Similarly, *xi* corresponds to a minimum if *yi*−1 *> yi < yi*+1. We can do this test for all "inner" points *i* = 1*,...,n* − 1 to find all local minima and maxima. In addition, we need to add an end point, *i* = 0 or *i* = *n*, if the corresponding *yi* is a global maximum or minimum.

**Implementation** The algorithm above can be translated to the following Python function (file brute\_force\_optimizer.py):

```
def brute_force_optimizer(f, a, b, n):
    from numpy import linspace
    x = linspace(a, b, n)
    y = f(x)
    # Let maxima and minima hold the indices corresponding
    # to (local) maxima and minima points
    minima = []
    maxima = []
    for i in range(1, n-1):
        if y[i-1] < y[i] > y[i+1]:
            maxima.append(i)
        if y[i-1] > y[i] < y[i+1]:
            minima.append(i)
    # What about the end points?
    y_max_inner = max([y[i] for i in maxima])
    y_min_inner = min([y[i] for i in minima])
    if y[0] > y_max_inner:
        maxima.append(0)
    if y[len(x)-1] > y_max_inner:
        maxima.append(len(x)-1)
    if y[0] < y_min_inner:
        minima.append(0)
    if y[len(x)-1] < y_min_inner:
        minima.append(len(x)-1)
```

```
# Return x and y values
return [(x[i], y[i]) for i in minima], \
       [(x[i], y[i]) for i in maxima]
```
The max and min functions are standard Python functions for finding the maximum and minimum element of a list or an object that one can iterate over with a for loop. An application to *f (x)* = *e*−*x*2 cos*(*4*x)* looks like

```
def demo():
    from numpy import exp, cos
    minima, maxima = brute_force_optimizer(
        lambda x: exp(-x**2)*cos(4*x), 0, 4, 1001)
    print('Minima:\n', minima)
    print('Maxima:\n', maxima)
```
Running the program gives

```
Minima:
 [(0.70000000000000007, -0.5772302750838405), (2.1520000000000001,
Maxima:
 [(1.4159999999999999, 0.10965467991643564), (2.8999999999999999,
0.00012651823896373234), (0.0, 1.0)]
```
#### **7.1.3 Model Problem for Algebraic Equations**

We shall consider the very simple problem of finding the square root of 9. That is, we want to solve *x*2 = 9, but will (for simplicity) seek only the positive solution. Knowing the solution beforehand, allows us to easily investigate how the numerical method (and the implementation of it) performs in the search for a solution. The *f (x)* function that corresponds to the equation *x*2 = 9 is

$$f(x) = x^2 - 9\dots$$

Our interval of interest for solutions will be [0*,* 1000] (the upper limit here is chosen somewhat arbitrarily).

In the following, we will present several efficient and accurate methods for solving nonlinear algebraic equations, both single equation and systems of equations. The methods all have in common that they search for *approximate* solutions. The methods differ, however, in the way they perform the search for solutions. The idea for the search influences the efficiency of the search and the reliability of actually finding a solution. For example, Newton's method is very fast, but not reliable, while the bisection method is the slowest, but absolutely reliable. No method is best at all problems, so we need different methods for different problems.

#### **What Is the Difference Between Linear and Nonlinear Equations?**

You know how to solve linear equations *ax* + *b* = 0: *x* = −*b/a*. All other types of equations *f (x)* = 0, i.e., when *f (x)* is not a linear function of *x*, are called nonlinear. A typical way of recognizing a nonlinear equation is to observe that *x* is "not alone" as in *ax*, but involved in a product with itself, such as in *x*3 + 2*x*2 − 9 = 0. We say that *x*3 and 2*x*2 are nonlinear terms. An equation like sin *x* +*ex* cos *x* = 0 is also nonlinear although *x* is not explicitly multiplied by itself, but the Taylor series of sin *x*, *ex*, and cos *x* all involve polynomials of *x* where *x* is multiplied by itself.

#### **7.2 Newton's Method**

*Newton's method*, also known as *Newton-Raphson's method*, is a very famous and widely used method for solving nonlinear algebraic equations.1 Compared to the other methods presented in this chapter, i.e., secant and bisection, it is generally the fastest one (although computational speed rarely is an issue with a single equation on modern laptops). However, it does not guarantee that an existing solution will be found.

A fundamental idea of numerical methods for nonlinear equations is to construct a series of linear equations (since we know how to solve linear equations) and hope that the solutions of these linear equations bring us closer and closer to the solution of the nonlinear equation. The idea will be clearer when we present Newton's method and the secant method.

#### **7.2.1 Deriving and Implementing Newton's Method**

Figure 7.1 shows the *f (x)* function in our model equation *x*2 − 9 = 0. Numerical methods for algebraic equations require us to guess at a solution first. Here, this guess is called *x*0. The fundamental idea of Newton's method is to approximate the original function *f (x)* by a straight line, i.e., a linear function, since it is straightforward to solve linear equations. There are infinitely many choices of how to approximate *f (x)* by a straight line. Newton's method applies the tangent of *f (x)* at *x*0, see the rightmost tangent in Fig. 7.1. This linear tangent function crosses the *x* axis at a point we call *x*1. This is (hopefully) a better approximation to the solution of *f (x)* = 0 than *x*0. The next fundamental idea is to repeat this process. We find the tangent of *f* at *x*1, compute where it crosses the *x* axis, at a point called *x*2, and repeat the process again. Figure 7.1 shows that the process brings us closer and closer to the left. It remains, however, to see if we hit *x* = 3 or come sufficiently close to this solution.

1 Read more about Newton's method, e.g., on https://en.wikipedia.org/wiki/Newton%27s\_method.

**Fig. 7.1** Illustrates the idea of Newton's method with *f (x)* = *x*2 − 9, repeatedly solving for crossing of tangent lines with the *x* axis

How do we compute the tangent of a function *f (x)* at a point *x*0? The tangent function, here called *f (x)* ˜ , is linear and has two properties:

So, if we write the tangent function as *f (x)* ˜ = *ax* + *b*, we must require *f*˜ *(x*0*)* = *f (x*0*)* and *f (x* ˜ 0*)* = *f (x*0*)*, resulting in

$$\bar{f}(\mathbf{x}) = f(\mathbf{x}\_0) + f'(\mathbf{x}\_0)(\mathbf{x} - \mathbf{x}\_0) \dots$$

The key step in Newton's method is to find where the tangent crosses the *x* axis, which means solving *f (x)* ˜ = 0:

$$
\tilde{f}(\mathbf{x}) = \mathbf{0} \quad \Rightarrow \quad \mathbf{x} = \mathbf{x}\_0 - \frac{f(\mathbf{x}\_0)}{f'(\mathbf{x}\_0)}.
$$

This is our new candidate point, which we call *x*1:

$$x\_1 = x\_0 - \frac{f(x\_0)}{f'(x\_0)} \dots$$

With *x*0 = 1000, we get *x*1 ≈ 500, which is in accordance with the graph in Fig. 7.1. Repeating the process, we get

$$\mathbf{x}\_2 = \mathbf{x}\_1 - \frac{f(\mathbf{x}\_1)}{f'(\mathbf{x}\_1)} \approx \mathbf{250}\dots$$

The general scheme2 of Newton's method may be written as

$$\mathbf{x}\_{n+1} = \mathbf{x}\_n - \frac{f(\mathbf{x}\_n)}{f'(\mathbf{x}\_n)}, \quad n = 0, 1, 2, \dots \tag{7.1}$$

The computation in (7.1) is repeated until *f (xn)* is close enough to zero. More precisely, we test if |*f (xn)*| *<* , with being a small number.

We moved from 1000 to 250 in two iterations, so it is exciting to see how fast we can approach the solution *x* = 3. A computer program can automate the calculations. Our first try at implementing Newton's method is in a function naive\_Newton (found in naive\_Newton.py):

```
def naive_Newton(f, dfdx, x, eps):
    while abs(f(x)) > eps:
        x = x - (f(x))/dfdx(x)
    return x
```
The argument x is the starting value, called *x*0 in our previous mathematical description.

To solve the problem *x*2 = 9 we also need to implement

```
def f(x):
    return x**2 - 9
def dfdx(x):
    return 2*x
```

```
print(naive_Newton(f, dfdx, 1000, 0.001))
```
which in naive\_Newton.py is included by use of an extra function and a test block.

2 The term *scheme* is often used as a synonym for method or computational recipe.

#### **Why Not Use an Array for the** *x* **Approximations?**

Newton's method is normally formulated with an *iteration index n*,

$$\mathbf{x}\_{n+1} = \mathbf{x}\_n - \frac{f(\mathbf{x}\_n)}{f'(\mathbf{x}\_n)}.$$

Seeing such an index, many would implement this as

x[n+1] = x[n] - f(x[n])/dfdx(x[n])

Such an array is fine, but requires storage of all the approximations. In large industrial applications, where Newton's method solves millions of equations at once, one cannot afford to store all the intermediate approximations in memory, so then it is important to understand that the algorithm in Newton's method has no more need for *xn* when *xn*+1 is computed. Therefore, we can work with one variable x and overwrite the previous value:

x = x - f(x)/dfdx(x)

Running naive\_Newton(f, dfdx, 1000, eps=0.001) results in the approximate solution 3.000027639. A smaller value of eps will produce a more accurate solution. Unfortunately, the plain naive\_Newton function does not return how many iterations it used, nor does it print out all the approximations *x*0*, x*1*, x*2*,...*, which would indeed be a nice feature. If we insert such a printout (print(x) in the while loop), a rerun results in

```
500.0045
250.011249919
125.02362415
62.5478052723
31.3458476066
15.816483488
8.1927550496
4.64564330569
3.2914711388
3.01290538807
3.00002763928
```
We clearly see that the iterations approach the solution quickly. This speed of the search for the solution is the primary strength of Newton's method compared to other methods.

#### **7.2.2 Making a More Efficient and Robust Implementation**

The naive\_Newton function works fine for the example we are considering here. However, for more general use, there are some pitfalls that should be fixed in an improved version of the code. An example may illustrate what the problem is.

Let us use naive\_Newton to solve tanh*(x)* = 0, which has solution *x* = 0 (interactively, you may define *f (x)* = *tanh(x)* and *f (x)* = 1−*tanh*2*(x)* as Python functions and import the naive\_Newton function from naive\_Newton.py). With |*x*0| ≤ 1*.*08 everything works fine. For example, *x*0 = 1*.*08 leads to six iterations if = 0*.*001:

Adjusting *x*0 slightly to 1.09 gives division by zero! The approximations computed by Newton's method become

```
1.10490354324
1.30303261823
13.4731428006
```
The division by zero is caused by *x*7 = −1*.*26055913647× 1011, because tanh*(x*7*)* is 1.0 to machine precision, and then *f (x)* = 1 − tanh*(x)*2 becomes zero in the denominator in Newton's method.

The underlying problem, leading to the division by zero in the above example, is that Newton's method *diverges*: the approximations move further and further away from *x* = 0. If it had not been for the division by zero, the condition in the while loop would always be true and the loop would run forever. Divergence of Newton's method occasionally happens, and the remedy is to abort the method when a maximum number of iterations is reached.

Another disadvantage of the naive\_Newton function is that it calls the *f (x)* function twice as many times as necessary. This extra work is of no concern when *f (x)* is fast to evaluate, but in large-scale industrial software, one call to *f (x)* might take hours or days, and then removing unnecessary calls is important. The solution in our function is to store the call f(x) in a variable (f\_value) and reuse the value instead of making a new call f(x).

To summarize, we want to write an improved function for implementing Newton's method where we

A more robust and efficient version of the function, inserted in a complete program (Newtons\_method.py) for solving *x*2 − 9 = 0, is listed below.

```
import sys
def Newton(f, dfdx, x, eps):
    f_value = f(x)
```

```
iteration_counter = 0
    while abs(f_value) > eps and iteration_counter < 100:
        try:
            x = x - f_value/dfdx(x)
        except ZeroDivisionError:
            print('Error! - derivative zero forx=', x)
            sys.exit(1) # Abort with error
        f_value = f(x)
        iteration_counter = iteration_counter + 1
    # Here, either a solution is found, or too many iterations
    if abs(f_value) > eps:
        iteration_counter = -1
    return x, iteration_counter
if __name__ == '__main__':
    def f(x):
        return x**2 - 9
    def dfdx(x):
        return 2*x
    solution, no_iterations = Newton(f, dfdx, x=1000, eps=1.0e-6)
    if no_iterations > 0: # Solution found
        print('Number of function calls: {:d}'.format(1+2*no_iterations))
        print('A solution is: {:f}'.format(solution))
    else:
        print('Solution not found!')
```
Handling of the potential division by zero is done by a try-except construction.3

The division by zero will always be detected and the program will be stopped. The main purpose of our way of treating the division by zero is to give the user a more informative error message and stop the program in a gentler way.

Calling sys.exit with an argument different from zero (here 1) signifies that the program stopped because of an error. It is a good habit to supply the value 1, because tools in the operating system can then be used by other programs to detect that our program failed.

To prevent an infinite loop because of divergent iterations, we have introduced the integer variable iteration\_counter to count the number of iterations in Newton's method. With iteration\_counterwe can easily extend the condition in the while loop such that no more iterations take place when the number of iterations reaches 100. We could easily let this limit be an argument to the function rather than a fixed constant.

The Newton function returns the approximate solution and the number of iterations. The latter equals −1 if the convergence criterion |*f (x)*| *<*  was not reached within the maximum number of iterations. In the calling code, we print out

3 Professional programmers would avoid calling sys.exit inside a function. Instead, they would raise a new exception with an informative error message, and let the calling code have another try-except construction to stop the program.

the solution and the number of function calls. The main cost of a method for solving *f (x)* = 0 equations is usually the evaluation of *f (x)* and *f (x)*, so the total number of calls to these functions is an interesting measure of the computational work. Note that in function Newton there is an initial call to *f (x)* and then one call to *f* and one to *f* in each iteration.

Running Newtons\_method.py, we get the following printout on the screen:

```
Number of function calls: 25
A solution is: 3.000000
```
The Newton scheme will work better if the starting value is close to the solution. A good starting value may often make the difference as to whether the code actually *finds* a solution or not. Because of its speed (and when speed matters), Newton's method is often the method of first choice for solving nonlinear algebraic equations, even if the scheme is not guaranteed to work. In cases where the initial guess may be far from the solution, a good strategy is to run a few iterations with the bisection method (see Sect. 7.4) to narrow down the region where *f* is close to zero and then switch to Newton's method for fast convergence to the solution.

**Using sympy to Find the Derivative** Newton's method requires the analytical expression for the derivative *f (x)*. Derivation of *f (x)* is not always a reliable process by hand if *f (x)* is a complicated function. However, Python has the symbolic package SymPy, which we may use to create the required dfdx function. With our sample problem, we get:

```
import sympy as sym
x = sym.symbols('x')
f_expr = x**2 - 9 # symbolic expression for f(x)
dfdx_expr = sym.diff(f_expr, x) # compute f'(x) symbolically
# turn f_expr and dfdx_expr into plain Python functions
f = sym.lambdify([x], # argument to f
            f_expr) # symbolic expression to be evaluated
dfdx = sym.lambdify([x], dfdx_expr)
print(f(3), dfdx(3)) # will print 0 and 6
```
The nice feature of this code snippet is that dfdx\_expr is the exact analytical expression for the derivative, 2\*x (seen if you print it out). This is a symbolic expression, so we cannot do numerical computing with it. However, with lambdify, such symbolic expression are turned into callable Python functions, as seen here with f and dfdx.

The next method is the secant method, which is usually slower than Newton's method, but it does not require an expression for *f (x)*, and it has only one function call per iteration.

#### **7.3 The Secant Method**

When finding the derivative *f (x)* in Newton's method is problematic, or when function evaluations take too long; we may adjust the method slightly. Instead of using tangent lines to the graph we may use secants. 4 The approach is referred to as the *secant method*, and the idea is illustrated graphically in Fig. 7.2 for our example problem *x*2 − 9 = 0.

The idea of the secant method is to think as in Newton's method, but instead of using *f (xn)*, we approximate this derivative by a finite difference or the *secant*, i.e., the slope of the straight line that goes through the points *(xn, f (xn))* and *(xn*−1*, f (xn*−1*))* on the graph, given by the two most recent approximations *xn* and *xn*−1. This slope reads

$$\frac{f(\mathbf{x}\_n) - f(\mathbf{x}\_{n-1})}{\mathbf{x}\_n - \mathbf{x}\_{n-1}}.\tag{7.2}$$

Inserting this expression for *f (xn)* in Newton's method simply gives us the secant method:

$$
\alpha\_{n+1} = \alpha\_n - \frac{f\left(\chi\_n\right)}{\frac{f\left(\chi\_n\right) - f\left(\chi\_{n-1}\right)}{\chi\_n - \chi\_{n-1}}},
$$

or

$$\mathbf{x}\_{n+1} = \mathbf{x}\_n - f(\mathbf{x}\_n) \frac{\mathbf{x}\_n - \mathbf{x}\_{n-1}}{f(\mathbf{x}\_n) - f(\mathbf{x}\_{n-1})} \,. \tag{7.3}$$

Comparing (7.3) to the graph in Fig. 7.2, we see how *two* chosen starting points (*x*0 = 1000, *x*1 = 700, and corresponding function values) are used to compute *x*2. Once we have *x*2, we similarly use *x*1 and *x*2 to compute *x*3. As with Newton's method, the procedure is repeated until *f (xn)* is below some chosen limit value, or some limit on the number of iterations has been reached. We use an iteration counter here too, based on the same thinking as in the implementation of Newton's method.

We can store the approximations *xn* in an array, but as in Newton's method, we notice that the computation of *xn*+1 only needs knowledge of *xn* and *xn*−1, not "older" approximations. Therefore, we can make use of only three variables: x for *xn*+1, x1 for *xn*, and x0 for *xn*−1. Note that x0 and x1 must be given (guessed) for the algorithm to start.

A program secant\_method.py that solves our example problem may be written as:

```
import sys
def secant(f, x0, x1, eps):
    f_x0 = f(x0)
    f_x1 = f(x1)
```
4 https://en.wikipedia.org/wiki/Secant\_line.

**Fig. 7.2** Illustrates the use of secants in the secant method when solving *x*2−9 = 0*, x* ∈ [0*,* 1000]. From two chosen starting values, *x*0 = 1000 and *x*1 = 700 the crossing *x*2 of the corresponding secant with the *x* axis is computed, followed by a similar computation of *x*3 from *x*1 and *x*2

```
iteration_counter = 0
    while abs(f_x1) > eps and iteration_counter < 100:
        try:
            denominator = (f_x1 - f_x0)/(x1 - x0)
            x = x1 - f_x1/denominator
        except ZeroDivisionError:
            print('Error! - denominator zero forx=', x)
            sys.exit(1) # Abort with error
        x0 = x1
        x1 = x
        f_x0 = f_x1
        f_x1 = f(x1)
        iteration_counter = iteration_counter + 1
    # Here, either a solution is found, or too many iterations
    if abs(f_x1) > eps:
        iteration_counter = -1
    return x, iteration_counter
if __name__ == '__main__':
    def f(x):
        return x**2 - 9
    x0 = 1000; x1 = x0 - 1
    solution, no_iterations = secant(f, x0, x1, eps=1.0e-6)
    if no_iterations > 0: # Solution found
```

```
print('Number of function calls: {:d}'.format(2+no_iterations))
    print('A solution is: {:f}'.format(solution))
else:
    print('Solution not found!')
```
The number of function calls is now related to no\_iterations, i.e., the number of iterations, as 2 + no\_iterations, since we need two function calls before entering the while loop, and then one function call per loop iteration. Note that, even though we need two points on the graph to compute each updated estimate, only a *single* function call (f(x1)) is required in each iteration since f(x0) becomes the "old" f(x1) and may simply be copied as f\_x0 = f\_x1 (the exception is the very first iteration where two function evaluations are needed).

Running secant\_method.py, gives the following printout on the screen:

```
Number of function calls: 19
A solution is: 3.000000
```
#### **7.4 The Bisection Method**

Neither Newton's method nor the secant method can guarantee that an existing solution will be found (see Exercises 7.1 and 7.2). The bisection method, however, does that. However, if there are several solutions present, it finds only one of them, just as Newton's method and the secant method. The bisection method is slower than the other two methods, so reliability comes with a cost of speed (but, again, for a single equation that is rarely an issue with laptops of today).

To solve *x*2 − 9 = 0, *x* ∈ [0*,* 1000], with the bisection method, we reason as follows. The first key idea is that if *f (x)* = *x*2 − 9 is *continuous* on the interval and the function values for the interval endpoints (*xL* = 0, *xR* = 1000) have *opposite signs*, *f (x) must* cross the *x* axis at least once on the interval. That is, we know there is at least one solution.

The second key idea comes from dividing the interval in two equal parts, one to the left and one to the right of the midpoint *xM* = 500. By evaluating the sign of *f (xM)*, we will immediately know whether a solution must exist to the left or right of *xM*. This is so, since if *f (xM)* ≥ 0, we know that *f (x)* has to cross the *x* axis between *xL* and *xM* at least once (using the same argument as for the original interval). Likewise, if instead *f (xM)* ≤ 0, we know that *f (x)* has to cross the *x* axis between *xM* and *xR* at least once.

In any case, we may proceed with half the interval only. The exception is if *f (xM)* ≈ 0, in which case a solution is found. Such interval halving can be continued until a solution is found. A "solution" in this case, is when |*f (xM)*| is sufficiently close to zero, more precisely (as before): |*f (xM)*| *<* , where is a small number specified by the user.

The sketched strategy seems reasonable, so let us write a reusable function that can solve a general algebraic equation *f (x)* = 0 (bisection\_method.py):

```
import sys
def bisection(f, x_L, x_R, eps):
    f_L = f(x_L)
    if f_L*f(x_R) > 0:
        print('Error! Function does not have opposite \
                 signs at interval endpoints!')
        sys.exit(1)
    x_M = (x_L + x_R)/2.0
    f_M = f(x_M)
    iteration_counter = 1
    while abs(f_M) > eps:
        if f_L*f_M > 0: # i.e. same sign
            x_L = x_M
            f_L = f_M
        else:
            x_R = x_M
        x_M = (x_L + x_R)/2
        f_M = f(x_M)
        iteration_counter = iteration_counter + 1
    return x_M, iteration_counter
if __name__ == '__main__':
    def f(x):
        return x**2 - 9
    a = 0; b = 1000
    solution, no_iterations = bisection(f, a, b, eps=1.0e-6)
    print('Number of function calls: {:d}'.format(1 + 2*no_iterations))
    print('A solution is: {:f}'.format(solution))
```
Note that we first check if *f* changes sign in [*a, b*], because that is a requirement for the algorithm to work. The algorithm also relies on a continuous *f (x)* function, but this is very challenging for a computer code to check.

We get the following printout to the screen when bisection\_method.py is run:

```
Number of function calls: 63
A solution is: 3.000000
```
We notice that the number of function calls is much higher than with the previous methods.

#### **Required Work in the Bisection Method**

If the starting interval of the bisection method is bounded by *a* and *b*, and the solution at step *n* is taken to be the middle value, the error is bounded as

$$\frac{|b-a|}{2^n},\tag{7.4}$$

because the initial interval has been halved *n* times. Therefore, to meet a tolerance , we need *n* iterations such that the length of the current interval equals :

$$\frac{|b-a|}{2^n} = \epsilon \quad \Rightarrow \quad n = \frac{\ln((b-a)/\epsilon)}{\ln 2}.$$

This is a great advantage of the bisection method: we know beforehand how many iterations *n* it takes to meet a certain accuracy in the solution.

#### **7.5 Rate of Convergence**

With the methods above, we noticed that the number of iterations or function calls could differ quite substantially. The number of iterations needed to find a solution is closely related to the *rate of convergence*, which dictates the speed of error reduction as we approach the root. More precisely, we introduce the error in iteration *n* as *en* = |*x* − *xn*|, and define the *convergence rate q* as

$$e\_{n+1} = Ce\_n^q,\tag{7.5}$$

where *C* is a constant. The exponent *q* measures how fast the error is reduced from one iteration to the next. The larger *q* is, the faster the error goes to zero (when *en <* 1), and the fewer iterations we need to meet the stopping criterion |*f (x)*| *<* .

#### **Convergence Rate and Iterations**

When we previously addressed numerical integration (Chap. 6), the approximation error *E* was related to the size *h* of the sub-intervals and the convergence rate *r* as *E* = *Khr*, *K* being some constant.

Observe that (7.5) gives a *different* definition of convergence rate. This makes sense, since numerical integration is based on a partitioning of the original integration interval into *n* sub-intervals, which is very different from the iterative procedures used here for solving nonlinear algebraic equations.

A single *q* in (7.5) is defined in the limit *n* → ∞. For finite *n*, and especially smaller *n*, *q* will vary with *n*. To estimate *q*, we can compute all the errors *en* and set up (7.5) for three consecutive experiments *n* − 1, *n*, and *n* + 1:

$$e\_n = Ce\_{n-1}^q,$$

$$e\_{n+1} = Ce\_n^q.$$

Dividing, e.g., the latter equation by the former, and solving with respect to *q*, we get that

$$q = \frac{\ln(e\_{n+1}/e\_n)}{\ln(e\_n/e\_{n-1})} \dots$$

Since this *q* will vary somewhat with *n*, we call it *qn*. As *n* grows, we expect *qn* to approach a limit (*qn* → *q*).

**Modifying Our Functions to Return All Approximations** To compute all the *qn* values, we need *all* the *xn* approximations. However, our previous implementations of Newton's method, the secant method, and the bisection method returned just the final approximation.

Therefore, we have modified our solvers5 accordingly, and placed them in nonlinear\_solvers.py. A user can choose whether the final value or the whole history of solutions is to be returned. Each of the extended implementations now takes an extra parameter return\_x\_list. This parameter is a boolean, set to True if the function is supposed to return all the root approximations, or False, if the function should only return the final approximation.

As an example, let us take a closer look at Newton:

```
def Newton(f, dfdx, x, eps, return_x_list=False):
    f_value = f(x)
    iteration_counter = 0
    if return_x_list:
        x_list = []
    while abs(f_value) > eps and iteration_counter < 100:
        try:
            x = x - float(f_value)/dfdx(x)
        except ZeroDivisionError:
            print('Error! - derivative zero for x = {:g}'.format(x))
            sys.exit(1) # Abort with error
        f_value = f(x)
        iteration_counter += 1
        if return_x_list:
            x_list.append(x)
    # Here, either a solution is found, or too many iterations
    if abs(f_value) > eps:
        iteration_counter = -1 # i.e., lack of convergence
```
5 An implemented numerical solution algorithm is often called a *solver*.

```
if return_x_list:
    return x_list, iteration_counter
else:
    return x, iteration_counter
```
We can now make a call

x, iter = Newton(f, dfdx, x=1000, eps=1e-6, return\_x\_list=True)

and get a list x returned. With knowledge of the exact solution *x* of *f (x)* = 0 we can compute all the errors *en* and all the associated *qn* values with the compact function (also found in nonlinear\_solvers.py)

```
def rate(x, x_exact):
    e = [abs(x_ - x_exact) for x_ in x]
    q = [log(e[n+1]/e[n])/log(e[n]/e[n-1])
         for n in range(1, len(e)-1, 1)]
    return q
```
The error model (7.5) works well for Newton's method and the secant method. For the bisection method, however, it works well in the beginning, but not when the solution is approached.

We can compute the rates *qn* and print them nicely (print\_rates.py),

```
def print_rates(method, x, x_exact):
    q = ['{:.2f}'.format(q_) for q_ in rate(x, x_exact)]
    print(method + ':')
    for q_ in q:
        print(q_, " ", end="") # end="" suppresses newline
```
The result for print\_rates('Newton', x, 3) is

Newton: 1.01 1.02 1.03 1.07 1.14 1.27 1.51 1.80 1.97 2.00

indicating that *q* = 2 is the rate for Newton's method. A similar computation using the secant method, gives the rates

```
secant:
1.26 0.93 1.05 1.01 1.04 1.05 1.08 1.13 1.20 1.30 1.43
1.54 1.60 1.62 1.62
```
Here it seems that *q* ≈ 1*.*6 is the limit.

**Remark** If we in the bisection method think of the length of the current interval containing the solution as the error *en*, then (7.5) works perfectly since *en*+1 = 1 2 *en*, i.e., *q* = 1 and *C* = 1 2 , but if *en* is the true error |*x*−*xn*|, it is easily seen from a sketch that this error can oscillate between the current interval length and a potentially very small value as we approach the exact solution. The corresponding rates *qn* fluctuate widely and are of no interest.

#### **7.6 Solving Multiple Nonlinear Algebraic Equations**

So far in this chapter, we have considered a single nonlinear algebraic equation. However, *systems* of such equations arise in a number of applications, foremost nonlinear ordinary and partial differential equations. Of the previous algorithms, only Newton's method is suitable for extension to systems of nonlinear equations.

#### **7.6.1 Abstract Notation**

Suppose we have *n* nonlinear equations, written in the following abstract form:

$$F0(\mathbf{x}\_0, \mathbf{x}\_1, \dots, \mathbf{x}\_n) = \mathbf{0},\tag{7.6}$$

$$F\_1(\mathbf{x}\_0, \mathbf{x}\_1, \dots, \mathbf{x}\_n) = \mathbf{0},\tag{7.7}$$

$$
\dot{\mathbf{i}} = \dot{\mathbf{i}}\tag{7.8}
$$

$$F\_n(\mathbf{x}\_0, \mathbf{x}\_1, \dots, \mathbf{x}\_n) = \mathbf{0} \,. \tag{7.9}$$

It will be convenient to introduce a *vector notation*

$$F = (F\_0, \dots, F\_{\mathbb{I}}), \quad \mathbf{x} = (\boldsymbol{\chi}\_0, \dots, \boldsymbol{\chi}\_n) \ .$$

The system can now be written as *F(x)* = **0**.

As a specific example on the notation above, the system

$$\mathbf{x}^2 = \mathbf{y} - \mathbf{x}\cos(\pi\mathbf{x})\tag{7.10}$$

$$
\mathbf{y}\mathbf{x} + e^{-\mathbf{y}} = \mathbf{x}^{-1} \tag{7.11}
$$

can be written in our abstract form by introducing *x*0 = *x* and *x*1 = *y*. Then

$$F\_0(\mathbf{x}\_0, \mathbf{x}\_1) = \mathbf{x}^2 - \mathbf{y} + \mathbf{x}\cos(\pi\mathbf{x}) = \mathbf{0},$$

$$F\_1(\mathbf{x}\_0, \mathbf{x}\_1) = \mathbf{y}\mathbf{x} + e^{-\mathbf{y}} - \mathbf{x}^{-1} = \mathbf{0}.$$

#### **7.6.2 Taylor Expansions for Multi-Variable Functions**

We follow the ideas of Newton's method for one equation in one variable: approximate the nonlinear *f* by a linear function and find the root of that function. When *n* variables are involved, we need to approximate a *vector function F(x)* by some linear function *F*˜ = *J x* + *c*, where *J* is an *n* × *n* matrix and *c* is some vector of length *n*.

The technique for approximating *F* by a linear function is to use the first two terms in a Taylor series expansion. Given the value of *F* and its partial derivatives with respect to *x* at some point *xi*, we can approximate the value at some point *xi*+1 by the two first term in a Taylor series expansion around *xi*:

$$F(\mathbf{x}\_{l+1}) \approx F(\mathbf{x}\_l) + \nabla F(\mathbf{x}\_l)(\mathbf{x}\_{l+1} - \mathbf{x}\_l) \ .$$

The next terms in the expansions are omitted here and of size ||*xi*+1 − *xi*||2, which are assumed to be small compared with the two terms above.

The expression ∇*F* is the matrix of all the partial derivatives of *F*. Component *(i, j )* in ∇*F* is

$$\frac{\partial F\_i}{\partial x\_j} \, .$$

For example, in our 2 × 2 system (7.10) and (7.11) we can use SymPy to compute the Jacobian:

```
In [1]: from sympy import *
In [2]: x0, x1 = symbols('x0 x1')
In [3]: F0 = x0**2 - x1 + x0*cos(pi*x0)
In [4]: F1 = x0*x1 + exp(-x1) - x0**(-1)
In [5]: diff(F0, x0)
Out[5]: -pi*x0*sin(pi*x0) + 2*x0 + cos(pi*x0)
In [6]: diff(F0, x1)
Out[6]: -1
In [7]: diff(F1, x0)
Out[7]: x1 + x0**(-2)
In [8]: diff(F1, x1)
Out[8]: x0 - exp(-x1)
```
We can then write

$$\nabla F = \begin{pmatrix} \frac{\partial F\_0}{\partial \mathbf{x}\_0} & \frac{\partial F\_0}{\partial \mathbf{x}\_1} \\ \frac{\partial F\_1}{\partial \mathbf{x}\_0} & \frac{\partial F\_1}{\partial \mathbf{x}\_1} \end{pmatrix} = \begin{pmatrix} 2\mathbf{x}\_0 + \cos(\pi \mathbf{x}\_0) - \pi \mathbf{x}\_0 \sin(\pi \mathbf{x}\_0) \ -1 \\ \mathbf{x}\_1 + \mathbf{x}\_0^{-2} & \mathbf{x}\_0 - e^{-\mathbf{x}\_1} \end{pmatrix}$$

The matrix ∇*F* is called the *Jacobian* of *F* and often denoted by *J*.

#### **7.6.3 Newton's Method**

The idea of Newton's method is that we have some approximation *xi* to the root and seek a new (and hopefully better) approximation *xi*+1 by approximating*F(xi*+1*)* by a linear function and solve the corresponding linear system of algebraic equations. We approximate the nonlinear problem *F(xi*+1*)* = 0 by the linear problem

$$F(\mathbf{x}\_{l}) + J(\mathbf{x}\_{l})(\mathbf{x}\_{l+1} - \mathbf{x}\_{l}) = \mathbf{0},\tag{7.12}$$

where *J(xi)* is just another notation for ∇*F(xi)*. The Eq. (7.12) is a linear system with coefficient matrix *J* and right-hand side vector *F(xi)*. We therefore write this system in the more familiar form

$$J(\mathbf{x}\_{i})\delta = -F(\mathbf{x}\_{i}),$$

where we have introduced a symbol *δ* for the unknown vector *xi*+1 − *xi* that multiplies the Jacobian *J*.

The *i*-th iteration of Newton's method for systems of algebraic equations consists of two steps:

Solving systems of linear equations must make use of appropriate software. Gaussian elimination is the most common, and in general the most robust, method for this purpose. Python's numpy package has a module linalg that interfaces the wellknown LAPACK package with high-quality and very well tested subroutines for linear algebra. The statement x = numpy.linalg.solve(A, b) solves a system *Ax* = *b* with a LAPACK method based on Gaussian elimination.

When nonlinear systems of algebraic equations arise from discretization of partial differential equations, the Jacobian is very often sparse, i.e., most of its elements are zero. In such cases it is important to use algorithms that can take advantage of the many zeros. Gaussian elimination is then a slow method, and (much) faster methods are based on iterative techniques.

#### **7.6.4 Implementation**

Here is a very simple implementation of Newton's method for systems of nonlinear algebraic equations:

```
import numpy as np
def Newton_system(F, J, x, eps):
    """
    Solve nonlinear system F=0 by Newton's method.
    J is the Jacobian of F. Both F and J must be functions of x.
    At input, x holds the start value. The iteration continues
    until ||F|| < eps.
    """
    F_value = F(x)
    F_norm = np.linalg.norm(F_value, ord=2) # l2 norm of vector
    iteration_counter = 0
    while abs(F_norm) > eps and iteration_counter < 100:
        delta = np.linalg.solve(J(x), -F_value)
        x = x + delta
        F_value = F(x)
        F_norm = np.linalg.norm(F_value, ord=2)
        iteration_counter = iteration_counter + 1
```

```
# Here, either a solution is found, or too many iterations
if abs(F_norm) > eps:
    iteration_counter = -1
return x, iteration_counter
```
We can test the function Newton\_systemwith the 2×2 system (7.10) and (7.11):

```
def test_Newton_system1():
    from numpy import cos, sin, pi, exp
    def F(x):
        return np.array(
            [x[0]**2 - x[1] + x[0]*cos(pi*x[0]),
             x[0]*x[1] + exp(-x[1]) - x[0]**(-1.)])
    def J(x):
        return np.array(
            [[2*x[0] + cos(pi*x[0]) - pi*x[0]*sin(pi*x[0]), -1],
             [x[1] + x[0]**(-2.), x[0] - exp(-x[1])]])
    expected = np.array([1, 0])
    tol = 1e-4
    x, n = Newton_system(F, J, x=np.array([2, -1]), eps=0.0001)
    print(n, x)
    error_norm = np.linalg.norm(expected - x, ord=2)
    assert error_norm < tol, 'norm of error ={:g}'.format(error_norm)
    print('norm of error ={:g}'.format(error_norm))
```
Here, the testing is based on the L2 norm6 of the error vector. Alternatively, we could test against the values of x that the algorithm finds, with appropriate tolerances. For example, as chosen for the error norm, if eps=0.0001, a tolerance of 10−4 can be used for x[0] and x[1].

#### **7.7 Exercises**

#### **Exercise 7.1: Understand Why Newton's Method Can Fail**

The purpose of this exercise is to understand when Newton's method works and fails. To this end, solve tanh *x* = 0 by Newton's method and study the intermediate details of the algorithm. Start with *x*0 = 1*.*08. Plot the tangent in each iteration of Newton's method. Then repeat the calculations and the plotting when *x*0 = 1*.*09. Explain what you observe.

Filename: Newton\_failure.\*.

#### **Exercise 7.2: See If the Secant Method Fails**

Does the secant method behave better than Newton's method in the problem described in Exercise 7.1? Try the initial guesses

1. *x*0 = 1*.*08 and *x*1 = 1*.*09

2. *x*0 = 1*.*09 and *x*1 = 1*.*1

6 https://en.wikipedia.org/wiki/Norm\_(mathematics)#Euclidean\_norm.

3. *x*0 = 1 and *x*1 = 2*.*3 4. *x*0 = 1 and *x*1 = 2*.*4

Filename: secant\_failure.\*.

#### **Exercise 7.3: Understand Why the Bisection Method Cannot Fail**

Solve the same problem as in Exercise 7.1, using the bisection method, but let the initial interval be [−5*,* 3]. Report how the interval containing the solution evolves during the iterations.

Filename: bisection\_nonfailure.\*.

#### **Exercise 7.4: Combine the Bisection Method with Newton's Method**

An attractive idea is to combine the reliability of the bisection method with the speed of Newton's method. Such a combination is implemented by running the bisection method until we have a narrow interval, and then switch to Newton's method for speed.

Write a function that implements this idea. Start with an interval [*a, b*] and switch to Newton's method when the current interval in the bisection method is a fraction *s* of the initial interval (i.e., when the interval has length *s(b*−*a)*). Potential divergence of Newton's method is still an issue, so if the approximate root jumps out of the narrowed interval (where the solution is known to lie), one can switch back to the bisection method. The value of *s* must be given as an argument to the function, but it may have a default value of 0.1.

Try the new method on tanh*(x)* = 0 with an initial interval [−10*,* 15]. Filename: bisection\_Newton.py.

#### **Exercise 7.5: Write a Test Function for Newton's Method**

The purpose of this function is to verify the implementation of Newton's method in the Newton function in the file nonlinear\_solvers.py. Construct an algebraic equation and perform two iterations of Newton's method by hand or with the aid of SymPy. Find the corresponding size of |*f (x)*| and use this as value for eps when calling Newton. The function should then also perform two iterations and return the same approximation to the root as you calculated manually. Implement this idea for a unit test as a test function test\_Newton().

Filename: test\_Newton.py.

#### **Exercise 7.6: Halley's Method and the Decimal Module**

A nonlinear algebraic equation *f (x)* = 0 may also be solved by Halley's method,7 given as:

$$\mathbf{x}\_{n+1} = \mathbf{x}\_n - \frac{2f(\mathbf{x}\_n)f'(\mathbf{x}\_n)}{2f'(\mathbf{x}\_n)^2 - f(\mathbf{x}\_n)f''(\mathbf{x}\_n)}, \quad n = 0, 1, \dots, \frac{1}{2}$$

with some starting value *x*0.

7 http://mathworld.wolfram.com/HalleysMethod.html.

The implementation of compute\_rates\_decimal should involve the decimal module (you search for the right documentation!), to better handle very small errors that may enter the rate computations. For comparison, you should also compute the rates without using the decimal module. Test and compare with several parameter values.

**Hint** The logarithms in the rate calculation might require some extra consideration when you use the decimal module. Filename: Halleys\_method.py.

#### **Exercise 7.7: Fixed Point Iteration**

A nonlinear algebraic equation *f (x)* = 0 may be solved in many different ways, and we have met some of these in this chapter. Another, very useful, solution approach is to first re-write the equation into *x* = *φ(x)* (this re-write is not unique), and then formulate the iteration

$$\alpha\_{n+1} = \phi(\alpha\_n), \quad n = 0, 1, \dots, \frac{1}{2}$$

with some starting value *x*0. If *φ(x)* is continuous, and if *φ(xn)* approaches *α* as *xn* approaches *α* (i.e., we get *α* = *φ(α)* as *n* → ∞), the iteration is called a *fixed point iteration* and *α* is referred to as a *fixed point* of the mapping *x* → *φ(x)*. Clearly, if a fixed point *α* is found, *α* will also be a solution to the original equation *f (x)* = 0.

In this exercise, we will briefly explore the fixed point iteration method by solving

$$x^3 + 2x = e^{-x}, \quad x \in [-2, 2] \text{ .} $$

For comparison, however, you will first be asked to solve the equation by Newton's method (which, in fact, can be seen as fixed point iteration8).

8 Check out https://en.wikipedia.org/wiki/Fixed\_point\_iteration.

When the program is executed, the original equation should be solved both with Newton's method and via fixed point iterations (as just described). Compare the output from the two methods.

Filename: fixed\_point\_iteration.py.

#### **Exercise 7.8: Solve Nonlinear Equation for a Vibrating Beam**

An important engineering problem that arises in a lot of applications is the vibrations of a clamped beam where the other end is free. This problem can be analyzed analytically, but the calculations boil down to solving the following nonlinear algebraic equation:

$$
\cosh \beta \cos \beta = -1,
$$

where *β* is related to important beam parameters through

$$
\beta^4 = \omega^2 \frac{\varrho A}{EI},
$$

where  is the density of the beam, *A* is the area of the cross section, *E* is Young's modulus, and *I* is the moment of the inertia of the cross section. The most important parameter of interest is *ω*, which is the frequency of the beam. We want to compute the frequencies of a vibrating steel beam with a rectangular cross section having width *b* = 25 mm and height *h* = 8 mm. The density of steel is 7850 kg/m3, and *E* = 2 × 1011 Pa. The moment of inertia of a rectangular cross section is *I* = *bh*3*/*12.

a) Plot the equation to be solved so that one can inspect where the zero crossings occur.

**Hint** When writing the equation as *f (β)* = 0, the *f* function increases its amplitude dramatically with *β*. It is therefore wise to look at an equation with damped amplitude, *g(β)* = *e*−*βf (β)* = 0. Plot *g* instead.

b) Compute the first three frequencies.

Filename: beam\_vib.py.

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **8 Solving Ordinary Differential Equations**

Differential equations constitute one of the most powerful mathematical tools to understand and predict the behavior of dynamical systems in nature, engineering, and society. A dynamical system is some system with some state, usually expressed by a set of variables, that evolves in time. For example, an oscillating pendulum, the spreading of a disease, and the weather are examples of dynamical systems. We can use basic laws of physics, or plain intuition, to express mathematical rules that govern the evolution of a system in time. These rules take the form of *differential equations*.

You are probably well experienced with equations, at least equations like *ax* + *b* = 0 and *ax*2 + *bx* + *c* = 0, where *a*, *b* and *c* are constants. Such equations are known as *algebraic equations*, and the unknowns are *numbers*. In a differential equation, the unknown is a *function*, and a differential equation will usually involve this function and one or more of its derivatives. When this function depends on a single independent variable, the equation is called an *ordinary differential equation* (ODE, plural: ODEs), which is different from a *partial differential equation* (PDE, plural: PDEs), in which the function depends on several independent variables (Chap. 9). As an example, *f (x)* = *f (x)* is a simple ODE (asking if there is any function *f* such that it equals its derivative—you might remember that *ex*

is a candidate).1 This is also an example of a *first-order* ODE, since the highest derivative appearing in the equation is a first derivative. When the highest derivative in an ODE is a second derivative, it is a *second-order* ODE, and so on (a similar terminology is used also for PDEs).

#### **Order of ODE versus order of numerical solution method**

Note that an ODE will have an *order* as just explained. This order, however, should not be confused with the *order* of a numerical solution method applied to solve that ODE. The latter refers to the convergence rate of the numerical solution method, addressed at the end of this chapter. We will present several such solution methods in this chapter, being first-order, second-order or fourth-order, for example (and a first-order ODE might, in principle, be solved by any of these methods).

The present chapter2 starts out preparing for ODEs and the Forward Euler method, which is a first-order method. Then we explain in detail how to solve ODEs numerically with the Forward Euler method, both *single (scalar)* first-order ODEs and *systems* of first-order ODEs. After the "warm-up" application—filling of a water tank—aimed at the less mathematically trained reader, we demonstrate all the mathematical and programming details through two specific applications: population growth and spreading of diseases. The first few programs we write, are deliberately made very simple and similar, while we focus the computational ideas.

Then we turn to oscillating mechanical systems, which arise in a wide range of engineering situations. The differential equation is now of second order, and the Forward Euler method does not perform too well. This observation motivates the need for other solution methods, and we derive the Euler-Cromer scheme, the second- and fourth-order Runge-Kutta schemes, as well as a finite difference scheme (the latter to handle the second-order differential equation directly without reformulating it as a first-order system). The presentation starts with undamped free oscillations and then treats general oscillatory systems with possibly nonlinear damping, nonlinear spring forces, and arbitrary external excitation. Besides developing programs from scratch, we also demonstrate how to access ready-made implementations of more advanced differential equation solvers in Python.

As we progress with more advanced methods, we develop more sophisticated and reusable programs. In particular, we incorporate good testing strategies, which allows us to bring solid evidence of correct computations. Consequently, the beginning—with water tank, population growth and disease modeling examples has a very gentle learning curve, while that curve gets significantly steeper towards the end of our section on oscillatory systems.

1 Note that the notation for the derivative may differ. For example, *f (x)* could equally well be written as just *f* (where the dependence on x is to be understood), or as *df dx* . 2 The reader should be aware of another excellent easy-to-read text by the late Prof. Langtangen,

"Finite Difference Computing with Exponential Decay Models" (Springer, open access, 2016, https://www.springer.com/gp/book/9783319294384). It fits right in with the material on ODEs and PDEs of the present book.

#### **8.1 Filling a Water Tank: Two Cases**

If "ordinary differential equation" is not among your favorite expressions, then this section is for you.

Consider a 25 L tank that will be filled with water in two different ways. In the first case, *the water volume that enters the tank per time* (rate of volume increase) is piecewise constant, while in the second case, it is continuously increasing.

For each of these two cases, we are asked to develop a code that can predict (i.e., compute) how the total water volume *V* in the tank will develop with time *t* over a period of 3 s. Our calculations must be based on the information given: the initial volume of water (1L in both cases), and the volume of water entering the tank per time.

#### **8.1.1 Case 1: Piecewise Constant Rate**

In this simpler case, there is initially 1 L of water in the tank, i.e.,

$$V(0) = 1\,\text{L},$$

while the rates of volume increase are given as:

$$r = 1\,\mathrm{L}\,\mathrm{s}^{-1}, \qquad 0\,\mathrm{s} < t < 1\,\mathrm{s},$$

$$r = 3\,\mathrm{L}\,\mathrm{s}^{-1}, \qquad 1\,\mathrm{s} \le t < 2\,\mathrm{s},$$

$$r = 7\,\mathrm{L}\,\mathrm{s}^{-1}, \qquad 2\,\mathrm{s} \le t \le 3\,\mathrm{s}.$$

Before turning to the programming, we should work out the exact solution by hand for this problem, since that is rather straight forward. Such a solution will of course be useful for verifying our implementation. In fact, comparing program output to these hand calculations should suffice for this particular problem.

**Exact Solution by Hand** Our reasoning goes like this: For each of the given subintervals (on the time axis), the total volume *V* of water in the tank will increase linearly. Thus, if we compute *V* after 1, 2 and 3 s, we will have what we need. We get

$$\begin{aligned} V(0) &= \text{l L}, \\ V(1) &= \text{l L} + (\text{l s})(\text{l L s}^{-1}) = 2\,\text{L}, \\ V(2) &= 2\,\text{L} + (\text{l s})(\text{3 L s}^{-1}) = 5\,\text{L}, \\ V(3) &= \text{S L} + (\text{l s})(7\,\text{L s}^{-1}) = 12\,\text{L}. \end{aligned}$$

We also have what is required for plotting the exact solution, since we can just tell Python to plot the computed *V* values against *t* for *t* = 0*,* 1*,* 2*,* 3, and let Python fill in the straight lines in between these points. Therefore, the exact solution is ready and we may proceed to bring our reasoning into code.

**Implementation and Performance** A simple implementation may read (rate\_piecewise\_constant.py):

```
import numpy as np
import matplotlib.pyplot as plt
a = 0.0; b = 3.0 # time interval
N = 3 # number of time steps
dt = (b - a)/N # time step (s)
V_exact = [1.0, 2.0, 5.0, 12.0] # exact volumes (L)
V = np.zeros(4) # numerically computed volume (L)
V[0] = 1 # initial volume
r = np.zeros(3) # rates of volume increase (L/s)
r[0] = 1; r[1] = 3; r[2] = 7
for i in [0, 1, 2]:
   V[i+1] = V[i] + dt*r[i]
time = [0, 1, 2, 3]
plt.plot(time, V, 'bo-', time, V_exact, 'r')
plt.title('Case 1')
plt.legend(['numerical','exact'], loc='upper left')
plt.xlabel('t (s)')
plt.ylabel('V (L)')
plt.show()
```
As you can see, we have included the exact (hand computed) solution in the code, so that it gets plotted together with the numerical solution found by the program. The time step dt will become 1 s here, and running the code, produces the plot seen in Fig. 8.1. We note that the numerical and exact solution can not be distinguished in the plot.

**Fig. 8.1** Water volume in a tank as it develops with piecewise constant rate of volume increase

#### **8.1.2 Case 2: Continuously Increasing Rate**

This case is more tricky. As in the previous case, there is 1 L of water in the tank from the start. When the tank fills up, however, the rate of volume increase, *r*, is *always* equal to the current volume of water, i.e., *r* = *V* (in units of L s−1). So, for example, at the time when there is 2 L in the tank, water enters the tank at 2 L s−1, when there is 2.1 L in the tank, water comes in at 2.1 L s−1, and so on. Thus, contrary to what we had in Case 1, *r* is not constant for any period of time within the 3 s interval, it increases continuously and gives us a steeper and steeper curve for *V (t)*. Writing this up as for Case 1, the information we have is

$$V(0) = 1\,\mathrm{L},$$

and, for the rate (in L s−1),

$$r(t) = V(t), \qquad 0 \le t \le 3 \,\text{s} \dots$$

Let us, for simplicity, assume that we also are given the exact solution in this case, which is *V (t)* = *et* . This allows us to easily check out the performance of any computational idea that we might try.

So, how can we compute the development of *V* , making it compare favorably to the given solution?

**An Idea** Clearly, we will be happy with an *approximately* correct solution, as long as the error can be made "small enough". In Case 1, we effectively computed connected straight line segments that matched the true development of *V* because of piecewise constant *r* values. Would it be possible to do something similar here in Case 2, i.e., compute straight line segments and use them as an *approximation* to the true solution curve? If so, it seems we could benefit from a very simple computational scheme! Let us pursue this idea further to see what comes out of it.

**The First Time Step** Considering the very first time step, we should realize that, since we are given the initial volume *V (*0*)* = 1 L, we do know the correct volume *and* correct rate at *t* = 0, since *r(t)* = *V (t)*. Thus, using this information and *pretending* that *r* stays constant as time increases, we will be able to compute a straight line segment for the very first time step (some *Δt* must be chosen). This straight line segment will then become tangent to the true solution curve when *t* = 0. The computed volume at the end of the first time step will have an error, but if our time step is not too large, the straight line segment will stay close to the true solution curve and the error in *V* should be "small".

**The Second Time Step** What about the second time step? Well, the volume we computed (with an error) at the end of the first time step, must now serve as the starting volume *and* ("constant") rate for the second time step. This allows us to compute an *approximation* to the volume also at the end of the second time step. If errors are "small", also the second straight line segment should be close to the true solution curve.

**What About the Errors?** We realize that, in this way, we can work our way all along the total time interval. Immediately, we suspect that the error may grow with the number of time steps, but since the total time interval is not too large, and since we may choose a *very* small time step on modern computers, this could still work!

**Implementation and Performance** Let us write down the code, which by choice gets very similar to the code in Case 1, and see how it performs. We realize that, for our strategy to work, the time steps should not be too large. However, during these initial investigations of ours, our aim is first and foremost to check out the computational idea. So, we pick a time step *Δt* = 0*.*1 s for a first try. A simple version of the code (rate\_exponential.py) may then read:

```
import numpy as np
import matplotlib.pyplot as plt
a = 0.0; b = 3.0 # time interval
N = 30 # number of time steps
dt = (b - a)/N # time step (s)
V = np.zeros(N+1) # numerically computed volume (L)
V[0] = 1 # initial volume
for i in range(0, N, 1):
   V[i+1] = V[i] + dt*V[i] # ...r is V now
time_exact = np.linspace(a, b, 1000)
V_exact = np.exp(time_exact) # make exact solution (for plotting)
time = np.linspace(0, 3, N+1)
plt.plot(time, V, 'bo-', time_exact, V_exact, 'r')
plt.title('Case 2')
plt.legend(['numerical','exact'], loc='upper left')
plt.xlabel('t (s)')
plt.ylabel('V (L)')
plt.show()
```
To plot the exact solution, we just picked 1000 points in time, which we consider "large enough" to get a representative curve. Compared to the code for Case 1, some more flexibility is introduced here, using range and N in the for loop header. Running the code gives the plot shown in Fig. 8.2.

This looks promising! Not surprisingly, the error grows with time, reaching about 2*.*64 L at the end. However, the time step is not particularly small, so we should expect much more accurate computations if *Δt* is reduced. We skip showing the plots,3 but if we increase N from 30 to 300, the maximum error drops to 0*.*30 L, while an N value of 3 · 106 gives an error of 3 · 10−5 L. It seems we are on to something!

3 With smaller time steps, it becomes inappropriate to use filled circles on the graph for the numerical values. Thus, in the plot command, one should change bo- to, e.g., only b.

**Fig. 8.2** Water volume in a tank as it develops with constantly increasing rate of volume change (*r* = *V* ). The numerical computations use piecewise constant rates as an approximation to the true rate

The similar code structures used for Case 1 and Case 2 suggests that a more general (and improved) code may be written, applicable to both cases. As we move on with our next example, population growth, we will see how this can be done.

#### **8.1.3 Reformulating the Problems as ODEs**

Typically, the problems we solved in Case 1 and Case 2, would rather have been presented in "more proper mathematical language" as ODEs.

**Case 1** When the rates were piecewise constant, we could have been requested to solve

$$V'(t) = 1\,\mathrm{L}\,\mathrm{s}^{-1}, \qquad 0\,\mathrm{s} < t < 1\,\mathrm{s},$$

$$V'(t) = 3\,\mathrm{L}\,\mathrm{s}^{-1}, \qquad 1\,\mathrm{s} \le t < 2\,\mathrm{s},$$

$$V'(t) = 7\,\mathrm{L}\,\mathrm{s}^{-1}, \qquad 2\,\mathrm{s} \le t \le 3\,\mathrm{s},$$

with

$$V(0) = 1\,\text{L},$$

where *V (*0*)* is known as an *initial condition*.

**Case 2** With a continuously increasing rate, we could have been asked to solve

$$V'(t) = V(t), \qquad V(0) = 1\,\mathrm{L}, \qquad 0\,\mathrm{s} < t \le \Im\,\mathrm{s} \,\mathrm{L}$$

This particular ODE is very similar to the ODE we will address when we next turn to population growth (in fact, it may be seen as a special case of the latter).

**The Forward Euler Method: A Brief Encounter** If we had proceeded to solve these ODEs by something called the *Forward Euler method* (or *Euler's method*), we could (if we wanted) have written the solution codes exactly as they were developed above! Thus, we have already used the essentials of Euler's method without stating it.

In the following sections, the Forward Euler method will be thoroughly explained and elaborated on, while we demonstrate how the approach may be used to solve different ODEs numerically.

#### **8.2 Population Growth: A First Order ODE**

Our first real taste of differential equations regards modeling the growth of some population, such as a cell culture, an animal population, or a human population. The ideas even extend trivially to growth of money in a bank.

Let *N(t)* be the number of individuals in the population at time *t*. How can we predict the evolution of *N* with time? Below we shall derive a differential equation whose solution is *N(t)*. The equation we will derive reads

$$N'(t) = rN(t),\tag{8.1}$$

where *r* is a number. Note that although *N* obviously is an integer in real life, we model *N* as a real-valued function. We choose to do this, because the solutions of differential equations are (normally continuous) real-valued functions. An integervalued *N(t)* in the model would lead to a lot of mathematical difficulties. Also, talking about, e.g., 2.5 individuals is no problem in mathematics, even though we must be a bit careful when applying this in a practical setting!

You may know, or find out, that the solution *N(t)* = *Cert*, where *C* is any number. To make this solution unique, we need to fix *C*, which is done by prescribing the value of *N* at some time, usually at *t* = 0. If *N(*0*)* is given as *N*0, we get *N(t)* = *N*0*ert* .

In general, a *differential equation model* consists of a *differential equation*, such as (8.1) *and* an *initial condition*, such as *N(*0*)* = *N*0. With a known initial condition, the differential equation can be solved for the unknown function and the solution is unique.

It is very rare that we can find the solution of a differential equation as easy as the ODE in this example allows. Normally, one has to apply certain mathematical methods. Still, these methods can only handle some of the simplest differential equations. However, with numerical methods and a bit of programming, we can easily deal with almost any differential equation! This is exactly the topic of the present chapter.

#### **8.2.1 Derivation of the Model**

It can be instructive to show how an equation like (8.1) arises. Consider some population of an animal species and let *N(t)* be the number of individuals in a certain spatial region, e.g. an island. We are not concerned with the spatial distribution of the animals, just the number of them in some region where there is no exchange of individuals with other regions. During a time interval *Δt*, some animals will die and some will be born. The numbers of deaths and births are expected to be proportional to *N*. For example, if there are twice as many individuals, we expect them to get twice as many newborns. In a time interval *Δt*, the net growth of the population will then be

$$N(t + \Delta t) - N(t) = bN(t) - dN(t),$$

where *bN(t)* ¯ is the number of newborns and *dN(t)* ¯ is the number of deaths. If we double *Δt*, we expect the proportionality constants *b*¯ and *d*¯ to double too, so it makes sense to think of *b*¯ and *d*¯ as proportional to *Δt* and "factor out" *Δt*. That is, we introduce *b* = *b/Δt* ¯ and *d* = *d/Δt* ¯ to be proportionality constants for newborns and deaths independent of *Δt*. Also, we introduce *r* = *b* − *d*, which is the net rate of growth of the population per time unit. Our model then becomes

$$N(t + \Delta t) - N(t) = \Delta t \, rN(t) \,. \tag{8.2}$$

Equation (8.2) is actually a computational model. Given *N(t)*, we can advance the population size by

$$N(t + \Delta t) = N(t) + \Delta t \, rN(t) \, .$$

This is called a *difference equation*. If we know *N(t)* for some *t*, e.g., *N(*0*)* = *N*0, we can compute

$$N(\Delta t) = N\_0 + \Delta t \, rN\_0,$$

$$N(2\Delta t) = N(\Delta t) + \Delta t \, rN(\Delta t),$$

$$N(3\Delta t) = N(2\Delta t) + \Delta t \, rN(2\Delta t),$$

$$\vdots$$

$$N((k+1)\Delta t) = N(k\Delta t) + \Delta t \, rN(k\Delta t),$$

where *k* is some arbitrary integer. A computer program can easily compute *N((k* + 1*)Δt)* for us with the aid of a little loop.

#### **The initial condition**

Observe that the computational formula cannot be started unless we have an initial condition!

The solution of *N* = *rN* is *N* = *Cert* for any constant *C*, and the initial condition is needed to fix *C* so the solution becomes unique. However, from a mathematical point of view, knowing *N(t)* at any point *t* is sufficient as initial condition. Numerically, we more literally need an initial condition: we need to know a starting value at the left end of the interval in order to get the computational formula going.

In fact, we do not really need a computer in this particular case, since we see a repetitive pattern when doing hand calculations. This leads us to a mathematical formula for *N((k* + 1*)Δt)*:

$$N((k+1)\Delta t) = N(k\Delta t) + \Delta t \, rN(k\Delta t) = N(k\Delta t)(1 + \Delta t \, r)$$

$$= N((k-1)\Delta t)(1 + \Delta t \, r)^2$$

$$\vdots$$

$$= N\_0(1 + \Delta t \, r)^{k+1}.$$

Rather than using (8.2) as a computational model directly, there is a strong tradition for deriving a differential equation from this difference equation. The idea is to consider a very small time interval *Δt* and look at the instantaneous growth as this time interval is shrunk to an infinitesimally small size. In mathematical terms, it means that we let *Δt* → 0. As (8.2) stands, letting *Δt* → 0 will just produce an equation 0 = 0, so we have to divide by *Δt* and then take the limit:

$$\lim\_{\Delta t \to 0} \frac{N(t + \Delta t) - N(t)}{\Delta t} = rN(t) \dots$$

The term on the left-hand side is actually the definition of the derivative *N (t)*, so we have

$$N'(t) = rN(t),$$

which is the corresponding differential equation.

There is nothing in our derivation that forces the parameter *r* to be constant—it can change with time due to, e.g., seasonal changes or more permanent environmental changes.

#### **Detour: Exact mathematical solution**

If you have taken a course on mathematical solution methods for differential equations, you may want to recap how an equation like *N* = *rN* or *N* = *r(t)N* is solved. The *method of separation of variables* is the most convenient solution strategy in this case:

$$N' = rN$$

$$\frac{dN}{dt} = rN$$

$$\frac{dN}{N} = rdt$$

$$\int\_{N\_0}^{N} \frac{dN}{N} = \int\_0^t rdt$$

$$\ln N - \ln N\_0 = \int\_0^t r(t)dt$$

$$N = N\_0 \exp\left(\int\_0^t r(t)dt\right),$$

which for constant *r* results in *N* = *N*0*ert*. Note that exp *(t)* is the same as *et* .

As will be described later, *r* must in more realistic models depend on *N*. The method of separation of variables then requires to integrate - *N N*0 *N/r(N)dN*, which quickly becomes non-trivial for many choices of *r(N)*. The only generally applicable solution approach is therefore a numerical method.

#### **8.2.2 Numerical Solution: The Forward Euler (FE) Method**

There is a huge collection of numerical methods for problems like (8.1), and in general any equation of the form *u* = *f (u, t)*, where *u(t)* is the unknown function in the problem, and *f* is some known formula of *u* and optionally *t*. In our case with population growth, i.e., (8.1), *u (t)* corresponds to *N (t)*, while *f (u, t)* corresponds to *rN(t)*.

We will first present a simple *finite difference method* solving *u* = *f (u, t)*. The idea is fourfold:

**Fig. 8.3** Mesh in time with corresponding discrete values (unknowns)

**Fig. 8.4** Illustration of a forward difference approximation to the derivative

The collection of points *t*0*, t*1*,...,tNt* is called a *mesh* (or *grid*), while the corresponding *un* values4 collectively are referred to as a *mesh function*.

Let us use an example to illustrate the given steps. First, we introduce the mesh, which often is *uniform*, meaning that the spacing between mesh points *tn* and *tn*+1

4 To get an estimate of *u* in between the mesh points, one often assumes a straight line relationship between computed *un* values. Check out *linear interpolation*, e.g., on Wikipedia (https://en.wikipedia.org/wiki/Linear\_interpolation).

is constant. This property implies that

$$t\_n = n\Delta t, \quad n = 0, 1, \ldots, N\_I \dots$$

Second, the differential equation is supposed to hold at the mesh points. Note that this is an approximation, because the differential equation is originally valid for all real values of *t*. We can express this property mathematically as

$$\mu'(t\_n) = f(\mu'', t\_n), \quad n = 0, 1, \dots, N\_l \dots$$

For example, with our model equation *u* = *ru*, we have the special case

$$u'(t\_n) = ru^n, \quad n = 0, 1, \ldots, N\_l, 1$$

or

$$
\mu'(t\_n) = r(t\_n)\mu'', \quad n = 0, 1, \dots, N\_{\mathcal{I}}, \\
$$

if *r* depends explicitly on *t*.

Third, derivatives are to be replaced by finite differences. To this end, we need to know specific formulas for how derivatives can be approximated by finite differences. One simple possibility is to use the definition of the derivative from any calculus book,

$$\mu'(t) = \lim\_{\Delta t \to 0} \frac{\mu(t + \Delta t) - \mu(t)}{\Delta t}.$$

At an arbitrary mesh point *tn* this definition can be written as

$$\mu'(t\_n) = \lim\_{\Delta t \to 0} \frac{\mu^{n+1} - \mu^n}{\Delta t}.$$

Instead of going to the limit *Δt* → 0 we can use a small *Δt*, which yields a computable approximation to *u (tn)*:

$$
\mu'(t\_n) \approx \frac{\mu^{n+1} - \mu^n}{\Delta t} \dots
$$

This is known as a *forward difference* since we go forward in time (*un*+1) to collect information in *u* to estimate the derivative. Figure 8.4 illustrates the idea. The error of the forward difference is proportional to *Δt* (often written as *O(Δt)*, but we will not use this notation in the present book).

We can now plug in the forward difference in our differential equation sampled at the arbitrary mesh point *tn*:

$$\frac{\mu^{n+1} - \mu^n}{\Delta t} = f(\mu^n, t\_n), \tag{8.3}$$

or with *f (u, t)* = *ru* in our special model problem for population growth,

$$\frac{u^{n+1} - u^n}{\Delta t} = ru^n. \tag{8.4}$$

If *r* depends on time, we insert *r(tn)* = *rn* for *r* in this latter equation.

The fourth step is to derive a computational algorithm. Looking at (8.3), we realize that if *un* should be known, we can easily solve with respect to *un*+1 to get a formula for *u* at the next time level *tn*+1:

$$
\mu^{n+1} = \mu^n + \Delta t f(\mu^n, t\_n) \,. \tag{8.5}
$$

Provided we have a known starting value, *u*0 = *U*0, we can use (8.5) to advance the solution by first computing *u*1 from *u*0, then *u*2 from *u*1, *u*3 from *u*2, and so forth.

Such an algorithm is called a *numerical scheme* for the differential equation. It is often written compactly as

$$u^{n+1} = u^n + \Delta t f(u^n, t\_n), \quad u^0 = U\_0, \quad n = 0, 1, \ldots, N\_l - 1. \tag{8.6}$$

This scheme is known as the *Forward Euler scheme*, also called *Euler's method*.

In our special population growth model, we have

$$u^{n+1} = u^n + \Delta t \, ru^n, \quad u^0 = U\_0, \quad n = 0, 1, \ldots, N\_l - 1. \tag{8.7}$$

We may also write this model using the problem-specific symbol *N* instead of the generic *u* function:

$$N^{n+1} = N^n + \Delta t \, rN^n, \quad N^0 = N\_0, \quad n = 0, 1, \ldots, N\_l - 1. \tag{8.8}$$

The observant reader will realize that (8.8) is nothing but the computational model (8.2) arising directly in the model derivation. The formula (8.8) arises, however, from a detour via a differential equation and a numerical method for the differential equation. This looks rather unnecessary! The reason why we bother to derive the differential equation model and then discretize it by a numerical method is simply that the discretization can be done in many ways, and we can create (much) more accurate and more computationally efficient methods than (8.8) or (8.6). This can be useful in many problems! Nevertheless, the Forward Euler scheme is intuitive and widely applicable, at least when *Δt* is chosen to be small.

#### **The numerical solution between the mesh points**

Our numerical method computes the unknown function *u* at discrete mesh points *t*1*, t*2*,...,tNt* . What if we want to evaluate the numerical solution between the mesh points? The most natural choice is to *assume* a linear variation between the mesh points, see Fig. 8.5. This is compatible with the fact that when we plot the array *u*0*, u*1*,...* versus *t*0*, t*1*,...*, a straight line is automatically drawn between the discrete points (unless we specify otherwise in the plot command).

**Fig. 8.5** The numerical solution at points can be extended by linear segments between the mesh points

#### **8.2.3 Programming the FE Scheme; the Special Case**

Let us compute (8.8) in a program. The input variables are *N*0, *Δt*, *r*, and *Nt* . Note that we need to compute *Nt* new values *N*1*,...,NNt* . A total of *Nt* + 1 values are needed in an array representation of *Nn*, *n* = 0*,...,Nt* .

Our first version of this program (growth1.py) is as simple as possible, and very similar to the codes we wrote previously for the water tank example:

```
import numpy as np
import matplotlib.pyplot as plt
N_0 = int(input('Give initial population size N_0: '))
r = float(input('Give net growth rate r: '))
dt = float(input('Give time step size: '))
N_t = int(input('Give number of steps: '))
t = np.linspace(0, N_t*dt, N_t+1)
N = np.zeros(N_t+1)
N[0] = N_0
for n in range(N_t):
    N[n+1] = N[n] + r*dt*N[n]
numerical_sol = 'bo' if N_t < 70 else 'b-'
plt.plot(t, N, numerical_sol, t, N_0*np.exp(r*t), 'r-')
plt.legend(['numerical', 'exact'], loc='upper left')
plt.xlabel('t'); plt.ylabel('N(t)')
filestem = 'growth1_{:d}steps'.format(N_t)
plt.savefig('{:s}.png'.format(filestem))
plt.savefig('{:s}.pdf'.format(filestem))
```
Note the compact assignment statement

numerical\_sol = 'bo' if N\_t < 70 else 'b-'

This is a one-line alternative to, e.g.,

```
if N_t < 70:
    numerical_sol = 'bo'
else:
    numerical_sol = 'b-'
```
Let us demonstrate a simulation where we start with 100 animals, a net growth rate of 10% (0.1) per time unit, which can be 1 month, and *t* ∈ [0*,* 20] months. We may first try *Δt* of half a month (0.5), which implies *Nt* = 40. Figure 8.6 shows the results. The solid line is the exact solution, while the circles are the computed numerical solution. The discrepancy is clearly visible. What if we make *Δt* 10 times smaller? The result is displayed in Fig. 8.7, where we now use a solid line also for the numerical solution (otherwise, 400 circles would look very cluttered, so the program has a test on how to display the numerical solution, either as circles or a solid line). We can hardly distinguish the exact and the numerical solution. The computing time is also a fraction of a second on a laptop, so it appears that the Forward Euler method is sufficiently accurate for practical purposes. (This is not always true for large, complicated simulation models in engineering, so more sophisticated methods may be needed.)

It is also of interest to see what happens if we increase *Δt* to 2 months. The results in Fig. 8.8 indicate that this is an inaccurate computation.

**Fig. 8.6** Evolution of a population computed with time step 0.5 month

**Fig. 8.7** Evolution of a population computed with time step 0.05 month

**Fig. 8.8** Evolution of a population computed with time step 2 months

#### **8.2.4 Understanding the Forward Euler Method**

The good thing about the Forward Euler method is that it gives an understanding of what a differential equation is and a geometrical picture of how to construct the solution. The first idea is that we have already computed the solution up to some time point *tn*. The second idea is that we want to progress the solution from *tn* to *tn*+1 as a straight line.

We know that the line must go through the solution at *tn*, i.e., the point *(tn, un)*. The differential equation tells us the slope of the line: *u (tn)* = *f (un, tn)* = *run*. That is, the differential equation gives a direct formula for the further *direction* of the solution curve. We can say that the differential equation expresses how the system (*u*) undergoes changes at a point.

There is a general formula for a straight line *y* = *ax* + *b* with slope *a* that goes through the point *(x*0*, y*0*)*: *y* = *a(x* − *x*0*)* + *y*0. Using this formula adapted to the present case, and evaluating the formula for *tn*+1, results in

$$u^{n+1} = ru^n(t\_{n+1} - t\_n) + u^n = u^n + \Delta t \, ru^n,$$

which is nothing but the Forward Euler formula. You are now encouraged to do Exercise 8.2 to become more familiar with the geometric interpretation of the Forward Euler method.

#### **8.2.5 Programming the FE Scheme; the General Case**

Our previous program was just a simple main program without function definitions, tailored to a special differential equation. When programming mathematics, it is always good to consider a (large) class of problems and making a Python function to solve any problem that fits into the class. More specifically, we will make software for the class of differential equation problems of the form

$$u'(t) = f(u, t), \quad u = U\_0, \ t \in [0, T], t$$

for some given function *f* , and numbers *U*0 and *T* . We also take the opportunity to illustrate what is commonly called a demo function. As the name implies, the purpose of such a function is solely to demonstrate how the function works (not to be confused with a test function, which does verification by use of assert). The Python function calculating the solution must take *f* , *U*0, *Δt*, and *T* as input, find the corresponding *Nt* , compute the solution, and return an array with *u*0*, u*1*,...,uNt* and an array with *t*0*, t*1*,...,tNt* . The Forward Euler scheme reads

$$
\mu^{n+1} = \mu^n + \Delta t f(\mu^n, t\_n), \quad n = 0, \dots, N\_l - 1 \dots
$$

The corresponding program may now take the form (file ode\_FE.py):

```
import numpy as np
import matplotlib.pyplot as plt
def ode_FE(f, U_0, dt, T):
    N_t = int(round(T/dt))
    u = np.zeros(N_t+1)
    t = np.linspace(0, N_t*dt, len(u))
    u[0] = U_0
```

```
for n in range(N_t):
        u[n+1] = u[n] + dt*f(u[n], t[n])
    return u, t
def demo_population_growth():
    """Test case: u'=r*u, u(0)=100."""
    def f(u, t):
        return 0.1*u
    u, t = ode_FE(f=f, U_0=100, dt=0.5, T=20)
    plt.plot(t, u, t, 100*np.exp(0.1*t))
    plt.show()
if __name__ == '__main__':
    demo_population_growth()
```
This program file, called ode\_FE.py, is a reusable piece of code with a general ode\_FE function that can solve any differential equation *u* = *f (u, t)* and a demo function for the special case *u* = 0*.*1*u*, *u(*0*)* = 100. Observe that the call to the demo function is placed in a test block. This implies that the call is not active if ode\_FE is imported as a module in another program, but active if ode\_FE.py is run as a program.

The solution should be identical to what the growth1.py program produces with the same parameter settings (*r* = 0*.*1, *N*0 = 100). This feature can easily be tested by inserting a print command, but a much better, automated verification is suggested in Exercise 8.2. You are strongly encouraged to take a "break" and do that exercise now.

#### **Remark on the Use of** u **as Variable**

In the ode\_FE program, the variable u is used in different contexts. Inside the ode\_FE function, u is an array, but in the f(u,t) function, as exemplified in the demo\_population\_growth function, the argument u is a number. Typically, we call f (in ode\_FE) with the u argument as one element of the array u in the ode\_FE function: u[n].

#### **8.2.6 A More Realistic Population Growth Model**

Exponential growth of a population according the model *N* = *rN*, with exponential solution *N* = *N*0*ert*, is unrealistic in the long run because the resources needed to feed the population are finite. At some point there will not be enough resources and the growth will decline. A common model taking this effect into account assumes that *r* depends on the size of the population, *N*:

$$N(t + \Delta t) - N(t) = r(N(t))N(t)\dots$$

The corresponding differential equation becomes

$$N' = r(N)N\dots$$

The reader is strongly encouraged to repeat the steps in the derivation of the Forward Euler scheme and establish that we get

$$N^{n+1} = N^n + \Delta t \, r(N^n) N^n,$$

which computes as easy as for a constant *r*, since *r(Nn)* is known when computing *Nn*+1. Alternatively, one can use the Forward Euler formula for the general problem *u* = *f (u, t)* and use *f (u, t)* = *r(u)u* and replace *u* by *N*.

The simplest choice of *r(N)* is a linear function, starting with some growth value *r*¯ and declining until the population has reached its maximum, *M*, according to the available resources:

$$r(N) = \bar{r}(1 - N/M)\dots$$

In the beginning, *N M* and we will have exponential growth *ert*¯ , but as *N* increases, *r(N)* decreases, and when *N* reaches *M*, *r(N)* = 0 so there is no more growth and the population remains at *N(t)* = *M*. This linear choice of *r(N)* gives rise to a model that is called the *logistic model*. The parameter *M* is known as the *carrying capacity* of the population.

Let us run the logistic model with aid of the ode\_FE function. We choose *N(*0*)* = 100, *Δt* = 0*.*5 month, *T* = 60 months, *r* = 0*.*1, and *M* = 500. The complete program, called logistic.py, is basically a call to ode\_FE:

```
from ode_FE import ode_FE
import matplotlib.pyplot as plt
for dt, T in zip((0.5, 20), (60, 100)):
    u, t = ode_FE(f=lambda u, t: 0.1*(1 - u/500.)*u, \
                               U_0=100, dt=dt, T=T)
    plt.figure() # Make separate figures for each pass in the loop
    plt.plot(t, u, 'b-')
    plt.xlabel('t'); plt.ylabel('N(t)')
    plt.savefig('tmp_{:g}.png'.format(dt))
    plt.savefig('tmp_{:g}.pdf'.format(dt))
```
Figure 8.9 shows the resulting curve. We see that the population stabilizes around *M* = 500 individuals. A corresponding exponential growth would reach *N*0*ert* = 100*e*0*.*1·60 ≈ 40*,* 300 individuals!

What happens if we use "large" *Δt* values here? We may set *Δt* = 20 and *T* = 100. Now the solution, seen in Fig. 8.10, oscillates and is hence qualitatively wrong, because one can prove that the exact solution of the differential equation is monotone.

**Fig. 8.9** Logistic growth of a population

**Fig. 8.10** Logistic growth with large time step

#### **Remark on the world population**

The number of people on the planet (http://en.wikipedia.org/wiki/Population\_ growth) follows the model *N* = *r(t)N*, where the net reproduction *r(t)* varies with time and has decreased since its top in 1990. The current world value of *r* is 1.2%, and it is difficult to predict future values. At the moment, the predictions of the world population point to a growth to 9.6 billion before declining.

This example shows the limitation of a differential equation model: we need to know all input parameters, including *r(t)*, in order to predict the future. It is seldom the case that we know all input parameters. Sometimes knowledge of the solution from measurements can help estimate missing input parameters.

#### **8.2.7 Verification: Exact Linear Solution of the Discrete Equations**

How can we verify that the programming of an ODE model is correct? One way, is to compute convergence rates with a solver and confirm that the rates are according to expectations. We address convergence rates for ODE solvers later (in Sect. 8.5) and will then show how a corresponding test function for the ode\_FE solver may be written.

The best verification method, however, is to find a problem where there are no unknown numerical approximation errors, because we can then compare the exact solution of the problem with the result produced by our implementation and expect the difference to be within a very small tolerance. We shall base a unit test on this idea and implement a corresponding *test function* (see Sect. 6.6.4) for automatic verification of our implementation.

It appears that most numerical methods for ODEs will exactly reproduce a solution *u* that is linear in *t*. We may therefore set *u* = *at* + *b* and choose any *f* whose derivative is *a*. The choice *f (u, t)* = *a* is very simple, but we may add anything that is zero, e.g.,

$$f(u,t) = a + (u - (at+b))^m.$$

This is a valid *f (u, t)* for any *a*, *b*, and *m*. The corresponding ODE looks highly non-trivial, however:

$$
\mu' = a + (\mu - (at + b))^m .
$$

Using the general ode\_FE function in ode\_FE.py, we may write a proper test function as follows (in file test\_ode\_FE\_exact\_linear.py):

```
def test_ode_FE():
    """Test that a linear u(t)=a*t+b is exactly reproduced."""
```

```
def exact_solution(t):
    return a*t + b
def f(u, t): # ODE
    return a + (u - exact_solution(t))**m
a = 4
b = -1
m = 6
dt = 0.5
T = 20.0
u, t = ode_FE(f, exact_solution(0), dt, T)
diff = abs(exact_solution(t) - u).max()
tol = 1E-15 # Tolerance for float comparison
success = diff < tol
assert success
```
As a measure of the error, we have here simply used the maximum error (picked out by a call to max and assigned to diff).

Recall that test functions should start with the name test\_, have no arguments, and formulate the test as a boolean expression success that is True if the test passes and False if it fails. Test functions should make the test as assert success (here success can also be a boolean expression as in assert diff < tol).

Observe that we cannot compare diff to zero, which is what we mathematically expect, because diff is a floating-point variable that most likely contains small rounding errors. Therefore, we must compare diff to zero with a tolerance, here 10−15.

You are encouraged to do Exercise 8.3 where the goal is to make a test function for a verification based on comparison with hand-calculated results for a few time steps.

#### **8.3 Spreading of Disease: A System of First Order ODEs**

Our aim with this section is to show in detail how one can apply mathematics and programming to solve a system of first-order ODEs. We will do this as we investigate the spreading of disease. The mathematical model is now a system of three differential equations with three unknown functions. To derive such a model, we can use mainly intuition, so no specific background knowledge of diseases is required.

#### **8.3.1 Spreading of Flu**

Imagine a boarding school out in the country side. This school is a small and closed society. Suddenly, one or more of the pupils get the flu. We expect that the flu may spread quite effectively or die out. The question is how many of the pupils and the school's staff will be affected. Some quite simple mathematics can help us to achieve insight into the dynamics of how the disease spreads.

Let the mathematical function *S(t)* count how many individuals, at time *t*, that have the possibility to get infected. Here, *t* may count hours or days, for instance. These individuals make up a category called susceptibles, labeled as S. Another category, I, consists of the individuals that are infected. Let *I (t)* count how many there are in category I at time *t*. An individual having recovered from the disease is assumed to gain immunity.  In either case, the individual is moved from the I category to a category we call the removed category, labeled with R. We let *R(t)* count the number of individuals in the *R* category at time *t*. Those who enter the *R* category, cannot leave this category.

To summarize, the spreading of this disease is essentially the dynamics of moving individuals from the S to the I and then to the R category:

We can use mathematics to more precisely describe the exchange between the categories. The fundamental idea is to describe the changes that take place during a small time interval, denoted by *Δt*.

Our disease model is often referred to as a *compartment model*, where quantities are shuffled between compartments (here a synonym for categories) according to some rules. The rules express *changes* in a small time interval *Δt*, and from these changes we can let *Δt* go to zero and obtain derivatives. The resulting equations then go from difference equations (with finite *Δt*) to differential equations (*Δt* → 0).

We introduce a uniform mesh in time, *tn* = *nΔt*, *n* = 0*,...,Nt* , and seek *S* at the mesh points. The numerical approximation to *S* at time *tn* is denoted by *Sn*. Similarly, we seek the unknown values of *I (t)* and *R(t)* at the mesh points and introduce a similar notation *I n* and *Rn* for the approximations to the exact values *I (tn)* and *R(tn)*.

In the time interval *Δt* we know that some people will be infected, so *S* will decrease. We shall soon argue by mathematics that there will be *βΔtSI* new infected individuals in this time interval, where *β* is a parameter reflecting how easy people get infected during a time interval of unit length. If the loss in *S* is *βΔtSI* , we have that the change in *S* is

$$S^{n+1} - S^n = -\beta \,\Delta t \, S^n \, I^n \,. \tag{8.9}$$

Dividing by *Δt* and letting *Δt* → 0, makes the left-hand side approach *S (tn)* such that we obtain a differential equation

$$S' = -\beta SI.\tag{8.10}$$

The reasoning in going from the difference equation (8.9) to the differential equation (8.10) follows exactly the steps explained in Sect. 8.2.1.

Before proceeding with how *I* and *R* develops in time, let us explain the formula *βΔtSI* . We have *S* susceptibles and *I* infected people. These can make up *SI* pairs. Now, suppose that during a time interval *T* we measure that *m* actual pairwise meetings do occur among *n* theoretically possible pairings of people from the S and I categories. The probability that people meet in pairs during a time *T* is (by the empirical frequency definition of probability) equal to *m/n*, i.e., the number of successes divided by the number of possible outcomes. From such statistics we normally derive quantities expressed per unit time, i.e., here we want the probability per unit time, *μ*, which is found from dividing by *T* : *μ* = *m/(nT )*.

Given the probability *μ*, the expected number of meetings per time interval of *SI* possible pairs of people is (from basic statistics) *μSI* . During a time interval *Δt*, there will be *μSIΔt* expected number of meetings between susceptibles and infected people such that the virus may spread. Only a fraction of the *μΔtSI* meetings are effective in the sense that the susceptible actually becomes infected. Counting that *m* people get infected in *n* such pairwise meetings (say 5 are infected from 1000 meetings), we can estimate the probability of being infected as *p* = *m/n*. The expected number of individuals in the S category that in a time interval *Δt* catch the virus and get infected is then *pμΔtSI* . Introducing a new constant *β* = *pμ* to save some writing, we arrive at the formula *βΔtSI* .

The value of *β* must be known in order to predict the future with the disease model. One possibility is to estimate *p* and *μ* from their meanings in the derivation above. Alternatively, we can observe an "experiment" where there are *S*0 susceptibles and *I*0 infected at some point in time. During a time interval *T* we count that *N* susceptibles have become infected. Using (8.9) as a rough approximation of how *S* has developed during time *T* (and now *T* is not necessarily small, but we use (8.9) anyway), we get

$$N = \beta^t S\_0 I\_0 \quad \Rightarrow \quad \beta = \frac{N}{T S\_0 I\_0} \,. \tag{8.11}$$

We need an additional equation to describe the evolution of *I (t)*. Such an equation is easy to establish by noting that the loss in the S category is a corresponding gain in the I category. More precisely,

$$I^{n+1} - I^n = \beta \Delta t S^n I^n \,. \tag{8.12}$$

However, there is also a loss in the I category because people recover from the disease. Suppose that we can measure that *m* out of *n* individuals recover in a time period *T* (say 10 of 40 sick people recover during a day: *m* = 10, *n* = 40, *T* = 24 h). Now, *γ* = *m/(nT )* is the probability that one individual recovers in a unit time interval. Then (on average) *γ ΔtI* infected will recover in a time interval *Δt*. This quantity represents a loss in the I category and a gain in the R category. We can therefore write the total change in the I category as

$$I^{n+1} - I^n = \beta \Delta t S^n I^n - \chi \Delta t I^n. \tag{8.13}$$

The change in the R category is simple: there is always an increase got from the I category:

$$R^{n+1} - R^n = \chi \, \Delta t I^n. \tag{8.14}$$

Since there is no loss in the R category (people are either recovered and immune, or dead), we are done with the modeling of this category. In fact, we do not strictly need Eq. (8.14) for *R*, but extensions of the model later will need an equation for *R*.

Dividing by *Δt* in (8.13) and (8.14) and letting *Δt* → 0, results in the corresponding differential equations

$$I' = \beta SI - \gamma I,\tag{8.15}$$

and

$$
\mathbb{R}' = \gamma I \,. \tag{8.16}
$$

To summarize, we have derived three difference equations and three differential equations, which we list here for easy reference. The difference equations are:

$$S^{n+1} = S^n - \beta \,\Delta t \, S^n I^n,\tag{8.17}$$

$$I^{n+1} = I^n + \beta \Delta t S^n I^n - \chi \Delta t I^n,\tag{8.18}$$

$$R^{n+1} = R^n + \chi \,\Delta t \, I^n \,. \tag{8.19}$$

Note that we have isolated the new unknown quantities *Sn*+1, *I n*+1, and *Rn*+1 on the left-hand side, such that these can readily be computed if *Sn*, *I n*, and *Rn* are known. To get such a procedure started, we need to know *S*0, *I* 0, *R*0. Obviously, we also need to have values for the parameters *β* and *γ* .

The three differential equations are:

$$S' = -\beta SI,\tag{8.20}$$

$$I' = \beta SI - \gamma I,\tag{8.21}$$

$$R' = \chi I \,. \tag{8.22}$$

This differential equation model (and also its discrete counterpart above) is known as an *SIR model*. The input data to the differential equation model consist of the parameter values for *β* and *γ* , as well as the initial conditions *S(*0*)* = *S*0, *I (*0*)* = *I*0, and *R(*0*)* = *R*0.

#### **8.3.2 A FE Method for the System of ODEs**

Let us apply the same principles as we did in Sect. 8.2.2 to discretize the differential equation system by the Forward Euler method. We already have a time mesh and time-discrete quantities *Sn*, *I n*, *Rn*, *n* = 0*,...,Nt* . The three differential equations are assumed to be valid at the mesh points. At the point *tn* we then have

$$S'(t\_n) = -\beta S(t\_n) I(t\_n),\tag{8.23}$$

$$I'(t\_n) = \beta S(t\_n) I(t\_n) - \gamma I(t\_n),\tag{8.24}$$

$$R'(t\_n) = \chi I(t\_n),\tag{8.25}$$

for *n* = 0*,* 1*,...,Nt* . This is an approximation since the differential equations are originally valid at all times *t* (usually in some finite interval [0*, T* ]). Using forward finite differences for the derivatives results in an additional approximation,

$$\frac{S^{n+1} - S^n}{\Delta t} = -\beta S^n I^n,\tag{8.26}$$

$$\frac{I^{n+1} - I^n}{\Delta t} = \beta S^n I^n - \gamma I^n,\tag{8.27}$$

$$\frac{\mathcal{R}^{n+1} - \mathcal{R}^n}{\Delta t} = \mathcal{Y}I^n. \tag{8.28}$$

As we can see, these equations are identical to the difference equations that naturally arise in the derivation of the model. However, other numerical methods than the Forward Euler scheme will result in slightly different difference equations.

#### **8.3.3 Programming the FE Scheme; the Special Case**

The computation of (8.26)–(8.28) can be readily made in a computer program SIR1.py:

```
import numpy as np
import matplotlib.pyplot as plt
# Time unit: 1 h
beta = 10./(40*8*24)
gamma = 3./(15*24)
dt = 0.1 # 6 min
D = 30 # Simulate for D days
N_t = int(D*24/dt) # Corresponding no of time steps
t = np.linspace(0, N_t*dt, N_t+1)
S = np.zeros(N_t+1)
I = np.zeros(N_t+1)
R = np.zeros(N_t+1)
# Initial condition
S[0] = 50
I[0] = 1
R[0] = 0
# Step equations forward in time
for n in range(N_t):
   S[n+1] = S[n] - dt*beta*S[n]*I[n]
   I[n+1] = I[n] + dt*beta*S[n]*I[n] - dt*gamma*I[n]
   R[n+1] = R[n] + dt*gamma*I[n]
fig = plt.figure()
l1, l2, l3 = plt.plot(t, S, t, I, t, R)
fig.legend((l1, l2, l3), ('S', 'I', 'R'), 'center right')
plt.xlabel('hours')
```

```
plt.savefig('tmp.pdf'); plt.savefig('tmp.png')
plt.show()
```
This program was written to investigate the spreading of flu at the mentioned boarding school, and the reasoning for the specific choices *β* and *γ* goes as follows. At some other school where the disease has already spread, it was observed that in the beginning of a day there were 40 susceptibles and 8 infected, while the numbers were 30 and 18, respectively, 24 h later. Using 1 h as time unit, we then have from (8.11) that *β* = 10*/(*40 · 8 · 24*)*. Among 15 infected, it was observed that 3 recovered during a day, giving *γ* = 3*/(*15 · 24*)*. Applying these parameters to a new case where there is one infected initially and 50 susceptibles, gives the graphs in Fig. 8.11. These graphs are just straight lines between the values at times *ti* = *iΔt* as computed by the program. We observe that *S* reduces as *I* and *R* grows. After about 30 days everyone has become ill and recovered again.

We can experiment with *β* and *γ* to see whether we get an outbreak of the disease or not. Imagine that a "wash your hands" campaign was successful and that the other school in this case experienced a reduction of *β* by a factor of 5. With this lower *β* the disease spreads very slowly so we simulate for 60 days. The curves appear in Fig. 8.12.

#### **8.3.4 Outbreak or Not**

Looking at the equation for *I* , it is clear that we must have *βSI* − *γI >* 0 for *I* to increase. When we start the simulation it means that

$$
\beta S(0)I(0) - \gamma I(0) > 0,
$$

**Fig. 8.11** Natural evolution of flu at a boarding school

**Fig. 8.12** Small outbreak of flu at a boarding school (*β* is much smaller than in Fig. 8.11)

or simpler

$$\frac{\beta S(0)}{\chi} > 1\tag{8.29}$$

to increase the number of infected people and accelerate the spreading of the disease. You can run the SIR1.py program with a smaller *β* such that (8.29) is violated and observe that there is no outbreak.

#### **The power of mathematical modeling**

The reader should notice our careful use of words in the previous paragraphs. We started out with modeling a very specific case, namely the spreading of flu among pupils and staff at a boarding school. With purpose we exchanged words like pupils and flu with more neutral and general words like *individuals* and *disease*, respectively. Phrased equivalently, we raised the abstraction level by moving from a specific case (flu at a boarding school) to a more general case (disease in a closed society). Very often, when developing mathematical models, we start with a specific example and see, through the modeling, that what is going on of essence in this example also will take place in many similar problem settings. We try to incorporate this generalization in the model so that the model has a much wider application area than what we aimed at in the beginning. This is the very power of mathematical modeling: by solving one specific case we have often developed more generic tools that can readily be applied to solve seemingly different problems. The next sections will give substance to this assertion.

#### **8.3.5 Abstract Problem and Notation**

When we had a specific differential equation with one unknown, we quickly turned to an abstract differential equation written in the generic form *u* = *f (u, t)*. We refer to such a problem as a *scalar ODE*. A specific equation corresponds to a specific choice of the formula *f (u, t)* involving *u* and (optionally) *t*.

It is advantageous to also write a system of differential equations in the same abstract notation,

$$u' = f(u, t),$$

but this time it is understood that *u* is a vector of functions and *f* is also vector. We say that *u* = *f (u, t)* is a *vector ODE* or *system of ODEs* in this case. For the SIR model we introduce the two 3-vectors, one for the unknowns,

*u* = *(S(t), I (t), R(t)),*

and one for the right-hand side functions,

$$f(\mu, t) = (-\beta SI, \beta SI - \gamma I, \gamma I) \dots$$

The equation *u* = *f (u, t)* means setting the two vectors equal, i.e., the components must be pairwise equal. Since *u* = *(S , I , R )*, we get that *u* = *f* implies

$$\begin{aligned} S' &= -\beta SI, \\ I' &= \beta SI - \gamma I, \\ R' &= \gamma I. \end{aligned}$$

The generalized short notation *u* = *f (u, t)* is very handy since we can derive numerical methods and implement software for this abstract system and in a particular application just identify the formulas in the *f* vector, implement these, and call functionality that solves the differential equation system.

#### **8.3.6 Programming the FE Scheme; the General Case**

In Python code, the Forward Euler step

$$
\mu^{n+1} = \mu^n + \Delta t f(\mu^n, t\_n),
$$

being a scalar or a vector equation, can be coded as

```
u[n+1] = u[n] + dt*f(u[n], t[n])
```
both in the scalar and vector case. In the vector case, u[n] is a one-dimensional numpy array of length *m* + 1 holding the mathematical quantity *un*, and the Python function f must return a numpy array of length *m* + 1. Then the expression u[n] + dt\*f(u[n], t[n]) is an array plus a scalar times an array.

For all this to work, the complete numerical solution must be represented by a two-dimensional array, created by u = zeros((N\_t+1, m+1)). The first index counts the time points and the second the components of the solution vector at one time point. That is, u[n,i] corresponds to the mathematical quantity *un i* . When we use only one index, as in u[n], this is the same as u[n,:] and picks out all the components in the solution at the time point with index n. Then the assignment u[n+1] = ... becomes correct because it is actually an in-place assignment u[n+1, :] = .... The nice feature of these facts is that the same piece of Python code works for both a scalar ODE and a system of ODEs!

The ode\_FE function for the vector ODE is placed in the file ode\_system\_FE.py and was written as follows:

```
import numpy as np
import matplotlib.pyplot as plt
def ode_FE(f, U_0, dt, T):
    N_t = int(round(T/dt))
    # Ensure that any list/tuple returned from f_ is wrapped as array
    f_ = lambda u, t: np.asarray(f(u, t))
    u = np.zeros((N_t+1, len(U_0)))
    t = np.linspace(0, N_t*dt, len(u))
    u[0] = U_0
    for n in range(N_t):
        u[n+1] = u[n] + dt*f_(u[n], t[n])
    return u, t
```
The line f\_ = lambda ... needs an explanation. For a user, who just needs to define the *f* in the ODE system, it is convenient to insert the various mathematical expressions on the right-hand sides in a list and return that list. Obviously, we could demand the user to convert the list to a numpy array, but it is so easy to do a general such conversion in the ode\_FE function as well. To make sure that the result from f is indeed an array that can be used for array computing in the formula u[n] + dt\*f(u[n], t[n]), we introduce a new function f\_ that calls the user's f and sends the results through the numpy function asarray, which ensures that its argument is converted to a numpy array (if it is not already an array).

Note also the extra parenthesis required when calling zeros with two indices.

Let us show how the previous SIR model can be solved using the new general ode\_FE that can solve *any* vector ODE. The user's f(u, t) function takes a vector u, with three components corresponding to *S*, *I* , and *R* as argument, along with the current time point t[n], and must return the values of the formulas of the right-hand sides in the vector ODE. An appropriate implementation is

```
def f(u, t):
    S, I, R = u
    return [-beta*S*I, beta*S*I - gamma*I, gamma*I]
```
Note that the S, I, and R values correspond to *Sn*, *I n*, and *Rn*. These values are then just inserted in the various formulas in the vector ODE. Here we collect the values in a list since the ode\_FE function will anyway wrap this list in an array. We could, of course, returned an array instead:

```
def f(u, t):
    S, I, R = u
    return array([-beta*S*I, beta*S*I - gamma*I, gamma*I])
```
The list version looks a bit nicer, so that is why we prefer a list and rather introduce f\_ = lambda u, t: asarray(f(u,t)) in the general ode\_FE function.

We can now show a function that runs the previous SIR example, while using the generic ode\_FE function:

```
def demo_SIR():
   """Test case using a SIR model."""
   def f(u, t):
       S, I, R = u
       return [-beta*S*I, beta*S*I - gamma*I, gamma*I]
   beta = 10./(40*8*24)
   gamma = 3./(15*24)
   dt = 0.1 # 6 min
   D = 30 # Simulate for D days
   N_t = int(D*24/dt) # Corresponding no of time steps
   T = dt*N_t # End time
   U_0 = [50, 1, 0]
   u, t = ode_FE(f, U_0, dt, T)
   S = u[:,0]
   I = u[:,1]
   R = u[:,2]
   fig = plt.figure()
   l1, l2, l3 = plt.plot(t, S, t, I, t, R)
   fig.legend((l1, l2, l3), ('S', 'I', 'R'), 'center right')
   plt.xlabel('hours')
   plt.show()
   # Consistency check:
   N = S[0] + I[0] + R[0]
   eps = 1E-12 # Tolerance for comparing real numbers
   for n in range(len(S)):
       SIR_sum = S[n] + I[n] + R[n]
       if abs(SIR_sum - N) > eps:
           print('*** consistency check failed: S+I+R={:g} != {:g}'\
                 .format(SIR_sum, N))
if __name__ == '__main__':
   demo_SIR()
```
Recall that the u returned from ode\_FE contains all components (*S*, *I* , *R*) in the solution vector at all time points. We therefore need to extract the *S*, *I* , and *R* values in separate arrays for further analysis and easy plotting.

Another key feature of this higher-quality code is the consistency check. By adding the three differential equations in the SIR model, we realize that *S* + *I* + *R* = 0, which means that *S* + *I* + *R* = const. We can check that this relation holds by comparing *Sn* + *I n* + *Rn* to the sum of the initial conditions. Exercise 8.6 suggests another method for controlling the quality of the numerical solution.

#### **8.3.7 Time-Restricted Immunity**

Let us now assume that immunity after the disease only lasts for some certain time period. This means that there is transport from the R state to the S state:

Modeling the loss of immunity is very similar to modeling recovery from the disease: the amount of people losing immunity is proportional to the amount of recovered patients and the length of the time interval *Δt*. We can therefore write the loss in the R category as −*νΔtR* in time *Δt*, where *ν*−1 is the typical time it takes to lose immunity. The loss in *R(t)* is a gain in *S(t)*. The "budgets" for the categories therefore become

$$S^{n+1} = S^n - \beta \Delta t S^n I^n + \nu \Delta t R^n,\tag{8.30}$$

$$I^{n+1} = I^n + \beta \,\Delta t S^n I^n - \chi \,\Delta t I^n,\tag{8.31}$$

$$R^{n+1} = R^n + \nu \Delta t I^n - \nu \Delta t R^n. \tag{8.32}$$

Dividing by *Δt* and letting *Δt* → 0 gives the differential equation system

$$S' = -\beta SI + \nu R,\tag{8.33}$$

$$I' = \beta SI - \gamma I,\tag{8.34}$$

$$R' = \chi I - \upsilon R\,. \tag{8.35}$$

This system can be solved by the same methods as we demonstrated for the original SIR model. Only one modification in the program is necessary: adding dt\*nu\*R[n] to the S[n+1] update and subtracting the same quantity in the R[n+1] update:

```
for n in range(N_t):
    S[n+1] = S[n] - dt*beta*S[n]*I[n] + dt*nu*R[n]
    I[n+1] = I[n] + dt*beta*S[n]*I[n] - dt*gamma*I[n]
    R[n+1] = R[n] + dt*gamma*I[n] - dt*nu*R[n]
```
The modified code is found in the file SIR2.py.

Setting *ν*−1 to 50 days, reducing *β* by a factor of 4 compared to the previous example (*β* = 0*.*00033), and simulating for 300 days gives an oscillatory behavior in the categories, as depicted in Fig. 8.13. It is easy now to play around and study how the parameters affect the spreading of the disease. For example, making the disease slightly more effective (increase *β* to 0.00043) and increasing the average time to loss of immunity to 90 days lead to other oscillations, see Fig. 8.14.

**Fig. 8.13** Including loss of immunity

**Fig. 8.14** Increasing *β* and reducing *ν* compared to Fig. 8.13

#### **8.3.8 Incorporating Vaccination**

We can extend the model to also include vaccination. To this end, it can be useful to track those who are vaccinated and those who are not. So, we introduce a fourth category, V, for those who have taken a successful vaccination. Furthermore, we assume that in a time interval *Δt*, a fraction *pΔt* of the S category is subject to a successful vaccination. This means that in the time *Δt*, *pΔtS* people leave from the S to the V category. Since the vaccinated ones cannot get the disease, there is no impact on the I or R categories. We can visualize the categories, and the movement between them, as

The new, extended differential equations with the *V* quantity become

$$S' = -\beta SI + \nu R - pS,\tag{8.36}$$

$$V' = p\mathcal{S},\tag{8.37}$$

$$I' = \beta SI - \gamma I,\tag{8.38}$$

$$
\boldsymbol{R}' = \boldsymbol{\gamma}\,\boldsymbol{I} - \boldsymbol{\nu}\boldsymbol{R}\,.\tag{8.39}
$$

We shall refer to this model as the SIRV model.

The new equation for *V* poses no difficulties when it comes to the numerical method. In a Forward Euler scheme we simply add an update

$$V^{n+1} = V'' + p\Delta t S'' \dots$$

The program needs to store *V (t)* in an additional array V, and the plotting command must be extended with more arguments to plot V versus t as well. The complete code is found in the file SIRV1.py.

Using *p* = 0*.*0005 and *p* = 0*.*0001 as values for the vaccine efficiency parameter, the effect of vaccination is seen in Figs. 8.15 and 8.16, respectively. (other parameters are as in Fig. 8.13).

#### **8.3.9 Discontinuous Coefficients: A Vaccination Campaign**

What about modeling a vaccination campaign? Imagine that 6 days after the outbreak of the disease, the local health station launches a vaccination campaign. They reach out to many people, say 10 times as efficiently as in the previous (constant vaccination) case. If the campaign lasts for 10 days we can write

$$p(t) = \begin{cases} 0.005, & 6 \cdot 24 \le t \le 15 \cdot 24, \\ 0, & \text{otherwise} \end{cases}$$

Note that we must multiply the *t* value by 24 because *t* is measured in hours, not days. In the differential equation system, *pS(t)* must be replaced by *p(t)S(t)*, and in this case we get a differential equation system with a term that is *discontinuous*. This

**Fig. 8.15** The effect of vaccination: *p* = 0*.*0005

**Fig. 8.16** The effect of vaccination: *p* = 0*.*0001

is usually quite a challenge in mathematics, but as long as we solve the equations numerically in a program, a discontinuous coefficient is easy to treat.

There are two ways to implement the discontinuous coefficient *p(t)*: through a function and through an array. The function approach is perhaps the easiest:

```
def p(t):
    return 0.005 if (6*24 <= t <= 15*24) else 0
```
Note the handy if-else construction in the return statement here. It is a one-line alternative to, for example,

```
if (6*24 <= t <= 15*24):
    return 0.005
```
**Fig. 8.17** The effect of a vaccination campaign

```
else:
    return 0
```
In the code for updating the arrays S and V, we then get a term p(t[n])\*S[n].

Alternatively, we can instead let *p(t)* be an array filled with correct values prior to the simulation. Then we need to allocate an array p of length N\_t+1 and find the indices corresponding to the time period between 6 and 15 days. These indices are found from the time point divided by *Δt*. That is,

```
p = zeros(N_t+1)
start_index = 6*24/dt
stop_index = 15*24/dt
p[start_index:stop_index] = 0.005
```
The *p(t)S(t)* term in the updating formulas for *S* and *V* simply becomes p[n]\*S[n]. The file SIRV2.py contains a program based on filling an array p.

The effect of a vaccination campaign is illustrated in Fig. 8.17. All the data are as in Fig. 8.15, except that *p* is ten times stronger for a period of 10 days and *p* = 0 elsewhere.

#### **8.4 Oscillating 1D Systems: A Second Order ODE**

Numerous engineering constructions and devices contain materials that act like springs. Such springs give rise to oscillations, and controlling oscillations is a key engineering task. We shall now learn to simulate oscillating systems.

As always, we start with the simplest meaningful mathematical model, which for oscillations is a second-order differential equation:

$$
\mu''(t) + \omega^2 \mu(t) = 0,\tag{8.40}
$$

where *ω* is a given physical parameter. Equation (8.40) models a one-dimensional system oscillating without damping (i.e., with negligible damping). Onedimensional here means that some motion takes place along one dimension only in some coordinate system. Along with (8.40) we need the two *initial conditions u(*0*)* and *u (*0*)*.

#### **8.4.1 Derivation of a Simple Model**

Many engineering systems undergo oscillations, and differential equations constitute the key tool to understand, predict, and control the oscillations. We start with the simplest possible model that captures the essential dynamics of an oscillating system. Some body with mass *m* is attached to a spring and moves along a line without friction, see Fig. 8.18 for a sketch (rolling wheels indicate "no friction"). When the spring is stretched (or compressed), the spring force pulls (or pushes) the body back and work "against" the motion. More precisely, let *x(t)* be the position of the body on the *x* axis, along which the body moves. The spring is not stretched when *x* = 0, so the force is zero, and *x* = 0 is hence the equilibrium position of the body. The spring force is −*kx*, where *k* is a constant to be measured. We assume that there are no other forces (e.g., no friction). Newton's second law of motion *F* = *ma* then has *F* = −*kx* and *a* = ¨*x*,

$$-k\mathbf{x} = m\ddot{\mathbf{x}},\tag{8.41}$$

which can be rewritten as

$$
\ddot{\mathbf{x}} + \omega^2 \mathbf{x} = \mathbf{0},
\tag{8.42}
$$

by introducing *ω* = √*k/m* (which is very common).

Equation (8.42) is a *second-order* differential equation, and therefore we need *two* initial conditions, one on the position *x(*0*)* and one on the velocity *x (*0*)*. Here

**Fig. 8.18** Sketch of a one-dimensional, oscillating dynamic system (without friction)

we choose the body to be at rest, but moved away from its equilibrium position:

$$\mathbf{x}(0) = X\_0, \quad \mathbf{x}'(0) = 0 \dots$$

The exact solution of (8.42) with these initial conditions is *x(t)* = *X*0 cos *ωt*. This can easily be verified by substituting into (8.42) and checking the initial conditions. The solution tells us that such a spring-mass system oscillates back and forth as described by a cosine curve.

The differential equation (8.42) appears in numerous other contexts. A classical example is a simple pendulum that oscillates back and forth. Physics books derive, from Newton's second law of motion, that

$$mL\theta'' + mg\sin\theta = 0,$$

where *m* is the mass of the body at the end of a pendulum with length *L*, *g* is the acceleration of gravity, and *θ* is the angle the pendulum makes with the vertical. Considering small angles *θ*, sin *θ* ≈ *θ*, and we get (8.42) with *x* = *θ*, *ω* = √*g/L*, *x(*0*)* = *Θ*, and *x (*0*)* = 0, if *Θ* is the initial angle and the pendulum is at rest at *t* = 0.

#### **8.4.2 Numerical Solution**

We have not looked at numerical methods for handling second-order derivatives, and such methods are an option, but we know how to solve first-order differential equations and even systems of first-order equations. With a little, yet very common, trick we can rewrite (8.42) as a first-order system of two differential equations. We introduce *u* = *x* and *v* = *x* = *u* as *two* new unknown functions. The two corresponding equations arise from the definition *v* = *u* and the original equation (8.42):

$$u' = v,\tag{8.43}$$

$$v' = -\omega^2 \mu \,. \tag{8.44}$$

(Notice that we can use *u* = *v* to remove the second-order derivative from Newton's second law.)

We can now apply the Forward Euler method to (8.43)–(8.44), exactly as we did in Sect. 8.3.2:

$$\frac{\mu^{n+1} - \mu^n}{\Delta t} = v^n,\tag{8.45}$$

$$\frac{v^{n+1} - v^n}{\Delta t} = -\omega^2 u^n,\tag{8.46}$$

resulting in the computational scheme

$$
\mu^{n+1} = \mu^n + \Delta t \, v^n,\tag{8.47}
$$

$$v^{n+1} = v^n - \Delta t \,\omega^2 u^n \,. \tag{8.48}$$

#### **8.4.3 Programming the FE Scheme; the Special Case**

A simple program for (8.47)–(8.48) follows the same ideas as in Sect. 8.3.3:

```
import numpy as np
import matplotlib.pyplot as plt
omega = 2
P = 2*np.pi/omega
dt = P/20
T = 3*P
N_t = int(round(T/dt))
t = np.linspace(0, N_t*dt, N_t+1)
u = np.zeros(N_t+1)
v = np.zeros(N_t+1)
# Initial condition
X_0 = 2
u[0] = X_0
v[0] = 0
# Step equations forward in time
for n in range(N_t):
    u[n+1] = u[n] + dt*v[n]
    v[n+1] = v[n] - dt*omega**2*u[n]
fig = plt.figure()
l1, l2 = plt.plot(t, u, 'b-', t, X_0*np.cos(omega*t), 'r--')
fig.legend((l1, l2), ('numerical', 'exact'), 'upper right')
plt.xlabel('t')
plt.savefig('tmp.pdf'); plt.savefig('tmp.png')
plt.show()
```
(See file osc\_FE.py.)

Since we already know the exact solution as *u(t)* = *X*0 cos *ωt*, we have reasoned as follows to find an appropriate simulation interval [0*, T* ] and also how many points we should choose. The solution has a period *P* = 2*π/ω*. (The period *P* is the time difference between two peaks of the *u(t)* ∼ cos*ωt* curve.) Simulating for three periods of the cosine function, *T* = 3*P*, and choosing *Δt* such that there are 20 intervals per period gives *Δt* = *P /*20 and a total of *Nt* = *T /Δt* intervals. The rest of the program is a straightforward coding of the Forward Euler scheme.

Figure 8.19 shows a comparison between the numerical solution and the exact solution of the differential equation. To our surprise, the numerical solution looks wrong. Is this discrepancy due to a programming error or a problem with the Forward Euler method?

**Fig. 8.19** Simulation of an oscillating system

First of all, even before trying to run the program, you should sit down and compute two steps in the time loop with a calculator so you have some intermediate results to compare with. Using *X*0 = 2, *dt* = 0*.*157079632679, and *ω* = 2, we get *u*1 = 2, *v*1 = −1*.*25663706, *u*2 = 1*.*80260791, and *v*2 = −2*.*51327412. Such calculations show that the program is seemingly correct. (Later, we can use such values to construct a unit test and a corresponding test function.)

The next step is to reduce the discretization parameter *Δt* and see if the results become more accurate. Figure 8.20 shows the numerical and exact solution for the cases *Δt* = *P /*40*,P/*160*,P/*2000. The results clearly become better, and the finest resolution gives graphs that cannot be visually distinguished. Nevertheless, the finest resolution involves 6000 computational intervals in total, which is considered quite much. This is no problem on a modern laptop, however, as the computations take just a fraction of a second.

Although 2000 intervals per oscillation period seem sufficient for an accurate numerical solution, the lower right graph in Fig. 8.20 shows that if we increase the simulation time, here to 20 periods, there is a little growth of the amplitude, which becomes significant over time. The conclusion is that the Forward Euler method has a fundamental problem with its growing amplitudes, and that a very small *Δt* is required to achieve satisfactory results. The longer the simulation is, the smaller *Δt* has to be. It is certainly time to look for more effective numerical methods!

#### **8.4.4 A Magic Fix of the Numerical Method**

In the Forward Euler scheme,

$$\begin{aligned} u^{n+1} &= u^n + \Delta t \, v^n, \\ v^{n+1} &= v^n - \Delta t \, w^2 u^n, \end{aligned}$$

**Fig. 8.20** Simulation of an oscillating system with different time steps. Upper left: 40 steps per oscillation period. Upper right: 160 steps per period. Lower left: 2000 steps per period. Lower right: 2000 steps per period, but longer simulation

we can replace *un* in the last equation by the recently computed value *un*+1 from the first equation:

$$
\mu^{n+1} = \mu^n + \Delta t \, v^n,\tag{8.49}
$$

$$v^{n+1} = v^n - \Delta t \, \alpha^2 u^{n+1} \,. \tag{8.50}$$

Before justifying this fix more mathematically, let us try it on the previous example. The results appear in Fig. 8.21. We see that the amplitude *does not grow*, but the phase is not entirely correct. After 40 periods (Fig. 8.21 right) we see a significant difference between the numerical and the exact solution. Decreasing *Δt* decreases the error. For example, with 2000 intervals per period, we only see a small phase error even after 50,000 periods (!). We can safely conclude that the fix results in an excellent numerical method!

Let us interpret the adjusted scheme mathematically. First we order (8.49)–(8.50) such that the difference approximations to derivatives become transparent:

$$\frac{\mu^{n+1} - \mu^n}{\Delta t} = v^n,\tag{8.51}$$

$$\frac{v^{n+1} - v^n}{\Delta t} = -\omega^2 u^{n+1}.\tag{8.52}$$

**Fig. 8.21** Adjusted method: first three periods (left) and period 36–40 (right)

**Fig. 8.22** Illustration of a backward difference approximation to the derivative

We interpret (8.51) as the differential equation sampled at mesh point *tn*, because we have *vn* on the right-hand side. The left-hand side is then a *forward difference* or Forward Euler approximation to the derivative *u* , see Fig. 8.4. On the other hand, we interpret (8.52) as the differential equation sampled at mesh point *tn*+1, since we have *un*+1 on the right-hand side. In this case, the difference approximation on the left-hand side is a *backward difference*,

$$v'(t\_{n+1}) \approx \frac{v^{n+1} - v^n}{\Delta t} \quad \text{or} \quad v'(t\_n) \approx \frac{v^n - v^{n-1}}{\Delta t}.$$

Figure 8.22 illustrates the backward difference. The error in the backward difference is proportional to *Δt*, the same as for the forward difference (but the proportionality constant in the error term has different sign). The resulting discretization method, seen in (8.52), is often referred to as a Backward Euler scheme (a first-order scheme, just like Forward Euler).

To summarize, using a forward difference for the first equation and a backward difference for the second equation results in a much better method than just using forward differences in both equations.

The standard way of expressing this scheme in physics is to change the order of the equations,

$$v' = -\omega^2 u,\tag{8.53}$$

$$u' = v,\tag{8.54}$$

and apply a forward difference to (8.53) and a backward difference to (8.54):

$$v^{n+1} = v^n - \Delta t \,\omega^2 u^n,\tag{8.55}$$

$$
\mu^{n+1} = \mu^n + \Delta t \,\, v^{n+1} \,. \tag{8.56}
$$

That is, first the velocity *v* is updated and then the position *u*, using the most recently computed velocity. There is no difference between (8.55)–(8.56) and (8.49)–(8.50) with respect to accuracy, so how you order the original differential equations does not matter. The scheme (8.55)–(8.56) goes by the name Semi-implicit Euler, 5 or Euler-Cromer (a first-order method). The implementation of (8.55)–(8.56) is found in the file osc\_EC.py. The core of the code goes like

```
u = zeros(N_t+1)
v = zeros(N_t+1)
# Initial condition
u[0] = 2
v[0] = 0
# Step equations forward in time
for n in range(N_t):
    v[n+1] = v[n] - dt*omega**2*u[n]
    u[n+1] = u[n] + dt*v[n+1]
```
#### **Explicit and implicit methods**