Solving systems of differential equations

Differential equations sometimes occur in systems consisting of two or more interlinked differential equations. A classical example is a simple model of the populations of competing species. This is a simple model of competing species labeled P (the prey)and W (the predators) given by the following equations:

The first equation dictates the growth of the prey species P, which, without any predators, would be exponential growth. The second equation dictates the growth of the predator species W, which, without any prey, would be exponential decay. Of course, these two equations are coupled; each population change depends on both populations. The predators consume the prey at a rate proportional to the product of their two populations, and the predators grow at a rate proportional to the relative abundance of prey (again the product of the two populations).

In this recipe, we will will analyze a simple system of differential equations and use the SciPy integrate module to obtain approximate solutions.

Getting ready

The tools for solving a system of differential equations using Python are the same as those for solving a single equation. We again use the solve_ivp routine from the integrate module in SciPy. However, this will only give us a predicted evolution over time with given starting populations. For this reason, we will also employ some plotting tools from Matplotlib to better understand the evolution.

How to do it...

The following steps walk through how to analyze a simple system of differential equations:

1. Our first task is to define a function that holds the system of equations. This function needs to take two arguments as for a single equation, except the dependent variable y (in the notation from the Solving simple differential equations numerically recipe) will now be an array with as many elements as there are equations. Here, there will be two elements. The function we need for the example system in this recipe is as follows:

def predator_prey_system(t, y):
    return np.array([5*y[0] - 0.1*y[0]*y[1], 0.1*y[1]*y[0] -
       6*y[1]])

1. Now we have defined the system in Python, we can use the quiver routine from Matplotlib to produce a plot that will describe how the populations will evolve—given by the equations—at numerous starting populations. We first set up a grid of points on which we will plot this evolution. It is a good idea to choose a relatively small number of points for the quiver routine, otherwise it becomes difficult to see details in the plot. For this example, we plot the population values between 0 and 100:

p = np.linspace(0, 100, 25)
w = np.linspace(0, 100, 25)
P, W = np.meshgrid(p, w)

1. Now, we compute the values of the system at each of these pairs. Notice that neither equation in the system is time-dependent (they are autonomous); the time variable t is unimportant in the calculation. We supply the value 0 for the t argument:

dp, dw = predator_prey_system(0, np.array([P, W]))

1. The variables dp and dw now hold the "direction" in which the population of P and W will evolve, respectively, if we started at each point in our grid. We can plot these directions together using the quiver routine from matplotlib.pyplot:

fig, ax = plt.subplots()
ax.quiver(P, W, dp, dw)
ax.set_title("Population dynamics for two competing species")
ax.set_xlabel("P")
ax.set_ylabel("W")

Plotting the result of these commands now gives us Figure 3.2, which gives a "global" picture of how solutions evolve:

To understand a solution more specifically, we need some initial conditions so we can use the solve_ivp routine described in the previous recipe.

1. Since we have two equations, our initial conditions will have two values. (Recall in the Solving simple differential equations numerically recipe, we saw that the initial condition provided to solve_ivp needs to be a NumPy array.) Let's consider the the initial values P(0) = 85 and W(0) = 40. We define these in a NumPy array, being careful to place them in the correct order:

initial_conditions = np.array([85, 40])

1. Now we can use solve_ivp from the scipy.integrate module. We need to provide the max_step keyword argument to make sure that we have enough points in the solution to give a smooth solution curve:

from scipy import integrate
sol = integrate.solve_ivp(predator_prey_system, (0., 5.),
   initial_conditions, max_step=0.01)

1. Let's plot this solution on our existing figure to show how this specific solution relates to the direction plot we have already produced. We also plot the initial condition at the same time:

ax.plot(initial_conditions[0], initial_conditions[1], "ko")
ax.plot(sol.y[0, :], sol.y[1, :], "k", linewidth=0.5)

The result of this is shown in Figure 3.3:

How it works...

The method used for a system of ordinary differential equations is exactly the same as for a single ordinary differential equation. We start by writing the system of equations as a single vector differential equation,

that can then be solved using a time-stepping method as though y were a simple scalar value.

The technique of plotting the directional arrows on a plane using the quiver routine is a quick and easy way of learning how a system might evolve from a given state. The derivative of a function represents the gradient of the curve (x, u(x)), and so a differential equation describes the gradient of the solution function at position y and time t. A system of equations describes the gradient of separate solution functions at a given position y and time t. Of course, the position is now a two-dimensional point, so when we plot the gradient at a point, we represent this as an arrow that starts at the point, in the direction of the gradient. The length of the arrow represents the size of the gradient; the longer the arrow, the "faster" the solution curve will move in that direction.

When we plot the solution trajectory on top of this direction field, we can see that the curve (starting at the point) follows the direction indicated by the arrows. The behavior shown by the solution trajectory is a limit cycle, where the solution for each variable is periodic as the two species populations grow or decline. This description of the behavior is perhaps more clear if we plot each population against time, as seen in Figure 3.4. What is not immediately obvious from Figure 3.3 is that the solution trajectory loops around several times, but this is clearly shown in Figure 3.4:

Figure 3.4: Plots of populations P and W against time. Both populations exhibit periodic behavior

There's more...

The technique of analyzing a system of ordinary differential equations by plotting the variables against one another, starting at various initial conditions, is called phase space (plane) analysis. In this recipe, we used the quiver plotting routine to quickly generate an approximation of the phase plane for the system of differential equations. By analyzing the phase plane of a system of differential equations, we can identify different local and global characteristics of the solution, such as limit cycles.