Characterization of Numerical Error in the Simulation of Translunar Trajectories
Using the Method of Nearby Problems
by
Ashish Ashok Jagat
A thesis submitted to the Graduate Faculty of
Auburn University
in partial fulfillment of the
requirements for the Degree of
Master of Science
Auburn, Alabama
May 9, 2011
Copyright 2011 by Ashish Ashok Jagat
Approved by
Andrew J. Sinclair, Chair, Associate Professor of Aerospace Engineering
David Cicci, Professor of Aerospace Engineering
Andrew Shelton, Assistant Professor of Aerospace Engineering
ii
Abstract
This thesis focuses on analyzing the effect numerical error has on the accurate
simulation of translunar trajectories. The method of nearby problems is employed to
estimate the numerical error. A simulation is developed to generate translunar
trajectories. Analytical curve fit is generated to this numerical solution and this curve fit
is used to compute analytical source terms. The addition of these source terms to the
governing equations defines a nearby problem, for which the curve fit serves as an exact
solution. By solving the nearby problem numerically, the numerical error in it can be
calculated. This facilitates the estimation of numerical error in the original problem.
iii
Acknowledgments
This work could not have been completed without the help and support of many
individuals. I am especially grateful to my advisor, Dr. Andrew Sinclair, for his constant
encouragement, guidance, and support. He went above and beyond to help me succeed.
His suggestions from time to time have been imperative to this work. Working with him
has been a tremendous learning experience for me. I am very thankful for the insights
and valuable time of my committee members Dr. David Cicci and Dr. Andrew Shelton. I
would also like to thank Mr. Dave Patrick of the Department of Physics at Auburn
University for giving me the opportunity to work as a teaching assistant. I can?t imagine
life without friends and have no words for their tremendous support both during good and
tough times. Finally, the support of my parents, my sister Anagha, and the rest of my
family in whatever I do is invaluable to me.
iv
Table of Contents
Abstract ............................................................................................................................... ii
Acknowledgments.............................................................................................................. iii
List of Tables ..................................................................................................................... vi
List of Figures ................................................................................................................... vii
1. Introduction .....................................................................................................................1
2. Original Problem .............................................................................................................7
3. Numerical Solution to the Original Problem ................................................................14
4. Curve Fitting Methods ..................................................................................................19
Least Squares .........................................................................................................19
Cubic Splines .........................................................................................................28
Fifth-degree Hermite Splines .................................................................................38
5. Generating Analytical Source Terms ............................................................................51
Using Least Squares ...............................................................................................52
Using Cubic Splines ...............................................................................................56
Using Fifth-degree Hermite Splines ......................................................................61
6. Nearby Problem to the Original Problem .....................................................................66
7. Conclusion and Recommendations ...............................................................................76
References .........................................................................................................................78
v
Appendices
A. MATLAB Code to Numerically Solve the Original Problem ........................81
B. Subroutines Used in the Code Given in Appendix A .....................................84
C. MATLAB Code to Calculate the Coefficients of Cubic Splines Curve Fit .....85
D. MATLAB Code to Calculate the Coefficients of
Fifth-degree Hermite Splines Curve Fit ..........................................................87
E. MATLAB Code to Calculate the Coefficients of
Multi-resolution Fifth-degree Hermite Splines Curve Fit ..............................90
F. Subroutines Used in the Code Given in Appendix E ......................................92
G. Tridiagonal Solver ..........................................................................................93
H. MATLAB Code to Numerically Solve the Nearby Problem ..........................94
I. Subroutines Used in the Code Given in Appendix H .....................................97
vi
List of Tables
1. LEO Parameters ...........................................................................................................13
2. Region Wise Number of Spline Zones Using Multi-resolution
Fifth-degree Hermite Splines Curve Fit to X ...............................................................40
3. Region Wise Number of Spline Zones Using Multi-resolution
Fifth-degree Hermite Splines Curve Fit to Y ...............................................................41
vii
List of Figures
1. n-body Problem ................................................................................................................8
2. Three-body Problem ........................................................................................................9
3. System Model ...............................................................................................................11
4. The Earth, LEO, Point of ?v and Initial Trajectory .......................................................13
5. Spacecraft and Moon Trajectories ................................................................................15
6. Time vs. X Position .......................................................................................................16
7. Time vs. Y Position .......................................................................................................16
8. Time vs. X Velocity ......................................................................................................17
9. Time vs. Y Velocity .......................................................................................................17
10. Time vs. X Acceleration ..............................................................................................18
11. Time vs. Y Acceleration ..............................................................................................18
12. Curve Fit to X - Least Squares Using 3rd Degree Polynomial ....................................22
13. Curve Fit Residuals for X - Least Squares Using 3rd Degree Polynomial ...................22
14. Curve Fit to X - Least Squares Using 5th Degree Polynomial ....................................23
15. Curve Fit Residuals for X - Least Squares Using 5th Degree Polynomial ...................23
16. Curve Fit to X - Least Squares Using 20th Degree Polynomial ..................................24
17. Curve Fit Residuals for X - Least Squares Using 20th Degree Polynomial .................24
18. Curve Fit to Y - Least Squares Using 3rd Degree Polynomial .....................................25
19. Curve Fit Residuals for Y - Least Squares Using 3rd Degree Polynomial ...................25
viii
20. Curve Fit to Y - Least Squares Using 5th Degree Polynomial .....................................26
21. Curve Fit Residuals for Y - Least Squares Using 5th Degree Polynomial ....................26
22. Curve Fit to Y - Least Squares Using 20th Degree Polynomial ...................................27
23. Curve Fit Residuals for Y - Least Squares Using 20th Degree Polynomial ..................27
24. Schematic of Cubic Splines Interpolation ...................................................................28
25. Curve Fit to X - Cubic Splines Using 8 Spline Zones ...............................................30
26. Curve Fit Residuals for X - Cubic Splines Using 8 Spline Zones ..............................30
27. Curve Fit to X - Cubic Splines Using 63 Spline Zones .............................................31
28. Curve Fit Residuals for X - Cubic Splines Using 63 Spline Zones ............................31
29. Curve Fit to X - Cubic Splines Using 504 Spline Zones ............................................32
30. Curve Fit Residuals for X - Cubic Splines Using 504 Spline Zones ..........................32
31. Curve Fit to X - Cubic Splines Using 5040 Spline Zones .........................................33
32. Curve Fit Residuals for X - Cubic Splines Using 5040 Spline Zones ........................33
33. Curve Fit to Y - Cubic Splines Using 8 Spline Zones ................................................34
34. Curve Fit Residuals for Y - Cubic Splines Using 8 Spline Zones ...............................34
35. Curve Fit to Y - Cubic Splines Using 63 Spline Zones ..............................................35
36. Curve Fit Residuals for Y - Cubic Splines Using 63 Spline Zones .............................35
37. Curve Fit to Y - Cubic Splines Using 504 Spline Zones ............................................36
38. Curve Fit Residuals for Y - Cubic Splines Using 504 Spline Zones ...........................36
39. Curve Fit to Y - Cubic Splines Using 5040 Spline Zones ..........................................37
40. Curve Fit Residuals for Y - Cubic Splines Using 5040 Spline Zones .........................37
41. Schematic of Hermite Splines Interpolation ...............................................................38
42. Curve Fit to X - Fifth-degree Hermite Splines Using 8 Spline Zones .......................42
ix
43. Curve Fit Residuals for X - Fifth-degree Hermite Splines Using 8 Spline Zones ......42
44. Curve Fit to X - Fifth-degree Hermite Splines Using 63 Spline Zones .....................43
45. Curve Fit Residuals for X - Fifth-degree Hermite Splines Using 63 Spline Zones ....43
46. Curve Fit to X - Fifth-degree Hermite Splines Using 504 Spline Zones ...................44
47. Curve Fit Residuals for X - Fifth-degree Hermite Splines
Using 504 Spline Zones ..............................................................................................44
48. Curve Fit to X - Fifth-degree Hermite Splines Using 5040 Spline Zones .................45
49. Curve Fit Residuals for X - Fifth-degree Hermite Splines Using
5040 Spline Zones.......................................................................................................45
50. Curve Fit to Y - Fifth-degree Hermite Splines Using 8 Spline Zones .......................46
51. Curve Fit Residuals for Y - Fifth-degree Hermite Splines Using 8 Spline Zones ......46
52. Curve Fit to Y - Fifth-degree Hermite Splines Using 63 Spline Zones .....................47
53. Curve Fit Residuals for Y - Fifth-degree Hermite Splines Using 63 Spline Zones ....47
54. Curve Fit to Y - Fifth-degree Hermite Splines Using 504 Spline Zones ...................48
55. Curve Fit Residuals for Y - Fifth-degree Hermite Splines
Using 504 Spline Zones ..............................................................................................48
56. Curve Fit to Y - Fifth-degree Hermite Splines Using 5040 Spline Zones .................49
57. Curve Fit Residuals for Y - Fifth-degree Hermite Splines
Using 5040 Spline Zones ............................................................................................49
58. Curve Fit Residuals for X - Multi-resolution Fifth-degree Hermite Splines Fit .........50
59. Curve Fit Residuals for Y - Multi-resolution Fifth-degree Hermite Splines Fit .........50
60. Source Term Histories for X - Least Squares Using 3rd Degree Polynomial ..............53
61. Source Term Histories for X - Least Squares Using 5th Degree Polynomial ..............53
62. Source Term Histories for X - Least Squares Using 20th Degree Polynomial ............54
63. Source Term Histories for Y - Least Squares Using 3rd Degree Polynomial ..............54
x
64. Source Term Histories for Y - Least Squares Using 5th Degree Polynomial ..............55
65. Source Term Histories for Y - Least Squares Using 20th Degree Polynomial ............55
66. Source Term Histories for X - Cubic Splines Using 8 Spline Zones ..........................57
67. Source Term Histories for X - Cubic Splines Using 63 Spline Zones ........................57
68. Source Term Histories for X - Cubic Splines Using 504 Spline Zones ......................58
69. Source Term Histories for X - Cubic Splines Using 5040 Spline Zones ....................58
70. Source Term Histories for Y - Cubic Splines Using 8 Spline Zones ..........................59
71. Source Term Histories for Y - Cubic Splines Using 63 Spline Zones ........................59
72. Source Term Histories for Y - Cubic Splines Using 504 Spline Zones ......................60
73. Source Term Histories for Y - Cubic Splines Using 5040 Spline Zones ....................60
74. Source Term Histories for X - Fifth-degree Hermite Splines
Using 8 Spline Zones ..................................................................................................62
75. Source Term Histories for X - Fifth-degree Hermite Splines
Using 63 Spline Zones ................................................................................................62
76. Source Term Histories for X - Fifth-degree Hermite Splines
Using 504 Spline Zones ..............................................................................................63
77. Source Term Histories for X - Fifth-degree Hermite Splines
Using 5040 Spline Zones ............................................................................................63
78. Source Term Histories for Y - Fifth-degree Hermite Splines
Using 8 Spline Zones ..................................................................................................64
79. Source Term Histories for Y - Fifth-degree Hermite Splines
Using 63 Spline Zones ................................................................................................64
80. Source Term Histories for Y - Fifth-degree Hermite Splines
Using 504 Spline Zones ..............................................................................................65
81. Source Term Histories for Y - Fifth-degree Hermite Splines
Using 5040 Spline Zones ............................................................................................65
82. Nearby Problem to X - Using Fifth-degree Hermite Splines
with 8 Spline Zones ....................................................................................................67
xi
83. Nearby Problem to X - Using Fifth-degree Hermite Splines
with 63 Spline Zones ..................................................................................................67
84. Nearby Problem to X - Using Fifth-degree Hermite Splines
with 504 Spline Zones ................................................................................................68
85. Nearby Problem to X - Using Fifth-degree Hermite Splines
with 5040 Spline Zones ..............................................................................................68
86. Nearby Problem to Y - Using Fifth-degree Hermite Splines
with 8 Spline Zones ....................................................................................................69
87. Nearby Problem to Y - Using Fifth-degree Hermite Splines
with 63 Spline Zones ..................................................................................................69
88. Nearby Problem to Y - Using Fifth-degree Hermite Splines
with 504 Spline Zones ................................................................................................70
89. Nearby Problem to Y - Using Fifth-degree Hermite Splines
with 5040 Spline Zones ..............................................................................................70
90. Exact Error in Nearby Problem to X Using Fifth-degree Hermite Splines .................71
91. Exact Error in Nearby Problem to Y Using Fifth-degree Hermite Splines .................72
92. Observed Order of Accuracies for Various
Numerical Solutions to the Nearby Problem Using Various Step Sizes ....................73
93. Numerical Error Estimates for X Using Various Methods .........................................75
94. Numerical Error Estimates for Y Using Various Methods ..........................................75
1
CHAPTER 1
INTRODUCTION
Numerical simulations play an important role in mathematical modeling of many
systems in engineering. They are of utmost importance when it comes to estimating the
performance of systems too complex for analytical solutions. Numerical simulations are
imperative to the field of orbital mechanics as there are many differential equations and
dynamical systems which cannot be solved analytically. Numerical simulations in this
field have been used since the era of the Apollo missions.
In numerical simulations, however, one must account for the numerical error as it
affects the efficacy of the simulation scheme. Previous studies on the numerical error in
orbital mechanics simulations have included three types of numerical error - iteration
error, round-off error, and discretization error. Iteration error is the difference between
the current iterative solution and the exact solution. In Runge-Kutta simulation of
translunar trajectories, however, iterative solutions are not used, and hence iteration error
is not considered in this study. Round-off error is caused by the fact that digital
computers can store numbers with only a finite precision. Discretization error is the
difference between the solution of the discretized equation and the exact solution of the
original differential equation.
Many researchers have studied numerical error in orbital mechanics simulations.
Some of them have focused on the development of accuracy assessment techniques.
2
Huang and Innanen [1] showed that the traditional ways of checking the accuracy of the
numerical solutions to the dynamical systems by the use of known integrals or the
integral invariant relations are neither exact nor reliable due to the tendency of these
numerical solutions to keep the integrals constant. They suggested a revised technique to
use the integral invariant relations for checking the accuracy of the numerical solutions to
the dynamical systems.
Other researchers have focused on surveying the accuracy of various numerical
integration schemes, often assessing the accuracy by comparison of the numerical
solutions to known exact solutions. Fox [2] did an accuracy based comparative study of
various categories of numerical integration methods applied to the solution of two-body
problem. Berry and Healy [3] compared the efficacy of various accuracy assessment
techniques of numerical integrators using two-body problem with and without
perturbations. Montenbruck [4] assessed the usefulness of various methods of numerical
integration such as Runge-Kutta, multi-step and extrapolation based methods for
generating numerical solutions to the problems involving solar system bodies or artificial
satellites. Hadjifotinou and Gousidou-Koutita [5] proposed a new method called the
recurrent power series (RPS) method for the integration of the system of n satellites
orbiting a point-mass planet.
One approach taken to check the numerical accuracy of a numerical solution is to
construct a similar problem to the original problem of interest. In this approach, the
similar problem has an exact solution. Therefore, error in the numerical solution to the
similar problem is exactly known which is then used to estimate the error in the
numerical solution to the original problem of interest. Researchers studying orbital
3
mechanics as well as fluid dynamics have taken this approach. Roach [6] proposed
?method of manufactured solutions (MMS)?. MMS involves manufacturing an exact
solution to a set of equations which are a modified form of the original differential
equations. The solution obtained to this set of modified equations may not have physical
significance. Therefore, MMS is used only to verify the mathematics involved in solving
the original equations, and does not verify the solution obtained by solving the original
equations.
Other researchers have used their similar problem to actually validate the
numerical solution to the original problem. All these researchers use curve fit to the
numerical solution of the original problem to construct the similar problem. The way in
which these curve fits are calculated and are then used to construct the similar problem
might differ slightly in each case. Zadunaisky [7, 8, 9, 10] suggested a technique in
which he called his similar problem as a ?pseudo-system? and applied it to the problems
in orbital mechanics. Junkins and Lee [11] constructed ?benchmark problem? for hybrid
coordinate systems of ordinary/partial differential equations. Hopkins and Roy [12, 13]
referred to their similar problem as ?nearby problem? and the approach was called as ?the
method of nearby problems (MNP)?. They applied this method to the problems in fluid
dynamics.
MNP is based on constructing a problem near the original problem of interest.
This nearby problem is constructed in such a way that it is both representative of the
original problem and also has an exact known solution. This nearby problem is then
solved numerically using the same numerical solution scheme that was used to
numerically solve the original problem. Because the exact solution to the nearby problem
4
is known, error in its numerical solution can be calculated. This information is then used
to estimate the error in the numerical solution of the original problem. MNP involves
five steps. These steps are explained as follows:
Establishing an Accurate Numerical Solution to the Original Problem
Once the problem of interest is identified, the first step is to discretize this
problem and produce an accurate numerical solution.
Generating an Analytical Curve Fit to the Above Numerical Solution
Once the accurate numerical solution is computed in previous step, this step
involves generating an analytical curve fit to this numerical solution. One of the many
curve fitting techniques is used to generate this curve fit. It should be kept in mind that
the technique used for curve fitting should provide a particular order of continuity which
is problem dependent. Once the curve fit is generated, it should be examined to see how
good the fit approximates the numerical solution. This analytical curve fit will serve as
the exact solution to the nearby problem.
Generating Analytical Source Terms
The nearby problem differs from the original problem by (hopefully) small source
terms. These source terms are obtained by operating the original equation on the analytic
curve fit obtained from the previous step. In the limit, as the magnitude of the source
terms approaches zero, the nearby problem approaches the original problem. The
nearness of the nearby problem to the original problem can be judged by the magnitude
of the source terms.
5
Numerically Solving the Nearby Problem
The nearby problem consists of original equations plus the analytical source
terms. This step involves solving the nearby problem numerically using the same
numerical solution scheme that was used to solve the original problem.
Estimating the Numerical Error in the Original Problem
Because both the exact and the numerical solution to the nearby problem are
known, error in the numerical solution can be calculated for the nearby problem. This
information can then be used to estimate the error in the numerical solution to the original
problem.
In this thesis, an effort is made to extend the application of MNP to the problems
in orbital mechanics. The objective is to demonstrate the usefulness of MNP in
validating the accuracy of the numerical solutions to the problems in orbital mechanics
by constructing a nearby problem to the Earth-spacecraft-Moon three-body problem.
While this work also uses the curve fit to the numerical solution of the original problem
to construct the nearby problem, unlike Zadunaisky, various curve fitting techniques are
examined first and then the technique satisfying certain criteria is used to construct the
nearby problem. Moreover, the way in which this curve fit is used to construct the
nearby problem differs from that of Zadunaisky?s.
Chapter 2 discusses the original problem studied in this work. In this chapter,
equations of motion of the n-body problem are derived first and these equations are then
specialized for the case of the three-body problem. System model for the Earth-
spacecraft-Moon three-body problem is also discussed and the equations of motion
governing the motion of the spacecraft in the cis-lunar space are formed. Chapter 3
6
discusses the numerical solution to the Earth-spacecraft-Moon three-body problem.
Chapter 4 discusses various curve fitting techniques and their feasibility to construct the
nearby problem. Chapter 5 discusses the calculation of analytical source terms using the
various curve fitting techniques studied in chapter 4. In chapter 6, the nearby problem is
constructed and the nearness of this nearby problem to the original problem is
established. This allows the exact error in the numerical solution to the nearby problem
to be considered as a good estimation of the error in the numerical solution to the original
problem. The numerical error estimated by MNP is compared to that estimated by
Richardson extrapolation using both global and local order of accuracy. In chapter 7,
a conclusion to this thesis is presented.
7
CHAPTER 2
ORIGINAL PROBLEM
The original problem studied in this work is the Earth-spacecraft-Moon three-
body problem. One of the most fundamental problems of orbital mechanics is to
accurately describe the motion of n gravitationally interacting massive particles also
known as the n-body problem. In this chapter, we will begin with the general n-body
problem and then we will specialize to our three-body problem. Newton introduced his
three laws of motion in The Mathematical Principles of General Philosophy, or, more
simply, the Principia, in 1687. Newton?s law of universal gravitation along with the
second law of motion can be used to describe the n-body problem. The second law can
be expressed mathematically in vector notation as follows:
? (1)
In Eq. (1), is the vector sum of all forces acting on the mass and ? is the vector
acceleration of the mass measured relative to an inertial reference frame. Similarly, the
universal law of gravitation can be expressed mathematically in vector notation as
follows:
(2)
In Eq. (2), is the gravitational force exerted on mass by mass is the vector
from to and is the universal gravitational constant.
8
Z
mj
m1 F1j
rj F2j Fnj
r1
Y
r2 rn
mn
X m2
Figure 1. n-body Problem [14]
The n-body problem is illustrated in figure 1. It shows a system of n
bodies . Let us assume an inertial reference frame (X, Y, Z) in which
the position vectors of these n masses are respectively. We wish to study
the motion of one of these bodies. Let us call this body as the jth body, . At any given
time in its journey, this body is being acted upon by several gravitational masses and may
be experiencing other forces such as drag, thrust, and solar radiation pressure. Let us not
consider all these other forces for the time being. Applying Newton?s law of universal
gravitation to the above system, the force exerted on by is:
(3)
In Eq. (3), .
9
The vector sum, , of all such gravitational forces acting on the jth body may be written
as:
? (
)
(4)
Applying Newton?s second law of motion,
? (5)
Dividing both sides of Eq. (5) by we get,
? ? (
)
(6)
Eq. (6) describes the inertial motion of the jth body. For our purposes, we will be
interested in describing the motion of a body relative to another body. Consider figure 2.
Figure 2. Three-body Problem [16]
10
Equation of motion of m1 with respect to the inertial reference frame can be given by:
? (
) (
) (7)
Also the equation of motion of m2 with respect to the same reference frame is:
? (
) (
) (8)
By subtracting Eq. (7) from (8), equation of motion of relative to is given as:
? ? (
)
(9)
Let as shown in figure 2.
Thus, Eq. (9) can be simplified as:
? ( ) ( ) (10)
Similarly, equation of motion of relative to is given as:
? ( ) ( ) (11)
Considering and above equations can be rewritten as:
? ( ) (12)
? ( ) (13)
In both the above equations, the first term on the right hand side is the central
acceleration due to the primary body, while the second term on the right hand side is the
11
disturbing acceleration from the perturbing body. As mentioned earlier, several external
forces such as drag, thrust, and solar radiation pressure are not considered while deriving
these equations of motion. Also, Newton's law of gravitation applies only if the bodies
are spherical and the mass is evenly distributed in spherical shells. Thus, the non-
spherical shape of the bodies also results in a perturbing force. Therefore, the governing
equations derived above are approximations and inclusion of all the perturbing forces
could facilitate a more accurate modeling of the trajectory.
S/C
r d
E M
?
Figure 3. System Model
For the Earth-spacecraft-Moon three-body problem, let m1 be the Earth, m2 be the
spacecraft and m3 be the moon. Figure 3 shows the system model for this problem. The
spacecraft motion is described either in the Earth-Centered (EC) reference frame or in the
Moon-Centered (MC) reference frame. Both EC and MC are non-rotating reference
frames. The EC frame has its origin at the center of mass of the Earth while the MC
frame has its origin at the center of mass of the Moon. The position vector of the
spacecraft in the EC frame is given by r while in the MC frame it is given by d. The
equation of motion of the spacecraft in vector form in the EC frame is given by Eq. (12)
rewritten below. Here, is the gravitational parameter for the Earth with a value of
3.98600436 ? 1014 m3/s2 and is the gravitational parameter for the Moon with a value
of 4.90266 ? 1012 m3/s2.
12
???????? ???? 33322dd ??? ?drr drt me (14)
In this equation, the Earth is the central body and the Moon causes the perturbative
acceleration. The Cartesian components of r in a non-rotating frame are written as
X and Y, and the scalar components of the equation of motion are as follows.
???????? ???? 33322dd ???? xmme dXXrtX
(15)
???????? ???? 33322dd ???? ymme dYYrt Y
(16)
Similarly, the equation of motion of the spacecraft in vector form in the MC frame is
given by Eq. (13) rewritten below.
???????? ???? 33322dd ??? ?rdd rdt em (17)
In this equation, the Moon is the central body and the Earth causes the perturbative
acceleration. The Cartesian components of d in a non-rotating frame are written as
Xm and Ym, and the scalar components of the equation of motion are as follows.
???????? ???? 33322dd ???? xemmm rXXdtX (18)
???????? ???? 33322dd ???? yemmm rYYdtY
(19)
13
These two sets of equations in the EC and the MC reference frames represent analytically
equivalent description of the three-body problem. The numerical solution described in
the next chapter will use switching between the two reference frames.
The problem description is completed by defining various initial conditions. The
starting location of the Moon is out from the Earth along the positive X-axis as shown in
figure 3. Because the mean eccentricity of the Moon?s orbit is only about 0.0549, it is
considered to be a circular orbit with a radius (?) of 384,400 km. The spacecraft
trajectory is assumed to be coplanar with the Moon?s orbit. The initial conditions for the
translunar trajectory are specified by the injection criteria relative to the low Earth orbit
(LEO): altitude (alt), ?v, and angle (?). These initial conditions are given in table 1 and
are illustrated in figure 4. [16]
Parameter Value(Units)
re 6372.797(km)
alt 359750(m)
?v 3102.13(m/s)
? -36.890(degrees)
Table 1. LEO Conditions
Figure 4. The Earth, LEO, Point of ?v and Initial Trajectory
14
CHAPTER 3
NUMERICAL SOLUTION TO THE ORIGINAL PROBLEM
To numerically solve the Earth-spacecraft-Moon three-body problem, a
MATLAB code (appendix A) implementing the fourth-order Runge-Kutta (RK4)
numerical integration scheme is developed to propagate the initial conditions forward in
time. [16] A trip time of 3.5 days is considered. An integration time step of 20 seconds is
used. Consistent with the RK4 algorithm, the lunar ephemeris are updated every half
time step. Choice of the reference frame for describing the spacecraft motion affects the
numerical error in the simulation. In order to avoid precision loss due to round-off error,
reference frame is switched from the EC to the MC at 2.905 days. [16] After the switch
point, equation of motion in the MC frame is integrated for d. This solution is then
converted to output a solution for r at each instant in time.
Figure 5 illustrates the resulting numerical solution of the original problem for the
spacecraft?s trajectory along with the Moon?s trajectory. Figures 6 and 7 show the plots
of X and Y with respect to time. Figures 8 and 9 show the plots of ? and ? with respect
to time. Figures 10 and 11 show the plots of ? and ? with respect to time. It can be seen
from the acceleration plots that the spacecraft experiences a higher acceleration when in
proximity to the Earth and the Moon. This is due to the strong nature of the gravitational
forces exerted by the Earth and the Moon on the spacecraft. The higher acceleration is
also reflected in the position and the velocity plots. The velocity of the spacecraft when
15
in proximity to the Earth and the Moon is higher, and it also changes rapidly. High
curvature of the trajectory during the initial and the final phases of the journey which
indicates a rapid change in position during these phases can also be seen by a close
examination of figures 5, 6 and 7. This behavior will affect the accuracy of the curve fits
calculated in the next chapter.
Figure 5. Spacecraft and Moon Trajectories
16
Figure 6. Time vs. X Position
Figure 7. Time vs. Y Position
17
Figure 8. Time vs. X Velocity
Figure 9. Time vs. Y Velocity
18
Figure 10. Time vs. X Acceleration
Figure 11. Time vs. Y Acceleration
19
CHAPTER 4
CURVE FITTING METHODS
Generating an accurate curve fit to the numerical solution to the original problem
is a critical step in MNP as subsequent analysis is based on this curve fit. The curve-fit
accuracy controls the nearness of the nearby problem. An analytical curve fit has to
satisfy two criteria in order to be considered accurate. The magnitude of the source terms
generated should be small, and the curve fit should maintain a particular order of
continuity in order to maintain slope continuity of the source terms. The order of
continuity to be maintained is problem dependent. The equations of motion of the three-
body dynamics are second order differential equations. Therefore, C3 continuity is
needed in the curve fit for this problem. Various curve fitting methods are available. The
methods used for generating curve fits in this study are least squares, cubic splines, and
fifth-degree Hermite splines. The accuracy of these curve fits is indicated by both the
residuals with respect to the numerical solution and the magnitude of the analytical
source terms calculated. The governing equations contain ? ?. However,
independent curve fits are only needed for and
Least Squares
The method of least squares is a standard approach to approximate the solution
of over determined systems, i.e., sets of equations in which there are more equations than
unknowns. "Least squares" means that the overall solution minimizes the sum of the
20
squares of the errors made in solving every single equation. Assume that we want to
approximate a function y(t), t being the independent variable. Assume that there are m
observations, i.e., values of y measured at specific values of t. [17]
),( ii tyy ? mi ,...,1? (20)
The idea is to model y(t) by a linear combination of n basis functions:
)(...)()( 11 txtxty nn?? ??? (21)
For example, a second degree polynomial can be written as 231201 txtxtx ?? , i.e., as
a linear combination of the basis functions t0, t1, and t2. The design matrix H is a
rectangular matrix of order m ? n with elements
)(, ijji th ?? (22)
The design matrix generally has more rows than columns. In matrix-vector notation:
Hxy? (23)
H is not invertible, but a pseudo inverse can be calculated as follows:
HxHyH TT ? (24)
yHHHx TT 1)( ?? (25)
In this study, the pseudo inverse is not calculated but instead the polyfit function in
MATLAB is used. The polyfit MATLAB file forms the Vandermonde matrix, H, whose
elements are powers of t:
jniji tH ??, (26)
21
It then uses the backslash operator to solve the least squares problem shown in Eq. (25).
The backslash operator selects from a variety of algorithms depending upon the structure
of the matrix H. Function ?polyfit (t, y, n)? finds the coefficients of a polynomial p(t) of
degree n that fits the data y best in a least squares sense. Various curve fits and their
residuals with respect to the numerical solution are shown in the following figures.
The residuals of these curve fits show that as the degree of the polynomial
increases, the curve fits are more accurate. However, as mentioned earlier, in MNP, in
order to be considered accurate a curve fit has to satisfy one more criterion, i.e., the
analytical source terms generated by these curve fits should be small in magnitude.
Source terms are discussed in detail in chapter 5.
22
Figure 12. Curve Fit to X ? Least Squares Using 3rd Degree Polynomial
Figure 13. Curve Fit Residuals for X - Least Squares Using 3rd Degree Polynomial
23
Figure 14. Curve Fit to X ? Least Squares Using 5th Degree Polynomial
Figure 15. Curve Fit Residuals for X - Least Squares Using 5th Degree Polynomial
24
Figure 16. Curve Fit to X ? Least Squares Using 20th Degree Polynomial
Figure 17. Curve Fit Residuals for X - Least Squares Using 20th Degree Polynomial
25
Figure 18. Curve Fit to Y ? Least Squares Using 3rd Degree Polynomial
Figure 19. Curve Fit Residuals for Y - Least Squares Using 3rd Degree Polynomial
26
Figure 20. Curve Fit to Y ? Least Squares Using 5th Degree Polynomial
Figure 21. Curve Fit Residuals for Y - Least Squares Using 5th Degree Polynomial
27
Figure 22. Curve Fit to Y ? Least Squares Using 20th Degree Polynomial
Figure 23. Curve Fit Residuals for Y - Least Squares Using 20th Degree Polynomial
28
Cubic Splines
Spline interpolation is a form of interpolation where piecewise functions called
?splines? are used. This has the advantage relative to a global least squares fit that
similar accuracy can be achieved using local fits that are each lower order.
S1 S2 S3 Sn
t1 t2 t3 tn tn+1
Figure 24. Schematic of Cubic Splines Interpolation
The basic idea of cubic splines is to fit a cubic polynomial on each interval
between points ti and ti+1 for i = 1,?, n. [18]
32 )()()()( iiiiiiii ttdttcttbatS ??????? (27)
This system has n+1 spline points (ti for i = 1,?, n+1) and n spline zones
(Si for i = 1,?, n) as shown in figure 24. The conditions that are used to construct these
polynomials are explained below for given , ?, ? data.
? ? iii XtS ? i = 1,?., n (28)
? ? 11 ?? ? nnn XtS
? ? ? ?iiii tStS 1?? i = 2,?., n (29)
? ? ? ?iiii tStS 1?? ?? i = 2,?., n (30)
? ? ? ?iiii tStS 1?? ???? i =2,?., n (31)
? ? 111 XtS ???? ? (32)
? ? 11 ?? ? nnn XtS ????
Conditions (28) set the value at each node point. Condition (29) makes the solution 0th
order continuous. Condition (30) makes the solution 1st order continuous, while
29
condition (31) makes the solution 2nd order continuous. Conditions (32) are the two
additional end point conditions.
Curve fits to X and Y using cubic splines are calculated using 8, 63, 504 and 5040
spline zones. The number 5040 corresponds to one third of the total time instances at
which the numerical solution is calculated. It was seen that further increase in the
number of spline zones did not improve the accuracy of the curve fit. A MATLAB
algorithm (appendix C) is developed to calculate the coefficients of these cubic splines.
[18] Various curve fits and their residuals with respect to the numerical solution are
shown in following figures. The residuals of these curve fits show that as the number of
spline zones increases, the curve fits are more accurate.
It is also seen that the residuals are higher in the initial and the final phases of the
trajectory. As discussed in chapter 2, the spacecraft experiences higher acceleration
during the initial and the final phases of the journey which results in the rapid change in
the position during these phases. The rapid change in the position results in high
curvature of the trajectory during these phases. Since the above cubic splines fit is
calculated using fixed resolution throughout the entire trajectory, it fails to capture these
high curvature regions precisely resulting in a less accurate curve fit in these regions.
The source terms calculated using cubic splines fit (discussed in chapter 5) indicate lesser
overall accuracy of cubic splines fit as compared to other methods discussed later in this
chapter. Therefore, no effort has been made to use a multi-resolution cubic splines fit in
order to capture the high curvature regions more precisely.
30
Figure 25. Curve Fit to X - Cubic Splines Using 8 Spline Zones
Figure 26. Curve Fit Residuals for X - Cubic Splines Using 8 Spline Zones
31
Figure 27. Curve Fit to X - Cubic Splines Using 63 Spline Zones
Figure 28. Curve Fit Residuals for X - Cubic Splines Using 63 Spline Zones
32
Figure 29. Curve Fit to X - Cubic Splines Using 504 Spline Zones
Figure 30. Curve Fit Residuals for X - Cubic Splines Using 504 Spline Zones
33
Figure 31. Curve Fit to X - Cubic Splines Using 5040 Spline Zones
Figure 32. Curve Fit Residuals for X - Cubic Splines Using 5040 Spline Zones
34
Figure 33. Curve Fit to Y - Cubic Splines Using 8 Spline Zones
Figure 34. Curve Fit Residuals for Y - Cubic Splines Using 8 Spline Zones
35
Figure 35. Curve Fit to Y - Cubic Splines Using 63 Spline Zones
Figure 36. Curve Fit Residuals for Y - Cubic Splines Using 63 Spline Zones
36
Figure 37. Curve Fit to Y - Cubic Splines Using 504 Spline Zones
Figure 38. Curve Fit Residuals for Y - Cubic Splines Using 504 Spline Zones
37
Figure 39. Curve Fit to Y - Cubic Splines Using 5040 Spline Zones
Figure 40. Curve Fit Residuals for Y - Cubic Splines Using 5040 Spline Zones
38
Fifth-degree Hermite Splines
An alternative approach to spline interpolation is to use ?fifth-degree Hermite
splines? instead of ?cubic splines?. The schematic of spline interpolation using fifth-
degree Hermite splines is as shown in figure 41.
S1 S2 S3 Sn
t1 t2 t3 tn tn+1
Figure 41. Schematic of Hermite Splines Interpolation
Here, the basic idea is to fit a fifth degree polynomial on each interval between
points ti and ti+1 for i = 1,?, n. [18]
5432 )()()()()()( iiiiiiiiiiii ttfttettdttcttbatS ??????????? (33)
This system has n+1 spline points (ti for i = 1,?, n+1) and n spline zones
(Si for i = 1,?, n). The conditions that are used to construct these polynomials are
explained below for given , ?, ? data.
? ? iii XtS ? i = 1,?., n (34)
? ? 11 ?? ? nnn XtS
? ? iii XtS ?? ? i = 1,?., n (35)
? ? 11 ?? ? nnn XtS ??
? ? ? ?iiii tStS 1?? i = 2,?., n (36)
? ? ? ?iiii tStS 1?? ?? i = 2,?., n (37)
? ? ? ?iiii tStS 1?? ???? i = 2,?., n (38)
? ? ? ?iiii tStS 1?? ?????? i = 2,?., n (39)
39
? ? 111 XtS ???? ? (40)
? ? 11 ?? ? nnn XtS ????
Conditions (34) set the value at each node point. Conditions (35) set the first derivative
at each node point. Typically, the cubic Hermite spline form consists of two control
points and two control tangents at the boundaries for each polynomial. Conditions (34)
and (35) make each polynomial of the spline fit to be in ?Hermite? form. Here, the
additional degrees of freedom in the fifth degree polynomial are used to enforce
additional continuity requirements. Condition (36) makes the solution 0th order
continuous. Condition (37) makes the solution 1st order continuous. Condition (38)
makes the solution 2nd order continuous, while condition (39) makes the solution 3rd order
continuous. Conditions (40) are the two additional end point conditions.
Curve fits to X and Y using fifth-degree Hermite splines are calculated using 8, 63,
504 and 5040 spline zones. A MATLAB algorithm (appendix D) is developed to
calculate the coefficients of these fifth-degree Hermite splines. [18] Various curve fits
and their residuals with respect to the numerical solution are shown in figures 42 to 57.
The residuals of these curve fits show that as the number of spline zones increases, the
curve fits are more accurate.
The phenomenon of higher residuals in the initial and the final phases of the
trajectory is also seen in the case of fifth-degree Hermite splines curve fit. An effort has
been made to address this phenomenon by using an iterative multi-resolution fifth-degree
Hermite splines fit in order to capture the high curvature regions more precisely. A
MATLAB code (appendix E) is developed for this purpose. In this code, the trajectory is
divided into four regions of equal amounts of time with region 1 and 4 covering the initial
40
and the final phases of the trajectory respectively. To start with, each region is divided
into two spline zones (total of eight spline zones). A curve fit is calculated using these
spline zones and each region is checked for the maximum value of residual which should
be below a selected threshold. If the maximum value of residual in a particular region is
below a selected threshold, addition of spline points to that particular region is
terminated. One additional spline point is added to each region where the threshold is not
met. For each iteration, the spline points within each region are evenly distributed.
However, spline points falling in the middle of a time step are rounded to the nearest
integer number of time steps. In this manner, the maximum value of residual in each
region is made to be somewhat consistent indicating that the curve fit exhibits nearly the
same level of accuracy throughout the trajectory.
It is seen that, if the selected threshold is below 10-3 meters, the maximum value
of residual fails to converge. Therefore, the threshold was selected to be 10-3 meters. A
reason for this is that the number of spline zones in certain regions began to approach the
number of data points. Tables 2 and 3 list, for the curve fits to X and Y, the number of
spline zones required in each region so that the maximum value of residual in that region
is below the selected threshold. It can be seen that the number of spline zones required
by regions 1 and 4 is much higher than that required by regions 2 and 3.
Region No. of Spline Zones
1 1708
2 24
3 17
4 933
Total 2682
Table 2. Region Wise Number of Spline Zones Using Multi-resolution Fifth-degree Hermite Splines
Curve Fit to X
41
Region No. of Spline Zones
1 1779
2 20
3 15
4 1043
Total 2857
Table 3. Region Wise Number of Spline Zones Using Multi-resolution Fifth-degree Hermite Splines
Curve Fit to Y
Figures 58 and 59 show the plot of the residuals with respect to time for X and Y
using the above discussed multi-resolution curve fit. It can be seen that, though this
curve fit requires less number of total spline zones and also maintains the level of
accuracy throughout the trajectory, the overall accuracy of this curve fit is fairly similar
to that of the fixed-resolution curve fit discussed earlier.
42
Figure 42. Curve Fit to X ? Fifth-degree Hermite Splines Using 8 Spline Zones
Figure 43. Curve Fit Residuals for X - Fifth-degree Hermite Splines Using 8 Spline Zones
43
Figure 44. Curve Fit to X ? Fifth-degree Hermite Splines Using 63 Spline Zones
Figure 45. Curve Fit Residuals for X - Fifth-degree Hermite Splines Using 63 Spline Zones
44
Figure 46. Curve Fit to X ? Fifth-degree Hermite Splines Using 540 Spline Zones
Figure 47. Curve Fit Residuals for X - Fifth-degree Hermite Splines Using 540 Spline Zones
45
Figure 48. Curve Fit to X ? Fifth-degree Hermite Splines Using 5040 Spline Zones
Figure 49. Curve Fit Residuals for X - Fifth-degree Hermite Splines Using 5040 Spline Zones
46
Figure 50. Curve Fit to Y ? Fifth-degree Hermite Splines Using 8 Spline Zones
Figure 51. Curve Fit Residuals for Y - Fifth-degree Hermite Splines Using 8 Spline Zones
47
Figure 52. Curve Fit to Y ? Fifth-degree Hermite Splines Using 63 Spline Zones
Figure 53. Curve Fit Residuals for Y - Fifth-degree Hermite Splines Using 63 Spline Zones
48
Figure 54. Curve Fit to Y ? Fifth-degree Hermite Splines Using 504 Spline Zones
Figure 55. Curve Fit Residuals for Y - Fifth-degree Hermite Splines Using 504 Spline Zones
49
Figure 56. Curve Fit to Y ? Fifth-degree Hermite Splines Using 5040 Spline Zones
Figure 57. Curve Fit Residuals for Y - Fifth-degree Hermite Splines Using 5040 Spline Zone
50
Figure 58. Curve Fit Residuals for X ? Multi-resolution Fifth-degree Hermite Splines Fit
Figure 59. Curve Fit Residuals for Y - Multi-resolution Fifth-degree Hermite Splines Fit
51
CHAPTER 5
GENERATING ANALYTICAL SOURCE TERMS
As mentioned earlier, in MNP, for a curve fit to be considered accurate the
analytical source terms generated by it should be small in magnitude. The smaller the
source terms, the nearer the nearby problem is to the original problem. Source terms are
calculated for all the curve fits discussed in the previous chapter. Consider rewriting Eq.
(14) in an operator form.
? ?? ? 0?drrrL ????????? ???? 33322dd, ??? drttt me (41)
The curve fits for X and Y generated in the previous chapter can be assembled into a
curve fit for the position vector, labeled ??tr Operating the original problem on this
curve fit defines the source terms (appendix H), ??ts .
? ?? ? ? ?ttt srL ?, (42)
By construction, the curve fit is the exact solution of a modified equation, which is the
nearby problem.
? ?? ? ? ? 0srL ?? ttt , (43)
Therefore, the source terms are equivalent to a set of time-dependent perturbing
accelerations that result in a trajectory following the curve fit. In Eq. (43), as ??ts
approaches zero, the nearby problem approaches the original problem.
52
We will now see the magnitude of the analytical source terms generated using
various curve fits discussed in the previous chapter.
Analytical Source Terms Using Least Squares
Following figures show the source-term histories for both X and Y for the least-
squares curve fits using varying degrees of polynomials. It can be seen that the
magnitude of the source terms increases as the degree of the polynomial increases. Note
that, though the curve fit using a higher-degree polynomial is more accurate than that
using a lower-degree polynomial; the source terms generated by using the higher-degree
polynomial curve fit are larger in magnitude than those generated by using the lower-
degree polynomial curve fit. This can be attributed to the fact that the least-squares
solution is trying to minimize errors only in position. The derivatives of the position do
not play any role in calculating the curve fit. The larger magnitude of the source terms
indicates that the least squares is not a feasible option for the construction of the nearby
problem.
53
Figure 60. Source Term Histories for X - Least Squares Using 3rd Degree Polynomial
Figure 61. Source Term Histories for X - Least Squares Using 5th Degree Polynomial
54
Figure 62. Source Term Histories for X - Least Squares Using 20th Degree Polynomial
Figure 63. Source Term Histories for Y - Least Squares Using 3rd Degree Polynomial
55
Figure 64. Source Term Histories for Y - Least Squares Using 5th Degree Polynomial
Figure 65. Source Term Histories for Y - Least Squares Using 20th Degree Polynomial
56
Analytical Source Terms Using Cubic Splines
Following figures show the source-term histories for both X and Y for the cubic
splines curve fits using varying numbers of spline zones. It can be seen that the
magnitude of the source terms decreases as the number of spline zones increases. Using
moderate number of spline zones, the magnitude of the source terms generated is small.
However, close examination of following figures shows that the source terms exhibit
slope discontinuities at the spline points. This is because cubic splines are only C2
continuous. As mentioned earlier, in MNP, the curve fit should maintain a particular
order of continuity in order to maintain the slope continuity of the source terms. The
equations of motion of three-body dynamics are second order differential equations and
demand C3 continuity in the curve fit. Since the continuity criterion is not satisfied, cubic
splines cannot be used to construct the nearby problem.
57
Figure 66. Source Term Histories for X - Cubic Splines Using 8 Spline Zones
Figure 67. Source Term Histories for X - Cubic Splines Using 63 Spline Zones
58
Figure 68. Source Term Histories for X - Cubic Splines Using 504 Spline Zones
Figure 69. Source Term Histories for X - Cubic Splines Using 5040 Spline Zones
59
Figure 70. Source Term Histories for Y - Cubic Splines Using 8 Spline Zones
Figure 71. Source Term Histories for Y - Cubic Splines Using 63 Spline Zones
60
Figure 72. Source Term Histories for Y - Cubic Splines Using 504 Spline Zones
Figure 73. Source Term Histories for Y - Cubic Splines Using 5040 Spline Zones
61
Analytical Source Terms Using Fifth-degree Hermite Splines
Following figures show the source-term histories for both X and Y for the fifth-
degree Hermite splines curve fits using varying numbers of spline zones. It can be seen
that the magnitude of the source terms decreases as the number of spline zones increases.
Using moderate number of spline zones, the magnitude of the source terms generated is
extremely small. Since the fifth-degree Hermite splines are C3 continuous, the source
terms are found to be slope continuous; thus satisfying the continuity criterion required
by MNP. Extremely small magnitude of source terms indicates that the nearby problem
constructed using these source terms can be considered to be a good representation of our
original problem. Therefore, the fifth-degree Hermite splines are a feasible option for
constructing the nearby problem. Because small source terms were achieved using a
reasonable number of fixed-resolution spline zones, for simplicity this approach was used
instead of the multi-resolution approach. Construction of the nearby problem using these
fifth-degree Hermite splines is discussed in chapter 6.
62
Figure 74. Source Term Histories for X ? Fifth-degree Hermite Splines Using 8 Spline Zones
Figure 75. Source Term Histories for X ? Fifth-degree Hermite Splines Using 63 Spline Zones
63
Figure 76. Source Term Histories for X ? Fifth-degree Hermite Splines Using 504 Spline Zones
Figure 77. Source Term Histories for X ? Fifth-degree Hermite Splines Using 5040 Spline Zones
64
Figure 78. Source Term Histories for Y ? Fifth-degree Hermite Splines Using 8 Spline Zones
Figure 79. Source Term Histories for Y ? Fifth-degree Hermite Splines Using 63 Spline Zones
65
Figure 80. Source Term Histories for Y ? Fifth-degree Hermite Splines Using 504 Spline Zones
Figure 81. Source Term Histories for Y ? Fifth-degree Hermite Splines Using 5040 Spline Zones
66
CHAPTER 6
NEARBY PROBLEM TO THE ORIGINAL PROBLEM
A nearby problem to our original problem of Earth-spacecraft-Moon three-body
dynamics is constructed by adding the source terms obtained from the curve fit using the
Fifth-degree Hermite splines to the governing equations of the original problem as shown
in Eq. (43). As discussed earlier, the curve fit serves as an exact solution to this nearby
problem. In order to be able to calculate the error in the numerical solution to the nearby
problem, it is solved numerically using the same numerical scheme that was used for
solving the original problem (appendix H). Note that solving the nearby problem using
RK4 with a time step of 20 seconds requires evaluating the source terms every 10
seconds.
Following figures show the construction of various nearby problems using Fifth-
degree Hermite splines with varying numbers of spline zones. It is seen that as the
number of spline zones increases, both the exact and the numerical solution to the nearby
problem start nearing the solution of the original problem.
The nearness of the nearby problem constructed using the Fifth-degree Hermite
splines with moderate number of spline zones has already been established through
previous chapters. Therefore, this nearby problem can be considered to be a good
representation of the original problem of interest. Thus, the exact numerical error in the
nearby problem can be used as an estimate of the numerical error in the original problem.
67
Figure 82. Nearby Problem to X - Using Fifth-degree Hermite Splines with 8 Spline Zones
Figure 83. Nearby Problem to X - Using Fifth-degree Hermite Splines with 63 Spline Zones
68
Figure 84. Nearby Problem to X - Using Fifth-degree Hermite Splines with 504 Spline Zones
Figure 85. Nearby Problem to X - Using Fifth-degree Hermite Splines with 5040 Spline Zones
69
Figure 86. Nearby Problem to Y - Using Fifth-degree Hermite Splines with 8 Spline Zones
Figure 87. Nearby Problem to Y - Using Fifth-degree Hermite Splines with 63 Spline Zones
70
Figure 88. Nearby Problem to Y - Using Fifth-degree Hermite Splines with 504 Spline Zones
Figure 89. Nearby Problem to Y - Using Fifth-degree Hermite Splines with 5040 Spline Zones
71
Figures 90 and 91 show, for X and Y, the exact error in numerical solutions to the
various nearby problems discussed above. It is seen that using varying numbers of spline
zones the exact solution to the nearby problem differs, but the error in the numerical
solution to the nearby problem remains fairly consistent. In general, the error tends to
grow as the numerical solution is propagated forward in time. In particular, the error
tends to change quickly as the spacecraft reaches the Moon in the terminal portion of the
trajectory. From the combined results for X and Y, the accumulated error is
approximately of the order of 10 meters.
Figure 90. Exact Error in Nearby Problem to X Using Fifth-degree Hermite Splines
72
Figure 91. Exact Error in Nearby Problem to Y Using Fifth-degree Hermite Splines
It is imperative to examine the reliability of these error estimates. This could be
done by comparing these estimates to the estimates given by some other methods.
Checking for the reliability of the numerical scheme used is also a good idea.
One of the best ways to verify the numerical scheme is to calculate the observed
order of accuracy and see how well it matches the formal order of accuracy. The RK4
numerical integration scheme is a 4th order accurate method i.e. the formal order of
accuracy of RK4 method is four. The observed order of accuracy, p, is computed for the
various numerical solutions of the nearby problem (constructed using fifth-degree
Hermite splines with 5040 spline zones) using different meshes. A mesh is indicated by
the step size used. Thus, the fine mesh uses a smaller step size than the coarse mesh.
Since the exact solution is known, the observed order of accuracy can be computed by the
following relation [19]:
73
( )
(44)
In Eq. (44), r is the grid refinement factor (the ratio between the coarse mesh and the fine
mesh). One of the tests used by Berry and Healy [3] to verify the accuracy of the
numerical integrators is the step-size halving test. For the step-size halving test, the
reference integration is produced with the same integrator but with the step-size cut in
half. On this basis, here, a value of r=2 is used. E2 is the L2 norm of the errors between
the exact solution and the numerical solution calculated at each instant in time using the
coarse mesh while E1 is the L2 norm of the errors between the exact solution and the
numerical solution calculated at each instant in time using the fine mesh. Figure 92
shows the observed order of accuracies for various numerical solutions to the nearby
problem (constructed using fifth-degree Hermite splines with 5040 spline zones) using
various step sizes. It can be seen that as the step size is reduced the observed order of
accuracy approaches four (which is the formal order of accuracy). This suggests that the
numerical solution is reliable.
Figure 92. Observed Order of Accuracies for Various Numerical Solutions to the Nearby
Problem using Different Step Sizes
74
Richardson extrapolation (RDE) can also be used to get an estimate of the exact
solution. [20, 21] RDE involves computation of numerical solutions on two or more
meshes. Solutions on these different meshes are then used to compute a higher-order
estimate of the exact solution. This estimate of the exact solution can then be used to
estimate the numerical error. For a pth order accurate scheme with solutions on a fine
mesh ( ) and a coarse mesh ( ) with refinement factor r, can be approximated
as: [19]
can be calculated using both global and formal order of accuracies. The error on
the fine mesh can then be given as:
Figures 93 and 94 show the comparison of error estimates using MNP, RDE with
formal order of accuracy, and RDE with observed order of accuracy. For MNP, a time
step of 20 seconds is used and therefore for RDE the fine mesh uses a time step of 20
seconds while the coarse mesh uses a time step of 40 seconds. It can be seen that the
error estimates using MNP and RDE match in the order of magnitude.
75
Figure 93. Numerical Error Estimates for X Using Various Methods
Figure 94. Numerical Error Estimates for Y Using Various Methods
76
CHAPTER 7
CONCLUSION AND RECOMMENDATIONS
In orbital mechanics, analytical solutions to many systems do not exist and thus,
an accurate numerical solution is imperative to mission planning. One such system is the
Earth-spacecraft-Moon three-body dynamics or the translunar trajectories. The
translunar trajectories are chaotic in nature. They are extremely sensitive when the
spacecraft is in proximity to the Earth and the Moon. So the numerical solution needs to
be as accurate as possible otherwise a small error can get escalated later on. MNP is
employed to validate the accuracy of this numerical solution by estimating the numerical
error in it. A nearby problem (having an exact solution to it) to the Earth-spacecraft-
Moon three-body problem is constructed. Fifth-degree Hermite splines are found to be a
feasible option to construct this nearby problem as it satisfies all the conditions necessary
to demonstrate its nearness to the original problem. This allows the exact error in the
numerical solution to this nearby problem to be considered as a good estimation of the
error in the numerical solution to the original problem.
The MNP estimate of the magnitude of the numerical error in the simulation of
translunar trajectories is of the order of 10 meters. This accuracy may be sufficient for
many aspects of mission planning; however, for critical mission phases higher accuracy
may be desired. Also, this error is accumulated when the propagation time is 3.5 days.
The propagation time for other missions could range from a few days to even years.
77
During such long propagation times, substantial amount of error can be accumulated.
This could lead to inaccurate results. Of course, accurate simulation also depends on
accurate modeling of the system dynamics. Mission planning requires both the
development of accurate models and the accurate numerical solution of those models.
This thesis demonstrates the usefulness of MNP in providing reliable estimates of
the error in the numerical solutions to the problems in orbital mechanics. These error
estimates, however, are dependent on the numerical scheme used and the type of problem
studied. The reliability of MNP can further be verified by using various numerical
integration schemes and/or by studying problems involving forces of a more complex
nature.
78
REFERENCES
1. Huang, T.-Y. and Innanen, K. A., ?The accuracy check in numerical integration of
dynamical systems,? Astronomical Journal, Vol. 88, No. 6, 1983, pp. 870?876, June
1983.
2. Fox, K., ?Numerical integration of the equations of motion of celestial mechanics,?
Celestial Mechanics, Vol. 33, No. 2, 1984, pp. 127-142, June 1984.
3. M. Berry, L. Healy, ?Comparison of accuracy assessment techniques for numerical
integration,? 13th AAS/AIAA Space Flight Mechanics Meeting, Ponce, Puerto Rico,
9-13 February 2003.
4. Montenbruck, O., ?Numerical Integration Methods for Orbital Motion,? Celestial
Mechanics and Dynamical Astronomy, Vol. 53, pp. 59?69, 1992.
5. Hadjifotinou, K. G. and Gousidou-Koutita, M., ?Comparison of Numerical Methods
for the Integration of Natural Satellite Systems,? Celestial Mechanics and Dynamical
Astronomy, Vol. 70, No. 2, 1998, pp. 99?113, 1998.
6. Roache, P.J., ?Code Verification by the Method of Manufactured Solutions,? ASME
Journal of Fluids Engineering, Vol. 124, March 2002, pp. 4-10.
7. Zadunaisky, P. E., ?A Method for the Estimation of Errors Propagated in the
Numerical Solution of a System of Ordinary Di?erential Equations,? In Contopoulos,
79
G., editor, The Theory of Orbits in the Solar System and in Stellar Systems, pp. 281?
287, New York, 1966. International Astronomical Union, Academic Press.
8. Zadunaisky, P. E., ?On the Estimation of Errors Propagated in the Numerical
Integration of Ordinary Differential Equations,? Numerische Mathematik, Vol. 27,
No. 1, 1976, pp. 21?39, 1976.
9. Zadunaisky, P. E., ?On the Accuracy in the Numerical Solution of the N-Body
Problem,? Celestial Mechanics, Vol. 20, pp. 209?230, 1979.
10. Zadunaisky, P. E., ?On the Accuracy in the Numerical Computation of Orbits,? In
Giacaglia, G.E. O., editor, Periodic Orbits, Stability and Resonances, pp. 216?227,
Dordrecht, Holland, 1970.D. Reidel Publishing Company.
11. Junkins, J.L., and Lee. S., ?Validation of Finite-Dimensional Approximate Solutions
for Dynamics of Distributed-Parameter Systems,? Journal of Guidance, Control, and
Dynamics, Vol. 18, No. 1, 1995, pp. 87-95.
12. M.M. Hopkins and C.J. Roy, "Introducing the Method of Nearby Problems."
European Congress on Computational Methods in Applied Sciences and Engineering,
ECCOMAS 2004, P. Neittaanmaki, T. Rossi, S. Korotov, E. Onate, J. Periaux, and
D. Knorzer (eds.), July 24-28, 2004.
13. Roy, C. J., and Hopkins, M. M., ?Discretization Error Estimates using Exact
Solutions to Nearby Problems,? AIAA Paper 2003-0629, January 2003
14. Bate R., Mueller D., White J., Fundamentals of Astrodynamics, Dover Publications,
New York, 1971.
15. V. Chobotov, "Orbital Mechanics," AIAA Education Series, Virginia 2002.
80
16. M.P. Vautier and A.J. Sinclair, ?Coordinate Switching for Accurate Simulation of
Chaotic Trajectories?, AAS 08-271, F. Landis Markley Astronautics Symposium,
Cambridge, Maryland, 29 June ? 2 July 2008.
17. C.B. Moler, "Numerical Computing with MATLAB." SIAM, Philadelphia 2004.
18. G.E. Mulleges and F. Uhlig, "Numerical Algorithms with Fortran," Springer-Verlag,
Berlin 1996
19. Roy, C. J., ?Review of Code and Solution Verification Procedures for Computational
Simulation? Journal of Computational Physics, Vol. 205, 2005, pp. 131-156
20. L.F. Richardson, The approximate arithmetical solution by finite differences of
physical problems involving differential equations with an application to the stresses
in a masonry dam, Trans. Royal Society London, Ser. A 210 (1910) 307?357.
21. L.F. Richardson, The deferred approach to the limit, Trans. Royal Society London,
Ser. A 226 (1927) 229?361.
81
APPENDICES
A. MATLAB CODE TO NUMERICALLY SOLVE THE ORIGINAL PROBLEM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% NUMRICAL SOLUTION TO THE ORIGINAL PROBLEM %%%%%%%%%
%%%%%%%%%%% PROPOGATION IN THE EC->MC FRAME %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
close all;
clc;
format long;
format compact;
%%%%%%%%%%%%%%%%%%%%%%%%% SETTING UP THE PROBLEM %%%%%%%%%%%%%%%%%%%%%%
%%%%% USER INPUTS %%%%%
alt = 359750; %meters
angle = -36.890; %degrees
dt = 20; %secs
fprintf('alt = %6.0f, ang = %6.3f, dt = %6.3f \n\n', alt, angle, dt);
deltav = 3102.13; %m/s
t0 = 0; %secs
tf = 3.5*86400; %302400 secs (1 day = 86400 secs)
SP = 2.905*86400; % Switch Point at 250992 secs
folder = char(['alt' int2str(alt) 'ang' int2str(abs(angle*1000)) 'dt'
int2str(dt)]);
fid = fopen([folder 'State_SP' '.txt'], 'w');
f = fopen([folder 'Moon_State_SP' '.txt'], 'w');
%%%%% CONSTANTS %%%%%
GMm = 4.90266e12; %m3/s2
GMe = 3.98600436e14; %m3/s2
radE = 6372797; %meters
EMdist = 384400000; %meters
omega = sqrt((GMe+GMm)/EMdist^3); %rad/sec
%Tm = 2*pi/omega = 27d 6h 49m 50.34879957310977s
alpha = angle*pi/180; %radians
h = radE + alt; %meters
Vorbit = sqrt(GMe/h); %m/s
V = Vorbit + deltav; %m/s
%%%%% INITIAL CONDITIONS %%%%%
%%%%% For the Spacecraft %%%%%
rxi = h*sin(alpha); %meters
82
ryi = -h*cos(alpha); %meters
vxi = V*cos(alpha); %m/s
vyi = V*sin(alpha); %m/s
x = [rxi;ryi;vxi;vyi];
fprintf(fid, '% 6.0f, % .16e, % .16e, % .16e, % .16e nn', t0, x);
%%%%% For the Moon %%%%%
xmt = EMdist*cos(omega*t0);
ymt = EMdist*sin(omega*t0);
fprintf(f, '% 6.0f, % .16f, % .16f nn', t0, xmt, ymt);
%%%%%%%%%%%%%%%%%%%%%%%%% PERFORM THE PROPOGATION %%%%%%%%%%%%%%%%%%%%%
%%%%% EARTH CENTERED FRAME %%%%%
i = 1;
for t = t0:dt:SP-dt
%Use the moon position at the end of the previous interval as the
%position at the start of this interval.
xm = xmt;
ym = ymt;
%Calculate the moon position at the middle and end of this
%interval.
xmh = EMdist*cos(omega*(t + 0.5*dt));
ymh = EMdist*sin(omega*(t + 0.5*dt));
xmt = EMdist*cos(omega*(t + dt));
ymt = EMdist*sin(omega*(t + dt));
t1 = t+dt;
fprintf(f, '% 6.0f, % .16f, % .16f nn', t1, xmt, ymt);
k1 = RK4_EC(x, xm, ym, GMm, GMe);
k2 = RK4_EC(x+(dt/2)*k1', xmh, ymh, GMm, GMe);
k3 = RK4_EC(x+(dt/2)*k2', xmh, ymh, GMm, GMe);
k4 = RK4_EC(x+(dt)*k3', xmt, ymt, GMm, GMe);
x = x + dt/6*(k1'+2*k2'+2*k3'+k4');
fprintf(fid, '% 6.0f, % .16e, % .16e, % .16e, % .16enn', t1, x);
i = i+1;
end
%%%%% SWITCH POINT %%%%%
%Convert state vector from EC to MC
rho = [xmt; ymt];
xmt_dot = -EMdist*omega*sin(omega*(t+dt));
ymt_dot = EMdist*omega*cos(omega*(t+dt));
rho_dot = [xmt_dot; ymt_dot];
d_state = x - [rho; rho_dot];
%%%%% MOON CENTERED FRAME %%%%%
j = 1;
for t = SP:dt:tf-dt
%Use the moon position at the end of the previous interval as the
%position at the start of this interval.
xm = xmt;
ym = ymt;
%Calculate the moon position at the middle and end of this
%interval.
xmh = EMdist*cos(omega*(t + 0.5*dt));
ymh = EMdist*sin(omega*(t + 0.5*dt));
xmt = EMdist*cos(omega*(t + dt));
83
ymt = EMdist*sin(omega*(t + dt));
t2 = t+dt;
fprintf(f, '% 6.0f, % .16f, % .16f nn', t2, xmt, ymt);
%Propogating the states forward in time
k1 = RK4_MC(d_state, xm, ym, GMm, GMe);
k2 = RK4_MC(d_state+(dt/2)*k1', xmh, ymh, GMm, GMe);
k3 = RK4_MC(d_state+(dt/2)*k2', xmh, ymh, GMm, GMe);
k4 = RK4_MC(d_state+(dt)*k3', xmt, ymt, GMm, GMe);
d_state = d_state + dt/6*(k1'+2*k2'+2*k3'+k4');
%Compute the moon state at the end of interval
rho = [xmt; ymt];
rho_dot = EMdist*omega*[-sin(omega*(t + dt)); cos(omega*(t + dt))];
%Convert state vector from MC to EC
x = d_state + [rho; rho_dot];
fprintf(fid, '% 6.0f, % .16e, % .16e, % .16e, % .16enn',t2,x);
j = j+1;
end
fclose(fid);
fclose(f);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
84
B. SUBROUTINES USED IN THE CODE GIVEN IN APPENDIX A
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% SUBROUTINE TO INTEGRATE %%%%%%%%%%%
%%%%%%%%%% THE EQUATION OF MOTION IN EC FRAME %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x_dot = RK4_EC(x,xmoon,ymoon,GMm,GMe)
r = x(1:2);
v = x(3:4);
rho = [xmoon; ymoon];
d = r-rho;
x_dot(1:2) = v;
x_dot(3:4) = -GMe*r/norm(r)^3-
GMm*(d/norm(d)^3+rho/norm(rho)^3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% SUBROUTINE TO INTEGRATE %%%%%%%%%%%
%%%%%%%%%% THE EQUATION OF MOTION IN MC FRAME %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d_dot = RK4_MC(d_state,xmoon,ymoon,GMm,GMe)
d = d_state(1:2);
v = d_state(3:4);
rho = [xmoon; ymoon];
r = d+rho;
d_dot(1:2) = v;
d_dot(3:4) = -GMm*d/norm(d)^3-GMe*(r/norm(r)^3-
rho/norm(rho)^3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
85
C. MATLAB CODE TO CALCULATE THE COEFFICIENTS OF
CUBIC SPLINES CURVE FIT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% MATLAB CODE TO CALCULATE THE COEFFICIENTS %%%%%%%%%%%
%%%%%%%%%%%%%%% OF CUBIC SPLINES CURVE FIT %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [X_fit Xd_fit Xdd_fit] = fCurveFit_Cu(time10,time20,X,Xdd)
% time10 => Time matrix with time step of 10 sec
% time20 => Time matrix with time step of 20 sec
% X => X position obtained from the simulation given in
Appendix A
% Xdd => X acceleration obtained from the simulation given in
Appendix A
S = size(time20); m = S(2); j = 1;
%Calculating Coefficients a
for i = 1:3:m
x(j) = time20(i);
y(j) = X(i);
a(j) = y(j);
j = j+1;
end
S = size(x); n = S(2);
%Calculating Coefficient c
c = zeros(1,n); A = zeros(n-2);
alpha = Xdd(1); beta = Xdd(m);
c(1) = 0.5*alpha; c(n) = 0.5*beta;
for i = 1:n-2
h1 = x(i+1)-x(i);
h2 = x(i+2)-x(i+1);
A(i,i) = 2*(h1+h2);
end
for i = 2:n-2
h2 = x(i+1)-x(i);
A(i-1,i) = h2;
A(i,i-1) = h2;
end
for i = 1:n-2
h2 = x(i+2)-x(i+1); h1 = x(i+1)-x(i);
a3 = a(i+2); a2 = a(i+1); a1 = a(i);
if i == 1
g(i) = (3/h2*(a3-a2))-(3/h1*(a2-a1))-h1*(alpha/2);
elseif i == n-2
g(i) = (3/h2*(a3-a2))-(3/h1*(a2-a1))-h1*(beta/2);
else
g(i) = (3/h2*(a3-a2))-(3/h1*(a2-a1));
end
end
C = fTRIDIAG(Amd,Ad,g,n-2);
for i = 2:n-1
c(i) = C(i-1);
86
end
%Calculating Coefficients b and d
for i = 1:n-1
h1 = x(i+1)-x(i);
a2 = a(i+1); a1 = a(i);
c2 = c(i+1); c1 = c(i);
b(i) = (1/h1*(a2-a1))-(h1/3*(c2+2*c1));
d(i) = (1/(3*h1))*(c2-c1);
end
%Evaluating Spline Fit
S = size(time10); m = S(2);
for i = 1:m
for j = 1:n-1
if time10(i) >= x(j) && time10(i) <= x(j+1)
X_fit(i) = d(j)*(time10(i)-x(j))^3 +
c(j)*(time10(i)-x(j))^2 +
b(j)*(time10(i)-x(j)) + a(j);
Xd_fit(i) = 3*d(j)*(time10(i)-x(j))^2 +
2*c(j)*(time10(i)-x(j)) + b(j);
Xdd_fit(i) = 6*d(j)*(time10(i)-x(j)) + 2*c(j);
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
87
D. MATLAB CODE TO CALCULATE THE COEFFICIENTS OF
FIFTH-DEGREE HERMITE SPLINES CURVE FIT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% MATLAB CODE TO CALCULATE THE COEFFICIENTS %%%%%%%%%%
%%%%%%%%%% OF FIFTH-DEGREE HERMITE SPLINES CURVE FIT %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [X_fit Xd_fit Xdd_fit] = fCurveFit_HS(time10,time20,X,Xd,Xdd)
% time10 => Time matrix with time step of 10 sec
% time20 => Time matrix with time step of 20 sec
% X => X position obtained from the simulation given in
Appendix A
% Xd => X velocity obtained from the simulation given in
Appendix A
% Xdd => X acceleration obtained from the simulation given in
Appendix A
S = size(time20); m = S(2); j = 1;
%Calculating Coefficients a and b
for i = 1:3:m
x(j) = time20(i);
y(j) = X(i);
yPR(j) = Xd(i);
a(j) = y(j);
b(j) = yPR(j);
j = j+1;
end
S = size(x); n = S(2);
%Calculating Coefficient c
c = zeros(1,n); A = zeros(n-2);
c(1) = 0.5*Xdd(1); c(n) = 0.5*Xdd(m);
alpha = 1; Beta1 = Xdd(1)/(2*(x(2)-x(1)));
Beta2 = Xdd(m)/(2*(x(n)-x(n-1)));
for i = 1:n-2
h1 = x(i+1)-x(i); h2 = x(i+2)-x(i+1);
if i == 1
Amd(i) = 3*(alpha/h1+1/h2);
elseif i == n-2
Amd(i) = 3*(1/h1+alpha/h2);
else
Amd(i) = 3*(1/h1+1/h2);
end
A(i,i) = Amd(i);
end
for i = 2:n-2
h2 = x(i+1)-x(i);
Ad(i-1) = -1/h2;
A(i-1,i) = Ad(i-1);
A(i,i-1) = Ad(i-1);
end
for i = 1:n-2
a3 = a(i+2); a2 = a(i+1); a1 = a(i);
88
b3 = b(i+2); b2 = a(i+1); b1 = b(i);
h2 = x(i+2)-x(i+1); h1 = x(i+1)-x(i);
if i == 1;
g(i) = 10*((a3-a2)/h2^3-(a2-a1)/h1^3)+4*(b1/h1^2-
1.5*b2*(1/h2^2-1/h1^2)-b3/h2^2)+Beta1;
elseif i == n-2
g(i) = 10*((a3-a2)/h2^3-(a2-a1)/h1^3)+4*(b1/h1^2-
1.5*b2*(1/h2^2-1/h1^2)-b3/h2^2)+Beta2;
else
g(i) = 10*((a3-a2)/h2^3-(a2-a1)/h1^3)+4*(b1/h1^2-
1.5*b2*(1/h2^2-1/h1^2)-b3/h2^2);
end
end
C = fTRIDIAG(Amd,Ad,g,n-2);
for i = 2:n-1
c(i) = C(i-1);
end
%Calculating Coefficients d
for i = 1:n
if i < n
h1 = x(i+1)-x(i);
a2 = a(i+1); a1 = a(i);
b2 = b(i+1); b1 = b(i);
c2 = c(i+1); c1 = c(i);
d(i) = 10/h1^3*(a2-a1)-2/h1^2*(2*b2+3*b1)+1/h1*(c2-
3*c1);
else
d4 = d(i-1); h4 = x(n)-x(n-1); b5 = b(i); b4 = b(i-1);
c5 = c(i); c4 = c(i-1);
d(i) = d4-2/h4^2*(b5-b4)+2/h4*(c5+c4);
end
end
%Calculating Coefficients e and f
for i = 1:n-1
h1 = x(i+1)-x(i); b2 = b(i+1); b1 = b(i);
c2 = c(i+1); c1 = c(i); d2 = d(i+1); d1 = d(i);
e(i) = 0.5/h1^3*(b2-b1)-1/h1^2*c1-0.25/h1*(d2+5*d1);
e1 = e(i);
f(i) = 0.1/(h1^3)*(c2-c1-3*d1*h1-6*e1*h1^2);
end
%Evaluating Spline Fit
S = size(time10); m = S(2);
for i = 1:m
for j = 1:n-1
if time10(i) >= x(j) && time10(i) <= x(j+1)
X_fit(i) = f(j)*(time10(i)-x(j))^5 +
e(j)*(time10(i)-x(j))^4 +
d(j)*(time10(i)-x(j))^3 +
c(j)*(time10(i)-x(j))^2 +
b(j)*(time10(i)-x(j)) + a(j);
Xd_fit(i) = 5*f(j)*(time10(i)-x(j))^4 +
4*e(j)*(time10(i)-x(j))^3 +
3*d(j)*(time10(i)-x(j))^2 +
2*c(j)*(time10(i)-x(j)) + b(j);
89
Xdd_fit(i) = 20*f(j)*(time10(i)-x(j))^3 +
12*e(j)*(time10(i)-x(j))^2 +
6*d(j)*(time10(i)-x(j)) + 2*c(j);
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
90
E. MATLAB CODE TO CALCULATE THE COEFFICIENTS OF
MULTI-RESOLUTION
FIFTH-DEGREE HERMITE SPLINES CURVE FIT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% MATLAB CODE TO CALCULATE THE COEFFICIENTS %%%%%%
%%%%%%%%%%%% OF FIFTH-DEGREE HERMITE SPLINES CURVE FIT %%%%%%
%%%%%%%%%%%% WITH ADAPTIVE ALGORITHM %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;
clear all;
close all;
format long E;
format compact;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dt = 20;
[time10 time20 X Y Xd Yd Xdd Ydd] = fSimulation(dt);
%Subroutine fSimulation is the same as Appendix A with dt as the
parameter.
S = size(time20); m20 = S(2);
S = size(time10); m10 = S(2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
REGIONS = [0 75600 151200 226800 302400]; %4 regions
%Initial No. of zones in each region
nzones1 = 2; nzones2 = 2; nzones3 = 2; nzones4 = 2;
iREGIONS = 1/dt*REGIONS; S = size(iREGIONS); n = S(2);
for i = 1:n-1
diff(i) = iREGIONS(i+1)-iREGIONS(i);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Max_accept = 10^-3;
p = Max_accept + 1;
Max1 = p; Max2 = p; Max3 = p; Max4 = p;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while Max1 >= Max_accept || Max2 >= Max_accept || Max3 >= Max_accept ||
Max4 >= Max_accept
%%%% SPLINE POINTS %%%%%
[x1 y1 yPR1] =
fSplinePoints(time20,diff(1),iREGIONS(1),iREGIONS(2),nzones1,X,Xd,1);
[x2 y2 yPR2] =
fSplinePoints(time20,diff(2),iREGIONS(2),iREGIONS(3),nzones2,X,Xd,2);
[x3 y3 yPR3] =
fSplinePoints(time20,diff(3),iREGIONS(3),iREGIONS(4),nzones3,X,Xd,2);
[x4 y4 yPR4] =
fSplinePoints(time20,diff(4),iREGIONS(4),iREGIONS(5),nzones4,X,Xd,2);
x = [x1 x2 x3 x4];
y = [y1 y2 y3 y4];
yPR = [yPR1 yPR2 yPR3 yPR4];
%%%%% CALCULATING CURVE FIT %%%%%
[iX_fit iXd_fit iXdd_fit] =
fCurveFit_HS(time10,time20,x,y,yPR,Xdd);
91
k = 1;
for i = 1:2:m10
X_fit(k) = iX_fit(i);
k = k+1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X1 = X(iREGIONS(1)+1:iREGIONS(2)+1);
X_fit1 = X_fit(iREGIONS(1)+1:iREGIONS(2)+1);
X2 = X(iREGIONS(2)+2:iREGIONS(3)+1);
X_fit2 = X_fit(iREGIONS(2)+2:iREGIONS(3)+1);
X3 = X(iREGIONS(3)+2:iREGIONS(4)+1);
X_fit3 = X_fit(iREGIONS(3)+2:iREGIONS(4)+1);
X4 = X(iREGIONS(4)+2:iREGIONS(5)+1);
X_fit4 = X_fit(iREGIONS(4)+2:iREGIONS(5)+1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S = size(X1); m1 = S(2); S = size(X2); m2 = S(2);
S = size(X3); m3 = S(2); S = size(X4); m4 = S(2);
%%%%% CALCULATING RESIDUALS %%%%%
%%%%% Entire Curve Fit %%%%%
[Residualx_CF RMSEx_CF Max_CF] = fResiduals(X,X_fit,m20);
%%%%% In Each Region %%%%%
[Residualx1 RMSEx1 Max1] = fResiduals(X1,X_fit1,m1);
[Residualx2 RMSEx2 Max2] = fResiduals(X2,X_fit2,m2);
[Residualx3 RMSEx3 Max3] = fResiduals(X3,X_fit3,m3);
[Residualx4 RMSEx4 Max4] = fResiduals(X4,X_fit4,m4);
%%%%% CHECKING FOR MAXIMUM RESIDUAL IN EACH REGION %%%%%
nzones1 = fCheck(Max1,Max_accept,nzones1);
nzones2 = fCheck(Max2,Max_accept,nzones2);
nzones3 = fCheck(Max3,Max_accept,nzones3);
nzones4 = fCheck(Max4,Max_accept,nzones4);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
92
F. SUBROUTINES USED IN THE CODE GIVEN IN APPENDIX E
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x y yPR] = fSplinePoints(time20,diff,m1,m2,nzones,X,Xd,k)
div = diff/nzones;
i = 1;
for j = m1:div:m2
j = round(j);
ix(i) = time20(j+1);
iy(i) = X(j+1);
iyPR(i) = Xd(j+1);
i = i+1;
end
if k == 1
for i = 1:1:nzones+1
x(i) = ix(i);
y(i) = iy(i);
yPR(i) = iyPR(i);
end
else
for i = 2:1:nzones+1
x(i-1) = ix(i);
y(i-1) = iy(i);
yPR(i-1) = iyPR(i);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Residual RMSE Max] = fResiduals(X,X_fit,m)
k = 1;
for i = 1:m
Residual(k) = abs(X(k)-X_fit(k));
k = k+1;
end
RMSE = sqrt(sum((X(:)-X_fit(:)).^2)/(m));
Max = max(Residual);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function nzones = fCheck(Max,Max_accept,nzones)
if Max <= Max_accept
nzones = nzones + 0;
else
nzones = nzones + 1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
93
G. TRIDIAGONAL SOLVER
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% TRIDIAGONAL SOLVER %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x = fTRIDIAG(d,f,b,n)
alpha(1) = d(1);
gamma(1) = f(1)/alpha(1);
for i = 2:1:n-1
alpha(i) = d(i)-f(i-1)*gamma(i-1);
gamma(i) = f(i)/alpha(i);
end
alpha(n) = d(n)-f(n-1)*gamma(n-1);
z(1) = b(1);
for i = 2:1:n
z(i) = b(i)-gamma(i-1)*z(i-1);
end
for i = 1:1:n
c(i) = z(i)/alpha(i);
end
x(n) = c(n);
for i = n-1:-1:1
i;
x(i) = c(i)-gamma(i)*x(i+1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
94
H. MATLAB CODE TO NUMERICALLY SOLVE THE NEARBY PROBLEM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% NUMERICAL SOLUTION TO THE NEARBY PROBLEM %%%%%%
%%%%%%%%%% USING FIFTH-DEGREE HERMITE SPLINES %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;
clear all;
close all;
format long E;
format compact;
%%%%%%%%%%%%%% NUMERICAL SOLUTION TO THE ORIGINAL PROBLEM %%%%%%%%%%%%%
dt = 20;
[time10 time20 SP X Y Xd Yd Xdd Ydd Xmoon Ymoon Xdd_moon Ydd_moon] =
fSimulation(dt);
%Subroutine fSimulation is the same as Appendix A with dt as the
parameter.
%%%%%%%%%%% GENERATING CURVE FIT TO THE NUMERICAL SIMULATION %%%%%%%%%%
%%%%%% USING FIFTH-DEGREE HERMITE SPLINES %%%%%%
[X_fit Xd_fit Xdd_fit] = fCurveFit_HS(time10,time20,X,Xd,Xdd);
[Y_fit Yd_fit Ydd_fit] = fCurveFit_HS(time10,time20,Y,Yd,Ydd);
S10 = size(time10); m10 = S10(2);
S20 = size(time20); m20 = S20(2);
%%%%%%%%%%%%%%%%%%%%%%% CALCULATING SOURCETERMS %%%%%%%%%%%%%%%%%%%%%%%
for k = 1:1:m10
GMm = 4.90266e12; GMe = 3.98600436e14;
rho = [Xmoon(k) Ymoon(k)];
r = [X_fit(k) Y_fit(k)];
d = r-rho;
if k <= SP
acc = -GMe*r/norm(r)^3-GMm*(d/norm(d)^3+rho/norm(rho)^3);
acc1(k) = acc(1); acc2(k) = acc(2);
Sourcetermx(k) = Xdd_fit(k)-acc1(k);
Sourcetermy(k) = Ydd_fit(k)-acc2(k);
else
acc = -GMm*d/norm(d)^3-GMe*(r/norm(r)^3-rho/norm(rho)^3);
acc3(k) = acc(1); acc4(k) = acc(2);
Sourcetermx(k) = Xdd_fit(k)-Xdd_moon(k)-acc3(k);
Sourcetermy(k) = Ydd_fit(k)-Ydd_moon(k)-acc4(k);
end
end
%%%%%%%%%%%%%%% NUMERICAL SOLUTION TO THE NEARBY PROBLEM %%%%%%%%%%%
%%%%% PROPOGATION IN THE EC->MC FRAME %%%%%
%%%%% SETTING UP THE PROBLEM %%%%%
95
%%%%% USER INPUTS %%%%%
alt = 359750; %meters
angle = -36.890; %degrees
fprintf('alt = %6.0f, ang = %6.3f, dt = %6.3f \n\n', alt, angle, dt);
deltav = 3102.13; %m/s
t0 = 0; %secs
tf = 3.5*86400; %302400 secs (1 day = 86400 secs)
SP = 2.905*86400; %Switch Point at 250992 secs
folder = char(['alt' int2str(alt) 'ang' int2str(abs(angle*1000)) 'dt'
int2str(dt)]);
fid = fopen([folder 'State_SP' '.txt'], 'w');
f = fopen([folder 'Moon_State_SP' '.txt'], 'w');
%%%%% CONSTANTS %%%%%
radE = 6372797; %meters
EMdist = 384400000; %meters
omega = sqrt((GMe+GMm)/EMdist^3); %rad/sec
%Tm = 2*pi/omega = 27d 6h 49m 50.34879957310977s
alpha = angle*pi/180; %radians
h = radE + alt; %meters
Vorbit = sqrt(GMe/h); %m/s
V = Vorbit + deltav; %m/s
%%%%% INITIAL CONDITIONS %%%%%
%%%%% For Spacecraft %%%%%
rxi = h*sin(alpha); %meters
ryi = -h*cos(alpha); %meters
vxi = V*cos(alpha); %m/s
vyi = V*sin(alpha); %m/s
x = [rxi;ryi;vxi;vyi];
fprintf(fid, '% 6.0f, % .16e, % .16e, % .16e, % .16e nn', t0, x);
%%%%% For the Moon %%%%%
xmt = EMdist*cos(omega*t0);
ymt = EMdist*sin(omega*t0);
fprintf(f, '% 6.0f, % .16f, % .16f nn', t0, xmt, ymt);
%%%%% PERFORM THE PROPOGATION %%%%%
%%%%% EARTH CENTERED FRAME %%%%%
i = 1; m = 1;
for t = t0:dt:SP-dt
%Use the moon position at the end of the previous interval as the
%position at the start of this interval.
xm = xmt;
ym = ymt;
%Calculate the moon position at the middle and end of this
interval.
xmh = EMdist*cos(omega*(t + 0.5*dt));
ymh = EMdist*sin(omega*(t + 0.5*dt));
xmt = EMdist*cos(omega*(t + dt));
ymt = EMdist*sin(omega*(t + dt));
t1 = t + dt;
fprintf(f, '% 6.0f, % .16f, % .16f nn', t1, xmt, ymt);
%Propogating the states forward in time
k1 = RK4_EC_NP(x, xm, ym, GMm, GMe, Sourcetermx(m),
96
Sourcetermy(m));
k2 = RK4_EC_NP(x+(dt/2)*k1', xmh, ymh, GMm, GMe, Sourcetermx(m+1),
Sourcetermy(m+1));
k3 = RK4_EC_NP(x+(dt/2)*k2', xmh, ymh, GMm, GMe, Sourcetermx(m+1),
Sourcetermy(m+1));
k4 = RK4_EC_NP(x+(dt)*k3', xmt, ymt, GMm, GMe, Sourcetermx(m+2),
Sourcetermy(m+2));
x = x + dt/6*(k1'+2*k2'+2*k3'+k4');
fprintf(fid, '% 6.0f, % .16e, % .16e, % .16e, % .16enn', t1, x);
i = i+1;
m = m+2;
end
%%%%% SWITCH POINT %%%%%
%Convert state vector from EC to MC
rho = [xmt; ymt];
xmt_dot = -EMdist*omega*sin(omega*(t+dt));
ymt_dot = EMdist*omega*cos(omega*(t+dt));
rho_dot = [xmt_dot; ymt_dot];
d_state = x - [rho; rho_dot];
%%%%% MOON CENTERED FRAME %%%%%
j = 1; n = m;
for t = SPt:dt:tf-dt
%Use the moon position at the end of the previous interval as the
%position at the start of this interval.
xm = xmt;
ym = ymt;
%Calculate the moon position at the middle and end of this
interval.
xmh = EMdist*cos(omega*(t + 0.5*dt));
ymh = EMdist*sin(omega*(t + 0.5*dt));
xmt = EMdist*cos(omega*(t + dt));
ymt = EMdist*sin(omega*(t + dt));
t2 = t + dt;
fprintf(f, '% 6.0f, % .16f, % .16f nn', t2, xmt, ymt);
%Propogating the states forward in time
k1 = RK4_MC_MNP(d_state, xm, ym, GMm, GMe, Sourcetermx(n),
Sourcetermy(n));
k2 = RK4_MC_MNP(d_state+(dt/2)*k1', xmh, ymh, GMm, GMe,
Sourcetermx(n+1), Sourcetermy(n+1));
k3 = RK4_MC_MNP(d_state+(dt/2)*k2', xmh, ymh, GMm, GMe,
Sourcetermx(n+1), Sourcetermy(n+1));
k4 = RK4_MC_MNP(d_state+(dt)*k3', xmt, ymt, GMm, GMe,
Sourcetermx(n+2), Sourcetermy(n+2));
d_state = d_state + dt/6*(k1'+2*k2'+2*k3'+k4');
%Compute moon state at the end of interval
rho = [xmt; ymt];
rho_dot = EMdist*omega*[-sin(omega*(t + dt)); cos(omega*(t + dt))];
%Convert state vector from MC to EC
x = d_state + [rho; rho_dot];
fprintf(fid, '% 6.0f, % .16e, % .16e, % .16e, % .16enn',t2,x);
j = j+1; n = n+2;
end
fclose(fid); fclose(f);
97
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I. SUBROUTINES USED IN THE CODE GIVEN IN APPENDIX H
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% SUBROUTINE TO INTEGRATE THE EQUATION OF MOTION IN EC FRAME %%%%
%%%%%%% FOR THE NEARBY PROBLEM %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x_dot =
RK4_EC_NP(x,xmoon,ymoon,GMm,GMe,Sourcetermx,Sourcetermy)
r = x(1:2);
v = x(3:4);
rho = [xmoon; ymoon];
d = r-rho;
Sourceterm = [Sourcetermx; Sourcetermy];
x_dot(1:2) = v;
x_dot(3:4) = -GMe*r/norm(r)^3-
GMm*(d/norm(d)^3+rho/norm(rho)^3)+Sourceterm;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% SUBROUTINE TO INTEGRATE THE EQUATION OF MOTION IN MC FRAME %%%%
%%%%%%% FOR THE NEARBY PROBLEM %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d_dot =
RK4_MC_NP(d_state,xmoon,ymoon,GMm,GMe,Sourcetermx,Sourcetermy)
d = d_state(1:2);
v = d_state(3:4);
rho = [xmoon; ymoon];
r = d+rho;
Sourceterm = [Sourcetermx; Sourcetermy];
d_dot(1:2) = v;
d_dot(3:4) = -GMm*d/norm(d)^3-GMe*(r/norm(r)^3-
rho/norm(rho)^3)+Sourceterm;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%