Stiff VS Non-stiff methods: Differing solutions

[Sh_166]

Apologies if there has been a trivial mistake or otherwise in what follows.

MWE:

using Symbolics, Plots, DifferentialEquations, ModelingToolkit, 
LaTeXStrings, Measures, DynamicalSystems

function forced_osc!(du, u, p, t) #the system
    ω = p[1]; F = p[2]; 
    du[1] = u[2]
    du[2] =  - sin(u[1])*ω^2 + F*cos(ω*t)*ω^2 
    return nothing
end

u0 = [0, 1e-4]
diffeq = (alg = Rodas4P(), abstol = 1e-8, reltol = 1e-8, maxiters = 1e8) #Solver specifications

forcedoscillator = ContinuousDynamicalSystem(forced_osc!, u0, [1,0.5]; diffeq)
ψ, t = trajectory(forcedoscillator6, 800; Δt = 0.1) 

PS = scatter(ψ6[:, 1], ψ6[:, 2], lw = 0, markersize = 0.5, size = (700, 400), markerstrokewidth = 0, framestyle = :box, label = L"F = 0.5", legendposition = :topright)

Now, this runs perfectly and so and so. The issue is that the solution differs based on what algorithm I use.

Examples:

With Rodas4P():

With RK4():

With Tsit5() (default):

With AutoVern7(Rodas4()) (recommended in the docs for lower tolerances):

All of these have differing long-time behaviours - is this expected instability associated with the solvers? I started with a stiff solver because the stiffness was unknown (and I had used this Rodas4P() combination in a previous problem) - as far as I know, stiff solvers can be used for non-stiff problems but the efficiency drops. Since my system wasn’t taking much time, I didn’t reflect on it much.

Is there a reason/solution for this? I noted that the Rodas4P() solution goes more towards the rest if the tolerances are further lowered to 1e-9. Since the rest of the solvers branch towards the right, I am guessing the Rodas4P() solution is the least credible here?

I have also verified that this same thing also happens if one only uses these methods in DifferentialEquations.jl only (i.e., without using DynamicalSystems.jl).

Any help will be appreciated!

Note: It does seem like an inbuilt solver error - lowering the tolerances to 1e-9 gives proper match between these solutions for around tspan = (0, 850).

[Datseris]

@Sh_166

du[2] = - sin(u[1])*ω^2 + F*cos(ω*t)*ω^2

While I have no real expertise with stiff solvers, I would argue that what you see is completely normal. Your dynamical system seems to be chaotic, non-autonomous, multistable, and metastable. This type of dynamical system will be very sensitive to the exact solver you use. In fact, even your operating system will have impact on the exact trajectories you will see. This isn’t a problem, this is numerical chaos.

What you want to make sure of is that dynamic invariants such as lyapunov exponents do not depend on the solver much (provided you use trajectories that go to the same attractor). If they do, you need to lower the tolerances more.

My hunch is further corroborated by what you stated:

@Sh_166

lowering the tolerances to 1e-9 gives proper match between these solutions for around tspan = (0, 850).

[Datseris]

May I also recommend that you study the PoincareMap of this system instead of the continuous time trajectories? I believe it will simplify a lot your analysis. The poincare map is “equivalent” due to the time-periodic forcing. (see chapter 9 of Nonlinear Dynanimics for details: Nonlinear Dynamics: A Concise Introduction Interlaced with Code | SpringerLink)

[Sh_166]

@Datseris I think so too - in fact, a claim we were trying to make/verify is that this system is susceptible to chaos for some range of the drive amplitude.

I observed that for F = 0.25 to F = 0.33, we see perfect overlap between all the solvers for an integration up to t = 10,000 secs. For F = 0.34, the overlap vanishes for t \sim 135 secs. It is not possible, of course, to check every single drive value, but I have done checks like this before as well (usually for small drive values) and I have never observed a disagreement between the solvers/packages for F<0.34. However, I am not an expert in either chaos theory or numerical methods.

Now, the Lyapunov exponents are somewhat important to our current study - I was using the inbuilt lyapunov() function to get the maximum Lyapunov exponent. Would there be an issue here? I was using a total time of integration of 800 secs - so, I presume I would need to lower the tolerances for this range?

We were actually studying the Poincare Maps of this system as well - for various drive values. The code was somewhat same, except that we had,

ψ, t = trajectory(forcedoscillator6, 25000; Δt = 2π) 

I believe that the time of integration and tolerances have to be modified in this case as well for relevant forcing values.

Thanks a lot for the prompt response!

[Datseris]

just use low enough tolerances so that you always get the same lyapunov exponent (for the same parameter configuration). I think this is already good enough.

Note: for chaotic systems I never go above 1e-9 tolerances, and I typically use high order solvers like Vern9.

[Sh_166]

Thanks for the heads up (and also, the very nice reference)!

[Tarny_GG_Channie]

Hmm? What’s happening? I thought stiff equation solvers were reliable at solving chaotic systems.

[Oscar_Smith]

This isn’t the case. The definition of chaotic systems is that after some time interval t_i epsilon differences in your initial conditions will result in solutions that differ by O(1). No method on a computer can “solve” this exactly because any rounding error, no matter how small will result in a completely different solution.

The one thing that saves you is the shadowing lemma, but that applies equally to stiff and non-stiff solvers.

[vettert]

As mentioned before chaotic systems like this will never give you the same trajectories when changing solvers/OS/hardware.

This kind of thing appears in many areas of science. For example, when performing molecular dynamics simulations, it is well known within the community that the position of a specific molecule over time will differ between repetitions of the same run unless extreme care is taken (that includes things like running the simulation on the same node in an HPC cluster (hardware!), same RNG state, etc etc.).

[Datseris]

@Oscar_Smith

The one thing that saves you is the shadowing lemma , but that applies equally to stiff and non-stiff solvers.

Unfortunately the shadowing lemma does not apply to the overwhelming majority of dynamical systems people use in practice to model the real world, as they are not hyperbolic (which is a rather special condition).