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). But then, I found out that this does not work for some cases, and I went to try some of the solvers proposed for stiff problems here ODE Solvers · DifferentialEquations.jl. I found that Kvaerno5 and KenCarp4 seemed more reliable, but also lead to some strange features in the dynamics in some cases.
For example, I consider a system of 3x3 particles with an interaction for direct neighbors of J_ij + i Gamma_ij ~ 23 + i and perform a second-order approximation (i.e., I only consider equations C3a - C3c and approximate three-particle terms). This corresponds to a system of 117 equations, and the dynamics of the central atom’s population look like what is shown on the screenshot below. Here, I would say that the automatic solver choice works the best (KenCarp4 and Kvaerno lie on top of each other).
Another example is when I consider 3x1 particles with an interaction for direct neighbors of J_ij + i Gamma_ij ~ 3000 + i and perform a third-order approximation (i.e., I do not approximate third-order terms but evolve them according to C3d-C3e). This corresponds to a system of 19 equations.
I know that these two examples are very different in terms of parameters and number of equations. Still, I would like to get an understanding of when to choose which solver - trying all of them for each parameter set seems extensive.
Optimally, I would like to find one reliable solver that works for all situations. Or should I lower my expectations?
Note: When I call the solve function, I use
sol = solve(
prob,
alg,
alg_hints = :stiff,
dtmax = 0.01,
maxiters = Int(1e10),
saveat = 0.001,
)
Maybe some other parameter is not set correctly?