Differential Equations in Python. Part 1: Initial Value Problems

This article starts a series of posts on Differential Equations and how to solve then in Python.

This Part 1 covers how to solve Ordinary Differential Equations (ODE), more specifically, Initial Value Problems (IVP). In following parts we will cover Boundary Value Problems (BVP) and Partial Differential Equations (PDE). You can check them out here:

Ordinary Differential Equations

An Ordinary Differential Equation (ODE) is a differential equation containing one or more functions of one independent variable and the derivatives of those functions.

An initial value problem (IVP) is an ODE together with an initial condition that specifies the value of the unknown function at a given point in the domain. It is a common problem in the modeling of physical and chemical systems where a differential equation describes how that system evolves with time for a known initial state. It can be generally stated as:

$$\frac{dy}{dt} = f(t, y(t))$$
$$t = 0: \quad y = y_0$$

Reduction of order

Differential equations containing second-order derivatives or higher order derivatives, can usually be transformed into a system of first-order ODEs. This approach greatly simplifies the process of solving the problem.

A ODE of order \(n\) can be transformed into a system of \(n\) first-order ODEs. In more general terms, an \(n\)-order ODE of the form:

$$\frac{d^ny}{dx^n} = F(x, y, \frac{dy}{dx}, \frac{d^2y}{dx^2}, ..., \frac{d^{n-1}y}{dx^{n-1}})$$

can be written as a system of \(n\) first-order ODEs by defining a new set of variables such that:

$$\frac{dy_1}{dx} = y_2$$
$$\frac{dy_2}{dx} = y_3$$
$$\frac{dy_{n-1}}{dx} = y_n$$
$$\frac{dy_n}{dx} = F(x, y_1, ..., y_n)$$

Solving an IVP in Python

The example below represents a chemical reaction occurring from an initial know condition until equilibrium is reached. It is a IVP consisting of a system of four ODEs.

The problem

The following three chemical reactions (two reversible and one irreversible) take place in liquid phase in a closed isothermal reactor.

$$A \underset{k_1'}{\stackrel{k_1}{\rightleftharpoons}} B \underset{}{\stackrel{k_2}{\rightarrow}} C \underset{k_3'}{\stackrel{k_3}{\rightleftharpoons}} D$$

All reaction steps are elementary and the initial concentrations of A, B, C, and D are 1, 0, 0, and 1 \(\mathrm{mol \ dm^{-3}}\), respectively. The kinetic constants are: \(k_1 = 1 \ \mathrm{min^{-1}}\), \(k_1' = 2 \ \mathrm{min^{-1}}\), \(k_2 = 1 \ \mathrm{min^{-1}}\), \(k_3 = 1 \ \mathrm{min^{-1}}\), and \(k_3' = 2 \ \mathrm{min^{-1}}\).

We want to calculate the concentrations of all species over time.

From the mass balance to all species we reach the following system of ODEs, in time, that represent our system:

$$\frac{dC_A}{dt} = -k_1 C_A + k_2 C_B $$
$$\frac{dC_B}{dt} = k_1 C_A - k_1' C_B - k_2 C_B$$
$$\frac{dC_C}{dt} = k_2 C_B - k_3 C_C + k_3' C_D$$
$$\frac{dC_D}{dt} = k_3 C_C - k_3' C_D$$

To solve this system of four ODEs we require four initial conditions, which are provided:

$$t = 0 \ \mathrm{min}: \quad C_A = 1, \quad C_B = 0, \quad C_C = 0, \quad C_D = 1 \quad \mathrm{mol \ dm^{-3}}$$


Python's scipy package has all necessary tools to solve this problem. We start by importing the necessary packages.

# Imports
import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt

Next, we input the know process data and set a time limit for integration. This is the time until which we want to obtain a solution. As in most actual processes the initial time corresponds to 0.

# Data
k1 = k2 = k3 = 1  # min-1, kinetic constants
k1_ = k3_ = 2  # min-1, kinetic constants
CA0 = CD0 = 1  # mol/dm3, initial concentrations
CB0 = CC0 = 0  # mol/dm3, initial concentrations
tfinal = 20  # min, final integration time

We need to create a function (here called sode) where we define the system of differential equations we want to solve. We will later pass this function to the scipy algorithm that will perform the integration (it is called solve_ivp). As such, we must provide our equations in the format expected by solve_ivp.

The function provided to solve_ivp must return a list of expressions where each list element represents the right hand side of a differential equation when placed in the format:

$$\frac{dC}{dt} = f(t, C(t))$$

The independent (\(t\)) and dependent (C) variables are arguments of the sode function.

# Define the system of ordinary differential equations (sode)
def sode(t, y):
    return [
        -k1*y[0] + k1_*y[1],
        k1*y[0] - k1_*y[1] - k2*y[1],
        k2*y[1] - k3*y[2] + k3_*y[3],
        k3*y[2] - k3_*y[3]

The hardest part is over. To solve this system of equations we simply call the scipy.integrate.solve_ivp algorithm. We must provide the function that defines our problem (sode), a vector containing the range of times we want a solution for (here called tspan), and a vector containing our initial conditions (y0).

# Solve the sedo using scipy.integrate.solve_ivp
tspan = [0, tfinal]  # integration interval
y0 = [CA0, CB0, CC0, CD0]  # initial state (initial concentrations)
sol = integrate.solve_ivp(sode, tspan, y0, dense_output=True)

We can obtain our final solution for specific time as follows.

t = np.linspace(0, tfinal, 1000)
C = sol.sol(t)

Finally, we plot the concentrations over time.

# Display results graphically
fig1, ax1 = plt.subplots()
ax1.plot(t, C.T)
plt.xlabel('$t$ ($\mathrm{min}$)')
plt.ylabel('$C_i$ ($\mathrm{mol \ dm^{-3}}$)')
plt.legend(['$C_A$', '$C_B$', '$C_C$', '$C_D$']);


Plotting the results of concentrations of each species over time we get.


As expected, for \(t = 0 \ \mathrm{min}\) we have only A and D. B is initially produced but immediately converts into C. As time progresses, A and B are completely consumed to produce C and D until these two reach an equilibrium.

Source code

You can find the complete source code of this post here.