Get Complete Project Material File(s) Now! »

**Numerical integrators**

There are two general ways to solve system (1.2.1) of second-order ordinary diﬀerential equations. One way is to solve the second-order ordinary diﬀerential equations directly by using numerical integrators such as Runge-Kutta-Nystr¨om, St¨ormer, or extrapolation in-tegrators. The other way is to convert system (1.2.1) into a system of first-order ordinary diﬀerential equations, and then use numerical integrators such as Runge-Kutta, Adams, or extrapolation integrators. In this chapter, in particular, we will survey N-body inte-grators that use the first approach. We will also measure the accuracy and the eﬃciency of diﬀerent types of integrators. We perform numerical experiments for the diﬀerent in-tegrators applied to the Jovian problem over a long interval, as long as 100 million years, with the local error tolerance ranging from 10^{−16} to 10^{−08}.

Reference solution and errors

Unlike the Kepler problem, an analytical solution for the Jovian problem is not available. Therefore, numerical experiments using the Jovian problem require a reference solution y_{ref} in order to obtain an estimate of the error in the position and velocity. The reference solution has to be more accurate than the numerical solution. Since we plan to test the numerical integrators near the limit of double-precision arithmetic (2.2 × 10^{−16}), it is essential to use quadruple-precision arithmetic for the reference solution.

Diﬀerent types of errors are discussed throughout this chapter. The global error is of major importance in the measurement of the quality of the numerical solution. We measure this global error in position and velocity. In addition we measure the relative error in energy and angular momentum, and the phase error.

For the total error in the system there are two main sources of error when an interpolation scheme is used: the integration error, which consists of truncation and round-oﬀ error, and the interpolation error. When performing accurate simulations, the round-oﬀ error can contribute significantly to the global error. For fixed-step-size schemes, Brouwer [4] showed that, if the step-size is smaller than a prescribed value, the round-oﬀ error for conserved quantities, such as total energy and angular momentum, grows as t 2 __3 __and for other dynamical variables, such as the coordinates of particles, as t2 . This error growth is known as Brouwer’s law in the literature; see, for example, [28, 34]. In contrast, when the round-oﬀ error is systematic, the power laws become t and t^{2}, respectively.

First we define the types of errors used in this thesis. We use the notation y_{true}(t) to denote the true solution of (1.2.1) and y_{num}(t) to denote its approximate solution. The diﬀerence is monitored time wise by considering the global error in y as The norm of the global error in y at time t is then ||y_{num}(t) − y_{true}(t)||_{2}.

The calculation of the global error is discussed later. A wide range of integrators, for example, Runge-Kutta [37, 43], linear multistep [42], Runge-Kutta-Nystr¨om [50], and St¨ormer [59] can be used to find y_{num}(t) for t ≥ 0. This leads to numerical solutions y_{n} = y_{num}(t_{n}) and y_{n}^{′} = y_{num}^{′}(t_{n}) at times t_{n} = t_{} + nh, n = 1, 2, …, where the time-step h can depend on n. The approximation over the continuous time interval can then be computed using (local) interpolation.

From one time-step to the next, the local problem is solved, which is defined as u_{n}^{′′} = f (t, u_{n}), u_{n}(t_{n−1}) = y_{n−1}, u_{n}^{′}(t_{n−1}) = y_{n}^{′}_{−1}, t ∈ [t_{n−1}, t_{n}],

where the function u_{n}(t) is the true local solution on the n^{th} interval. The initial conditions y_{n−1} and y_{n}^{′}_{−1} are the values of the approximate solution at the end of the (n − 1)^{st} step. Hence, for the local problem the diﬀerential equations are the same as for the original problem, but the initial conditions depend upon the numerical solution. Since, in general, the approximation y_{n−1} is not equal to the true solution y_{true}(t_{n−1}), there is a diﬀerence between the error made in this local step and the global error. The norm of the local error is defined as ||y_{num}(t) − u_{n}(t)||_{2}, where t ∈ [t_{n−1}, t_{n}]. In general, the global error cannot be calculated because the true solution is not known.

The estimated global error in the position over each time-step, n, is defined in terms of an accurate numerical solution (computed in quadruple precision), y_{ref} (t), and an underlying polynomial approximation introduced on each step, P_{n}(t). We then estimate global error for t ∈ [t_{n−1}, t_{n}] by evaluating P_{n}(t) − y_{ref} (t) at k evenly spaced sample points over the step. A more detailed description of y_{ref} (t), k, P_{n}(t), and the definition of the estimated global error in the position is presented later in this section.

The global error in the position for all t ∈ [t_{}, t_{f} ], where t_{f} is the pre-specified final value of t in the integration, is more diﬃcult to define and estimate, because we need a continuous approximation to y(t), and not just the discrete values y_{1}, y_{2}, … . We obtain the continuous approximation by first forming a polynomial P_{n}(t) that approximates the local solution on the interval [t_{n−1}, t_{n}]. The numerical approximation is then the piecewise-defined function accurate and have suﬃcient continuity. These requirements along with the types of local interpolants we used in this thesis are discussed in detail in Chapter 3.

We store the information, such as positions and times, in separate files. In a post-processing program, we estimate the maximum global error by sampling the norm at evenly spaced data points on each sub-interval [t_{n−1}, t_{n}] with respect to the reference solution, which we obtain at these stored values of time by forcing the integrator to hit the sample points, and then taking the maximum of these norms. To see the eﬀect of diﬀerent values of k on the estimated maximum global error in position we performed experiments using the combination of ERKN101217 and the 23-stage interpolant, which are described in Sections 2.2.1.1 and 3.1.2. The integration was performed using the Jovian problem over one million years for the local error tolerances 10^{−08}, 10^{−10}, 10^{−12}, 10^{−14}, and 10^{−16}. We found that k = 10 led to an estimated maximum global error that was suﬃciently accurate for our work; this is illustrated by the results in Table 2.1. The rows labelled k = 50 and k = 100 list the percentage changes in the estimated maximum global error when compared with the values in row labelled k = 10.

We observe from Table 2.1 that the data for T OL = 10^{−08} shows a percentage change of zero. In this case, the estimated maximum global error occurs at t_{f} and this is possible with N-body simulations; if this occurs, the percentage diﬀerence between diﬀerent values of k would be zero. We repeated the same set of experiments with other combinations. For example, when using the integrator ODEX2 with its interpolant, which are described in Sections 2.2.3.1 and 3.1.3, an increase of k from 10 to 100 changed the estimated maximum global error by not more than 1%. Throughout the remainder of this thesis, all maximum global errors are estimated by sampling at 10 evenly spaced data points on every time-step.

Appropriate combinations of integrators and interpolation schemes are essential to main-tain accuracy; for more details on appropriate combinations, see Chapter 3. The ex-periments, for example, as reported in Table 2.3 as well as figures 2.1, 2.2, 2.4, and 2.5 are performed with local interpolation schemes based on definition of the maximum global error. Here, we have used the 12-stage, 23-stage, built-in interpolants, and quintic Hermite interpolation with integrators ERKN689, ERKN101217, ODEX2, and St¨ormer, respectively.

Physical systems often have conserved quantities, for example, the total energy H or the total angular momentum L. Usually, these quantities will not be conserved exactly by the numerical solution and this deviation provides assessment about the accuracy of the solution. The total energy for a system of N bodies interacting with one another through Newtonian forces is defined as where G is the gravitational constant, m_{i} the mass of the i^{th} body, v_{i,num}(t) the numerical approximation to velocity of the i^{th} body (components 3i − 2 to 3i of y_{num}^{′}(t)), and d_{ij} (t) = ||r_{i,num}(t) − r_{j,num}(t)||_{2} is the distance between the i^{th} and j^{th} bodies. Here, r_{i,num}(t) is formed by the components 3i−2 to 3i of y_{num}(t) and r_{j,num}(t) by the components 3j − 2 to 3j of y_{num}(t), which represent the numerical approximations to the positions of the i^{th} and j^{th} bodies, respectively.

Abstract

Acknowledgments

Symbols and abbreviations

1 Introduction

1.1 Background

1.2 Overview

1.3 Imaginary sphere

1.4 Test problems

1.5 Structure of the Thesis

2 Numerical integrators

2.1 Reference solution and errors

2.2 Numerical methods and integrators

2.3 Step-size variation

2.4 Numerical experiments for long-term simulation

3 Accuracy and efficiency of interpolation schemes

3.1 Continuous approximations methods

3.2 Numerical experiments for long-term simulations .

4 Algorithm and initial conditions

4.1 Core Algorithm

4.2 Other integration details

4.3 Initial conditions of asteroids

5 N-Body simulations

5.1 Selection of quality asteroids

5.2 Close-encounter error

5.3 Close-encounters at a given numerical accuracy

5.4 Best observed numerical accuracy

6 Summary

References

GET THE COMPLETE PROJECT