在各种初始条件下求解微分方程组

在这个例子中，我们考虑在不同初始条件下求解微分方程的两种方法：1。
使用循环 for 随着n.o的变化。
2. 使用 EnsembleProblem 对于多线程解决方案

任务说明

作为一个例子，考虑[Lotka-Volterra]方程（https://ru.wikipedia.org/wiki/Model_lotki_-_Volterra），描述捕食者-猎物模型:

哪里
和 -猎物和捕食者的种群，分别

-时间

-描述物种代表之间相互作用的恒定参数。

为了确定起见，让我们把它们视为平等 .

在给定的任务中，n.o.是物种种群的大小。通过解决问题，我们将能够观察人口随时间的变化。

让我们来看看需要哪些软件包来解决我们的问题。:

OrdinaryDiffEq -一个包，允许您解决普通问题，无需连接包即可工作 DifferentialEquations
Plots 它能够在一行中显示微分方程的解。
Distributed 允许您将函数放置在不同的处理器内核上，以提高计算效率。

Pkg.add(["DifferentialEquations", "OrdinaryDiffEq", "Distributed", "BenchmarkTools"])

using OrdinaryDiffEq, Plots

一个普通的d.u.有这样的形式

要解决问题，有必要确定以下初始数据:

f -功能本身（设备的结构）

u0 -初始条件

tspan -寻求解决方案的时间段

p -也决定其外观的系统参数

kwargs -传递给解决方案的关键字

有两种方法可以设置函数 f:

f(du,u,p,t) -内存效率高，但可能不支持[软件包]（https://fluxml.ai/Zygote ...jl/dev/limitations/#Array-mutation-1），其中[mutation]被禁止（https://docs.julialang.org/en/v1/manual/variables/#man-assignment-expressions ）
f(u,p,t) -回报 du. 适用于禁止突变的包装。

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)

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，向函数指示保存解决方案的步骤。

sol = solve(prob, Tsit5(), saveat=0.01)

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

plot(sol)

现在让我们显示相位轨迹-相位坐标的变化（x 和 y）不时。在图中，初始位置是红色的色调，并且在溶液的末端是蓝色的。

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

现在我们已经能够描绘一对车辆的相位轨迹，我们将构建这些轨迹的总体相位肖像。

为此，我们将使用循环 for 逐渐增加捕食者种群的初始大小。

tspan=(0., 15.)
initial_conditions = range(10, stop=400, length=40)

10.0:10.0:400.0

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个线程（进程）中。

using Distributed
N_PROCS = 4  # зададим количество процессов
addprocs(N_PROCS)

4-element Vector{Int64}:
 2
 3
 4
 5

通过启用 DifferentialEquations -解决各种类型问题的软件包，我们将使用该功能 EnsembleProb，它允许您通过调用函数来设置初始条件的更改 prob_func.

例如，可以随机地改变初始条件。:
``'茱莉亚
函数prob_func(prob,i,repeat)
@. prob.u0=randn()*prob.u0prob
结束

或者使用外部数据容器:
``'茱莉亚
函数prob_func(prob,i,repeat)

    重拍（prob，u0=[50，initial_conditions[i]]）

结束

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)

EnsembleProblem with problem ODEProblem

设置任务后，我们得到了一个解决方案，使用 solve. 但这一次，我们将得到的不仅仅是一个解决方案，而是一整套模拟。 sim. 请注意，在 solve 添加参数:

EnsembleThreads（）-意味着分布将跨处理器线程，而不是，例如，跨不同的群集计算机
trajectories=N_TRAJECTS -我们要运行的模拟数量

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

让我们绘制一个相图，通过模拟集中的每个解决方案。

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

display(plt)

执行时间比较

现在让我们比较循环和使用多线程的代码执行时间。

为此，请连接库 BenchmarkTools. 让我们设置我们要比较的函数，并使用宏 @btime 我们的函数将运行多次，并计算它们的平均执行时间。

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.