Comparison of integrators DifferentialEquations.jl

Specific Domains Modelling & Simulations
1

I get two different trajectoriesf or identical parametres, time, and initial condition.
for comparison i using integrator from Julia and python. I have used several integrators from Julia and python and get completely different results. Even when I imported an integrator with scipy to Julia, the result did not change.

Code from python

Gamma = lambda x_i: 1 / (1 + np.exp(-lambd * (x_i - Theta)))
def two_HRs(t, X, k1, k2):
    x1, y1, z1, x2, y2, z2  = X
    dx1 = y1 - a*x1**3 + b*x1**2 -z1 + I - k1*(x1 - v_s) * Gamma(x2)
    dy1 = c - d*x1**2 - y1
    dz1 = r*(s*(x1 - x_r) - z1)
    
    dx2 = y2 - a*x2**3 + b*x2**2 -z2 + I - k2*(x2 - v_s) * Gamma(x1)
    dy2 = c - d*x2**2 - y2
    dz2 = r*(s*(x2 - x_r) - z2)
    return [dx1, dy1, dz1, dx2, dy2, dz2]

a = 1
b = 3
c = 1
d = 5
x_r = -1.6
r = 0.01
s = 5
I = 4
v_s = 2
lambd = 10
Theta = -0.25

k1 = -0.17
k2 = -0.17

initials = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1]

sol = solve_ivp(two_HRs, [0, 100000], initials, rtol = 1e-11, atol = 1e-11, args = (-0.17, -0.17), method="LSODA")
ts = sol.t
x1s, y1s, z1s, x2s, y2s, z2s = sol.y

i try use RK45, LSODA.
Code from Julia:

function sigma(x)
    return 1.0 / ( 1.0 + exp( -10.0 * ( x  - ( - 0.25 ) ) ) )
end

function HR(u, p, t)
        
    a, b, c, d, s, xr, r,  I, vs, k1, k2, el_link  = p
    x1, y1, z1, x2, y2, z2 = u
    
    du1 = y1 + b * x1 ^ 2 - a * x1 ^3 - z1 + I - k1 * ( x1 - vs ) * sigma(x2) + el_link * ( x2 - x1 )
    du2 = c - d * x1 ^2 - y1
    du3 = r * ( s * ( x1 - xr ) - z1 )
    
    du4 = y2 + b * x2 ^ 2 - a * x2 ^3 - z2 + I - k2 * ( x2 - vs ) * sigma(x1) + el_link * ( x1 - x2 )
    du5 = c - d * x2 ^2 - y2
    du6 = r * ( s * ( x2 - xr ) - z2 )
    return SVector(du1, du2, du3,
                    du4, du5, du6)
end
a = 1.0; b = 3.0; c = 1.0; d = 5.0
xr = -1.6; r = 0.01; s = 5.0; I = 4.0; xv = 2.0
k12 = -0.17
k1= k12; k2 = k12

k = 0.0

condition = SA [0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
tspan = (0.0, 100000.0)
p = SA[a, b, c, d, s, xr, r, I, xv, k1, k2, k];

probsa = ODEProblem(HR, condition, tspan, p)
sol = solve(probsa, AutoVern9(Rodas5()), abstol = 1e-11, reltol = 1e-11, dense=false, maxiters = 50000000);

Result from python

image
image2268×588 57 KB

and Julia

image
image1000×400 55.4 KB

in time series i using sum x1 + x2.
Also i try use dense = false and dense = true

2

Is the system chaotic?

3

Yes, system with such parameters has multistability, that is, there is a regular chaotic attractor.

Using completely identical initial conditions Julia show regular attractor and python chaotic.
I also calculated the Lyapunov spectrum using two different algorithms jitcode and DynamicalSystems.jl, and the result indicates hyperhaos.

In addition, I identify the attractors and their basins. In result i have two attractor this regular and chaotic.

4

python integrators have too low accuracy and therefore fly into chaotic attractor?

5

No, this is just literally a property of chaos. Changes of 1e-16 will give O(1) differences in the trajectory. The difference of using muladd vs x*y + z is then enough to, over time after a Lyopunov time, give O(1) differences. You can only get exactly the same result out of chaotic systems if you have exactly, floating point exact, the same implementation.

Take a look at https://diffeqcallbacks.sciml.ai/stable/uncertainty_quantification/ to explore the uncertainty of such systems.

Everything that approximates a chaotic attractor is incorrect. You only have the shadowing lemma: the fact that there’s a trajectory with an epsilon different initial condition on the attractor that you happen to have pulled, but it can be O(1) different from the real trajectory. Numerical approximations to chaotic systems are only trustworthy in a probabilistic sense.

Here’s a nice post that shows some of the details:

1 Like
6

Thank you. Is it possible to somehow accurately determine and say that this trajectory with an initial condition is a regular mode? it exists only in one single case if all the initial conditions are the same

7

Can you define the regularity (i.e. the rules) of the regular mode?

8

The regular regime can be determined by a period or lyapunov spectrum or largest lyapunov exponent, also power spectrum

9

So you are saying a regular mode in your case means a periodic orbit

10

Yes, in articles people often write a regular mode, either a periodic orbit or a periodic trajectory

11

I would personally use ApproxFun.jl combined with an optimization package to find periodic solutions. (Of course it’s quite hard to show non-existence of a periodic solution this way.)

EDIT: Especially since strange attractors can have unstable periodic orbits which would show up in this approach.

12

Thank you, I’ll try to use this packages. I started writing code on Julia relatively recently and don’t know many packages

13

The gist of this approach is explained in section 2.3.4 of this article (it’s definitely not original; just a methods section)

14

Thank you. Do you happen to have a code example for using this with a continuous system?

15

Yes, have a look at this—it’s doesn’t have much documentation, unfortunately. And since it’s a few years old the code might not work completely with current versions of the packages.

1 Like
16

Thank you