# **Simula SpringerBriefs on Computing 15**

**Joakim Sundnes**

# **Solving Ordinary Differential Equations in Python**

# **Simula SpringerBriefs on Computing**

### Volume 15

### **Editor-in-Chief**

Joakim Sundnes, Simula Research Laboratory, Oslo, Norway

### **Series Editors**

Shaukat Ali, Simula Research Laboratory, Oslo, Norway Evrim Acar Ataman, Simula Metropolitan Centre for Digital Engineering, Oslo, Norway Are Magnus Bruaset, Simula Research Laboratory, Oslo, Norway Xing Cai, Simula Research Laboratory & University, Oslo, Norway Kimberly Claffy, San Diego Supercomputer Center, CAIDA, University of California, San Diego, San Diego, CA, USA Andrew Edwards, Simula Research Laboratory, Oslo, Norway Arnaud Gotlieb, Simula Research Laboratory, Oslo, Norway Magne Jørgensen, Software Engineering, Simula Research Laboratory, Oslo, Norway Olav Lysne, Simula Research Laboratory, Oslo, Norway Kent-Andre Mardal, University of Oslo & Simula Research Lab, Oslo, Norway Kimberly McCabe, Simula Research Laboratory, Oslo, Norway Andrew McCulloch, Bioengineering 0412, University of California, San Diego, La Jolla, CA, USA Leon Moonen, Simula Research Laboratory, Oslo, Norway Michael Riegler, Simula Metropolitan Centre for Digital Engineering & UiT The Arctic University of Norway, Oslo, Norway Marie Rognes, Simula Research Laboratory & University, Oslo, Norway Fabian Theis, Institute of Computational Biology, Helmholtz Zentrum München, Neuherberg, Germany Aslak Tveito, Simula Research Laboratory, Oslo, Norway Karen Willcox, Oden Institute for Computational Engineering & Science, The University of Texas at Austin, Austin, MA, USA Tao Yue, Nanjing University of Aeronautics and Astronautics & Simula Research Laboratory, Oslo, Norway Andreas Zeller, Saarland University, Saarbrücken, Germany

Yan Zhang, University of Oslo & Simula Research Laboratory, Oslo, Norway

In 2016, Springer and Simula launched the book series *Simula SpringerBriefs on Computing*, which aims to provide introductions to selected research topics in computing. The series provides compact introductions for students and researchers entering a new field, brief disciplinary overviews of the state-of-the-art of select fields, and raises essential critical questions and open challenges in the field of computing. Published by SpringerOpen, all *Simula SpringerBriefs on Computing* are open access, allowing for faster sharing and wider dissemination of knowledge.

Simula Research Laboratory is a leading Norwegian research organization which specializes in computing. Going forward, the book series will provide introductory volumes on the main topics within Simula's expertise, including communications technology, software engineering and scientific computing.

By publishing the *Simula SpringerBriefs on Computing,* Simula Research Laboratory acts on its mandate of emphasizing research education. Books in this series are published by invitation from one of the series editors. Authors interested in publishing in the series are encouraged to contact any member of the editorial board.

Joakim Sundnes

# Solving Ordinary Differential Equations in Python

Joakim Sundnes Simula Research Laboratory AS Oslo, Norway

ISSN 2512-1677 ISSN 2512-1685 (electronic) Simula SpringerBriefs on Computing ISBN 978-3-031-46767-7 ISBN 978-3-031-46768-4 (eBook) https://doi.org/10.1007/978-3-031-46768-4

Mathematics Subject Classification: 34-01, 34-04, 65-L05, 65-L04, 65-L06, 39-08

© The Author(s) 2024. This book is an open access publication.

**Open Access** This book 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 license and indicate if changes were made.

The images or other third party material in this book are included in the book's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the book's Creative Commons license 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.

The use of general descriptive names, registered names, trademarks, service marks, etc. in this publication does not imply, even in the absence of a specific statement, that such names are exempt from the relevant protective laws and regulations and therefore free for general use.

The publisher, the authors, and the editors are safe to assume that the advice and information in this book are believed to be true and accurate at the date of publication. Neither the publisher nor the authors or the editors give a warranty, expressed or implied, with respect to the material contained herein or for any errors or omissions that may have been made. The publisher remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

This Springer imprint is published by the registered company Springer Nature Switzerland AG The registered company address is: Gewerbestrasse 11, 6330 Cham, Switzerland

Paper in this product is recyclable.

# **Series Foreword**

### Dear reader,

Scientific research is increasingly interdisciplinary, and both students and experienced researchers often face the need to learn the foundations, tools, and methods of a new research field. This process can be quite demanding, and typically involves extensive literature searches and reading dozens of scientific papers in which the notation and style of presentation varies considerably. Since the establishment of this series in 2016 by founding editor-in-chief Aslak Tveito, the briefs in this series have aimed to ease the process by introducing and explaining important concepts and theories in a relatively narrow field, and to outline open research challenges and pose critical questions on the fundamentals of that field. The goal is to provide the necessary understanding and background knowledge and to motivate further studies of the relevant scientific literature. A typical brief in this series should be around 100 pages and should be well suited as material for a research seminar in a well-defined and limited area of computing.

We publish all items in this series under the SpringerOpen framework, as this allows authors to use the series to publish an initial version of their manuscript that could subsequently evolve into a full-scale book on a broader theme. Since the briefs are freely available online, the authors do not receive any direct income from the sales; however, remuneration is provided for every completed manuscript. Briefs are written on the basis of an invitation from a member of the editorial board. Suggestions for possible topics are most welcome and can be sent to<EMAIL_ADDRESS>
March 2023 *Dr. Joakim Sundnes* Editor-in-Chief Simula Research Laboratory<EMAIL_ADDRESS>
> *Dr. Martin Peters* Executive Editor Mathematics Springer Heidelberg, Germany<EMAIL_ADDRESS>
Aslak Tveito, Simula Research Laboratory, Oslo, Norway

# **Preface**

This book was based on a set of lecture notes originally written for the book *A Primer on Scientific Programming with Python* by Hans Petter Langtangen [14], mainly covering topics from Appendices A, C, and E. To provide a more comprehensive overview of state-of-the art solvers for ordinary differential equations (ODEs), the notes have been extended with additional material on implicit solvers and automatic time-stepping methods. The main purpose of the notes is to serve as a concise and gentle introduction to solving differential equations in Python, specifically for the course *Introduction to programming for scientific applications* (IN1900, 10 ETCS credits) at the University of Oslo. These notes will be most useful for readers with a basic knowledge of Python and NumPy, see for instance [16], and it is also useful to have a fundamental understanding of ODEs.

One may question the usefulness of learning how to write your own ODE solvers in Python when there are already multiple solvers available, such as those in the SciPy library. However, no single ODE solver is universally optimal and efficient for all ODE problems, and the choice of solver should always be based on the specific characteristics of the problem at hand. To make the right choice, it is extremely beneficial to understand the strengths and weaknesses of different solvers, and the best way to gain this knowledge is by programming your own collection of ODE solvers. Different ODE solvers are conveniently grouped into families and hierarchies, offering an excellent example of how object-oriented programming (OOP) can maximize code reuse and minimize duplication.

The book's presentation style is compact and pragmatic, incorporating numerous code examples to illustrate how various ODE solvers can be implemented and applied in practice. The complete source code for all examples, as well as Jupyter notebooks for each chapter, are provided in the accompanying online resources. The programs and code examples are written in a simple and compact Python style, avoiding the use of advanced tools and features. Experienced Python programmers may find more elegant and modern solutions to many of the examples, utilizing abstract base classes, type hints, data classes, and other advanced features. However, the book's main goal is to introduce the fundamentals of ODE solvers and OOP as part of an introductory programming course, and we believe this purpose is best served by focusing on the basics.

Readers familiar with scientific computing or numerical software may also miss a discussion of computational performance. While performance is certainly relevant when solving ODEs, optimizing the performance of a Pythonbased solver easily becomes quite technical, and requires features like just-intime compilers (e.g., Numba) or mixed-language programming. The solvers in this book use fairly basic features of Python and NumPy, sacrificing some performance in favor of enhancing understanding of solver properties and implementation.1

The book is organized as follows: Chapter 1 introduces the forward Euler method, serving as a foundation for understanding the principles underlying all the methods covered later. It introduces the notation and mathematical formulation used throughout the book for scalar ODEs and systems of ODEs, and is essential reading for those with limited prior experience with ODEs and ODE solvers. Additionally, it briefly explains how to use the ODE solvers from the SciPy library. Readers already familiar with the fundamentals of the forward Euler method and its implementation may consider proceeding straight to Chapter 2, which presents explicit Runge-Kutta methods. The chapter introduces the fundamental ideas of these methods, but the main focus is on the implementation and how a collection of ODE solvers is conveniently implemented as a class hierarchy. Chapter 3 introduces *stiff* ODEs, presents techniques for performing simple stability analysis of Runge-Kutta methods, and introduces implicit Runge-Kutta methods. The majority of the chapter is dedicated to the programming of these solvers, which exhibit better stability properties than explicit methods and are therefore more suitable for solving stiff ODEs. Chapter 4 concludes the presentation of ODE solvers by introducing methods for adaptive time step control, which is an essential component of all modern ODE software. Chapter 5 takes a different approach from the preceding chapters, as it focuses on a specific class of ODE models rather than a set of solvers. While the simpler ODE problems discussed in earlier chapters serve the purpose of introducing and testing the solvers, it is valuable to explore more complex ODE models in order to appreciate both the potential and the challenges of modeling with ODEs. As an example, the chapter examines the famous Kermack-McKendrick SIR (Susceptible-Infected-Recovered) model from epidemiology. These classic models were developed in the early 1900s (see [12]) and remain fundamental for predicting and understanding the spread of infectious diseases. We describe the derivation of the models from a set of fundamental assumptions, and discuss the implications and limitations resulting from these assumptions. The main focus of the chapter is then on modifying and extending the models to capture new phenomena, and

1Complete source code for all the solvers and examples in the book can be found here: https://sundnes.github.io/solving\_odes\_in\_python/

demonstrating how these changes can be implemented and explored using the solvers developed in preceding chapters.

Finally, while the main focus of the text is on differential equations, Appendix A is dedicated to the related topic of *difference equations*. Difference equations have important applications on their own and may serve as a stepping stone towards understanding and solving ODEs, since numerical methods for ODEs essentially involve transforming differential equations into difference equations. The standard formulation of difference equations found in mathematical textbooks is already well-suited for computer implementation, using for-loops and arrays. Some students find difference equations easier to grasp than differential equations, making Appendix A a useful resource to begin with. However, others may prefer to dive straight into ODEs and explore Appendix A at a later stage.

*July 2023 Joakim Sundnes*

# **Contents**

# **Chapter 1 Programming a Simple ODE Solver**

Ordinary differential equations (ODEs) are widely used in science and engineering, particularly when it comes to modeling dynamic processes. Although analytical methods can be employed to solve simple ODEs, nonlinear ODEs typically require numerical methods for solutions. In this chapter we demonstrate how to program general numerical solvers capable of handling any ODE. Initially we will focus on scalar ODEs, which consist of a single equation and a single unknown. Subsequently, in Section 1.3, we will extend these concepts to systems of coupled ODEs. Acquiring a solid grasp of the concepts presented in this chapter will not only help you with programming your own ODE solvers but also in using a diverse range of readily available, general-purpose ODE solvers in Python or other programming languages.

### **1.1 Creating a General-Purpose ODE Solver**

When solving ODEs analytically, one typically considers a specific ODE or a class of ODEs and attempts to derive a formula for the solution. However, in this chapter, our goal is to implement numerical solvers that can be applied to any ODE, without being limited to a single example or a specific class of equations. To achieve this, we need a general abstract notation for an arbitrary ODE. We will write the ODEs on the following form:

$$u'(t) = f(t, u(t)),\tag{1.1}$$

which means that the ODE is fully specified by the definition of the righthand side function *f*(*t,u*). Examples of this function may be:

$$\begin{aligned} f(t, u) &= \alpha u, \quad \text{exponential growth} \\ f(t, u) &= \alpha u \left( 1 - \frac{u}{R} \right), \quad \text{logistic growth} \\ f(t, u) &= -b |u| u + g, \quad \text{fallsing body in a fluid} \end{aligned}$$

Notice that, for the sake of generality, we write all the right-hand sides of the ODEs as functions of both *t* and *u*, even though the mathematical formulations only involve *u*. This general formulation is not strictly necessary in the mathematical equations, but it proves to be highly convenient when we start programming and want to use the same solver for a diverse range of ODE models. We will delve into this topic in greater detail later. Now, our objective is to write functions and classes that accept the function *f* as input and solve the corresponding ODE to generate the output *u*.

To ensure a unique solution for (1.1), it is necessary to specify the *initial condition* for *u*. This initial condition corresponds to the value of the solution at a specific time *t* = *t*0. The resulting mathematical problem can be expressed as

$$\begin{aligned} u'(t) &= f(t, u(t)), \\ u(t\_0) &= u\_0, \end{aligned}$$

and is commonly referred to as an *initial value problem*, or simply an IVP. Every ODE problem discussed in this book is an initial value problem. To illustrate, let us consider the very simple ODE

$$u' = u.$$

This general solution of this equation is given by *u*(*t*) = *Cet* for any constant *C*, implying that there exist an infinite number of solutions. However, by specifying an initial condition *u*(*t*0) = *u*0, we get *C* = *u*0 and the unique solution *u*(*t*) = *u*0*e t* . When solving the equation numerically, it is necessary to define the initial condition *u*0 in order to start our method and compute a solution at all.

**A Simple and General Solver: the Forward Euler Method.** A numerical method for (1.1) can be derived by using a finite difference approximation for the derivative in the equation *u* 0 = *f*(*t,u*). To introduce this idea, let us assume that we have already computed *u* at discrete time points *t*0*, t*1*,..., tn*. At time *tn* we have the ODE

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

and we can now approximate *u* 0 (*tn*) with a forward finite difference:

$$u'(t\_n) \approx \frac{u(t\_{n+1}) - u(t\_n)}{\Delta t}.$$

By inserting this approximation into the ODE at *t* = *tn*, we obtain the following equation

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

and we can rearrange the terms to obtain an explicit formula for *u*(*tn*+1):

$$u(t\_{n+1}) = u(t\_n) + \Delta t f(t\_n, u(t\_n)).$$

This method, known as the *Forward Euler (FE) method* or the *Explicit Euler method*, is the simplest numerical method for solving an ODE. The terms *forward* and *explicit* refer to the fact that we have an explicit update formula for *u*(*tn*+1) that only involves known quantities at time *tn*. In contrast, an *implicit* ODE solver would have an update formula that includes terms like *f*(*tn*+1*,u*(*tn*+1)), requiring the solution of a generally nonlinear equation to determine the unknown *u*(*tn*+1). We will explore other explicit ODE solvers in Chapter 2 and implicit solvers in Chapter 3.

To simplify the formula, we introduce the notation *un* = *u*(*tn*), i.e., we let *un* represent the numerical approximation to the exact solution *u*(*t*) at *t* = *tn*. With this notation, the update formula reads

$$u\_{n+1} = u\_n + \Delta t f(t\_n, u\_n),\tag{1.2}$$

which, if we know the *u*0 at time *t*0, can be applied repeatedly to *u*1, *u*2, *u*3 and so forth. If we again consider the very simple ODE given by *u* 0 = *u* (i.e., *f*(*t,u*) = *u*), we have

$$\begin{aligned} u\_1 &= u\_0 + \Delta t u\_0, \\ u\_2 &= u\_1 + \Delta t u\_1, \\ u\_3 &= u\_2 + \dots, \end{aligned}$$

and the general update formula

$$u\_{n+1} = u\_n + \Delta t u\_n = (1 + \Delta t) u\_n \dots$$

In a Python program, the repeated application of the same formula can be conveniently implemented using a for-loop, and the solution can be stored in a list or a NumPy array. If you are unfamiliar with NumPy arrays and their usage, we recommend referring to [16], which provides an introduction to NumPy arrays and tools that will be extensively used through this book. To solve the ODE numerically, given a final time *T* and the number of time steps *N*, we can follow these steps:

1For *N* time steps, the length of the arrays needs to be *N* + 1 since we need to store both end points, i.e., *t*0*, t*1*,..., tn* and *u*0*,u*1*,...,un*.

	- t[n + 1] = t[n] + dt • u[n + 1] = (1 + dt) \* u[n]

A complete Python implementation of this algorithm may look like

```
import numpy as np
import matplotlib.pyplot as plt
N = 20
T = 4
dt = T/N
u0 = 1
t = np.zeros(N + 1)
u = np.zeros(N + 1)
u[0] = u0
for n in range(N):
    t[n + 1] = t[n] + dt
    u[n + 1] = (1 + dt) * u[n]
plt.plot(t, u)
plt.show()
```
Notice that there is no need to set t[0]= 0 when t is created in this way, but it is important to update u[0]. Forgetting to do so is a common error in ODE programming, so it is worth taking note of the line u[0] = u0. The solution is shown in Figure 1.1 for two different choices of the time step *∆t*. As observed, the approximate solution improves as *∆t* is reduced, although both the solutions deviate from the exact solution. However, reducing the time step further would easily yield a solution that is indistinguishable from the exact solution.

The for-loop in the aforementioned example could also be implemented differently, for instance

**for** n **in range**(1, N+1): t[n] = t[n - 1] + dt u[n] = (1 + dt) \* u[n - 1]

Here, the index n runs from 1 to N, and all the indices inside the loop have been decreased by one to achieve the same outcome. In this simple case, it is easy to verify that both loop formulations give the same result. However, mixing up the two formulations can easily lead to errors, such as a loop that exceeds the array bounds (resulting in an IndexError) or a loop where the last elements of t and u are not computed. Although these errors may appear trivial, they are common pitfalls when working with for-loops and it is good practice to always examine the loop formulation to ensure consistent use of indices and bounds.

**Fig. 1.1** Solution of *u* 0 = *u,u*(0) = 1 with *∆t* = 0*.*4 (*N* = 10) and *∆t* = 0*.*2 (*N* = 20).

**Extending the Solver to the General ODE.** As mentioned earlier, the goal of this chapter is to develop general-purpose ODE solvers capable of solving any ODE expressed in the form *u* 0 = *f*(*t,u*). Achieving this requires only a slight modification of the algorithm presented above:

	- t[n + 1] = t[n] + dt
	- u[n + 1] = u[n] + dt \* f(t[n], u[n])

The modified version of the algorithm only requires a small change in the formula for computing u[n+1] from u[n]. In the previous case we had *f*(*t,u*) = *u*, and to create a general-purpose ODE solver we simply replace u[n] with the more general f(t[n],u[n]). The following Python function implements this generic version of the FE method:2

```
import numpy as np
def forward_euler(f, u0, T, N):
    """Solve u'=f(t, u), u(0)=u0, with n steps until t=T."""
    t = np.zeros(N + 1)
    u = np.zeros(N + 1) # u[n] is the solution at time t[n]
    u[0] = u0
    dt = T / N
    for n in range(N):
        t[n + 1] = t[n] + dt
        u[n + 1] = u[n] + dt * f(t[n], u[n])
    return t, u
```
This simple function can solve any ODE expressed in the form (1.1). The right-hand side function *f*(*t,u*) must be implemented as a Python function, which is then passed as an argument to forward\_euler, along with the initial condition u0, the stop time T and the number of time steps N. Inside the function, the time step dt is calculated using T and N.

To illustrate the usage of the forward\_euler function, let us apply it to solve the same problem as before: *u* 0 = *u*, with the initial condition *u*(0) = 1, for *t* ∈ [0*,*4]. The following code uses the forward\_euler function to solve this problem:

```
def f(t, u):
    return u
u0 = 1
T = 4
N = 30
t, u = forward_euler(f, u0, T, N)
```
The forward\_euler function returns two arrays, t and u, which can be further processed or plotted as desired. An important aspect to note in this code is the definition of the right-hand side function f. As mentioned earlier, this function should always be written with two arguments, t and u, although in this case only u is used inside the function. The inclusion of both arguments is necessary because we want our solver to be applicable for all ODEs in the form *u* 0 = *f*(*t,u*). Therefore, inside the forward\_euler function, the f function is called as f(t[n], u[n]). If the right-hand side function were defined as a function of u only, i.e., using def f(u):, an error would occur when

2The source code for this function, as well as all subsequent solvers and examples, can be found here: https://sundnes.github.io/solving\_odes\_in\_python/

calling the function inside forward\_euler. To avoid this issue, we simply write def f(t,u): even if t is not used inside the function.3

For being only 15 lines of code, the capabilities of the forward\_euler function are quite remarkable. Using this function, we can solve any kind of linear or nonlinear ODE, most of which would be impossible to solve using analytical techniques. To use this function, follow these general steps:

It is important to note that the FE method is the simplest of all ODE solvers, and many will argue that it is not very good. This is partly true, there exist other methods that offer greater accuracy and stability when applied to challenging ODEs. As we will see later, numerical solutions obtained using the FE method may not only be inaccurate, as depicted in Figure 1.1, but can also diverge or exhibit instability. In Chapter 3, we will explore solvers that address such issues by providing improved stability. However, the FE method is quite suitable for solving a wide range of interesting ODEs. If we are not happy with the accuracy we can simply reduce the time step, and in most cases this will give the accuracy we need with a negligible increase in computational time.

### **1.2 The ODE Solver Implemented as a Class**

We can increase the flexibility of the forward\_euler solver function by implementing it as a class. There are many ways to implement such a class, but one possible usage can be as follows:

```
method = ForwardEuler_v0(f)
method.set_initial_condition(u0)
t, u = method.solve(t_span=(0, 10), N=100)
plot(t, u)
```
The benefits of using a class instead of a function may not be obvious at this point, but it will become clear when we introduce different ODE solvers later. For now, let us just look at how such a solver class can be implemented to support the specified use case:

3This way of defining the right-hand side is a standard used by most available ODE solver libraries, both in Python and other languages. The right-hand side function always takes two arguments t and u, but, annoyingly, the order of the two arguments varies between different solver libraries. Some expect the t argument first, while others expect u first.

In addition to the mentioned methods, it can be convenient to implement a separate method, for instance called advance, for advancing the solution one time step. This approach simplifies the implementation of new numerical methods, as we often only need to modify the advance method. A first version of the solver class can be implemented as follows:

```
import numpy as np
class ForwardEuler_v0:
    def __init__(self, f):
        self.f = f
    def set_initial_condition(self, u0):
        self.u0 = u0
    def solve(self, t_span, N):
        """Compute solution for t_span[0] <= t <= t_span[1],
        using N steps."""
        t0, T = t_span
        self.dt = T / N
        self.t = np.zeros(N + 1) # N steps ~ N+1 time points
        self.u = np.zeros(N + 1)
        msg = "Please set initial condition before calling solve"
        assert hasattr(self, "u0"), msg
        self.t[0] = t0
        self.u[0] = self.u0
        for n in range(N):
            self.n = n
            self.t[n + 1] = self.t[n] + self.dt
            self.u[n + 1] = self.advance()
        return self.t, self.u
    def advance(self):
        """Advance the solution one time step."""
        # Create local variables to get rid of "self." in
        # the numerical formula
```

```
u, dt, f, n, t = self.u, self.dt, self.f, self.n, self.t
return u[n] + dt * f(t[n], u[n])
```
This class performs the same tasks as the forward\_euler function mentioned earlier, with the main advantage of the class implementation being the enhanced flexibility provided by the advance method. As we shall see later, implementing a different numerical method typically only requires implementing a new version of this method, leaving the rest of the code unchanged. An additional improvement in the class implementation is the inclusion of an assert statement within the solve method. This statement verifies that the user has called set\_initial\_condition before calling solve. Forgetting to do so is a common mistake, and the assert statement ensures that a useful error message is raised rather than a less informative AttributeError.

We can also use a class to represent the right-hand side function *f*(*t,u*), which is particularly convenient for functions with parameters. Consider, for instance, the model for logistic growth:

$$u'(t) = \alpha u(t) \left(1 - \frac{u(t)}{R}\right), \quad u(0) = u\_0, \quad t \in [0, 40],$$

which is typically used to model self-limiting growth of a biological population, i.e., growth that is constrained by limited resources. Initially, the growth follows an approximately exponential pattern with growth rate *α*. As the population size approaches the *carrying capacity R*, the population curve flattens out, see Figure 1.2 for an example solution. The right-hand side function includes two parameters *α* and *R*, but if we want to solve this model using the FE function or class, the function must be implemented as a function of *t* and *u* only. There are several ways to achieve this in Python, but a convenient approach is to implement the function as a class with a call method.4 We can then define the parameters as attributes in the constructor and use them within the \_\_call\_\_ method:

```
class Logistic:
    def __init__(self, alpha, R):
        self.alpha, self.R = alpha, float(R)
    def __call__(self, t, u):
        return self.alpha * u * (1 - u / self.R)
```
The main program for solving the logistic growth problem may now look like:

```
problem = Logistic(alpha=0.2, R=1.0)
solver = ForwardEuler_v0(problem)
u0 = 0.1
solver.set_initial_condition(u0)
t, u = solver.solve(t_span=(0, 40), N=400)
```
4Recall that if we equip a class with a special method named \_\_call\_\_, instances of the class will be callable and will behave like regular Python functions. See, for instance, Chapter 8 of [16] for a brief introduction to \_\_call\_\_ and other special methods.

**Fig. 1.2** Solution of the logistic growth model.

### **1.3 Systems of ODEs**

Up until now, our focus has been on solving ODEs with a single solution component, commonly know as scalar ODEs. However, many interesting processes can be described by systems of ODEs, which consist of multiple ODEs where the right-hand side of one equation depends on the solution of the others. Such equation systems are also referred to as vector ODEs. One simple example is

$$\begin{aligned} u' &= v, & u(0) &= 1\\ v' &= -u, & v(0) &= 0. \end{aligned}$$

The solution of this particular system is *u* = cos*t,v* = sin*t*, which can be easily verified by inserting the solution into the equations and the initial conditions. For more general cases, it is usually even more difficult to find analytical solutions of ODE systems than of scalar ODEs, and numerical methods are usually necessary. In this section we will extend the solvers introduced in sections 1.1-1.2 to handle systems of ODEs. We shall see that such an extension requires relatively small modifications of the code.

Our goal is to develop general software capable of solving any vector ODE or scalar ODE. To achieve this, it is helpful to introduce some general mathematical notation. We have *m* unknowns

$$u^{(0)}(t), u^{(1)}(t), \dots, u^{(m-1)}(t)$$

in a system of *m* ODEs:

$$\begin{aligned} \frac{d}{dt}u^{(0)} &= f^{(0)}(t, u^{(0)}, u^{(1)}, \dots, u^{(m-1)}),\\ \frac{d}{dt}u^{(1)} &= f^{(1)}(t, u^{(0)}, u^{(1)}, \dots, u^{(m-1)}),\\ \vdots &= \vdots\\ \frac{d}{dt}u^{(m-1)} &= f^{(m-1)}(t, u^{(0)}, u^{(1)}, \dots, u^{(m-1)}).\end{aligned}$$

To simplify the notation (and later the implementation), we can collect both the solutions *u* (*i*) (*t*) and right-hand side functions *f* (*i*) into vectors;

$$u = (u^{(0)}, u^{(1)}, \dots, u^{(m-1)}),$$

and

$$f = (f^{(0)}, f^{(1)}, \dots, f^{(m-1)}).$$

Note that *f* is now a vector-valued function. It takes *m*+ 1 input arguments (*t* and the *m* components of *u*) and returns a vector of *m* values. Using this notation, the ODE system can be written

$$u' = f(t, u), \quad u(t\_0) = u\_0,$$

where *u* and *f* are now vectors and *u*0 is a vector of initial conditions. We observe that the notation used for scalar ODEs remains the same, and whether we are solving a scalar or system of ODEs is determined by how we define *f* and the initial condition *u*0. This general notation is commonly employed in ODE textbooks, and we can easily make the Python implementation just as general. The use of NumPy arrays and vectorized computations greatly simplifies the generalization process and enhances the efficiency of our ODE solvers.

### **1.4 A ForwardEuler Class for Systems of ODEs**

The ForwardEuler\_v0 class above was written for scalar ODEs, and we now want to modify it to handle a system of equations: *u* 0 = *f*, *u*(0) = *u*0, where *u*, *f* and *u*0 are vectors (arrays). To identify how the code needs to be changed, let us first revisit the underlying numerical method. Using the general notation introduced earlier, when we apply the FE method to a system of ODEs, the update formula looks exactly the same as in the scalar case, but with all the terms being vectors:

$$
\underbrace{u\_{k+1}}\_{\text{vector}} = \underbrace{u\_k}\_{\text{vector}} + \Delta t \underbrace{f(u\_k, t\_k)}\_{\text{vector}}.
$$

We could also write this formula in terms of the individual components, as in

$$u\_{k+1}^{(i)} = u\_k^{(i)} + \Delta t f^{(i)}(t\_k, u\_k), \text{ for } i = 0, \dots, m-1, \dots$$

but the compact vector notation is more readable. Fortunately, the way we write the vector version of the formula is also how NumPy arrays are used in calculations. The Python code for the formula above may therefore look identical to the version for scalar ODEs:

u[k + 1] = u[k] + dt \* f(t[k], u[k])

with the crucial difference that both u[k], u[k+1], and f(t[k], u[k]) are now arrays.5 Since these are arrays, the solution u must be a two-dimensional array, and u[k],u[k+1], etc. are the rows of this array. The function f expects an array as its second argument, and must return a one-dimensional array containing all the right-hand sides *f* (0)*,...,f*(*n*−1). To gain a better feel for how these arrays look and how they are used, let us compare the array holding the solution of a scalar ODE with that of a system of two ODEs. For the scalar equation, both t and u are one-dimensional NumPy arrays, and indexing into u gives us numbers representing the solution at each time step. For instance, in an interactive Python session we may have arrays t and u with the following contents:

>>> t array([0. , 0.4, 0.8, 1.2, ... ]) >>> u array([1. , 1.4, 1.96, 2.744, ... ])

and indexing into u then gives

>>> u[0] 1.0

5This compact notation requires that the solution vector u is represented by a NumPy array. We could, in principle, use lists to hold the solution components, but the resulting code would need to loop over the components and would be far less elegant and readable.

>>> u[1] 1.4

In the case of a system of two ODEs, the array t remains one-dimensional, but the solution array u becomes two-dimensional, with one column for each solution component. We can index it in the same way as described earlier, and the result is a one-dimensional array of length two, containing the two solution components at a specific time step:

```
>>> u
array([[1.0, 0.8],
        [1.4, 1.1],
        [1.9, 2.7],
        ... ])
>>> u[0]
array([1.0, 0.8])
>>> u[1]
array([1.4, 1.1])
```
Equivalently, we could write

```
>>> u[0,:]
array([1.0, 0.8])
>>> u[1,:]
array([1.4, 1.1])
```
to make explicit which of the two array dimensions (or *axes*) that we are indexing into.

The similarity between the generic mathematical notation for vector and scalar ODEs, as well as the convenient algebra of NumPy arrays, suggests that the implementation of the solver for scalar and system ODEs can be very similar. Indeed, this is true, and the ForwardEuler\_v0 class introduced earlier can be modified with a few minor adjustments to work for ODE systems:

If these two aspects are handled and initialized correctly, the remaining code from Section 1.2 will work without any modifications.

The extended class implementation may look like:

```
import numpy as np
class ForwardEuler:
    def __init__(self, f):
        self.f = lambda t, u: np.asarray(f(t, u), float)
```
6This step is not strictly needed, since we could use a two-dimensional array with shape (N + 1, 1) for scalar ODEs. However, using a one-dimensional array for scalar ODEs gives simpler and more intuitive indexing.

```
def set_initial_condition(self, u0):
   if np.isscalar(u0): # scalar ODE
       self.neq = 1 # no of equations
       u0 = float(u0)
   else: # system of ODEs
       self.neq = u0.size # no of equations
       u0 = np.asarray(u0)
   self.u0 = u0
def solve(self, t_span, N):
   """Compute solution for
   t_span[0] <= t <= t_span[1],
   using N steps."""
   t0, T = t_span
   self.dt = (T - t0) / N
   self.t = np.zeros(N + 1)
   if self.neq == 1:
       self.u = np.zeros(N + 1)
   else:
       self.u = np.zeros((N + 1, self.neq))
   msg = "Please set initial condition before calling solve"
   assert hasattr(self, "u0"), msg
   self.t[0] = t0
   self.u[0] = self.u0
   for n in range(N):
       self.n = n
       self.t[n + 1] = self.t[n] + self.dt
       self.u[n + 1] = self.advance()
   return self.t, self.u
def advance(self):
   """Advance the solution one time step."""
   u, dt, f, n, t = self.u, self.dt, self.f, self.n, self.t
   return u[n] + dt * f(t[n], u[n])
```
It is worth commenting on certain parts of this code. First, the constructor is almost identical to the scalar case, but we use a lambda function and the convenient np.asarray function to convert any f that returns a list or tuple into a function that returns a NumPy array. If f already returns an array, np.asarray will simply return this array with no changes. This modification is not strictly necessary, since we could just assume that the user implements f to return an array. However, it enhances the robustness and flexibility of the class. We have also used the function isscalar from NumPy in the set\_initial\_condition method, to check if u0 is a single number or a NumPy array. This allows us to determine the number of equations self.neq and ensures the class can handle both scalar and system ODEs. The final modification can be observed in the solve method, where the self.neq attribute is inspected. Depending on its value, u is initialized as a one- or two-dimensional array with the appropriate size. The actual for-loop and the advance method remain unchanged from the previous version of the class.

**Example: ODE Model for a Pendulum.** As an example, let us consider a system of ODEs that models the motion of a simple pendulum, as illustrated in Figure 1.3. This nonlinear system is a classic physics problem, and despite its simplicity, it is not possible to find an exact analytical solution. The system is formulated in terms of two main variables; the angle *θ* and the angular velocity *ω*, see Figure 1.3. For a simple pendulum with no friction, the dynamics of these variables are governed by

$$\frac{d\theta}{dt} = \omega,\tag{1.3}$$

$$\frac{d\omega}{dt} = -\frac{g}{L}\sin(\theta),\tag{1.4}$$

where *L* denotes the length of the pendulum and *g* represents the gravitational constant. Eq. (1.3) follows directly from the definition of angular velocity, while (1.4) follows from Newton's second law, where *dω/dt* is the acceleration and the right-hand side is the tangential component of the gravitational force acting on the pendulum, divided by its mass. To solve the system we need to define initial conditions for *θ* and *ω*, i.e., we need to know the initial position and velocity of the pendulum.

**Fig. 1.3** Illustration of the pendulum problem. The main variables of interest are the angle *θ* and its derivative *ω* (the angular velocity).

Since the right-hand side defined by (1.3)-(1.4) includes the parameters *L* and *g*, it is convenient to implement it as a class, similar to the logistic growth model discussed earlier. A possible implementation may look like this:

```
from math import sin
class Pendulum:
    def __init__(self, L, g=9.81):
        self.L = L
```

```
self.g = g
def __call__(self, t, u):
    theta, omega = u
    dtheta = omega
    domega = -self.g / self.L * sin(theta)
    return [dtheta, domega]
```
We observe that the function returns a list. However, this list will be automatically wrapped into a function returning an array by the constructor of the solver class, as mentioned above. The main program remains quite similar to the examples presented earlier, with the exception that we now need to define an initial condition with two components. Assuming that this class definition as well as the ForwardEuler exist in the same file, the code to solve the pendulum problem can look like this:

```
import matplotlib.pyplot as plt
problem = Pendulum(L=1)
solver = ForwardEuler(problem)
solver.set_initial_condition([np.pi / 4, 0])
T = 10
N = 1000
t, u = solver.solve(t_span=(0, T), N=N)
plt.plot(t, u[:, 0], label=r'$\theta$')
plt.plot(t, u[:, 1], label=r'$\omega$')
plt.xlabel('t')
plt.ylabel(r'Angle ($\theta$) and angular velocity ($\omega$)')
plt.legend()
plt.show()
```
Notice that in order to extract and plot each solution component, we need to index into the second dimension of u, using array slicing. If we were to use the first index, such as u[0] or u[0,:], it would return an array of length two containing the solution components at the first time point. In this specific example, a call like plt.plot(t, u) would also work and would plot both solution components. However, there are cases where we are interested in plotting specific components of the solution, and in such cases, array slicing becomes necessary. The resulting plot is shown in Figure 1.4. Additionally, it is worth mentioning the use of Python's raw string format for the labels, indicated by the r in front of the string. Raw strings treat the backslash (\) as a regular character and are often needed when using LaTeX encoding for mathematical symbols. Furthermore, an observant reader may notice that the amplitude of the pendulum motion appears to increase over time, which is clearly not physically accurate. In reality, for an undamped pendulum problem defined by equations (1.3)-(1.4), the energy is conserved, and the amplitude should remain constant. The increasing amplitude is a numerical artifact introduced by the FE method, and the solution may be improved by reducing the time step or using a different numerical method.

**Fig. 1.4** Solution of the simple pendulum problem, computed with the forward Euler method.

### **1.5 Checking the Error in the Numerical Solution**

Recall from Section 1.1 that we derived the FE method by approximating the derivative using a finite difference formula:

$$u'(t\_n) \approx \frac{u(t\_{n+1}) - u(t\_n)}{\Delta t}.\tag{1.5}$$

This approximation obviously introduces an error, and since we approach the true derivative as *∆t* → 0, it is intuitive that the error depends on the size of *∆t*. We visually demonstrated this relationship in Figure 1.1, but it would be valuable to have a way of more precisely quantifying how the error depends on the time step. Analyzing the error in numerical methods is a broad field within applied mathematics, which we will not cover in detail here, and the interested reader is referred to, for instance, [8]. However, when implementing a numerical method it is very useful to know its theoretical accuracy, and in particular to be able to compute the error and verify that the method performs as expected.

The Taylor expansion, which is also discussed briefly in Appendix A.4, is an essential tool for estimating the error in numerical methods for ODEs. For a smooth function *u*ˆ(*t*), if we can compute the function value and its derivatives at time *tn*, we can approximate the value at *tn* + *∆t* using the following series:

$$
\hat{u}(t\_n + \Delta t) = \hat{u}(t\_n) + \Delta t \hat{u}'(t\_n) + \frac{\Delta t^2}{2} \hat{u}''(t\_n) + \frac{\Delta t^3}{6} \hat{u}'''(t\_n) + O(\Delta t^4).
$$

We can include as many terms as we like, and since *∆t* is small we always have *∆t*(*n*+1) *∆tn*, so the error in the approximation is dominated by the first neglected term. The update formula for the FE method, derived from (1.5), is

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

We can recognize this formula as a Taylor series truncated after the first order term, and we expect the error |*un*+1 −*u*ˆ*n*+1| to be proportional to *∆x*2 . Since this is the error for a single time step, the accumulated error after *N* ∼ 1*/∆t* steps is proportional to *∆t*, and the FE method is hence a *first order* method. As we will see in Chapter 2, more accurate methods can be constructed by deriving update formulas that make more terms in the Taylor expansion of the error cancel. This process is fairly straightforward for low-order methods, e.g., of second or third order, but it quickly gets complicated for high order solvers, see, for instance [8] for details.

Knowing the theoretical accuracy of an ODE solver is important for a number of reasons, including the verification of the solver implementation. If we can solve a given problem and demonstrate that the error behaves as predicted by the theory, we gain confidence in the correctness of our solver. To illustrate this procedure, let us consider the simple initial value problem introduced earlier:

$$u' = u, \quad u(0) = 1.$$

As stated above, the analytical solution to this problem is *u* = *e t* , and we can use this to compute the error in our numerical solution. But how should the error be defined? There is no unique answer to this question. For practical applications, common error measures include the root-mean-square (RMS) or relative-root-mean-square (RRMS), which are defined by

$$RMSE = \sqrt{\frac{1}{N} \sum\_{n=0}^{N} (u\_n - \hat{u}(t\_n))^2},$$

$$RRMSE = \sqrt{\frac{1}{N} \sum\_{n=0}^{N} \frac{(u\_n - \hat{u}(t\_n))^2}{\hat{u}(t\_n)^2}},$$

respectively. Here, *un* is the numerical solution at time step *n* and *u*ˆ(*tn*) the corresponding exact solution. In more mathematically oriented texts, the errors are usually defined in terms of *norms*, for instance the discrete *l*1*,l*2*,* and *l*∞ norms:

$$\begin{aligned} e\_{l\_1} &= \sum\_{i=0}^{N} (|u\_i - \hat{u}(t\_i)|), \\ e\_{l\_2} &= \sum\_{i=0}^{N} (u\_i - \hat{u}(t\_i)^2), \\ e\_{l\_\infty} &= \max\_{i=0}^{N} (\hat{u}\_i - u(t\_i)). \end{aligned}$$

While the choice of error norm may be important for certain cases, it is usually not crucial for practical applications, and all the different error measures can generally be expected to behave as predicted by the theory. For simplicity, we will use an even simpler error measure in our example, where we compute the error at the final time *T*, given by *e* = |*uN* −*u*ˆ(*tN* )|. Using the ForwardEuler class introduced above, the complete code for checking the convergence can be written as follows:

```
from forward_euler_class_v1 import ForwardEuler
import numpy as np
def rhs(t, u):
    return u
def exact(t):
    return np.exp(t)
solver = ForwardEuler(rhs)
solver.set_initial_condition(1.0)
T = 3.0
t_span = (0,T)
N = 30
print('Time step (dt) Error (e) e/dt')
for _ in range(10):
    t, u = solver.solve(t_span, N)
    dt = T / N
    e = abs(u[-1] - exact(T))
    print(f'{dt:<14.7f} {e:<12.7f} {e/dt:5.4f}')
    N = N * 2
```
Most of the lines in the code are identical to the previous programs. However, we have enclosed the call to the solve method within a for loop, and the last line ensures that the number of time steps N is doubled for each iteration of the loop. Also, note the f-string format specifiers used, such as {dt:<14.7f}, which specifies that the output should be a left-aligned decimal number with seven decimals, occupying a total of 14 characters. These format specifiers ensure that the numbers are displayed as vertically aligned columns, improving readability, which may be important for visually inspecting the convergence. See, for instance, [16] for a brief introduction to f-strings and format specifiers. The program will produce the following output:

In the rightmost column we observe that the ratio of the error to the time step remains approximately constant. This observation supports the theoretical result that the error is proportional to *∆t*. In the upcoming chapters, we will perform similar calculations for higher order methods to verify that the error is proportional to *∆tr* , where *r* is the theoretical order of convergence for the method.

To compute the error in our numerical solution, we need to determine the true solution to our initial value problem. This task was straightforward for the simple example above because we knew the analytical solution to the equation. However, for more complex ODE problems, estimating the error and the order of convergence requires a different approach. Several alternatives exist, including methods like the *method of manufactured solutions*, where we choose a solution function *u*(*t*) and compute its derivative analytically to determine the right-hand side of the ODE. An even simpler approach, which usually yields good results, involves computing a highly accurate numerical solution using a high-order solver and small time steps. This solution then serves as the reference for computing the error. To obtain accurate error estimates it is essential that the reference solution is significantly more accurate than the numerical solution we want to evaluate. Generating the reference solution typically requires very small time steps and can take some time to compute, but in most cases the computation time for the reference solution is not a significant issue.

### **1.6 Using ODE Solvers from SciPy**

As mentioned in the book's preface, there exist many ODE solvers available for direct use, and it can be argued that there is no need to implement our own solvers. While there is some truth to this, as we have emphasized, it can be beneficial to understand the inner workings of these solvers, in order to apply and use them correctly, and the best way to obtain this knowledge is to implement the solvers ourselves. However, when we have a specific ODE model that we need to solve as efficiently as possible, there are several existing solvers to choose from. For Python programmers, a natural choice may be the solvers provided by *SciPy*, which have evolved into a robust and fairly efficient suite of ODE solvers. SciPy is a comprehensive scientific software package for Python that includes various tools for tasks such as linear algebra, optimization, integration, and more.7 To solve initial value problems, the recommended tool is the solve\_ivp function from the integrate module in SciPy. The following code demonstrates the use of solve\_ivp with the Pendulum class presented above to solve the simple pendulum problem defined by (1.3)-(1.4). We here assume that the Pendulum class is stored in a separate file named pendulum.py.

```
from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
from pendulum import Pendulum
problem = Pendulum(L = 1)
t_span = (0, 10.0)
u0 = (np.pi/4, 0)
solution = solve_ivp(problem, t_span, u0)
plt.plot(solution.t, solution.y[0,:])
plt.plot(solution.t, solution.y[1,:])
plt.legend([r'$\theta$',r'$\omega$'])
plt.show()
```
Running this code will generate a plot similar to Figure 1.5, and we observe that the solution does not appear as smooth as the one obtained from the ForwardEuler solver introduced earlier. This seeming discrepency is due to the nature of the solve\_ivp solver, which is an *adaptive* solver that automatically selects the time step to meet a specified error tolerance. The default value of this tolerance is relatively large, leading to the solver using very few time steps and resulting in jagged-looking solution plots. Comparing the plot with a highly accurate numerical solution, represented by the two dotted curves in Figure 1.5, we notice that the solution at the specified time points *tn* is fairly accurate. However, the visual appearance is compromised by the linear interpolation between these time points. To obtain a more visually appealing solution, there are several approaches we can take. We may, for instance, pass the function an additional argument t\_eval, which is a NumPy array containing the desired time points for evaluating the solution:

t\_eval = np.linspace(0, 10.0, 1001)

7See https://scipy.org/

**Fig. 1.5** Solution of the simple pendulum problem, computed with the SciPy solve\_ivp function and the default tolerance.

```
solution = solve_ivp(problem, t_span, u0, t_eval=t_eval)
```
Alternatively, we can reduce the error tolerance of the solver, for instance, by setting

```
rtol = 1e-6
solution = solve_ivp(problem, t_span, u0, rtol=rtol)
```
This latter call will reduce the relative tolerance rtol from its default value of 1e-3 (0.001). We could also adjust the absolute tolerance using the parameter atol. We will not cover all the possible arguments and options to the solve\_ivp function here, but it is worth mentioning that we can also change the numerical method used by the function, by passing in a parameter named method. For instance, a call like

```
rtol = 1e-6
solution = solve_ivp(problem, t_span, u0, method='Radau')
```
will replace the default solver (called rk45) with an implicit Radau ODE solver, which we will cover in Chapter 3. For a complete description of parameters accepted by the solve\_ivp function we recommend referring to the online documentation available on the SciPy website.

**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 license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license 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.

# **Chapter 2 Improving the Accuracy**

As mentioned previously, the FE method derived in Chapter 1 is not the most sophisticated ODE solver. Although it provides sufficient accuracy for most of the applications covered in this book, there are alternative methods available that offer improved accuracy and stability, making them better suited for solving challenging ODE systems. In this chapter, we will focus on enhancing accuracy, and thus we will primarily explore explicit methods. Implicit methods, which exhibit superior stability properties and are better suited for solving *stiff* ODEs, will be discussed in Chapter 3.

In Chapter 1 we demonstrated that the FE method is a first-order accurate method, meaning that the error in the numerical solution is proportional to the size of the time step *∆t*. In this chapter, we will derive solvers of higher order, where the numerical error is proportional to a higher power of *∆t*. To explain the derivation of such higher order ODE solvers, we will revisit the general formulation of the ODE system:

$$u' = f(t, u), \quad u(t\_0) = u\_0.$$

Instead of simply replacing the derivative with a finite difference approximation, as done in the derivation of the FE method, we will adopt a slightly different approach in this chapter. Assuming that we know the solution *un* at time *tn*, we can find the solution at time *tn*+1 by integrating both sides of the equation. We have

$$\int\_{t\_n}^{t\_{n+1}} \frac{du}{dt} dt = \int\_{t\_n}^{t\_{n+1}} f(t, u(t)) dt,$$

which gives us the exact solution at time *tn*+1 as

$$u(t\_{n+1}) = u(t\_n) + \int\_{t\_n}^{t\_{n+1}} f(t, u(t))dt. \tag{2.1}$$

© The Author(s) 2024 23 J. Sundnes, *Solving Ordinary Differential Equations in Python*, Simula SpringerBriefs on Computing 15, https://doi.org/10.1007/978-3-031-46768-4\_2

In general, it is not feasible to compute the integral on the right-hand side of the equation analytically, as the function *f* is often nonlinear and the function *u*(*t*) is unknown. However, we can approximate the integral using various numerical integration techniques. The simplest approximation is to set

$$f(t, u(t)) \approx f(t\_n, u\_n), \text{ for } t\_n < t < t\_{n+1},$$

meaning that we approximate the integrand as a constant over the interval from *tn* to *tn* + *∆t*. Substituting this choice into (2.1) yields the update formula:

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

which is recognized as the FE method introduced in Chapter 1. Approximating the function *f*(*tn,un*) as constant over the interval *tn < t < tn*+1 is obviously not the most accurate choice, and we will explore more accurate approximations that lead to higher-order ODE solvers.

The classical approach for approximating the integral of a general nonlinear function is to approximate the function with a polynomial and then integrate this polynomial analytically. This technique forms the basis for standard quadrature rules in numerical integration, and has also been used to derive accurate ODE solvers. Two main ideas have been explored for constructing the polynomial approximation of *f*(*t,u*), resulting in two important classes of ODE solvers. The first approach is to approximate *f*(*t,u*) with a polynomial that interpolates *f* at previous time points, such as *f*(*tn*−1*,un*−1)*,f*(*tn*−2*,un*−2)*,...*. This method gives rise to *multistep* methods, which are widely used for solving ODEs. We will not cover multistep methods in this book, but interested readers can refer to references [1, 8] for further details. The second approach entails computing a number of intermediate approximations of *f*(*t,u*) on the interval *tn < t < tn*+1, and using these values to define the polynomial approximation of *f*. This approach is analogous to the derivation of classical quadrature rules in numerical integration, and leads to a class of ODE solvers known as Runge-Kutta methods. These methods come in many forms, exhibiting different accuracy and stability properties, and will be the main focus of Chapters 2-4.

### **2.1 Explicit Runge-Kutta Methods**

As mentioned earlier, an intuitive way to improve the accuracy of the approximate integral in (2.1) is to calculate several intermediate approximations of *f*(*t*∗*,u*∗) for *tn* ≤ *t*∗ ≤ *tn*+1, and calculate the integral as a weighted sum of these values. This approach builds upon standard numerical integration techniques and gives rise to a widely used class of ODE solvers called *Runge-Kutta (RK) methods*. The simplest example of an RK method is in fact the FE method discussed earlier, which is a one-stage, first-order, explicit RK method. An alternative formulation of the FE method is

$$\begin{aligned} k\_1 &= f(t\_n, u\_n), \\ u\_{n+1} &= u\_n + \Delta t k\_1. \end{aligned}$$

It can be observed that this is the same formula as introduced earlier, and there is no real advantage in writing the formula in two lines instead of one. However, this alternative formulation aligns with the typical representation of RK methods and facilitates understanding the relationship between the FE method and more advanced solvers. The intermediate value *k*1 is often referred to as a *stage derivative* in the ODE literature.

To enhance the accuracy of the FE method to second order, i.e., with error proportional to *∆t*2 , we can employ more accurate approximations of the integral in (2.1). One option is to maintain the assumption that *f*(*t,u*(*t*)) is constant over *tn* ≤ *t*∗ ≤ *tn*+1, but to approximate it at the midpoint of the interval instead of the left end. This approach requires one additional stage:

$$k\_1 = f(t\_n, u\_n),\tag{2.2}$$

$$k\_2 = f(t\_n + \frac{\Delta t}{2}, u\_n + \frac{\Delta t}{2}k\_1),\tag{2.3}$$

$$
u\_{n+1} =
u\_n +
\Delta t \, k\_2.\tag{2.4}$$

This method is known as the explicit midpoint method or the modified Euler method. The first step is identical to that of the FE method, but instead of using the stage derivative *k*1 to advance the solution to the next step, we use it to calculate an intermediate midpoint solution

$$u\_{n+1/2} = u\_n + \frac{\Delta t}{2}k\_1.$$

This solution is then used to compute the corresponding stage derivative *k*2, which serves as an approximation to the derivative of *u* at time *tn* + *∆t/*2. Finally, we use this midpoint derivative to advance the solution to *tn*+1.

Another second-order method is Heun's method, also known as the explicit trapezoidal method, which can be derived by approximating the integral in equation (2.1) using the trapezoidal rule:

$$k\_1 = f(t\_n, u\_n),\tag{2.5}$$

$$k\_2 = f(t\_n + \Delta t, u\_n + \Delta t k\_1),\tag{2.6}$$

$$u\_{n+1} = u\_n + \frac{\Delta t}{2}(k\_1 + k\_2). \tag{2.7}$$

This method also computes two stage derivatives *k*1 and *k*2. However, note that the formula for *k*2 approximates the derivative at *tn*+1 rather than at the midpoint *tn* +*∆t/*2. The solution is then advanced from *tn* to *tn*+1 using the mean value of *k*1 and *k*2.

All RK methods follow the same recipe as the two second-order methods considered above; we calculate one or more intermediate values (i.e., stage derivatives) and then advance the solution using a combination of these stage derivatives. The method's accuracy can be improved by adding more stages. A general RK method with *s* stages can be written as

$$k\_i = f(t\_n + c\_i \Delta t, u\_n + \Delta t \sum\_{j=1}^s a\_{ij} k\_j), \text{ for } i = 1, \dots, s \tag{2.8}$$

$$u\_{n+1} = u\_n + \Delta t \sum\_{i=1}^{s} b\_i k\_i. \tag{2.9}$$

Here *ci , bi , aij* , for *i, j,*= 1*,...,s* are coefficients specific to the method. Every RK method can be written in this form, and a method is uniquely determined by the number of stages, *s*, and the values of the coefficients.

As mentioned earlier, there exists a wide variety of RK methods, where the coefficients are typically chosen to optimize the accuracy for a given number of stages. While we will not delve into the details of how the methods are constructed here, it can be useful to mention some of the underlying principles. For instance, it can be shown that the P *bi* coefficients must satisfy *s i*=1 *bi* = 1 in order to ensure convergence. This condition naturally arises from the motivation for RK methods as numerical integrators applied to equation (2.1). When approximating the integral as a weighted sum, the sum of the weights must be equal to one. Another common constraint on the coefficients is to set *ci* = P*s j*=1 *aij* . While not strictly necessary, this constraint can simplify the derivation of the methods and aligns with our interpretation of the stage derivative *ki* as approximations of the right-hand side *f*(*t,u*) at time *tn* + *ci∆t*. When implementing a new solver it is easy to introduce errors in the coefficient values, and it may be useful to include tests to verify that the most fundamental conditions on the coefficients are satisfied. It is possible to derive general *order conditions* that the coefficients must satisfy for a method to achieve a given order, see, for instance, [1,8] for details. In this chapter we exclusively consider explicit Runge-Kutta (ERK) methods, which means that *aij* = 0 for *j* ≥ *i*. It can be shown that the order *p* of an ERK method with *s* stages satisfies *p* ≤ *s*, and for *p* ≥ 5, the bound is *p* ≤ *s*−1. However, it remains unknown whether this latter bound is sharp, and it may be even stricter for methods of very high order. For instance, all known methods with *p* = 8 have at least eleven stages, and it is not known whether eight-order methods with nine or ten stages exist.

In the ODE literature, method coefficients are often specified in the form of a *Butcher tableau*, which provides a concise representation of any RK method. The Butcher tableau is simply a specification of all the method coefficients, and for a general RK method it is written as

$$
\begin{array}{c}
 c\_i \left| a\_{11} \cdot \cdots \cdot a\_{1s} \right. \\
 \vdots & \vdots \\
 c\_s \left| a\_{s1} \cdot \cdots \cdot a\_{ss} \right. \\
 \hline
 & b\_1 \cdot \cdots \cdot b\_s \end{array}
$$

The Butcher tableaus of the three methods discussed above: FE, explicit midpoint, and Heun's method, are

$$
\begin{array}{c c c c}
\begin{array}{c c c}
0 & 0 \\
\hline
1
\end{array} & \begin{array}{c c c}
0 & 0 \\
\hline
1/2 & 1/2 \\
\hline
0 & 1
\end{array} & \begin{array}{c c c}
0 & 0 \\
1 & 1 & 0 \\
\hline
1/2 & 1/2
\end{array} \\
\end{array},
$$

respectively. To grasp the concept of Butcher tableaus, it is beneficial to practice inserting the coefficients from these three tableaus into equations (2.8)-(2.9) and verifying that the correct formulas for the respective methods are obtained. As an example of a higher order method, we may consider the "original" Runge-Kutta method, which is a fourth-order, four-stage method defined by

$$
\begin{array}{c|cccc}
0 & 0 & 0 & 0 & 0 \\
1/2 & 1/2 & 0 & 0 & 0 \\
1/2 & 0 & 1/2 & 0 & 0 \\
1 & 0 & 0 & 1 & 0 \\
\hline
\end{array},
$$

which gives the formulas

$$k\_1 = f(t\_n, u\_n),\tag{2.10}$$

$$k\_2 = f(t\_n + \frac{\Delta t}{2}, u\_n + \frac{\Delta t}{2}k\_1),\tag{2.11}$$

$$k\_3 = f(t\_n + \frac{\Delta t}{2}, u\_n + \frac{\Delta t}{2}k\_2),\tag{2.12}$$

$$k\_4 = f(t\_n + \Delta t, u\_n + \Delta t k\_3),\tag{2.13}$$

$$u\_{n+1} = u\_n + \frac{\Delta t}{6} \left(k\_1 + 2k\_2 + 2k\_3 + k\_4\right). \tag{2.14}$$

As mentioned earlier, all the methods discussed in this chapter are explicit methods, meaning that *aij* = 0 for *j* ≥ *i*. Examining equations (2.10)-(2.14) or the general formula (2.8) more closely, we observe that this conditions implies that each stage derivative *ki* only depends on previously computed stage derivatives. Consequently, all *ki* can be computed sequentially using explicit formulas. In contrast, for implicit RK methods, *aij* 6= 0 for some *j* ≥ *i*. As seen in equation (2.8), the formula for computing *ki* will then include *ki* on the right-hand side, as part of the argument to the function *f*. Therefore, equations need to be solved to compute the stage derivatives, and since *f* is typically nonlinear, we need to solve these equations with an iterative solver such as Newton's method. These steps make implicit RK methods more complex to implement and more computationally expensive per time step, but they are also more stable than explicit methods and perform much better for certain classes of ODEs. We will cover implicit RK methods in Chapter 3.

### **2.2 A Class Hierarchy of Runge-Kutta Methods**

We now want to implement RK methods as classes, similar to the FE classes introduced above. Upon examining the ForwardEuler class, we may notice that most of the code is common to all ODE solvers and is not specific to the FE method. For instance, we always need to create an array to store the solution, and the general solution method using a for-loop remains the same for all methods. The only difference among the methods lies in how the solution is advanced from one time step to the next. Recalling the ideas of Object-Oriented Programming, it becomes apparent that a class hierarchy is a suitable structure for implementing such a collection of ODE solvers. This approach allows us to consolidate common code in a superclass (base class), and use inheritance to avoid code duplication. The superclass can handle most of the administrative steps of the ODE solver, such as

To address these aspects, we can introduce a superclass called ODESolver and implement the method-specific details in subclasses. This is precisely why we isolated the code to perform a single step in the advance method, as it becomes the only method that needs to be implemented in the subclasses. The implementation of the superclass can closely resemble the ForwardEuler class introduced earlier:

```
import numpy as np
class ODESolver:
   def __init__(self, f):
       # Wrap user's f in a new function that always
       # converts list/tuple to array (or let array be array)
       self.model = f
       self.f = lambda t, u: np.asarray(f(t, u), float)
   def set_initial_condition(self, u0):
       if np.isscalar(u0): # scalar ODE
          self.neq = 1 # no of equations
          u0 = float(u0)
       else: # system of ODEs
          u0 = np.asarray(u0)
```

```
self.neq = u0.size # no of equations
    self.u0 = u0
def solve(self, t_span, N):
    """Compute solution for t_span[0] <= t <= t_span[1],
    using N steps."""
    t0, T = t_span
    self.dt = (T - t0) / N
    self.t = np.zeros(N + 1) # N steps ~ N+1 time points
    if self.neq == 1:
        self.u = np.zeros(N + 1)
    else:
        self.u = np.zeros((N + 1, self.neq))
    msg = "Please set initial condition before calling solve"
    assert hasattr(self, "u0"), msg
    self.t[0] = t0
    self.u[0] = self.u0
    for n in range(N):
        self.n = n
        self.t[n + 1] = self.t[n] + self.dt
        self.u[n + 1] = self.advance()
    return self.t, self.u
    def advance(self):
        raise NotImplementedError(
            "Advance method is not implemented in the base class")
```
It is important to note that the ODESolver is designed to be a pure superclass, and the implementation of the advance method is left for subclasses. To clearly convey this abstract nature of the class, we have included an advance method that raises a NotImplementedError when it is called. If an attempt is made to create an instance of ODESolver and use it as a standalone solver, an error will occur in the line self.u[n + 1] = self.advance(). Omitting the definition of advance entirely would lead to an error in the same line, but in this case it would be a less informative AttributeError. By raising the NotImplementedError, it explicitly indicates to anyone reading or using the code that this behavior is intentional and that the specific functionality is meant to be implemented in subclasses. It is worth mentioning that there are alternative approaches in Python to make explicit the abstract nature of the ODESolver class, for instance using the module abc, for "Abstract Base Class". However, while this solution may be considered more modern, we have decided to not use it here, in the interest of keeping the code simple and compact.

With the superclass at hand, the implementation of a ForwardEuler subclass becomes very simple:

```
class ForwardEuler(ODESolver):
    def advance(self):
```

```
u, f, n, t = self.u, self.f, self.n, self.t
dt = self.dt
return u[n] + dt * f(t[n], u[n])
```
Similarly, the explicit midpoint method and the fourth-order RK method can be subclasses, each implementing a single method:

```
class ExplicitMidpoint(ODESolver):
    def advance(self):
        u, f, n, t = self.u, self.f, self.n, self.t
        dt = self.dt
        dt2 = dt / 2.0
        k1 = f(t[n], u[n])
        k2 = f(t[n] + dt2, u[n] + dt2 * k1)
        return u[n] + dt * k2
class RungeKutta4(ODESolver):
    def advance(self):
        u, f, n, t = self.u, self.f, self.n, self.t
        dt = self.dt
        dt2 = dt / 2.0
        k1 = f(t[n], u[n],)
        k2 = f(t[n] + dt2, u[n] + dt2 * k1, )
        k3 = f(t[n] + dt2, u[n] + dt2 * k2, )
        k4 = f(t[n] + dt, u[n] + dt * k3, )
        return u[n] + (dt / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
```
The use of these classes is nearly identical to the FE class introduced in Section 1.3. Considering the same simple ODE used above; *u* 0 = *u*, *u*(0) = 1, *t* ∈ [0*,*3], *∆t* = 0*.*5, the code looks like:

```
import numpy as np
import matplotlib.pyplot as plt
from ODESolver import ForwardEuler, ExplicitMidpoint, RungeKutta4
def f(t, u):
    return u
t_span = (0, 3)
N = 6
fe = ForwardEuler(f)
fe.set_initial_condition(u0=1)
t1, u1 = fe.solve(t_span, N)
plt.plot(t1, u1, label='Forward Euler')
em = ExplicitMidpoint(f)
em.set_initial_condition(u0=1)
t2, u2 = em.solve(t_span, N)
plt.plot(t2, u2, label='Explicit Midpoint')
rk4 = RungeKutta4(f)
rk4.set_initial_condition(u0=1)
t3, u3 = rk4.solve(t_span, N)
```

```
plt.plot(t3, u3, label='Runge-Kutta 4')
# plot the exact solution in the same plot
time_exact = np.linspace(0, 3, 301)
plt.plot(time_exact, np.exp(time_exact), label='Exact')
plt.title('RK solvers for exponential growth, $\\Delta t = 0.5$')
plt.xlabel('$t$')
plt.ylabel('$u(t)$')
plt.legend()
plt.show()
```
This code will solve the same simple equation using three different methods, and plot the solutions in the same window, as shown in Figure 2.1. To emphasize the disparity in accuracy between the methods, we have set N = 6, resulting in a very large time step (*∆t* = 0*.*5).

**Fig. 2.1** Numerical solutions of the exponential growth problem, computed with ForwardEuler, ImplicitMidpoint and RungeKutta4. All the solvers use *∆t* = 0*.*5, to highlight the difference in accuracy.

### **2.3 Testing the Solvers**

In Chapter 1 we demonstrated how to compute the error in the numerical solution, and in particular how we could verify that the error behaved as predicted by the theoretical convergence of the applied solvers. Such tests are extremely valuable for verifying that we have implemented the ODE solvers correctly, and can easily be extended to the higher order solvers. As an example, the following code defines a dictionary containing three different solver classes and their theoretical order, and solves the simple exponential ODE using all three solvers.

```
from ODESolver import *
import numpy as np
def rhs(t, u):
    return u
def exact(t):
    return np.exp(t)
solver_classes = [(ForwardEuler,1), (Heun,2),
                  (ExplicitMidpoint,2), (RungeKutta4,4)]
for solver_class, order in solver_classes:
    solver = solver_class(rhs)
    solver.set_initial_condition(1.0)
    T = 3.0
    t_span = (0, T)
    N = 30
    print(f'{solver_class.__name__}, order = {order}')
    print(f'Time step (dt) Error (e) e/dt**{order}')
    for _ in range(10):
        t, u = solver.solve(t_span, N)
        dt = T / N
        e = abs(u[-1] - exact(T))
        if e < 1e-13: # break if error is close to machine precision
            break
        print(f'{dt:<14.7f} {e:<12.7f} {e/dt**order:5.4f}')
        N = N * 2
```
The code is nearly identical to the FE convergence test in Section 1.5, with the only difference being that we loop over a list of tuples containing the four method classes and their corresponding orders. The output is also similar to the previous version, but now repeated for all four solvers. The built-in class attribute \_\_name\_\_ is used to extract and print the name of each solver. Three columns are displayed, representing the time step *∆t*, the error *e* at time *t* = 3*.*0, and finally *e/∆tp* , where *p* is the order of the method. The output matches the expected values for the first two methods, as the numbers in the rightmost column are approximately constant, confirming that the error is in fact proportional to *∆tp* . However, the last part of the output, for the forth order RK method, looks like this:

We see that the *e/∆tp* numbers remain approximately constant for a while, consistent with the convergence order of the method. However, for very small values of *∆t*, these values start to increase. This behavior is not uncommon in convergence tests, especially for high-order methods, and is caused by the finite accuracy of number representation on computers. As the numerical errors become smaller and approach the machine precision (≈ 10−16), roundoff errors start to dominate the overall error, leading to a loss of convergence.

There are many alternative ways to check the implementation of ODE solvers. One approach is to consider an even simpler ODE where the righthand side is a constant, i.e., *u* 0 (*t*) = *f*(*t,u*) = *C*. The solution to this simple ODE is given by *u*(*t*) = *Ct*+*u*0, where *u*0 is the initial condition. All the numerical methods discussed in this book will capture this solution to machine precision, and we can write a general test function that takes advantage of this:

```
def test_exact_numerical_solution():
    solver_classes = [ForwardEuler, Heun,
                       ExplicitMidpoint, RungeKutta4]
    a = 0.2
    b = 3
    def f(t, u):
        return a
    def u_exact(t):
        """Exact u(t) corresponding to f above."""
        return a * t + b
    u0 = u_exact(0)
    T = 8
    N = 10
    tol = 1E-14
    t_span = (0, T)
    for solver_class in solver_classes:
        solver = solver_class(f)
```

```
solver.set_initial_condition(u0)
t, u = solver.solve(t_span, N)
u_e = u_exact(t)
max_error = abs((u_e - u)).max()
msg = f'{solver_class.__name__} failed, error={max_error}'
assert max_error < tol, msg
```
Similar to the convergence check illustrated below, this code will loop through all the solver classes, solve the simple ODE, and check that the resulting error falls within the specified tolerance.

Both of the methods shown here for verifying the implementation of our solvers have certain limitations. The most important one is that they both solve very simple ODEs, and it is possible to introduce errors in the code that may only manifest themselves when dealing with more complex problems. However, the methods presented here offer the advantages of simplicity and generality, and they can be applied to any newly implemented ODE solver class. Many common implementation errors, such as incorrectly specifying a single parameter in an RK method, will often become apparent even when solving these simple problems. Therefore, these methods can provide an initial indication of whether the implementation is correct, which can be followed by more extensive tests if needed.

**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 license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license 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.

# **Chapter 3 Stable Solvers for Stiff ODE Systems**

In the previous chapter, we introduced explicit Runge-Kutta (ERK) methods and demonstrated how they can be implemented as a hierarchy of Python classes. For most ODE systems, replacing the simple forward Euler method with a higher-order ERK method will significantly reduce the number of time steps needed to reach a specified accuracy. Furthermore, it often leads to reduced computation time, since the additional cost per time step is outweighed by the reduced number of steps. However, there exists a class of ODEs known as *stiff* systems, where all the ERK methods require very small time steps, and any attempt to increase the time step leads to spurious oscillations and possible divergence of the solution. Stiff ODE systems pose a challenge for explicit methods, and they are better addressed by implicit solvers such as implicit Runge-Kutta (IRK) methods. IRK methods are well-suited for stiff problems and can offer substantial reductions in computation time when tackling challenging problems.

### **3.1 Stiff ODE Systems and Stability**

One well-known example of a stiff ODE system is the Van der Pol equation, which can be written as an initial value problem on the form

$$y\_1' = y\_2,\qquad\qquad y\_1(0) = 1,\tag{3.1}$$

$$y\_2' = \mu (1 - y\_1^2) y\_2 - y\_1, \quad y\_2(0) = 0. \tag{3.2}$$

Here, the parameter *µ* represents a constant that determines the properties of the system, including its "stiffness". When *µ* = 0 the problem is a simple oscillator with analytical solution *y*1 = cos(*t*)*,y*2 = sin(*t*). However, for non-zero values of *µ*, the solution exhibits far more complex behavior. The following code implements this system and solves it with the ForwardEuler subclass from the ODESolver class hierarchy.

```
from ODESolver import *
import numpy as np
import matplotlib.pyplot as plt
class VanderPol:
    def __init__(self, mu):
        self.mu = mu
    def __call__(self, t, u):
        du0 = u[1]
        du1 = self.mu * (1 - u[0]**2) * u[1] - u[0]
        return du0, du1
model = VanderPol(mu=1)
solver = ForwardEuler(model)
solver.set_initial_condition([1, 0])
t,u = solver.solve(t_span=(0, 20), N=1000)
plt.plot(t, u)
plt.show()
```
Figure 3.1 shows the solutions of the Van der Pol equation for *µ* = 0*,*1 and 5. When the parameter *µ* is set even higher, such as *µ* = 50, the solution diverges (becomes unstable) with the given time step (*∆t* = 0*.*02). Although using a more accurate ERK method instead of the FE method may provide some improvement, it does not resolve the issue significantly. It does help to reduce the time step considerably, but the resulting computation time may be substantial. In this problem, the time step is determined by stability requirements rather than the desired accuracy, and opting for a solver that is more stable than the previously discussed ERK methods may yield significant benefits.

Before introducing more stable solvers, it is useful to examine the observed stability problems in more detail. Why does the solution of the Van der Pol model deteriorate significantly for large values of *µ*? More generally, what are the properties of an ODE system that make it stiff? To address these questions, it is useful to start with a simpler problem than the Van der Pol model. Consider, for instance, a simple IVP known as the Dahlquist test equation:

$$u' = \lambda u, \quad u(0) = 1,\tag{3.3}$$

where *λ* can be a complex number.1 When *λ* = 1, it corresponds to the simple exponential growth problem we discussed earlier. However, in this chapter

1Note that the implementation of the solvers in this book does not support solving this ODE for complex *λ*. However, considering complex values in the stability analysis is still important because, for systems of ODEs, the relevant values are the *eigenvalues* of the right-hand side, and these may be complex.

**Fig. 3.1** Solutions of the Van der Pol model for different values of *µ*.

we primarily focus on *λ* values with a negative real part, i.e., either real or complex *λ* values that satisfy <(*λ*) *<* 0. In such cases, the solution of equation (3.3) decays over time and remains stable. However, we will discover that the numerical solutions may not always preserve this stability.

Following the definition in [1], we classify problem (3.3) as stiff for an interval [0*, b*] if the real part of *λ* satisfies

$$b\mathfrak{R}(\lambda) \ll -1.$$

For more general nonlinear problems, such as the Van der Pol model in (3.1)- (3.2), the stiffness of the system is determined by the eigenvalues *λi* of the local Jacobian matrix *J*, which is the matrix of partial derivatives of the right-hand side function *f*. The Jacobian is defined by

$$J\_{ij} = \frac{\partial f\_i(t, y)}{\partial y\_j},$$

and the problem is considered stiff for an interval [0*, b*] if

$$b \min\_i \Re(\lambda\_i) \ll -1.$$

These definitions highlight that the stiffness of a problem depends not only on the ODE itself, but also on the length of the solution interval (*b*), which may seem somewhat surprising. To understand why the interval of interest is important, let us consider the equation (3.3). If *λ* is large and negative, we need to choose a small *∆t* to maintain stability of explicit solvers, as we will discuss in more detail later. However, if our goal is to solve the equation over a very small time interval, i.e., *b* is small, using a small *∆t* is not a problem, and according to the definition above the problem will no longer be considered stiff. In addition to these definitions, the ODE literature also provides more pragmatic definitions of stiffness. For example, an equation is often classified as stiff if the time step needed to maintain stability of an explicit method is much smaller than the time step dictated by the accuracy requirements [1,2]. For a more comprehensive discussion of stiff ODE systems, refer to [1, 9].

Equation (3.3) serves as the foundation for *linear stability analysis*, a valuable technique for analyzing and understanding the stability of ODE solvers. The solution to this equation is given by *u*(*t*) = *e λt*, which grows rapidly if *λ* has a positive real part. Therefore, our primary interest lies in the case where <(*λ*) *<* 0, for which the analytical solution is stable, but our choice of solver may introduce *numerical instabilities*. When the FE method is applied to (3.3), we obtain the update formula

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

and for the first step, with the initial condition *u*(0) = 1, we have

$$u\_1 = 1 + \Delta t \lambda.\tag{3.4}$$

The analytical solution decays exponentially for <(*λ*) *<* 0, and it is natural to require that the numerical solution decreases monotonically. This leads to the requirement |1 + *∆tλ*| ≤ 1. When *λ* is a negative real number, the time step must satisfy *∆t* ≤ −2*/λ* to ensure stability. It is important to note that meeting this stability criterion does not necessarily guarantee a highly accurate solution; the numerical solution may exhibit oscillate and differ substantially from the exact solution. Nevertheless, by selecting *∆t* to satisfy the stability criterion, we ensure that the solution, along with any spurious oscillations or other numerical artifacts, decays with time.

We have observed that the right-hand side of (3.4) contains critical information about the stability of the FE method. This expression is commonly referred to as the *stability function* or *amplification factor* of the method, and is often written as

$$R(z) = 1 + z.$$

For the FE method to be stable, all values of *λ∆t* must satisfy |*R*(*λ∆t*)| *<* 1. This region of *λ∆t* values in the complex plane is referred to as the method's *region of absolute stability* or its *stability region*. The stability region for the FE method is shown in the left panel of Figure 3.2, taking the form of a circle with center at (−1*,*1) and a radius of one. It is evident that if *λ* 0, the requirement for *λ∆t* to lie within this circle is quite restrictive for the choice of *∆t*.

**Fig. 3.2** Stability regions for explicit Runge-Kutta methods. From left: forward Euler, explicit midpoint, and the fourth order method given by (2.10)-(2.14).

We can easily extend the linear stability analysis to the other explicit RK methods introduced in Chapter 2. For instance, applying a single step of the explicit midpoint method given by (2.2)-(2.4) to (3.3) gives

$$u(\Delta t) = 1 + \lambda \Delta t + \frac{(\Delta t \lambda)^2}{2},$$

and we identify the stability function for this method as

$$R(z) = 1 + z + \frac{z^2}{2}.$$

The corresponding stability region is shown in the middle panel of Figure 3.2. For the fourth-order RK method defined in (2.10)-(2.14), the same steps yield the stability function

$$R(z) = 1 + z + \frac{z^2}{2} + \frac{z^3}{6} + \frac{z^4}{24},$$

and the stability region is shown in the right panel of Figure 3.2. We observe that the stability regions of these higher-order RK methods are slightly larger than that of the FE method. However, the difference is not very large, and when also considering the computational cost of each time step, the FE method is usually superior for problems where the time step is governed by stability.

It can be shown that the stability function for an *s*-stage explicit RK method is always a polynomial of degree ≤ *s*, and it can easily be verified that the stability region defined by such a polynomial will never be very large. To obtain a significant improvement of this situation, we need to replace the explicit methods discussed so far with implicit RK methods.

### **3.2 Implicit methods for stability**

Given that equation (3.3) is stable for all values of *λ* with a negative real part, it is natural to look for numerical methods with the same stability property. This implies that the stability domain of the method covers the entire left half of the complex plane, or in other words, that its stability function |*R*(*z*)| ≤ 1 whenever <(*z*) ≤ 0. This property is called *A-stability*. As noted above, the stability function of an explicit RK method is always a polynomial, and no polynomial can satisfy |*R*(*z*)| *<* 1 for all *z <* 0. Hence, there are no Astable explicit RK methods. An even stronger stability requirement can be motivated by the fact that for *λ* 0, the solution to (3.3) decays very rapidly. It is reasonable to expect the same behavior from the numerical solution, which leads to the requirement that |*R*(*z*)| → 0 as *z* → −∞. This property is referred to as *stiff decay*, and an A-stable method that exhibits stiff decay is known as an *L-stable* method. For further detail, readers can refer to [1, 9].

The simplest implicit RK method is the backward Euler (BE) method, which can be derived in exactly the same way as the FE method, by approximating the derivative with a simple finite difference. The only difference from the FE method is that the right-hand side is evaluated at step *n*+ 1 rather than step *n*. For a general ODE, we have

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

and if we rearrange the terms we get

$$u\_{n+1} - \Delta t f(t\_{n+1}, u\_{n+1}) = u\_n. \tag{3.5}$$

Although the derivation is similar to the FE method, there is a fundamental difference. In the BE method, the unknown variable *un*+1 occurs as an argument in the right-hand side function *f*(*t,u*). Therefore, for nonlinear *f*, (3.5) becomes a nonlinear algebraic equation that must be solved to determine *un*+1, instead of having an explicit update formula as in the FE method. This requirement makes implicit methods more complex to implement than explicit methods, and they tend to require more computations per time step. However, as we will demonstrate later, the superior stability properties of implicit solvers still make them better suited for stiff problems.

We will consider the implementation of implicit solvers in Section 3.3 below, but let us first examine the stability of the BE method and other implicit RK solvers using the linear stability analysis introduced earlier. Applying the BE method to (3.3) yields

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

and for the first time step, with *u*(0) = 1, we get

$$u\_1 = \frac{1}{1 - \Delta t \lambda}.$$

The stability function of the BE method is therefore *R*(*z*) = 1*/*(1 − *z*), and its corresponding stability domain is shown in the left panel of Figure 3.3. The method is stable for all choices of *λ∆t outside* the circle with a radius of one and centered at (1,0) in the complex plane. This confirms that the BE method is a highly stable method. It is both A-stable, as its stability domain covers the entire left half of the complex plane, and L-stable, as the stability function satisfies *R*(*z*) → 0 as <(*z*) → −∞.

The BE method fits into the general RK framework defined by (2.8)-(2.9) in Chapter 2, with a single stage (*s* = 1), and *a*11 = *b*1 = *c*1 = 1. Similar to the FE method discussed in Chapter 2, we can reformulate the method slightly to introduce a stage derivative and emphasize its connection to the RK family:

$$k\_1 = f(t\_n + \Delta t, u\_n + \Delta t k\_1),\tag{3.6}$$

$$u\_{n+1} = u\_n + \Delta t k\_1. \tag{3.7}$$

The explicit midpoint and trapezoidal methods mentioned earlier also have their implicit counterparts. The implicit midpoint method is given by

$$k\_1 = f(t\_n + \Delta t/2, u\_n + k\_1 \Delta t/2),\tag{3.8}$$

$$u\_{n+1} = u\_n + \Delta t k\_1,\tag{3.9}$$

while the implicit trapezoidal rule, or Crank-Nicolson method, is given by

$$k\_1 = f(t\_n, u\_n),\tag{3.10}$$

$$k\_2 = f(t\_n + \Delta t, u\_n + \Delta t k\_2),\tag{3.11}$$

$$u\_{n+1} = u\_n + \frac{\Delta t}{2}(k\_1 + k\_2). \tag{3.12}$$

**Fig. 3.3** Stability regions for the backward Euler method (left) and the implicit midpoint method and trapezoidal method (right).

Note that this formulation of the Crank-Nicolson is not very common, and it can be simplified by eliminating the stage derivatives and defining the method in terms of *un* and *un*+1. However, the given formulation in (3.10)- (3.12) highlights its implicit RK nature. The implicit nature of these methods is apparent from the formulas, as one of the stage derivatives must be found by solving an equation involving the nonlinear function *f* instead of using an explicit update formula. The Butcher tableaus of the three methods are given by

$$
\begin{array}{c|cc}
\frac{1\mid 1}{\mid 1}, & \frac{1/2\mid 1/2}{\mid 1}, & \frac{0}{1}\mid 0 & 1 \\
\end{array}
\tag{3.13}
\\
\begin{array}{c|cc}
0\\\hline 1/2 \ 1/2\end{array}
\tag{3.13}
$$

from left to right for backward Euler, implicit midpoint and the implicit trapezoidal method.

The implicit midpoint method and the implicit trapezoidal method share the same stability function, given by *R*(*z*) = (2+*z*)*/*(2−*z*). The corresponding stability domain covers the entire left half-plane of the complex plane, as shown in the right panel of Figure 3.3. Both the implicit midpoint method and the trapezoidal method are therefore A-stable methods. However, since *R*(*z*) → 1 as *z* → −∞, these methods lack stiff decay and are therefore not L- stable. In general, the stability functions of implicit RK methods are rational functions, i.e., given by

$$R(z) = \frac{P(z)}{Q(z)},$$

where *P,Q* are polynomials of degree at most *s*. This is in contrast to the stability functions of explicit methods, which are always polynomials of degree at most *s*, as mentioned in Section 3.1 above.

The accuracy of the implicit methods mentioned can be easily determined using a Taylor series expansion, as outlined in Section 1.5, to confirm that the backward Euler method is a first-order accurate method, while the implicit midpoint and trapezoidial methods are both second-order accurate. It is worth noting that while explicit Runge-Kutta methods with *s* stages have an order of accuracy *p* ≤ *s*, implicit methods offer more flexibility in choosing the coefficients *aij* , potentially leading to higher accuracy for a given number of stages. In fact, the maximum order achievable for an implicit RK method is *p* = 2*s*, which is the case for the implicit midpoint method with *s* = 1 and *p* = 2. More advanced implicit RK methods will be explored later, but first, let us examine the implementation of the methods introduced so far.

### **3.3 Implementing Implicit Runge-Kutta Methods**

In the previous section, we highlighted the superior stability and generally higher accuracy of the implicit methods. Following this discussion, one might question why IRK solvers are not the default choice for all ODE problems. The answer lies in the fact that they are implicit, so the stage derivatives are defined in terms of nonlinear equations rather than explicit formulae. This fact complicates the implementation of the methods and increases the computational cost of each time step. Due to this computational overhead, explicit solvers are usually more efficient for non-stiff problems, while implicit solvers are primarily suited for stiff ODEs.

For scalar ODEs, solving equations such as (3.5) or (3.8) using Newton's method is usually not overly challenging. However, when dealing with systems of ODEs, the task becomes more complex as we need to solve a system of coupled nonlinear equations. Applying Newton's method to a system of equations involves solving a system of linear equations at each iteration, which opens up a range of possible methods for solving these systems, as well as other solver choices and parameters that can be tuned to optimize the performance. In this text, our focus is primarily on understanding the fundamental concepts of IRK solvers, and we will not delve into the details of performance optimization. Therefore, we will base our implementation on built-in equation solvers from SciPy. We will start with the backward Euler method, which is the simplest implicit method, but we will keep the implementation sufficiently general to be easily extendable to more advanced implicit methods. For a more detailed discussion on solver optimization and choices to enhance computational performance, interested readers can refer to references [1, 9].

When examining the ODESolver class introduced in Chapter 2, we can observe that many administrative tasks involved in RK methods are the same for both implicit and explicit methods. Specifically, the initialization of solution arrays and the for-loop that advances the solution remain unchanged. However, advancing the solution from one step to the next differs significantly. Therefore, it is convenient to implement implicit solvers within the existing class hierarchy and let the ODESolver superclass handle the tasks of initializing the solver and the main solver loop. The different explicit methods introduced in Chapter 2 were realized through different implementations of the advance method. We can use the same approach for implicit methods, but since each step in implicit methods involves a few more operations it is useful to introduce a couple of additional methods. For instance, a concise implementation of the backward Euler method could appear as follows:

```
from ODESolver import *
from scipy.optimize import root
class BackwardEuler(ODESolver):
    def stage_eq(self, k):
        u, f, n, t = self.u, self.f, self.n, self.t
        dt = self.dt
        return k - f(t[n] + dt, u[n] + dt * k)
    def solve_stage(self):
        u, f, n, t = self.u, self.f, self.n, self.t
        k0 = f(t[n], u[n])
        sol = root(self.stage_eq, k0)
        return sol.x
    def advance(self):
        u, f, n, t = self.u, self.f, self.n, self.t
        dt = self.dt
        k1 = self.solve_stage()
        return u[n] + dt * k1
```
Compared with the explicit solvers presented in Chapter 2, we have introduced two additional methods in our BackwardEuler class. The first method, stage\_eq(self, k), is a Python implementation of (3.6), which defines the nonlinear equation for the stage derivative. The method takes the stage derivative k as input and returns the residual of (3.6). This formulation allows us to use SciPy's nonlinear equation solvers effectively. The actual solution of the stage derivative equation is handled in the solve\_stage method. This method first computes an initial guess k0 for the stage derivative, and then passes this guess and the function stage\_eq to SciPy's root function to solve the equation. The root function is a general tool for solving nonlinear equations of the form *g*(*x*) = 0, and we apply it to solve the stage equation *k*1 − *f*(*tn* + *∆t,un* + *∆tk*1) = 0. The function returns an object of the OptimizeResult class, which includes the solution as an attribute x, along with numerous other attributes containing information about the solution process. For further details on the OptimizeResult and the root function, we refer to the SciPy documentation.

**Fig. 3.4** Solutions of the Van der Pol model for *µ* = 10, using the forward and backward Euler methods with *∆t* = 0*.*04.

We can demonstrate the superior stability of the BE method by revisiting the Van der Pol equation discussed earlier. Setting, for instance, *µ* = 10, and solving the model using the FE and BE methods gives the plots shown in Figure 3.4. The top panel shows a reference solution computed with the SciPy solve\_ivp solver using very low tolerance (rtol=1e-10). The middle panel shows the solution produced by FE with *∆t* = 0*.*04, showing visible oscillations in one of the solution components. Attempting to increase the time step further for this method leads to a divergent solution. On the other hand, the lower panel presents the solution obtained with BE, which is substantially more stable but still deviates from the reference solution in the top panel. With the BE method, increasing the time step further will still yield a stable solution, although it will deviate further from the exact solution. This simple experiment highlights the importance of considering both accuracy and stability when solving challenging systems of ODEs.2

Just as we did for the explicit methods in Chapter 2, it is possible to reuse code from the BackwardEuler class to implement other solvers. To facilitate extensive code reuse for a large range of implicit solvers, a slight modification of the code is required, which will be discussed in the next section. However, it is worth noting that a simple solver like the Crank-Nicolson method can be implemented with minimal changes to the BackwardEuler class. The class implementation could resemble the following:

```
class CrankNicolson(BackwardEuler):
    def advance(self):
        u, f, n, t = self.u, self.f, self.n, self.t
        dt = self.dt
        k1 = f(t[n], u[n])
        k2 = self.solve_stage()
        return u[n] + dt / 2 * (k1 + k2)
```
In this implementation, we leverage the fact that the stage *k*1 in the Crank-Nicolson is explicit and does not require solving an equation. On the other hand, while the definition of *k*2 is identical to that of *k*1 in the backward Euler method. Consequently, We can directly reuse both the stage\_eq and solve\_stage methods, with only the advance method needing to be reimplemented. While this compact implementation of the Crank-Nicolson method allows for code reuse, it can be argued that it violates a common principle of object-oriented programming. Subclassing and inheritance represent an "is-a" relationship, implying that an instance of the Crank-Nicolson class is also an instance of the BackwardEuler class. While this works fine in the program, and is convenient for code reuse, it is not a correct representation of the relationship between the two numerical methods. Both methods belong to the group of implicit RK solvers, but the Crank-Nicolson method is not a special case of BackwardEuler. In the following sections, we will introduce an alternative class hierarchy that reflects this relationship and enables a compact implementation of RK methods using the general formulation in (2.8)-(2.9).

### **3.4 Implicit Methods of Higher Order**

Similar to the ERK methods discussed in Chapter 2, the accuracy of IRK methods can be enhanced by increasing the number of stages. However, for implicit methods, we have more flexibility in selecting the parameters *aij* ,

2Note that the accompanying source code includes a script to generate Figure 3.4, as well as many other figures featured in the book. It is recommended to run these scripts independently and experiment with different time steps and parameters. Doing so will provide a better understanding of how the solvers work.

and this choice affects both the accuracy and computational complexity of the methods. In this section, we will explore two main branches of IRK methods: *fully implicit* methods and *diagonally implicit* methods. Both classes of methods are widely used and both have their advantages and drawbacks.

### **3.4.1 Fully Implicit RK Methods**

The most general form of RK methods is known as fully implicit methods or FIRK methods. These solvers are defined by (2.8)-(2.9), with all coefficients *aij* (potentially) non-zero. When a method has more than one stage, this formulation implies that all stage derivatives depend on all other stage derivatives, so we need to determine them all at once by solving a single system of nonlinear equations. This operation is quite expensive, but the reward is that the FIRK methods have superior stability and accuracy for a given number of stages. A FIRK method with *s* stages can have order at most 2*s*, which was the case for the implicit midpoint method in (3.8)-(3.9).

Many popular FIRK methods are based on combining standard numerical integration quadrature methods with the idea of *collocation*. Here, we present a brief overview of the derivation to illustrate the foundation shared by many important methods. For a comprehensive understanding, we recommend referring to references such as [9]. Recall from Chapter 2 that all RK methods can be viewed as approximations of equation (2.1), where the integral is approximated by a weighted sum. We set

$$u(t\_{n+1}) = u(t\_n) + \int\_{t\_n}^{t\_{n+1}} f(t, u(t)) \approx u(t\_n) + \sum\_{i=1}^{s} b\_i k\_i,\tag{3.14}$$

where *bi* are the weights and *ki* are the stage derivatives, which could be interpreted as approximations of the right-hand side function *f*(*t,u*) at distinct time points *tn* + *∆tci* .

Numerical integration is a well-established field in numerical analysis, and it is natural to choose the integration points *ci* and weights *bi* in (3.14) based on standard quadrature rules with known properties. Such quadrature rules are often derived by approximating the integrand with a polynomial which interpolates the function *f* at distinct points, and then integrating the polynomial exactly. A similar approach can be employed in deriving implicit RK methods. We approximate the solution *u* on the interval *tn < t* ≤ *tn*+1 using a polynomial *P*(*t*) of degree up to *s*, and require that *P*(*t*) satisfies the ODE exactly at distinct points *tn* +*ci∆t*. This requirement, expressed as

$$P'(t\_i) = f(t\_i, P(t\_i)), \quad t\_i = t\_n + c\_i \Delta t, i = 1, \dots, s. \tag{3.15}$$

is known as collocation, and is a widely used concept in numerical analysis. It can be shown that, given a choice of quadrature points *ci* , the collocation equations (3.15) uniquely determine the remaining coefficients *aij* and *bi* of the method, see [1] for details.

A convenient approach to deriving FIRK methods is to choose a set of collocation points *ci* , typically chosen from standard quadrature rules, and solve (3.15) to determine the remaining parameters. This strategy has led to families of FIRK methods based on common numerical integration rules. For instance, choosing *ci* as Gauss points gives rise to the Gauss methods, which are the most accurate methods for a given number of stages, achieving order 2*s*. The single-stage Gauss method corresponds to the implicit midpoint method introduced earlier, while the fourth-order Gauss method with *s* = 2 is defined by the Butcher tableau

The Gauss methods are A-stable but not L-stable. Since FIRK methods are primarily used for challenging stiff problems where stability is crucial, another family of FIRK methods, known as Radau IIA methods, is more commonly employed. These methods are based on Radau quadrature points, which include the right end of the integration interval (i.e., *cs* = 1). The onestage Radau IIA method is the backward Euler method, while the two- and three-stage versions are given by

$$\begin{array}{c|cccccc} 1/3 & 5/12 - 1/12 & & & & & & & & & \\ \hline 1 & 3/4 & 1/4 & & & & & & & & & & \\ \hline 1 & 3/4 & 1/4 & & & & & & & & & & \\ \hline 1 & 2/3 & 1/4 & & & & & & & & & & \\ \hline & 2/3 & 1/4 & & & & & & & & & & \\ \hline & & & & & & & & & & & & & & \\ \hline & & & & & & & & & & & & & \\ \hline & & & & & & & & & & & & & \\ \hline & & & & & & & & & & & & & \\ \hline & & & & & & & & & & & & \\ \hline & & & & & & & & & & & & \\ \hline & & & & & & & & & & & & \\ \hline & & & & & & & & & & & & \\ \hline & & & & & & & & & & & & \\ \hline & & & & & & & & & & & & \\ \hline & & & & & & & & & & & & \\ \hline & & & & & & & & & & & & \\ \hline & & & & & & & & & & & & \\ \hline & & & & & & & & & & & \\ \hline & & & & & & & & & & & \\ \end{array}$$

The Radau IIA methods exhibit order 2*s* − 1, and their stability functions are (*s*−1*,s*) Padé approximations of the exponential function, as described in [9]. For the two- and three-stage methods mentioned earlier, the stability functions are given by

$$\begin{aligned} R(z) &= \frac{1+z/3}{1-2z/3+z^2/6}, \\ R(z) &= \frac{1+2z/5+z^2/20}{1-3z/5+3z^2/20-z^2/60}, \end{aligned}$$

respectively. The stability domains of these methods are depicted in Figure 3.5. Due to their L-stability, the Radau IIA methods are commonly used for solving stiff ODE systems. However, as noted above, the fact that all *aij* 6= 0 complicates the implementation of the methods and makes each time step computationally expensive. All the *s* equations of (2.8) become fully coupled and need to be solved simultaneously. In the case of an ODE system comprising *m* equations, we must solve a system of *ms* nonlinear equations for each time step. We will come back the implementation of FIRK methods in Section 3.5, but let us first introduce a slightly simpler class of implicit RK solvers.

**Fig. 3.5** The shaded area represents the stability region for two of the RadauIIA methods, with *s* = 2 (left) and *s* = 3 (right).

### **3.4.2 Diagonally Implicit RK Methods**

Diagonally implicit RK (DIRK) methods, also known as semi-explicit methods, are a subclass of implicit RK methods. For DIRK methods, we have *aij* = 0 for all *j > i*. (Notice the small but important difference from the explicit methods, where we have *aij* = 0 for *j* ≥ *i*.) The consequence of this choice is that the equation for a single stage derivative *ki* does not involve stages *ki*+1*, ki*+2, and so on, and we can sequentially solve for the stage derivatives one by one. We still need to solve nonlinear equations to determine each *ki* , but we can solve *s* systems of *m* equations rather than solving one large system to compute all stages simultaneously. This simplifies the implementation and reduces the computational cost per time step. However, the restriction on the method coefficients also reduces the accuracy and stability compared with FIRK methods. A general DIRK method with *s* stages has a maximum order of *s*+ 1, and methods optimized for stability typically have even lower order.

It is worth nothing that the implicit midpoint method discussed earlier technically falls under the category of DIRK methods. However, it is also a fully implicit Gauss method, and is not commonly referred to as a DIRK method. The distinction between FIRK and DIRK methods is meaningful only for *s >* 1. The Crank-Nicolson (implicit trapezoidal) method given by (3.10)-(3.12) is another example of a DIRK method, evident from the rightmost Butcher tableau in (3.13). These methods are, however, only A-stable, and it is possible to derive DIRK methods with better stability properties. An example of an L-stable, two-stage DIRK method of order two is given by

$$\frac{\begin{array}{c} \gamma\\1 \end{array} \Big| \begin{array}{c} \gamma \quad 0\\1 - \gamma \ \gamma \end{array}}{\begin{array}{c} 1 - \gamma \ \gamma \end{array}},\tag{3.16}$$

with stability function

$$R(z) = \frac{1 + z(1 - 2\gamma)}{(1 - z\gamma)^2}.$$

This method is A-stable for *γ >* 1*/*4, and for *γ* = 1 ± √ 2*/*2 the method is L-stable and second-order accurate. Note that choosing *γ >* 1 means that we estimate the stage derivatives outside the interval (*tn, tn*+1), and for the last step beyond the time interval of interest. While this does not affect the stability or accuracy of the method, it may not be suitable for all ODE problems, and the most popular choice is therefore *γ* = 1− √ 2*/*2 (≈ 0*.*293). Notice also that in this method the two diagonal entries of *aij* are identical, with *a*11 = *a*22 = *γ*. This choice is very common in DIRK methods, and methods of this kind are known as *singly diagonally implicit* RK (SDIRK) methods. The main benefit of this structure is that the nonlinear equations for each stage derivative become very similar, which can be utilized when solving the equations using quasi-Newton methods. This benefit may not be very obvious for the examples in this book, since we rely on the generic root function from scipy.optimize to solve the nonlinear equations. However, if we wanted to improve the computational performance of the solvers, a natural place to start would be to implement a custom quasi-Newton solver that exploits the particular structure of the nonlinear equations. We will not go into the details of such an implementation here, but it is worth commenting on some aspects in order to understand why SDIRK methods are so widely used.

The central point is that when applying Newton's method to solve a general nonlinear system *g*(*u*) = 0, each iteration involves solving linear systems of the form *Jg∆u* = −*g*(*u k* ), where *∆u* is the increment to the solution, *u k* is the solution value at the previous iteration, and *Jg* is the Jacobian matrix of *g*, defined by

$$J\_g = \frac{\partial g\_i}{\partial u\_j}.$$

For a general DIRK method, the nonlinear equation to compute stage derivative *ki* is given by

$$k\_i = f(t\_n + c\_i \Delta t, u\_n + \Delta \sum\_{j=1}^i a\_{ij} k\_j),$$

which can be written in the form *g*(*ki*) = 0, with

$$g(k\_i) = k\_i - f(t\_n + c\_i \Delta t, u\_n + \Delta \left(\sum\_{j=1}^{i-1} a\_{ij} k\_j + a\_{ii} k\_i\right).$$

Note that we have split the sum over the stage derivatives, highlighting that when solving for *ki* , the values *kj* for *j < i* are already known. The Jacobian matrix *Jg* is found by differentiating *g* with respect to *ki* , resulting in

$$J\_g = I - \Delta t a\_{ii} J\_f,$$

where *Jf* is the Jacobian of the right-hand side function *f*. If we have *aii* = *γ* for all stages, the Jacobian matrices *Jg* will also be identical, which can be exploited to optimize the solution of the linear systems. For more detailed information on solving nonlinear equations arising in SDIRK methods, refer to references such as [9].

Although we do not aim to present a complete overview of all subcategories of RK methods, one additional method class is worth mentioning. These are the so called ESDIRK (explicit singly diagonally implicit RK) methods, which are simply SDIRK methods where the first stage is explicit. The motivation behind these methods is that the nonlinear algebraic equations involved in the implicit methods are always solved with iterative methods, requiring an initial guess for the solution. For SDIRK methods, it is convenient to use the previous stage derivate as initial guess for the next one, which will usually provide a good initial guess. This approach is obviously not possible for the first stage, but an explicit formula for the first stage solves this problem. The simplest ESDIRK method is the implicit trapezoidal (Crank-Nicolson) method introduced above. A popular extension of this method is given by the following Butcher tableau:

$$\begin{array}{c} 0 \\ 2\gamma \Big| \gamma \text{ } 0 \\ 1 \Big| \beta \text{ } \beta \text{ } \gamma \text{ } \\ \hline \beta \text{ } \beta \text{ } \gamma \end{array} \tag{3.17}$$

with *γ* = 1− √ 2*/*2 and *β* = √ 2*/*4. The resulting equations for each time step are

$$\begin{aligned} k\_1 &= f(t\_n, u\_n), \\ k\_2 &= f(t\_n + 2\gamma \Delta t, u\_n + \Delta t(\gamma k\_1 + \gamma k\_2)), \\ k\_3 &= f(t\_n + \Delta t, u\_n + \Delta t(\beta k\_1 + \beta k\_2 + \gamma k\_3)), \\ u\_{n+1} &= u\_n + \Delta t(\beta k\_1 + \beta k\_2 + \gamma k\_3). \end{aligned}$$

This method can be interpreted as the sequential application of the trapezoidal method and a popular multistep solver called BDF2 (*backward differentiation formula* of order 2), and it is commonly known as the TR-BDF2 method. It is second order accurate, like the trapezoidal rule, but it is also L-stable, making it suitable for stiff problems.

### **3.5 Implementing Higher Order IRK Methods**

In Section 3.3 we implemented two of the simplest implicit RK methods by a relatively small extension of the ODEsolver class hierarchy. We could easily continue this idea for the more complex IRK methods, and all the different methods could be realized by separate implementations of the three methods solve\_stage, stage\_eq, and advance. However, these three methods essentially implement the equations given by (2.8)-(2.9), which are common for all RK solvers. It is natural to look for an implementation that allows even more code reuse between the various methods, and we shall see that such a general implementation is indeed possible. However, it is still useful to treat the fully implicit methods and SDIRK methods separately, since the stage calculations of these two method classes are fundamentally different.

### **3.5.1 A Base Class for Fully Implicit Methods**

One approach to implementing the fully implicit RK methods is to rewrite the solve\_stage, stage\_eq, and advance methods of the BackwardEuler class in a more general manner that can handle any number of stages and method parameters *aij , bi* , and *ci* . By adopting this approach, new methods can easily be implemented by specifying the number of stages and defining the parameter values. In the methods we have discussed so far, the method coefficients have been hard-coded within the mathematical expressions, often inside the advance methods. However, with the generic approach, it is more natural to define these coefficients as class attributes in the constructor. Following this general approach, a base class for implicit RK methods can be defined as follows:

```
from ODESolver import *
from scipy.optimize import root
class ImplicitRK(ODESolver):
    def solve_stages(self):
        u, f, n, t = self.u, self.f, self.n, self.t
        s = self.stages
        k0 = f(t[n], u[n])
        k0 = np.tile(k0,s)
        sol = root(self.stage_eq, k0)
        return np.split(sol.x, s)
    def stage_eq(self, k_all):
        a, c = self.a, self.c
        s, neq = self.stages, self.neq
        u, f, n, t = self.u, self.f, self.n, self.t
        dt = self.dt
        res = np.zeros_like(k_all)
        k = np.split(k_all, s)
        for i in range(s):
            fi = f(t[n] + c[i] * dt, u[n] + dt *
                   sum([a[i, j] * k[j] for j in range(s)]))
            res[i * neq:(i + 1) * neq] = k[i] - fi
        return res
    def advance(self):
        b = self.b
        u, n, t = self.u, self.n, self.t
        dt = self.dt
        k = self.solve_stages()
        return u[n] + dt * sum(b_ * k_ for b_, k_ in zip(b, k))
```
Note that we assume that the method parameters are stored in NumPy arrays self.a, self.b, self.c, which need to be defined in subclasses. It is important to note that, just as the ODESolver class discussed earlier, the ImplicitRK class is intended as a pure base class for holding common code. It is not meant to be used as a standalone solver class. In accordance with the principles described in Section 2.2, we could make the abstract nature of this class explicit by using the abc module, but for the present text we focus on the fundamentals of the solvers and the class structure, keeping the code as simple and compact as possible.

The three methods of the ImplicitRK class are generalizations of the corresponding methods in the BackwardEuler class. They perform the same tasks but at a higher abstraction level and they rely on a bit of NumPy magic:

With the general base class in place, it becomes straightforward to implement new solvers by writing constructors that define the method coefficients. The following code implements the implicit midpoint and the two- and three-stage Radau methods:

```
class ImplicitMidpoint(ImplicitRK):
    def __init__(self, f):
        super().__init__(f)
        self.stages = 1
        self.a = np.array([[1 / 2]])
        self.c = np.array([1 / 2])
        self.b = np.array([1])
class Radau2(ImplicitRK):
    def __init__(self, f):
        super().__init__(f)
        self.stages = 2
```

```
self.a = np.array([[5 / 12, -1 / 12], [3 / 4, 1 / 4]])
        self.c = np.array([1 / 3, 1])
        self.b = np.array([3 / 4, 1 / 4])
class Radau3(ImplicitRK):
    def __init__(self, f):
        super().__init__(f)
        self.stages = 3
        sq6 = np.sqrt(6)
        self.a = np.array([[(88 - 7 * sq6) / 360,
                            (296 - 169 * sq6) / 1800,
                            (-2 + 3 * sq6) / (225)],
                           [(296 + 169 * sq6) / 1800,
                            (88 + 7 * sq6) / 360,
                            (-2 - 3 * sq6) / (225)],
                           [(16 - sq6) / 36, (16 + sq6) / 36, 1 / 9]])
        self.c = np.array([(4 - sq6) / 10, (4 + sq6) / 10, 1])
        self.b = np.array([(16 - sq6) / 36, (16 + sq6) / 36, 1 / 9])
```
Notice that we always define the method coefficients as NumPy arrays, even for the implicit midpoint method where they only contain a single number. This definition is necessary for the generic methods of the ImplicitRK class to work.

### **3.5.2 Base Classes for SDIRK and ESDIRK Methods**

We could, in principle, implement both the SDIRK and ESDIRK methods in the same manner as the FIRK methods, by defining the method coefficients in the constructor and using the generic methods from the ImplicitRK base class. These generic methods would handle the cases where *aij* = 0 for *j > i*. However, the motivation behind the diagonally implicit methods is to avoid solving large systems of nonlinear equations, so it does not make much sense to implement them in this way. Instead, we should take advantage of the specific structure of the method coefficients and solve for the stage variables sequentially. This requires rewriting solve\_stages and stage\_eq methods from the base class. However, the advance method, which advances the solution to the next step, can remain unchanged as it is common to all RK methods.

Considering first the SDIRK methods, we can implement these as subclasses of the ImplicitRK class, which allows for some moderate code reuse and reflects the fact that SDIRK methods are special cases of implicit RK methods. To illustrate this implementation, let us first consider the two-stage SDIRK method defined by (3.16), and write out the equations for the stage derivatives to get a better view of the tasks involved. Inserting the coefficients from (3.16) into (2.8)-(2.9) gives

$$k\_1 = f(t\_n + \gamma \Delta t, u\_n + \gamma \Delta t k\_1),\tag{3.18}$$

$$k\_2 = f(t\_n + \Delta t, u\_n + \Delta t((1 - \gamma)k\_1 + \gamma k\_2)),\tag{3.19}$$

$$u\_{n+1} = u\_n + \Delta t((1 - \gamma)k\_1 + \gamma k\_2). \tag{3.20}$$

Here, (3.18) is nearly identical to the equation defining the stage derivative in the backward Euler method, with the only difference being that *∆t* is replaced with *γ∆t*. Similarly, the only difference between (3.18) and (3.19) is the additional term *∆t*(1 − *γ*)*k*1 inside the function call. In general, any stage equation for any DIRK method can be written as

$$k\_i = f(t\_n + c\_i \Delta t, u\_n + \Delta t(\sum\_{j=0}^{i-1} a\_{ij} k\_j + \gamma k\_i)),\tag{3.21}$$

where the sum inside the function call only includes previously computed stages.

Given the similarity of (3.21) with the stage equation from the backward Euler method, it is natural to implement the SDIRK stage equation as a generalization of the stage\_eq method from the BackwardEuler class. To achieve this, we can create an SDIRK base class that contains the general versions of both the stage\_eq and solve\_stages methods. This base class can then be used as a foundation for deriving specific SDIRK solver classes. By writing the stage equations in this general form, it becomes straightforward to generalize the algorithm for looping through the stages and computing the individual stage derivatives. The complete base class implementation may appear as follows.

```
class SDIRK(ImplicitRK):
    def stage_eq(self, k, c_i, k_sum):
        u, f, n, t = self.u, self.f, self.n, self.t
        dt = self.dt
        gamma = self.gamma
        return k - f(t[n] + c_i * dt, u[n] + dt * (k_sum + gamma * k))
    def solve_stages(self):
        u, f, n, t = self.u, self.f, self.n, self.t
        a, c = self.a, self.c
        s = self.stages
        k = f(t[n], u[n]) # initial guess for first stage
        k_sum = np.zeros_like(k)
        k_all = []
        for i in range(s):
            k_sum = sum(a_ * k_ for a_, k_ in zip(a[i, :i], k_all))
            k = root(self.stage_eq, k, args=(c[i], k_sum)).x
            k_all.append(k)
        return k_all
```
The modified stage\_eq method takes two additional parameters: the coefficient c\_i, corresponding to the current stage, and the array k\_sum, which holds the sum P*i*−1 *j*=1 *aijkj* . These arguments need to be initialized correctly for each stage and passed as additional arguments to the SciPy root function. For convenience, we also assume that the method parameter *γ* has been stored as a separate class attribute. With the stage\_eq method implemented in this general way, the solve\_stages method simply needs to update the weighted sum of previous stages (k\_sum), and pass this and the correct *c* value as additional arguments to the SciPy root function. The implementation uses a for loop to compute the stage derivatives sequentially and returns them as a list k\_all.

As for the FIRK method classes, the only method we need to implement specifically for each solver class is the constructor, in which we define the number of stages and the method coefficients. A class implementation of the method in (3.16) may look as follows.

```
class SDIRK2(SDIRK):
    def __init__(self, f):
        super().__init__(f)