hope you're all fine.
That's my first question, so I'm sorry if something's not right about it.
I'm studying numerical stability and chaoticity of some dynamical systems, more specifically, about the Circular Restricted Three Body Problem (CR3BP), defined by the set of 3 second order differential equations. After transforming those three second order differential equations in six first order differential equations as seen here i can finally finally work with them numerically using scipy's ODEINT. Here's an example of an orbit integrated for T = 2^10 with n = 2^18 points (np.linspace(1, 2^10, 2^18)) and here's a bit of my code for this, the main function to be integrated:
def func(init, t, mu):
#x0, y0, z0, vx0, vy0, vz0 = init
x0, y0, z0, vx0, vy0, vz0 = init[0], init[1], init[2], init[3], init[4], init[5]
dx1dt1 = vx0
dy1dt1 = vy0
dz1dt1 = vz0
dx2dt2 = 2*dy1dt1+d_omega_x(mu, x0, y0, z0)
dy2dt2 = -2*dx1dt1+d_omega_y(mu, x0, y0, z0)
dz2dt2 = d_omega_z(mu, x0, y0, z0)
return np.array([dx1dt1, dy1dt1, dz1dt1, dx2dt2, dy2dt2, dz2dt2])#, dx2dt2, dy2dt2, dz2dt2])
in which x, y, z, vx, vy, vz = (0.848, 0, 0, 0, 0.0423, 0) and mu = 0.01215.
Then comes the stabillity part. I'm using a chaos detection tool called Fast Lyapunov Indicator. It's basically defined by v'(t)=Df(x)v(t), where Df(x) is the Jacobian matrix of the system of equations (in this case, a 6x6 matrix) and v(t) a tangent vector to be evolved in time alongside the six original equations of the CR3BP, then I take the log10 of the norm of the six components of the integrated v(t) and pick the highest value.
One might note the 6 "auxiliary" vectors obtained from v'(t)=Df(x)v(t) depend on the original 6 equations (more specifically on the positions x, y, z of the particle), but the six original equations do not depend on the six new equations related to the tangent mapping defined by v'(t) and the six initial conditions for v(0).
Thus, we end having to integrate 12 first order differential equations. What happens is that, each time i set a non-null initial vector for v(0), for some initial conditions of the CR3BP (just like the one used to generate the above graphics) the results obtained are different than those obtained by the integration of only the six original equations since the system "collapses" going to x = y = 0 and the integrator tells me some errors rather than "integration successful", differently than what should happen. Here, v(0) = (1, 0, 0, 1, 0, 0).
The only circunstance in which i have the same results for both the integrations is when v(0) = (0, 0, 0, 0, 0, 0), but then I can't get values for the Fast Lyapunov Indicator.
Here's the code snippet of the new function to include the six new equations for the Lyapunov Indicator:
def func2(init, t, mu):
#x0, y0, z0, vx0, vy0, vz0 = init
x0, y0, z0, vx0, vy0, vz0 = init[0], init[1], init[2], init[3], init[4], init[5]
v0, v1, v2, v3, v4, v5 = init[6], init[7], init[8], init[9], init[10], init[11]
#print(init)
dx1dt1 = vx0
dy1dt1 = vy0
dz1dt1 = vz0
dx2dt2 = 2*dy1dt1+d_omega_x(mu, x0, y0, z0)
dy2dt2 = -2*dx1dt1+d_omega_y(mu, x0, y0, z0)
dz2dt2 = d_omega_z(mu, x0, y0, z0)
r1 = r11(mu, x0, y0, z0)
r2 = r22(mu, x0, y0, z0)
jacobiana = [v3,
v4,
v5,
(v0*(mu*(-3*mu - 3*x0 + 3)*(-mu - x0 + 1)/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) -
mu/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(3/2) +
(1 - mu)*(-3*mu - 3*x0)*(-mu - x0)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2) -
(1 - mu)/(y0**2 + z0**2 + (mu + x0)**2)**(3/2) + 1.0) +
v1*(-3*mu*y0*(-mu - x0 + 1)/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) -
3*y0*(1 - mu)*(-mu - x0)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2)) +
v2*(-3*mu*z0*(-mu - x0 + 1)/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) -
3*z0*(1 - mu)*(-mu - x0)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2)) + 2*v4),
(v0*(-mu*y0*(-3*mu - 3*x0 + 3)/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) -
y0*(1 - mu)*(-3*mu - 3*x0)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2)) +
v1*(3*mu*y0**2/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) -
mu/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(3/2) +
3*y0**2*(1 - mu)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2) -
(1 - mu)/(y0**2 + z0**2 + (mu + x0)**2)**(3/2) + 1.0) +
v2*(3*mu*y0*z0/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) +
3*y0*z0*(1 - mu)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2)) - 2*v3),
(v0*(-mu*z0*(-3*mu - 3*x0 + 3)/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) -
z0*(1 - mu)*(-3*mu - 3*x0)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2)) +
v1*(3*mu*y0*z0/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) +
3*y0*z0*(1 - mu)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2)) +
v2*(3*mu*z0**2/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) -
mu/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(3/2) +
3*z0**2*(1 - mu)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2) -
(1 - mu)/(y0**2 + z0**2 + (mu + x0)**2)**(3/2) + 1.0))]
fli = jacobiana
dv1 = fli[0]
dv2 = fli[1]
dv3 = fli[2]
dv4 = fli[3]
dv5 = fli[4]
dv6 = fli[5]
return [dx1dt1, dy1dt1, dz1dt1, dx2dt2, dy2dt2, dz2dt2, dv1, dv2, dv3, dv4, dv5, dv6]
What to do? It seems to clearly be a problem of floating point precision, since i get some different results every new time I run the code. When I increase the number of points in np.linspace (in this case to 2^20 points), i tend to get correct results, but I can't aways be handling more than a million points when, for the other situation, I can get correct results with 4 times less data. I need to apply a Continuous Wavelet Transform for the data and so it becomes really consuming.
Once more, i'm sorry if the question is rather long or confusing, if needed I can provide more informatio about it.
Anyway, thank you very much for the attention, stay safe.
EDIT:
As pointed by Lutz and following the natural dynamics of the system, by which Chaotic orbits are defined by exponentially increasing value of the Lyapunov Indicator -- which is actually defined by log10 of the norm, not only the norm --, it turns out the initial vector v(0) must be reasonably small so the result won't overflow. Trying (1e-10, 0, 0, 0, 0, 0) returns the correct integration. The curve profile for the Lyapunov indicator is also correct, just shifted by a factor log10(1e-10).
Thank you very much once again!
This is probably due to the step size control being also influenced by the rapidly growing v vector. Either by regulating step sizes rapidly down due to stiffness, or more likely, due to increasing the step size to match the dominant components, thus becoming unsuitable for an exact integration of the original trajectory. This rapid growth is the reason that Lyapunov exponents were introduced, as they capture this growth in nicely bounded numbers.
What you can do is to split up the integration into smaller chunks and normalize the v vector at the start of each chunk. One would have to experiment on how long it takes until the v component unduly dominates the step size control. As the coupling is purely multiplicative, the dynamic theoretically is linear. So it could also help if you scale the initial v to have norm 1e-100.
First however check what error tolerances you use. Setting them narrower also tends to stabilize the computation. You might also get some progress be setting the maximal step size hmax to half or so of the external step.
Or you could do the Lyapunov exponent computation like I explored in https://scicomp.stackexchange.com/questions/36013/numerical-computation-of-lyapunov-exponent. This approach however increases a system of dimension n by the n x n matrix of eigen/singular vectors and the n products of exponents times time.
