Computation of coefficients of the normal form in a singular case

SADOV, Sergey Yu. Institute of Applied Mathematics, Moscow, Russia

Keywords: Systems of ordinary differential equations, normal form, singular point, satellite motion.

1. Consider a system of ordinary differential equations

\frac{dX}{dt} = εF (X, t), X = (x_{1}, …, x_{n}),  (1)

where ε is a small parameter, and the vector-function F = (f_{1}, …, f_{n}) is analytic in X in some region O and 2π-periodic in t. Its normal form is an autonomous system

Ẏ = ∑_{k=1}^{∞} ε^{k}G_{k}(Y) ≡ G(Y, ε).  (2)

Variables X in (1) and Y in (2) are connected by the normalizing transformation

X = Y + ∑_{k=1}^{∞} ε^{k}H_{k}(Y, t).  (3)

The vector-functions H_{k}(Y, t) and G_{k}(Y) are analytic in Y ∈ O, but the series in (2) and (3) diverge, as a rule. The functions H_{k} are 2π-periodic in t.

A computation of the coefficients H_{k} and G_{k} represents the well-known recurrent procedure of averaging. For example,

G_{1} = M[F],  H_{1} = L[F],
G_{2} = M [\frac{∂F}{∂Y} · H_{1} − \frac{∂H_{1}}{∂Y} · G_{1}] , H_{2} = L[ the same ],

where the operators M and L are given by

M [f] = f_{0} = \frac{1}{2π} ∫_{−π}^{π} f(t) dt  (4)

(L [f] )(t) = ∑_{j≠0} \frac{f_{j}}{ij} exp(ijt)  (5)

for a 2π-periodic function

f(t) = ∑_{j=−∞}^{∞} f_{j} exp(ijt).

Hence, a numerical computation of the normal form coefficients amounts to recursive quadratures.

2. We consider a system of the form (1) depending on an additional parameter e. We have n = 2 and right-hand sides in (1)

f_{1}(X, t) = x_{2}, f_{2}(X, t) = − (\frac{1 + e cos ν(t, e)}{1 − e^{2}})^{3} sin(mt − 2ν(t, e) + x_{1}),

(6)

where ν(·, e) is the inverse function for

t(ν, e) = (1 − e^{2})^{3/2} ∫_{0}^{ν} (1 + e cos ν)^{−2} dν.

The system under consideration arises from the equation

(1 + e cos ν)\frac{d^{2}δ}{dν^{2}} − 2e sin ν \frac{dδ}{dν} + μ sinδ = 4e sin ν  (7)

that describes a satellite motion about its center of mass in the plane of its elliptic orbit [1]. Here e is the orbit’s eccentricity, μ the inertial parameter which is small for the satellite being almost dynamically symmetric. The independent variable ν is the so called true anomaly, and the dependent variable δ represents the doubled angle of the satellite’s rotation with respect to the radius-vector.

The system (1), (6) follows ¿from (7) by the change of variables [3]

μ = ε^{2},  ν = ν(t, e),  δ = mt − 2ν + x_{1},  x_{2} = ε^{−1} dx_{1}/dt,

m being an auxiliary integer parameter.

A numerical study of the equation (7) for e close to 1 is very difficult (see, e.g. [2]), because the equation is singular at the point e = 1, ν = π. Also the straight computation of the normal form coefficients meets with a precision loss because of singularities of the integrands. We propose a regularization of the computational procedure. For the coefficients, which are bounded as e → 1, it enables us to compute them numerically without the precision loss, and for the unbounded coefficients it gives their asymptotics.

Our method is based on integrating by a path in the complex plane and factorizing integrands into functions with close but diverse poles. Apart from the poles we have a regular expression, which one can calculate numerically.

It is convenient to use a parameter

η = \sqrt{1 − e^{2}}

instead of (or together with) e as e close to one.

3. Before further discussions, we introduce some notations. (a). We define functions γ_{±}(t; e, m) by

γ_{±}(t; e, m) = η^{−6}(1 + e cos ν(t, e))^{3}) exp {±i(mt − 2ν(t, e))} . (8)

(b) Let {a_{±}, …, z_{±}} be two sets of nonnegative integer indices, one of them possible empty. Let γ̃(t, e, m) be the product of functions L^{a_{+}} γ_{+} … L^{z_{+}} γ_{+} and L^{a_{−}} γ_{−} … L^{z_{−}} γ_{−} (with the operator L defined by (5)). We set

Ξ_{a_{−},...,z_{−}}^{a_{+},...,z_{+}}(e, m) = M(γ̃).  (9)

For example, M γ_{+} = Ξ^{0}, M γ_{−} = Ξ_{0}.

Proposition The expressions (9) with the even sum of indices are invariant under the permutation of the lower and upper sets of indices.

Formulas for some coefficients G_{j} of the normal form (2) in our case are written out in [4, Ch. V]. It appears that

G(y_{1}, y_{2}, ε, e, m) = (εy_{2}, ε sin y_{1} ∑_{j} ε^{j}p_{j}(y_{1}; e, m) + y_{2}q(y_{1}, y_{2}, ε, e, m))^{†},

and p_{2l+1} ≡ 0. Formulas from [4] for the first three nonzero p_{k} ’s take the following form in our notations.

Theorem 1 We have

p_{0} = −Ξ_{0} = −Ξ^{0}, p_{2} = −Ξ_{1,1} · cos y_{1} = −Ξ^{1,1} · cos y_{1},

p_{4} = −\frac{1}{4}Ξ_{1,1}^{2} + (\frac{3}{8} − \frac{3}{2} cos^{2} y_{1}) Ξ_{2,2}^{0} + Ξ_{0} · (−\frac{5}{4}Ξ_{2}^{2} + (\frac{11}{8} − 3 cos^{2} y_{1}) Ξ_{2,2}) .

We introduce a complex variable

ζ = \frac{exp(iν) + ζ_{−}}{1 + ζ_{−} exp(iν)} = \frac{1 + ζ_{+} exp(iν)}{exp(iν) + ζ_{+}},

where

ζ_{−} = \frac{1 − η}{e} , ζ_{+} = \frac{1 + η}{e} = 1/ζ_{−}.

Functions (8) have the explicit expressions in terms of ζ:

γ_{±}(ζ, e, m) = (\frac{2ζ}{1 ∓ η} )^{3} (ζ − ζ_{±})^{−5}(1 − ζ_{±}ζ)^{−1}ζ^{±m} exp \{±\frac{em}{2}(ζ^{−1} − ζ)\} .

Now we replace the integration over [−π, π] in (9) by the integration over the unit circle |ζ| = 1 with the weight function

ρ(ζ, e) = \frac{dt}{dζ} = \frac{e(ζ − ζ_{−})(ζ − ζ_{+})}{2iζ^{2}}.

Analytic functions ργ_{+} and L^{k} γ_{+} (k ≥ 1) have no singularities outside the unit circle, and the functions ργ_{−} and L^{k} γ_{−} (k ≥ 1) are regular inside the unit circle except of the point ζ = 0. This fact allows us to compute the expressions Xi_{...} (or Xi^{...}) as integrals over circles not close to the singular points of integrands ζ = ζ_{+} (or ζ = ζ_{−}). A Turbo-Pascal program (see Appendix) is designed for computing p_{0}(m, e) and p_{1}(m, e). This program gives results that are compatible with analytical asymptotics (e.g. (11) below); on the other hand, the old numerical results by Lutze and Abbit [5] proved to be incorrect as m ≤ 3 and e ≥ 0.8.

For Ξ ’s with both upper and lower indices the above regularization failes because the original path of integration lies between the singular points ζ_{+} and ζ_{−} of an integrand, and one cannot move it aside these points without crusial change of the value of the integral. However, in this case we can explicitly write out singular terms of the integrands. The formula (10) below is derived in this way.

So, we formulate the main analytical results.

Theorem 2 (a) For a fixed integer m and e → 1 the functions Ξ_{0}(e, m), Ξ_{1,1}(e, m), Ξ_{2,2}^{0}(e, m), Ξ_{2}^{2}(e, m) and Ξ_{2,2}^{0}(e, m) have finite limits, and the asymptotics of the function Ξ_{1,1}^{2} as e → 1 is the following:

Ξ_{1,1}^{2}(e, m) = \frac{1}{162}η^{−3} − \frac{7}{54}η^{−1} − \frac{32}{9}m ln^{2} η + O(ln η).  (10)

(b) The limit value of Ξ_{0}(e, m) is

Ξ_{0}(1, m) = \frac{m}{3} (−1 + J_{0}(m) + ∫_{0}^{m} J_{0}(x) dx − 2 ∑_{k=0}^{m−1} J_{k}(m) − J_{m−1}(m)) ,

(11)

J_{k} being the Bessel functions. In particular, Ξ_{0}(1, 2) = −1.018863 ….

It is interesting that for m = 0 one can obtain expressions for the quantities Ξ in terms of elementary functions due to rationality of the functions γ_{±}(ζ).

Theorem 3 For m = 0 we have

Ξ_{0}(e, 0) = Ξ_{1,1}(e, 0) = Ξ_{2,2}^{0}(e, 0) ≡ 0,

Ξ_{1,1}^{2}(e, 0) = −\frac{e^{−6}}{27} (\frac{(1 − η)(3η^{3} + 23η^{2} − η − 1)}{8η^{3}} + 8 log (\frac{2η}{1 + η} ) ) .

These results are proven in [6, 7].

4. The motivation of this work and our goal consists in describing periodic solutions of (7) for the parameters μ ≪ 1 and e → 1. Periodic solutions of (7) are the periodic solutions of (1), (refrhs) with m = 2. By Bruno’s theorem about analytical integral sets [4, Ch. III,§ 3], the periodic solutions lie at the analytical set A, given by the following equations in the normal form coordinates Y:

A = {Y, ε, e | G(Y, ε, e) = 0}.

There is always a branch of the set A of the form sin y_{1} = 0. Other branches, if they exist, correspond to so called critical periodic solutions. Such solutions in the region μ ≪ 1, e → 1 were recently discovered numerically by Bruno and Petrovich [8]. Our analysis doesn’t give a rigorous description of this branch because our results concern only a few beginning terms of the normal form, but if we suppose that the following terms are minor, we come to the following

Conjecture A branch of the set A, which corresponds to critical periodic solutions, exists for parameters e → 1, μ → 0, in the neighbourhood of the curve

μ = ±36 · 2^{1/4} · (−Ξ_{0}(1, 2))^{1/2} · (1 − e)^{3/4} + O((1 − e)^{5/4}.  (12)

For 0.9999 ≤ e ≤ 0.99999 the curve (12) differs from the numerically computed curve [8] by less than 5%.

References

[1] V.V. Beletskii. Motion of an Artificial Satellite about Its Center of Mass, NASA TT F-429, pp. 257–261.

[2] V.A. Sarychev. Topics in satellite orientation. Itogi nauki i tehniki, ser. ”Issledovanije kosmicheskogo prostranstva”, Vol. 11. Moscow: VINITI, 1978. (Russian)

[3] F.L. Chernous’ko. Resonance Phenomena in the Motion of a Satellite Relative to Its Mass Center. USSR Computational Mathematics and Mathematical Phys., 3 (1964), 699–713.

[4] A.D. Bruno. Local Methods in Nonlinear Differential Equations. Springer-Verlag: Berlin-Heidelberg-New York-London-Paris-Tokyo, 1989.

[5] F.H. Lutze, Jr., M.W. Abbit, Jr. Rotational locks for near-symmetric sattelites. Selest. Mech. 1:1 (1969), 31–35.

[6] S.Yu Sadov. Analysis of a function that defines rotational stability of almost symmetric satellite. Preprint of the Keldysh Institute of Applied Mathematics of Rus. Ac. Sci., Moscow, 1994, N 84 (Russian).

[7] S.Yu. Sadov. Coefficients of an averaged equation of satellite oscillations. Preprint of the Keldysh Institute of Applied Mathematics of Rus. Ac. Sci., Moscow, 1995, N 27 (Russian).

[8] A.D. Bruno and V.Yu. Petrovich. Regularization of satellite’s oscillations on a very stretched orbit. Preprint of the Keldysh Institute of Applied Mathematics of Rus. Ac. Sci., Moscow, 1994, N 4. (Russian)

APPENDIX. A program for computing Ξ_{0}(m, e) and Ξ_{1,1}(m, e).

(In this program there are notations ϕ(m, e) and −s(m, e) instead of Ξ_{0}(m, e) and Ξ_{1,1} respectively. References relate to the prepint [7].)

Program Coefficients;
{ @  S.Yu. Sadov, 1995.              }
{    In: Coefficients of an averaged equation of satellite
     oscillations. Preprint Inst.Appl.Math., Rus.Ac.Sci.,
     1995, N 27                      }
{    e-mail: sadovs@applmat.msk.su   }
Uses Complex_Ariphmetics ;
 (*  In the Unit Complex_Ariphmetics there are defined:
    type cmpl=record Re,Im:real end;
    procedure cAdd(a,b:cmpl; var c:cmpl); {c:=a+b}
    procedure cSub(a,b:cmpl; var c:cmpl); {c:=a-b}
    procedure cMul(a,b:cmpl; var c:cmpl); {c:=a*b}
    procedure cSqr(var c:cmpl); {c:=c^2}
    procedure cInv(var c:cmpl); {c:=1/c}
    procedure cExp(theta:real; var c:cmpl); {c:=exp(i*theta) }
  *)
Const         {manually defined constants}
      m=2   ; {parameter m }
      e_start=0; e_end=1;
      delta_e=0.1; {step of a table}

N=512; {Number of integration points}
Z=0.5;  {radius of integration contour}
Const         {automatically defined constants}
pi2=6.283185307179586477;
cPi2:cmpl=(Re:pi2; Im:0.0);
circ_step=pi2/N;
Zinv=1/Z;
eps=delta_e/2;
Type  fun=array[0..N-1] of cmpl;
Var  e, eta: real;
rho, gamma_minus: fun;
Procedure Calculation_Rho;
{here theta is the eccentric anomaly, rho=1/N * dt/d(theta),
calculated on the contour |zeta|=Z }
var     i: integer;   theta: real;
begin
  for i:=0 to (N-1) do begin
     theta:=i*circ_step;
     rho[i].Re:=(1-cos(theta)*(Z+Zinv)*(e/2))/N;
     rho[i].Im:=(sin(theta)*(Zinv-Z)*e/2)/N
  end
end;
Procedure Calculation_Gamma_minus; {see (32)}
var   theta : real;
 Procedure exp_imt (var expimt:cmpl);
     {expimt:=exp(-imt(zeta)) }
   var  c1,c2 : cmpl;
       arg, zm : real;
             i : integer;
   begin  {exp_imt}
      c1.Re:=(Z-Zinv)*cos(theta)* m *e/2;
      c1.Im:=(Z+Zinv)*sin(theta)* m *e/2;

 arg:=-theta*m + c1.Im;
 cExp(arg, c2);
 zm:=1;  {we want now to compute zm=1/Z^m}
 if m>0 then for i:=1 to m do zm:=zm*Zinv
     else for i:=1 to (-m) do zm:=zm*Z;
 expimt.Re:=c2.Re*exp(c1.Re)*zm;
 expimt.Im:=c2.Im*exp(c1.Re)*zm;
end;  {exp_imt}
var c1,c2, expimt : cmpl;
     i      : integer;
begin {Calculation_Gamma_minus}
for i:=0 to (N-1) do begin
theta:=circ_step*i;
c1.Re:=((1-eta)*Z+(1+eta)*Zinv)*cos(theta)/2-e;
c1.Im:=((1-eta)*Z-(1+eta)*Zinv)*sin(theta)/2;
cSqr(c1);
c2.Re:=-(Z+Zinv)*cos(theta)*e/2+1;
c2.Im:=(Zinv-Z)*sin(theta)*e/2;
cMul(c1,c2,c1);
cInv(c1);
exp_imt(expimt);
cMul(expimt,c1,gamma_minus[i]);
end
end; {Calculation_Gamma_minus}
Procedure op_M(var f:fun; var  Mf:cmpl); {see(39)}
var i: integer;
c: cmpl;
begin
Mf.Re:=0; Mf.Im:=0;
for i:=0 to (N-1) do begin
{computation of an integral by
  the method of rectangles}
cMul(f[i],rho[i],c);  cAdd(Mf,c,Mf);
end;
{Mf.Re:=Mf.Re/pi2;  Mf.Im:=Mf.Im/pi2}

end;  { M }
Procedure L(var f:fun; var Lf:fun); {see(41),(42)}
var  c1,c2,Mf: cmpl;   i: integer;
begin
op_M(f,Mf);
Lf[0].Re:=0; Lf[0].Im:=0;
cSub(f[0],Mf,c1);
cMul(c1,rho[0],c1);
for i:=1 to (N-1) do begin
  {integration by the method of triangles}
  c2:=c1;
  cSub(f[i],Mf,c1);
  cMul(c1,rho[i],c1);
  Lf[i].Re:=Lf[i-1].Re+(c1.Re+c2.Re)/2;
  Lf[i].Im:=Lf[i-1].Im+(c1.Im+c2.Im)/2;
end; {Now array Lf corresponds to funct.(41)}
op_M(Lf,Mf);
for i:=0 to (N-1) do begin
   cSub(Lf[i],Mf,Lf[i]);
   cMul(cPi2,Lf[i],Lf[i])
end
end;  { L }
Function phi: real;  {see (14)}
var c: cmpl;
begin
 op_M(gamma_minus, c);
 phi:=c.Re
end;  {phi}
Function s: real;    {see (15)}
var c: cmpl; Lf: fun; i:integer;
begin
 L(gamma_minus,Lf);
 for i:=0 to (N-1) do cSqr(Lf[i]);
 op_M(Lf,c);
 s:=-c.Re

end; {s}
BEGIN {main}
e:=e_start;
Repeat
eta:=sqrt(1.0-sqr(e)); {see(25)}
Calculation_rho;
Calculation_Gamma_minus;
writeln('phi(',m,' ,',e:5:2,')=',phi);
(* or:    writeln('s(',m,' ,',e:5:2,')=',s); *)
e:=e+delta_e;
Until (e>e_end+eps)
END.