Research Article | Open Access
The Asymptotic Expansion Method via Symbolic Computation
This paper describes an algorithm for implementing a perturbation method based on an asymptotic expansion of the solution to a second-order differential equation. We also introduce a new symbolic computation system which works with the so-called modified quasipolynomials, as well as an implementation of the algorithm on it.
The origin of symbolic manipulation derives from the sheer magnitude of the work involved in the building of perturbation theories, which made inevitable that scientific community became interested in the possibility of constructing those theories with the help of computers.
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é (1892) . 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 seem to be a very useful tool.
As explained in , 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 equations problems, in the field of Celestial Mechanics.
One of the first applications of these processors 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  pushed it to order 25. This solution was improved by iteration by Chapront-Touzé , and planetary perturbations were also introduced by Chapront-Touzé . At present, the most complete solution, ELP (Ephemeride Lunaire Parisien) contains more than 50000 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 (e.g., ), and artificial satellite theories (AST). Abad et al.  have analysed the convenience of developing specific computer algebra systems in order to deal with AST. As it is explained in detail in their work, the series involved in the computation of the solution through the application of a Lie transformation have a total amount of almost 2000000 terms.
Nowadays, many general purpose computer algebra packages—as, for example, Maple, Mathematica, and Matlab—contain tools for the calculation of the solution of certain classes of ODEs. All these packages have the advantage of being very general, so they can deal with a lot of problems of different nature. However, if one is interested in higher order solutions, the most common perturbation methods tend to produce expressions containing thousands of terms, and their treatment with those general processors becomes inefficient.
To achieve better accuracies in the applications of analytical theories, high orders of the approximate solution must be computed, making a continuous maintenance and revision of the existing symbolic manipulation systems necessary, as well as the development of new packages adapted to the peculiarities of the problem to be treated.
In order to contribute to the solution of this problem, we have developed a symbolic computation package (called MQSP) based on the object-oriented philosophy which manipulates objects of the form where , , and are real constants with unknown value. We will refer to those elements as modified quasipolynomials . The kernel of this symbolic processor has been developed in C++. The operations on series implemented in the manipulator are the usual operations of the Algebra of quasipolynomials: addition, subtraction, multiplication, multiplication by a scalar, differentiation with respect to , and substitution of a quasipolynomial into an undetermined coefficient.
We have also constructed a set of subroutines to deal with the solution to the perturbed differential equation
where , is given by (1.2), and
In a previous contribution, the author employed the kernel of this processor to compute periodic solutions in equations of type (1.3) via the Poincaré-Lindstedt method [11, 12]. If the unperturbed equation () has periodic solutions and is a measure of the size of the perturbing terms, then the trajectories for the full system will remain pretty close to those of the nonperturbed system, for any finite period of time () with an error not larger than . In general, even a small perturbation is enough to destroy periodicity, that is, nonlinearity will end with most of the periodic orbits of the unperturbed system, but some of them may persist. To calculate those periodic orbits, the solution and the modified frequency are expanded with respect to the small parameter, allowing to kill secular terms which appear in the recursive scheme.
The aim of this paper is to construct an algorithm to implement the asymptotic expansion method . This new implementation is general and does not depend on the function as given in (1.4), that is, the user does not need to programme the algorithm described in this paper, and only has to introduce the adequate parameters when calling the corresponding routine of the package. The code of this specific symbolic system is not available on internet but it can be provided by contacting its author.
2. Data Structure
The algorithms that can be implemented to perform the basic manipulation on a series and their efficiency depend on the way a series is coded. An overcoded structure that makes good use of memory generally requires complex algorithms, which increase the computational cost in terms of time. On the other hand, an undercoded computational representation of the terms generates simple algorithms, because the location of all the coefficients can be obtained directly. However, this scheme presents the inconvenience of being very wasteful in the memory resources required for the storage of the series .
In this section, we follow San-Juan and Abad  to introduce the representation of a mathematical object in a computer. To that purpose, let us introduce the concepts of normal and canonical functions. Let be a set of symbolic objects, and let ~ be an equivalence relation in , defined as follows: if , with . Here, the operator is considered as the equality on the mathematical object level. Moreover, if and are identical as symbolic objects. A function is said to be normal in if for all , and is said to be canonical in if it is normal and for all . Thus, a canonical function provides identical objects when objects are equivalent, that is, when they represent the same mathematical object.
For the sake of simplicity, let us focus our attention on the set of quasipolynomials in the independent variable . A quasipolynomial is a map , defined by where , , and . Let us now consider the set of quasipolynomials in the independent variable , .
2.1. QuasiPolynomials as Symbolic Objects
We look for a canonical representation for each equivalence class defined in . For that purpose, the following operations must be performed over each quasipolynomial. Let us consider a quasipolynomial . If , the following rules must be applied: The terms of a quasipolynomial will be ordered as follows: let us consider two term of a quasipolynomial: and . We say that if () or ( and ) or (, , and ) or (, , , and ) or (, , , , and ). The terms of a quasipolynomial
with identical values of , and must be grouped together.
2.2. QuasiPolynomials as Computational Objects
In this section, we will consider the basic information which characterizes a quasipolynomial, as well as the data structure to store it in the computer. This must be done preserving the canonical representation we have chosen.
Each quasipolynomial is collected in a sorted dynamical list: a sorted list is one in which the order of the items is defined by some collating sequence. The codification of each term of the list contained in a quasipolynomial is statical, and given by the following elements. are the coefficients of the term. is the degree of the monomial . is the exponent of the exponential part of the term. is the frequency of the periodic part.
In Scheme 1 we represent the way a quasipolynomial is stored in memory by using a simple list.
As pointed out in , most of the operations involving a series are based on navigating and searching through the structure that represents the series. For example, the addition of two quasipolynomials is equivalent to inserting each term of one series into the other one. So, a good choice of the data structure results in simple and efficient algorithms. The binary tree resulting seems to be a very useful data structure for rapidly storing sorted data and rapidly retrieving saved data. A binary tree is composed of parent nodes, or leaves, each of which stores data and also links to up to two other child nodes (leaves), which can be visualized spatially as below the first node with one placed to the left and the other to the right. In this structure, the relationship between the linked leaves and the linking leaf makes the binary tree an efficient data structure: the leaf on the left has a lesser key value, and the leaf on the right has an equal or greater key value.
A special type of tree is the red-black tree. In a red-black tree, each node has a color attribute, the value of which is either red or black. In addition to the ordinary requirements imposed on binary search trees, the following additional requirements of any valid red-black tree apply: A node is either red or black. The root is black. All leaves are black, even when the parent is black. Both children of every red node are black. Every simple path from a node to a descendant leaf contains the same number of black nodes. A critical property of red-black trees is enforced by these constraints: the longest path from the root to a leaf is no more than twice as long as the shortest path from the root to a leaf in that tree. The result is that the tree is roughly balanced. Since operations such as inserting, deleting, and finding values requires worst-case time proportional to the height of the tree, this fact makes the red-black tree efficient, for instance, the search-time results to be .
With the use of this structure, the complexity of the algorithms for addition, multiplication, derivation, and integration of quasipolynomials is significantly reduced. Unfortunately, this structure, which results to be ideal for Poisson series, cannot be applied directly in our case due to the fact that the numbers which identify each term of a quasipolynomial are not indexed arrays. However, an alternative aggrupation of terms in a quasipolynomial can be performed in order to introduce this balanced structure. To do that, let us express a quasipolynomial as follows:
where and are polynomials in the variable with constant coefficients, of degree and respectively, being : with . If we aggrupate terms of a quasipolynomial in such a manner, we can use a tree structure to store it, saving only significant terms. In Scheme 2 we show the tree structure in which the quasipolynomial is stored.
In Scheme 2 we show how we store in memory the quasipolynomial where
with and being sets of indices for any , and for any , and for any , and for any , and and for any . Moreover, and for any , 2 and belonging to the corresponding set of indices.
To illustrate the way a quasipolynomial is stored by means of the tree structure described above, we show in Scheme 3 the structure associated to the quasipolynomial given by
2.3. Modified QuasiPolynomials as Computational Objects
In our symbolic system, we represent a modified quasipolynomial by an ordered and dynamical list of terms, keeping in memory only significant terms. The codification of a term is statical, and given by the following elements.(1) are the coefficients of the term (for and , resp.).(2). For each , represents the exponent of the undetermined coefficient .(3) is the degree of the monomial .(4) is the exponent of the exponential part of the term.(5) is the frequency of the periodic part.
We have included the two following additional vectors, which are common to all the terms of a modified quasipolynomial:. For each , indicates that is a coefficient with unknown value, while implies that has a real value assigned,. If the coefficient has a real value assigned, its content is given by , for each .
It is also absolutely essential to store the number of undetermined coefficients that are currently in use, to generate correctly new undetermined constants if needed. In Table 1 we illustrate the way a term is coded with a few examples, where it has been assumed that .
3. Design of the Symbolic System in C++
As pointed out in , object-oriented programming using C++ provides many advantages in the design of computer algebra systems, as this programming technique combines both the data and the functions that operate on that data into a single unit (called class). The main reasons given by Hardy et al. to use C++ to implement a symbolic system are as follows.C++ allows the introduction of abstract data types. Thus, we can define a modified quasipolynomial as an abstract data type. The language C++ supports encapsulation, inheritance, polymorphism, and operator overloading. Consequently, we can overload the operators +, −, and for modified quasipolynomials, as well as andfor multiplication and division of modified quasipolynomials by real numbers.
Some symbolic computation systems have been constructed in C++. MuPAD is a computer algebra system developed by the MuPAD research group at the University of Paderborn, in Germany. This symbolic system manipulates formulas symbolically and provides packages for linear algebra, differential equations, number theory, statistics, and functional programming, as well as an interactive graphic system that supports animations and transparent areas in 3D. MuPAD also offers a programming language that supports object-oriented programming and functional programming .
Symbolic C++ also uses C++ to develop a computer algebra system. This package introduces the Symbolic class which is used for all symbolic computation, and provides almost all of the features required for symbolic computation including symbolic terms, substitution, noncommutative multiplication, and vectors and matrices.
Both symbolic systems could have been used to implement the algorithm over a general symbolic class. However, as the goal of the symbolic system we are developing is to handle modified quasipolynomial to apply perturbation methods to solve some types of differential equations, we have constructed it directly over C++, instead of using some other symbolic processor.
The specific symbolic processor designed is written in clean C++ code, is very portable, it can compile stand-alone, and is easily embeddable. It implements a new data type, called MQ, which represents a series of the form given by (1.2). The class MQ is defined as an ordered and dynamical list of terms. To that end, we have also defined a class associated to the structure of a term of a modified quasipolynomial, called Term. The definition of these two classes in C++ code is shown in Algorithms 1 and 2, respectively.
The set of routines developed includes: addition, subtraction and product of series, multiplication of a series by a-real number, differentiation with respect to , and computation of the solution to a linear second-order differential equation of type where is a modified quasipolynomial presenting undetermined coefficients, and computation of the solution to (1.3) via the asymptotic expansion method. These two algorithms are described in detail in Sections 4 and 5, respectively.
The basic algebra associated to quasipolynomials is easily implemented because of the undercoded data scheme chosen for their computational representation. Thus, for example, the addition of two quasipolynomials is performed by directly juxtaposing both quasipolynomials, arranging the resulting series, and joining terms with equal elements. In a similar way, the rest of algebraic operations are simply implemented.
4. Solution of a Linear Second-Order ODE
The general solution to a nonhomogeneous differential equation can be expressed as the sum of general solutions to the corresponding homogenous, linear differential equation and any solution to the complete equation . The solution to the homogeneous ODE is expressed in terms of the roots of the characteristic equation, , and it is well-known. We will resume now the formulae that are required to construct a particular solution to a complete ordinary differential equation of second-order (3.1). Without loss of generality, we will assume that is written as follows: where , and , and are th and th degree polynomials in with undetermined coefficients respectively, of the form being
with , , , for , and for .
The principle of superposition is applied to calculate the particular solution, so we can focus our attention on the computation of a particular solution to the equation where being with , , , for , and for .
At this point, we will distinguish two cases depending on if or .
Subcase 1.1. If , the particular solution to the complete differential equation is expressed as Thus, substituting and its derivatives in (4.7), we get that and, for any ,
Subcase 1.2. If and (i.e., ), the particular solution can be written as Now, the substitution of and its derivatives into (4.7) leads to the following formula for , and, for any ,
Subcase 1.3. If and (i.e., ), the particular solution can be written as Now, the substitution of and its derivatives into (4.7) leads to the following formula for ,
Case 2. Let us consider the second-order differential equation (4.4), where now . There are two possible situations:is not a root of the characteristic equation, that is, or , is a root of the characteristic equation, that is, and .
Subcase 2.1. Let us call . Then, Now, if we assume and to be and substitute and its derivatives into (4.4), we get that where . Note that , because or . Now, we can compute and by solving the system As before, this system can be solved by applying the Cramer’s rule. Finally, for any , we have to solve the system in and , as follows:
Subcase 2.2. Now we consider the case where and . This implies that and . In this case, the particular solution adopts the form and substituting , , and into (4.4), we obtain for any .
5. Computation of the Solution to the Perturbed Problem
The symbolic package is thought to compute the solution to (1.3). The standard approach  is to try a power series solution of the formThis series is inserted into the governing equation and initial conditions, and coefficients of same powers of are then grouped to obtain a collection of equations for the coefficient functions , which are then solved in a sequential manner. The resulting series need not converge for any value of ; nevertheless, the solution can be useful in approximating the function when is small.
Considering the zero-order term in yields the so-called nonperturbed problem. The symbolic manipulation system calculates the solution to a differential equation of the form (5.2), and arranges it as a quasipolynomial, as it has been described in detail in the previous section.
The coefficient of the solution to the order is computed by solving the equation where the notation refers to the th order term of the series .
At each order of the solution, the series , , and must be computed once the order has been solved, for each , following the formulae
According to this, the algorithm to apply the asymptotic expansion method to solve the initial value problem given by (1.3), consists of the following steps.Define a three-dimensional array of quasipolynomials , , , where , , , and , being the order of the asymptotic expansion.Define an array of quasipolynomials , where .Initialize , and proceed as follows: Compute as the solution to (5.2). Compute . Calculate, for each such that , For each , such that , , determine the quasipolynomial For each such that , do the following. Compute the quasipolynomial Calculate as the solution to (5.3), Compute . Calculate, for each such that , For each such that , compute For each such that ,
The input arguments of the algorithm consist of the order of the asymptotic approximation, the coefficients , of the differential equation, the real values , , and which define the initial conditions of the problem, the quasipolynomial , and the perturbed part of the equation . As can be written as given in (1.4),
it can be specified by a real matrix,
The parameter is also required as input.
The output argument of the algorithm is the array containing the coefficients of the asymptotic expansion. The asymptotic approximation to the order is given by
6. On the Implementation of the Algorithm
In Algorithms 3, 4, 5, 6, and 7 we show the C++ code for the implementation of the algorithm as described in Section 5, omitting error control sentences for the sake of simplicity. The head of the definition of the function that implements the algorithm is shown below.
MQ *solveAM (double a1, double a0, double t0, double x0, double dx0, MQ u, double **f, int m, int order)
Let us consider the initial value problem given by where . Here, , , , , , , and is defined by the following matrix:
In Table 2 we show the output generated from the symbolic manipulator, just to order nine for the sake of simplicity, as the number of terms of the asymptotic expansion increases very quickly with the order of the expansion.