Help with numerical accuracy when solving a system of 300 ODE's using DifferentialEquations.jl

General Usage
,
1

Hello,

I have difficulties to accurately solve my stiff ODE system which contains 300 nonlinear differential equations representing the motion (position, velocity and state) of 100 inter-connected blocks. A basic working version of the julia code is provided in git/test-ode.

I have tried different solver algorithms and combinations of relative and absolute tolerances. I have represented below the maximal differences in position over time for all blocks for different tolerances compared to a reference tolerance of abstol=1e-12 and reltol=1e-12 using the Rodas5() solver. The difference can be reduced at lower tolerances but nevertheless increases over time and no tolerance choice can effectively force convergence among the results.

EqSim_Benchmark
EqSim_Benchmark1129Ă—916 92.3 KB

I am able to run simulations over the desired timespan of tspan=(0.0,10000.0) without any ODE instability warnings/errors using my reference tolerance + solver algorithm. The simulation reproduces the expected stick & slip motion (Blocks are initially stuck until the driving forces exceed the friction forces and blocks begin to slip, over multiple repetitions) but the details (e.g. maximum block velocity, motion start) significantly depend on the tolerances. Given the larger timespans, small initial differences in the motion become significant and alter the system behaviour.

I have read through (Stiffness?) Problem with Differential Equation System and wonder if there are nevertheless unnoted instabilities (parameter depend) related to my ODE system (e.g. during the “slip phase”)?

Considering efficiency, I determined that my implementation of the ode system is faster using the @ ode_def macro provided by ParameterizedFunctions.jl than using the matrix du=A*u+f(u) formulation.

I am thankful for any suggestions and help of how I can reach numerical stability and/or improve the efficiency when solving my ode system? I am available for further clarifications.

2

Is the system chaotic? Some systems just have the behavior that differences will grow exponentially over time and thus over longer timespans their result is simply sensitive. The tolerances are tracking the solver’s accuracy quite well here so I’m not sure it can really do much better: at this point it’s starting to run into floating point accuracy issues (especially since it has to solve a 300x300 linear system). So it could just be a property of the system, like the Lorenz system.

1 Like