On simulations of the classical harmonic oscillator equation by difference equations

We show that any second order linear ordinary diffrential equation with constant coefficients (including the damped and undumped harmonic oscillator equation) admits an exact discretization, i.e., there exists a difference equation whose solutions exactly coincide with solutions of the corresponding differential equation evaluated at a discrete sequence of points (a lattice). Such exact discretization is found for an arbitrary lattice spacing.


Introduction
The motivation for writing this paper is an observation that small and apparently not very important changes in the discretization of a differential equation lead to difference equations with completely different properties.
By the discretization we mean a simulation of the differential equation by a difference equation [1].
In this paper we consider the damped harmonic oscillator equation where x = x(t) and the dot means the t-derivative. This is a linear equation and its general solution is well known. Therefore discretization procedures are not so important (but sometimes are applied, see [2]). However, this example allows us to show and illustrate some more general ideas. The most natural discretization, known as the Euler method (Appendix B, compare [1,3]) consists in replacing x by x n ,ẋ by the difference ratio (x n+1 − x n )/ε,ẍ by the difference ratio of difference ratios, i.e., and so on. This possibility is not unique. We can replace, for instance, x by x n+1 , orẋ by (x n − x n−1 )/ε, orẍ by (x n+1 − 2x n + x n−1 )/ε 2 . Actually the last formula, due to its symmetry, seems to be more natural than (2) (and it works better indeed, see Section 2). In any case we demand that the continuum limit, i.e., x n = x(t n ) , t n = εn , ε → 0 , applied to any discretization of a differential equation yields this differential equation. The continuum limit consists in replacing x n by x(t n ) = x(t) and the neighbouring values are computed from the Taylor expansion of the function x(t) at t = t n x n+k = x(t n + kε) = x(t n ) +ẋ(t n )kε + 1 2ẍ (t n )k 2 ε 2 + . . .
Substituting these expansions into the difference equation and leaving only the leading term we should obtain the considered differential equation.
In this paper we compare various discretizations of the damped (and undamped) harmonic oscillator equation. The main result of the present paper consists in finding the exact discretization of the damped harmonic oscillator equation (1). By exact discretization we mean that x n = x(t n ) holds for any ε and not only in the limit (3).

Simplest discretizations of the harmonic oscillator
Let us consider the following three discrete equations where ε is a constant. The continuum limit (3) yields, in any of these cases, the harmonic oscillator equation To fix our attention, in this paper we consider only solutions corresponding to the initial conditions x(0) = 0,ẋ(0) = 1. The initial data for the discretizations are chosen in the simplest form: we assume that x 0 and x 1 belong to the graph of the exact continuous solution. For small t n and small ε the discrete solutions of any of these equations approximate the corresponding continuous solution quite well (see Fig. 1). However, the global behaviours of the solutions (still for small ε!) are different (see Fig. 2). The solution of the equation (6) vanishes at t → ∞ while the solution of (4) oscillates with rapidly incrising amplitude (all black points are outside the range of Fig. 2). Qualitatively, only the discretization (5) resembles the continuous case. However, for very large t the discrete solution becomes increasingly different from the exact continuous solution even in the case (5) (compare Fig. 2 and Fig. 3).
The natural question arises how to find a discretization which is the best as far as global properties of solutions are concerned.
In this paper we will show how to find the "exact" discretization of the damped harmonic oscillator equation. In particular, we will present the discretization of (7) which is better than (5) and, in fact, seems to be the best possible. We begin, however, with a very simple example which illustrates the general idea of this paper quite well.  3 The exact discretization of the exponential growth equation We consider the discretization of the equationẋ = x. Its general solution reads The simplest discretization is given by This discrete equation can be solved immediately. Actually this is just the geometric sequence x n+1 = (1 + ε)x n . Therefore To compare with the continuous case we write (1 + ε) n in the form where t n := εn and κ := ε −1 ln(1+ε). Thus the solution (10) can be rewritten as Therefore we see that for κ = 1 the continuous solution (8), evaluated at t n , i.e.
Although the qualitative behaviour of the "naive" simulation (9) is in good agreement with the smooth solutions (exponential growth in both cases) but quantitatively the discrepancy is very large for t → ∞ because the exponents are different.

Discretizations of the harmonic oscillator: exact solutions
The general solution of the harmonic oscillator equation (7) is well known In Section 2 we compared graphically this solution with the simplest discrete simulations: (4), (5), (6). Now we are going to present exact solutions of these discrete equations.
Because the discrete case is usually less known than the continuous one we recall shortly that the simplest approach consists in searching solutions in the form x n = Λ n (this is an analogue of the ansatz x(t) = exp(λt) made in the continuous case, for more details see Appendix A). As a result we get the characteristic equation for Λ.
We illustrate this approach on the example of the equation (4) resulting from the Euler method. Substituting x n = Λ n we get the following characteristic equation with solutions Λ 1 = 1 + iε, Λ 2 = 1 − iε. The general solution of (4) reads x n = c 1 Λ n 1 + c 2 Λ n 2 , and, expressing c 1 , c 2 by the initial conditions x 0 , x 1 , we have where ρ = √ 1 + ε 2 and α = arctan ε, we obtain after elementary calculations It is convenient to denote ρ = e κε and t n = nε , κ := 1 2ε and then One can check that κ > 0 and ω < 1 for any ε > 0. For ε → 0 we have κ → 0, ω → 1. Therefore the discrete solution (20) is characterized by the exponential growth of the envelope amplitude and a smaller frequency of oscillations than the corresponding continuous solution (15). A similar situation is in the case (6), with only one (but very important) difference: instead of the growth we have the exponential decay. The formulas (19) and (20) need only one correction to be valid in this case. Namely, κ has to be changed to −κ.
The third case, (5), is characterized by ρ = 1, and, therefore, the amplitude of the oscillations is constant (this case will be discussed below in more detail).
These results are in perfect agreement with the behaviour of the solutions of discrete equations illustrated at Fig. 1 and Fig. 2.
Let us consider the following family of discrete equations (parameterized by real parameters p, q): The continuum limit (3) applied to (21) yields the harmonic oscillator (7) for any p, q. The family (21) contains all three examples of Section 2 and (for p = q = 1/4) the equation resulting from the Gauss-Legendre-Runge-Kutta method (see Appendix B): Substituting x n = Λ n into (21) we get the following characteristic equation: We formulate the following problem: find a discrete equation in the family (21) with the global behaviour of solutions as much similar to the continuous case as possible.
At least two conditions seem to be very natural in order to get a "good" discretization of the harmonic oscillator: oscillatory character and a constant amplitude of the solutions (i.e., ρ = 1, κ = 0). These conditions can be easily  expressed in terms of roots (Λ 1 , Λ 2 ) of the quadratic equation (23). First, the roots should be imaginary (i.e., ∆ < 0), second, their modulus should be equal to 1 i.e., Λ 1 = e iα , Λ 2 = e −iα . Therefore 1 + pε 2 = 1 + qε 2 , i.e., q = p. In the case q = p the discriminant ∆ of the quadratic equation (23) is given by There are two possiblities: if p ≥ 1/4, then ∆ < 0 for any ε = 0, and if p < 1/4, then ∆ < 0 for sufficiently small ε, namely ε 2 < 4(1 − 4p) −1 . In any case, these requirements are not very restrictive and we obtained p-family of good discretizations of the harmonic oscilltor. If Λ 1 = e iα and Λ 2 = e −iα , then the solution of (21) is given by where ω = α/ε, i.e., Note that the formula (24) is invariant with respect to the transformation α → −α which means that we can choose as Λ 1 any of the two roots of (23). The equation (5) is a special case of (21) for q = p = 0. As we have seen in Section 2, for small ε this discretization simulates the harmonic oscillator (7) much better than (4) or (6). However, for sufficiently large ε (namely, ε > 2) the properties of this discretization change dramatically. Its generic solution grows exponentially without oscillations.
Expanding (25) in the Maclaurin series with respect to ε we get Therefore the best approximation of (7) from among the family (21) is characterized by p = 1/12: In this case ω ≈ 1 + ε 4 /480 + . . . is closest to the exact value ω = 1.
The standard numerical methods give similar results (in all cases presented in Appendix B the discretization of the second derivative is the simplest one, the same as described in our Introduction). The corresponding discrete equations do not simulate (7) better than the discretizations presented in Section 2.

Damped harmonic oscillator and its discretization
.
To obtain some simple discretization of (1) we should replace the first derivative and the second derivative by discrete analogues. The results of Section 2 suggest that the best way to discretize the second derivative is the symmetric one, like in Eq. (5). There are at least 3 possibilities for the discretization of the first derivative leading to the following simulations of the damped harmonic oscillator equation: As one could expect, the best simulation is given by the most symmetric equation, i.e., Eq. (31), see  6 The exact discretization of the damped harmonic oscillator equation In order to find the exact discretization of (1) we consider the general linear discrete equation of second order The general solution of (33) has the following closed form (compare Appendix A): where Λ 1 , Λ 2 are roots of the characteristic equation Λ 2 − 2AΛ − B = 0, i.e., The formula (34) is valid for Λ 1 = Λ 2 , which is equivalent to Is it possible to identify x n given by (34) with x(t n ) where x(t) is given by (28)?
Yes! It is sufficient to express in an appropriate way λ 1 and λ 2 by Λ 1 and Λ 2 and also the initial conditions x(0),ẋ(0) by x 0 , x 1 . It is quite surprising that the above identification can be done for any ε.
The crucial point consists in setting Λ n k = exp(n ln Λ k )) = exp(t n λ k ) , where, as usual, t n := nε. It means that (note that for imaginary Λ k , say Λ k = ρ k e iα k , we have ln Λ k = ln ρ k + iα k ). Then (34) assumes the form Comparing (28) with (39) we get x n = x(t n ) provided that The degenerate case, Λ 1 = Λ 2 (which is equivalent to λ 1 = λ 2 ) can be considered analogically (compare Appendix A). The formula (36) is obtained from (34) in the limit Λ 2 → Λ 1 . Therefore all formulas for the degenerate case can be derived simply by taking the limit λ 2 → λ 1 .
Thus we have a one-to-one correspondence between second order differential equations with constant coefficients and second order discrete equations with constant coefficients. This correspondence, referred to as the exact discretization, is induced by the relation (38) between the eigenvalues of the associated characteristic equations.
The damped harmonic oscillator (1) corresponds to the discrete equation (33) such that In the case of the weakly damped harmonic oscillator (ω 0 > γ > 0) the exact discretization is given by where ω := ω 2 0 − γ 2 . In other words, the exact discretization of (1) reads The initial data are related as follows (see (40)): x(0) = x 0 ,ẋ(0) = x 1 ωe γε − (γ sin(ωε) + ω cos(ωε))x 0 sin(ωε) , Fig. 6 and Fig. 7 compare our exact discretization with two other good discretizations of the weakly damped harmonic oscillator. The exact discretization is really exact, i.e., the discrete points belong to the graph of the exact continuous solution (for any ε and any n). Similarly as in the undumped case, the fully symmetric simple discretization (31) is better than the equation resulting from GLRK method. The exact discretization of the harmonic oscillator equationẍ + x = 0 is a special case of (43) and is given by It is easy to verify that the formula (45) can be rewritten as which reminds the "symmetric" version of Euler's discretization scheme (see (2) and (5)) but ε appearing in the discretization of the second derivative is replaced by 2 sin(ε/2). For small ε we have 2 sin(ε/2) ≈ ε.
The comparison of the exact discretization (45) with three other discrete equations simulating the harmonic oscillator is done in Fig. 3, and in Fig. 4. We point out that the considered simulations are very good indeed (although t in Fig 3 is very large) but they cannot be better than the exact discretization. The discretization (27) is also excellent. The coefficient by −2x n in Eq. (27) 12 − 5ε 2 12 + ε 2 ≈ 1 − approximates cos ε up to 4th order. Actually, for the choice of parameteres made in Fig. 3 and Fig. 4, the discretization (27) practically cannot be discerned from the exact one.

Conclusions
In this paper we have shown that for linear ordinary differential equations of second order with constant coefficients there exists a discretization which is "exact" and simulates properly all features of the differential equation. The solutions of this discrete equation exactly coincide with solutions of the corresponding differential equation evaluated at a discrete lattice. Such exact discretization can be found for an arbitrary lattice spacing ε. Therefore we conclude that in this case differential and difference equations are in one-to-one correspondence: to any linear differential equation with constant coefficients there corresponds a difference equation which we call the exact discretization.
Analogical considerations can be made for linear ordinary differential equations (with constant coefficients) of any order (the details will be presented elsewhere).
We point out that to achieve our goal we had to assume an essential dependence of the discretization on the considered equation, in contrast to the standard numerical approach to ordinary differential equations where practically no assumptions are imposed on the considered system (i.e., universal methods, applicable for any equation, are considered) [4].
In the last years one can observe the development of the numerical methods which are applicable for some clases of equations (e.g., admitting Hamiltonian formulation) but are much more powerful (especially when global or qualitative aspects are considered) [5,6].
"Recent years have seen a shift in paradigm away from classical considerations which motivated the construction of numerical methods for ordinary differential equations. Traditionally the focus has been on the stability of difference schemes for dissipative systems on compact time intervals. Modern research is instead shifting in emphasis towards the preservation of invariants and the reproduction of correct qualitative features" [7].
In our paper we have a kind of extremal situation: the method is applicable for a very narrow class of equations but as a result we obtain the discretization which seems to be surprisingly good.
Similar situation occurs for the integrable (soliton) nonlinear systems. It is believed (and proved for a very large class of equations) that integrable equations admit integrable discretizations which preserve the unique features of these equations (infinite number of conservation laws, solitons, transformations generating explicit solutions etc.) [8,9].
The exact discretization considered in this paper is the best possible simulation of a differential equation. Linear ordinary differential equations always admit the unique exact discretization. An open problem is to find such discretization for some other classes of differential equations.
just the eigenvectors of M, and the diagonal diagonal coefficients of D are eigenvalues of M.
The characteristic equation (det(M − λI) = 0) for m = 2 (i.e., for (47)) has the form Its roots will be denoted by Λ 1 , Λ 2 (see (35)). If Λ 1 = Λ 2 , then the diagonalization procedure yields where the columns of N are the eigenvectors of M, i.e., Therefore and performing the multiplication we get (34). The case of multiple eigenvalues of M is technically more complicated. In order to compute M n we can, for instance, transform M to the Jordan canonical form (see, for instance, [10]). Here we suggest a method which is very efficient for 2 × 2 matrices. By the Cayley-Hamilton theorem ( [10]) any matrix satisfies its characteristic equation. In the case of (47) it means that M 2 = 2AλM + B. In the case of the double root (B = −A 2 ) one can easily prove by induction M n = (1 − n)A n M + nA n−1 . (56) Substituting it to (48) we get immediately (36).

B Numerical methods for ordinary differential equations
In this short note we give basic informations about some numerical methods for ordinary differential equations and we appply them to case of harmonic oscillator equation (7).
A system of linear ordinary differential equations (of any order) can always be rewritten as a single matrix equation of the first order: where the unknown y is a vector and S is a given matrix (in general tdependent). Numerical methods are almost always (see [4]) constructed for a large class of ordinary differential equations (including nonlinear ones): We denote by y n a numerical approximant to the exact solution y(t n ).
Euler's method In this case the discretization ofẍ + x = 0 is given by (4).