Functional form for BAOAB solver

[jebej]

From PR #397, it is stated that the functional form assumed for the SDE with the BAOAB solver is

dq = v*dt
dv = f1*dt - gamma*v*dt + g*sqrt(2gamma)*dW

This is not what I would expect, since it has an extra -gamma*v*dt term. Is there a reason for this difference?

Comparing with the ODE case in the following script, the solutions are different:

using OrdinaryDiffEq, StochasticDiffEq, PyPlot

jj!(ϕ′,ϕ,τ,β,I) = (I - ϕ′ - sin(ϕ))/β
jj!(ϕ′,ϕ,p::Tuple,τ) = jj!(ϕ′,ϕ,τ,p[1],p[2])
jj!(ϕ′::Vector,ϕ::Vector,p::Tuple,τ) = jj!(ϕ′[],ϕ[],τ,p[1],p[2]) # special version to workaround #517 (fixed on master)
identity_f(v, u, p, t) = v # needed to form second order dynamical ODE

setup_jj(β,I,ϕ_init=0.0,ϕ′_init=0.0,tspan=(0.0,50.0)) = DynamicalODEProblem{false}(jj!,identity_f,ϕ′_init,ϕ_init,tspan,(β,I))

g(u,p,t) = 0.01 # small noise
setup_noisy_jj(β,I,ϕ_init=[0.0],ϕ′_init=[0.0],tspan=(0.0,50.0)) = DynamicalSDEProblem{false}(jj!,identity_f,g,ϕ′_init,ϕ_init,tspan,(β,I))

function analyze_junction_dc(res, clear=true)
    noisy = res isa RODESolution
    V = getindex.(res.u,1)
    I = sin.(getindex.(res.u,2))
    clear && clf()
    plot(res.t, V, label="\$V/φ₀\$"*ifelse(noisy," (Langevin)"," (noise-free)"))
    plot(res.t, I, label="\$I/I_c\$"*ifelse(noisy," (Langevin)"," (noise-free)"))
    xlabel("Time τ")
    legend()
    grid(true)
    display(gcf())
end

prob = setup_jj(1.25,2.0)
res = solve(prob,VelocityVerlet(),dt=1/100)
analyze_junction_dc(res)

prob = setup_noisy_jj(1.25,2.0)
res = solve(prob,BAOAB(1),dt=1/100)
analyze_junction_dc(res, false)

[jebej]

After reading the paper a little bit, I now understand that the point of the BAOAB method is to include this term in the middle "O" operation for better accuracy. As a compromise for applications different than molecular sampling, I'd suggest going back to @jamesgardner1421's second proposition, where gamma is only included in the gamma*v*dt term, and not in the noise term. That way, one can choose gamma=0 to disable this extra operation, or alternatively, it lets people choose exactly how to define the noise term without the current forced gamma dependence.

For my application, for example, where I need:

jj!(ϕ′,ϕ,τ,β,I) = (I - ϕ′ - sin(ϕ))/β

I can use

jj!(ϕ′,ϕ,τ,β,I) = (I - sin(ϕ))/β

with gamma = 1/β. However, I do not want the noise term to include this β.

I can make a PR, if that is acceptable.

[ChrisRackauckas]

What about a gamma and a beta, separate terms?

[jebej]

Fixed by #520

[ChrisRackauckas]

need to update docs?

[jebej]

There weren't any docs for this algorithm. I actually found it by chance, going through the code. As far as I understand, this is the only algo for second-order SDEs, and there isn't even a section in the docs for that.