julia - 模拟反射边界 SDEProblem

Question

我正在尝试模拟反射边界。根据此处找到的建议：我尝试了 Julia 中带有回调的随机微分方程

using DifferentialEquations
using Plots
using Random

m(x,p,t) -> 0
s(x,p,t) -> 1
x0 = 0.1
tspan = (0.0, 2.5)

prob = SDEProblem(m, s, x0, tspan)

condition(u,t,integrator) = true

function affect!(integrator)
    if integrator.u < 0    
        integrator.u = -integrator.u
    end
end
        
cb = DiscreteCallback(condition,affect!;save_positions=(false,false))

Random.seed!(2001)
sol1 = solve(prob, EM(), dt = 0.001, callback = cb)
plot(sol1)
Random.seed!(2021)
sol2 = solve(prob, EM(), dt = 0.01, callback = cb)
plot(sol2)

分别产生和。我注意到的一件事是，当 dt 较小时，反射的“质量”要好得多。我怀疑这是因为求解器只检查插值中的节点，而不是每个点。

这对选择自己的时间步长的自适应求解器有影响。例如，如果我现在使用求解器运行相同的问题，这是https://diffeq.sciml.ai/stable/solvers/sde_solve/SOSRI上的第一个建议，我得到：

Random.seed!(2021)
sol3 = solve(prob, SOSRI(), callback = cb)
plot(sol3)

可以说反射质量较差。

鉴于问题似乎是仅在结处评估条件，这是 a 的想法DiscreteCallback，我尝试了使用的最后一种方法ContinuousCallback：

condition(u,t,integrator) = u<0

function affect!(integrator)
    if integrator.u < 0    
        integrator.u = -integrator.u
    end
end
        
cb = ContinuousCallback(condition,affect!;save_positions=(false,false))

Random.seed!(2001)
sol4= solve(prob, EM(), dt = 0.001, callback = cb)
plot(sol4)
Random.seed!(2021)
sol5 = solve(prob, SOSRI(), callback = cb)
plot(sol5)

但这不起作用（我想我可能没有正确使用 ContinuousCallback。结果是and ，可以说没有反射

模拟这些过程的推荐方法是什么，以及显式时间步求解器是唯一受支持的方法吗？

Answer 1

这真的只是节省。您拥有它的方式将保存每一步，这意味着它将“保存、反映、保存”。您真正想要的只是反思后的保存：

using StochasticDiffEq
using Plots
using Random

m(x,p,t) = 0
s(x,p,t) = 1
x0 = 0.1
tspan = (0.0, 2.5)

prob = SDEProblem(m, s, x0, tspan)

condition(u,t,integrator) = true

function affect!(integrator)
    if integrator.u < 0
        integrator.u = -integrator.u
    end
end

cb = DiscreteCallback(condition,affect!;save_positions=(false,true))

Random.seed!(2001)
sol1 = solve(prob, EM(), dt = 0.001, callback = cb, save_everystep=false)
plot(sol1)
Random.seed!(2021)
sol2 = solve(prob, EM(), dt = 0.01, callback = cb, save_everystep=false)
plot(sol2)
Random.seed!(2021)
sol2 = solve(prob, SOSRI(), callback = cb, save_everystep=false)
plot(sol2)