Документация Engee
Notebook

Решение системы диффренциальных уравнений при различных начальных условиях

В этом примере рассматриваются два подхода к решению дифференциального уравнения (д.у.) при различных начальных условиях (н.у.):

  1. используя цикл for с изменением н.у.
  2. используя EnsembleProblem для решения в многопоточном режиме

Описание задачи

В качестве д.у.  рассмотрим уравнение Лотки-Вольтерры, описывающего модель "хищник-жертва":

где
и - популяции жертвы и хищника, соответственно

- время

- постоянные параметры, описывающие взаимодействие между представителями видов.

Для определённости возьмём их равными .

В поставленной задаче н.у. являются размеры популяции видов. Решая д.у. мы сможем наблюдать изменение популяции во времени.

Рассмотрим, какие пакеты нужны для решения нашей задачи:

  • OrdinaryDiffEq - пакет, позволяющий решать обыкновенные д.у., который может работать без подключения пакета DifferentialEquations
  • Plots имеет возможность отобразить решение дифференциального уравнения буквально одной строчкой
  • Distributed позволит разместить функции на разных ядрах процессора, чтобы повысить эффективность вычисления
In [ ]: 
Pkg.add(["DifferentialEquations", "OrdinaryDiffEq", "Distributed", "BenchmarkTools"])
In [ ]: 
using OrdinaryDiffEq, Plots

Обыкновенное д.у. имеет вид

Для решения задачи необходимо определить следующие исходные данные:

f - саму функцию (структуру д.у.)

u0 - начальные условия

tspan - промежуток времени, на котором отыскивается решение

p - параметры системы, которые также определяют её вид

kwargs - ключевые слова, передаваемые в решение

Есть два способа задать функцию f:

  1. f(du,u,p,t) - эффективный по памяти, но может не поддерживаться пакетами, в которых запрещена мутация
  2. f(u,p,t) - возвращает du. Применима в пакетах, где мутация запрещена.
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

После того, как мы поставили задачу, перейдём к её решению. Функция solve позволяет решать д.у. с применением различных решателей (solver). В нашем случае имеем ОДУ, которое можно решать с постоянным шагом. В качестве параметров будет передаваться поставленная проблема, решатель Tsit5() и параметр saveat=0.01, указывающий функции шаг, с которым необходимо сохранять решение.

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]

Теперь отобразим наше решение, с автоматически указанными легендой и подписью оси t

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

Теперь отобразим фазовую траекторию - изменение фазовых координат (x и y) от времени. На рисунке начальные положения в оттенках красного, а в конце решения - синие.

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]:

Теперь, когда мы смогли изобразить фазовую траекторию для одной пары н.у., построим фазовый портрет - совокупность этих тракторий.

Для этого будем использовать цикл for, постепенно увеличивая начальный размер популяции хищника.

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)

Многопоточное построение фазового портрета

Теперь рассмотрим другой вариант решения задачи построения фазового портрета. Для этого подключим стандартный пакет Distributed, позволяющий выполнять задачи распределённо. И будем распределять нашу задачу на 4 потока (процесса).

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

Подключив DifferentialEquations - пакет для решения д.у. различных типов, Воспользуемся функцией EnsembleProb, позволяющей задать изменение начальных условий посредством вызова функции prob_func.

Например, изменять начальные условия можно случайно:

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

либо же использовать внешний контейнер данных:

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

Поставив задачу, получим решение, используя solve. Но в этот раз нам вернёnся не одно решение д.у., а целый набор симуляции sim. Заметим, что в solve добавились аргументы:

EnsembleThreads() - означающий, что распределение будет по потокам процессора, а не, например, по разным компьютерам кластера
trajectories=N_TRAJECTS - количество симуляций, которые мы хотим провести

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

Отрисуем фазовый портрет, пройдясь по каждому решению в наборе симуляции

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

display(plt)

Сравнение времени выполнения

Теперь сравним время выполнения кода в цикле и с использованием многопоточности.

Для этого подключим библиотеку BenchmarkTools, зададим наши функции, которые мы хоти сравнить, и при помощи макроса @btime наши функции прогонятся несколько раз и будет вычисленно их среднее время выполнения.

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)

Заключение

Исходя из результатов оценки времени выполнения, ensemble_prob_solution оказался быстрее в 4 раза, что соответствует количеству ядер, выделенных под задачу.

Из чего можно сделать вывод, что при решении системы дифференциальных уравнений с различными начальными условиями циклу for стоит предпочитать EnsemleProblem и распределённое вычисление, с использованием пакета Distributed.