# Understanding abstol and reltol in DifferentialEquations.jl

**URL:** <https://discourse.julialang.org/t/understanding-abstol-and-reltol-in-differentialequations-jl/129250>\
**Category:** General Usage\
**Tags:** differentialequation\
**Created:** [May 22, 2025, 4:47pm UTC](https://discourse.julialang.org/t/understanding-abstol-and-reltol-in-differentialequations-jl/129250 "2025-05-22T16:47:53Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![diffeqslvr](https://avatars.discourse-cdn.com/v4/letter/d/f475e1/32.png) [@diffeqslvr](https://discourse.julialang.org/u/diffeqslvr)\
**Post date:** [May 22, 2025, 4:47pm UTC](https://discourse.julialang.org/t/understanding-abstol-and-reltol-in-differentialequations-jl/129250/1 "2025-05-22T16:47:53Z")

</div>

I am solving a set of differential equations, see eqs C3a-C3e from [https://journals.aps.org/prresearch/pdf/10.1103/PhysRevResearch.5.013091](https://journals.aps.org/prresearch/pdf/10.1103/PhysRevResearch.5.013091)

 ![Screenshot 2025-05-22 at 18.30.45](https://global.discourse-cdn.com/julialang/original/3X/9/a/9ae726c564fa4a42b28ab7c49f3638038c416a2a.png)

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:

```julia
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.

 ![Screenshot 2025-05-22 at 18.41.59](https://global.discourse-cdn.com/julialang/original/3X/0/e/0ec11713246ebcee03d288fed55988a48ee3766a.png)

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.

 ![Screenshot 2025-05-22 at 18.46.16](https://global.discourse-cdn.com/julialang/original/3X/e/e/ee19e13a271b8fbaa9023edc8782a38ee2bde9b0.png)

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?

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [May 22, 2025, 5:12pm UTC](https://discourse.julialang.org/t/understanding-abstol-and-reltol-in-differentialequations-jl/129250/2 "2025-05-22T17:12:27Z")

</div>

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](https://discourse.julialang.org/t/setting-abstol-and-reltol-when-solving-the-schrodinger-equation-with-ordinarydiffeq/125534) (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.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [May 22, 2025, 5:45pm UTC](https://discourse.julialang.org/t/understanding-abstol-and-reltol-in-differentialequations-jl/129250/3 "2025-05-22T17:45:00Z")

</div>

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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [May 22, 2025, 7:16pm UTC](https://discourse.julialang.org/t/understanding-abstol-and-reltol-in-differentialequations-jl/129250/4 "2025-05-22T19:16:12Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 22, 2025, 8:27pm UTC](https://discourse.julialang.org/t/understanding-abstol-and-reltol-in-differentialequations-jl/129250/5 "2025-05-22T20:27:59Z")

</div>

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

[![](https://global.discourse-cdn.com/julialang/original/3X/5/d/5d8be9be2d861eecc910a9bdda15b5b546dd839f.jpeg "How to Debug Julia Simulation Codes (ODEs, Optimization, etc.!) | Chris Rackauckas | JuliaCon 2022") ](https://www.youtube.com/watch?v=g-iOOhh2U6o)

---

<div class="post-metadata">

**Author:** ![diffeqslvr](https://avatars.discourse-cdn.com/v4/letter/d/f475e1/32.png) [@diffeqslvr](https://discourse.julialang.org/u/diffeqslvr)\
**Post date:** [May 23, 2025, 12:34pm UTC](https://discourse.julialang.org/t/understanding-abstol-and-reltol-in-differentialequations-jl/129250/6 "2025-05-23T12:34:41Z")

</div>

Thanks for all the input!

> You really need to manually choose the solver that works best for your problem.

I have actually tried a few different solvers. So far, I have not found any solver that worked for all configurations I have considered. Also, sometimes, high tolerances work better than low tolerances, see screenshot for a 3x3 lattice.

 ![Screenshot 2025-05-23 at 14.27.30](https://global.discourse-cdn.com/julialang/original/3X/e/f/ef3f7baabfbe3692e30ee4915c196eff0b431060.png)

> 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.

No I am not modifying `u`, just computing `du`. It is entirely possible that I made a mistake somewhere, but I compared my results to some cases shown in the paper, where I got the equations from, and my code matched their results.

It is important to mention that we expect there to be some instabilities in the equations (as also mentioned in the paper). I am working on removing those by performing a projection between propagation steps. In that case, I would really modify `u` between time steps - but I am not at this point yet.

What troubles me is that I don’t know when the instabilities or the unphysical wiggles come from my solver choice, and when the equations are in fact unstable.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 23, 2025, 12:39pm UTC](https://discourse.julialang.org/t/understanding-abstol-and-reltol-in-differentialequations-jl/129250/7 "2025-05-23T12:39:10Z")

</div>

> [@diffeqslvr](#):
>
> What troubles me is that I don’t know when the instabilities or the unphysical wiggles come from my solver choice, and when the equations are in fact unstable.

That doesn’t look like a possible solver instability. That looks like something you’d get if you violated one of the principles of the ODE solver’s assumptions. Did you check the whole list?

---

<div class="post-metadata">

**Author:** ![diffeqslvr](https://avatars.discourse-cdn.com/v4/letter/d/f475e1/32.png) [@diffeqslvr](https://discourse.julialang.org/u/diffeqslvr)\
**Post date:** [May 23, 2025, 1:30pm UTC](https://discourse.julialang.org/t/understanding-abstol-and-reltol-in-differentialequations-jl/129250/8 "2025-05-23T13:30:19Z")

</div>

Yes, as far as I can tell I did.

- I do not modify u
- There is no stochastic term in the equation
- I do not use values from “previous” time steps
- There is not unsmoothness like an if condition
- There are some constraints like that the population of the ground state `ng` and the population of the excited state `ne` need to sum up to 1. However, I completely remove `ng` by replacing it `ng = 1 - ne`.

Here is the my function

```julia
function eoms!(du, u, p, t)
    slv, sys, params = p
    # Extract the atomic operators from u
    σee, σeg_σge, σee_σee, σee_σeg_σge, σee_σee_σee = to_sigma(slv, sys, u)
    # Use multithreading to loop over all complex equations
    @threads for i = 1:slv.n_cmplx_eqs
        if i <= sys.n1_cmplx
            # Get indices pertaining to the ith equation 
            # E.g., the first equation pertains to the population of the
            # first atom σ_ee^1, whose indices are (e, e, 1)
            inds = System.inds_from_equationind(sys, i)
            res = eom1(σee, σeg_σge, inds, params)
        elseif i <= (sys.n1_cmplx + sys.n2_cmplx)
            # Equations for second-order operators, e.g., σ_eg^1 σ_ge^2
            inds = System.inds_from_equationind(sys, i)
            res = eom2(
                σee,
                σeg_σge,
                σee_σ_ee,
                σee_σeg_σge,
                inds,
                params
            )
        else
            # Equations for third-order operators
            inds = System.inds_from_equationind(sys, i)
            res = eom3(
                σee,
                σeg_σge,
                σee_σee,
                σee_σeg_σge,
                σee_σee_σee,
                inds,
                params
            )
        end
        # Get the indices of u that pertain to the ith equation
        uinds = slv.equationind_to_uind[i]
        # If the equation is real, it pertains to 1 index of u, 
        # if it is complex then it pertains to 2
        if length(uinds) == 1
            du[uinds[1]] = real(res)
        else
            du[uinds[1]] = real(res)
            du[uinds[2]] = imag(res)
        end
    end
end

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 23, 2025, 2:14pm UTC](https://discourse.julialang.org/t/understanding-abstol-and-reltol-in-differentialequations-jl/129250/9 "2025-05-23T14:14:48Z")

</div>

Is your function holomorphic?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 23, 2025, 2:16pm UTC](https://discourse.julialang.org/t/understanding-abstol-and-reltol-in-differentialequations-jl/129250/10 "2025-05-23T14:16:17Z")

</div>

I think to help more I’ll need something I can run.

---

<div class="post-metadata">

**Author:** ![Brad\_Carman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/brad_carman/32/17631_2.png) [@Brad\_Carman](https://discourse.julialang.org/u/Brad_Carman)\
**Post date:** [August 5, 2025, 9:31am UTC](https://discourse.julialang.org/t/understanding-abstol-and-reltol-in-differentialequations-jl/129250/11 "2025-08-05T09:31:28Z")

</div>

Please check out [bradcarman/ModelingToolkitTolerances.jl: A package to help determine `abstol` and `reltol` when using ModelingToolkit.jl and DifferentialEquations.jl](https://github.com/bradcarman/ModelingToolkitTolerances.jl). Something I’ve been working on that I hope can help. Should be in the `General` registry in a couple days.