Instability after reducing number of equations in DifferentialEquations.jl

General Usage
1

I have a large set of differential equations (O(10000)) which I try to solve using Tsit5() and DifferentialEquations.jl.

sol = solve(
prob,
Tsit5(),
dtmax = 0.1,
maxiters = 100000,
progress = true,
saveat = 0.01,
)

I found a symmetry in the system, which allowed me to reduce the amount of equations. However, after some time, the symmetrized set of equations becomes unstable (see screenshot). Before that the dynamics looks the same as in the non-symmetrized case.

I do not understand this behavior - is the redundancy making the equations more stable?

I tried choosing alg_hints = :auto to automatically choose a solver, but the execution time became very large and the instabilities remained.

comparison
comparison700×500 19.9 KB

2

You can try the most accurate solvers(I like radauIIA5 but rosenbrock are also great) if that fails then you can try Sundial Fortran ones and if this fail too, it’s the model that may have something weird going one.

3

Did you try to set abstol and reltol to a low value, e.g. 1e-6?

4

Thanks for the quick reply! My equations are complex, and the solvers you suggested only work with real equations as far as I can see…
Do you think it is worthwile to restructure my code such that I propagate real and imaginary part separately?

5

Thanks for the quick reply! I set them 1e-6 but the instability is still there…

6

how is Rodas5(autodiff=false) doing ? maybe later on for performance you can split depends on the solver you will need

7

It gets stuck at the beginning… Maybe there are too many equations?

8

Can you share the equations? This will be impossible to debug at a distance.

9

Hi sorry for the late reply I was on a Christmas holiday!

I have over 10 000 equations and quite a huge code base, so it is not easy to create a minimal working example… Are there no diagnostics that I can do?

10

It seems to be a pretty problem specific issue, nonetheless a few questions that could maybe lead you to the solution:

  • What do you mean by unstable? Does the solver terminate with an error or is the result merely different?
  • Is there any way to reduce the number of equations without using symmetries? I.e. a numerical parameter, problem size etc that you can tune to make checking the problem a bit simpler?
  • If you solve the unsymmetrized ODE, does your solution maintain the symmetry? I.e. at each time step you could check this explicitly: If the symmetry is not fully maintained, does this happen immediately at the beginning or is there perhaps a gradual increase in the error?
  • If the symmetry is not kept, this could be due to numerical errors, in computing the derivative. An error of about 1e-14 would likely be simple floating point accuracy which may be harmless. If the derivative you compute strictly obeys symmetries, then I would not expect the choice of the ODE integrator to matter.
  • If the symmetry is maintained also in your numerics, you can try to gradually move from the non-symmetric to the symmetric solution:
  1. In the nonsymmetric part, compute the neccessary derivatives. Then set the symmetry equivalent ones explicitly to be equal to the nonsymmetric parts. Note that you keep the size of the ODE state array the same. If your numerics keep the symmetry this should still produce the same result.
  2. If 1 is successful, you could try actually reducing the number of equations.
11

This still needs code to do any debugging.