JuliaDynamics / ChaosTools.jl

Tools for the exploration of chaos and nonlinear dynamics
https://juliadynamics.github.io/DynamicalSystemsDocs.jl/chaostools/stable/
MIT License
188 stars 36 forks source link

callbacks on ds #324

Closed gacarita closed 7 months ago

gacarita commented 7 months ago

Hi there,

I'm facing some trouble with the callback routines when running chaos tools. So far running my program directly with the DifferentialEquations library and got no errors. The energy is maintained as expected.

image

However when passing everything to the Chaostools library the energy deviates every time that the callback is called out

image

My problem is quite simple I need to solve elastic collisions between ellipsoids. I'm doing something wrong or it could be a bug ?

Below follows a piece of my code, I'm using the same parameters that I've used for the normal solution with DifferentialEquations. In the end I want to use a chaotic indicator such as GALI.


function condition(u, t, integrator) # Event when condition(u,t,integrator) == 0

    x1 = a1*cos(u[3])
    y1 = b1*sin(u[3])

    x2 = a2*cos(u[4])
    y2 = b2*sin(u[4])

    sep1 = sqrt(x1^2 + y1^2)
    sep2 = sqrt(x2^2 + y2^2)

    u[1] - sep1 - sep2

end

function affect!(integrator)
    integrator.u[5] = -1.0*(0.5*integrator.u[5] + 0.5*integrator.u[5])
end

cb = ContinuousCallback(condition, affect!;reltol=1e-3)

diffeq = (abstol=1e-12, retol=1e-12, alg = Vern9(), callback=cb)

ds = CoupledODEs(absolute_motion, y_0; diffeq)

y_0 = [1.0,0.0,0.0,0.0,0.2,1.0,0.0,0.0]# test vector

tr = trajectory(ds, 100.0, y_0;)[1]
Datseris commented 7 months ago

Hi, thanks for reporting this. First of all, can you provide a minimal working example so I can actually run the code and be able to examine it?

Now for the callback. I think the trajectory function is not really meant to be used with callbacks but rest of the library should respect them normally. We can test this. You can do:

energies = map(1:1000) do i # iterate over i steps
step!(ds)
u = current_state(ds)
e = energy(u)
end

where energy the function that computies the energy. Can you please do that and see if the vector energies stays constant?

gacarita commented 7 months ago

Hi Datseris,

I could share a minimal working example, but I prefer to do it by email instead of leaving it online, as this is related to my research. If it is ok for you I can sent a example by email.

Thank you for the reply. Actually, the trajectory function does not work properly. I have followed what you have mentioned, and I was able to reproduce the results from the differential equations, and everything is working properly. Now I wonder if it is possible to use the GALI function in this case.

When I was running tests with the trajectory function, I noticed that changing between diffeq = (alg = Vern9(), reltol = 1e-12, abstol=1e-12, callback=cb) and diffeq = (reltol = 1e-12, abstol=1e-12, alg = Vern9(), callback=cb) caused differences in the results, especially in energy conservation.

Datseris commented 7 months ago

I have followed what you have mentioned, and I was able to reproduce the results from the differential equations, and everything is working properly. Now I wonder if it is possible to use the GALI function in this case.

GALI will work fine. Everything in DynamicalSystems.jl will work fine because everyting uses only step!(ds).

trajectory is the only special function that doesn't use step!. Hence, you do not need to share your code for now.

When I was running tests with the trajectory function, I noticed that changing between diffeq = (alg = Vern9(), reltol = 1e-12, abstol=1e-12, callback=cb) and diffeq = (reltol = 1e-12, abstol=1e-12, alg = Vern9(), callback=cb) caused differences in the results, especially in energy conservation.

I am confused here, is there any typo in this? because both diffeq are identical.

Datseris commented 7 months ago

Oh, actually, you have to be careful. The dynamical system used by GALI is a TangentDynamicalSystem. The state of it is different, and you may need to adjust the callback accordingly. Try initializing a TangentDYnamicalSystem with all deviation vectors being the 0 vector. Then evolve this TangentDYnamicalSystem it with the callback with the same for ... step!() energy() code as above, and make sure that the callback conserves energy.

gacarita commented 7 months ago

Hi there,

Thank you.

Regarding the diffeq, looks like even changing the args positions caused an difference in the results when I was using trajectory. Since I'm not using trajectory anymore I can't tell you if this is correct or not. I was doing many tests and might this was the cause to my wrong value observation.

So far, I've tested your suggestion to evolve the TangentDynamicalSystem with Q0=zeros(8,8) using step! and determining the energy using the ds with the callbacks and the energy conserves, at this point I would say that everything is ok and I may use GALI without problems.

Below follows a piece of the code just in case that I'm wrong.

ds = CoupledODEs(absolute_motion, y_0;diffeq)

tands = TangentDynamicalSystem(ds;Q0=zeros(8,8))

checked with

current_deviations(tands)

returned zero matrix(8,8). Evolved over time and the energy was conserved.

Datseris commented 7 months ago

Perfect then I will assume tat everything is okay and closing this. trajectory isn't planned to work with callbacks at the moment so we'll leave it at that. I added a note to the trajectory docstring about this.