An Application Of Modified Equation To Reduce Some Numerical Instabilities

The presence of errors in the initial conditions and the different physical parameters affect to numerical integration of equations of motion necessarily. It is interesting to determine which the causes of the possible phase errors are and/or the dumping/dissipative terms that can appear. The modified equation can be useful to reduce numerical instability errors which may be confused by geometrical or physical errors. This is a necessary first step to detect another kind of errors related to the mathematical/physical model and not to the integration method.


I NTRODUCTI ON
Celestial Mechanics deals with motions of celestial bodies such as dynamics of planetary systems or orbital and rotational motion of planets and stars, but the dynamical system is not integrable and so only numerical methods provide accurate and efficient solutions. This fact explains why the astronomical community has shown great interest in the numerical methods of ordinary differential equation. The classification and description of ordinary methods is widely known (see, for example, [1]) while a review of symplectic methods may be seen in [3]. Numerical stability depends only to the integrator and is independent of the problem that we want to solve. A first approach on this question has already been attempted in [4].
Symplectic integrators are widely used in different fields of Celestial Mechanics, such as the study of resonances or satellite integration. In a discrete level, the conservation of a quantity is not possible without implying any error propagation for any variables of the problem at the same time. For example, let us remember Yoshida's work [7], where no explanation is found for the fact that the symplectic integrator SI4 provides errors which are much larger in magnitude for the argument of the periaster than the ones that would have been obtained if we had used the classic RK4, which is not a symplectic method. The application of the method makes impossible the integration step to tend to zero and this fact, together with the impossibility to adjust the step to the desired local error while the symplecticity is preserved, implies that the practical application depends on each particular problem. A remarkable example is provided by the "leapfrog" scheme in the P.D.E. frame, which is unconditionally stable, but a more careful study reveals that a significant error in phase is produced in time when the method is applied to the linear wave equation. Both examples illustrate the problem of a practical instability linked to a particular problem. In addition, it also shows us how to analyze the errors in more detail and how to modify the scheme in such a way that we can obtain another better scheme. In this sense, let us remember that it is possible to obtain the differential equation that, under some regularity hypothesis and for a given scheme, really verifies it. Yet, this equation is not the initial one but instead is the so called modified equation. The study of this equation makes it possible to modify the initial scheme that is behaving in an unstable way due to some parameter.
Finally, let us think of the numerical integration of a set of orbital elements of a body in the solar system: the metric factors are subjected to amplification/diffusion and the angular factors may suffer from errors in phase. All of these may happen while the theoretical stability of the methods is not affected.
It would be interesting to find out if the observed errors come from the numerical treatment or from the mathematical/physical model (local reference system or physical parameters) Practical examples may be seen, for example, in [2] and [5] (when he cites [6]).
Our objective is to briefly study different symplectic integrators in terms of its application to the cases of the harmonic oscillator and the two-body problem, in what refers to amplitude and phase errors. These integrators will be modified to improve the problems in amplification and phase errors that make its symplecticity less effective.

PHASE ERROR I N SEVERAL METHODS APPLI ED TO THE HARM ONI C OSCI LLATOR
Let us consider the harmonic oscillator equation A first step will be to check that the numerical solution presents in most cases a phase error with respect the true solution. We know that the numerical solution is the exact solution of the equivalent equation and the application of the scheme to the modified equation allows us to see that the phase error increases in a similar magnitude. This fact allows an artificial parameter w to eliminate the former error. As an example, it is immediate to check that the equations of the SI2 method applied to this problem are: being h the time step. The differential equation for which the exact solution is the exact solution of the numerical scheme is called the differential modified equation and takes the form:   (5) In order to illustrate the procedure, in Figure 1a we see the results for the SI2 taking n=4000, h=0.01, m=0.00025, being w=1/m, q 0 =0.5 and p 0 =0. In Fig 1b we

CORRECTI ON TO THE PERI HELI ON ADVANCE PRODUCED BY THE DI FFERENT METHODS
The problems concerning the existence of phase errors in the application of any numerical method to the simple problem of the harmonic oscillator are especially dramatic if the two-body problem is to be integrated, in particular, in which concerns to the precession of the perihelion. To study this, let us consider the Kepler problem: We take the initial conditions q 1 =0.4, q 2 =0, p 1 =0, p 2 =2, =1, and the step h=0.1 with 1200 steps of integrations. Bearing in mind the non-existence of the corresponding first integral for the Runge-Lenz vector, which is expressed by (e 1 ,e 2 ), with the problem increases when considering the modified equation of the method as applied to the two-body problem [8]. In practice, the main difficulty lies in the magnitude of the inverse of the radius vector in the perihelion. This is because in the perihelion the radius is small enough to imply a large magnitude for its inverse, which requires a dramatic reduction of the step h. Hence, a special treatment would be required, with a reduction of the size of step in the perihelion, which would destroy the symplecticity of the method. It should be remembered that the stability of a method and numerical stability are different concepts. The former refers to the case when h tends to 0 (in practice, a 'small' stepsize) and the latter refers to the behaviour of the numerical solution versus the real one, when the number of steps increases with a constant stepsize h. It is obvious that stability is a necessary condition for numerical stability, but it is not enough on its own. In this regard, a first analysis of the modified equati q 2 and the p 1 components. So, a simple solution consists in including a diffusive term in the computation of p 1 , without modifying the symplectic character of the method and also weighting the vector radius obtained from SI2 and its modification.
After these improvements of the method, we can see its effect on the Runge-Lenz vector, shown in the following Figure 2a. We see the precession effect on the Runge-Lenz vector using the SI2 method (crossing the figure) and its rough correction (the errors are bounded on the right top of the figure, inside the little square). In the same way, a improvement in the results concerning the whole orbit is also observed, as can be seen in the following figure  (Figure 2b), where the second order modified integration (lower section of Fig. 2b) may be compared with the provided by the classic four order RK4 (upper part of Fig. 2b).

CONCLUSI ONS
We have seen that the property of being symplectic does not prevent the numerical instability of some methods in either the harmonic oscillator problem or in the two-body problem. We have conducted a study of the different symplectic lower order methods with the aim of adapting them in a specific way for each of the problem to which it is applied. Obviously, this is a first step that must be followed by a more in-depth study which allows for a more systematic treatment of the issues mentioned above, focusing special attention in the non linear problems that appear in celestial mechanics, such as the mathematical pendulum (autonomous and non autonomous). When applying these models to real proble have to be especially careful with the possible existence of errors in phase and in the amplitude for the different coordinates used in the problems of evolution, as shown in the examples. Each problem may require specific coordinates and the numerical schemes have to be adapted to them so that the physical properties are conserved as much as possible and, at the same time, the sensitivity to changes remains small (which means having 'physical' stability through the 'toughness' of the method). In the case of problems that allow both metric and angular coordinates to be used, the kind of errors we are dealing with, amplitude errors or phase errors is fundamental for these purposes. Bearing in mind the above issues, our main objective is to carry out further research.