dan hammer

              

About that which is me

Projects that I work on

Academics that I am part of

Posts that I wrote

Life that I live

Press that mentions

Differential equations in Python


There is often no analytical solution to systems with nonlinear, interacting dynamics. We can, however, examine the dynamics using numerical methods. Consider the predator-prey system of equations, where there are fish () and fishing boats ():

We use the built-in SciPy function odeint to solve the system of ordinary differential equations, which relies on lsoda from the FORTRAN library odepack. First, we define a callable function to compute the time derivatives for a given state, indexed by the time period. We also load libraries that we’ll use later to animate the results.

import matplotlib.animation as animation
from scipy.integrate import odeint
from numpy import arange
from pylab import *

def BoatFishSystem(state, t):
    fish, boat = state
    d_fish = fish * (2 - boat - fish)
    d_boat = -boat * (1 - 1.5 * fish)
    return [d_fish, d_boat]

Then, we define the state-space and intital conditions, so that we can solve the system of linear equations. The result is animated below. (The code for some of the graphical bells and whistles is omitted for the sake of exposition.)

t = arange(0, 20, 0.1)
init_state = [1, 1]
state = odeint(BoatFishSystem, init_state, t)

fig = figure()
xlabel('number of fish')
ylabel('number of boats')
plot(state[:, 0], state[:, 1], 'b-', alpha=0.2)

def animate(i):
    plot(state[0:i, 0], state[0:i, 1], 'b-')

ani = animation.FuncAnimation(fig, animate, interval=1)
show()

The red, dashed lines indicate the nullclines, derived from the first-order conditions of the equation system. These lines delineate the phase space of the top graph; and the lines intersect at the equilibrium levels of fish and boats.

It is easy to break this result by messing with the solver parameters or the size of the time steps (relative to the total time), demonstrating the fragility of the result for real-world applications. If, for example, we increase the step size from 0.1 to 5, we lose most of the dynamics that characterize the system. The same goes for fiddling with the iteration parameters of the ODE solver.

Suppose we wanted to figure out the behavior of this system near the equilibrium before going through the numerical estimation. First, we linearize the system near the equilibrium, yielding the Jacobian. Let and , then the nonlinear system can be approximated by the following linear system near the equilibrim :

Noting that by definition of the equilibrium, the nonlinear system is approximated by the linear system defined by the Jacobian, evaluated at the equilibrium:

Then the eigenvalues are given by:

The real part is negative and there is an imaginary component, such that the system will oscillate around the equilibrium tending inward. This behavior is reflected in the animation. I guess we didn’t really have to go through all this work; but whatever, it’s useful for other problems.


Schaefer model

Bjorndal and Conrad (1987) modelled open-access exploitation of North Sea herring between 1963 - 1977. Their model is similar to the one above, except slightly more complicated. Let fish stock () and fishing effort () be modelled by the following system:

where is a catchability constant, is the intrinsic growth rate of the fish stock, is the carrying capacity, is the fish price, and is the marginal cost of one unit of effort. Then, through the same process as above, we find that the equilibrium point (at the intersection of the nullclines) is:

Using the constants in Bjorndal and Conrad (1987) we model the system similarly:

price = 735
effort_scale = 2e-6
marginal_cost = 556380
carrying_capacity = 3.2e6
intrinsic_growth = 0.08
catchability = marginal_cost / (price * 0.25e6)

def BoatFishSystem(state, t, time_scale=0.1):
    stock, effort = state
    net_growth = intrinsic_growth * stock * (1 - (stock/carrying_capacity))
    harvest = catchability * stock * effort
    d_stock = net_growth - harvest
    d_effort = effort_scale * (price * catchability * stock - marginal_cost)
    return [d_stock * time_scale, d_effort * time_scale]

The Jacobian for this system evaluated at the equilibrium is:

I don’t want to solve this by hand, so I plug it into Python.

from numpy import linalg
values, vectors = linalg.eig(J_equil)
print values
>>> [-0.003125+41.01820462j -0.003125-41.01820462j]

The eigenvalues are . Once again, the behavior seen in the numerical approximation is confirmed by math. The system oscillates and tends inward.

Shit can get weird. The numerical approximations, while very good in Python, can misrepresent the system of equations, given certain parameters. Specifically, the system is solved through an iterative process of calculating the linear change at each interval, approximating the continuous system. Choosing certain step sizes and tolerances will send Python or Matlab into a tailspin. Although, the checks and balances within odeint are really quite good, such that it’s way easier to break the numerical approximation if you try to write it explicitly in a for-loop.