Stopping Simulation - Question

[gacarita]

Hello there,

My name is Gabriel and I work in astronomy and space engineering and am new to julia. My tests involve collisions and escapes, so I need to check each step. I wonder if there is any way to stop numerical integration without using differential equation compatibility callbacks ?

I would also like to know if there is any way to control the number of points in the output?

[lbenet]

Hi Gabriel! Thanks for asking and using our package.

Regarding the question about checking some condition (which I guess you define from the distance of two objects) at each step of the integration, there is no specific way of doing that. Yet, there are two ways I can think of doing that:

• You can define a function g that defines the condition you are interested, which is analogous to the condition defining a Poincare surface of section. Have a look here in the documentation, where an example for a surface of section is defined. The interface is a bit more general, so I think it may be the solution you are asking for.
• Based on one of the methods of taylorinteg, you could define your own function which includes checking the condition after each step of the integration (after taylorstep!) but before updating all needed quantities for the next step of the integration.

Regarding the second question, whether we have a way of controlling the number of points of the output, yes we do. We have methods that return all time steps, and you can fix the maximum number of integration steps through the keyword argument maxsteps, and also have methods that return the quantities at given time lapses, by using ranges to define the integration limits. You can also pass a vector of times where you want the solution to be evaluated. Have a look on the docstrings of taylorinteg.

I hope these comments are helpful!

[gacarita]

Thank you Luis, totally helpful.

I've runned few tests comparing Vern9() with TaylorIntregration and for very long tspans TaylorIntegration sounds better since the energy was conserved better than Vern9() for a million years and 10 million years and the running time was faster.

I'm trying to implement a simple problem like Planar CR3BP and see how poincare surfaces was going on, my sections are doing well i'm obtaining the expected results and the code seems fine. However i still faces some implementation difficulties.

My problem is :

1) I would like to add a new routine running with the poincare (g function below) to checkout for collision and escapes. Instead of just entering g I should enter g and esc!(escape function) in the format [g,esc!] ?

example: tv_i, xv_i, tvS_i, xvS_i, gvS_i = taylorinteg(f, [g,esc!], u0, 0.0, 15000.0, 25, 1e-17,p ,maxsteps=90000,parse_eqs=false);

2) To stop the simulation at exact collision or escape time what should i return in esc! function ? (i imagine should be something like [true, u] but i'm no sure how to stop the simulation.

My routines to call out the simulation are defined below.

# poincare section routine
function g(dx, x, p, t)
    px_ = constant_term(x[4])

    # crossing variable >0 or <0
    if px_ > zero(px_)
        return (true, x[3])
    else
        #otherwise, discard the crossing
        return (false, x[3])
    end
end

@inline @inbounds function pmap(f,g,u0,p)
        tv_i, xv_i, tvS_i, xvS_i, gvS_i = taylorinteg(f, g, u0, 0.0, 15000.0,
                25, 1e-17,p ,maxsteps=90000,parse_eqs=false);
        tvSv = vcat(0.0, tvS_i)
        xvSv = vcat(transpose(u0), xvS_i)
        gvSv = vcat(0.0, gvS_i)
        return tvSv,xvSv[:,1],xvSv[:,2],xvSv[:,3],xvSv[:,4]
end

function esc!(du, u, p, t)
    local mu2 = p
        r=sqrt(u[1]^2+u[2]^2)
    # escape
.
.
.
    if r > 10
        return (what should i return here to stop the numerical integration and return the escape pos,vel and time ?)
    end
end

Thank you

[lbenet]

I've runned few tests comparing Vern9() with TaylorIntregration and for very long tspans TaylorIntegration sounds better since the energy was conserved better than Vern9() for a million years and 10 million years and the running time was faster.

Regarding this let me first ask, what is the version of the package you are using? The reason I ask this is because the last version (v0.9.0) should be faster than the previous one (v0.8.x), if you use @taylorize. (I hope you are using 0.8.x or you are not using yet @taylorize, so when you will have an even faster implementation!) Note that using @taylorize imposes a bunch of constrains in the way you code the differential equations...

I'm trying to implement a simple problem like Planar CR3BP and see how poincare surfaces was going on, my sections are doing well i'm obtaining the expected results and the code seems fine.

This problem is actually close to my heart; see this and this papers :-) Once I have promoted my own (old!) papers (I hope you knew them already), let me say that on that time I chose quite uncommon coordinates for the Poincare section, precisely because I was interested in escapes. Collisions eventually appeared, depending on the Jacobi constant, and they are difficult to handle...

However i still faces some implementation difficulties.

What I would do is to write another function based on the method of the method taylorinteg you are using, as I proposed before.

The idea is to introduce the "escape switch" (your esc! function) once the integration step is finished, and perhaps even after the integration coordinates have been updated, so esc! uses the last calculated values. If you read the code of the taylorinteg method you are using, the check involving the g function (for the Poincare section) involves the function findroot!, then a bunch of things are updated, and finally a check related to the number of steps is performed here; the break that appears within that if-block exits the while loop, which is the actual integration cycle. I think you could either check right there if the particle has reached the escaping region, e.g., if (nsteps > maxsteps) || esc!(...), or if you prefer with its own if-block. Then, esc! could simply return true if that is the case and then breaking the while loop, or false, so the integration proceeds as before.

If you follow this scheme, you don't need to change anything related to the returned output, though you certainly may, because you can always evaluate the esc! function, at the final integration point obtained, and know if the particle reached the escaping region, or not.

I hope these comments are helpful!

[gacarita]

Gotcha Luis.

Yes i'm using the last version and the @taylorize. I'm really enjoing the precision of this numerical integrator, i'm probably gonna use it on my future works.

Unfortunately, I didn't know yours published works. I've have finished my masters 3months ago and i worked a half of my thesis on the PCR3BP with surface sections in retrograde orbits. I'm an enthusiast about anything related with chaotic/periodic orbits analysis methods and tools. I'm probably gonna read your work with pleasure since is a very related area to me. I've started my PhD and i'm gonna search something related with chaos 3BP and perturbations somehow. Handling with close-encounters and collisions was the problem i was facing using anothers languages and numerical integrators. I'm very new to Julia (i have some experience in python) and i'm trying to understand and starting to make my tools to start my work.

I've forked the files and coppied to my computer explicit everything as functions and made some changes taylorinteg functions to always check for my collision/escape conditions following the steps you shared. It was trick at first because i didn't know about the evalute! function but now seems it's working now i just need to keep update the library when new versions come out. Now i probably may use it to any Differential Equations i have.

A suggestion is an implementation of a generic function to break simulation when this conditions becomes true.

Thank you for the attention. Feel free to close this topic.

[lbenet]

Thanks Gabriel for all comments and suggestions. I just opened #145 to not forget your suggestion.

Regarding the RC3BP, there is a lot of recent work carried out by E. Zottos and collaborators. Also, there is some work by (applied) mathematicians which may interest you, specially related to collisions. Collisions are numerically nasty, in the sense that the manifolds yielding the collision are parabolic; yet, the numerical instabilities are apparent. A nice problem all together!

[PerezHz]

I've runned few tests comparing Vern9() with TaylorIntregration and for very long tspans TaylorIntegration sounds better since the energy was conserved better than Vern9() for a million years and 10 million years and the running time was faster.

I'm a bit late to the party, but just wanted to mention that you can also use TaylorIntegration via the common interface with DifferentialEquations.jl, and that way you can use callbacks to stop an integration given a certain condition, both with Vern9 or TaylorMethod. Hope this helps.