The Asymptotic Expansion Method via Symbolic Computation

Abstract

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.

1. Introduction

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) [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 seem to be a very useful tool.

As explained in [2], the first symbolic processors were developed to work with Poisson series, that is, multivariate Fourier series whose coefficients are multivariate Laurent series,_{𝑖1,…,𝑖𝑛𝑗1,…,𝑗𝑚}𝐶^{𝑗1,…,𝑗𝑚}_{𝑖1,…,𝑖𝑛}𝑥^𝑖1₁𝑥^𝑖2₂⋯𝑥^𝑖𝑛_𝑛cos𝑗sin¹𝜙₁+⋯+𝑗_𝑚𝜙_𝑚, (1.1) where 𝐶^{𝑗1,…,𝑗𝑚}_{𝑖1,…,𝑖𝑛}∈ℝ, 𝑖₁,…,𝑖_𝑛,𝑗₁,…,𝑗_𝑚∈ℤ, 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 [5] pushed it to order 25. This solution was improved by iteration by Chapront-Touzé [6], and planetary perturbations were also introduced by Chapront-Touzé [7]. 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., [8]), and artificial satellite theories (AST). Abad et al. [9] 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_𝜈∑_≥0𝜏^𝜎1₁×⋯×𝜏^𝜎𝑄_𝑄𝑡^𝑛𝜈𝑒^𝛼𝜈𝑡𝜆_𝜈cos𝜔_𝜈𝑡+𝜇_𝜈sin𝜔_𝜈𝑡, (1.2) where 𝑛_𝜈∈ℕ,𝛼_𝜈,𝜔_𝜈,𝜆_𝜈,𝜇_𝜈∈ℝ, 𝜎_𝜈∈ℤ, and 𝜏₁,…,𝜏_𝑄 are real constants with unknown value. We will refer to those elements as modified quasipolynomials [10]. 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𝑥̈+𝑎₁𝑥̇+𝑎₀𝑥=𝑢(𝑡)+𝜖𝑓(𝑥,𝑥̇),𝑥𝑡₀=𝑥₀,𝑥̇𝑡₀=𝑥̇₀, (1.3)

where 𝑎₀,𝑎₁,𝑡₀,𝑥₀,𝑥̇₀∈ℝ, 𝑢(𝑡) is given by (1.2), and𝑓(𝑥,𝑥̇)=^𝑀_{𝜅=00≤𝜈≤𝜅}𝑓_{𝜈,𝜅−𝜈}𝑥^𝜈𝑥̇^𝜅−𝜈,𝑓_{𝜈,𝜅−𝜈}∈ℝ. (1.4)

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 (𝜖=0) 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 𝑡₀<𝑡<𝑡₀+𝛼 (𝛼>0) 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 [13]. 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 [2].

In this section, we follow San-Juan and Abad [14] 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𝑢(𝑡)=_𝜈≥0𝑡^𝑛𝜈𝑒^𝛼𝜈𝑡𝜆_𝜈cos𝜔_𝜈𝑡+𝜇_𝜈sin𝜔_𝜈𝑡, (2.1) 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. (1)Let us consider a quasipolynomial 𝑢(𝑡)=∑_𝜈≥0𝑡^𝑛𝜈𝑒^𝛼𝜈𝑡(𝜆_𝜈cos(𝜔_𝜈𝑡)+𝜇_𝜈sin(𝜔_𝜈𝑡)). If 𝜔_𝜈<0, the following rules must be applied:sin−𝜔_𝜈𝑡=−sin𝜔_𝜈𝑡,cos−𝜔_𝜈𝑡=cos𝜔_𝜈𝑡. (2.2)(2) The terms of a quasipolynomial will be ordered as follows: let us consider two term of a quasipolynomial: 𝜏₁=𝑡^𝑛1𝑒^𝛼1𝑡(𝜆₁cos𝜔₁𝑡+𝜇₂sin𝜔₁𝑡) and 𝜏₂=𝑡^𝑛2𝑒^𝛼2𝑡(𝜆₂cos𝜔₂𝑡+𝜇₂cos𝜔₂𝑡). We say that 𝜏₁<𝜏₂ if (𝑛₁<𝑛₂) or (𝑛₁=𝑛₂ and 𝛼₁<𝛼₂) or (𝑛₁=𝑛₂, 𝛼₁=𝛼₂, and 𝜔₁<𝜔₂) or (𝑛₁=𝑛₂, 𝛼₁=𝛼₂, 𝜔₁=𝜔₂, and 𝜆₁<𝜆₂) or (𝑛₁=𝑛₂, 𝛼₁=𝛼₂, 𝜔₁=𝜔₂, 𝜆₁<𝜆₂, and 𝜇₁<𝜇₂).(3) The terms of a quasipolynomial𝑢(𝑡)=_𝜈≥0𝑡^𝑛𝜈𝑒^𝛼𝜈𝑡𝜆_𝜈cos𝜔_𝜈𝑡+𝜇_𝜈sin𝜔_𝜈𝑡 (2.3)

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.(1)𝜆,𝜇∈ℝ are the coefficients of the term.(2)𝑛∈ℕ is the degree of the monomial 𝑡.(3)𝛼∈ℝ is the exponent of the exponential part of the term.(4)𝜔≥0 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.

Scheme 1 

As pointed out in [14], 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 𝑂(log𝑛).

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:𝑢(𝑡)=_𝜈≥0𝐶_𝑝,𝜈(𝑡)𝑒^𝛼𝜈𝑡cos𝜔_𝜈𝑡+𝑆_𝑞,𝜈(𝑡)𝑒^𝛼𝜈𝑡sin𝜔_𝜈𝑡, (2.4)

where 𝐶_𝑝,𝜈(𝑡) and 𝑆_𝑞,𝜈(𝑡) are polynomials in the variable 𝑡 with constant coefficients, of degree 𝑝 and 𝑞 respectively, being 𝑝,𝑞∈ℕ:𝐶_𝑝,𝜈(𝑡)=𝜆₀+𝜆₁𝑡+𝜆₂𝑡²+⋯+𝜆_𝑝𝑡^𝑝,𝑆_𝑞,𝜈(𝑡)=𝜇₀+𝜇₁𝑡+𝜇₂𝑡²+⋯+𝜇_𝑞𝑡^𝑞, (2.5) 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.

Scheme 2 

In Scheme 2 we show how we store in memory the quasipolynomial𝑢(𝑡)=𝐶₁(𝑡)𝑒^𝛼1𝑡cos𝜔₁𝑡+𝑆₁(𝑡)𝑒^𝛼1𝑡sin𝜔₁𝑡+𝐶₂(𝑡)𝑒^𝛼2𝑡cos𝜔₂𝑡+𝑆₂(𝑡)𝑒^𝛼2𝑡sin𝜔₂𝑡, (2.6) where𝐶₁(𝑡)=_{𝑖∈𝕀𝐶,1}𝜆_1,𝑖𝑡^𝜈1,𝑖,𝑆₁(𝑡)=_{𝑖∈𝕀𝑆,1}𝜇_1,𝑖𝑡^𝜎1,𝑖,𝐶₂(𝑡)=_{𝑖∈𝕀𝐶,2}𝜆_2,𝑖𝑡^𝜈2,𝑖,𝑆₂(𝑡)=_{𝑖∈𝕀𝑆,2}𝜇_2,𝑖𝑡^𝜎2,𝑖, (2.7)

with 𝕀_𝐶,𝑗 and 𝕀_𝑆,𝑗 being sets of indices for any 𝑗=1,2, 𝜈_1,𝑖≥0 and 𝜆_1,𝑖≠0 for any 𝑖∈𝕀_𝐶,1, 𝜈_2,𝑖≥0 and 𝜆_2,𝑖≠0 for any 𝑖∈𝕀_𝐶,2, 𝜎_1,𝑖≥0 and 𝜇_1,𝑖≠0 for any 𝑖∈𝕀_𝑆,1, and 𝜎_2,𝑖≥0 and 𝜇_2,𝑖≠0 for any 𝑖∈𝕀_𝑆,2. Moreover, 𝜈_{𝑗,𝑖−1}<𝜈_𝑗,𝑖<𝜈_𝑗,𝑖+1 and 𝜎_{𝑗,𝑖−1}<𝜎_𝑗,𝑖<𝜎_𝑗,𝑖+1 for any 𝑗=1, 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𝑢(𝑡)=3𝑡+8𝑡³−𝑡⁷𝑒^2𝑡cos𝑡+𝑒^2𝑡sin𝑡𝑡−𝑡⁴+3𝑡⁷. (2.8)

Scheme 3 

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 cos and sin, resp.).(2)(𝜎₁,…,𝜎_𝑄)∈ℤ^𝑄. For each 0≤𝜈≤𝑄, 𝜎_𝜈 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)𝜔≥0 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:(1)(𝑚₁,…,𝑚_𝑄)∈{0,1}^𝑄. For each 0≤𝜈≤𝑄, 𝑚_𝜈=0 indicates that 𝜏_𝜈 is a coefficient with unknown value, while 𝑚_𝜈=1 implies that 𝜏_𝜈 has a real value assigned,(2)𝜏=(𝜏₁,…,𝜏_𝑄)∈ℝ^𝑄. If the coefficient 𝜏_𝜈 has a real value assigned, its content is given by 𝜏_𝜈, for each 0≤𝜈≤𝑄.

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 𝑄=4.

Table 1 

Some examples of representation of terms.

3. Design of the Symbolic System in C++

As pointed out in [15], 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.(1)C++ allows the introduction of abstract data types. Thus, we can define a modified quasipolynomial as an abstract data type.(2) The language C++ supports encapsulation, inheritance, polymorphism, and operator overloading. Consequently, we can overload the operators +, −, and ∗ for modified quasipolynomials, as well as ∗ and/for 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 [16].

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.

typedef Term *serie; class MQ {     public:      [definition of functions]       private:     serie first; };

Algorithm 1 

C++ code for the definition of a modified quasipolynomial.

class Term {      public:      Term ();      Term (double L, double M,         int *sigm, int m, double alph, double omeg,         Term *nxt);   ~Term ();      private:    double lambda, mu;    int *sigma, n;    double alpha, omega;    Term *next;      friend class MQ; };

Algorithm 2 

C++ code for the definition of a term of a modified quasipolynomial.

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𝑥̈+𝑎₁𝑥̇+𝑎₀𝑥=𝑢(𝑡), (3.1) 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 [17]. The solution to the homogeneous ODE is expressed in terms of the roots of the characteristic equation, 𝛼²+𝑎₁𝛼+𝑎₀=0, 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:𝑢(𝑡)=_𝜈≥0𝑒^𝛼𝜈𝑡𝑝_𝜈,𝑛(𝑡,𝜏)cos𝜔_𝜈𝑡+𝑞_𝜈,𝑚(𝑡,𝜏)sin𝜔_𝜈𝑡, (4.1) where 𝛼,𝜔∈ℝ, and 𝑝_𝜈,𝑛(𝑡,𝜏), and 𝑞_𝜈,𝑚(𝑡,𝜏) are 𝑛th and 𝑚th degree polynomials in 𝑡 with undetermined coefficients respectively, of the form𝑝_𝜈,𝑛(𝑡,𝜏)=𝑢_𝜈,0+𝑢_𝜈,1𝑡+⋯+𝑢_𝜈,𝑛𝑡^𝑛,𝑞_𝜈,𝑚(𝑡,𝜏)=𝑣_𝜈,0+𝑣_𝜈,1𝑡+⋯+𝑣_𝜈,𝑚𝑡^𝑚, (4.2) being𝑢_𝜈,𝜌=𝑢^∗_𝜈,𝜌𝜏^{𝜎𝜈,𝜌,1}₁×⋯×𝜏^{𝜎𝜈,𝜌,𝑄}_𝑄,𝑣_𝜈,𝜌=𝑣^∗_𝜈,𝜌𝜏^{𝜎∗𝜈,𝜌,1}₁×⋯×𝜏^{𝜎∗𝜈,𝜌,𝑄}_𝑄, (4.3)

with 𝑢^∗_𝜈,𝜌,𝑣^∗_𝜈,𝜌∈ℝ, 𝜎_{𝜈,𝜌,𝑖},𝜎^∗_{𝜈,𝜌,𝑖}∈ℤ, 1≤𝑖≤𝑄, 0≤𝜌≤𝑛 for 𝑢_𝜈,𝜌, and 0≤𝜌≤𝑚 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𝑥̈+𝑎₁𝑥̇+𝑎₀𝑥=𝑒^𝛼𝑡𝑝_𝑛(𝑡,𝜏)cos(𝜔𝑡)+𝑞_𝑚(𝑡,𝜏)sin(𝜔𝑡), (4.4) where𝑝_𝑛(𝑡,𝜏)=𝑢₀+𝑢₁𝑡+⋯+𝑢_𝑛𝑡^𝑛,𝑞_𝑚(𝑡,𝜏)=𝑣₀+𝑣₁𝑡+⋯+𝑣_𝑚𝑡^𝑚, (4.5) being 𝑢_𝜌=𝑢^∗_𝜌𝜏^𝜎𝜌,1₁×⋯×𝜏^{𝜎𝜌,𝑄}_𝑄,𝑣_𝜌=𝑣^∗_𝜌𝜏^{𝜎∗𝜌,1}₁×⋯×𝜏^{𝜎∗𝜌,𝑄}_𝑄, (4.6) with 𝑢^∗_𝜌,𝑣^∗_𝜌∈ℝ, 𝜎_𝜌,𝑖,𝜎^∗_𝜌,𝑖∈ℤ, 1≤𝑖≤𝑄, 0≤𝜌≤𝑛 for 𝑢_𝜌, and 0≤𝜌≤𝑚 for 𝑣_𝜌.

At this point, we will distinguish two cases depending on if 𝜔=0 or 𝜔≠0.

Case 1. Let us consider the second-order ODE (4.4), with 𝜔=0. Then, the equation is written as 𝑥̈+𝑎₁𝑥̇+𝑎₀𝑥=𝑝_𝑛(𝑡,𝜏)𝑒^𝛼𝑡, (4.7) where 𝑝_𝑛(𝑡,𝜏) is given by (4.5).

Subcase 1.1. If 𝛼²+𝑎₁𝛼+𝑎₀≠0, the particular solution to the complete differential equation is expressed as𝑥(𝑡)=𝛼₀(𝜏)+𝛼₁(𝜏)𝑡+𝛼₂(𝜏)𝑡²+⋯+𝛼_𝑛(𝜏)𝑡^𝑛𝑒^𝛼𝑡. (4.8) Thus, substituting 𝑥(𝑡) and its derivatives in (4.7), we get that𝛼_𝑛(𝜏)=𝑢_𝑛(𝜏)𝛼²+𝑎₁𝛼+𝑎₀,𝛼_𝑛−1(𝜏)=𝑢_𝑛−1(𝜏)−𝑎₁+2𝛼𝑛𝛼_𝑛(𝜏)𝛼²+𝑎₁𝛼+𝑎₀, (4.9) and, for any 𝑝<𝑛−1,𝛼_𝑝(𝜏)=𝑢_𝑝(𝜏)−𝑎₁+2𝛼𝑝+1𝛼_𝑝+1(𝜏)−𝑝+1𝑝+2𝛼_𝑝+2(𝜏)𝛼²+𝑎₁𝛼+𝑎₀. (4.10)

Subcase 1.2. If 𝛼²+𝑎₁𝛼+𝑎₀=0 and 𝑎²₁≠4𝑎₀ (i.e., 𝛼≠−𝑎₁/2), the particular solution can be written as𝑥(𝑡)=𝛼₀(𝜏)𝑡+𝛼₁(𝜏)𝑡²+𝛼₂(𝜏)𝑡³+⋯+𝛼_𝑛(𝜏)𝑡^𝑛+1𝑒^𝛼𝑡. (4.11) Now, the substitution of 𝑥(𝑡) and its derivatives into (4.7) leads to the following formula for 𝛼_𝑛(𝜏),𝛼_𝑛(𝜏)=𝑢_𝑛(𝜏)2𝛼+𝑎₁(𝑛+1), (4.12) and, for any 𝑝≤𝑛−1,𝛼_𝑝(𝜏)=𝑢_𝑝(𝜏)−𝑝+1𝑝+2𝛼_𝑝+1(𝜏)2𝛼+𝑎₁𝑝+1. (4.13)

Subcase 1.3. If 𝛼²+𝑎₁𝛼+𝑎₀=0 and 𝑎²₁=4𝑎₀ (i.e., 𝛼=−𝑎₁/2), the particular solution can be written as𝑥(𝑡)=𝛼₀(𝜏)𝑡²+𝛼₁(𝜏)𝑡³+𝛼₂(𝜏)𝑡⁴+⋯+𝛼_𝑛(𝜏)𝑡^𝑛+2𝑒^𝛼𝑡. (4.14) Now, the substitution of 𝑥(𝑡) and its derivatives into (4.7) leads to the following formula for 𝛼_𝑝(𝜏),𝛼_𝑝(𝜏)=𝑢_𝑝(𝜏)𝑝+1𝑝+2,𝑝=0,…,𝑛. (4.15)

Case 2. Let us consider the second-order differential equation (4.4), where now 𝜔≠0. There are two possible situations:

(1)𝛼±𝑖𝜔 is not a root of the characteristic equation, that is, 𝛼≠−𝑎₁/2 or 𝜔≠4𝑎₀−𝑎²₁/2,(2)𝛼±𝑖𝜔 is a root of the characteristic equation, that is, 𝛼=−𝑎₁/2 and 𝜔=4𝑎₀−𝑎²₁/2.

Subcase 2.1. Let us call 𝑁=max{𝑛,𝑚}. Then,𝑥(𝑡)=𝛼₀(𝜏)+𝛼₁(𝜏)𝑡+⋯+𝛼_𝑁(𝜏)𝑡^𝑁𝑒^𝛼𝑡()+𝜌cos𝜔𝑡₀(𝜏)+𝜌₁(𝜏)𝑡+⋯+𝜌_𝑁(𝜏)𝑡^𝑁𝑒^𝛼𝑡sin(𝜔𝑡). (4.16) Now, if we assume 𝑝_𝑁(𝑡,𝜏) and 𝑞_𝑁(𝑡,𝜏) to be𝑝_𝑁(𝑡,𝜏)=𝑢₀(𝜏)+𝑢₁(𝜏)𝑡+⋯+𝑢_𝑁(𝜏)𝑡^𝑁,𝑞_𝑁(𝑡,𝜏)=𝑣₀(𝜏)+𝑣₁(𝜏)𝑡+⋯+𝑣_𝑁(𝜏)𝑡^𝑁, (4.17) and substitute 𝑥(𝑡) and its derivatives into (4.4), we get that𝛼_𝑁(𝜏)=1Δ𝑢_𝑁(𝜏)𝜔𝑎₁𝑣+2𝛼_𝑁(𝜏)𝑎₀+𝑎₁𝛼+𝛼²−𝜔²,𝜌_𝑁(𝜏)=1Δ𝑎₀+𝑎₁𝛼+𝛼²−𝜔²𝑢_𝑁(𝜏)−𝜔𝑎₁+2𝛼𝑣_𝑁(𝜏), (4.18) where Δ=(𝑎₀+𝑎₁𝛼+𝛼²−𝜔²)²+𝜔²(𝑎₁+2𝛼)². Note that Δ≠0, because 𝑎₀+𝑎₁𝛼+𝛼²−𝜔²≠0 or 𝜔(2𝛼+𝑎₁)≠0. Now, we can compute 𝛼_𝑁−1(𝜏) and 𝜌_𝑁−1(𝜏) by solving the system𝑢_𝑁−1(𝜏)−𝑁𝑎₁+2𝛼𝛼_𝑁(𝜏)−2𝑁𝜔𝜌_𝑁(𝜏)=𝑎₀+𝑎₁𝛼+𝛼²−𝜔²𝛼_𝑁−1(𝜏)+𝑎₁𝜔+2𝛼𝜔𝜌_𝑁−1(𝜏),𝑣_𝑁−1(𝜏)−𝑁𝑎₁+2𝛼𝜌_𝑁(𝜏)+2𝑁𝜔𝛼_𝑁(𝜏)=−𝜔𝑎₁+2𝛼𝛼_𝑁−1(𝜏)+𝑎₀+𝑎₁𝛼+𝛼²−𝜔²𝜌_𝑁−1(𝜏). (4.19) As before, this system can be solved by applying the Cramer’s rule. Finally, for any 𝑝<𝑁−1, we have to solve the system in 𝛼_𝑝(𝜏) and 𝜌_𝑝(𝜏), as follows:𝑢_𝑝(𝜏)−𝑝+2𝑝+1𝛼_𝑝+2(𝜏)−𝑝+1𝑎₁+2𝛼𝛼_𝑝+1(𝜏)−2𝑝+1𝜔𝜌_𝑝+1(𝜏)=𝑎₀+𝑎₁𝛼+𝛼²−𝜔²𝛼_𝑝(𝜏)+𝜔𝑎₁+2𝛼𝜌_𝑝(𝜏),𝑣_𝑝(𝜏)−𝑝+2𝑝+1𝜌_𝑝+2(𝜏)−𝑝+1𝑎₁+2𝛼𝜌_𝑝+1(𝜏)+2𝑝+1𝜔𝛼_𝑝+1(𝜏)=−𝜔𝑎₁+2𝛼𝛼_𝑝(𝜏)+𝑎₀+𝑎₁𝛼+𝛼²−𝜔²𝜌_𝑝(𝜏). (4.20)

Subcase 2.2. Now we consider the case where 𝛼=−𝑎₁/2 and 𝜔=4𝑎₀−𝑎²₁/2. This implies that 𝑎₁+2𝛼=0 and 𝛼²−𝜔²+𝑎₁𝛼+𝑎₀=0. In this case, the particular solution adopts the form𝑥(𝑡)=𝛼₀(𝜏)𝑡+𝛼₁(𝜏)𝑡²+𝛼₂(𝜏)𝑡³+⋯+𝛼_𝑁(𝜏)𝑡^𝑁+1𝑒^𝛼𝑡()+𝜌cos𝜔𝑡₀(𝜏)𝑡+𝜌₁(𝜏)𝑡²+𝜌₂(𝜏)𝑡³+⋯+𝜌_𝑁(𝜏)𝑡^𝑁+1𝑒^𝛼𝑡sin(𝜔𝑡), (4.21) and substituting 𝑥(𝑡), 𝑥̇(𝑡), and 𝑥̈(𝑡) into (4.4), we obtain𝜌_𝑁(𝜏)=𝑢_𝑁(𝜏)2(𝑁+1)𝜔,𝛼_𝑁(𝜏)𝑣=−^𝑁(𝜏)2()𝜔,𝜌𝑁+1_𝑝(𝜏)=𝑢_𝑝(𝜏)−𝑝+1𝑝+2𝛼_𝑝+12𝑝+1𝜔,𝛼_𝑝(𝜏)=𝑣_𝑝(𝜏)−𝑝+1𝑝+2𝜌_𝑝+1𝜔,−2𝑝+1 (4.22) for any 𝑝<𝑁.

From (4.8)–(4.22), it is a straightforward task; the derivation of an algorithm for the computation of the solution to (4.4).

5. Computation of the Solution to the Perturbed Problem

The symbolic package is thought to compute the solution to (1.3). The standard approach [13] is to try a power series solution of the form𝑥(𝑡,𝜖)=𝑥₀(𝑡)+𝑥₁(𝑡)𝜖+𝑥₂(𝑡)𝜖²+⋯. (5.1)This 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𝑥̈₀+𝑎₁𝑥̇₀+𝑎₀𝑥₀=𝑢(𝑡),𝑥𝑡₀=𝑥₀,𝑥̇𝑡₀=𝑥̇₀, (5.2) 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 𝑞≥1 is computed by solving the equation𝑥̈_𝑞+𝑎₁𝑥̇_𝑞+𝑎₀𝑥_𝑞=^𝑀_{𝜅=00≤𝜈≤𝜅}𝑓_{𝜈,𝜅−𝜈}𝑥^𝜈𝑥̇^𝜅−𝜈_𝑞−1,𝑥_𝑞𝑡₀=0,𝑥̇_𝑞𝑡₀=0, (5.3) 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 𝑞≥0, following the formulae𝑥^𝜈_𝑞=_{0≤𝑝≤𝑞}𝑥^𝜈−1_𝑝(𝑥)_𝑞−𝑝,𝑥̇^𝜈_𝑞=_{0≤𝑝≤𝑞}𝑥̇^𝜈−1_𝑝(𝑥̇)_𝑞−𝑝,𝑥^𝜈𝑥̇^𝜅−𝜈_𝑞=_{0≤𝑝≤𝑞}𝑥^𝜈_𝑝𝑥̇^𝜅−𝜈_𝑞−𝑝. (5.4)

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.(1)Define a three-dimensional array of quasipolynomials 𝑋(𝜌₁, 𝜌₂, 𝑞), where 𝜌₁,𝜌₂,𝑞∈ℕ, 0≤𝜌₁, 𝜌₂≤𝑀, and 0≤𝑞≤𝑄, 𝑄 being the order of the asymptotic expansion.(2)Define an array of quasipolynomials 𝑥(𝜌), where 0≤𝜌≤𝑄.(3)Initialize 𝑋(0,0,0)=1, and proceed as follows:(3.1) Compute 𝑋(1,0,0) as the solution to (5.2).(3.2) Compute 𝑋(0,1,0)=𝑑/𝑑𝑡(𝑋(1,0,0)).(3.3) Calculate, for each 𝜌 such that 2≤𝜌≤𝑀, 𝑋(),𝑋().𝜌,0,0=𝑋1,0,0×𝑋𝜌−1,0,00,𝜌,0=𝑋0,1,0×𝑋0,𝜌−1,0 (5.5)(3.4) For each 𝜌₁, 𝜌₂ such that 1≤𝜌₁, 𝜌₂≤𝑀, determine the quasipolynomial 𝑋𝜌₁,𝜌₂,0=𝑋𝜌₁,0,0×𝑋0,𝜌₂,0. (5.6)(4) For each 𝑞 such that 1≤𝑞≤𝑄, do the following.(4.1) Compute the quasipolynomial 𝑈=^𝑀_{𝜅=00≤𝜌≤𝜅}𝑓_{𝜌,𝜅−𝜌}×𝑋𝜌,𝜅−𝜌,𝑞−1. (5.7)(4.2) Calculate 𝑋(1,0,𝑞) as the solution to (5.3), 𝑥̈_𝑞+𝑎₁𝑥̇_𝑞+𝑎₀𝑥_𝑞=𝑈,𝑥_𝑞𝑡₀=0,𝑥̇_𝑞𝑡₀=0. (5.8)(4.3) Compute 𝑋(0,1,𝑞)=𝑑/𝑑𝑡(𝑋(1,0,𝑞)).(4.4) Calculate, for each 𝜌 such that 2≤𝜌≤𝑀, 𝑋𝜌,0,𝑞=_{0≤𝑝≤𝑞}𝑋,𝑋𝜌−1,0,𝑝×𝑋1,0,𝑞−𝑝0,𝜌,𝑞=_{0≤𝑝≤𝑞}𝑋0,𝜌−1,𝑝×𝑋0,1,𝑞−𝑝. (5.9)(4.5) For each 𝜌₁,𝜌₂ such that 1≤𝜌₁,𝜌₂≤𝑀, compute 𝑋𝜌₁,𝜌₂,𝑞=_{0≤𝑝≤𝑞}𝑋𝜌₁,0,𝑝×𝑋0,𝜌₂,𝑞−𝑝. (5.10)(5)For each 𝜌 such that 0≤𝜌≤𝑄, 𝑥𝜌=𝑋1,0,𝜌,𝑥(𝑡)=_{0≤𝜌≤𝑄}𝜖^𝜌𝑥𝜌. (5.11)

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),𝑓(𝑥,𝑥̇)=^𝑀_{𝜅=00≤𝜈≤𝜅}𝑓_{𝜈,𝜅−𝜈}𝑥^𝜈𝑥̇^𝜅−𝜈, (5.12)

it can be specified by a real (𝑀+1)×(𝑀+1) matrix,𝑓(𝑥,𝑥̇)=1𝑥⋯𝑥^𝑀𝑓₀₀𝑓₀₁⋯𝑓_0𝑀𝑓₁₀𝑓₁₁⋯𝑓_1𝑀𝑓⋮⋮⋱⋮_𝑀0𝑓_𝑀1⋯𝑓_𝑀𝑀1⋮𝑥̇𝑥̇^𝑀. (5.13)

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𝑥(𝑡)=_{0≤𝜌≤𝑄}𝜖^𝜌𝑥𝜌. (5.14)

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 xdx [𝙼+𝟷][𝙼+𝟷][𝙾𝚁𝙳𝙴𝚁+𝟷];          Define a three-dimensional array of quasi                  polynomials𝑋(𝜌1,𝜌2,𝑞), where 𝜌1,𝜌2,𝑞∈                  𝑁, 0≤𝜌1,𝜌2≤𝑀 and 0≤𝑞≤𝑄, being                  𝑄 the order of the asymptotic expansion. MQ *s;                 Define an array of quasipolynomials 𝑥(𝜌), s=new MQ [𝙾𝚁𝙳𝙴𝚁+𝟷];          where 0≤𝜌≤𝑄. int 𝚒,𝚓,𝚔,𝚙,𝚚,𝚗𝚞;         Definition of auxiliary variables. MQ U;

Algorithm 3 

Implementation of the algorithm. Steps (1) and (2).

xdx [𝟶][𝟶][𝟶].add            Initialize 𝑋(0,0,0)=1. (𝟷.,𝟶.,𝟶,𝟶.,𝟶.); xdx [𝟷][𝟶][𝟶] = u.solvePVI       Compute 𝑋(1,0,0) as the solution to (5.2), (a1, a0, t0, x0, dx0);       and 𝑋(0,1,0)=𝑑/𝑑𝑡(𝑋(1,0,0)). xdx [𝟶][𝟷][𝟶] = xdx [𝟷][𝟶][𝟶].der ();

Algorithm 4 

Implementation of the algorithm. Steps (3.1) and (3.2).

for (𝚒=2;𝚒<=𝚖;𝚒++){         Calculate, for each 𝜌 such that 2≤𝜌≤𝑀,     xdx [i][𝟶][𝟶] =       xdx [𝟷][𝟶][𝟶] *xdx [i−𝟷][𝟶][𝟶];      𝑋(𝜌,0,0)=𝑋(1,0,0)×𝑋(𝜌−1,0,0),     xdx [𝟶][i][𝟶] =       xdx [𝟶][𝟷][𝟶] *xdx [𝟶][i−𝟷][0];     𝑋(0,𝜌,0)=𝑋(0,1,0)×𝑋(0,𝜌−1,0). } for (i=1; i<=m; i++)          For each 𝜌1,𝜌2 such that 1≤𝜌1,𝜌2≤𝑀,   for (j=1; j<=m; j++)            determine the quasipolynomial     xdx [i][j][𝟶] =      xdx [i][𝟶][𝟶]*xdx [𝟶][𝚓][𝟶];   𝑋(𝜌1,𝜌2,0)=𝑋(𝜌1,0,0)×𝑋(0,𝜌2,0).

Algorithm 5 

Implementation of the algorithm. Steps (3.3) and (3.4).

for (q=1; q<=order; q++)         For each 𝑞 such that 1≤𝑞≤𝑄, do: {   U=U*0.;                  Compute the quasipolynomial      for (k=0; k<=m; k++)         for (nu=0; nu<=k; nu++)       𝑈=∑𝑀𝜅=0∑0≤𝜌≤𝜅𝑓𝜌,𝜅−𝜌×𝑋(𝜌,𝜅−𝜌,𝑞−1).         U = U + xdx [nu] [k-nu] [q-1]            *f [nu] [k-nu];  xdx [𝟷][𝟶][𝚚] = U.solvePVI         Calculate 𝑋(1,0,𝑞) as the solution to (5.3),  (a1, a0, t0, 0., 0.);  xdx [𝟶][𝟷][𝚚]= xdx [𝟷][𝟶][𝚚].der ();              𝑥̈𝑞+𝑎1𝑥̇𝑞+𝑎0𝑥𝑞=𝑈,                          𝑥𝑞(𝑡0)=0,𝑥̇𝑞(𝑡0)=0,                       and 𝑋(0,1,𝑞)=𝑑/𝑑𝑡(𝑋(1,0,𝑞)).  for (i=2; i<=m; i++)             Calculate, for each 𝜌 such that 2≤𝜌≤𝑀,    for (nu=0; nu<=q; nu++) {    xdx [i][𝟶][𝚚] = xdx[i][𝟶][𝚚]       𝑋(𝜌,0,𝑞)=∑0≤𝑝≤𝑞𝑋(𝜌−1,0,𝑝)         + xdx [i−𝟷][𝟶][𝚗𝚞]                         *xdx [𝟷][𝟶][𝚚−𝚗𝚞];           ×𝑋(1,0,𝑞−𝑝),    xdx [𝟶][i][𝚚] = xdx [𝟶][i][𝚚]      𝑋(0,𝜌,𝑞)=∑0≤𝑝≤𝑞𝑋(0,𝜌−1,𝑝)      + xdx [𝟶][i−𝟷][𝚗𝚞]         *xdx [𝟶][𝟷][𝚚−𝚗𝚞];           ×𝑋(0,1,𝑞−𝑝).   }  for (i=1; i<=m; i++)           For each 𝜌1,𝜌2 such that 1≤𝜌1,𝜌2≤𝑀,   for (j=1; j<=m; j++)         compute    for (nu=0; nu<=q; nu++)     xdx [i][𝚓][𝚚] =            𝑋(𝜌1,𝜌2,𝑞)      xdx [i][𝚓][𝚚]      + xdx [i][𝟶][𝚗𝚞]            =∑0≤𝑝≤𝑞𝑋(𝜌1,0,𝑝)×𝑋(0,𝜌2,𝑞−𝑝).          *xdx [𝟶][𝚓][𝚚−nu]; }

Algorithm 6 

Implementation of the algorithm. Step (4).

for (i=0; i<=order; i++){ For each 𝜌 such that 0≤𝜌≤𝑄,   s [i] = xdx [𝟷][𝟶][i];    s [i].normalice ();       𝑥(𝜌)=𝑋(1,0,𝜌),and𝑥(𝑡)=∑0≤𝜌≤𝑄𝜖𝜌𝑥(𝜌).   s [i].order ();   s [i].join ();   s [i].neglect (PREC);  }

Algorithm 7 

Implementation of the algorithm. Step (5).

MQ *solveAM (double a1, double a0, double t0, double x0, double dx0, MQ u, double **f, int m, int order)

7. Example

Let us consider the initial value problem given by𝑥̈+𝑥=𝜖−𝑥³−𝑥̇=𝜖𝑓(𝑥,𝑥̇),𝑥(0)=1,𝑥̇(0)=0, (7.1) where 𝜖=10⁻¹. Here, 𝑎₁=0, 𝑎₀=1, 𝑡₀=0, 𝑥₀=1, 𝑥̇₀=0, 𝑀=3, and 𝑓 is defined by the following 4×4 matrix:𝐹=0−10000000000−1000. (7.2)

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.

Table 2 

Some coefficients of the asymptotic expansion of the solution to (7.1).

In Figure 1, we show a comparison of the solution computed through the symbolic method (to the fourth-order of the solution) presented in this paper, and the numeric solution to the problem calculated by a Runge-Kutta fourth-order method with a step of ℎ=0.001. Let us stress that the difference between both solutions at any time is smaller than 10⁻⁷. The numerical solution has been calculated with the only purpose of validating the symbolic solution, as it is not the goal of this paper to develop a numerical tool.

Figure 1 

Comparison between a fourth order symbolic and a numerical solution.

In Table 3, we compare the CPU time of the algorithm implemented with MQSP and Maple running on an iMac 2.8 GHz Intel Core 2 Duo. For the test problem described by (7.1), we have calculated the solution up to order 14. In Table 3, we show the CPU time for the different orders of the solution from 𝑥₀(𝑡) (nonperturbed problem) to 𝑥₁₄(𝑡). The computation time of 𝑥₄(𝑡) with the use of the dsolve command of Maple in the implementation of the algorithm is larger than 100000 seconds. This is due to the fact that dsolve is a general solver which handles different types of ordinary differential equations. We have also computed the symbolic solution taking into account the linearity of (7.1) and implementing a routine following the algorithm detailed in Section 4. These CPU times are given in the fourth column of Table 3. For orders of the solution larger than 11, Maple does not give any response. We can see that the difference between MQSP and Maple in terms of CPU time is significant.

Table 3 

CPU time (in seconds) for the implementation with the MQSP, dsolve, and Maple.

8. Conclusions

In this paper, we have described an algorithm for the computation of the solution to a perturbed second-order differential equation through the asymptotic expansion technique. This algorithm has been implemented via a symbolic computation system which handles quasipolynomials.