Simulating a reflecting boundary SDEProblem

Question

I'm trying to simulate a reflecting boundary. Based on the suggestions found here: Stochastic differential equation with callback in Julia I tried

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)

Which yield and , respectively. One thing that I noticed is that the "quality" of the reflection is much better when dt is smaller. I suspect that this is because the solver only checks at knots in the interpolation, rather than at every point.

This has implications for adaptive solvers, which choose their own time steps. For example, if I now run the same problem using the SOSRI solver, which is the first recommendation on https://diffeq.sciml.ai/stable/solvers/sde_solve/, I get:

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

which features arguably poorer quality reflection.

Given that the problem seems to be that the condition is only evaluated at the knots, which is the idea of a DiscreteCallback, I tried one final approach of using 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)

but this didn't work (I think I'm likely not using ContinuousCallback correctly. The results were and , which feature arguably no reflection

What is the recommended way to simulate these processes, and are explicit timestepping solvers the only supported ones?

Answer 1

It's really just saving. The way you had it would save every step, which means it would "save, reflect, save". What you really want are just the post-reflection saves:

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)