Yes that is correct.
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.
Yeah, makes total sense. The surprising part was that reltol applies elementwise rather than to the overall norm.
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)?
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.
I assume you meant explicit here.
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.
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 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:
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):