Journal of Applied Mathematics

Volume 2013 (2013), Article ID 518194, 15 pages

http://dx.doi.org/10.1155/2013/518194

Research Article

Symbolic Solution to Complete Ordinary Differential Equations with Constant Coefficients

Juan F. Navarro and Antonio Pérez-Carrió

Department of Mathematics, University of Alicante, Carretera San Vicente del Raspeig s/n, 03690 San Vicente del Raspeig, Alicante, Spain

Received 7 May 2013; Accepted 17 July 2013

Academic Editor: Debasish Roy

Copyright © 2013 Juan F. Navarro and Antonio Pérez-Carrió. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 .

Figure 1: 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 .

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 .

Figure 2: The first window shown by the symbolic processor allows us to introduce the parameters which define the differential equation to be solved, including the initial conditions.

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

Figure 3: Matrix .

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.

Figure 4: Matrix .

Figure 5: Matrix .

Figure 6: Matrix .

Figure 7: Matrix .

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.

Figure 8: Parameters related to the accuracy of the solution to the differential equation.

Figure 9: Matrix .

Figure 10: Comparison of real solution (black) and symbolic solutions for (light blue line) and (yellow line).

Figure 11: Comparison of real solution (black) and symbolic solutions for (dark blue line) and (green line).

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 .

Table 1: Numerical evolution of the symbolic approximation , for .

Step

1
0.1
0.0951665336810

2
0.2
0.1813305224920

3
0.3
0.2594801299955

4
0.4
0.3305856479131

5
0.5
0.3955920028686

6
0.6
0.4554122210115

7
0.7
0.5109218189711

8
0.8
0.5629540836788

9
0.9
0.6122961986812

10
1.0
0.6596861707192

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.

Table 2: Numerical error in the symbolic () and numeric () approximations, with respect to the exact solution .

0.1
0.0000047853239
0.0000001672031

0.4
0.0000712453497
0.0000358751267

0.7
0.0002074165163
0.0003544605868

1.0
0.0003710415717
0.0014307873617

1.3
0.0004547161744
0.0038253117836

1.6
0.0003284030320
0.0080436899562

1.9
0.0000622406880
0.0144197157174

2.2
0.0006216534547
0.0230329551851

2.5
0.0011277539714
0.0336675456541

2.8
0.0013572764337
0.0458134712942

3.1
0.0012347982225
0.0587076431208

3.4
0.0008953272672
0.0714078077594

3.7
0.0006002231245
0.0828896894497

4.0
0.0005492334199
0.0921556872924

4.3
0.0007154938930
0.0983431460014

4.6
0.0008300921248
0.1008204330678

4.9
0.0005487250652
0.0992608551151

5.2
0.0002963137375
0.0936868279861

5.5
0.0015298326576
0.0844798593631

5.8
0.0026838757810
0.0723553936143

6.1
0.0032571351779
0.0583046260627

6.4
0.0030289419613
0.0435091879952

6.7
0.0022189180928
0.0292364988851

7.0
0.0013691340100
0.0167257096706

7.3
0.0010043530807
0.0070750555622

7.6
0.0012821439144
0.0011409176473

7.9
0.0018689726696
0.0005417339524

8.2
0.0021436005756
0.0021897330323

8.5
0.0016148602301
0.0091082291225

8.8
0.0002869277587
0.0196147078451

9.1
0.0012796703448
0.0327897203523

9.4
0.0022771057968
0.0474740017777

9.7
0.0021518593816
0.0623710756305

10.0
0.0009803998394
0.0761625070468

Figure 12: Exact solution (black), symbolic solution (blue), and numeric solution (green).

Figure 13: Error in the symbolic (black) and numeric (blue) approximations.

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.

References

H. Poincare, Les Methodes Nouvelles de la Mecanique Celeste I, Gauthiers-Villars, Paris, France, 1892.

J. Henrard, “A survey of Poisson series processors,” Celestial Mechanics, vol. 45, no. 1–3, pp. 245–253, 1989.

A. Deprit, J. Henrard, and A. Rom, “La theorie de la lune de Delaunay et son prolongement,” Comptes Rendus de l'Académie des Sciences, vol. 271, pp. 519–520, 1970.

A. Deprit, J. Henrard, and A. Rom, “Analytical Lunar Ephemeris: Delaunay's theory,” The Astronomical Journal, vol. 76, pp. 269–272, 1971.

J. Henrard, “A new solution to the Main Problem of Lunar Theory,” Celestial Mechanics, vol. 19, no. 4, pp. 337–355, 1979.

M. Chapront-Touze, “La solution ELP du probleme central de la Lune,” Astronomy & Astrophysics, vol. 83, pp. 86–94, 1980.

M. Chapront-Touze, “Progress in the analytical theories for the orbital motion of the Moon,” Celestial Mechanics, vol. 26, no. 1, pp. 53–62, 1982.

J. F. Navarro and J. M. Ferrándiz, “A new symbolic processor for the Earth rotation theory,” Celestial Mechanics & Dynamical Astronomy, vol. 82, no. 3, pp. 243–263, 2002.

J. F. Navarro, “On the implementation of the Poincaré-Lindstedt technique,” Applied Mathematics and Computation, vol. 195, no. 1, pp. 183–189, 2008.

J. F. Navarro, “Computation of periodic solutions in perturbed second-order ODEs,” Applied Mathematics and Computation, vol. 202, no. 1, pp. 171–177, 2008.

J. F. Navarro, “On the symbolic computation of the solution to a differential equation,” Advances and Applications in Mathematical Sciences, vol. 1, no. 1, pp. 1–21, 2009.

J. F. Navarro and A. Pérez, “Principal matrix of a linear system symbolically computed,” in Proceedings of the International Conference on Numerical Analysis and Applied Mathematics (ICNAAM '08), pp. 400–402, September 2008.

C. Moler and C. van Loan, “Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later,” SIAM Review, vol. 45, no. 1, pp. 3–49, 2003.

J. F. Navarro and A. Pérez, “Symbolic computation of the solution to an homogeneous ODE with constant coefficients,” in Proceedings of the International Conference on Numerical Analysis and Applied Mathematics (ICNAAM '09), pp. 400–402, 2009.