模拟反射边界 SDEProblem

Question

我正在尝试模拟反射边界。根据此处找到的建议：Stochastic differential equation with callback in 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 较小时，反射的“质量”要好得多。我怀疑这是因为求解器只检查插值中的节点，而不是每个点。

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

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

可以说反射质量较差。

鉴于问题似乎是只在结处评估条件，这是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。结果是 和 ，可以说没有反射

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

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)