# The Asymptotic Expansion Method via Symbolic Computation

Authors: Juan F. Navarro
Journal: Journal of Applied Mathematics
 (2012)

---

## 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.


---

## Body

## 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)∑i1,…,in∑j1,…,jmCi1,…,inj1,…,jmx1i1x2i2⋯xnincos⁡sin⁡(j1ϕ1+⋯+jmϕm),
where Ci1,…,inj1,…,jm∈ℝ, i1,…,in,j1,…,jm∈ℤ, and x1,…,xn and ϕ1,…,ϕm 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 (calledMQSP) based on the object-oriented philosophy which manipulates objects of the form(1.2)∑ν≥0τ1σ1×⋯×τQσQtnνeανt(λνcos⁡⁡(ωνt)+μνsin⁡⁡(ωνt)),
where nν∈ℕ,αν,ων,λν,μν∈ℝ, σν∈ℤ, and τ1,…,τQ 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 t, 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)ẍ+a1ẋ+a0x=u(t)+ϵf(x,ẋ),x(t0)=x0,ẋ(t0)=ẋ0,wherea0,a1,t0,x0,ẋ0∈ℝ, u(t) is given by (1.2), and(1.4)f(x,ẋ)=∑κ=0M∑0≤ν≤κfν,κ-νxνẋκ-ν,fν,κ-ν∈R.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 t0<t<t0+α (α>0) with an error not larger than O(α). 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 f(x,ẋ) 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 E be a set of symbolic objects, and let ~ be an equivalence relation in E, defined as follows: a~b if a=b, with a,b∈E. Here, the operator = is considered as the equality on the mathematical object level. Moreover, a≡b if a and b are identical as symbolic objects. A function f:E→E is said to be normal in (E,~) if f(a)~a for all a∈E, and f is said to be canonical in (E,~) if it is normal and a~b⇒f(a)=f(b) for all a,b∈E. 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 variableQ[t]. A quasipolynomial is a map u:ℝ→ℝ, defined by(2.1)u(t)=∑ν≥0tnνeανt(λνcos⁡(ωνt)+μνsin⁡(ωνt)),
where nν∈ℕ, αν,ων,λν, and μν∈ℝ. Let us now consider the set of quasipolynomials in the independent variable t, Q[t].


### 2.1. QuasiPolynomials as Symbolic Objects

We look for a canonical representation for each equivalence class defined inQ[t]. For that purpose, the following operations must be performed over each quasipolynomial. (1)
Let us consider a quasipolynomialu(t)=∑ν≥0tnνeανt(λνcos(ωνt)+μνsin(ωνt)). If ων<0, the following rules must be applied:(2.2)sin⁡(-ωνt)=-sin⁡(ωνt),cos⁡(-ωνt)=cos⁡(ωνt).(2)
The terms of a quasipolynomial will be ordered as follows: let us consider two term of a quasipolynomial:τ1=tn1eα1t(λ1cosω1t+μ2sinω1t) and τ2=tn2eα2t(λ2cosω2t+μ2cosω2t). We say that τ1<τ2 if (n1<n2) or (n1=n2 and α1<α2) or (n1=n2, α1=α2, and ω1<ω2) or (n1=n2, α1=α2, ω1=ω2, and λ1<λ2) or (n1=n2, α1=α2, ω1=ω2, λ1<λ2, and μ1<μ2).(3)
The terms of a quasipolynomial(2.3)u(t)=∑ν≥0tnνeανt(λνcos⁡(ωνt)+μνsin⁡(ωνt))with identical values ofnν,αν, 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)
n∈ℕ is the degree of the monomial t.(3)
α∈ℝ is the exponent of the exponential part of the term.(4)
ω≥0 is the frequency of the periodic part.In Scheme1 we represent the way a quasipolynomial is stored in memory by using a simple list.Scheme 1As 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 beO(logn).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:(2.4)u(t)=∑ν≥0Cp,ν(t)eανtcos⁡(ωνt)+Sq,ν(t)eανtsin⁡(ωνt),whereCp,ν(t) and Sq,ν(t) are polynomials in the variable t with constant coefficients, of degree p and q respectively, being p,q∈ℕ:(2.5)Cp,ν(t)=λ0+λ1t+λ2t2+⋯+λptp,Sq,ν(t)=μ0+μ1t+μ2t2+⋯+μqtq,
with λ0,λ1,…,λp,μ0,μ1,…,μq∈ℝ. 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 2In Scheme2 we show how we store in memory the quasipolynomial(2.6)u(t)=C1(t)eα1tcos⁡(ω1t)+S1(t)eα1tsin⁡(ω1t)+C2(t)eα2tcos⁡(ω2t)+S2(t)eα2tsin⁡(ω2t),
where(2.7)C1(t)=∑i∈IC,1λ1,itν1,i,S1(t)=∑i∈IS,1μ1,itσ1,i,C2(t)=∑i∈IC,2λ2,itν2,i,S2(t)=∑i∈IS,2μ2,itσ2,i,with𝕀C,j and 𝕀S,j being sets of indices for any j=1,2, ν1,i≥0 and λ1,i≠0 for any i∈𝕀C,1, ν2,i≥0 and λ2,i≠0 for any i∈𝕀C,2, σ1,i≥0 and μ1,i≠0 for any i∈𝕀S,1, and σ2,i≥0 and μ2,i≠0 for any i∈𝕀S,2. Moreover, νj,i-1<νj,i<νj,i+1 and σj,i-1<σj,i<σj,i+1 for any j=1, 2 and i 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 Scheme3 the structure associated to the quasipolynomial given by(2.8)u(t)=(3t+8t3-t7)e2tcos⁡⁡t+e2tsin⁡⁡t(t-t4+3t7).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)
(σ1,…,σQ)∈ℤQ. For each 0≤ν≤Q, σν represents the exponent of the undetermined coefficient τν.(3)
n∈ℕ is the degree of the monomial t.(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)
(m1,…,mQ)∈{0,1}Q. For each 0≤ν≤Q, mν=0 indicates that τν is a coefficient with unknown value, while mν=1 implies that τν has a real value assigned,(2)
τ=(τ1,…,τQ)∈ℝQ. If the coefficient τν has a real value assigned, its content is given by τν, for each 0≤ν≤Q.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 Table1 we illustrate the way a term is coded with a few examples, where it has been assumed that Q=4.Table 1
Some examples of representation of terms.
Termλμσ1σ2σ3σ4nαω3τ12τ4cos(5t)302001005-τ1τ2τ32t7e-tsin(3t)0−111207−132τ2τ32te2tcos(3t)200120123


## 2.1. QuasiPolynomials as Symbolic Objects

We look for a canonical representation for each equivalence class defined inQ[t]. For that purpose, the following operations must be performed over each quasipolynomial. (1)
Let us consider a quasipolynomialu(t)=∑ν≥0tnνeανt(λνcos(ωνt)+μνsin(ωνt)). If ων<0, the following rules must be applied:(2.2)sin⁡(-ωνt)=-sin⁡(ωνt),cos⁡(-ωνt)=cos⁡(ωνt).(2)
The terms of a quasipolynomial will be ordered as follows: let us consider two term of a quasipolynomial:τ1=tn1eα1t(λ1cosω1t+μ2sinω1t) and τ2=tn2eα2t(λ2cosω2t+μ2cosω2t). We say that τ1<τ2 if (n1<n2) or (n1=n2 and α1<α2) or (n1=n2, α1=α2, and ω1<ω2) or (n1=n2, α1=α2, ω1=ω2, and λ1<λ2) or (n1=n2, α1=α2, ω1=ω2, λ1<λ2, and μ1<μ2).(3)
The terms of a quasipolynomial(2.3)u(t)=∑ν≥0tnνeανt(λνcos⁡(ωνt)+μνsin⁡(ωνt))with identical values ofnν,αν, 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)
n∈ℕ is the degree of the monomial t.(3)
α∈ℝ is the exponent of the exponential part of the term.(4)
ω≥0 is the frequency of the periodic part.In Scheme1 we represent the way a quasipolynomial is stored in memory by using a simple list.Scheme 1As 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 beO(logn).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:(2.4)u(t)=∑ν≥0Cp,ν(t)eανtcos⁡(ωνt)+Sq,ν(t)eανtsin⁡(ωνt),whereCp,ν(t) and Sq,ν(t) are polynomials in the variable t with constant coefficients, of degree p and q respectively, being p,q∈ℕ:(2.5)Cp,ν(t)=λ0+λ1t+λ2t2+⋯+λptp,Sq,ν(t)=μ0+μ1t+μ2t2+⋯+μqtq,
with λ0,λ1,…,λp,μ0,μ1,…,μq∈ℝ. 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 2In Scheme2 we show how we store in memory the quasipolynomial(2.6)u(t)=C1(t)eα1tcos⁡(ω1t)+S1(t)eα1tsin⁡(ω1t)+C2(t)eα2tcos⁡(ω2t)+S2(t)eα2tsin⁡(ω2t),
where(2.7)C1(t)=∑i∈IC,1λ1,itν1,i,S1(t)=∑i∈IS,1μ1,itσ1,i,C2(t)=∑i∈IC,2λ2,itν2,i,S2(t)=∑i∈IS,2μ2,itσ2,i,with𝕀C,j and 𝕀S,j being sets of indices for any j=1,2, ν1,i≥0 and λ1,i≠0 for any i∈𝕀C,1, ν2,i≥0 and λ2,i≠0 for any i∈𝕀C,2, σ1,i≥0 and μ1,i≠0 for any i∈𝕀S,1, and σ2,i≥0 and μ2,i≠0 for any i∈𝕀S,2. Moreover, νj,i-1<νj,i<νj,i+1 and σj,i-1<σj,i<σj,i+1 for any j=1, 2 and i 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 Scheme3 the structure associated to the quasipolynomial given by(2.8)u(t)=(3t+8t3-t7)e2tcos⁡⁡t+e2tsin⁡⁡t(t-t4+3t7).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)
(σ1,…,σQ)∈ℤQ. For each 0≤ν≤Q, σν represents the exponent of the undetermined coefficient τν.(3)
n∈ℕ is the degree of the monomial t.(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)
(m1,…,mQ)∈{0,1}Q. For each 0≤ν≤Q, mν=0 indicates that τν is a coefficient with unknown value, while mν=1 implies that τν has a real value assigned,(2)
τ=(τ1,…,τQ)∈ℝQ. If the coefficient τν has a real value assigned, its content is given by τν, for each 0≤ν≤Q.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 Table1 we illustrate the way a term is coded with a few examples, where it has been assumed that Q=4.Table 1
Some examples of representation of terms.
Termλμσ1σ2σ3σ4nαω3τ12τ4cos(5t)302001005-τ1τ2τ32t7e-tsin(3t)0−111207−132τ2τ32te2tcos(3t)200120123


## 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, calledMQ, 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.Algorithm 1:C++ code for the definition of a modified quasipolynomial.
typedef Term *serie;class MQ{public:[definition of functions]private:serie first;};Algorithm 2:C++ code for the definition of a term 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;};The set of routines developed includes: addition, subtraction and product of series, multiplication of a series by a-real number, differentiation with respect tot, and computation of the solution to a linear second-order differential equation of type(3.1)ẍ+a1ẋ+a0x=u(t),
where u(t) 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, α2+a1α+a0=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 u(t) is written as follows:(4.1)u(t)=∑ν≥0eανt(pν,n(t,τ)cos⁡(ωνt)+qν,m(t,τ)sin⁡(ωνt)),
where α,ω∈ℝ, and pν,n(t,τ), and qν,m(t,τ) are nth and mth degree polynomials in t with undetermined coefficients respectively, of the form(4.2)pν,n(t,τ)=uν,0+uν,1t+⋯+uν,ntn,qν,m(t,τ)=vν,0+vν,1t+⋯+vν,mtm,
being(4.3)uν,ρ=uν,ρ*τ1σν,ρ,1×⋯×τQσν,ρ,Q,vν,ρ=vν,ρ*τ1σν,ρ,1*×⋯×τQσν,ρ,Q*,withuν,ρ*,vν,ρ*∈ℝ, σν,ρ,i,σν,ρ,i*∈ℤ, 1≤i≤Q, 0≤ρ≤n for uν,ρ, and 0≤ρ≤m for vν,ρ.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(4.4)ẍ+a1ẋ+a0x=eαt(pn(t,τ)cos⁡(ωt)+qm(t,τ)sin⁡(ωt)),
where(4.5)pn(t,τ)=u0+u1t+⋯+untn,qm(t,τ)=v0+v1t+⋯+vmtm,
being (4.6)uρ=uρ*τ1σρ,1×⋯×τQσρ,Q,vρ=vρ*τ1σρ,1*×⋯×τQσρ,Q*,
with uρ*,vρ*∈ℝ, σρ,i,σρ,i*∈ℤ, 1≤i≤Q, 0≤ρ≤n for uρ, and 0≤ρ≤m for vρ.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)ẍ+a1ẋ+a0x=pn(t,τ)eαt,
where pn(t,τ) is given by (4.5).Subcase 1.1.
Ifα2+a1α+a0≠0, the particular solution to the complete differential equation is expressed as(4.8)x(t)=(α0(τ)+α1(τ)t+α2(τ)t2+⋯+αn(τ)tn)eαt.
Thus, substituting x(t) and its derivatives in (4.7), we get that(4.9)αn(τ)=un(τ)α2+a1α+a0,αn-1(τ)=un-1(τ)-(a1+2α)nαn(τ)α2+a1α+a0,
and, for any p<n-1,(4.10)αp(τ)=up(τ)-(a1+2α)(p+1)αp+1(τ)-(p+1)(p+2)αp+2(τ)α2+a1α+a0.Subcase 1.2.
Ifα2+a1α+a0=0 and a12≠4a0 (i.e., α≠-a1/2), the particular solution can be written as(4.11)x(t)=(α0(τ)t+α1(τ)t2+α2(τ)t3+⋯+αn(τ)tn+1)eαt.
Now, the substitution of x(t) and its derivatives into (4.7) leads to the following formula for αn(τ),(4.12)αn(τ)=un(τ)(2α+a1)(n+1),
and, for any p≤n-1,(4.13)αp(τ)=up(τ)-(p+1)(p+2)αp+1(τ)(2α+a1)(p+1).Subcase 1.3.
Ifα2+a1α+a0=0 and a12=4a0 (i.e., α=-a1/2), the particular solution can be written as(4.14)x(t)=(α0(τ)t2+α1(τ)t3+α2(τ)t4+⋯+αn(τ)tn+2)eαt.
Now, the substitution of x(t) and its derivatives into (4.7) leads to the following formula for αp(τ),(4.15)αp(τ)=up(τ)(p+1)(p+2),p=0,…,n.Case 2.
Let us consider the second-order differential equation (4.4), where now ω≠0. There are two possible situations:(1)
α±iω is not a root of the characteristic equation, that is, α≠-a1/2 or ω≠4a0-a12/2,(2)
α±iω  is a root of the characteristic equation, that is, α=-a1/2 and ω=4a0-a12/2.Subcase 2.1.
Let us callN=max{n,m}. Then,(4.16)x(t)=(α0(τ)+α1(τ)t+⋯+αN(τ)tN)eαtcos⁡(ωt)+(ρ0(τ)+ρ1(τ)t+⋯+ρN(τ)tN)eαtsin⁡(ωt).
Now, if we assume pN(t,τ) and qN(t,τ) to be(4.17)pN(t,τ)=u0(τ)+u1(τ)t+⋯+uN(τ)tN,qN(t,τ)=v0(τ)+v1(τ)t+⋯+vN(τ)tN,
and substitute x(t) and its derivatives into (4.4), we get that(4.18)αN(τ)=1Δ|uN(τ)ω(a1+2α)vN(τ)a0+a1α+α2-ω2|,ρN(τ)=1Δ|a0+a1α+α2-ω2uN(τ)-ω(a1+2α)vN(τ)|,
where Δ=(a0+a1α+α2-ω2)2+ω2(a1+2α)2. Note that Δ≠0, because a0+a1α+α2-ω2≠0 or ω(2α+a1)≠0. Now, we can compute αN-1(τ) and ρN-1(τ) by solving the system(4.19)uN-1(τ)-N(a1+2α)αN(τ)-2NωρN(τ)=(a0+a1α+α2-ω2)αN-1(τ)+(a1ω+2αω)ρN-1(τ),vN-1(τ)-N(a1+2α)ρN(τ)+2NωαN(τ)=-ω(a1+2α)αN-1(τ)+(a0+a1α+α2-ω2)ρN-1(τ).
As before, this system can be solved by applying the Cramer’s rule. Finally, for any p<N-1, we have to solve the system in αp(τ) and ρp(τ), as follows:(4.20)up(τ)-(p+2)(p+1)αp+2(τ)-(p+1)(a1+2α)αp+1(τ)-2(p+1)ωρp+1(τ)=(a0+a1α+α2-ω2)αp(τ)+ω(a1+2α)ρp(τ),vp(τ)-(p+2)(p+1)ρp+2(τ)-(p+1)(a1+2α)ρp+1(τ)+2(p+1)ωαp+1(τ)=-ω(a1+2α)αp(τ)+(a0+a1α+α2-ω2)ρp(τ).Subcase 2.2.
Now we consider the case whereα=-a1/2 and ω=4a0-a12/2. This implies that a1+2α=0 and α2-ω2+a1α+a0=0. In this case, the particular solution adopts the form(4.21)x(t)=(α0(τ)t+α1(τ)t2+α2(τ)t3+⋯+αN(τ)tN+1)eαtcos⁡(ωt)+(ρ0(τ)t+ρ1(τ)t2+ρ2(τ)t3+⋯+ρN(τ)tN+1)eαtsin⁡(ωt),
and substituting x(t), ẋ(t), and ẍ(t) into (4.4), we obtain(4.22)ρN(τ)=uN(τ)2(N+1)ω,αN(τ)=-vN(τ)2(N+1)ω,ρp(τ)=up(τ)-(p+1)(p+2)αp+12(p+1)ω,αp(τ)=vp(τ)-(p+1)(p+2)ρp+1-2(p+1)ω,
for any p<N.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)x(t,ϵ)=x0(t)+x1(t)ϵ+x2(t)ϵ2+⋯.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 xi(t), which are then solved in a sequential manner. The resulting series need not converge for any value of ϵ; nevertheless, the solution x(t,ϵ) can be useful in approximating the function when ϵ is small.Considering the zero-order term inϵ yields(5.2)ẍ0+a1ẋ0+a0x0=u(t),x(t0)=x0,ẋ(t0)=ẋ0,
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 coefficientxq(t) of the solution to the order q≥1 is computed by solving the equation(5.3)ẍq+a1ẋq+a0xq=∑κ=0M∑0≤ν≤κfν,κ-ν(xνẋκ-ν)q-1,xq(t0)=0,ẋq(t0)=0,
where the notation (xνẋκ-ν)q refers to the qth order term of the series xνẋκ-ν.At each order of the solution, the series(xν)q, (ẋν)q, and (xνẋκ-ν)q must be computed once the order q has been solved, for each q≥0, following the formulae(5.4)(xν)q=∑0≤p≤q(xν-1)p(x)q-p,(ẋν)q=∑0≤p≤q(ẋν-1)p(ẋ)q-p,(xνẋκ-ν)q=∑0≤p≤q(xν)p(ẋκ-ν)q-p.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 quasipolynomialsX(ρ1, ρ2, q), where ρ1,ρ2,q∈ℕ, 0≤ρ1, ρ2≤M, and 0≤q≤Q, Q being the order of the asymptotic expansion.(2)
Define an array of quasipolynomialsx(ρ), where 0≤ρ≤Q.(3)
InitializeX(0,0,0)=1, and proceed as follows:
(3.1)
ComputeX(1,0,0) as the solution to (5.2).(3.2)
ComputeX(0,1,0)=d/dt(X(1,0,0)).(3.3)
Calculate, for eachρ such that 2≤ρ≤M,
(5.5)X(ρ,0,0)=X(1,0,0)×X(ρ-1,0,0),X(0,ρ,0)=X(0,1,0)×X(0,ρ-1,0).(3.4)
For eachρ1, ρ2 such that 1≤ρ1, ρ2≤M, determine the quasipolynomial
(5.6)X(ρ1,ρ2,0)=X(ρ1,0,0)×X(0,ρ2,0).(4)
For eachq such that 1≤q≤Q, do the following.
(4.1)
Compute the quasipolynomial(5.7)U=∑κ=0M∑0≤ρ≤κfρ,κ-ρ×X(ρ,κ-ρ,q-1).(4.2)
CalculateX(1,0,q) as the solution to (5.3),
(5.8)ẍq+a1ẋq+a0xq=U,xq(t0)=0,ẋq(t0)=0.(4.3)
ComputeX(0,1,q)=d/dt(X(1,0,q)).(4.4)
Calculate, for eachρ such that 2≤ρ≤M,
(5.9)X(ρ,0,q)=∑0≤p≤qX(ρ-1,0,p)×X(1,0,q-p),X(0,ρ,q)=∑0≤p≤qX(0,ρ-1,p)×X(0,1,q-p).(4.5)
For eachρ1,ρ2 such that 1≤ρ1,ρ2≤M, compute
(5.10)X(ρ1,ρ2,q)=∑0≤p≤qX(ρ1,0,p)×X(0,ρ2,q-p).(5)
For eachρ such that 0≤ρ≤Q,
(5.11)x(ρ)=X(1,0,ρ),x(t)=∑0≤ρ≤Qϵρx(ρ).The input arguments of the algorithm consist of the orderQ of the asymptotic approximation, the coefficients a1, a0 of the differential equation, the real values t0, x0, and ẋ0 which define the initial conditions of the problem, the quasipolynomial u(t), and the perturbed part of the equation f(x,ẋ). As f can be written as given in (1.4),(5.12)f(x,ẋ)=∑κ=0M∑0≤ν≤κfν,κ-νxνẋκ-ν,it can be specified by a real(M+1)×(M+1) matrix,(5.13)f(x,ẋ)=(1x⋯xM)(f00f01⋯f0Mf10f11⋯f1M⋮⋮⋱⋮fM0fM1⋯fMM)(1ẋ⋮ẋM).The parameterM is also required as input.The output argument of the algorithm is the arrayx(ρ) containing the coefficients of the asymptotic expansion. The asymptotic approximation to the order Q is given by(5.14)x(t)=∑0≤ρ≤Qϵρx(ρ).


## 6. On the Implementation of the Algorithm

In Algorithms3, 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.Algorithm 3:Implementation of the algorithm. Steps (1) and (2).
MQ xdx[M+1][M+1][ORDER+1];          Define a three-dimensional array of quasipolynomialsX(ρ1,ρ2,q), where ρ1,ρ2,q∈N, 0≤ρ1,ρ2≤M and 0≤q≤Q, beingQ  the order of the asymptotic expansion.MQ *s;                 Define an array of quasipolynomials x(ρ),s=new MQ[ORDER+1];          where 0≤ρ≤Q.inti,j,k,p,q,nu;         Definition of auxiliary variables.MQ U;Algorithm 4:Implementation of the algorithm. Steps (3.1) and (3.2).
xdx[0][0][0].add            Initialize X(0,0,0)=1.(1.,0.,0,0.,0.);xdx[1][0][0]= u.solvePVI       Compute X(1,0,0) as the solution to (5.2),(a1, a0, t0, x0, dx0);       and X(0,1,0)=d/dt(X(1,0,0)).xdx[0][1][0]= xdx[1][0][0].der ();Algorithm 5:Implementation of the algorithm. Steps (3.3) and (3.4).
for(i=2;i<=m;i++){         Calculate, for each ρ such that 2≤ρ≤M,xdx[i][0][0]=xdx[1][0][0]*xdx[i-1][0][0];      X(ρ,0,0)=X(1,0,0)×X(ρ-1,0,0),xdx[0][i][0]=xdx[0][1][0]*xdx[0][i-1][0];     X(0,ρ,0)=X(0,1,0)×X(0,ρ-1,0).}for (i=1; i<=m; i++)          For each ρ1,ρ2 such that 1≤ρ1,ρ2≤M,for (j=1; j<=m; j++)            determine the quasipolynomialxdx[i][j][0]=xdx[i][0][0]*xdx[0][j][0];   X(ρ1,ρ2,0)=X(ρ1,0,0)×X(0,ρ2,0).Algorithm 6:Implementation of the algorithm. Step (4).
for (q=1; q<=order; q++)         For each q such that 1≤q≤Q, do:{U=U*0.;                  Compute the quasipolynomialfor (k=0; k<=m; k++)for (nu=0; nu<=k; nu++)       U=∑κ=0M∑0≤ρ≤κfρ,κ-ρ×X(ρ,κ-ρ,q-1).U = U + xdx [nu] [k-nu] [q-1]*f [nu] [k-nu];xdx[1][0][q]= U.solvePVI         Calculate X(1,0,q) as the solution to (5.3),(a1, a0, t0, 0., 0.);xdx[0][1][q]= xdx[1][0][q].der ();              ẍq+a1ẋq+a0xq=U,xq(t0)=0,ẋq(t0)=0,andX(0,1,q)=d/dt(X(1,0,q)).for (i=2; i<=m; i++)             Calculate, for each ρ such that 2≤ρ≤M,for (nu=0; nu<=q; nu++){xdx[i][0][q]= xdx[i][0][q]       X(ρ,0,q)=∑0≤p≤qX(ρ-1,0,p)+ xdx[i-1][0][nu]*xdx[1][0][q-nu];           ×X(1,0,q-p),xdx[0][i][q]= xdx[0][i][q]      X(0,ρ,q)=∑0≤p≤qX(0,ρ-1,p)+ xdx[0][i-1][nu]*xdx[0][1][q-nu];           ×X(0,1,q-p).}for (i=1; i<=m; i++)           For each ρ1,ρ2 such that 1≤ρ1,ρ2≤M,for (j=1; j<=m; j++)         computefor (nu=0; nu<=q; nu++)xdx[i][j][q]=            X(ρ1,ρ2,q)xdx[i][j][q]+ xdx[i][0][nu]            =∑0≤p≤qX(ρ1,0,p)×X(0,ρ2,q-p).*xdx[0][j][q-nu];}Algorithm 7:Implementation of the algorithm. Step (5).
for (i=0; i<=order; i++){ For each ρ such that 0≤ρ≤Q,s [i] = xdx[1][0][i];s [i].normalice ();       x(ρ)=X(1,0,ρ),andx(t)=∑0≤ρ≤Qϵρx(ρ).s [i].order ();s [i].join ();s [i].neglect (PREC);}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(7.1)ẍ+x=ϵ(-x3-ẋ)=ϵf(x,ẋ),x(0)=1,ẋ(0)=0,
where ϵ=10-1. Here, a1=0, a0=1, t0=0, x0=1, ẋ0=0, M=3, and f is defined by the following 4×4 matrix:(7.2)F=(0-10000000000-1000).In Table2 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).
x0=+1.000000cos(t)x1=-0.031250cos(t)+0.031250cos(3t)+0.500000sin(t)-0.500000tcos(t)-0.375000tsin(t)x2=+0.022461cos(t)-0.023438cos(3t)+0.000977cos(5t)+0.226563sin(t)+0.070313sin(3t)-0.390625tcos(t)-0.046875tcos(3t)-0.031250tsin(t)-0.035156tsin(3t)+0.054688t2cos(t)+0.375000t2sin(t)x3=+0.047760cos(t)-0.046326cos(3t)-0.001465cos(5t)+0.000031cos(7t)+0.122640sin(t)-0.005859sin(3t)+0.003988sin(5t)-0.201660tcos(t)+0.079102tcos(3t)-0.002441tcos(5t)-0.133057tsin(t)-0.059692tsin(3t)-0.001831tsin(5t)+0.103760t2cos(t)+0.015381t2cos(3t)+0.140625t2sin(t)+0.070313t2sin(3t)+0.084635t3cos(t)-0.194336t3sin(t)x4=+0.031557cos(t)-0.025802cos(3t)-0.005688cos(5t)-0.000069cos(7t)+0.000001cos(9t)+0.314784sin(t)+0.001236sin(3t)-0.003563sin(5t)+0.000181sin(7t)-0.329605tcos(t)+0.018219tcos(3t)+0.009552tcos(5t)-0.000107tcos(7t)-0.441519tsin(t)+0.023666tsin(3t)-0.006212tsin(5t)-0.000080tsin(7t)+0.340805t2cos(t)-0.113068t2cos(3t)+0.001335t2cos(5t)+0.323486t2sin(t)-0.007141t2sin(3t)+0.005493t2sin(5t)-0.113566t3cos(t)+0.031860t3cos(3t)-0.168498t3sin(t)-0.064362t3sin(3t)-0.081533t4cos(t)+0.060547t4sin(t)x5=+0.101090cos(t)-0.103018cos(3t)+0.002327cos(5t)-0.000396cos(7t)-0.000003cos(9t)+0.000000cos(11t)+0.567049sin(t)+0.011443sin(3t)-0.003304sin(5t)-0.000300sin(7t)+0.000007sin(9t)-0.630124tcos(t)+0.042347tcos(3t)+0.004305tcos(5t)+0.000649tcos(7t)-0.000004tcos(9t)-0.861960tsin(t)-0.136819tsin(3t)+0.014968tsin(5t)-0.000407tsin(7t)-0.000003tsin(9t)+0.769419t2cos(t)+0.056821t2cos(3t)-0.019625t2cos(5t)+0.000082t2cos(7t)+0.653250t2sin(t)-0.009047t2sin(3t)-0.002311t2sin(5t)+0.000320t2sin(7t)-0.271820t3cos(t)+0.077614t3cos(3t)+0.003465t3cos(5t)-0.512211t3sin(t)+0.066607t3sin(3t)-0.007243t3sin(5t)+0.084625t4cos(t)-0.054769t4cos(3t)+0.110530t4sin(t)+0.027557t4sin(3t)+0.046019t5cos(t)-0.005018t5sin(t)x6=+0.041699cos(t)-0.035635cos(3t)-0.006540cos(5t)+0.000498cos(7t)-0.000022cos(9t)-0.000000cos(11t)+0.000000cos(13t)+1.753506sin(t)-0.124680sin(3t)+0.000293sin(5t)-0.000409sin(7t)-0.000018sin(9t)+0.000000sin(11t)-1.732879tcos(t)+0.361830tcos(3t)-0.007317tcos(5t)+0.000416tcos(7t)+0.000035tcos(9t)-0.000000tcos(11t)-1.950362tsin(t)+0.090323tsin(3t)-0.009679tsin(5t)+0.001762tsin(7t)-0.000022tsin(9t)-0.000000tsin(11t)+1.753296t2cos(t)-0.260570t2cos(3t)+0.014763t2cos(5t)-0.001785t2cos(7t)+0.000004t2cos(9t)+1.911049t2sin(t)+0.270222t2sin(3t)-0.022385t2sin(5t)-0.000220t2sin(7t)+0.000016t2sin(9t)-1.003710t3cos(t)-0.074659t3cos(3t)+0.018743t3cos(5t)+0.000255t3cos(7t)-1.212436t3sin(t)-0.124235t3sin(3t)+0.019212t3sin(5t)-0.000552t3sin(7t)+0.305366t4cos(t)-0.009847t4cos(3t)-0.008565t4cos(5t)+0.511480t4sin(t)-0.077788t4sin(3t)+0.004177t4sin(5t)-0.061540t5cos(t)+0.046341t5cos(3t)-0.067702t5sin(t)+0.007443t5sin(3t)-0.018889t6cos(t)-0.007806t6sin(t)x7=-0.030141cos(t)+0.026721cos(3t)+0.003535cos(5t)-0.000160cos(7t)+0.000046cos(9t)-0.000001cos(11t)-0.000000cos(13t)+0.000000cos(15t)+5.481497sin(t)+0.136811sin(3t)-0.022934sin(5t)+0.000391sin(7t)-0.000032sin(9t)-0.000001sin(11t)+0.000000sin(13t)-5.726072tcos(t)-0.099462tcos(3t)+0.047957tcos(5t)-0.002159tcos(7t)+0.000029tcos(9t)+0.000002tcos(11t)-0.000000tcos(13t)-5.240757tsin(t)+0.347041tsin(3t)+0.020306tsin(5t)-0.000669tsin(7t)+0.000136tsin(9t)-0.000001tsin(11t)-0.000000tsin(13t)+4.898568t2cos(t)-0.633927t2cos(3t)-0.017751t2cos(5t)+0.002044t2cos(7t)-0.000122t2cos(9t)+0.000000t2cos(11t)+5.910751t2sin(t)-0.423093t2sin(3t)+0.042183t2sin(5t)-0.003593t2sin(7t)-0.000015t2sin(9t)+0.000001t2sin(11t)-3.289085t3cos(t)+0.621827t3cos(3t)-0.037683t3cos(5t)+0.002208t3cos(7t)+0.000015t3cos(9t)-3.616263t3sin(t)-0.288847t3sin(3t)+0.002661t3sin(5t)+0.002305t3sin(7t)-0.000034t3sin(9t)+1.238058t4cos(t)-0.002983t4cos(3t)-0.002782t4cos(5t)-0.000822t4cos(7t)+1.664171t4sin(t)+0.252574t4sin(3t)-0.030580t4sin(5t)+0.000408t4sin(7t)-0.311598t5cos(t)-0.037610t5cos(3t)+0.009709t5cos(5t)-0.443419t5sin(t)+0.045286t5sin(3t)+0.002005t5sin(5t)+0.041874t6cos(t)-0.023862t6cos(3t)+0.045263t6sin(t)-0.022910t6sin(3t)+0.005561t7cos(t)+0.006555t7sin(t)x8=+0.707697cos(t)-0.737652cos(3t)+0.029647cos(5t)+0.000299cos(7t)+0.000006cos(9t)+0.000003cos(11t)-0.000000cos(13t)-0.000000cos(15t)+17.859287sin(t)+0.537553sin(3t)+0.020622sin(5t)-0.001769sin(7t)+0.000063sin(9t)-0.000002sin(11t)-0.000000sin(13t)+0.000000sin(15t)-19.388600tcos(t)-0.105739tcos(3t)-0.072738tcos(5t)+0.004113tcos(7t)-0.000258tcos(9t)+0.000002tcos(11t)+0.000000tcos(13t)-0.000000tcos(15t)-16.767762tsin(t)-0.906868tsin(3t)+0.089594tsin(5t)+0.000388tsin(7t)-0.000052tsin(9t)+0.000008tsin(11t)-0.000000tsin(13t)-0.000000tsin(15t)+16.120089t2cos(t)+0.421931t2cos(3t)-0.092803t2cos(5t)+0.001418t2cos(7t)+0.000192t2cos(9t)-0.000007t2cos(11t)+0.000000t2cos(13t)+19.608050t2sin(t)-0.420458t2sin(3t)-0.128984t2sin(5t)+0.005901t2sin(7t)-0.000351t2sin(9t)-0.000001t2sin(11t)+0.000000t2sin(13t)-11.042356t3cos(t)+0.875871t3cos(3t)+0.101426t3cos(5t)-0.006986t3cos(7t)+0.000187t3cos(9t)+0.000001t3cos(11t)-12.224352t3sin(t)+0.801292t3sin(3t)-0.047084t3sin(5t)+0.002053t3sin(7t)+0.000194t3sin(9t)-0.000002t3sin(11t)+4.715957t4cos(t)-0.906645t4cos(3t)+0.033974t4cos(5t)-0.000385t4cos(7t)-0.000061t4cos(9t)+5.809016t4sin(t)+0.172731t4sin(3t)+0.039748t4sin(5t)-0.004671t4sin(7t)+0.000031t4sin(9t)-1.468495t5cos(t)+0.132440t5cos(3t)-0.017574t5cos(5t)+0.001171t5cos(7t)-1.830821t5sin(t)-0.280024t5sin(3t)+0.025814t5sin(5t)+0.000249t5sin(7t)+0.278969t6cos(t)+0.046180t6cos(3t)-0.006166t6cos(5t)+0.377937t6sin(t)-0.001896t6sin(3t)-0.006983t6sin(5t)-0.028889t7cos(t)+0.004692t7cos(3t)-0.029867t7sin(t)+0.021274t7sin(3t)-0.000816t8cos(t)-0.003373t8sin(t)x9=+2.175721cos(t)-2.094928cos(3t)-0.084502cos(5t)+0.003726cos(7t)-0.000018cos(9t)+0.000001cos(11t)+0.000000cos(13t)-0.000000cos(15t)-0.000000cos(17t)+64.332742sin(t)+0.333348sin(3t)+0.044040sin(5t)+0.002384sin(7t)-0.000104sin(9t)+0.000006sin(11t)-0.000000sin(13t)-0.000000sin(15t)+0.000000sin(17t)-68.359632tcos(t)+2.878461tcos(3t)-0.079121tcos(5t)-0.008828tcos(7t)+0.000335tcos(9t)-0.000022tcos(11t)+0.000000tcos(13t)+0.000000tcos(15t)-0.000000tcos(17t)-64.047243tsin(t)-1.871954tsin(3t)-0.188447tsin(5t)+0.009868tsin(7t)-0.000169tsin(9t)-0.000004tsin(11t)+0.000000tsin(13t)-0.000000tsin(15t)-0.000000tsin(17t)+61.716104t2cos(t)-0.768222t2cos(3t)+0.295936t2cos(5t)-0.011360t2cos(7t)+0.000383t2cos(9t)+0.000014t2cos(11t)-0.000000t2cos(13t)+0.000000t2cos(15t)+69.633323t2sin(t)+2.498056t2sin(3t)-0.113525t2sin(5t)-0.012762t2sin(7t)+0.000674t2sin(9t)-0.000026t2sin(11t)-0.000000t2sin(13t)+0.000000t2sin(15t)-40.085044t3cos(t)-0.028980t3cos(3t)+0.047373t3cos(5t)+0.009129t3cos(7t)-0.000811t3cos(9t)+0.000013t3cos(11t)+0.000000t3cos(13t)-46.345112t3sin(t)+0.014214t3sin(3t)+0.328854t3sin(5t)-0.010719t3sin(7t)+0.000291t3sin(9t)+0.000013t3sin(11t)-0.000000t3sin(13t)+18.758768t4cos(t)-1.045383t4cos(3t)-0.199353t4cos(5t)+0.009593t4cos(7t)-0.000040t4cos(9t)-0.000004t4cos(11t)+21.493321t4sin(t)-0.827863t4sin(3t)-0.028074t4sin(5t)+0.005286t4sin(7t)-0.000478t4sin(9t)+0.000002t4sin(11t)-6.218301t5cos(t)+0.995799t5cos(3t)+0.007797t5cos(5t)-0.003375t5cos(7t)+0.000105t5cos(9t)-7.652998t5sin(t)+0.017785t5sin(3t)-0.070424t5sin(5t)+0.004954t5sin(7t)+0.000022t5sin(9t)+1.549574t6cos(t)-0.229644t6cos(3t)+0.028038t6cos(5t)-0.000898t6cos(7t)+1.852341t6sin(t)+0.213697t6sin(3t)-0.008050t6sin(5t)-0.001065t6sin(7t)-0.230798t7cos(t)-0.028635t7cos(3t)+0.000743t7cos(5t)-0.309074t7sin(t)-0.025324t7sin(3t)+0.008028t7sin(5t)+0.019906t8cos(t)+0.004745t8cos(3t)+0.017987t8sin(t)-0.012705t8sin(3t)-0.000299t9cos(t)+0.001329t9sin(t)In Figure1, 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 h=0.001. Let us stress that the difference between both solutions at any time is smaller than 10-7. 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 Table3, 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 x0(t) (nonperturbed problem) to x14(t). The computation time of x4(t) 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 theMQSP, dsolve, and Maple.
xnMQSPMaple dsolve(using dsolve)Maplex00.0001870.0590.042x10.0008470.1310.078x20.0046370.1960.147x30.0205430.3000.299x40.070672—0.622x50.209268—1.330x60.538049—2.694x71.255629—5.426x82.711981—10.715x95.421418—21.052x1010.367602—39.304x1118.494896—73.471x1231.738276——x1352.468991——x1483.748298——


## 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.



---

*Source: 201340-2012-06-09.xml*