I am solving a set of differential equations, see eqs C3a-C3e from https://journals.aps.org/prresearch/pdf/10.1103/PhysRevResearch.5.013091
I separated the equations into real and imaginary parts and propagate them using different solvers. In the beginning, I set the algorithm to nothing and chose alg_hints = :stiff (which corresponds to choosing the solver automatically). Here is how I call the solver:
sol = solve(
prob,
alg,
alg_hints = :stiff,
dtmax = 0.01,
maxiters = Int(1e10),
saveat = 0.01,
reltol = reltol,
abstol = abstol
)
Depending on how I choose abstol and reltol, the dynamics can have weird wiggles or even become unstable. For new configurations, I always need to check which tolerances work best.
For example, I consider 3 particles (~19 real equations) that interact strongly (J_ij ~ 3000) than it is best to choose abstol = 1e-9 and reltol = 1e-8, see screenshot.
However, if I consider 5 particles (~105 real equations) that interact less strongly (J_ij ~ 375), then it is best to choose abstol = 1e-6 and reltol = 1e-5, see screenshot.
My questions are: Is it expected that one has to play with abstol and reltol for different configurations? Is there not an automatic way of determining them?
Can one understand why in one case it is good to have low tolerances and in the other it is not?
Just linking for reference an in-depth discussion about abstol and reltol from a few months ago: Setting `abstol` and `reltol` when solving the Schrödinger equation with OrdinaryDiffEq (note that the answer that’s marked as correct isn’t 100 % correct, so don’t stop reading there).
I have no idea why lowering the tolerance makes your solver unstable, though. Probably a good idea to experiment with different solvers, but it would also be good to know why the default choice screws up like this.
1 Like
You really need to manually choose the solver that works best for your problem. Might take a couple of days because there are ~150 solvers to choose from. And many of them also have extra options. We need to use for example:
solver = FBDF(nlsolve=OrdinaryDiffEqNonlinearSolve.NLNewton(relax=0.4, max_iter=1000))
integrator = OrdinaryDiffEqCore.init(s.prob, solver; dt, abstol=s.set.abs_tol, reltol=s.set.rel_tol, save_on=false, save_everystep=false)
for some problems.
Trying different solvers is often necessary, but even an unsuitable solver usually shouldn’t go unstable when adaptive time stepping is used. It should either be painfully slow or error out with dt < dt_min. It would be good to understand why stability is lost in this case, and especially why it’s triggered by tightening the tolerances.
1 Like
Are you modifying u in f? What is plotted doesn’t look natural, it looks like it’s violating something. And the instabilities around low tolerances is a hallmark that you’re doing something that violates step rejection.
See
