Trajectories of $L_4$ and Lyapunov Characteristic Exponents in the Generalized Photogravitational Chermnykh-Like Problem

0.32 MB·English

by  Badam Singh Kushvah

Astrophysics & Space Science, doi:10.1007/s10509-011-0632-y
Trajectories of L and Lyapunov Characteristic Exponents in the
4
Generalized Photogravitational Chermnykh-Like Problem
1
1
Badam Singh Kushvah
0
2
Department of Applied Mathematics, Indian School of Mines, Dhanbad - 826004
n
a Jharkhand(India)
J

8
[email protected],[email protected]
2

]
P
ABSTRACT
E
.
h
p
The dynamical behaviour of near by trajectories is being estimated by Lya-
-
o
punov Characteristic Exponents(LCEs) in the Generalized Photogravitational
r
t Chermnykh-Like problem. It is found that the trajectories of the Lagrangian
s
a
point L move along the epicycloid path, and spirally depart from the vicinity
[ 4

of the point. The LCEs remain positive for all the cases and depend on the
1
v initial deviation vector as well as renormalization time step. It is noticed that
4
3 the trajectories are chaotic in nature and the L4 is asymptotically stable. The
4
effects of radiation pressure, oblateness and mass of the belt are also examined
5
. in the present model.
1
0
1
Subjectheadings: Trajectory:LagrangianPoint:LCEs:Photograviational:Chermnykh-
1
: Like Problem:RTBP
v
i
X
r
a
1. Introduction
In present paper our aim is to obtain trajectories of L and is to estimate the rate of
4
deviation for initially closely related trajectories in the modified restricted three body prob-
lem model(as in Kushvah (2008, 2009a)) with radiation from Sun, oblateness of the second
primary(massive body) and influence of the belt. It is supposed that the primary bodies and
a belt are moving in a circular orbits about the common center of mass of both primaries.
First time such problem was discussed by Chermnykh (1987) and its importance in astron-
omy has been addressed by Jiang and Yeh (2004a). More generalized cases of the problem
were studied by many scientists such as Jiang and Yeh (2004b), Papadakis (2004),Papadakis
– 2 –
(2005) and Jiang and Yeh (2006); Yeh and Jiang (2006). The effect of radiation pressure,
Poynting-Robertson(P-R)drag and oblateness on the linear stability and nonlinear stability
of the L have been discussed by Kushvah and Ishwar (2006) ; Kushvah et al. (2007a,b,c).
4(5)
In our article Kushvah (2010), we have described the design of the trajectory and analysis
of the stability of collinear point L in the Sun-Earth system.
2
The first fundamental article about LCN’s was written by Oseledec (1968)in their study
of the ergodictheory ofdynamical system andBenettin et al. (1980)presented explicit meth-
ods for computing all LCEs of a dynamical system. Then Jefferys and Yi (1983) examined
stability in the restricted problem of three bodies with Liapunov Characteristic number.
First time Wolf et al. (1985) presented an algorithm with FORTRAN code that allows to
estimate non-negative Lyapunov Exponents(LEs) from an experimental time series. San-
dri (1996) have presented method for numerical calculation of Lyapunov Exponents for a
smooth dynamical system with Mathematica[Wolfram (2003)] code. Tancredi et al. (2001)
compared two different methods to compute Lyapunov Exponents(LEs). They have shown
that since the errors are introduced in the renormalization procedure, it is natural to expect
a dependency of the estimated LCEs with the number of renormalization performed in the
sense that the smaller the step the worse the estimation. In his study they made conclusion
that the two-particle method is not recommended to calculate LCEs in these cases where
the solution can fall in a region of regular or quasi regular solution of the phase space. For
a region of strong stochastically the LCEs calculated with the two-particle method gives
acceptable value.
This paper is organized as follows: In section 2, we state the model of the dynamical
system and compute the trajectories of L . Section 3 gives method to compute the LCEs,
4
where subsection 3.1 presents the first order LCEs for various set values of parameters, time
ranges and renormalization time steps. Section 4 presents comment about stability using
trajectories of L . Lastly, section 5 concludes the paper.
4
2. Trajectory of L
4
It is supposed that the motion of an infinitesimal mass particle be influenced by the
gravitational force from the two primaries(massive bodies) and a belt of mass M . We also
b
assume thatinfinitesimal massdoesnot influence themotionof thetwo massive bodieswhich
move in circular orbit under their mutual gravitational attraction. Let us assume that m
1
and m be the masses of the bigger and smaller primary respectively, m be the mass of the
2
infinitesimal body. The units are normalized by supposing that the sum of the masses to be
unity, the distance between both massive bodies to be unity. The rotating frame normalized
– 3 –
torotatewithunitangularvelocity andthetimeisnormalizedinsuch awaythatthetimefor
one period as a unit so that, the Gaussian constant of gravitational k2 = 1. For the present
model, perturbed mean motion n of the primaries is given by n2 = 1+3A2 + 2Mbrc , where
2 (r2+T2)3/2
c
T = a + b, a,b are flatness and core parameters respectively[as in Yeh and Jiang (2006)]
which determine the density profile of the belt; where r2 = (1 µ)q2/3+µ2, A = re2−rp2 is the
c − 1 2 5r2
oblateness coefficient of m ; r , r are the equatorial and polar radii of m respectively, r is
2 e p 2
the distance between primaries and the radius of the belt; µ = m2 is a mass parameter;
m1+m2
q = 1 Fp is a mass reduction factor and F is the solar radiation pressure force which is
1 − Fg p
exactly apposite to the gravitational attraction force F . In a rotating reference frame the
g
coordinates of m and m are ( µ,0) and (1 µ,0) respectively. We consider the model
1 2
− −
proposed by Miyamoto and Nagai (1975), and equations of motion are given as in Kushvah
(2008) and Kushvah (2009a):
x¨ 2ny˙ = Ω , (1)
x
−
y¨+2nx˙ = Ω , (2)
y
where
(1 µ)q (x+µ) µ(x+µ 1) 3µA (x+µ 1)
Ω = n2x − 1 − 2 −
x − r3 − r3 − 2 r5
1 2 2
M x W (x+µ)
b 1
(x+µ)x˙ +yy˙ +x˙ ny ,
−(r2 +T2)3/2 − r2 (cid:20) r2 { } − (cid:21)
1 1
(1 µ)q y µy 3µA y
Ω = n2y − 1 2
y − r3 − r3 − 2 r5
1 2 2
M y W y
b 1
(x+µ)x˙ +yy˙ +y˙ +n(x+µ) ,
−(r2 +T2)3/2 − r2 (cid:20)r2{ } (cid:21)
1 1
n2(x2 +y2) (1 µ)q µ µA M
1 2 b
Ω = + − + + +
2 r1 r2 2r23 (r2 +T2)1/2
(x+µ)x˙ +yy˙ y
+W narctan ,
1(cid:26) 2r2 − (cid:18)x+µ(cid:19)(cid:27)
1
(1 µ)(1 q )
W = − − 1 , r2 = (x+µ)2 +y2, r2 = (x+µ 1)2 +y2.
1 c 1 2 −
d
TheparameterW isconsideredduetoP-Rdrag[morereviewinPoynting(1903),Robert-
1
son (1937), Chernikov (1970), Murray (1994) and Kushvah (2009b)]. Where r , r are the
1 2
distances of m from first and second primary respectively. The dimensionless velocity of the
– 4 –
light is supposed to be c = 299792458. Then from equations (1) and (2) energy integral is
d
given as:
1
E = x˙2 +y˙2 Ω(x,y,x˙,y˙) = (Constant) (3)
2 −
(cid:0) (cid:1)
where the quantity E is an energy integral related to the Jacobi’s constant C(= 2E).
−
For numerical computation of equilibrium points, we divide the orbital plane Oxy into
three parts x µ, µ < x < 1 µ and 1 µ x with respect to the primaries.
≤ − − − − ≤
For the simplicity, we set µ = 9.537 10 4,T = 0.01. The equilibrium points are given by
−
×
substitutingΩ = Ω = 0,andpresentedinfigure1whenq = 0.75,A = 0.25,M = 0.25. In
x y 1 2 b
this figure the dark blue dotes present the position of L (5) : (x = 0.347988,y = 0.70645),
4
±
the light blue represent the collinear equilibrium points L : x = 0.753578,L : x = 1.14795
1 2
and L : x = 0.788385 for which y = 0.
3
−
The equations (1-2) with initial conditions x(0) = 1 2µ, y(0) = √3, x(0) = y (0) = 0
−2 2 ′ ′
are used to determine the trajectories of L for different possible cases. At at time t = 0,
4
the origin of coordinate axes is supposed at the equilibrium point.
In the present model all the computed trajectories of the L follow approximately the
4
same path described by an epitrochoid whose parametric equations are given as:
a +b
1 1
x(t) = (a +b )cost d cos t (4)
1 1 1
− (cid:18) b (cid:19)
1
a +b
1 1
y(t) = (a +b )sint d sin t (5)
1 1 1
− (cid:18) b (cid:19)
1
where a is radius of a fixed circle, b is radius of rolling circle and d is distance form
1 1 1
center of rolling circle to to the point(x(t),y(t)) which forms a trajectory. It is evident from
above equations that if d depends on time then orbit is unstable and trajectory moves
1
spirally outward the vicinity of the initial point.
Whenq = 0.75,A = 0, thetrajectoryisshowninfigure2withpanels(a-d)forM = 0.0
1 2 b
and panels(e-f) for M = 0.0, where frames(a& d) 0 t 50, (b&e) 50 t 100 and (c&f)
b
≤ ≤ ≤ ≤
0 t 200. It is clear from figure that if 0 t 50 and M = 0.0 the trajectory of L is
b 4
≤ ≤ ≤ ≤
similar to the curve described by epitrochoid (4, 5) for a = 1/7,b = 1 = d . When t > 50
1 1 1
then d becomes function of time t, and the trajectory moves spirally outward. When M =
1 b
0.25, the trajectory follows the path correspond to parameters a = 1,b = √5(irrational),
1 1
d = 3. Herethevalueofb isirrationalnumber which shows thatthemotionisnonperiodic.
1 1
When q = 0.75,M = 0.25, figure 3 depicts the trajectory for L with frames(a-c) for
1 b 4
A = 0.25 and frames(d-f) for A = 0.50. In frame(a) 0 t 50, (b) 50 t 75 and
2 2
≤ ≤ ≤ ≤
– 5 –
1
L4
0.5
L3 L1 L2
y 0 S
J
-0.5
L5
-1
-1.5 -1 -0.5 0 0.5 1 1.5
x
Fig. 1.— The position of equilibrium points.
y@tD y@tD
30
7.5 10
20
5
2.5 5 10
0 x@tD 0 x@tD
0
-2.5 -10
-5 -5 -20
-7.5
-7.5-5-2.50 2.5 5 7.5 -10 -5 0 5 -20-10 0 10 20 30
HaL HbL HcL
y@tD y@tD
7.5 15
20
5 10
2.5 5 10
x@tD
0 0 0 x@tD
-2.5
-5 -5 -10
-7.5 -10 -20
-10 -15
-10 -5 0 5 -15-10-5 0 5 10 15 -30-20-10 0 10 20
HdL HeL HfL
Fig. 2.— Trajectory of L when q = 0.75,A = 0.0 in frame(a) blue curve for 0 t 50,
4 1 2
≤ ≤
(b) red for 50 t 100 and (c) green for 0 t 200. Frames(a-c) M = 0.0, and frames
b
≤ ≤ ≤ ≤
(d-f) M = 0.25.
b
– 6 –
(c) 0 t 77 while (d) 0 t 8.2, (e) 8.2 t 8.3 and (f)0 t 9. It is clear from
≤ ≤ ≤ ≤ ≤ ≤ ≤ ≤
frames(a-c) that the trajectory moves along approximately epicycloid path, when t increases
it departs form the vicinity of L . The region of stability shrinks and trajectory moves along
4
a single cusped epicycloid, then it departs far from the initial point. Hence oblateness effect
is significant factor to reducing the stability region.
3. Lyapunov Characteristic Exponents(LCEs)
Itiswellknownthat, ifLCE> 0forsomeinitialconditionswhichindicatesthetrajectory
of initial condition is unstable. If LCE= 0 for some values of initial conditions the orbit is
neutrally stable and which corresponds to regular motion. If LCE< 0, the corresponding
orbit is asymptotically stable. Now suppose S be a 4-dimensional phase space such that
S = X : X = [x(t),y(t),p (t),p (t)]Tran , then the time evaluation of the orbit is governed
x y
{ }
by the equation
tran
∂H ∂H ∂H ∂H
X˙ = f(X) = = J DHX (6)
4
(cid:20) ∂x ∂y − ∂p − ∂p (cid:21)
x y
where D = ∂ and
∂X
0 0 1 0
J =  1 0 0 1. (7)
4
−
0 1 0 0
 − 
The dynamical system is described by the Hamiltonian H which depends on Jacobian
constant and given by
1
H = p2 +p2 +n(yp xp ) U(x,y) (8)
2 x y x − y −
(cid:0) (cid:1)
where p ,p are the momenta coordinates given by
x y
∂H ∂H
p˙ = , p˙ = ,
x y
−∂p −∂p
x y
n2(x2 +y2)
U(x,y) = Ω ,
− 2
Consider v = (δx,δy,δp ,δp ) be a deviation vector from initial condition X(0) such
x y
that v = 1. Then the variational equation is given
0
|| ||
v˙ = Df(X).v (9)
– 7 –
δx˙ 0 n 1 0 δx
δy˙   n 0 0 1δy
or = − , (10)
δp˙ Ut Ut 0 n δp
 x  xx xy  x
    
δp˙  Ut Ut n 0δp 
 y  yx yy −  y
where superscript t over partial derivatives of U indicates their respective values at t etc.
Then the Lyapunov Characteristic Exponent is given by
v(t)
λ(v(t)) = lim log || ||. (11)
t > v(0
− ∞ || ||
For numerical computation of LCEs we use method presented in Skokos and Gerlach
(2010) and Skokos (2010). To avoid overflow in numerical computation, we partition the
closed interval I = [t ,Tmax] into n sub intervals with time step ∆t and the time to run
0 1
from 0 to Tmax i.e.
P(I) = 0 = t ,t ,t ,t ,...t ,t ,...,t = Tmax , (12)
{ 0 1 2 3 k−1 k n1 }
then equation (11) can be written as
n1
λ(v(t)) = lim logα(t ), (13)
k
n1t−>∞Xk=0
where α(t ) = v(t ) . To determine first order LCEs in next section, we will use initial
k k
|| ||
vector X(0) = (0.499046,0.866025, 0.866025,0.499046) for classical RTBP(q = 1,A =
1 2
−
0,M = 0) and X(0) = (0.337957,0.81415, 0.954676,0.39629) for modified RTBP( q =
b 1
−
0.75,A = 0.25,M = 0.25). As in figure 4, at each step v(t ) will be evaluated from (6,10)
2 b k
using X(t ) and unit deviation vector vˆ(t ) = V(t )(say).
k 1 k 1 k 1
− − −
3.1. First Order LCEs
Now consider R4 = LD , R3 = LD , R2 = LD and R1 = LD spaces such that
1 2 3 4
LD LD LD LD . To find the first order LCE(λ ),(i = 1,2,3,4), we choose
1 2 3 4 i
⊃ ⊃ ⊃
initial unitdeviation vectors fromLD LD : v = (1/2,1/2,1/2,1/2),v = (0, 1 , 1 , 1 ),
1\ 2 11 12 √3 √3 √3
v = (0,0, 1 , 1 ), v = (0,0,0,1). The values of LCEs are presented in log-log plot figure
13 √2 √2 14
5 for t = 0 10000 when q = 0.75,A = 0.25,M = 0.25 left panel corresponding to
1 2 b
−
∆t = 1 and right for ∆t = 2. Initially the values of LCE(λ ) are different, they are shown by
1
curves(I)-(IV) correspond to four vectors respectively, but when t increases they merge into
a single curve. To obtain LCE(λ ), we choose initial unit deviation vectors from LD LD
2 2 3
\
such thatv = ( 1 , 1 , 1 ,0), v = (0, 1 , 1 ,0), v = (0,0,1,0). Figure6 shows LCE(λ )
21 √3 √3 √3 22 √2 √2 23 2
– 8 –
y@tD y@tD
3
1
2 2
1 x@tD 0 x@tD
0 -1 0
-1 -2 -2
-2
-3 -3 -4
-2-1 0 1 2 3 -3-2-1 0 1 2 3 -3-2-10 1 2 3 4
HaL HbL HcL
y@tD 0.95 0.9540.9580.962 y@tD
2 -0.08 2
x@tD
1 x@tD -0.09 0
-0.1
0 -2
-0.11
-1 -0.12 -4
-2 -0.13 -6
-2 -1 0 1 0.95 0.9540.9580.962 -6 -4 -2 0 2
HdL HeL HfL
Fig. 3.— Trajectory of L when q = 0.75,M = 0.25, where frames(a-c) for A = 0.25 and
4 1 b 2
frames(d-f) for A = 0.50.
2
Fig. 4.— In Plot (I) only one step is used to obtain LCEs while in (II) more that one
normalization steps are used and V(t ) denotes the unit deviation vector for all k.
k 1
−
– 9 –
6
LCE4
L
Hlog
2
I(cid:144)II
IV
0 III
0 2 4 6 8
logHTL
Fig. 5.— LCE(λ ) when q = 0.75,A = 0.25,M = 0.25 and 0 t 10000; curve (I)
1 1 2 b
≤ ≤
v = (1/2,1/2,1/2,1/2),(II):v = (0, 1 , 1 , 1 ), (III):v = (0,0, 1 , 1 ) and (IV):v =
11 12 √3 √3 √3 13 √2 √2 14
(0,0,0,1).
8
6 6
LLCE4 LLCE4
Hlog Hlog
2
I(cid:144)II 2
0 III
I(cid:144)II
0 III
0 2 4 6 8 0 2 4 6 8
logHTL logHTL
Fig. 6.— LCE(λ ) when q = 0.75,A = 0.25,M = 0.25 and 0 t 10000; curves (I)
2 1 2 b
≤ ≤
v = ( 1 , 1 , 1 ,0),(II):v = (0, 1 , 1 ,0) and (III):v = (0,0,1,0).
21 √3 √3 √3 22 √2 √2 23
8
6 6
HLlogLCE 24 HLlogLCE4
I 2
II
0 III I(cid:144)II
0 III
-2
0 2 4 6 8 0 2 4 6 8
logHTL logHTL
Fig. 7.— LCEs when 0 t 10000, q = 0.75,A = 0.25 and M = 0.25; where LCE(λ ):(I)
1 2 b 3
≤ ≤
v = ( 1 , 1 ,0,0),(II):v = (0,1,0,0) and LCE(λ )(III):v = (1,0,0,0).
31 √2 √2 32 4 41
– 10 –
when q = 0.75,A = 0.25,M = 0.25, where left panel corresponds to ∆t = 1 and right for
1 2 b
∆t = 2.
NowforcomputationofLCE(λ ),wechooseinitialunitdeviationvectorsfromLD LD .
3 3 4
\
The results are presented in figure 7 for q = 0.75,A = 0.25,M = 0.25 and 0 t 10000
1 2 b
≤ ≤
with left frame for ∆t = 1 and right for ∆t = 2. In figures 8 and 9, curves are plotted
when q = 1,A = 0,M = 0, where left panel corresponding to ∆t = 0.1 and right for
1 2 b
∆t = 1. In figure 8, curves are labeled as(I) v = ( 1 , 1 , 1 ,0),(II):v = (0, 1 , 1 ,0),
21 √3 √3 √3 22 √2 √2
(III):v = (0,0,1,0) and in figure 9, curves are plotted for 0 t 100, where (I)
23
≤ ≤
v = ( 1 , 1 ,0,0),(II):v = (0,1,0,0). The curves are in wave form with decreasing ampli-
31 √2 √2 32
tudes which tend to zero at infinity and curves become constant.
To determine LCE(λ ), we choose v = (1,0,0,0) from LD . The corresponding LCE
4 41 4
is shown by curve (III) in figure 7:q = 0.75,A = 0.25,M = 0.25 with left frame for ∆t = 1
1 2 b
and right for ∆t = 2. In figure 10, we consider q = 1,A = 0.0,M = 0.0 in which curve(I)
1 2 b
represents renormalization time step ∆t = 0.1 and (II) for ∆t = 1. It can be seen that (I)
is a smooth curve and (II) is initially stepped curve but both curves are initially increasing
in nature and after certain time they become constant. The LCEs(λ ),(i = 1,2,3,4) are
i
presented in Table 1 for initial point X(0) = (0.337957,0.81415, 0.954676,0.39629) and
−
q = 0.75,A = 0.25,M = 0.25. It is clear from figures and Table that all first order LCEs
1 2 b
are positive for various set values of parameters and renormalization time steps. This shows
that the present dynamical system is stochastic. It is also noticed that if Tmax is not very
large the LCEs depend on the choice of the renormalization time step as well as initial
deviation vectors while if Tmax is very large then LCEs depend on renormalization time
step only.
1
0.75 III
LE 0.5 I
C
L
Hg 0.25
o
l
0
-0.25
II
-0.5
1.5 2 2.5 3 3.5 4 4.5
logHTL
Fig. 8.— LCE(λ ) when q = 1,A = 0,M = 0 and 0 t 100; where (I):v =
2 1 2 b 21
≤ ≤
( 1 , 1 , 1 ,0),(II):v = (0, 1 , 1 ,0) and (III):v = (0,0,1,0).
√3 √3 √3 22 √2 √2 23