Stiff SDE and `DiffEq`

I’m trying to integrate an SDE with constant, additive, diagonal noise and a “stiff” drift component (the ODE component comes from a singularly perturbed DAE system…). I need to be able to integrate indefinitely (I’m interested in when the particle reaches some faraway region), so I think that adaptive time-stepping methods won’t work for this problem (but I may be wrong).

I’ve looked through the SDE solvers page and was thinking of TangXiaoSROCK2 but see that it’s under development. So, I’m wondering what integrators are recommended for problems like this (ideally any explicit integrators). Since the diffusion term is constant, Milstein methods are identical to Euler-Maruyama (with strong order 1 convergence), and I’m wondering if I can just use a “better” (explicit) ODE method for the drift term (to deal with the stiffness) while handling the diffusion in the same way as Euler-Maruyama.

This is precisely the case where adaptive methods perform better. The transition can be a big constraint to the time steps, so you want to adapt during transient behavior. Also, any implicit method needs adaptivity since Newton iterations can diverge. Anyways, any method can be set to fixed timesteps with adaptive=false.

I would give SKenCarp a try. There’s something I want to try with its error estimator so give it a try and if it’s not doing well I might want to play with it.

For this style, SROCK2 might be the best bet? Those methods are brand new so we haven’t benchmarked between them yet, so I wouldn’t be so caught up in TangXiaoSROCK2 (there seems to be an issue with the tableau in the paper, IIRC it doesn’t satisfy the order conditions).

And if you need to brute force it, SOSRA should get the job done, but not efficiently. Start giving these a try and let us know if you need more help.

1 Like

Thanks @ChrisRackauckas! By “indefinitely”, I mean that I can’t specify an end time in the problem, as I don’t know how long my simulation will take to reach this region. Am I able to handle such cases with SKenCarp? I tried looking through the documentation, but it appears that an end time is required for each SDEProblem.

You’re looking for the terminate! command in the callbacks. An example is shown here:

http://docs.juliadiffeq.org/latest/features/callback_functions.html#Example-2:-Terminating-an-Integration-1

Thanks, though this still requires an end time, right? My issue is that a priori I don’t know how long the simulation should run. Is there an effective way to set the end time in tspan as something like infinity? I feel that perhaps I’m missing something…

Yes, just use tspan = (0.0,Inf) or something like that.

1 Like