模拟反射边界 SDEProblem

Question

我正在尝试模拟反射边界。根据此处找到的建议：我试过了

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)

分别产生 and 。我注意到的一件事是，当 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。结果是 and ，可以说没有反射

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

Answer 1

真的只是省钱。你拥有它的方式会保存每一步，这意味着它会“保存、反映、保存”。你真正想要的只是 post-反射保存：

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)

模拟反射边界 SDEProblem

Simulating a reflecting boundary SDEProblem

julia

differentialequations.jl