Numerical integration, the search for solutions of differential equations, is a hallmark of scientific computing. In this lesson, we will see how we can apply multipe concepts to write our own routine for the second-order Runge-Kutta method. In practice, it is never recommended to write one’s own routine for numerical integration, as there are specific packages to handle this task. In Julia, this is DifferentialEquations.jl. This being said, writing a Runge-Kutta method is an interesting exercise.

Glossing over a few details, numerical integration requires a function $f$, giving the derivative at time $t$, a timestep $h$ at which to evaluate $f$, and an initial quantity $y$. There are multiple notations, but for this lesson we will adopt $y(t)$ for the value of $y$ at time $t$. The second-order Runge-Kutta method (RK2) works by evaluating the derivative at time $t+h/2$, and using this information to correct the initial assessment. This helps (very moderately) in problems where $y$ varies a lot over short periods of time.

In practice, RK2 is described by the following system of equations:

$$y(t+\frac{h}{2}) = y(t)+\frac{h}{2}f\left(t, y(t)\right)$$

$$y’(t+\frac{h}{2}) = f\left(t+\frac{h}{2}, y(t+\frac{h}{2})\right)$$

$$y(t+h) = y(t) + h y’(t+\frac{h}{2})$$

But before we start writing this in code, let’s think about the function we want to integrate for a minute or two. In this example, we will use the situation of two species, one of which ($x$) is a prey, the other ($y$) a predator. The equations we use to model this system are:

$$\frac{x’}{x} = r - γ y$$ $$\frac{y’}{y} = \beta x - δ$$

In other words, the prey grows at rate $r$. Every unit of predator removes $\gamma$ units of prey, which are converted into $\beta$ unit of predators; finally, predators decay at rate $\delta$.

Most of these parameters are things we might want to change, so our function should be able to accept them as parameters. We may also want to substitute another problem later, instead of working on the competition model. For this reason, it is better to write a fonction which is as general as possible.

We know that all of our problems (i.e. all of our models) will have two things in common: they will require the time $t$, and the value $u$; the rest can go in keyword arguments. This gives us a template to write the function for the competitive model:

function predation(t::Float64, u::Vector{Float64}; r=1.0, γ=1.0, β=1.0, δ=1.0)
x, y = u
dx = x*r-γ*x*y
dy = β*x*y-δ*y
return [dx, dy]
end

predation (generic function with 1 method)


There are many things that can go wrong with this function – so we may want to add a few checks before we call it. This is an interesting exercise to do, and you can re-read through the lesson on avoiding mistakes to refresh your memory. For now, we will use this basic (but unsafe) version.

Let’s try! If we have a single species (the preys), and we give it any initial population larger than 0, then it should grow (and the population of the predator should not change):

predation(0.0, [0.1, 0.0])

2-element Array{Float64,1}:
0.1
0.0


Nothing stops us from running this function a few times, to figure out if it gives the correct dynamics:

u0 = [0.1, 0.1]
population = [u0]
for i in 2:30
push!(population, population[i-1] + predation(0.0, population[i-1]))
end


We can plot this:

using StatsPlots
plot(hcat(population...)[1,:], c=:teal, lab="Prey",
xlim=(0,30), ylim=(0,5), frame=:origin,
xlab="Time", ylab="Abundance"
)
plot!(hcat(population...)[2,:], c=:purple, lab="Predator")


We obviously need to start thinking about the RK2 integrator! We want our RK2 integrator (a function named rk2) to return two values: an array of times, going from $t_0$ to $t_f$ by steps of $h$, and an array of population sizes at every point. This means that we can start thinking about the arguments: t=(t0,tf) is a tuple with the beginning and end time; h is the integration step. None of these have obvious defaults, but we will still put them as default, keyword arguments, to make sure that we can run our function with a minimum number of keystrokes. If this seems contradictory to the best practices we mentionned in some lessons, it’s because it really is! But best practices are not laws, and sometimes it is worth applying critical judgement.

Our function will therefore look something like

function rk2(;t=0.0=>100.0, h=0.1)
# something
# something else
T = collect(t.first:h:t.second)
return something, T
end


The notation we use for t is called a Pair (and is covered in the Julia documentation). The next step is to think about the “meat” of the function: for every timestep, we want to apply the RK2 equations that are given early in this lesson. Assuming that the current population is stored at position i of an object called U, this would look like:

y_th2 = U[i] .+ 0.5*h .* f(ti, U[i]; p...)
yprime_th2 = f(ti+0.5*h, y_th2; p...)
y_th = U[i] .+ h*yprime_th2
U[i+1] = y_th


There are a few things worth looking at here. First, note that we use the . notation, to make sure that operations are vectorized. Second, we call the f function (which here would be our predation function) using an argument called p.... This syntax allows a function to *capture keyword arguments, and pass them to another function. Let’s have a look:

function a(x; y=2.0, z=3.0)
return x + y + z
end
function b(x; k...)
return a(x; k...)
end

b(1.0) # Should be 6

6.0

b(1.0; y=0.0) # Should be 4

4.0


In short, this allows passing “all the keyword arguments” from a function to another. This is really useful in our context, because although we currently use the competition model, we may want to change it later, possibly to something that accepts different keyword arguments! This gives us the two extra parameters we need to add to rk2: the function f, and the term to catch arguments p....

function rk2(f; t=0.0=>100.0, h=0.1, p...)
# something
for some things
y_th2 = U[i] .+ 0.5*h .* f(ti, U[i]; p...)
yprime_th2 = f(ti+0.5*h, y_th2; p...)
y_th = U[i] .+ h*yprime_th2
U[i+1] = y_th
end
# something else
T = collect(t.first:h:t.second)
return U T
end


Believe it or not, this is almost finished! All we need is to set a value for U, and to iterate over all values of T:

function rk2(u0, f; t=0.0=>100.0, h=0.1, p...)
T = collect(t.first:h:t.second)
U = typeof(u0)[]
push!(U, u0)
for (i,ti) in enumerate(T[2:end])
y_th2 = U[i] + 0.5*h * f(ti, U[i]; p...)
yprime_th2 = f(ti+0.5*h, y_th2; p...)
y_th = U[i] + h*yprime_th2
push!(U, y_th)
end
return (T, hcat(U...)')
end

rk2 (generic function with 1 method)


Let’s try to run the same example as before:

t, u = rk2([0.1, 0.2], predation; t=1.0=>30.0, r=0.9, h=0.05)
plot(t, u[:,1], c=:teal, lab="Prey",
xlim=(0,30), ylim=(0,5), frame=:origin,
xlab="Time", ylab="Abundance"
)
plot!(t, u[:,2], c=:purple, lab="Predator")


Notice that we had previously not used $t$ in the integration – but because our function is general, we can add something like a seasonal effect on the intensity of predation – specifically, the effect of the predator on the prey is a sinusoidal function of time, with “high” and “low” seasons.

function seasonal_predation(t::Float64, u::Vector{Float64}; r=1.0, γ=1.0, β=1.0, δ=1.0)
x, y = u
dx = x*r-(γ+sin(2.0*t)*0.15)*x*y
dy = β*x*y-δ*y
return [dx, dy]
end

seasonal_predation (generic function with 1 method)


This require no changes to the code to produce the figure:

t, u = rk2([0.1, 0.2], seasonal_predation; t=1.0=>200.0, r=0.9, h=1e-2)
plot(t, u[:,1], c=:teal, lab="Prey",
xlim=(0,50), ylim=(0,5), frame=:origin,
xlab="Time", ylab="Abundance"
)
plot!(t, u[:,2], c=:purple, lab="Predator")