Add BAOAB algorithm #397
Add BAOAB algorithm #397
Conversation
src/alg_utils.jl
Outdated
| @@ -95,6 +95,8 @@ alg_order(alg::SMEB) = 1//1 | |||
| alg_order(alg::TauLeaping) = 1//1 | |||
| alg_order(alg::CaoTauLeaping) = 1//1 | |||
|
|
|||
| alg_order(alg::BAOAB) = 1//1 # Not sure what order it is | |||
There was a problem hiding this comment.
that would be good to look into
|
Needs a lower bound on SciMLBase 1.7 |
| prob = DynamicalSDEProblem(f1_harmonic_iip,f2_harmonic_iip,g_iip,v0,u0,(0.0,5.0); noise=NoiseWrapper(sol1.W)) | ||
|
|
||
| sol2 = solve(prob,BAOAB(gamma=γ);dt=1/10) | ||
| @test sol1[:] ≈ sol2[:] |
There was a problem hiding this comment.
this should get convergence tests
| dts = (1/2) .^ (9:-1:3) | ||
| # Can't use NoiseGrid as noise is not generated with the correct size in convergence.jl. We require noise with shape of v. | ||
| # I don't think there are any analytic solutions for problems of this type | ||
| sim1 = analyticless_test_convergence(dts,prob1,BAOAB(gamma=γ),(1/2)^10;trajectories=Int(1e2),use_noise_grid=false) |
There was a problem hiding this comment.
That's definitely not enough trajectories to produce a weak convergence estimate. Do a local plot with like 5e5 trajectories and show me the plot(sim1).
There was a problem hiding this comment.
Here's the plot with 5e4 trajectories with tspan=(0.0,0.5). I seemed to be running into memory issues if I try more trajectories or a longer tspan. However, it does seem to be converged now. I also tried with 1e4 trajectories and the plot is only marginally different.
Here's the estimate I obtain:
julia> sim1.𝒪est
Dict{Any,Any} with 2 entries:
:weak_final => 1.45594
:final => 0.802602There was a problem hiding this comment.
What if you make the dts a bit smaller? the strong convergence estimate seems to not entire its asymtopic regime until half way through.
There was a problem hiding this comment.
This is with only 1e4 trajectories as it takes a very long time with smaller timesteps. The test dt here is 0.5^12.
julia> sim1.𝒪est
Dict{Any,Any} with 2 entries:
:weak_final => 1.19034
:final => 0.777283This doesn't quite match up with the previous plot however so I guess it's not properly converged. I'm not sure I can run it for many more trajectories though.
There was a problem hiding this comment.
More trajectories won't help. You'd see that because it would flatline in a noisy fashion on the left end due to not having enough values for the sampling error. This doesn't show sampling eror. It shows 1st order convergence.
There was a problem hiding this comment.
Lower dts would be the check to get rid of that wobble though.
There was a problem hiding this comment.
Oh ok, does this mean the algorithm is 1st order then? I don't know much about SDEs, I was just going off the plots in the paper, but that's converging observables calculated from the distributions obtained from the integration. Is that different?
There was a problem hiding this comment.
A lot of the convergence diagrams aren't well-behaved in their paper either. Can you try and recreate one and see if it's around the same?
There was a problem hiding this comment.
It turns out I had made a mistake in the algorithm which I've now fixed. The convergence is a lot more well behaved now and I was able to roughly reproduce one of their plots.
I did 2000 trajectories with 1e6 steps and they did 2000 with 1e8 so I think it's safe to assume I'd converge to the same result given enough timesteps.
Seems like the algorithm is actually just 1st order. Order estimates both strong and weak are 1, even with small dt and many trajectories.
I've made a start on adding the BAOAB algorithm mentioned in issue #259. It relies upon the corresponding pull request I've submitted to SciMLBase that contains the new DynamicalSDEProblem/Function.
Currently I've been able to implement both the inplace and out-of-place versions and both give the same result.
I'm not sure that I've set up the caches/initialisation correctly with the prototypes and such but it works for now.
I've assumed a specific structure of the SDE with:
dq = v*dtdv = f1*dt - gamma*v*dt + g*sqrt(2gamma)*dWWhen f1 is the force, and g is sqrt(kT/m) this is equivalent to Langevin dynamics.
It's possible to make it more general and instead have
dv = f1*dt - gamma*v*dt + g*dW.If you think this form is better I can change it.
Let me know how I should modify things to fit in with everything else.