在 Julia 中使用 VectorContinuousCallback 将数组的值设为零
Use VectorContinuousCallback to take values of array to zero in Julia
我无法让 VectorContinuousCallback
按预期工作,我不确定自己做错了什么。我有一个大型方程组,基本上,任何时候任何值超过某个阈值(在我的系统中它是 10e-30 但在这个 reprex 0.05 中),我希望该值变为零。
也就是说,如果在任何时候 u 的值低于 0.05,我希望回调将值设为零,但现在,求解器似乎几乎忽略了回调?没有任何阈值交叉被识别。
一个代表:
using DifferentialEquations, Plots
function biomass_sim!(du, u, p, t)
# change = growth + gain from eating - loss from eating - loss
du[1] = 0.2*u[1] + (0.1*u[2] + 0.15*u[3]) - (0.2*u[4]) - 0.9*u[1]
du[2] = 0.2*u[2] + (0.1*u[1] + 0.05*u[3]) - (0.1*u[1] + 0.4*u[4]) - 0.5*u[2]
du[3] = 1.2*u[3] + 0 - (0.15*u[1] + 0.005*u[2]) - 1.3*u[3]
du[4] = 0.2*u[4] + (0.2*u[1] + 0.4*u[2]) - 1.9*u[1]
end
# set up extinction callback
function extinction_threshold(out,u,t,integrator)
# loop through all species to make the condition check all of them
for i in 1:4
out[i] = 0.05 - u[i]
end
end
function extinction_affect!(integrator, event_idx)
# loop again through all species
for i in 1:4
if event_idx == i
integrator.u[i] = 0
end
end
end
extinction_callback =
VectorContinuousCallback(extinction_threshold,
extinction_affect!,
4,
save_positions = (true, true),
interp_points = 1000
)
tspan = (0.0, 10.0)
u0 = [10, 10, 10, 10]
prob = ODEProblem(biomass_sim!,
u0,
tspan)
sol= solve(prob,
Tsit5(),
abstol = 1e-15,
reltol = 1e-10,
callback = extinction_callback,
progress = true,
progress_steps = 1)
plot(sol)
我想在这里看到的是,至少在视觉上超过阈值的两个值 < 0.05(u3 和 u4 明显低于零),我希望这些值变为零。
这个应用是对一个物种的灭绝,所以如果它们低于某个阈值,我想认为它们已经灭绝,因此不能被其他物种消耗。
我尝试过使用不同的求解器更改公差(我没有与 Tsit5() 结合),但尚未找到执行此操作的方法。
非常感谢任何帮助!!
完整输出:
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 165-element Vector{Float64}:
0.0
0.004242012928189618
0.01597661154828718
0.03189297583643294
0.050808376324350105
⋮
9.758563212772982
9.850431863368996
9.94240515017787
10.0
u: 165-element Vector{Vector{Float64}}:
[10.0, 10.0, 10.0, 10.0]
[9.972478301795496, 9.97248284719235, 9.98919421326857, 9.953394015118882]
[9.89687871881005, 9.896943262844019, 9.95942008072302, 9.825050883392004]
[9.795579619066798, 9.795837189358107, 9.919309541156514, 9.652323759097303]
[9.677029343447844, 9.67768414040866, 9.872046236050455, 9.449045554530718]
⋮
[43.86029800419986, 110.54225328286441, -12.173991695732434, -186.40702483057268]
[45.33660997599057, 114.35164541304869, -12.725800474246844, -192.72257104623995]
[46.86398454218351, 118.2922142830212, -13.295579652115606, -199.25572621838901]
[47.84633546050675, 120.82634905853745, -13.661479860003494, -203.45720035095707]
已在 https://github.com/SciML/DifferentialEquations.jl/issues/843 中回答。这是一个“用户错误”。当您检查回调时:
function extinction_affect!(integrator, event_idx)
# loop again through all species
@show integrator.u,event_idx
for i in 1:4
if event_idx == i
integrator.u[i] = 0
end
end
biomass_sim!(get_tmp_cache(integrator)[1], integrator.u, integrator.p, integrator.t)
@show get_tmp_cache(integrator)[1]
end
它确实被调用并且完全符合预期。
(integrator.u, event_idx) = ([5.347462662161639, 5.731062469090074, 7.64667777801325, 0.05000000000000008], 4)
(get_tmp_cache(integrator))[1] = [-2.0231159499021523, -1.3369848518263594, -1.5954424894710222, -6.798261538038756]
(integrator.u, event_idx) = ([12.968499097445866, 30.506371944743357, 0.050000000000001314, -53.521085634736835], 3)
(get_tmp_cache(integrator))[1] = [4.676904953209599, 12.256522670471728, -2.0978067243405967, -20.548116814707996]
但这也表明,即使 u[4] = 0
、du[4] < 0
等也很清楚为什么它变为负数:这是由于 ODE 的定义方式。如果你想在回调点之后保持零,你应该翻转一个参数或其他东西使导数 = 0。
我无法让 VectorContinuousCallback
按预期工作,我不确定自己做错了什么。我有一个大型方程组,基本上,任何时候任何值超过某个阈值(在我的系统中它是 10e-30 但在这个 reprex 0.05 中),我希望该值变为零。
也就是说,如果在任何时候 u 的值低于 0.05,我希望回调将值设为零,但现在,求解器似乎几乎忽略了回调?没有任何阈值交叉被识别。
一个代表:
using DifferentialEquations, Plots
function biomass_sim!(du, u, p, t)
# change = growth + gain from eating - loss from eating - loss
du[1] = 0.2*u[1] + (0.1*u[2] + 0.15*u[3]) - (0.2*u[4]) - 0.9*u[1]
du[2] = 0.2*u[2] + (0.1*u[1] + 0.05*u[3]) - (0.1*u[1] + 0.4*u[4]) - 0.5*u[2]
du[3] = 1.2*u[3] + 0 - (0.15*u[1] + 0.005*u[2]) - 1.3*u[3]
du[4] = 0.2*u[4] + (0.2*u[1] + 0.4*u[2]) - 1.9*u[1]
end
# set up extinction callback
function extinction_threshold(out,u,t,integrator)
# loop through all species to make the condition check all of them
for i in 1:4
out[i] = 0.05 - u[i]
end
end
function extinction_affect!(integrator, event_idx)
# loop again through all species
for i in 1:4
if event_idx == i
integrator.u[i] = 0
end
end
end
extinction_callback =
VectorContinuousCallback(extinction_threshold,
extinction_affect!,
4,
save_positions = (true, true),
interp_points = 1000
)
tspan = (0.0, 10.0)
u0 = [10, 10, 10, 10]
prob = ODEProblem(biomass_sim!,
u0,
tspan)
sol= solve(prob,
Tsit5(),
abstol = 1e-15,
reltol = 1e-10,
callback = extinction_callback,
progress = true,
progress_steps = 1)
plot(sol)
我想在这里看到的是,至少在视觉上超过阈值的两个值 < 0.05(u3 和 u4 明显低于零),我希望这些值变为零。
这个应用是对一个物种的灭绝,所以如果它们低于某个阈值,我想认为它们已经灭绝,因此不能被其他物种消耗。
我尝试过使用不同的求解器更改公差(我没有与 Tsit5() 结合),但尚未找到执行此操作的方法。
非常感谢任何帮助!!
完整输出:
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 165-element Vector{Float64}:
0.0
0.004242012928189618
0.01597661154828718
0.03189297583643294
0.050808376324350105
⋮
9.758563212772982
9.850431863368996
9.94240515017787
10.0
u: 165-element Vector{Vector{Float64}}:
[10.0, 10.0, 10.0, 10.0]
[9.972478301795496, 9.97248284719235, 9.98919421326857, 9.953394015118882]
[9.89687871881005, 9.896943262844019, 9.95942008072302, 9.825050883392004]
[9.795579619066798, 9.795837189358107, 9.919309541156514, 9.652323759097303]
[9.677029343447844, 9.67768414040866, 9.872046236050455, 9.449045554530718]
⋮
[43.86029800419986, 110.54225328286441, -12.173991695732434, -186.40702483057268]
[45.33660997599057, 114.35164541304869, -12.725800474246844, -192.72257104623995]
[46.86398454218351, 118.2922142830212, -13.295579652115606, -199.25572621838901]
[47.84633546050675, 120.82634905853745, -13.661479860003494, -203.45720035095707]
已在 https://github.com/SciML/DifferentialEquations.jl/issues/843 中回答。这是一个“用户错误”。当您检查回调时:
function extinction_affect!(integrator, event_idx)
# loop again through all species
@show integrator.u,event_idx
for i in 1:4
if event_idx == i
integrator.u[i] = 0
end
end
biomass_sim!(get_tmp_cache(integrator)[1], integrator.u, integrator.p, integrator.t)
@show get_tmp_cache(integrator)[1]
end
它确实被调用并且完全符合预期。
(integrator.u, event_idx) = ([5.347462662161639, 5.731062469090074, 7.64667777801325, 0.05000000000000008], 4)
(get_tmp_cache(integrator))[1] = [-2.0231159499021523, -1.3369848518263594, -1.5954424894710222, -6.798261538038756]
(integrator.u, event_idx) = ([12.968499097445866, 30.506371944743357, 0.050000000000001314, -53.521085634736835], 3)
(get_tmp_cache(integrator))[1] = [4.676904953209599, 12.256522670471728, -2.0978067243405967, -20.548116814707996]
但这也表明,即使 u[4] = 0
、du[4] < 0
等也很清楚为什么它变为负数:这是由于 ODE 的定义方式。如果你想在回调点之后保持零,你应该翻转一个参数或其他东西使导数 = 0。