Preliminary Orbit Determination of Artificial Satellites: A Vectorial Sixth-Order Approach

Abstract

A modified classical method for preliminary orbit determination is presented. In our proposal, the spread of the observations is considerably wider than in the original method, as well as the order of convergence of the iterative scheme involved. The numerical approach is made by using matricial weight functions, which will lead us to a class of iterative methods with a sixth local order of convergence. This is a process widely used in the design of iterative methods for solving nonlinear scalar equations, but rarely employed in vectorial cases. The numerical tests confirm the theoretical results, and the analysis of the dynamics of the problem shows the stability of the proposed schemes.

1. Introduction

The analysis of linear systems has a well-developed mathematical and computational theory. However, many of the applied problems in science and engineering are nonlinear. This situation is more complicated than the linear one, and the estimation of its solution needs a numerical treatment.

While computational engineering has achieved significant maturity, computational costs can be extremely large when high accuracy simulations are required. The development of a practical high-order solution method could diminish this problem by significantly decreasing the computational time required to achieve an acceptable error level (see, for instance, [1]).

The existence of an extensive literature on higher order methods reveals that they are only limited by the nature of the problem to be solved: in particular, the numerical solution of nonlinear equations and systems is needed in the study of dynamical models of chemical reactors [2] or in radioactive transfer [3]. Moreover, many of numerical applications use high precision in their computations; in [4], high-precision calculations are used to solve interpolation problems in astronomy; in [5], the authors describe the use of arbitrary precision computations to improve the results obtained in climate simulations; the results of these numerical experiments show that the high order methods associated with a multiprecision arithmetic floating point are very useful, because it yields a clear reduction in iterations. A motivation for an arbitrary precision in interval methods can be found in [6], in particular for the calculation of zeros of nonlinear functions.

In last decades, many researchers have proposed different iterative methods to improve Newton’s one, which is still the most used scheme in practice. These variants of Newton’s method have been designed by means of different techniques, providing in the most of cases multistep schemes. Some of them are extensions of one-dimensional schemes (see, e.g., [7, 8]), and others come from Adomian decomposition (see e.g., [9, 10]), specifically the methods proposed by Darvishi and Barati in [11, 12] with super cubic convergence and the schemes proposed by Cordero et al. in [13] with order of convergence 4 and 5. Another procedure to develop iterative methods for nonlinear systems is the replacement of the second derivative in Chebyshev-type methods by some approximation. In [14], Traub presented a family of multipoint methods based on approximating the second derivative that appears in the iterative formula of Chebyshev’s scheme, and Babajee et al. in [15] designed two Chebyshev-like methods free from second derivatives.

A common way to generate new schemes is the direct composition of known methods with a later treatment to reduce the number of functional evaluations (see e.g., [16–19]). A variant of this technique is the so called pseudocomposition, introduced in [20, 21]. Let us note that if the initial approximation or any of the successive estimations make the jacobian matrix almost singular, the convergence is not guaranteed. In some of these cases, the problem can be avoided by using some kind of pseudoinverse to solve the linear system involved in each step of the iterative process (see, for instance, [22, 23]).

Recently, for 𝑛=1, the weight-function procedure has been used to increase the order of convergence of known methods [7]. This technique can be also used, with some restrictions, in the development of high order iterative methods for systems: see, for example the papers of Sharma et al. [24, 25] and Abad et al. [26], where the authors apply the designed method to the software improvement of the Global Positioning System.

1.1. Preliminary Orbit Determination

A classical reference in preliminary orbit determination is F. Gauss (1777–1855), who deduced the orbit of the minor planet Ceres, discovered in 1801 and afterwards lost. The calculation of its trajectory by means of the procedure designed by Gauss marked the international recognition of Gauss and his work.

The first step in orbit determination methods is to obtain preliminary orbits, as the motion analyzed is under the premises of the two bodies problem. It is possible to set a two-dimensional coordinate system (see Figure 1), where the X axis points to the perigee of the orbit, the closest point of the elliptical orbit to the focus and center of the system, the Earth. In this picture, the true anomaly 𝜈 and the eccentric anomaly 𝐸 can be observed. In order to place this orbit in the celestial sphere and determine completely the position of a body in the orbit, some elements (called orbital or keplerian elements) must be determined. These orbital elements are as follows.(i)Ω (right ascension of the ascending node): defined as the equatorial angle between the Vernal point 𝛾 and the ascending node 𝑁; it orients the orbit in the equatorial plane.(ii)𝜔 (argument of the perigee): defined as the angle of the orbital plane, centered at the focus, between the ascending node 𝑁 and the perigee of the orbit; it orients the orbit in its plane.(iii)𝑖 (inclination): dihedral angle between the equatorial and the orbital planes.(iv)𝑎 (semimajor axis): which sets the size of the orbit.(v)𝑒 (eccentricity): which gives the shape of the ellipse.(vi)𝑇₀ (perigee epoch): time for the passing of the object over the perigee, to determine a reference origin in time. It can be denoted by a exact date, in Julian days, or by the amount of time ago the object was over the perigee.

Figure 1 

Size, shape, and anomalies in orbital plane 2-dimensional coordinate system.

The so-called Gauss’ method is based on the rate 𝑦 between the triangle and the ellipse sector defined by two position vectors, 𝑟₁ and 𝑟₂, from astronomical observations. This proportion is related with the geometry of the orbit and the observed position by 𝑦=1+𝑋(𝑙+𝑥), (1) where 𝑟𝑙=¹+𝑟₂4√𝑟₁𝑟₂cos𝜈₂−𝜈₁−1/22,𝑥=sin²𝐸₂−𝐸₁4,𝐸𝑋=²−𝐸₁−sin𝐸₂−𝐸₁sin³𝐸₂−𝐸₁./2 (2)

The angles 𝐸_𝑖, 𝜈_𝑖, 𝑖=1,2, are the eccentric and true anomalies, respectively, associated to the observed positions 𝑟₁ and 𝑟₂ (let us denote by 𝑟_𝑖 the modulus of vector 𝑟₁, 𝑖=1,2).

Equation (1) is, actually, the composition of the first Gauss equation 𝑦²=𝑚𝑙+𝑥 (3) and the second Gauss equation 𝑦²𝑦−1=𝑚𝑋, (4) where 𝑚=𝜇𝜏²/[2√𝑟₁𝑟₂cos((𝜈₂−𝜈₁)/2)]³, 𝜇 is the gravitational parameter of the motion, and 𝜏 is a modified time variable.

The original scheme used by Gauss (see [27]) was based on applying fixed point method on unified (1). By using the initial estimation 𝑦₀=1, the classical Gauss procedure gets to calculate the orbital elements only if the range of the true anomalies corresponding to the observed positions is lower than 𝜋/4. Our aim is to widen the admissible range of true anomalies up to 𝜋 by solving the same problem but using the nonlinear system formed by (3) and (4), being the unknowns 𝑦 and 𝐸₂−𝐸₁.

In the second section of this paper, we present a class of sixth-order Newton-type methods by using the composition technique and matricial weight functions. The convergence results of this family have been obtained by means of the 𝑛-dimensional Taylor expansion of the involved functions, by using the notation introduced in [28]. Section 3 is devoted to analyze the efficiency of the proposed methods applied on preliminary orbit determination and other nonlinear problems, compared with the classical Newton’, Traub’ [14], and Jarratt’s [29] methods. In Section 4, the preliminary orbit determination problem is revisited in order to analyze the stability of the mentioned schemes by means of dynamical planes, comparing the set of starting points that leads each of the methods to converge to the solution. We finish the paper with some conclusions and the references used in it.

2. The New Family and Its Convergence

In this section, we present our three-step iterative methods designed from Newton's one composing with itself twice, once with a “frozen” function and other one with a “frozen” Jacobian matrix. By using matricial weight functions in the second and third step, we proof that the methods of the proposed family have order of convergence six, under certain conditions of function 𝐹 and of weight functions.

Let us consider the following three-step iterative method which makes use of the weight functions: 𝑦^(𝑘)=𝑥^(𝑘)−𝐹^′𝑥^(𝑘)−1𝐹𝑥^(𝑘),𝑧^(𝑘)=𝑦^(𝑘)−𝐻𝜇^(𝑘)𝐹^′𝑦^(𝑘)−1𝐹𝑥^(𝑘),𝑥^(𝑘+1)=𝑧^(𝑘)−𝐺𝜇^(𝑘)𝐹^′𝑦^(𝑘)−1𝐹𝑧^(𝑘), (5) where 𝜇^(𝑘)=[𝐹^′(𝑦^(𝑘))]⁻¹𝐹^′(𝑥^(𝑘)) and 𝐻,𝐺∶ℝ^𝑛×𝑛→ℝ^𝑛×𝑛 are matricial functions.

Theorem 1. Let 𝛼∈𝐷 be a zero of a sufficiently differentiable function 𝐹∶𝐷⊆ℝ^𝑛→ℝ^𝑛 in a convex set 𝐷 with non-singular Jacobian in 𝛼. Let 𝐻 and 𝐺 be any function with the following conditions: 𝐻(𝐼)=0, 𝐻^′(𝐼)=(1/2)𝐼, 𝐻^′′(𝐼)=0 and 𝐺(𝐼)=𝐼, 𝐺^′(𝐼)=0, 𝐺^′′(𝐼)=(1/2)𝐼, being 𝐼 the identity matrix. Then, the scheme defined in (5) provides sixth order of convergence, whose error equation is given by 𝑒^(𝑘+1)=−32𝐶₃𝐶₂𝐶₃+14𝐶₂𝐶²₃+6𝐶₃𝐶³₂−𝐶₂𝐶₃𝐶²₂+𝐶³₂𝐶₃−4𝐶⁵₂𝑒^(𝑘)6+𝑂𝑒^(𝑘)7, (6) where 𝐶_𝑘=(1/𝑘!)[𝐹^′(𝛼)]⁻¹𝐹^(𝑘)(𝛼), 𝑘=2,3,4,…, and 𝑒^(𝑘)=𝑥^(𝑘)−𝛼.

Proof. If we use Taylor expansion of 𝐹(𝑥^(𝑘)) and 𝐹^′(𝑥^(𝑘)) around 𝛼, we obtain 𝐹𝑥^(𝑘)=𝐹^′(𝛼)𝑒^(𝑘)+𝐶₂𝑒^(𝑘)2+𝐶₃𝑒^(𝑘)3+𝐶₄𝑒^(𝑘)4+𝐶₅𝑒^(𝑘)5+𝐶₆𝑒^(𝑘)6+𝑂𝑒^(𝑘)7,𝐹^′𝑥^(𝑘)=𝐹^′(𝛼)𝐼+2𝐶₂𝑒^(𝑘)+3𝐶₃𝑒^(𝑘)2+4𝐶₄𝑒^(𝑘)3+5𝐶₅𝑒^(𝑘)4+6𝐶₆𝑒^(𝑘)5+𝑂𝑒^(𝑘)6. (7)
From (7), we calculate the expression of the inverse 𝐹^′𝑥^(𝑘)−1=𝐼+𝑋₂𝑒^(𝑘)+𝑋₃𝑒^(𝑘)2+𝑋₄𝑒^(𝑘)3+𝑋₅𝑒^(𝑘)4+𝑋₆𝑒^(𝑘)5×𝐹^′(𝛼)⁻¹+𝑂𝑒^(𝑘)6, (8) where 𝑋₂=−2𝐶₂, 𝑋₃=4𝐶²₂−3𝐶₃, 𝑋₄=6𝐶₃𝐶₂−6𝐶₂𝐶₃−8𝐶³₂−4𝐶₄, 𝑋₅=−5𝐶₅+8𝐶₂𝐶₄+8𝐶₄𝐶₂−12𝐶²₂𝐶₃−12𝐶₃𝐶²₂−12𝐶₂𝐶₃𝐶₂+9𝐶²₃+16𝐶⁴₂, and 𝑋₆=−6𝐶₆+10𝐶₂𝐶₅+10𝐶₅𝐶₂−16𝐶²₂𝐶₄−16𝐶₄𝐶²₂−16𝐶₂𝐶₄𝐶₂+12𝐶₃𝐶₄+12𝐶₄𝐶₃−18𝐶₃𝐶₂𝐶₃−18𝐶₂𝐶²₃−18𝐶²₃𝐶₂+24𝐶³₂𝐶₃+24𝐶²₂𝐶₃𝐶₂+24𝐶₂𝐶₃𝐶²₂+24𝐶₃𝐶³₂−32𝐶⁵₂.
These values have been obtained by imposing the conditions 𝐹^′𝑥^(𝑘)−1𝐹^′𝑥^(𝑘)=𝐹^′𝑥^(𝑘)𝐹^′𝑥^(𝑘)−1=𝐼. (9)
Then, the error expression in the first step of the method is 𝑦^(𝑘)−𝛼=𝐶₂𝑒^(𝑘)2+−2𝐶²₂+2𝐶₃𝑒^(𝑘)3+4𝐶³₂−4𝐶₂𝐶₃−3𝐶₃𝐶₂+4𝐶³₂𝑒^(𝑘)4+4𝐶₅−6𝐶₂𝐶₄−4𝐶₄𝐶₂+8𝐶²₂𝐶₃+6𝐶₃𝐶²₂+6𝐶₂𝐶₃𝐶₂−6𝐶²₃−8𝐶⁴₂𝑒^(𝑘)5+𝑂𝑒^(𝑘)6. (10)
Furthermore, we know that 𝐹𝑦^(𝑘)=𝐹^′(𝛼)𝑦^(𝑘)−𝛼+𝐶₂𝑦^(𝑘)−𝛼²+𝐶₃𝑦^(𝑘)−𝛼³+𝐶₄𝑦^(𝑘)−𝛼⁴+𝑂𝑦^(𝑘)−𝛼⁵ (11) and if we replace in (11) the powers of (𝑦^(𝑘)−𝛼), we obtain after some operations 𝐹𝑦^(𝑘)=𝐹^′(𝛼)𝐶₂𝑒^(𝑘)2+2−𝐶²₂+𝐶₃𝑒^(𝑘)3+3𝐶₄−4𝐶₂𝐶₃−3𝐶₃𝐶₂+5𝐶₂𝑒^(𝑘)4+4𝐶₅−6𝐶₂𝐶₄−4𝐶₄𝐶₂+8𝐶³₂𝐶₃+6𝐶₃𝐶²₂+6𝐶₂𝐶₃𝐶₂+2𝐶²₂𝐶₃+2𝐶₂𝐶₃𝐶₂−6𝐶²₃−12𝐶⁴₂𝑒^(𝑘)5+𝑂𝑒^(𝑘)6 (12) and also 𝐹^′𝑦^(𝑘)=𝐹^′(𝛼)𝐼+2𝐶²₂𝑒^(𝑘)2+4𝐶₂𝐶₃−𝐶³₂𝑒^(𝑘)3+6𝐶₂𝐶₄−8𝐶²₂𝐶₃−6𝐶₂𝐶₃𝐶₂+8𝐶⁴₂+3𝐶₃𝐶²₂𝑒^(𝑘)4+8𝐶₂𝐶₅−12𝐶²₂𝐶₄−8𝐶₂𝐶₄𝐶₂+16𝐶³₂𝐶₃+12𝐶₂𝐶₃𝐶²₂+12𝐶²₂𝐶₃𝐶₂−12𝐶₃𝐶³₂−12𝐶₂𝐶²₃+6𝐶₃𝐶₂𝐶₃+6𝐶²₃𝐶₂−16𝐶⁵₂𝑒^(𝑘)5+𝑂𝑒^(𝑘)6. (13)
In a similar way as in (8), 𝐹^′𝑦^(𝑘)−1=𝐼+𝑌₂𝑒^(𝑘)+𝑌₃𝑒^(𝑘)2+𝑌₄𝑒^(𝑘)3+𝑌₅𝑒^(𝑘)4+𝑌₆𝑒^(𝑘)5×𝐹^′(𝛼)⁻¹+𝑂𝑒^(𝑘)6, (14) where 𝑌₂=0,𝑌₃=−2𝐶²₂,𝑌₄=4(𝐶³₂−𝐶₂𝐶₃),𝑌₅=−6𝐶₂𝐶₄+8𝐶²₂𝐶₃+6𝐶₂𝐶₃𝐶₂−4𝐶⁴₂−3𝐶₃𝐶²₂, and 𝑌₆=−8𝐶₂𝐶₅+12𝐶²₂𝐶₄+8𝐶₂𝐶₄𝐶₂−16𝐶³₂𝐶₃−12𝐶₂𝐶₃𝐶²₂−12𝐶²₂𝐶₃𝐶₂+12𝐶₃𝐶³₂+8𝐶³₂𝐶₃+8𝐶₂𝐶₃𝐶²₂+12𝐶₂𝐶²₃−6𝐶₃𝐶₂𝐶₃−6𝐶²₃𝐶₂.
So, 𝐹^′𝑦^(𝑘)−1𝐹𝑥^(𝑘)=𝑒^(𝑘)+𝐶₂𝑒^(𝑘)2+−2𝐶²₂+𝐶₃𝑒^(𝑘)3+𝐶₄+2𝐶³₂−4𝐶₂𝐶₃𝑒^(𝑘)4+𝐶₅+6𝐶²₂𝐶₃+2𝐶₂𝐶₃𝐶₂−3𝐶₃𝐶²₂−6𝐶₂𝐶₄𝑒^(𝑘)5+𝑂𝑒^(𝑘)6. (15)
On the other hand, we get the Taylor expansion of 𝜇^(𝑘) by using (7) and (15): 𝜇^(𝑘)=𝐹^′𝑦^(𝑘)−1𝐹^′𝑥^(𝑘)=𝐼+2𝐶₂𝑒^(𝑘)+−2𝐶²₂+3𝐶₃𝑒^(𝑘)2+4𝐶₄−𝐶₂𝐶₃𝑒^(𝑘)3+5𝐶₅+2𝐶²₂𝐶₃−2𝐶₂𝐶₃𝐶₂−3𝐶₃𝐶²₂+4𝐶⁴₂−6𝐶₂𝐶₄𝑒^(𝑘)4+6𝐶₆+4𝐶²₂𝐶₄−4𝐶₂𝐶₄𝐶₂+4𝐶³₂𝐶₃+4𝐶²₂𝐶₃𝐶₂+8𝐶₂𝐶₃𝐶²₂+6𝐶₃𝐶³₂−6𝐶₃𝐶₂𝐶₃−6𝐶²₃𝐶₂−8𝐶⁵₂−8𝐶₂𝐶₅𝑒^(𝑘)5+𝑂𝑒^(𝑘)6. (16)
Let us note that 𝜇^(𝑘) tends to the identity matrix when 𝑥^(𝑘) and 𝑦^(𝑘) tend to 𝛼, thereby the second order polynomial approximation of the weight function, 𝐻(𝜇^(𝑘)), is 𝐻𝜇^(𝑘)≈𝐻(𝐼)+𝐻^′(𝐼)𝜇^(𝑘)+𝐻−𝐼^′′(𝐼)2𝜇^(𝑘)−𝐼². (17)
Let us consider now the following notation: 𝐻(𝐼)=𝐻₀,𝐻^′(𝐼)=𝐻₁,𝐻^′′(𝐼)2=𝐻₂. (18)
So, 𝑧^(𝑘)−𝛼=𝑦^(𝑘)−𝛼−𝐻𝜇^(𝑘)𝐹^′𝑦^(𝑘)−1𝐹𝑥^(𝑘)=−𝐻₀𝑒^(𝑘)+𝐼−𝐻₀−2𝐻₁𝐶₂𝑒^(𝑘)2+−2𝐼+2𝐻₀−4𝐻₂𝐶²₂+2𝐼−𝐻₀−3𝐻₁𝐶₃𝑒^(𝑘)3+3𝐼−𝐻₀−4𝐻₁𝐶₄+4𝐼−2𝐻₀+6𝐻₁+4𝐻₂𝐶³₂+−4𝐼+4𝐻₀+2𝐻₁−6𝐻₂𝐶₂𝐶₃+−3𝐼−3𝐻₁−6𝐻₂𝐶₃𝐶₂𝑒^(𝑘)4+4𝐼−𝐻₀−5𝐻₁𝐶₅+−6𝐼−8𝐻₂+6𝐻₀+4𝐻₁×𝐶₂𝐶₄+−4𝐼−8𝐻₂−4𝐻₁𝐶₄𝐶₂+8𝐼+10𝐻₂−6𝐻₀+8𝐻₁𝐶²₂𝐶₃+6𝐼+3𝐻₀+9𝐻₁𝐶₃𝐶²₂+6𝐼−2𝐻₀+2𝐻₂+6𝐻₁𝐶₂𝐶₃𝐶₂+−6𝐼−9𝐻₂−3𝐻₁𝐶²₃+−8𝐼+12𝐻₂−12𝐻₁𝐶⁴₂𝑒^(𝑘)5+𝑂𝑒^(𝑘)6. (19)
Then, 𝐹𝑧^(𝑘)=𝐹^′(𝛼)−𝐻₀𝑒^(𝑘)+𝐼−𝐻₀−2𝐻₁+𝐻²₀𝐶₂𝑒^(𝑘)2+−2𝐼+2𝐻²₀−4𝐻₂+4𝐻₀𝐻₁𝐶²₂+2𝐼−𝐻₀−3𝐻₁−𝐻³₀𝐶₃𝑒^(𝑘)3+3𝐼−𝐻₀−4𝐻₁+𝐻⁴₀𝐶₄+5𝐼+2𝐻₁+4𝐻₂−3𝐻²₀+8𝐻₀𝐻₂+4𝐻²₁+4𝐻₀𝐻₁𝐶³₂+−4𝐼+2𝐻₁−6𝐻₂+2𝐻²₀+6𝐻₀𝐻₁𝐶₂𝐶₃+−3𝐼−3𝐻₁−6𝐻₂+3𝐻²₀−6𝐻²₀𝐻₁−3𝐻³₀𝐶₃𝐶₂𝑒^(𝑘)4+4𝐼−𝐻₀−𝐻⁵₀−5𝐻₁𝐶₅+−6𝐼+4𝐻₁−8𝐻₂+2𝐻²₀+8𝐻₀𝐻₁𝐶₂𝐶₄+−4𝐼−4𝐻₁−8𝐻₂−4𝐻³₀+4𝐻⁴₀+8𝐻³₀𝐻₁×𝐶₄𝐶₂+10𝐼−𝐻₀+𝐻₁+10𝐻₂+𝐻₀𝐻₁+12𝐻₀𝐻₂−7𝐻²₀+6𝐻²₁𝐶²₂𝐶₃+6𝐼+9𝐻₁+12𝐻₀𝐻₁−12𝐻²₀𝐻₁−12𝐻₀𝐻²₁−12𝐻²₀𝐻₂+3𝐻³₀𝐶₃𝐶²₂+8𝐼+𝐻₀−𝐻₁+2𝐻₂+11𝐻₀𝐻₁+12𝐻₀𝐻₂+𝐻²₀+6𝐻²₁𝐶₂𝐶₃𝐶₂+−6𝐼+6𝐻²₀−3𝐻³₀−3𝐻₁−9𝐻²₀𝐻₁−9𝐻₂𝐶²₃+−12𝐼−4𝐻₁−20𝐻₀𝐻₁+4𝐻₂+16𝐻₁𝐻₂𝐶⁴₂𝑒^(𝑘)5+𝑂𝑒^(𝑘)6. (20)
Let us consider the truncated Taylor expansion of order two of the weight function 𝐺(𝜇^(𝑘)), 𝐺𝜇^(𝑘)≈𝐺(𝐼)+𝐺^′(𝐼)𝜇^(𝑘)+𝐺−𝐼^′′(𝐼)2𝜇^(𝑘)−𝐼², (21) and let us denote by 𝐺(𝐼)=𝐺₀,𝐺^′(𝐼)=𝐺₁,𝐺^′′(𝐼)2=𝐺₂. (22)
Finally, the error equation is expressed as 𝑒^(𝑘+1)=𝐻₀𝐼−𝐺₀𝑒^(𝑘)+𝐼−𝐻₀−𝐺₀+𝐺₀𝐻₀+2𝐺₁𝐻₀−𝐺₀𝐻²₀−2𝐻₁+2𝐺₀𝐻₁𝐶₂𝑒^(𝑘)2+−2𝐼+2𝐺₀−2𝐺₁+2𝐻₀−2𝐺₀𝐻₀+4𝐺₂𝐻₀−2𝐺₀𝐻²₀−2𝐺₁𝐻²₀+4𝐺₁𝐻₁−4𝐺₀𝐻₀𝐻₁−4𝐻₂+4𝐺₀𝐻₂𝐶²₂+2𝐼−𝐻₀−3𝐻₁−2𝐺₀+𝐺₀𝐻₀+3𝐺₀𝐻₁+𝐺₀𝐻³₀+3𝐺₁𝐻₀𝐶₃×𝑒^(𝑘)3+𝑂𝑒^(𝑘)4. (23)
If we choose 𝐻₀=0,𝐻₁=(1/2)𝐼, and 𝐻₂=0, we obtain 𝑒^(𝑘+1)=−2𝐼+2𝐺₀𝐶²₂+12𝐺𝐼−⁰2𝐶₃𝑒^(𝑘)3+𝑂𝑒^(𝑘)4. (24)
So, to increase the order of convergence up to four, the value of 𝐺₀ must be 𝐺₀=𝐼, and then, 𝑒^(𝑘+1)=𝐺₁−𝐶₃+4𝐶²₂𝑒^(𝑘)4+𝑂𝑒^(𝑘)5. (25)
Moreover, in order to reach fifth order of convergence, 𝐺₁ must be null. Therefore, the error equation is 𝑒^(𝑘+1)=𝐼−2𝐺₂𝐶²₂𝐶₃+−4𝐼+8𝐺₂𝐶⁴₂𝑒^(𝑘)5+𝑂𝑒^(𝑘)6. (26)
Finally, if 𝐺₂=(1/2)𝐼, we have 𝑒^(𝑘+1)=−32𝐶₃𝐶₂𝐶₃+14𝐶₂𝐶²₃+6𝐶₃𝐶³₂−𝐶₂𝐶₃𝐶²₂+𝐶³₂𝐶₃−4𝐶⁵₂𝑒^(𝑘)6+𝑂𝑒^(𝑘)7, (27) and the theorem is proved.

From the previous theorem, (5) defines a family of sixth-order methods. We can find different weight functions satisfying the conditions of the theorem. Specifically, we propose the following examples for the next sections.

Example 2. An element of our family of sixth-order is given by the weight functions 𝐻(𝑡)=12(),𝐺(𝑡)=[]𝑡−𝐼𝐼+𝑡⁻¹2𝐼−𝑡+𝑡², (28) which is called NAJC1.

Example 3. Another combination of weight functions that can be used is 𝐻(𝑡)=12(),𝐺(𝑡)1𝑡−𝐼=𝐼+2(𝑡−𝐼)². (29)
We will refer to this element of the class as NAJC2.

3. Numerical Results

In this section, we analyze the computational efficiency of our methods and compare them with other classical ones in the problem of preliminary orbit determination as well as in academic examples. The classical methods used are Newton’, Traub’, and Jarratt’s ones of convergence order 2, 3, and 4, respectively, whose iterative expressions are as follows.(i)Newton (N) 𝑥^(𝑘+1)=𝑥^(𝑘)−𝐹^′𝑥^(𝑘)−1𝐹𝑥^(𝑘). (30)(ii)Traub (T) 𝑦^(𝑘)=𝑥^(𝑘)−𝐹^′𝑥^(𝑘)−1𝐹𝑥^(𝑘),𝑥^(𝑘+1)=𝑦^(𝑘)−𝐹^′𝑥^(𝑘)−1𝐹𝑦^(𝑘). (31)(iii)Jarratt (J) 𝑧^(𝑘)=𝑥^(𝑘)−23𝐹^′𝑥^(𝑘)−1𝐹𝑥^(𝑘),𝑥^(𝑘+1)=𝑥^(𝑘)−123𝐹^′𝑧^(𝑘)−𝐹^′𝑥^(𝑘)−1×3𝐹^′𝑧^(𝑘)−𝐹^′𝑥^(𝑘)𝐹^′𝑥^(𝑘)−1𝐹𝑥^(𝑘). (32)

In our tests, we have used the following numerical settings: variable precision arithmetics of two hundred and fifty digits in Mathematica 8.0; moreover, in each iterative method, we have used the stopping criterium ‖𝐹(𝑥^(𝑘+1))‖+‖𝑥^(𝑘+1)−𝑥^(𝑘)‖<10⁻¹⁰⁰, and the approximated computational order of convergence 𝜌 (see [30]) has been calculated by using the formula: 𝑥𝑝≈𝜌=ln^(𝑘+1)−𝑥^(𝑘)/𝑥^(𝑘)−𝑥^(𝑘−1)ln𝑥^(𝑘)−𝑥^(𝑘−1)/𝑥^(𝑘−1)−𝑥^(𝑘−2). (33)

From this value, we have designed two practical indices to measure the computational efficiency: the approximated efficiency index 𝐼̃=𝜌^1/𝑑, (34) and the approximated computational index 𝐼_𝑐̃=𝜌^1/op, (35)𝑑 and op being the number of functional evaluations and the number of operations (products and quotients) per iteration, respectively.

In Tables 3, 4, 5, 6, and 7, we show the number of iterations, the previously defined indices, and the absolute errors committed between theoretical and practical values that we denote by 𝜖.

Two reference orbits have been used in the test for the preliminary orbit determination. The first can be found in [27], and the second one is a commercial real orbit called Tundra. As the orbital elements of each one are known, the vector positions (measured in Earth radius) at the instants 𝑡₁ and 𝑡₂ have been recalculated with 500 exact digits. These vector positions are 𝑟₁≈[],𝑟2.46080928705339,2.04052290636432,0.14381905768815₂≈[],1.98804155574820,2.50333354505224,0.31455350605251 (36) for Reference Orbit I and 𝑟₁≈[],𝑟−2.02862564034533,−0.74638890547506,−4.322222156844465₂≈[],4.24372000256074,−1.689387746496,6.79724893784587 (37) for Tundra Orbit. Then, our aim is to gain from these positions the orbital elements showed in Tables 1 and 2 with a precision as high as possible by means of proposed iterative schemes.

Table 1 

Orbital elements and temporal distance between 𝑟₁ and 𝑟₂. Reference Orbit I.

Table 2 

Orbital elements and temporal distance between 𝑟₁ and 𝑟₂. Tundra Orbit.

Table 3 

Results of reference Orbit I.

Table 4 

Results of tundra Orbit.

Table 5 

Results of system (a), 𝑥⁽⁰⁾ = (4,−3)^𝑇.

Table 6 

Results of system (b), 𝑥⁽⁰⁾ = (12,−2,−1)^𝑇.

Table 7 

Results of system (c), 𝑥⁽⁰⁾ = (5,5,5,−1)^𝑇.

If we look at the numerical results of Table 3, our methods NAJC1 and NAJC2 have the least number of iterations, the highest order of convergence with the highest efficiency index. For this case, the absolute error committed by Jarratt’s method is lower than in our procedures, which is due to use of an initial estimation very close to the solution and a very small temporal distance.

If we observe the numerical results of Tundra Orbit (Table 4), we note that the absolute error is stabilized, and we maintain the good results of the Reference Orbit I. We can conclude that working with the two equations provided by Gauss as a system, we improve the original procedure which has the restriction that the difference of true anomalies cannot be greater than 𝜋/4. We are able to increase the range of the difference of true anomalies associated to the observations until values close to 𝜋.

3.1. Other Nonlinear Problems

In order to continue checking the computational efficiency of the proposed schemes, NAJC1 and NAJC2 we use, in the following, some academic examples.(a)𝐹₁(𝑥)=(𝑓₁(𝑥),𝑓₂(𝑥))^𝑇: 𝑥=(𝑥₁,𝑥₂)^𝑇 and 𝑓_𝑖∶ℝ²→ℝ, 𝑖=1,2, 𝛼≈(3.47063096,−2.47063096), 𝑓₁(𝑥)=𝑒^𝑥1𝑒^𝑥2+𝑥₁cos𝑥₂,𝑓₂(𝑥)=𝑥₁+𝑥₂−1. (38)(b)𝐹₂(𝑥)=(𝑓₁(𝑥),𝑓₂(𝑥),𝑓₃(𝑥))^𝑇: 𝑥=(𝑥₁,𝑥₂,𝑥₃)^𝑇 and 𝑓_𝑖∶ℝ³→ℝ, 𝑖=1,2,3, 𝛼≈(2.14025,−2.09029,−0.223525), 𝑓₁(𝑥)=𝑥²₁+𝑥²₂+𝑥²₃𝑓−9,₂(𝑥)=𝑥₁𝑥₂𝑥₃𝑓−1,₃(𝑥)=𝑥₁+𝑥₂−𝑥²₃. (39)(c)𝐹₃(𝑥)=(𝑓₁(𝑥),𝑓₂(𝑥),𝑓₃(𝑥),𝑓₄(𝑥))^𝑇: 𝑥=(𝑥₁,𝑥₂,𝑥₃,𝑥₄)^𝑇 and 𝑓_𝑖∶ℝ⁴→ℝ, 𝑖=1,…,4, 𝛼=(±1/√3,±1/√3,±1/√3,±2/√3), 𝑓₁(𝑥)=𝑥₂𝑥₃+𝑥₄𝑥₂+𝑥₃,𝑓₂(𝑥)=𝑥₁𝑥₃+𝑥₄𝑥₁+𝑥₃,𝑓₃(𝑥)=𝑥₁𝑥₂+𝑥₄𝑥₁+𝑥₂,𝑓₄(𝑥)=𝑥₁𝑥₂+𝑥₁𝑥₃+𝑥₂𝑥₃−1. (40)

The results are showed in Tables 5, 6, and 7, where 𝜀_𝛼𝑖 denotes |𝑓_𝑖(𝛼)|. From these tables, the best schemes in terms of precision are NAJC1 and NAJC2, even when the initial estimation is far from the solution. In system (b), Traub’s method does not converge after five hundred iterations from this initial estimation.

4. Dynamical Planes

From the numerical results presented in Section 3, our proposed methods show to be competitive with respect to existing ones. Nevertheless, under the dynamical point of view, we have checked that they can have better behavior in terms of stability and wideness of the region of convergence.

For the representation of the convergence basins of our procedures and classical methods, we have used the software described in [31]. We draw a mesh with two thousand points per axis; each point of the mesh is a different initial estimation which we introduce in each procedure. If the method reaches the final solution in less than five hundred iterations, this point is drawn in orange. The color will be more intense when the number of iterations is lower. Otherwise, if the method arrives at the maximum of iterations, the point will be drawn in black. In each axis, we will represent each of the variables with which we work. The ratio sector triangle is represented in the abscissas and the difference of eccentric anomalies in the ordinates. In addition, we will use the reference orbit I, which is defined in Table 1, and the solution of the nonlinear system is in this case around (1, 0.1). For this reason, we choose [0,3]×[−1,1] as the region of representation.

In Figure 2, we show the dynamical planes of the classical methods. It can be observed that, in general, higher order means lower stability. If we focus our attention on the attraction basins of each plane, the method with the greatest stability is Newton, and the procedure with the lowest number of iterations is Jarratt.

(a) Newton

(b) Traub

(c) Jarratt

Figure 2 

Dynamical planes from classical methods and Reference Orbit I.

As we can see in Figure 3, both NAJC1 and NAJC2 have large areas of stability, similar to Newton’s one, but with order of convergence six. For the intensity of the orange in the attraction basins, the two schemes will have the least number of iterations. Moreover, if we compare both procedures, the attraction basins of NAJC1 are more disperse than the convergence basins of NAJC2 which makes the first one more unstable than the second one.

(a) NAJC1

(b) NAJC2

Figure 3 

Dynamical planes from new methods and Reference Orbit I.

5. Conclusions

The classical Gauss’ method for preliminary orbit determination has been improved, introducing a new performance by means of a nonlinear system. This fact increases the admissible spread of the observations (in order to ensure the convergence) from 𝜋/4 to 𝜋 and reduces the number of iterations of the process.

The new sixth-order methods NAJC1 and NAJC2 belonging to the class of methods designed by using matricial weight functions have good global properties of convergence and stability, even for initial estimations far from the solution.

It is a well-known fact that the size of the area of convergence is inversely proportional to the order of convergence. However, our methods hold a basin of attraction comparable with Newton’s one in spite of their sixth-order of convergence.

Acknowledgments

The authors thank the anonymous referees for their valuable comments and suggestions. This research was supported by Ministerio de Ciencia y Tecnología MTM2011-28636-C02-02.