A modified Chebyshev Picard iteration method is proposed for solving orbit propagation initial/boundary value problems. Cosine sampling techniques, known as Chebyshev-Gauss-Lobatto (CGL) nodes, are used to reduce Runge’s phenomenon that plagues many series approximations. The key benefit of using the CGL data sampling is that the nodal points are distributed nonuniformly, with dense sampling at the beginning and ending times. This problem can be addressed by a nonlinear time transformation and/or by utilizing multiple time segments over an orbit. This paper suggests a method, called a multisegment method, to obtain accurate solutions overall regardless of initial states and albeit eccentricity by dividing the given orbit into two or more segments based on the true anomaly.

1. Introduction

A modified Chebyshev Picard iteration (MCPI) is an iterative numerical method for approximating solutions of linear or nonlinear ordinary differential equations to obtain time histories of system state trajectories [1, 2]. In contrast to many step-by-step integrators, the MCPI algorithm approximates long arcs of the state trajectory with an iterative path approximation approach and is ideally suited to parallel computation [3]. It is well known that Picard iteration has theoretical guarantees for converging to the solution assuming the forces are continuous, once differentiable, and the solution of the differential equation is unique [4]. The rate of convergence of Picard iteration is geometric rather than quadratic for Jacobian based methods. However, given a good starting approximation, excellent efficiency is possible, and the case for parallelization provides a significant advantage [5, 6].

Orthogonal Chebyshev polynomials are used as basis functions during each path iteration, and the integrations of Picard iteration are then performed analytically. The orthogonality of the Chebyshev basis functions implies that the least-square approximations can be computed to arbitrary precision without a matrix inversion; the coefficients are conveniently and robustly computed from discrete inner products [7]. Similar approximation approaches that use Legendre polynomials can be utilized, but the authors obtain slightly better results because the starting and ending points of the fits are not sampled as densely as the MCPI algorithm, and importantly the location of the nodes for the Chebyshev basis functions is computed exactly without iterations. The MCPI algorithm utilizes a vector-matrix framework for computational efficiency. Additionally, all Chebyshev coefficients and integrand function evaluations are independent, meaning that they can be simultaneously computed in parallel for further decreased computational costs [3].

For the MCPI algorithm, the cosine sampling techniques, known as Chebyshev-Gauss-Lobatto (CGL) nodes [8], are utilized to reduce Runge’s phenomenon. The Runge phenomenon is a problem of oscillation at the edges of an interval that occurs when using polynomial interpolation with polynomials of high degree [9]. Since dense sample points are distributed at the beginning and ending locations, less accurate solutions are usually obtained where sample points are more uniformly distributed [10].

For the most extreme counterexample, let us consider an unperturbed two-body problem, where the initial position is not located near the periapsis. Obviously, large errors can be observed near the periapsis due to sparse sample points where the dynamics are most nonlinear, yet we waste the dense sample points at apoapsis when the problem is most linear as shown in Figure 1.

This problem is overcome by introducing a multisegment method and the results are compared with the basic MCPI algorithm. This paper only considers two and three segments per one orbit. The performance of the proposed approach is established by numerical examples of the two-body problem.

2. Modified Chebyshev Picard Iteration

The MCPI algorithm combines the discoveries of two great mathematicians: Émile Picard (Picard iteration) and Rafnuty Chebyshev (Chebyshev polynomials). Combing these techniques was first proposed by Clenshaw and Norton in 1963 [11].

Picard stated that any first-order differential equation with an initial condition can be rearranged without approximation as follows:

In the MCPI algorithm, orthogonal Chebyshev polynomials are used as basis functions to approximate the integrand in the Picard integral. Chebyshev polynomials reside in the domain and are defined recursively as

The system dynamics are normalized such that the time span of integration is projected onto the domain of the Chebyshev polynomials. Thus, the system states are approximated using the Chebyshev polynomial basis functions. The orthogonal nature of the basis functions means the coefficients that linearly scale the basis functions are computed independently as simple ratios of inner products without requiring matrix inversions.

A key feature of the MCPI algorithm is a nonuniform cosine density sampling of the domain of the Chebyshev basis functions, called CGL nodes, defined as follows:

This sampling scheme provides much higher density towards the edges (beginning and ending points), which enables high accuracy solutions near the boundaries of the state trajectory. This scheme eliminates the Runge phenomenon, a common issue in function approximations, whereby noisy estimates are returned near the edges due to lack of knowledge of the states on the other sides of the boundaries (see Figure 2). The coefficients multiplying the Chebyshev basis functions are approximated by the method of least squares, which generally requires a matrix inversion. A wonderful side effect of the cosine sampling scheme is that the matrix required to be inverted in the normal equations of least squares is diagonal; thus the inverse is trivial.

A full derivation of the MCPI algorithm is not included in this work (refer to Bai [3]). Instead, the authors present a flowchart in Figure 3 briefly summarizing the mathematics underlying the MCPI algorithm for solution of initial value problems.

3. Multisegment Approach for MCPI Algorithm

This work considers an unperturbed two-body problem, where the initial position is not located near the periapsis. As expected, large errors are observed near the periapsis where dense sample points are required, but sparse sample points are distributed. In addition, even though initial positions are located near the periapsis, accurate solutions cannot be obtained for highly elliptical orbits. To obtain accurate solutions for the above cases using the MCPI algorithm, the multisegment approach is proposed.

Given the initial true anomaly , two and three segmented orbits are considered as shown in Figure 4. These two cases require patch times to link the divided segments. To distribute dense sample points near the periapsis, several strategies are presented.

For two segmented orbits, the time for the patch point is selected as the time at the perigee, where the true anomaly is degrees. For three segmented orbits, the time for the first patch point is selected where degree for symmetry, and the time for the second patch point is selected where degrees. To find propagation times for each segment, the following calculation needs to be performed. First, given the initial position and velocity vectors, prescribe the break point and one orbit period time as follows [12]:

Second, calculate the initial mean anomaly as follows [13]: where is the eccentricity and the eccentric anomaly () is defined as [13]

Finally, the propagation times for each segment are calculated as follows: where is the number of segments for the orbit.

The sets of propagation times are determined as follows:

For more than three segments and the associated break points, the above logic is readily extended. The optimization of the break points to achieve efficiency and accuracy is not addressed in this paper but is an easy-to-pose optimization problem research for a future study.

4. Numerical Examples

A satellite motion integration problem, idealized for the case with only the inverse square gravitational force from the Earth, is considered. The three-dimensional dynamical equations are given by [12] where , , and are the three coordinates in Earth-centered inertial reference frame; is the distance of the satellite from the Earth; is the Earth gravitational constant and is chosen as .

To verify the results, the following normalized energy error check is utilized: where is the initial energy and the energy is calculated as follows: Note that the goal for the demonstration example in this paper is to obtain solutions where . Moreover, for the unperturbed two-body problem, the analytical solution [12] for the function can also be used to confirm the accuracy of the solution.

Two sets of initial position and velocity vectors are given as follows: The given initial states lead to and degrees, respectively, and the classical orbital elements are listed in Table 1.

For the MCPI algorithm implementation to solve this problem, various tuning parameters are determined in a prior calculation: maximum iteration number , error tolerance , degree of polynomial , and number of sample points . This work focuses on finding a methodology to improve MCPI accuracy and reduce computational burden given the described factors in Table 2 and initial conditions listed in Table 1.

Numerical simulations are performed, and the normalized energy error results are shown in Figures 58. Figure 5 shows that the normalized energy errors are much larger than the requirement (). Obviously, the largest error is observed at the periapsis when degrees because of sparse sample point distributions at the periapsis.

Figure 6 shows that the solution satisfies the requirement when degrees using the two-segment scheme. The same number of sample points is distributed for each segment, and the total number of the sample points is equal to the number of sample points for the basic (one-segment) MCPI algorithm. Note that only the two-segment orbit approach is used when the initial position is located at apoapsis for symmetry.

Figure 7 shows that the solution satisfies the requirement when degrees using the three-segmentation scheme. The same number of sample points is distributed for the second and third segments, and the total number of the sample points is equal to the number of sample points for the basic MCPI algorithm.

For the case where degrees, both approaches such as the two- and three-segment schemes are applicable. As shown in Figure 8, both approaches satisfy the requirement, but the three-segment orbit approach outperforms the other methods. The number of nodes for each approach is determined by a heuristic method for this paper (and tuned numerically); and a methodology to select optimal number of nodes is under development.

5. Conclusion

The modified Chebyshev Picard iteration (MCPI) algorithm uses Chebyshev-Gauss-Lobatto (CGL) nodes to reduce the Runge phenomenon. By using the CGL nodes, however, less accurate solutions may be obtained where sparse sample points are distributed. Physical insights indicate that the dense nodes should be located where the orbit is most nonlinear. However, the stating epoch state can be at a random point in the orbit. For the unperturbed two-body problem, where the initial state is not located near the periapsis and the eccentricity is high, the multisegment approach is utilized to obtain an accurate solution. The final perigee passage can be used to make all subsequent segment breaks symmetrical about the major axis. As a result, the multisegment approach provides much more accurate solutions when compared to the solution from the basic MCPI algorithm with random user-specified segmentation logic. Moreover, it is shown that the three-segment orbit approach outperforms others in terms of computational efficiency. To improve the performance of the MCPI algorithm, this approach will be very useful, especially when the initial position is not located near the periapsis and high eccentric orbits are given.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.


This work is supported by the Air Force Office of Scientific Research, USA.