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., [1619]). 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 , 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.

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 where

The angles , , , are the eccentric and true anomalies, respectively, associated to the observed positions and (let us denote by the modulus of vector , ).

Equation (1) is, actually, the composition of the first Gauss equation and the second Gauss equation where , 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 , 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 . 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: 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: , , and , , , being the identity matrix. Then, the scheme defined in (5) provides sixth order of convergence, whose error equation is given by where , , and .

Proof. If we use Taylor expansion of and around , we obtain
From (7), we calculate the expression of the inverse where , , , , and .
These values have been obtained by imposing the conditions
Then, the error expression in the first step of the method is
Furthermore, we know that and if we replace in (11) the powers of , we obtain after some operations and also
In a similar way as in (8), where , and .
On the other hand, we get the Taylor expansion of by using (7) and (15):
Let us note that tends to the identity matrix when and tend to , thereby the second order polynomial approximation of the weight function, , is
Let us consider now the following notation:
Let us consider the truncated Taylor expansion of order two of the weight function , and let us denote by
Finally, the error equation is expressed as
If we choose , and , we obtain
So, to increase the order of convergence up to four, the value of must be , and then,
Moreover, in order to reach fifth order of convergence, must be null. Therefore, the error equation is
Finally, if , we have 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 which is called NAJC1.

Example 3. Another combination of weight functions that can be used is
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) (ii)Traub (T) (iii)Jarratt (J)

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 , and the approximated computational order of convergence (see [30]) has been calculated by using the formula:

From this value, we have designed two practical indices to measure the computational efficiency: the approximated efficiency index and the approximated computational index and 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 for Reference Orbit I and 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.

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 , , , (b): and , , , (c): and , , ,

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 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.

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.

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 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.


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.