# Setting \`abstol\` and \`reltol\` when solving the Schrödinger equation with OrdinaryDiffEq

**URL:** <https://discourse.julialang.org/t/setting-abstol-and-reltol-when-solving-the-schrodinger-equation-with-ordinarydiffeq/125534>\
**Category:** Quantum\
**Tags:** question, ode\
**Created:** [February 4, 2025, 2:19pm UTC](https://discourse.julialang.org/t/setting-abstol-and-reltol-when-solving-the-schrodinger-equation-with-ordinarydiffeq/125534 "2025-02-04T14:19:04Z")\
**Posts on this page:** 8\
**Page:** 2

<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:** [February 5, 2025, 11:07pm UTC](https://discourse.julialang.org/t/setting-abstol-and-reltol-when-solving-the-schrodinger-equation-with-ordinarydiffeq/125534/21 "2025-02-05T23:07:41Z")

</div>

Yes that is correct.

> [@danielwe](#):
>
> (I think this is used because DifferentialEquations allows you to change the length of your state vector during the evolution)

This gives the property that if you have a diagonal ODE and you add more of the same ODE, then you get the same step sizes. It would be weird if bigger ODEs just naturally have smaller dt even if it’s the same equation just repeated a few times.

---

<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:** [February 5, 2025, 11:11pm UTC](https://discourse.julialang.org/t/setting-abstol-and-reltol-when-solving-the-schrodinger-equation-with-ordinarydiffeq/125534/22 "2025-02-05T23:11:17Z")

</div>

> [@ChrisRackauckas](#):
>
> It would be weird if bigger ODEs just naturally have smaller dt even if it’s the same equation just repeated a few times.

Yeah, makes total sense. The surprising part was that `reltol` applies elementwise rather than to the overall norm.

---

<div class="post-metadata">

**Author:** ![BdeKoning](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bdekoning/32/214557_2.png) [@BdeKoning](https://discourse.julialang.org/u/BdeKoning)\
**Post date:** [March 11, 2025, 8:29am UTC](https://discourse.julialang.org/t/setting-abstol-and-reltol-when-solving-the-schrodinger-equation-with-ordinarydiffeq/125534/23 "2025-03-11T08:29:42Z")

</div>

Just for my own understanding if nothing else: would it be possible to explicitly enforce the solution norm with a DAE (or even eliminate part of the state)?

---

<div class="post-metadata">

**Author:** ![cwittens](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cwittens/32/213036_2.png) [@cwittens](https://discourse.julialang.org/u/cwittens)\
**Post date:** [March 11, 2025, 12:25pm UTC](https://discourse.julialang.org/t/setting-abstol-and-reltol-when-solving-the-schrodinger-equation-with-ordinarydiffeq/125534/24 "2025-03-11T12:25:08Z")

</div>

If I am not mistaken, it would be technically possible, but in practice you shouldn’t do it, because if you are already using an implicit solver to solve the DAE, you can just use an implicit method which preserves quadratic invariants (like the norm), which would be cheaper and simpler to implement.

---

<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:** [March 11, 2025, 1:52pm UTC](https://discourse.julialang.org/t/setting-abstol-and-reltol-when-solving-the-schrodinger-equation-with-ordinarydiffeq/125534/25 "2025-03-11T13:52:14Z")

</div>

> [@cwittens](#):
>
> you can just use an implicit method which preserves quadratic invariants (like the norm)

I assume you meant explicit here.

---

<div class="post-metadata">

**Author:** ![cwittens](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cwittens/32/213036_2.png) [@cwittens](https://discourse.julialang.org/u/cwittens)\
**Post date:** [March 11, 2025, 3:55pm UTC](https://discourse.julialang.org/t/setting-abstol-and-reltol-when-solving-the-schrodinger-equation-with-ordinarydiffeq/125534/26 "2025-03-11T15:55:47Z")

</div>

If you want to preserve arbitrary quadratic invariants for any f, your RK method needs to fulfill \forall i,j: b\_i a\_{ij} + b\_j a\_{ji} = b\_i b\_j, (as e.g the Gauss–Legendre methods do). And if I am not completely mistaken, this condition implies an implicit method.  
But of course in many cases for special f you can also find explicit methods preserving the norm.

---

<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:** [March 11, 2025, 4:27pm UTC](https://discourse.julialang.org/t/setting-abstol-and-reltol-when-solving-the-schrodinger-equation-with-ordinarydiffeq/125534/27 "2025-03-11T16:27:48Z")

</div>

A few things.

First of all, Gauss-Legendre methods are not L-stable and do not achieve higher than stage order 1, since that’s impossible with symplectic methods, so they are not suitable for most DAEs. If you have a DAE you almost certainly need to use some other kind of method, like Radau or BDF methods. In that case, just adding the norm condition is a good way to make sure it’s satisfied and you won’t be paying pretty much any extra cost. You do need to be careful with how it’s written though to ensure the norm condition is index 1, likely this will need a dummy derivative index reduction treatment, see [Dummy Derivatives and Reordering of Equations · ModelingToolkit Course](https://sciml.github.io/ModelingToolkitCourse/dev/lectures/lecture8/) for copious amount of detail on this. But the good thing is, ModelingToolkit automates this transformation so I’d just suggest MTK for these kinds of problems where you have a DAE and also want other conditions.

Second of all, it should be mentioned that implicit symplectic integrators only preserve the norm when the nonlinear solver has converged. So explicit symplectic integrators, such as verlet methods, tend to show this norm preservation behavior of quadratic invariants much better than Gauss-Legendre methods. That isn’t to say it won’t, it’s just that you’ll see a numerical drift over time which is linear based on the tolerance of the nonlinear solver that is used.

But now in their favor, third of all implicit symplectic integrators, since they are not stable enough for stiff equations, that also means you don’t need to use Newton methods with them and can just use function iteration / predictor corrector forms. So while they are technically implicit their implementation still looks very explicit and avoids any linear algebra.

So in total:

1. If you can use an explicit symplectic integrator, for example a non-stiff Hamiltonian, use it. That’s just going to be the best at preserving norm over long integrations for the amount of work you’re putting in, and the work per step is small.
2. If you have a non-stiff equation and specifically want to better enforce quadratic invariants, consider using a Gauss-Legendre method.
3. If you already have a DAE or if your equation is stiff, add that conservation law as a constraint and use ModelingToolkit’s `structural_simplify` to generate a stable index-1 form to solve with Rodas, Radau, or BDF methods.

And just putting it here since people don’t know about it, SciML does have a really good 16th order Gauss-Legendre implementation for the case (2):

> **[GitHub - SciML/IRKGaussLegendre.jl: Implicit Runge-Kutta Gauss-Legendre 16th order...](https://github.com/SciML/IRKGaussLegendre.jl)**
>
> Implicit Runge-Kutta Gauss-Legendre 16th order (Julia)

---

<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:10pm UTC](https://discourse.julialang.org/t/setting-abstol-and-reltol-when-solving-the-schrodinger-equation-with-ordinarydiffeq/125534/28 "2025-08-05T21:10:57Z")

</div>

Note I’ve added this example to a new package I’m working on to help with abstol and reltol settings: [bradcarman/ModelingToolkitTolerances.jl: A package to help determine `abstol` and `reltol` when using ModelingToolkit.jl and DifferentialEquations.jl](https://github.com/bradcarman/ModelingToolkitTolerances.jl)

The plot below shows the analysis this package can provide for this problem. The red arrow is pointing to the residual that the default tolerances give, which is fairly large. To get a residual closer to zero, the `reltol` can be adjusted to 1e-6 (or lower), as shown by the green arrow.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/d/8/d87d2372d8003b60619c52fe3f3f553585c84977.png)

[Previous page](https://discourse.julialang.org/t/setting-abstol-and-reltol-when-solving-the-schrodinger-equation-with-ordinarydiffeq/125534.md?page=1)