This is taught in the intro class in basic numerics.
Methods are considered more robust and practical if they have larger stability regions. Von-Neumann is the OG analyst for PDEs
> Those kinds of models aren’t really realistic? Energy goes to infinity, angular momentum goes to infinity, the chemical concentration goes to infinity: whatever you’re modeling just goes crazy! If you’re in this scenario, then your model is probably wrong. Or if the model isn’t wrong, the numerical methods aren’t very good anyways.
Uh... that's not how it works. The whole point of making a model and using a simulator is you don't already know what's going to happen. If your system is going to blow up, your model needs to tell you that. If it does, that's good, because it caught your mistake. And if it doesn't, it's obviously wrong, possibly in a deadly way.
This aligns with what I've observed in computational physics.
Trying to handle instantaneous constraints and propagating modes with a single timestepper is often suboptimal.
I developed a framework that systematically separates these components: using direct elliptic solvers for constraints and explicit methods for flux evolution. The resulting algorithms are both more stable and more efficient than unified implicit approaches.
The key insight is that many systems (EM, GR, fluids, quantum mechanics) share the same pattern:
- An elliptic constraint equation (solved directly, not timestepped)
- A continuity law for charge/mass/probability flux
- Wave-like degrees of freedom (handled with explicit methods)
Witgh this structure, you can avoid the stiffness issues entirely rather than trying to power through them with implicit methods.
In practical terms it’s not unusual to reset the integrator when something instantaneous happens. When i did magnetic research, an application of instantaneous field for e.g. usually required this because otherwise the adaptive integrator spends a lot of time reducing the time step size
The orbital example where BDF loses momentum is really about the difference between a second-order method (BDF2) and a fourth-order method (RK), rather than explicit vs implicit (but: no method with order > 2 can be A-stable; since the whole point of implict methods is to achieve stability, the higher order BDF formulas are relatively niche).
There are whole families of _symplectic_ integrators that conserve physical quantities and are much more suitable for this sort of problem than either option discussed. Even a low-order symplectic method will conserve momentum on an example like this.
I always see RK(Fe) mentioned but almost never Bulirsch-Stoer (https://en.wikipedia.org/wiki/Bulirsch%E2%80%93Stoer_algorit...). Some years ago, I was looking to do high precision orbital motion (involving several bodies) simulations, and came across this in Numerical Recipes for C. I don't pretend to understand how B-S works, but it produced amazing precision with adaptive step size adjustment. For highly chaotic systems, it held out far longer than RK.
I was expecting a mention of symplectic ODE solvers, although perhaps that was beyond the scope of the blog post. For Hamiltonian ODEs, you can design methods that explicitly preserve energy, outperforming more generic methods.
In the vast majority of cases, people using ODE solvers are working with physical systems, where there likely aren't more than a dozen digits of accuracy in any of the parameters. So you set some fixed level of accuracy for the output, and seek a method that can reliably attain that accuracy at the lowest cost. (Either that, or the system is too chaotic, in which case you're told not to bother since the output would be meaningless anyway.)
Something I've had trouble finding any information about is, what if we have a mathematical system, where the initial configuration is known to infinite precision, and we want to solve the ODE to an arbitrary precision (say, 100 digits, or 1000 digits)? All the usual stepwise methods would seemingly require a very high order to avoid a galactic step count, but initializing all the coefficients, etc., would be a struggle in itself.
Is this the best that can be done, or are there less-usual methods that would win out at extremely high accuracy levels?
8 comments
[ 2.5 ms ] story [ 31.3 ms ] threadUh... that's not how it works. The whole point of making a model and using a simulator is you don't already know what's going to happen. If your system is going to blow up, your model needs to tell you that. If it does, that's good, because it caught your mistake. And if it doesn't, it's obviously wrong, possibly in a deadly way.
Trying to handle instantaneous constraints and propagating modes with a single timestepper is often suboptimal.
I developed a framework that systematically separates these components: using direct elliptic solvers for constraints and explicit methods for flux evolution. The resulting algorithms are both more stable and more efficient than unified implicit approaches.
The key insight is that many systems (EM, GR, fluids, quantum mechanics) share the same pattern: - An elliptic constraint equation (solved directly, not timestepped) - A continuity law for charge/mass/probability flux - Wave-like degrees of freedom (handled with explicit methods)
Witgh this structure, you can avoid the stiffness issues entirely rather than trying to power through them with implicit methods.
Paper: https://zenodo.org/records/16968225
There are whole families of _symplectic_ integrators that conserve physical quantities and are much more suitable for this sort of problem than either option discussed. Even a low-order symplectic method will conserve momentum on an example like this.
Something I've had trouble finding any information about is, what if we have a mathematical system, where the initial configuration is known to infinite precision, and we want to solve the ODE to an arbitrary precision (say, 100 digits, or 1000 digits)? All the usual stepwise methods would seemingly require a very high order to avoid a galactic step count, but initializing all the coefficients, etc., would be a struggle in itself.
Is this the best that can be done, or are there less-usual methods that would win out at extremely high accuracy levels?