Convergence of SDE solvers

General Usage
,
1

I am wondering how to make sure that an SDE solver has converged on a trajectory level. In other words, I want to make sure that I have chosen abstol and reltol small enough.

What I would do in the deterministic case is that I would just solve the differential equation twice for different tolerances. If the results lie on top of each other, I can be sure that I have chosen the tolerances low enough.

In the stochastic case, however, this is not possible. I am new to these types of equations, but as far as I understand, stochastic solvers only have strong convergence, which means that the difference between two trajectories with different step sizes dt1 and dt2 vanishes as ~ |dt1 - dt2|^p.

I find this a bit unintuitive, as I don’t want my solution to depend on my choice of the time step (or the tolerances). To test this, I took an example from the tutorial and executed it with the same seed but different tolerances:

import StochasticDiffEq as SDE
α = 1
β = 1
u₀ = 1 / 2
f(u, p, t) = α * u
g(u, p, t) = β * u
tspan = (0.0, 1.0)
prob = SDE.SDEProblem(f, g, u₀, tspan)

sol = SDE.solve(prob, SDE.SRIW1(), seed = 123, abstol = 1e-7, reltol = 1e-7, maxiters = Int(1e10), saveat = 0.01)
sol_lower_tol = SDE.solve(prob, SDE.SRIW1(), seed = 123, abstol = 1e-8, reltol = 1e-8, maxiters = Int(1e10), saveat = 0.01)

Plotting the two solutions yields

Comparison of the solutions with different tolerances.
Comparison of the solutions with different tolerances.919×489 37.4 KB

Therefore my questions:

  1. Is it better to play with abstol and reltol or with dt and dtmax to ensure convergence?
  2. How can I be sure that I have chosen the correct solver and that my solver has converged?
  3. Is it not a problem that my solution depends on the solver parameters? Or is it just that the stochastics differ (i.e., that different random numbers are chosen) when I modify the tolerances?
2

Depends on the problem but generally yes.

I would instead state it as, for a chosen brownian motion, it converges like order p.

Those aren’t using the same Brownian motion. You’d want to use things like NoiseWrapper in order to refine the Brownian motion of a different solve.

It’s that the random process is completely different. The random seed only makes sure that your sequence of random numbers is the same. But the Brownian motion is a pairing of dt’s with random numbers, and so if you change the dts you just have a completely different Brownian motion.

3

Thanks for the swift reply. Here is the updated code using DiffEqNoiseProcess:

import StochasticDiffEq as SDE
using DiffEqNoiseProcess
α = 1
β = 1
u₀ = 1 / 2
f(u, p, t) = α * u
g(u, p, t) = β * u
dt = 1 // 2^4
tspan = (0.0, 1.0)
prob1 = SDE.SDEProblem(f, g, u₀, tspan)

sol = SDE.solve(prob1, SDE.SRIW1(), seed = 123, abstol = 1e-6, reltol = 1e-6, maxiters = Int(1e10), saveat = 0.01, save_noise = true)
W2 = NoiseWrapper(sol.W)
prob2 = SDE.SDEProblem(f, g, u₀, tspan, noise = W2)

sol_lower_tol = SDE.solve(prob2, SDE.SRIW1(), seed = 123, abstol = 1e-7, reltol = 1e-7, maxiters = Int(1e10), saveat = 0.01)

Now the curves lie on top of each other:

919×489 29.7 KB

So, summing up, one could say the following:
To make sure that the solver has converged, one can perform the same analysis as in the deterministic case, i.e., solving the dynamics twice for different abstol and reltol. If the curves coincide, than one has chosen a sufficiently low tolerance. However, one must be careful to use the same noise process in both runs.

Is that correct? Or is it ok for the solutions to be different for different tolerances?

4

The simpler answer is to trust that the tolerances are working correctly and just set reltol=1e-6 without bothering to check, but a little paranoia definitely can be a good idea for complex systems like SDEs.

5

Because you never have convergence better than order 1.5 in a strong sense, hitting that kind of tolerance is almost impossible with an SDE solver :sweat_smile: more like 1e-3 is strict.

That is correct. The tolerances are local strong error, kind of like how the tolerances are local error of the ODE solver.

6

I am still trying to get a feeling for SDEs (or maybe I am fundamentally misunderstanding something).

My experience is also that choosing a tolerance of 1e-6 makes the simulation time prohibitively long. However, when I choose a tolerance of 1e-3 - 1e-4, which should be strict, I find the following dynamics for the example script from my previous posts:

919×489 31.7 KB

This looks fairly converged to me, but not quite - or do I have to lower my standards when working with SDEs?

I then did the same analysis for my set of differential equations. I am solving 82 nonlinear equations. Randomness enters through a phase phi which follows a Wiener process. Some of the remaining equations are coupled to the phase like d/dt u[i] = cos(phi) * u [j] + … Additionally, there is a ContinuousCallback, which applies discontinuous changes to u at random times (but the times also depend on the u vector).

When I now compare the dynamics for tolerances between 1e-3 and 1e-4, I find that the phase is quite similar:

919×489 30.2 KB

However, the remaining entries of u look quite different:

919×489 43.6 KB

To me it looks like that due to slight differences of the stochastics the ContinuousCallback is activated at a different times (i.e., the first time around t = 2.8 for the blue curve).

Would you consider this as “converged”?

7

You really have to lower your standards for SDEs: they converge very slowly.

But also, at that point your error is likely to be dominated by sampling error anyways.

That’s a very different model, since an SDE has a discontinuous change at every time point, it’s continuously discontinuous in the derivative (i.e. it has a derivative discontinuity almost everywhere in a measure theoretic sense). Your model has finitely many discontinuities, it’s infinitely less discontinuous.