Diffential Equations with Julia, updated.

The new release of the DifferentialEquations package introduces a number of breaking changes, which are all for the best. So I will revisit my [previous] post on this package, to show some of the new syntax in action.

[previous]: {% link _posts/2017-11-06-solving-differential-equations-for-ecological-problems.md %}

In a few words, the new release defines parameters of the differential equations system as part of the problem, so they can be passed directly to the function. The new way to declare a problem to solve is f(u, p, t), where u is the quantity to model, p are the parameters, and t is the time.

So let’s start:

# Logistic model with r=1.1, K=10.0
f(u,p,t) = u*1.1*(1.0-u/10.0)
# Now we solve the problem
problem = ODEProblem(f, 0.5, (0.0, 10.0))
solution = solve(problem)

Now, we want to modify the values of $r$ and $K$ – this can be done using the p parameter. First, let’s re-write f so that it calls values within p:

f(u,p,t) = u*p[1]*(1.0-u/p[2])

The next step is to provide ODESolver with an array containing the parameters:

problem = ODEProblem(f, 0.5, (0.0, 10.0), [1.4, 8.0])
solution = solve(problem)

With these ingredients, we have anything to reproduce the Allee effect model from the [previous post][previous], in a way which is much more concise:

allee(u,p,t) = u*p[1]*(u/p[3]-1.0)*(1.0-u/p[2])

This takes three parameters, $r$, $K$, and $A$. This allows calling the function directly without having to re-define the problem multiple times. As a consequence, the code is easier to read, to write, and to maintain. Of course, we might also want to use a dictionary to store the parameters. It works out of the box too.

p = Dict(:r => 1.1, :K => 3.0, :A => 1.0)
allee(u,p,t) = u*p[:r]*(u/p[:A]-1.0)*(1.0-u/p[:K])
pr = ODEProblem(allee, 1.2, (0.0, 10.0), p)
solve(pr)

Now, it’s time for a short example. Let’s write-up the logistic growth model with competition - we’ll use the version with a pre-allocated variable for the derivative:

p = Dict(
  :r1 => 1.1, :r2 => 1.4,
  :a11 => 1.0, :a22 => 1.0,
  :a21 => 1.2, :a12 => 1.1
)

function competition(du, u,p,t)
  du[1] = u[1]*p[:r1] - p[:a11]*u[1]^2 - p[:a12]*u[1]*u[2]
  du[2] = u[2]*p[:r2] - p[:a22]*u[2]^2 - p[:a21]*u[2]*u[1]
end

du = zeros(2)
u0 = [0.1, 0.1]

pro = ODEProblem(competition, u0, (0.0, 10.0), p)
sol = solve(pro)

Output of the initial runNow we can run this model, and as expected given the parameters values, the first species will gradually go extinct. The nice thing with this interface is that we don’t need to re-declare the function describing the differential equations to use different parameters: this can be done when declaring the ODEProblem to solve.

As a consequence, the real power of the new interface is that iterations over large collections of parameters becomes trivial. As an example, we can look at transgressive overyielding, which is defined as the coexistence equilibrium having a higher density than the best single population equilibrium, when changing the values of $\alpha_{12}$ and $\alpha_{21}$. In principle we don’t need to run the model for the two species alone (because we vary parameters that are not involved in the single species equilibrium), but computing time is cheap.

N = 150
output = zeros(N, N)
alpha = collect(linspace(0.2, 2.0, N))
for a1 in eachindex(alpha), a2 in eachindex(alpha)
    p[:a12] = alpha[a1]
    p[:a21] = alpha[a2]
    coexist = last(solve(ODEProblem(competition, [0.1, 0.1], (0.0, 100.0), p)))
    s1 = last(solve(ODEProblem(competition, [0.1, 0.0], (0.0, 100.0), p)))[1]
    s2 = last(solve(ODEProblem(competition, [0.0, 0.1], (0.0, 100.0), p)))[2]
    output[a1,a2] = sum(coexist)/(max(s1,s2))
end

Heatmap

And that’s it! We see more transgressive overyielding when both interspecific competition rates are low relative to intraspecific competition.

The DifferentialEquations package is amazing, and the code is both easy to write and fast to run. You should definitely check out some of the other outstanding functions it offers.