Abstract

The aim of this paper is to introduce a symbolic technique for the computation of the solution to a complete ordinary differential equation with constant coefficients. The symbolic solution is computed via the variation of parameters method and, thus, constructed over the exponential matrix of the linear system associated with the homogeneous equation. This matrix is also symbolically determined. The accuracy of the symbolic solution is tested by comparing it with the exact solution of a test problem.

1. Introduction

Perturbation theories for differential equations containing a small parameter are quite old. The small perturbation theory originated by Sir Isaac Newton has been highly developed by many others, and an extension of this theory to the asymptotic expansion, consisting of a power series expansion in the small parameter, was devised by Poincaré [1]. The main point is that, for the most of the differential equations, it is not possible to obtain an exact solution. In cases where equations contain a small parameter, we can consider it as a perturbation parameter to obtain an asymptotic expansion of the solution. In practice, the work involved in the application of this approach to compute the solution to a differential equation cannot be performed by hand, and algebraic processors result in being a very useful tool.

As explained in Henrard [2], the first symbolic processors were developed to work with Poisson series, that is, multivariate Fourier series whose coefficients are multivariate Laurent series, where , , and and are called polynomial and angular variables, respectively. These processors were applied to problems in nonlinear mechanics or nonlinear differential equation problems, in the field of celestial mechanics. One of the first applications of these processor was concerned with the theory of the Moon. Delaunay invented his perturbation method to treat it and spent 20 years doing algebraic manipulations by hand to apply it to the problem. Deprit et al. [3, 4] prolongated the solution of Delaunay's work with the help of a special purpose symbolic processor, and Henrard [5] pushed to order . This solution was improved by iteration by Chapront-Touze [6], and planetary perturbations were also introduced by Chapront-Touze [7]. At present, the most complete solution, Ephemeride Lunaire Parisien (ELP) contains more than 50 000 periodic terms. But the motion of the Moon is not the only application of algebraic processors. There are many problems where the facilities provided by Poisson series processors can lead rather quickly to very accurate results. As examples, we would like to mention planetary theories, the theory of the rotation of the Earth (see, e.g. Navarro and Ferrándiz [8]), and artificial satellite theories.

In order to achieve better accuracies in the applications of analytical theories, high orders of the approximate solution must be computed, making necessary a continuous maintenance and revision of the existing symbolic manipulation systems, as well as the development of new packages adapted to the peculiarities of the problem to be treated. Recently, Navarro [9, 10] developed a symbolic processor to deal with the solution to second order differential equations of the form with initial conditions where , , , and , is a quasipolynomial, that is, an object of the form where , , , , and , and admits the expansion To face this problem, the algebraic processor handles objects called quasipolynomials, and it has resulted in being a useful tool in the computation of the solution to (2) by the application of the asymptotic expansion method. A modification of this processor has been employed to compute periodic solutions in equations of type (2) via the Poincaré-Lindstedt method in Navarro [9, 10]. To that end, the idea is to expand both the solution and the modified frequency with respect to the small parameter, allowing to kill secular terms which appear in the recursive scheme. The elimination of secular terms is performed through a manipulation system which works with modified quasipolynomials, that is, quasipolynomials containing undetermined constants: where , , , , and , , and are real constants with unknown value.

One year later, Navarro [11] presented a symbolic computation package based on the object-oriented philosophy for manipulating matrices whose elements lie on the set of quasipolynomials of type (4). The kernel of the symbolic processor was developed in C++, defining a class for this new object as well as a set of functions that operate on the data structure: addition, substraction, differentiation, and integration with respect to , substitution of an undetermined coefficient by a series, and many others. The goal of this processor is to provide a tool to solve a perturbed -order differential equation of the class with initial conditions where is a small real parameter, , is a quasipolynomial, and is such that with , , and .

As a first step, Navarro and Pérez [12] have developed a symbolic technique for the computation of the principal matrix of the linear system associated with a homogeneous ordinary differential equation with constant coefficients of the form where , with initial conditions being . This method provides a final analytical solution which can be completely computed in a symbolic way.

The second step in the construction of the solution to (7) is to obtain a symbolic procedure for computing the solution to the nonhomogeneus equation where, as above, is a quasipolynomial.

This new symbolic package can be useful for research and educational purposes in the field of differential equations or dynamical systems. Nowadays there are many open problems which requires massive symbolic computation as, to cite one example, the analytical theory of the resonant motion of Mercury. We would like to stress that the aim of this work is to develop a symbolic tool, not a numeric method, as numeric solutions cannot be used in perturbation methods for differential equations. However, we have performed comparisons between the symbolic solution with a numeric one, just to show the efficiency of the technique.

In next section, we summarize the scheme proposed in [12] to calculate the solution to (10), which is needed to express the solution to the complete problem (12). The technique to construct this solution will be described in Section 3.

2. Solution to the Homogeneous Problem

As mentioned, Navarro and Pérez [12] have proposed a symbolic method to compute the solution to (10) where , with initial conditions being . This method provides a final analytical solution which can be completely computed in a symbolic way.

2.1. Description of the Method

With the aid of the substitutions (10) is transformed into the system of differential equations given by where is the companion matrix, To compute the exponential of , the matrix is splitted into , where Then, we use the approximation with . This approach to obtain is of potential interest when the exponentials of matrices and can be efficiently computed and requires to be nonequal to zero. In our case, being, for any , and, for any , As it is established in Moler and van Loan [13], for a general splitting , can be determine from Thus, with being the companion matrix, we get that where . Here, for any , and , for any matrix

2.2. Adaptation to a Symbolic Formalism

Navarro and Pérez [12, 14] have proposed an adaptation of the method described above in order to compute the matrix instead of . If we do so, we obtain the principal matrix of (16), whose elements lie on the set of quasipolynomials, and the symbolic processor results in being suitable to work with those matrices. The approach is to compute taking into account that being, for any , and, for any ,

Thus, is a matrix of quasipolynomials that can be completely computed through the symbolic processor developed by Navarro [11]. The procedure for the computation of the exponential matrix is substantially simplified by using the following equation, which avoids the symbolic multiplication of and : where This relation can be obtained directly from the multiplication of matrices and as given in (28).

2.3. Symbolic Expansion of the Exponential Matrix

As it has been stated in Section 2.2, the exponential matrix can be symbolically calculated through (27) and (31) This expression can be calculated symbolically. Let us express as with Thus, where is the so-called noncommutative parenthesis, defined as where and are two noncommutative matrices, , and is determined by the following properties:(1) for all ;(2) for all , with being the entire part of ; (3) for all ;(4)if and , then for all and .

In the following, we summarize some expressions which simplify the way in which the matrix is symbolically computed. First, let us introduce the following matrices:

Lemma 1. For any such that , with for each , where is given by (29).

Lemma 2. For any such that , where and are given by (39).

Lemma 3. For any such that , where and are given by (39), and

Lemma 4. For any such that , where and are given by (39), for any , and is given by (29).

For the sake of simplicity, we have omitted the dependence on and of .

3. Solution to the Nonhomogeneous Problem

The general solution to a nonhomogeneous linear differential equation of order can be expressed as the sum of the general solution to the corresponding homogeneous, linear differential equation and any solution to the complete equation. The symbolic manipulation system calculates the solution to a nonperturbed differential equation with initial conditions of the form (12) with initial conditions where , , and being , and . With the aid of the substitutions (12) is reduced to the system of differential equations given by where

The computation of the solution to the constant coefficients linear part requires the calculation of the exponential of the matrix , :

3.1. Symbolic Expansion of the Solution

In the following, we give a formula for the symbolic expansion of the solution to the complete ordinary differential equation. To that end, let us express the noncommutative parenthesis as where each is a polynomial of degree in the indeterminate with coefficients from , and .

Let us also express being Thus, taking into account (33), the exponential matrix of can be arranged as where

In order to develop a symbolic expression of the solution to the complete differential equation, let us express as with , , and and being real polynomials of degrees and , respectively, in the indeterminate .

The product can be arranged as follows: Here, the integrand of each element of the above matrix can be written in the form of a quasipolynomial. Hence, we arrive to the formulae where , , , and These functions are computed recursively through and, for any ,

Thus, (53) can be expressed via a quasipolynomial and, thus, obtained through the designed symbolic system.

4. Symbolic and Numeric Results

4.1. On the Form of the Exponential Matrix

The aim of this section is to illustrate the form that the approximation of the exponential matrix adopts. For that purpose, let us consider the differential equation given by with initial conditions This equation is transformed into the system of differential equations given by where is the companion matrix According to the technique described in Section 2, we split matrix into , where Equation (36), allows us to compute a symbolic approximation to the solution to the differential equation. In the problem we are discussing, , , , and can be written as with

Taking, for instance, , we get that For the sake of simplicity, we have introduced the following notation:

Matrices , , , , , , , and are computed with the help of the designed symbolic processor, with the result where where where where where where So, the noncommutative parentheses are given by where being

In Figure 1, we show a comparison of the solution computed through the symbolic method presented in this paper, and the numeric solution to the problem calculated by a Runge-Kutta fourth order method with a step of . Let us stress that the difference between both solutions at any time is smaller than .

4.2. Description of the Program

In order to describe the algebraic processor, let us introduce the following test problem: with initial conditions This equation describes the (small) angular position in radians of a forced damped pendulum with a periodic driving force. In this equation, represents the inertia, represents the friction, and represents a sinusoidal driving torque applied at the pivot of the pendulum. The exact solution to this problem is given by The program proceeds as follows. The first window (Figure 2) allows us introducing the definition of variables, the order of the ODE and its coefficients, and the function of the nonhomogeneous term .

The next window visualizes the expression of the companion matrix related to the ODE (Figure 3).

Then, matrices , , , and , which are used in the calculation of the exponential matrix, are computed by using a value of satisfying (23). These matrices are presented in Figures 4, 5, 6, and 7.

In Figures 8 and 9, we evaluate the exponential matrix in and for . The algebraic processor enables us to change these two parameters. We should emphasize that the symbolic solution to the problem is dependent on the value of parameter . We have computed this solution for different values of this parameter (, and ). In Figures 10 and 11, we compare the real solution with symbolic solutions for these values of . We observe that the symbolic solution gets closer to the real solution as the value of is increased.

Table 1 shows the numerical evolution of the symbolic approximation , computed for . These values have been obtained by substituting variable in the symbolic solution by values .

In Figure 12, we compare the symbolic approximation, computed via the symbolic processor, with the exact solution (87). We also compare these two functions with a numeric solution computed through a Runge-Kutta 4th order method with step . We see that the solution symbolically computed fits better to the exact solution to than the described numerical solution. In Table 2, we show the error on the solution computed using both methods. We observe that error obtained in the symbolic approximation is quite smaller that the numeric one. This fact can be better observed in Figure 13, where we show the difference in absolute value between the exact solution and the symbolic and numeric solutions. Nevertheless, it is not the goal of this paper to develop a numeric tool. The symbolic technique we have developed provides an analytical solution which can be used as a kernel to apply perturbation methods to compute the solution to a perturbed differential equation depending on a small parameter.

5. Conclusions

We have developed a symbolic processor as well as a symbolic technique in order to deal with the solution to a linear ordinary differential equation with constant coefficients of order . The solution to this equation is computed completely in a symbolic way that is, the solution is expressed as a function of . This function can be of great interest to face a perturbed differential equation by means of a perturbation method. The aim of this work is not to develop a numeric method. However, we have compared the symbolic solution obtained via our symbolic processor with a numeric solution computed with a Runge-Kutta method, in order to show how our solution is closer to the exact solution.