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.
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.