Engee documentation
Notebook

Solution of the system of diffrential equations at different initial conditions

This example considers two approaches to solving the differential equation (d.u.) under different initial conditions (n.u.):

  1. using the loop for with changing n.u.
  2. using EnsembleProblem to solve in multithreaded mode

Problem description

As a d.u., consider the Lotka-Volterra equation describing the predator-prey model:

$$\begin{cases} \frac{d x}{dt} = \alpha x - \beta x y \\ \frac{d y}{dt} = \delta xy - \gamma y \\ \end{cases}$$

where $x$ and $y$ are the victim and predator populations, respectively

$t$ - time

$\alpha, \beta, \delta, \gamma$ - constant parameters describing the interaction between representatives of species.

For definiteness, let us take them equal to $\alpha = \gamma = 1, \beta =0.01, \delta = 0.02$.

In this problem, the n.d. is the size of the species population. By solving the d.u., we can observe the population change over time.

Let us consider what packages are needed to solve our problem:

  • OrdinaryDiffEq is a package that allows us to solve ordinary d.u., which can work without a package connection DifferentialEquations
  • Plots has the ability to display the solution of a differential equation literally in one line
  • Distributed will allow you to place functions on different processor cores to increase the efficiency of calculation
In [ ]: 
Pkg.add(["DifferentialEquations", "OrdinaryDiffEq", "Distributed", "BenchmarkTools"])
In [ ]: 
using OrdinaryDiffEq, Plots

An ordinary d.u. has the form $$\frac{du}{dt}=f(u,,p,t)$$

To solve the problem it is necessary to determine the following initial data:

f - the function itself (d.u. structure)

u0 - initial conditions

tspan - time interval at which the solution is found

p - parameters of the system, which also determine its form

kwargs - keywords passed to the solution

There are two ways to specify the function f:

  1. f(du,u,p,t) - efficient from memory, but may not be supported by packages that do not allow mutation
  2. f(u,p,t) - Returns du. Applicable in packages where mutation is forbidden.
In [ ]: 
function lotka(du, u, p, t) 
    du[1] = p[1] * u[1] - p[2] * u[1] * u[2]
    du[2] = p[4] * u[1] * u[2] - p[3] * u[2]  
end
α = 1; β = 0.01; γ = 1; δ = 0.02;
p = [α, β, γ, δ]
tspan = (0.0, 6.5)
u0 = [50; 50]
prob = ODEProblem(lotka, u0, tspan, p)
Out[0]:
ODEProblem with uType Vector{Int64} and tType Float64. In-place: true
timespan: (0.0, 6.5)
u0: 2-element Vector{Int64}:
 50
 50

After we have set the problem, let's proceed to its solution. The function solve allows solving ODEs using different solvers (solver). In our case we have an ODE that can be solved with a constant step. As parameters we will pass the set problem, the solver Tsit5() and the parameter saveat=0.01, indicating the step with which the solution should be saved.

In [ ]: 
sol = solve(prob, Tsit5(), saveat=0.01)
Out[0]:
retcode: Success
Interpolation: 1st order linear
t: 651-element Vector{Float64}:
 0.0
 0.01
 0.02
 0.03
 0.04
 0.05
 0.06
 0.07
 0.08
 0.09
 0.1
 0.11
 0.12
 ⋮
 6.39
 6.4
 6.41
 6.42
 6.43
 6.44
 6.45
 6.46
 6.47
 6.48
 6.49
 6.5
u: 651-element Vector{Vector{Float64}}:
 [50.0, 50.0]
 [50.25062391161723, 50.00125211549466]
 [50.50249138370918, 50.00501692739517]
 [50.755595782113325, 50.01130738411668]
 [51.009930152636706, 50.020136767794746]
 [51.26548722105585, 50.031518694285374]
 [51.52225939311683, 50.04546711316499]
 [51.78023875453521, 50.06199630773043]
 [52.039417070996095, 50.081120894998996]
 [52.29978578815409, 50.10285582570839]
 [52.56133603163334, 50.127216384316746]
 [52.82405860702751, 50.15421818900265]
 [53.08794399989971, 50.18387719166507]
 ⋮
 [47.63881663670754, 50.10217268569618]
 [47.877146554422595, 50.07972805051839]
 [48.11676888713707, 50.05968794764736]
 [48.35768016469967, 50.04206191833404]
 [48.59987678709717, 50.026859716639514]
 [48.84335502445415, 50.01409130943499]
 [49.08811101703379, 50.00376687640185]
 [49.33414077523645, 49.99589681003162]
 [49.581440179601316, 49.99049171562596]
 [49.83000498080491, 49.98756241129669]
 [50.07983079966195, 49.98711992796579]
 [50.33091312712557, 49.98917550936536]

Now let's display our solution, with automatically specified legend and axis signature t

In [ ]: 
plot(sol)
Out[0]:

Now let's display the phase trajectory - the variation of the phase coordinates (x and y) with time. In the figure, the initial positions are in shades of red, and at the end of the solution they are in shades of blue.

In [ ]: 
scatter(
sol[1, :], sol[2, :],                           # x и y
legend=false, marker_z=sol.t, color=:RdYlBu,    # marker_z позволяет изменять цвет маркера от t
xlabel="Жертва x", ylabel="Хищник y",           # RdYlBu - от Красного до Синего через Жёлтый
title="Фазовая траектория системы")
Out[0]:

Now that we have been able to depict the phase trajectory for one pair of n.u., let us construct a phase portrait - the set of these trajectories.

For this purpose, we will use the cycle for, gradually increasing the initial size of the predator population.

In [ ]: 
tspan=(0., 15.)
initial_conditions = range(10, stop=400, length=40)
Out[0]:
10.0:10.0:400.0
In [ ]: 
plt = plot()
for y0 in initial_conditions
    u0 = [50; y0]
    prob = ODEProblem(lotka, u0, tspan, p)
    sol = solve(prob, Tsit5(), saveat=0.01)
    plot!(plt, sol[1, :], sol[2, :], title="Фазовый портрет системы", xlabel="Жертва", ylabel="Хищник" ,legend=false)
end

display(plt)

Multithreaded construction of phase portrait

Now let's consider another variant of the phase portrait construction task. For this purpose, let's connect the standard package Distributed, which allows to perform tasks in a distributed manner. And we will distribute our task into 4 threads (processes).

In [ ]: 
using Distributed
N_PROCS = 4  # зададим количество процессов
addprocs(N_PROCS)
Out[0]:
4-element Vector{Int64}:
 2
 3
 4
 5

Having connected DifferentialEquations - a package for solving various types of d.u., we will use the function EnsembleProb, which allows us to set the change of initial conditions by calling the function prob_func.

For example, it is possible to change the initial conditions randomly:

function prob_func(prob, i, repeat) 
@. prob.u0 = randn() * prob.u0 prob
end

Or use an external data container:

function prob_func(prob, i, repeat)

    remake(prob, u0=[50, initial_conditions[i]])

end
In [ ]: 
using DifferentialEquations
u0 = [50; 10];
prob = ODEProblem(lotka, u0, tspan, p);

function prob_func(prob, i, repeat)
    remake(prob, u0=[50, initial_conditions[i]])
end
ensemble_prob = EnsembleProblem(prob, prob_func=prob_func)
Out[0]:
EnsembleProblem with problem ODEProblem

Having posed the problem, we get the solution using solve. But this time we are returned not a single d.u. solution, but a whole set of simulations sim. Note that we have added arguments to solve:

EnsembleThreads() - meaning that the distribution will be over processor threads and not, for example, over different computers in the cluster trajectories=N_TRAJECTS - number of simulations we want to run

In [ ]: 
N_TRAJECTS = length(initial_conditions)
sim = solve(ensemble_prob, Tsit5(), EnsembleThreads(), trajectories=N_TRAJECTS, saveat=0.01);

Draw a phase portrait by traversing each solution in the simulation set

In [ ]: 
plt = plot()
for sol in sim.u
    plot!(plt, sol[1, :], sol[2, :],  title="Фазовый портрет системы", xlabel="Жертва", ylabel="Хищник" ,legend=false)
end

display(plt)

Execution time comparison

Now let's compare the code execution time in a loop and using multithreading.

To do this, let's connect the library BenchmarkTools, set our functions that we want to compare, and using the macro @btime our functions will be run several times and their average execution time will be calculated.

In [ ]: 
using BenchmarkTools

function for_loop_solution()
    for y0 in range(10, stop=400, length=20)
        u0 = [50; y0]
        prob = ODEProblem(lotka, u0, tspan, p)
        solve(prob, Tsit5(), saveat=0.01)
    end
end

function ensemble_prob_solution()
    prob = ODEProblem(lotka, u0, tspan, p)

    initial_conditions = range(10, stop=400, length=20)
     function prob_func(prob, i, repeat)
        remake(prob, u0=[50, initial_conditions[i]])
    end
    ensemble_prob = EnsembleProblem(prob, prob_func=prob_func)
    solve(ensemble_prob, Tsit5(), EnsembleThreads(), trajectories=20, saveat=0.01)
end

@btime for_loop_solution()
@btime ensemble_prob_solution()
[];  #заглушка чтобы не выводить EnsembleProblem в ячейку вывода

  8.418 ms (33260 allocations: 3.59 MiB)
  1.869 ms (31424 allocations: 3.53 MiB)

Conclusion

Based on the results of the runtime evaluation, ensemble_prob_solution was 4 times faster, which corresponds to the number of cores allocated to the task.

From what we can conclude that when solving the system of differential equations with different initial conditions for cycle should be preferred to EnsemleProblem and distributed computation using the package Distributed.