# PERFORMANCE OF ORDINARY DIFFERENTIAL EQUATION SOLVERS FOR THE SIMULATION OF SUN-MERCURY SYSTEM.

Byline: A. Amin, S. Rehman, A. Pervaiz and I. GhafoorABSTRACT

N-body simulations are frequently used to understand the evolution of the Solar System. A wide range of numerical integrators have been developed for performing such simulations; see, for example,[1, 2].The primary objective of this paper is to analyze and compare the error growth and efficiency of different ODE solvers applied to the Kepler's two-body problem. Throughout this paper, the error growth is examined in terms of the global error in the position and velocity, and the relative error in terms of total energy and angular momentum of the system. We performed numerical experiments for different ODE solvers applied to the Kepler problem for Sun-Mercury System with local error tolerances ranging from to .

Keywords: N-body simulations, Kepler problem, ODE Solvers

1. INTRODUCTION

Computational astronomers widely use N-body simulations to understand the development of the Solar System. The developments of the Solar System include the dynamics of the planets, small celestial bodies and asteroids. N-body simulations are performed by deriving a collection of second order ordinary differential equations, and specifying the initial positions and initial velocities of the bodies at time t t0 . The initial value problem (IVP) for N-body simulation is of the form Eqs. where Eqs. represent the initial positions and Eqs.velocities, the operator denotes the differentiation with respect to , and Eqs. is a smooth function,

where is the dimension of the IVP. In some cases of simulations, may be changed as bodies are added or removed from the simulation. In celestial mechanics, motion of two bodies under the influence of gravitational attraction is the simplest form of N-body problem. Large numbers of numerical integrators [3, 4] and ODE solvers for non stiff problems [5, 6, 7] are used to find numerical approximations of differential equations at Eqs. , with and time-step h, which can depend on .

1.1. Kepler Problem

The Kepler's two-body problem [10, 11] is the oldest and simplest problem in the dynamics of the Solar System, because the exact solution of two-body problem exists. Kepler's two-body problem involves the motion of one body about another under the influence of their mutual gravitational force. The difference of masses and the large distance between the planets allows the orbits of most planets and test particles to be approximated by two-body motion. In this research we used Sun-Mercury System as the test problem. Mercury [12, 13] is a planet in the Solar System that is closest to the Sun. The equations of motion of Kepler problem can be written as

Equations

where Eqs. and Eqs. are the x- and y-coordinates of one body relative to the other, Eqs.. The initial conditions are

Equations

The parameter is the orbital eccentricity and . The analytical solution of Kepler problem is available, so the Kepler problem is well suited for observing the accuracy of ODE solvers over a short interval of time. The analytical solution of the Kepler problem is

Equations

where the eccentric anomaly satisfies Kepler's equation Different types of errors have been discussed in this paper. The global error is of main importance in the measurement of the excellence of the numerical solution. We measure this global error in position and velocity and the relative error in energy and angular momentum.

Let Eqs. and Eqs. are the position vectors of the solutions obtained numerically and analytically, respectively, and and are the velocity vectors of numerical and analytical solutions, respectively. The -norm of the global error in position and velocity are given by

Equations

where, is the -norm. Physical systems usually have conserved quantities, like, the total energy and the angular momentum . Generally, these quantities will not be conserved accurately by the numerical solution and this digression provides evaluation about the accuracy of the numerical solution. The total energy of the two-body system is defined as

Equations

where Eqs. is the total energy at the initial time Eqs. . The total angular momentum Eqs. is defined as

Equation

The relative error in angular momentum is defined as

Equation

where Eqs. is the angular momentum at the initial time Eqs. . The analytical solution is not required to calculate . and Eqs. , unlike the global errors in the position and velocity. Hence, fewer computing resources are required to observe the performance of the ODE solvers here.

2. MATERIAL AND METHODS

Mathematicians have developed a wide range of numerical integration techniques for solving the ordinary differential equations (ODEs) that correspond to the continuous state of dynamic systems. A general set of fixed-step and variable- step solvers are provided, each of which implements a specific ODE solution method.

MATLAB is a software that is used to solve an extensive variety of problems. The MATLAB ODE suit is a set of codes for solving first order systems of ordinary differential equations for initial value problems and plotting mathematical results of that problem [5].ODE Solvers control the estimated local error for initial value problems. A local error tolerance is specified and, if the estimated error is too large comparative to this tolerance, the step is rejected and a new attempt is made with a smaller step size [8].The codes; ODE23, ODE45, ODE113 are designed to solve non-stiff problems [5].

2.1. ODE23

The ODE23 solver is based on an explicit Runge-Kutta (2, 3) pair of Bogacki and Shampine [6]. The ODE23 may be more efficient than ODE45 at crude tolerances and in the presence of mild stiffness. The ODE23 is a one step solver, used for non-stiff problems. The code ODE23 consist in a four stage pair of embedded explicit Runge-Kutta methods of order 2 and 3 with error control.

2.2. ODE45

The code ODE45 is a popular (4, 5) pair due to Dormand- Prince [7]. The ODE45 consist 6 stage pair of embedded runge-Kutta method of order 4 and 5. The ODE45 is a one step solver for non-stiff problems. In computing , it needs only the solution at the immediately preceding time point. In general ODE45 is the best solver to apply as a first try" for most problems.

2.3. ODE113

The ODE113 is a variable order and variable step solver which uses Adams-Bashforth-Moulton predictor-correctors of order 1 to 13 [9]. The ODE113 may be more efficient than ODE45 at stringent tolerances and when the ODE function is particularly expensive to evaluate. The ODE113 normally needs the solution at several preceding time point to compute the current solution values [5].

3. RESULTS AND DISCUSSION

Here, we observe the global error in the position and velocity, and the relative error in energy and angular momentum for the Sun-Mercury System. The eccentricity (e) of Mercury is 0.21 and the time of the orbital motion to complete one vibration of Mercury is 28p approximately. The experiments are performed using three ODE solvers applied to the Sun- Mercury System over the period of 28p.

The results in Table 1, Table 2 and Table 3 are for experiments performed with ODE23, ODE45 and ODE113, respectively, applied to the Kepler problem over the interval [0, 28p]. These experiments are performed with a sequence of time steps p/2, p/ , p/ and p/ over the interval [0, 28p] with the tolerances ranging from to . This particular selection of the time step is due to the fact that we wish to obtain best accuracy of the integrators ODE23, ODE45 and ODE113. First we evaluate the position and velocity on each of time step using ODE solvers. The values of positions, velocities and times are saved in separate files. Then we calculate the errors in positions and velocities with respect to the analytical solution that we obtain at the stored values of time. Table 1 shows the maximum global error in position ( ) using the integrator ODE23 with the local error tolerances in the range of to . We observe that at a combination of tolerance and p/2, the maximum global error

Table 1. Maximum global error in position for ODE23 at different tolerances and time steps

###/2###0.342912###0.0362790###0.003665###0.000367###3.6771E-05###3.6779E-06###3.6782E-07###3.6781E-08###3.6778E-09

###/4###0.342913###0.0362799###0.003666###0.000368###3.6763E-05###4.5376E-06###1.2275E-06###8.9652E-07###8.6342E-07

###/8###0.342913###0.0362799###0.003666###0.000368###3.6763E-05###4.5376E-06###1.2275E-06###8.9652E-07###8.6342E-07

###/16###0.363613###0.0572880###0.024676###0.021378###0.0210476###0.0210145###0.0210112###0.0210109###0.0210109

Table 2. Maximum global error in position for ODE45 at different tolerances and time steps

/2###1.52888###0.11872###0.00339###1.73E-05###1.1259E-05###1.6579E-06###1.8758E-07###1.9620E-08###1.9929E-09

/4###1.52888###0.11872###0.00339###1.64E-05###1.2119E-05###2.5177E-06###1.0473E-06###8.7936E-07###8.6173E-07

/8###1.52888###0.11872###0.00339###1.64E-05###1.2119E-05###2.5177E-06###1.0473E-06###8.7936E-07###8.6173E-07

/16###1.51627###0.09775###0.01761###0.020993###0.0210221###0.0210125###0.0210111###0.0210109###0.0210109

Table 3. Maximum global error in position for ODE113 at different tolerances and time steps

/2###0.13984###0.00358###0.00025###3.00E-05###6.8574E-06###2.0444E-06###1.8946E-07###1.0664E-08###5.8392E-09

/4###0.13984###0.00358###0.00025###2.91E-05###5.9977E-06###1.1847E-06###6.7028E-07###8.7040E-07###8.5390E-07

/8###0.13984###0.00358###0.00025###2.91E-05###5.9977E-06###1.1847E-06###6.7028E-07###8.7040E-07###8.5390E-07

/16###0.16075###0.02459###0.02125###0.020981###0.0210040###0.0210088###0.02101072###0.0210109###0.0210109

In the position is less than the errors obtained by other three combinations of different time steps with tolerance . Same is the case at other tolerances, the time step p/2 gives least error as compared to other three time steps for ODE23 integrator.

The same sets of experiments describe in Table 2 are performed to obtain the maximum global error in the position ( ), but now using the integrator ODE45 with the local error tolerances ranging from Eqs. to Eqs . We observe that at a combination of tolerance and p/2, the maximum global error in the position is less than the errors obtained by other three combinations of different time steps tolerance . Same is the case at other tolerances, the time step p/2 gives least error as compared to other three time steps for ODE45 integrator Like Table 1 and 2, the same sets of experiments are performed to obtain the maximum global error in position ( ) using ODE113 integrator with the local error tolerances ranging from to .

We observe that the previous conclusion holds, i.e., at a combination of tolerance and p/2 the maximum global error in the position is less than the errors obtained by other three combinations of different time steps with tolerance . Same is the case at other tolerances, the time step p/2 gives least error as compared to other three time steps for ODE113 integrator. From Table 1, Table2 and Table 3, we observe a clear pattern. When the tolerance is increased, the maximum global error in position is also increased. We also observe that accuracy of the given ODE solver improves if the time step is large at tolerance Eqs. We conclude that a combination of Tolerance =Eqs and time step = p/2 gives better results in terms of maximum global error in position for all three integrators.

Furthermore, we have performed experiments to observe the global error in velocity and relative error in energy and angular momentum for ODE23, ODE45 and ODE113 using the combination of tolerance and step size = p/2. We observe that using ODE45 solver the least maximum global error in position is approximately , which is the best observed accuracy. The ODE45 integrator has achieved approximately 84% and 192% better accuracy then ODE23 and ODE113, respectively.

Table 4. Least maximum global error in position and velocity

Solver

ODE45###2.574536075

ODE23###4.761516179

ODE113###5.839207460###7.617207406

We also computed the global error in velocity. From Table 4, we found that least maximum global error in velocity is approximately 1 digit larger than the error in the position for the integrators ODE23, ODE45 and ODE113. Now consider the accuracy of ODE integrators using relative error in the energy and angular momentum.

Figure1: shows the error behavior in total energy using three integrators ODE23, ODE45 and ODE113 applied to the Kepler problem over 14 periods of time for the planet Mercury. The tolerance and time step was selected to give the smallest maximum global error. From Figure 1, we observe that the best observed accuracy in terms of the relative error in energy is again achieved by the ODE45 integrator.

Figure 2: shows the error behavior in angular momentum using three integrators ODE23, ODE45 and ODE113 applied to the Kepler problem over 14 periods of time for the planet Mercury. The tolerance and time step was selected to give the smallest maximum global error. From Figure 2, we observe that the best observed accuracy in terms of the relative error in angular momentum is again achieved by the ODE45 integrator.

Let us now consider the efficiency of the integrators discussed in this paper, which is the amount of work to attain a given accuracy. One way of measuring the amount of work is to count the number of function evaluations. The number of function evaluations against the least maximum global error in position for the ODE45, ODE23 and ODE113 integrator with tolerance and time step p/2 are 39223, 995641 and 3966, respectively. We observe that ODE113 takes the least function evaluations. Whereas, the integrator ODE23 is approximately 25.4 times more expensive than ODE45. Whereas, ODE45 is approximately 9.9 times more expensive than ODE113.

4. CONCLUSION

The main purpose of this paper was to examine the error growth and efficiency for different ODE integrators applied to the real word problem involving the Sun and the planet Mercury. The simulations were performed over one period of Mercury. For these simulations, we performed experiments to observe the error growth in position and velocity using three integrators ODE23, ODE45 and ODE113 for the local error tolerances ranging from to and with variation in time steps. For the given range of tolerances from to , and time steps from p/16 to p/2, we observed that the integrator ODE45 achieves the best observed accuracy. The integrator ODE23 attained the 2nd best accuracy whereas the ODE113 was the least accurate at the same combination.

We then analyzed the efficiency of the integrators by counting the number of function evaluations against the least maximum global error in position. We observed that the best accuracy attained integrator ODE45 uses approximately 9.9 times more function evaluations than ODE113 and approximately 25.4 times less function evaluations than ode23.

REFERENCES

[1]. P.W. Sharp, Requirement of a package for N-body simulation of the solar system," Kluwer Academic Publishers, 31, 271-279(2002).

[2]. S. Rehman, Jovian problem: Performance of some high-order Numerical Integrators," American journal of computational Mathematics, 3(3), 195-204(2013).

[3]. E. S. P. NArsett, and G. Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems," Springer, second edition, (1993).

[4]. J.C. Butcher, A history of Runge Kutta methods," Applied Numerical Mathematics, 20, 247-260(1996).

[5]. R. Ashino, M. Nagase and R.Vaillancourt, Behind and Beyond the MATLAB ODE Suite," Computers and Mathematics with Applications, 40, 491-512(2000).

[6]. P. Bogacki and L. F. Shampine, A 3(2) pair of Runge- Kutta formulas," Appl. Math. Letters, 2(4), 1-9(1989).

[7]. J. R. Dormand and P. J. Prince, A family of embedded Rune-Kutta formulae," J. Comp. Appl. Math., 6(1), 19- 26(1980).

[8]. L.F. Shampine, I. Gladwell, S. Thompson, Solving ODEs with MATLAB," Cambridge university press, (2003).

[9]. L. F. Shampine and M. K. Gordon, Computer Solution of Ordinary Differential Equations: the Initial Value Problem," W. H. Freeman, San Francisco, (1975).

[10]. D. Lucic and M.Varga, Simulation of the two-body problem in geogebra,"12(3), 4750(2012).

[11]. H. Goldstein, Classical Mechanics," Addison Wesley, second edition, (1980).

[12]. Joseph J. Smulsky, Gravitational, field, and rotation of Mercury perihelion,"5(2), 254-260(2008).

[13]. G. Faure. Teresa M. Mensing, Introduction to planetary science: The geological perspective," published by Springer, (2007).

[14]. Muhammad Amer Qureshi, Efficient Numerical Integration for gravitational N-body Simulations," University of Auckland, (2012).

Printer friendly Cite/link Email Feedback | |

Publication: | Science International |
---|---|

Article Type: | Report |

Date: | Apr 30, 2015 |

Words: | 2720 |

Previous Article: | EFFECT OF MAGNETIC FIELD AND DOUBLE DISPERSION ON MIXED CONVECTION HEAT AND MASS TRANSFER IN NON-DARCY POROUS MEDIUM. |

Next Article: | UNSTEADY MHD FLOW AND HEAT TRANSFER FOR NEWTONIAN FLUIDS OVER AN EXPONENTIALLY STRETCHING SHEET. |

Topics: |