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.

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, 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)
```

Now 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
```

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.