Project topics and materials
Get Complete Project Material File(s) Now! »
Numerical integrators
There are two general ways to solve system (1.2.1) of second-order ordinary differential equations. One way is to solve the second-order ordinary differential 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 differential 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 efficiency of different types of integrators. We perform numerical experiments for the different 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 yref 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.
Different 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-off error, and the interpolation error. When performing accurate simulations, the round-off 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-off 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-off error is systematic, the power laws become t and t2, respectively.
First we define the types of errors used in this thesis. We use the notation ytrue(t) to denote the true solution of (1.2.1) and ynum(t) to denote its approximate solution. The difference 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 ||ynum(t) − ytrue(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 ynum(t) for t ≥ 0. This leads to numerical solutions yn = ynum(tn) and yn′ = ynum′(tn) at times tn = t0 + 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 un′′ = f (t, un), un(tn−1) = yn−1, un′(tn−1) = yn′−1, t ∈ [tn−1, tn],
where the function un(t) is the true local solution on the nth interval. The initial conditions yn−1 and yn′−1 are the values of the approximate solution at the end of the (n − 1)st step. Hence, for the local problem the differential equations are the same as for the original problem, but the initial conditions depend upon the numerical solution. Since, in general, the approximation yn−1 is not equal to the true solution ytrue(tn−1), there is a difference between the error made in this local step and the global error. The norm of the local error is defined as ||ynum(t) − un(t)||2, where t ∈ [tn−1, tn]. 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), yref (t), and an underlying polynomial approximation introduced on each step, Pn(t). We then estimate global error for t ∈ [tn−1, tn] by evaluating Pn(t) − yref (t) at k evenly spaced sample points over the step. A more detailed description of yref (t), k, Pn(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 ∈ [t0, tf ], where tf is the pre-specified final value of t in the integration, is more difficult to define and estimate, because we need a continuous approximation to y(t), and not just the discrete values y1, y2, … . We obtain the continuous approximation by first forming a polynomial Pn(t) that approximates the local solution on the interval [tn−1, tn]. The numerical approximation is then the piecewise-defined function accurate and have sufficient 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 [tn−1, tn] 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 effect of different 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 sufficiently 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 tf and this is possible with N-body simulations; if this occurs, the percentage difference between different 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, mi the mass of the ith body, vi,num(t) the numerical approximation to velocity of the ith body (components 3i − 2 to 3i of ynum′(t)), and dij (t) = ||ri,num(t) − rj,num(t)||2 is the distance between the ith and jth bodies. Here, ri,num(t) is formed by the components 3i−2 to 3i of ynum(t) and rj,num(t) by the components 3j − 2 to 3j of ynum(t), which represent the numerical approximations to the positions of the ith and jth 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
