# Stiff VS Non-stiff methods: Differing solutions
**URL:** https://discourse.julialang.org/t/stiff-vs-non-stiff-methods-differing-solutions/114664
**Category:** General Usage
**Tags:** package, differentialequation, dynamical-systems
**Created:** [May 24, 2024, 8:48am UTC](https://discourse.julialang.org/t/stiff-vs-non-stiff-methods-differing-solutions/114664 "2024-05-24T08:48:11Z")
**Posts on this page:** 10
**Page:** 1
### Author:  [@Sh\_166](https://discourse.julialang.org/u/Sh_166)
#### Post date: [May 24, 2024, 8:48am UTC](https://discourse.julialang.org/t/stiff-vs-non-stiff-methods-differing-solutions/114664/1 "2024-05-24T08:48:12Z")
Apologies if there has been a trivial mistake or otherwise in what follows.
MWE:
```julia
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)`.
---
### Author:  [@Datseris](https://discourse.julialang.org/u/Datseris)
#### Post date: [May 24, 2024, 9:20am UTC](https://discourse.julialang.org/t/stiff-vs-non-stiff-methods-differing-solutions/114664/2 "2024-05-24T09:20:42Z")
> [@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)`.
---
### Author:  [@Datseris](https://discourse.julialang.org/u/Datseris)
#### Post date: [May 24, 2024, 9:22am UTC](https://discourse.julialang.org/t/stiff-vs-non-stiff-methods-differing-solutions/114664/3 "2024-05-24T09:22:44Z")
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](https://link.springer.com/book/10.1007/978-3-030-91032-7))
---
### Author:  [@Sh\_166](https://discourse.julialang.org/u/Sh_166)
#### Post date: [May 24, 2024, 9:30am UTC](https://discourse.julialang.org/t/stiff-vs-non-stiff-methods-differing-solutions/114664/4 "2024-05-24T09:30:12Z")
@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,
```julia
ψ, 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!
---
### Author:  [@Datseris](https://discourse.julialang.org/u/Datseris)
#### Post date: [May 24, 2024, 9:39am UTC](https://discourse.julialang.org/t/stiff-vs-non-stiff-methods-differing-solutions/114664/5 "2024-05-24T09:39:31Z")
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.
---
### Author:  [@Sh\_166](https://discourse.julialang.org/u/Sh_166)
#### Post date: [May 24, 2024, 9:47am UTC](https://discourse.julialang.org/t/stiff-vs-non-stiff-methods-differing-solutions/114664/6 "2024-05-24T09:47:37Z")
Thanks for the heads up (and also, the very nice reference)!
---
### Author:  [@Tarny\_GG\_Channie](https://discourse.julialang.org/u/Tarny_GG_Channie)
#### Post date: [May 24, 2024, 1:34pm UTC](https://discourse.julialang.org/t/stiff-vs-non-stiff-methods-differing-solutions/114664/7 "2024-05-24T13:34:20Z")
Hmm? What’s happening? I thought stiff equation solvers were reliable at solving chaotic systems.
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [May 24, 2024, 1:50pm UTC](https://discourse.julialang.org/t/stiff-vs-non-stiff-methods-differing-solutions/114664/8 "2024-05-24T13:50:44Z")
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](https://en.wikipedia.org/wiki/Shadowing_lemma), but that applies equally to stiff and non-stiff solvers.
---
### Author:  [@vettert](https://discourse.julialang.org/u/vettert)
#### Post date: [May 24, 2024, 2:39pm UTC](https://discourse.julialang.org/t/stiff-vs-non-stiff-methods-differing-solutions/114664/9 "2024-05-24T14:39:12Z")
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.).
---
### Author:  [@Datseris](https://discourse.julialang.org/u/Datseris)
#### Post date: [May 24, 2024, 2:44pm UTC](https://discourse.julialang.org/t/stiff-vs-non-stiff-methods-differing-solutions/114664/10 "2024-05-24T14:44:25Z")
> [@Oscar\_Smith](#):
>
> The one thing that saves you is the [shadowing lemma](https://en.wikipedia.org/wiki/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).