不同初始条件下的差分方程组解法¶
本例考虑了在不同初始条件（n.u.）下求解微分方程（d.u.）的两种方法：
- 使用循环
for并改变 n.u. - 使用
EnsembleProblem以多线程模式求解
问题描述¶
考虑描述捕食者-猎物模型的 Lotka-Volterra 公式：
$$\begin{cases} \frac{d x}{dt} = \alpha x - \beta x y \\ \frac{d y}{dt} = \delta xy - \gamma y \\ \end{cases}$$
其中 $x$ 和$y$ 分别是受害种群和捕食种群
$t$ - 时间
$\alpha, \beta, \delta, \gamma$ - 常数参数，描述物种代表之间的相互作用。
为明确起见，让我们把它们取为$\alpha = \gamma = 1, \beta =0.01, \delta = 0.02$ 。
在这个问题中，n.d.是物种种群的规模。通过求解 d.u.，我们可以观察种群随时间的变化。
让我们考虑一下需要哪些软件包来解决我们的问题：
-OrdinaryDiffEq 是一个可以求解普通 d.u.的软件包，无需软件包插件即可运行DifferentialEquations
-Plots 可以在一行中显示微分方程的解
-Distributed 允许将函数放在不同的处理器内核上，以提高计算效率
Pkg.add(["DifferentialEquations", "OrdinaryDiffEq", "Distributed", "BenchmarkTools"])
using OrdinaryDiffEq, Plots
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)
设定问题后，我们开始求解。函数solve 允许使用不同的求解器（solver）求解 ODE。在我们的例子中，我们有一个可以用常数步求解的 ODE。作为参数，我们将传递问题集、求解器Tsit5() 和参数saveat=0.01 ，后者表示保存解的步长。
sol = solve(prob, Tsit5(), saveat=0.01)
现在让我们显示解决方案，并自动指定图例和坐标轴签名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="Фазовая траектория системы")
既然我们已经能够描绘出一对 n.u. 的相位轨迹，那么让我们构建一个相位图--这些轨迹的集合。
为此，我们将使用循环for ，逐渐增加捕食者种群的初始规模。
tspan=(0., 15.)
initial_conditions = range(10, stop=400, length=40)
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)
连接DifferentialEquations - 用于求解各种类型 d.u. 的软件包，然后使用函数EnsembleProb ，通过调用函数prob_func 来设置初始条件的变化。
例如，可以随机改变初始条件： ``朱利亚 函数 prob_func(prob, i, repeat) @.prob.u0 = randn() * prob.u0 prob 结束
或使用外部数据容器：
```julia
函数 prob_func(prob, i, repeat)
remake(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)
提出问题后，我们使用solve 得到解。但这次返回给我们的不是单一的 d.u. 解，而是一整套模拟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 в ячейку вывода
结论¶
根据运行时间评估的结果，ensemble_prob_solution 快 4 倍，这与分配给任务的内核数量相符。
由此我们可以得出结论，在求解具有不同初始条件的微分方程系统时，for 循环应优于EnsemleProblem 和使用软件包Distributed 进行的分布式计算。