# Symbolic Solution to Complete Ordinary Differential Equations with Constant Coefficients

Authors: Juan F. Navarro; Antonio Pérez-Carrió
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

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


---

## Body

## 1. Introduction

Perturbation theories for differential equations containing a small parameterϵ are quite old. The small perturbation theory originated by Sir Isaac Newton has been highly developed by many others, and an extension of this theory to the asymptotic expansion, consisting of a power series expansion in the small parameter, was devised by Poincaré [1]. The main point is that, for the most of the differential equations, it is not possible to obtain an exact solution. In cases where equations contain a small parameter, we can consider it as a perturbation parameter to obtain an asymptotic expansion of the solution. In practice, the work involved in the application of this approach to compute the solution to a differential equation cannot be performed by hand, and algebraic processors result in being a very useful tool.As explained in Henrard [2], the first symbolic processors were developed to work with Poisson series, that is, multivariate Fourier series whose coefficients are multivariate Laurent series,
(1)∑i1,…,in∑j1,…,jm‍Ci1,…,inj1,…,jmx1i1⋯xnincossin(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 equation problems, in the field of celestial mechanics. One of the first applications of these processor was concerned with the theory of the Moon. Delaunay invented his perturbation method to treat it and spent 20 years doing algebraic manipulations by hand to apply it to the problem. Deprit et al. [3, 4] prolongated the solution of Delaunay's work with the help of a special purpose symbolic processor, and Henrard [5] pushed to order 25. This solution was improved by iteration by Chapront-Touze [6], and planetary perturbations were also introduced by Chapront-Touze [7]. At present, the most complete solution, Ephemeride Lunaire Parisien (ELP) contains more than 50 000 periodic terms. But the motion of the Moon is not the only application of algebraic processors. There are many problems where the facilities provided by Poisson series processors can lead rather quickly to very accurate results. As examples, we would like to mention planetary theories, the theory of the rotation of the Earth (see, e.g. Navarro and Ferrándiz [8]), and artificial satellite theories.In order to achieve better accuracies in the applications of analytical theories, high orders of the approximate solution must be computed, making necessary a continuous maintenance and revision of the existing symbolic manipulation systems, as well as the development of new packages adapted to the peculiarities of the problem to be treated. Recently, Navarro [9, 10] developed a symbolic processor to deal with the solution to second order differential equations of the form
(2)x¨+a1x˙+a0x=u(t)+ϵf(x,x˙),
with initial conditions
(3)x(0)=x0,x˙(0)=x˙0,
where a0, a1, x0, and x˙0∈ℝ, u(t) is a quasipolynomial, that is, an object of the form
(4)u(t)=∑ν≥0‍tnνeανt(λνcos(ωνt)+μνsin(ωνt)),
where nν∈ℕ, αν, ων, λν, and μν∈ℝ, and f(x,x˙) admits the expansion
(5)f(x,x˙)=∑κ=0M‍∑0≤ν≤κ‍fν,κ-νxνx˙κ-ν,fν,κ-ν∈ℝ.
To face this problem, the algebraic processor handles objects called quasipolynomials, and it has resulted in being a useful tool in the computation of the solution to (2) by the application of the asymptotic expansion method. A modification of this processor has been employed to compute periodic solutions in equations of type (2) via the Poincaré-Lindstedt method in Navarro [9, 10]. To that end, the idea is to expand both the solution and the modified frequency with respect to the small parameter, allowing to kill secular terms which appear in the recursive scheme. The elimination of secular terms is performed through a manipulation system which works with modified quasipolynomials, that is, quasipolynomials containing undetermined constants:
(6)u(t)=∑ν≥0‍τ1σ1×⋯×τQσQtnνeανt×(λνcos(ωνt)+μνsin(ωνt)),
where nν∈ℕ, αν, ων, λν, and μν∈ℝ, σν∈ℤ, and τ1,…,τQ are real constants with unknown value.One year later, Navarro [11] presented a symbolic computation package based on the object-oriented philosophy for manipulating matrices whose elements lie on the set of quasipolynomials of type (4). The kernel of the symbolic processor was developed in C++, defining a class for this new object as well as a set of functions that operate on the data structure: addition, substraction, differentiation, and integration with respect to t, substitution of an undetermined coefficient by a series, and many others. The goal of this processor is to provide a tool to solve a perturbed n-order differential equation of the class
(7)x(n)+an-1x(n-1)+⋯+a0x=u(t)+ϵf(x,x˙,…,x(n-1)),
with initial conditions
(8)x(0)=x10,x˙(0)=x20,…,x(n-1)(0)=xn0,
where ϵ is a small real parameter, a0,a1,…,an-1∈ℝ, u(t) is a quasipolynomial, and f is such that
(9)f(x,x˙,…,x(n-1))=∑0≤ν1,…,νn≤M‍fν1,…,νnxν1⋯(x(n-1))νn,
with M∈ℕ, ν1,…,νn∈ℕ, and fν1,…,νn∈ℝ.As a first step, Navarro and Pérez [12] have developed a symbolic technique for the computation of the principal matrix of the linear system associated with a homogeneous ordinary differential equation with constant coefficients of the form
(10)x(n)+an-1x(n-1)+⋯+a1x˙+a0x=0,
where a0,a1,…,an-1∈ℝ, with initial conditions
(11)x(0)=x10,x˙(0)=x20,…,x(n-1)(0)=xn0
being x10,…,xn0∈ℝ. This method provides a final analytical solution which can be completely computed in a symbolic way.The second step in the construction of the solution to (7) is to obtain a symbolic procedure for computing the solution to the nonhomogeneus equation
(12)x(n)+an-1x(n-1)+⋯+a0x=u(t),
where, as above, u(t) is a quasipolynomial.This new symbolic package can be useful for research and educational purposes in the field of differential equations or dynamical systems. Nowadays there are many open problems which requires massive symbolic computation as, to cite one example, the analytical theory of the resonant motion of Mercury. We would like to stress that the aim of this work is to develop a symbolic tool, not a numeric method, as numeric solutions cannot be used in perturbation methods for differential equations. However, we have performed comparisons between the symbolic solution with a numeric one, just to show the efficiency of the technique.In next section, we summarize the scheme proposed in [12] to calculate the solution to (10), which is needed to express the solution to the complete problem (12). The technique to construct this solution will be described in Section 3.


## 2. Solution to the Homogeneous Problem

As mentioned, Navarro and Pérez [12] have proposed a symbolic method to compute the solution to (10)
(13)x(n)+an-1x(n-1)+⋯+a1x˙+a0x=0,
where a0,a1,…,an-1∈ℝ, with initial conditions
(14)x(0)=x10,x˙(0)=x20,…,x(n-1)(0)=xn0
being x10,…,xn0∈ℝ. This method provides a final analytical solution which can be completely computed in a symbolic way.


### 2.1. Description of the Method

With the aid of the substitutions(15)x1=x,x2=x˙,…,xn=x(n-1),
(10) is transformed into the system of differential equations given by
(16)X˙(t)=AX(t),X(0)=X0,
where A is the companion matrix,
(17)A=(010⋯0001⋯0⋮⋮⋮⋱⋮-a0-a1-a2⋯-an-1),X(t)=(x1(t)x2(t)⋮xn(t)),X0=(x10x20⋮xn0).
To compute the exponential of A, the matrix is splitted into B+C, where
(18)B=(010⋯0001⋯0⋮⋮⋮⋱⋮000⋯0),C=(000⋯0000⋯0⋮⋮⋮⋱⋮-a0-a1-a2⋯-an-1).
Then, we use the approximation
(19)eA≈(eB/meC/m)m,
with m∈ℕ. This approach to obtain eA is of potential interest when the exponentials of matrices B and C can be efficiently computed and requires an-1 to be nonequal to zero. In our case,
(20)eB/m=(g0(m)g1(m)g2(m)⋯gn-1(m)0g0(m)g1(m)⋯gn-2(m)00g0(m)⋯gn-3(m)⋮⋮⋮⋱⋮000⋯g0(m)),eC/m=I-1an-1(e-an-1/m-1)C=(100⋯0010⋯0001⋯0⋮⋮⋮⋱⋮ω0ω1ω2⋯ωn-1),
being, for any n∈ℕ,
(21)gn(m)=1n!mn,
and, for any i=0,…,n-2,
(22)ωi=aian-1(e-an-1/m-1),ωn-1=e-an-1/m.
As it is established in  Moler and van Loan [13], for a general splitting A=B+C, m can be determine from
(23)∥eA-(eB/meC/m)m∥≤12m∥[B,C]∥e∥B∥+∥C∥.
Thus, with A being the companion matrix, we get that
(24)∥eA-(eB/meC/m)m∥≤1me1+∥a∥∥a∥,
where a=(a0,a1,…,an-1). Here,
(25)∥x∥=∥x∥1=|x1|+⋯+|xn|,
for any x=(x1,…,xn)∈ℝn, and ∥A∥=max∥x∥1=1∥Ax∥1=maxν(|a1ν|+⋯+|anν|), for any matrix
(26)A=(a11a12⋯a1na21a22⋯a2n⋮⋮⋱⋮an1an2⋯ann).


### 2.2. Adaptation to a Symbolic Formalism

Navarro and Pérez [12, 14] have proposed an adaptation of the method described above in order to compute the matrix eAt instead of eA. If we do so, we obtain the principal matrix of (16), whose elements lie on the set of quasipolynomials, and the symbolic processor results in being suitable to work with those matrices. The approach is to compute
(27)eAt≃(eBt/meCt/m)m,
taking into account that
(28)eBt/m=(g0(t,m)g1(t,m)⋯gn-1(t,m)0g0(t,m)⋯gn-2(t,m)⋮⋮⋱⋮00⋯g0(t,m)),eCt/m=(100⋯0010⋯0001⋯0⋮⋮⋮⋱⋮ω0(t)ω1(t)ω2(t)⋯ωn-1(t)),
being, for any n∈ℕ,
(29)gn(t,m)=tnn!mn,
and, for any i=0,…,n-2,
(30)ωi(t)=aian-1(e-an-1t/m-1),ωn-1(t)=e-an-1t/m.Thus,eAt is a matrix of quasipolynomials that can be completely computed through the symbolic processor developed by Navarro [11]. The procedure for the computation of the exponential matrix is substantially simplified by using the following equation, which avoids the symbolic multiplication of eBt/m and eCt/m:
(31)eBt/meCt/m=eBt/m+H(t,m),
where
(32)H(t,m)=1an-1f(t,m)(gn-1(t,m)⋮g0(t,m))(a0a1⋯an-1),f(t,m)=e-an-1t/m-1.
This relation can be obtained directly from the multiplication of matrices eBt/m and eCt/m as given in (28).


### 2.3. Symbolic Expansion of the Exponential Matrix

As it has been stated in Section2.2, the exponential matrix eAt can be symbolically calculated through (27) and (31)
(33)eAt≃(eBt/meCt/m)m=(eBt/m+H(t,m))m.
This expression can be calculated symbolically. Let us express H(t,m) as
(34)H(t,m)=λ(t,m)J(t,m),
with
(35)λ(t,m)=1an-1f(t,m),J(t,m)=(gn-1(t,m)⋮g0(t,m))(a0a1⋯an-1).
Thus,
(36)(eBt/m+H(t,m))m=∑k=0m‍λk(t,m)(eBt/m,J(t,m))km,
where (A,B)km is the so-called noncommutative parenthesis, defined as
(37)(A,B)km=∑i=1(mk)∏j=1σ(m,k)(A(1-(-1)j)/2B(1+(-1)j)/2)cij‍,
where A and B are two noncommutative n×n matrices, m,k,cij∈ℤ+∪{0},
(38)σ(m,k)={2k+1,if2k≤m,2(m-k+1),if2k>m,
and cij is determined by the following properties:(1)
∑j=1σ(m,k)‍cij=m for all i=1,2,…,(mk);(2)
∑j=1[σ(m,k)/2]‍ci(2j)=m-k for all i=1,2,…,(mk),  with [σ(m,k)/2] being the entire part of σ(m,k)/2;(3)
ci2≠0 for all i=1,2,…,(mk);(4)
ifcij≠0 and ci(j+2)≠0, then ci(j+1)≠0 for all i=1,2,…,(mk) and j=1,2,…,σ(m,k)-2.In the following, we summarize some expressions which simplify the way in which the matrix(eBt/m+H(t,m))m is symbolically computed. First, let us introduce the following matrices:
(39)Sn-1(t,m)=(a0a1⋯an-1)F(t,m),F(t,m)=(tn-1(n-1)tn-2m(n-1)(n-2)tn-3m2⋮(n-1)!t0mn-1).Lemma 1.
For anyk∈ℤ such that k>1,
(40)(eBt/m)k=(G0,k(t,m)G1,k(t,m)⋯Gn-1,k(t,m)0G0,k(t,m)⋯Gn-2,k(t,m)⋮⋮⋱⋮00⋯G0,k(t,m)),
with
(41)Gν,k(t,m)=kνgν(t,m),
for each ν=0,…,n-1, where gν(t,m) is given by (29).Lemma 2.
For anyp∈ℤ such that p>1,
(42)(H(t,m))p=(f(t,m)an-1(n-1)!mn-1)p(Sn-1(t,m))p-1×F(t,m)(a0a1⋯an-1),
where Sn-1(t,m) and F(t,m) are given by (39).Lemma 3.
For anyp,k∈ℤ such that p,k≥1,
(43)(eBt/m)k(H(t,m))p=(f(t,m)an-1(n-1)!mn-1)p(Sn-1(t))p-1×ΩF(t,m)(a0a1⋯an-1),
where Sn-1(t,m) and F(t,m) are given by (39), and
(44)Ω=((k+1)n-10⋯00(k+1)n-2⋯0⋮⋮⋱⋮00⋯1).Lemma 4.
For anyp,k∈ℤ such that p,k≥1,
(45)Hp(t,m)(eBt/m)k=(f(t,m)an-1(n-1)!mn-1)p×(Sn-1(t))p-1F(t,m)A(t,m),
where Sn-1(t,m) and F(t,m) are given by (39),
(46)A(t,m)=(A0(t,m)A1(t,m)⋯An-1(t,m)),Aμ(t,m)=∑ν=0μ‍aνkμ-νgμ-ν,
for any μ=0,…,n-1, and gν is given by (29).For the sake of simplicity, we have omitted the dependence onm and t of gν.


## 2.1. Description of the Method

With the aid of the substitutions(15)x1=x,x2=x˙,…,xn=x(n-1),
(10) is transformed into the system of differential equations given by
(16)X˙(t)=AX(t),X(0)=X0,
where A is the companion matrix,
(17)A=(010⋯0001⋯0⋮⋮⋮⋱⋮-a0-a1-a2⋯-an-1),X(t)=(x1(t)x2(t)⋮xn(t)),X0=(x10x20⋮xn0).
To compute the exponential of A, the matrix is splitted into B+C, where
(18)B=(010⋯0001⋯0⋮⋮⋮⋱⋮000⋯0),C=(000⋯0000⋯0⋮⋮⋮⋱⋮-a0-a1-a2⋯-an-1).
Then, we use the approximation
(19)eA≈(eB/meC/m)m,
with m∈ℕ. This approach to obtain eA is of potential interest when the exponentials of matrices B and C can be efficiently computed and requires an-1 to be nonequal to zero. In our case,
(20)eB/m=(g0(m)g1(m)g2(m)⋯gn-1(m)0g0(m)g1(m)⋯gn-2(m)00g0(m)⋯gn-3(m)⋮⋮⋮⋱⋮000⋯g0(m)),eC/m=I-1an-1(e-an-1/m-1)C=(100⋯0010⋯0001⋯0⋮⋮⋮⋱⋮ω0ω1ω2⋯ωn-1),
being, for any n∈ℕ,
(21)gn(m)=1n!mn,
and, for any i=0,…,n-2,
(22)ωi=aian-1(e-an-1/m-1),ωn-1=e-an-1/m.
As it is established in  Moler and van Loan [13], for a general splitting A=B+C, m can be determine from
(23)∥eA-(eB/meC/m)m∥≤12m∥[B,C]∥e∥B∥+∥C∥.
Thus, with A being the companion matrix, we get that
(24)∥eA-(eB/meC/m)m∥≤1me1+∥a∥∥a∥,
where a=(a0,a1,…,an-1). Here,
(25)∥x∥=∥x∥1=|x1|+⋯+|xn|,
for any x=(x1,…,xn)∈ℝn, and ∥A∥=max∥x∥1=1∥Ax∥1=maxν(|a1ν|+⋯+|anν|), for any matrix
(26)A=(a11a12⋯a1na21a22⋯a2n⋮⋮⋱⋮an1an2⋯ann).


## 2.2. Adaptation to a Symbolic Formalism

Navarro and Pérez [12, 14] have proposed an adaptation of the method described above in order to compute the matrix eAt instead of eA. If we do so, we obtain the principal matrix of (16), whose elements lie on the set of quasipolynomials, and the symbolic processor results in being suitable to work with those matrices. The approach is to compute
(27)eAt≃(eBt/meCt/m)m,
taking into account that
(28)eBt/m=(g0(t,m)g1(t,m)⋯gn-1(t,m)0g0(t,m)⋯gn-2(t,m)⋮⋮⋱⋮00⋯g0(t,m)),eCt/m=(100⋯0010⋯0001⋯0⋮⋮⋮⋱⋮ω0(t)ω1(t)ω2(t)⋯ωn-1(t)),
being, for any n∈ℕ,
(29)gn(t,m)=tnn!mn,
and, for any i=0,…,n-2,
(30)ωi(t)=aian-1(e-an-1t/m-1),ωn-1(t)=e-an-1t/m.Thus,eAt is a matrix of quasipolynomials that can be completely computed through the symbolic processor developed by Navarro [11]. The procedure for the computation of the exponential matrix is substantially simplified by using the following equation, which avoids the symbolic multiplication of eBt/m and eCt/m:
(31)eBt/meCt/m=eBt/m+H(t,m),
where
(32)H(t,m)=1an-1f(t,m)(gn-1(t,m)⋮g0(t,m))(a0a1⋯an-1),f(t,m)=e-an-1t/m-1.
This relation can be obtained directly from the multiplication of matrices eBt/m and eCt/m as given in (28).


## 2.3. Symbolic Expansion of the Exponential Matrix

As it has been stated in Section2.2, the exponential matrix eAt can be symbolically calculated through (27) and (31)
(33)eAt≃(eBt/meCt/m)m=(eBt/m+H(t,m))m.
This expression can be calculated symbolically. Let us express H(t,m) as
(34)H(t,m)=λ(t,m)J(t,m),
with
(35)λ(t,m)=1an-1f(t,m),J(t,m)=(gn-1(t,m)⋮g0(t,m))(a0a1⋯an-1).
Thus,
(36)(eBt/m+H(t,m))m=∑k=0m‍λk(t,m)(eBt/m,J(t,m))km,
where (A,B)km is the so-called noncommutative parenthesis, defined as
(37)(A,B)km=∑i=1(mk)∏j=1σ(m,k)(A(1-(-1)j)/2B(1+(-1)j)/2)cij‍,
where A and B are two noncommutative n×n matrices, m,k,cij∈ℤ+∪{0},
(38)σ(m,k)={2k+1,if2k≤m,2(m-k+1),if2k>m,
and cij is determined by the following properties:(1)
∑j=1σ(m,k)‍cij=m for all i=1,2,…,(mk);(2)
∑j=1[σ(m,k)/2]‍ci(2j)=m-k for all i=1,2,…,(mk),  with [σ(m,k)/2] being the entire part of σ(m,k)/2;(3)
ci2≠0 for all i=1,2,…,(mk);(4)
ifcij≠0 and ci(j+2)≠0, then ci(j+1)≠0 for all i=1,2,…,(mk) and j=1,2,…,σ(m,k)-2.In the following, we summarize some expressions which simplify the way in which the matrix(eBt/m+H(t,m))m is symbolically computed. First, let us introduce the following matrices:
(39)Sn-1(t,m)=(a0a1⋯an-1)F(t,m),F(t,m)=(tn-1(n-1)tn-2m(n-1)(n-2)tn-3m2⋮(n-1)!t0mn-1).Lemma 1.
For anyk∈ℤ such that k>1,
(40)(eBt/m)k=(G0,k(t,m)G1,k(t,m)⋯Gn-1,k(t,m)0G0,k(t,m)⋯Gn-2,k(t,m)⋮⋮⋱⋮00⋯G0,k(t,m)),
with
(41)Gν,k(t,m)=kνgν(t,m),
for each ν=0,…,n-1, where gν(t,m) is given by (29).Lemma 2.
For anyp∈ℤ such that p>1,
(42)(H(t,m))p=(f(t,m)an-1(n-1)!mn-1)p(Sn-1(t,m))p-1×F(t,m)(a0a1⋯an-1),
where Sn-1(t,m) and F(t,m) are given by (39).Lemma 3.
For anyp,k∈ℤ such that p,k≥1,
(43)(eBt/m)k(H(t,m))p=(f(t,m)an-1(n-1)!mn-1)p(Sn-1(t))p-1×ΩF(t,m)(a0a1⋯an-1),
where Sn-1(t,m) and F(t,m) are given by (39), and
(44)Ω=((k+1)n-10⋯00(k+1)n-2⋯0⋮⋮⋱⋮00⋯1).Lemma 4.
For anyp,k∈ℤ such that p,k≥1,
(45)Hp(t,m)(eBt/m)k=(f(t,m)an-1(n-1)!mn-1)p×(Sn-1(t))p-1F(t,m)A(t,m),
where Sn-1(t,m) and F(t,m) are given by (39),
(46)A(t,m)=(A0(t,m)A1(t,m)⋯An-1(t,m)),Aμ(t,m)=∑ν=0μ‍aνkμ-νgμ-ν,
for any μ=0,…,n-1, and gν is given by (29).For the sake of simplicity, we have omitted the dependence onm and t of gν.


## 3. Solution to the Nonhomogeneous Problem

The general solution to a nonhomogeneous linear differential equation of ordern can be expressed as the sum of the general solution to the corresponding homogeneous, linear differential equation and any solution to the complete equation. The symbolic manipulation system calculates the solution to a nonperturbed differential equation with initial conditions of the form (12)
(47)x(n)+an-1x(n-1)+⋯+a1x˙+a0x=u(t),
with initial conditions
(48)x(0)=x10,…,x(n-1)(0)=xn0,
where a0,a1,…,an-1∈ℝ, x10,…,xn0∈ℝ, and
(49)u(t)=∑ν≥0‍tnνeανt(λνcos(ωνt)+μνsin(ωνt)),
being nν∈ℕ, and αν,ων,λν,μν∈ℝ. With the aid of the substitutions
(50)x1=x,x2=x˙,…,xn=x(n-1),
(12) is reduced to the system of differential equations given by
(51)X˙(t)=AX(t)+B(t),X(0)=X0,
where
(52)A=(010⋯0001⋯0⋮⋮⋮⋱⋮-a0-a1-a2⋯-an-1),B(t)=(00⋮u(t)),X(t)=(x1(t)x2(t)⋮xn(t)),X0=(x10x20⋮xn0).The computation of the solution to the constant coefficients linear part requires the calculation of the exponential of the matrixA, Φ(t)=eAt:
(53)X(t)=Φ(t)X0+Φ(t)∫0t‍exp(Aτ)B(τ)dτ.


### 3.1. Symbolic Expansion of the Solution

In the following, we give a formula for the symbolic expansion of the solution to the complete ordinary differential equation. To that end, let us express the noncommutative parenthesis(eBt/m,J(t,m))km as
(54)(eBt/m,J(t,m))km=(pr11(m,k)(t)⋯pr1n(m,k)(t)⋮⋱⋮prn1(m,k)(t)⋯prnn(m,k)(t))=(prij(m,k)(t))i,j=1n=P(m,k)(t),
where each pr1n(m,k)(t) is a polynomial of degree rij≤m(n-1) in the indeterminate t with coefficients from ℝ, and m≠0.Let us also express(55)λk(t,m)=(1an-1(ean-1t/m-1))k=∑ν=0k‍ανeβνt,
being
(56)αν=(1an-1)k(kν)(-1)ν,βν=-an-1m(k-ν).
Thus, taking into account (33), the exponential matrix of A can be arranged as
(57)eAt≃∑k=0m‍(ϕk(t)pr11(m,k)(t)⋯ϕk(t)pr1n(m,k)(t)⋮⋱⋮ϕk(t)prn1(m,k)(t)⋯ϕk(t)prnn(m,k)(t))=(∑k=0m‍ϕk(t)pr11(m,k)(t)⋯∑k=0m‍ϕk(t)pr1n(m,k)(t)⋮⋱⋮∑k=0m‍ϕk(t)prn1(m,k)(t)⋯∑k=0m‍ϕk(t)prnn(m,k)(t)),
where
(58)ϕk(t)=∑ν=0k‍ανeβνt.In order to develop a symbolic expression of the solution to the complete differential equation, let us expressu(t) as
(59)u(t)=eδt(Qa(t)cos(ωt)+Tb(t)sin(ωt)),
with a,b∈ℤ∪{0}, δ,ω∈ℝ, and Qa(t) and Tb(t) being real polynomials of degrees a and b, respectively, in the indeterminate t.The producte-AτB(τ) can be arranged as follows:
(60)e-AτB(τ)=(∫0t‍(∑k=0m‍(∑ν=0k‍ανeβντ)pr1n(m,k)(-τ))u(τ)dτ∫0t‍(∑k=0m‍(∑ν=0k‍ανeβντ)pr2n(m,k)(-τ))u(τ)dτ⋮∫0t‍(∑k=0m‍(∑ν=0k‍ανeβντ)prnn(m,k)(-τ))u(τ)dτ).
Here, the integrand of each element of the above matrix can be written in the form of a quasipolynomial. Hence, we arrive to the formulae
(61)∫0t‍(∑k=0m‍(∑ν=0k‍ανeβντ)prjn(m,k)(-τ))u(τ)dτ=∑αjn‍(ηαjnCαjn,βjn,γjn(t)+ναjnSαjn,βjn,γjn(t)),
where αjn∈ℤ+∪{0}, βjn,γjn∈ℝ, j=1,…,n, and
(62)Cm,β,γ(t)=∫‍tmeβtcos(γt)dt,Sm,β,γ(t)=∫‍tmeβtsin(γt)dt.
These functions are computed recursively through
(63)C0,β,γ(t)=γβ2+γ2eβtsin(γt)+ββ2+γ2eβtcos(γt),S0,β,γ(t)=ββ2+γ2eβtsin(γt)-γβ2+γ2eβtcos(γt),
and, for any m≥1,
(64)Cm,β,γ(t)=tmC0,β,γ(t)-mγβ2+γ2Sm-1,β,γ-mββ2+γ2Cm-1,β,γ(t),Sm,β,γ(t)=tmS0,β,γ(t)-mββ2+γ2Sm-1,β,γ+mγβ2+γ2Cm-1,β,γ(t).Thus, (53) can be expressed via a quasipolynomial and, thus, obtained through the designed symbolic system.


## 3.1. Symbolic Expansion of the Solution

In the following, we give a formula for the symbolic expansion of the solution to the complete ordinary differential equation. To that end, let us express the noncommutative parenthesis(eBt/m,J(t,m))km as
(54)(eBt/m,J(t,m))km=(pr11(m,k)(t)⋯pr1n(m,k)(t)⋮⋱⋮prn1(m,k)(t)⋯prnn(m,k)(t))=(prij(m,k)(t))i,j=1n=P(m,k)(t),
where each pr1n(m,k)(t) is a polynomial of degree rij≤m(n-1) in the indeterminate t with coefficients from ℝ, and m≠0.Let us also express(55)λk(t,m)=(1an-1(ean-1t/m-1))k=∑ν=0k‍ανeβνt,
being
(56)αν=(1an-1)k(kν)(-1)ν,βν=-an-1m(k-ν).
Thus, taking into account (33), the exponential matrix of A can be arranged as
(57)eAt≃∑k=0m‍(ϕk(t)pr11(m,k)(t)⋯ϕk(t)pr1n(m,k)(t)⋮⋱⋮ϕk(t)prn1(m,k)(t)⋯ϕk(t)prnn(m,k)(t))=(∑k=0m‍ϕk(t)pr11(m,k)(t)⋯∑k=0m‍ϕk(t)pr1n(m,k)(t)⋮⋱⋮∑k=0m‍ϕk(t)prn1(m,k)(t)⋯∑k=0m‍ϕk(t)prnn(m,k)(t)),
where
(58)ϕk(t)=∑ν=0k‍ανeβνt.In order to develop a symbolic expression of the solution to the complete differential equation, let us expressu(t) as
(59)u(t)=eδt(Qa(t)cos(ωt)+Tb(t)sin(ωt)),
with a,b∈ℤ∪{0}, δ,ω∈ℝ, and Qa(t) and Tb(t) being real polynomials of degrees a and b, respectively, in the indeterminate t.The producte-AτB(τ) can be arranged as follows:
(60)e-AτB(τ)=(∫0t‍(∑k=0m‍(∑ν=0k‍ανeβντ)pr1n(m,k)(-τ))u(τ)dτ∫0t‍(∑k=0m‍(∑ν=0k‍ανeβντ)pr2n(m,k)(-τ))u(τ)dτ⋮∫0t‍(∑k=0m‍(∑ν=0k‍ανeβντ)prnn(m,k)(-τ))u(τ)dτ).
Here, the integrand of each element of the above matrix can be written in the form of a quasipolynomial. Hence, we arrive to the formulae
(61)∫0t‍(∑k=0m‍(∑ν=0k‍ανeβντ)prjn(m,k)(-τ))u(τ)dτ=∑αjn‍(ηαjnCαjn,βjn,γjn(t)+ναjnSαjn,βjn,γjn(t)),
where αjn∈ℤ+∪{0}, βjn,γjn∈ℝ, j=1,…,n, and
(62)Cm,β,γ(t)=∫‍tmeβtcos(γt)dt,Sm,β,γ(t)=∫‍tmeβtsin(γt)dt.
These functions are computed recursively through
(63)C0,β,γ(t)=γβ2+γ2eβtsin(γt)+ββ2+γ2eβtcos(γt),S0,β,γ(t)=ββ2+γ2eβtsin(γt)-γβ2+γ2eβtcos(γt),
and, for any m≥1,
(64)Cm,β,γ(t)=tmC0,β,γ(t)-mγβ2+γ2Sm-1,β,γ-mββ2+γ2Cm-1,β,γ(t),Sm,β,γ(t)=tmS0,β,γ(t)-mββ2+γ2Sm-1,β,γ+mγβ2+γ2Cm-1,β,γ(t).Thus, (53) can be expressed via a quasipolynomial and, thus, obtained through the designed symbolic system.


## 4. Symbolic and Numeric Results

### 4.1. On the Form of the Exponential Matrix

The aim of this section is to illustrate the form that the approximation of the exponential matrix adopts. For that purpose, let us consider the differential equation given by(65)x(4)+4x⃛+3x¨+2x˙+x=0,
with initial conditions
(66)x(0)=1,x˙(0)=0,x¨(0)=0,x⃛(0)=0.
This equation is transformed into the system of differential equations given by
(67)X˙(t)=AX(t),
where A is the companion matrix
(68)A=(010000100001-1-2-3-4).
According to the technique described in Section 2, we split matrix A into B+C, where
(69)B=(0100001000010000),C=(000000000000-1-2-3-4).
Equation (36),
(70)(eBt/m+H(t,m))m=∑k=0m‍λk(t,m)(eBt/m,J(t,m))km,
allows us to compute a symbolic approximation to the solution to the differential equation. In the problem we are discussing, eBt/m, H(t,m), λ(t,m), and J(t,m) can be written as
(71)eBt/m=(1tmt2(2m)t3(6m)01tmt2(2m)001tm0001),H(t,m)=λ(t,m)J(t,m),
with
(72)λ(t,m)=14(e-4t/m-1),J(t,m)=(t3(6m3)2t3(6m3)3t3(6m3)4t3(6m3)t2(2m2)2t2(2m2)3t2(2m2)4t2(2m2)tm2tm3tm4tm1234).Taking, for instance,m=3, we get that
(73)eAt≃∑k=03‍λk(t,3)(eBt/3,J(t,3))k3=𝒜3+λ(𝒜2ℬ+ℬ𝒜2+𝒜ℬ𝒜)+λ2(𝒜ℬ2+ℬ2𝒜+ℬ𝒜ℬ)+λ3ℬ3.
For the sake of simplicity, we have introduced the following notation:
(74)𝒜=eBt/3,ℬ=J(t,3).Matrices𝒜3, 𝒜2ℬ, ℬ𝒜2, 𝒜ℬ𝒜, 𝒜ℬ2, ℬ2𝒜, ℬ𝒜ℬ, and ℬ3 are computed with the help of the designed symbolic processor, with the result
(75)𝒜3=(1tt22t3601tt22001t0001),𝒜2ℬ=(t36t33t322t33t22t23t222t2t2t3t4t1234),𝒜ℬ2=((𝒜ℬ2)1,1(𝒜ℬ2)1,2(𝒜ℬ2)1,3(𝒜ℬ2)1,4(𝒜ℬ2)2,1(𝒜ℬ2)2,2(𝒜ℬ2)2,3(𝒜ℬ2)2,4(𝒜ℬ2)3,1(𝒜ℬ2)3,2(𝒜ℬ2)3,3(𝒜ℬ2)3,4(𝒜ℬ2)4,1(𝒜ℬ2)4,2(𝒜ℬ2)4,3(𝒜ℬ2)4,4),
where
(76)(𝒜ℬ2)1,1=26561t6+4729t5+481t4+1681t3,(𝒜ℬ2)1,2=46561t6+8729t5+881t4+3281t3,(𝒜ℬ2)1,3=22187t6+4243t5+427t4+1627t3,(𝒜ℬ2)1,4=86561t6+16729t5+1681t4+6481t3,(𝒜ℬ2)2,1=1729t5+281t4+29t3+89t2,(𝒜ℬ2)2,2=2729t5+481t4+49t3+169t2,(𝒜ℬ2)2,3=1243t5+227t4+23t3+83t2,(𝒜ℬ2)2,4=4729t5+881t4+89t3+329t2,(𝒜ℬ2)3,1=1243t4+227t3+23t2+83t,(𝒜ℬ2)3,2=2243t4+427t3+43t2+163t,(𝒜ℬ2)3,3=181t4+29t3+2t2+8t,(𝒜ℬ2)3,4=4243t4+827t3+83t2+323t,(𝒜ℬ2)4,1=1162t3+19t2+t+4,(𝒜ℬ2)4,2=181t3+29t2+2t+8,(𝒜ℬ2)4,3=154t3+13t2+3t+12,(𝒜ℬ2)4,4=281t3+49t2+4t+16,𝒜ℬ𝒜=((𝒜ℬ𝒜)1,1(𝒜ℬ𝒜)1,2(𝒜ℬ𝒜)1,3(𝒜ℬ𝒜)1,4(𝒜ℬ𝒜)2,1(𝒜ℬ𝒜)2,2(𝒜ℬ𝒜)2,3(𝒜ℬ𝒜)2,4(𝒜ℬ𝒜)3,1(𝒜ℬ𝒜)3,2(𝒜ℬ𝒜)3,3(𝒜ℬ𝒜)3,4(𝒜ℬ𝒜)4,1(𝒜ℬ𝒜)4,2(𝒜ℬ𝒜)4,3(𝒜ℬ𝒜)4,4),
where
(77)(𝒜ℬ𝒜)1,1=481t3,(𝒜ℬ𝒜)1,2=4243t4+881t3,(𝒜ℬ𝒜)1,3=2729t5+8243t4+427t3,(𝒜ℬ𝒜)1,4=26561t6+4729t5+481t4+1681t3,(𝒜ℬ𝒜)2,1=29t2,(𝒜ℬ𝒜)2,2=227t3+49t2,(𝒜ℬ𝒜)2,3=181t4+427t3+23t2,(𝒜ℬ𝒜)2,4=1729t5+281t4+29t3+89t2,(𝒜ℬ𝒜)3,1=23t,(𝒜ℬ𝒜)3,2=29t2+43t,(𝒜ℬ𝒜)3,3=127t3+49t2+2t,(𝒜ℬ𝒜)3,4=1243t4+227t3+23t2+83t,(𝒜ℬ𝒜)4,1=1,(𝒜ℬ𝒜)4,2=13t+2,(𝒜ℬ𝒜)4,3=118t2+23t+3,(𝒜ℬ𝒜)4,4=1162t3+19t2+t+4,ℬ3=((ℬ3)1,1(ℬ3)1,2(ℬ3)1,3(ℬ3)1,4(ℬ3)2,1(ℬ3)2,2(ℬ3)2,3(ℬ3)2,4(ℬ3)3,1(ℬ3)3,2(ℬ3)3,3(ℬ3)3,4(ℬ3)4,1(ℬ3)4,2(ℬ3)4,3(ℬ3)4,4),
where
(78)(ℬ3)1,1=14251528t9+1118098t8+16561t7+116561t6+171458t5+481t4+881t3,(ℬ3)1,2=12125764t9+159049t8+26561t7+226561t6+17729t5+881t4+1681t3,(ℬ3)1,3=11417176t9+139366t8+12187t7+112187t6+17486t5+427t4+827t3,(ℬ3)1,4=11062882t9+259049t8+46561t7+446561t6+34729t5+1681t4+3281t3,(ℬ3)2,1=1472392t8+113122t7+1729t6+11729t5+17162t4+49t3+89t2,(ℬ3)2,2=1236196t8+16561t7+2729t6+22729t5+1781t4+89t3+169t2,(ℬ3)2,3=1157464t8+14374t7+1243t6+11243t5+1754t4+43t3+83t2,(ℬ3)2,4=1118098t8+26561t7+4729t6+44729t5+3481t4+169t3+329t2,(ℬ3)3,1=178732t7+12187t6+2243t5+22243t4+1727t3+83t2+163t,(ℬ3)3,2=139366t7+22187t6+4243t5+44243t4+3427t3+163t2+323t,(ℬ3)3,3=126244t7+1729t6+281t5+2281t4+179t3+8t2+16t,(ℬ3)3,4=119683t7+42187t6+8243t5+88243t4+6827t3+323t2+643t,(ℬ3)4,1=126244t6+1729t5+281t4+2281t3+179t2+8t+16,(ℬ3)4,2=113122t6+2729t5+481t4+4481t3+349t2+16t+32,(ℬ3)4,3=18748t6+1243t5+227t4+2227t3+173t2+24t+48,(ℬ3)4,4=16561t6+4729t5+881t4+8881t3+689t2+32t+64,ℬ2𝒜=((ℬ2𝒜)1,1(ℬ2𝒜)1,2(ℬ2𝒜)1,3(ℬ2𝒜)1,4(ℬ2𝒜)2,1(ℬ2𝒜)2,2(ℬ2𝒜)2,3(ℬ2𝒜)2,4(ℬ2𝒜)3,1(ℬ2𝒜)3,2(ℬ2𝒜)3,3(ℬ2𝒜)3,4(ℬ2𝒜)4,1(ℬ2𝒜)4,2(ℬ2𝒜)4,3(ℬ2𝒜)4,4),
where
(79)(ℬ2𝒜)1,1=126244t6+11458t5+1162t4+281t3,(ℬ2𝒜)1,2=178732t7+26561t6+51458t5+5243t4+481t3,(ℬ2𝒜)1,3=1472392t8+578732t7+22187t6+111458t5+17486t4+227t3,(ℬ2𝒜)1,4=14251528t9+1118098t8+16561t7+116561t6+171458t5+481t4+881t3,(ℬ2𝒜)2,1=12916t5+1162t4+118t3+29t2,(ℬ2𝒜)2,2=18748t6+2729t5+5162t4+527t3+49t2,(ℬ2𝒜)2,3=152488t7+58748t6+2243t5+11162t4+1754t3+23t2,(ℬ2𝒜)2,4=1472392t8+113122t7+1729t6+11729t5+17162t4+49t3+89t2,(ℬ2𝒜)3,1=1486t4+127t3+13t2+43t,(ℬ2𝒜)3,2=11458t5+4243t4+527t3+109t2+83t,(ℬ2𝒜)3,3=18748t6+51458t5+481t4+1127t3+179t2+4t,(ℬ2𝒜)3,4=178732t7+12187t6+2243t5+22243t4+1727t3+83t2+163t,(ℬ2𝒜)4,1=1162t3+19t2+t+4,(ℬ2𝒜)4,2=1486t4+481t3+59t2+103t+8,(ℬ2𝒜)4,3=12916t5+5486t4+427t3+119t2+173t+12,(ℬ2𝒜)4,4=126244t6+1729t5+281t4+2281t3+179t2+8t+16,ℬ𝒜2=((ℬ𝒜2)1,1(ℬ𝒜2)1,2(ℬ𝒜2)1,3(ℬ𝒜2)1,4(ℬ𝒜2)2,1(ℬ𝒜2)2,2(ℬ𝒜2)2,3(ℬ𝒜2)2,4(ℬ𝒜2)3,1(ℬ𝒜2)3,2(ℬ𝒜2)3,3(ℬ𝒜2)3,4(ℬ𝒜2)4,1(ℬ𝒜2)4,2(ℬ𝒜2)4,3(ℬ𝒜2)4,4),
where
(80)(ℬ𝒜2)1,1=1162t3,(ℬ𝒜2)1,2=1243t4+181t3,(ℬ𝒜2)1,3=1729t5+2243t4+154t3,(ℬ𝒜2)1,4=26561t6+2729t5+181t4+281t3,(ℬ𝒜2)2,1=118t2,(ℬ𝒜2)2,2=127t3+19t2,(ℬ𝒜2)2,3=181t4+227t3+16t2,(ℬ𝒜2)2,4=2729t5+281t4+19t3+29t2,(ℬ𝒜2)3,1=13t,(ℬ𝒜2)3,2=29t2+23t,(ℬ𝒜2)3,3=227t3+49t2+t,(ℬ𝒜2)3,4=4243t4+427t3+23t2+43t,(ℬ𝒜2)4,1=1,(ℬ𝒜2)4,2=23t+2,(ℬ𝒜2)4,3=29t2+43t+3,(ℬ𝒜2)4,4=481t3+49t2+2t+4,ℬ𝒜ℬ=((ℬ𝒜ℬ)1,1(ℬ𝒜ℬ)1,2(ℬ𝒜ℬ)1,3(ℬ𝒜ℬ)1,4(ℬ𝒜ℬ)2,1(ℬ𝒜ℬ)2,2(ℬ𝒜ℬ)2,3(ℬ𝒜ℬ)2,4(ℬ𝒜ℬ)3,1(ℬ𝒜ℬ)3,2(ℬ𝒜ℬ)3,3(ℬ𝒜ℬ)3,4(ℬ𝒜ℬ)4,1(ℬ𝒜ℬ)4,2(ℬ𝒜ℬ)4,3(ℬ𝒜ℬ)4,4),
where
(81)(ℬ𝒜ℬ)1,1=26561t6+2729t5+181t4+281t3,(ℬ𝒜ℬ)1,1=46561t6+4729t5+281t4+481t3,(ℬ𝒜ℬ)1,1=22187t6+2243t5+127t4+227t3,(ℬ𝒜ℬ)1,1=86561t6+8729t5+481t4+881t3,(ℬ𝒜ℬ)1,1=2729t5+281t4+19t3+29t2,(ℬ𝒜ℬ)1,1=4729t5+481t4+29t3+49t2,(ℬ𝒜ℬ)1,1=2243t5+227t4+13t3+23t2,(ℬ𝒜ℬ)1,1=8729t5+881t4+49t3+89t2,(ℬ𝒜ℬ)1,1=4243t4+427t3+23t2+43t,(ℬ𝒜ℬ)1,1=8243t4+827t3+43t2+83t,(ℬ𝒜ℬ)1,1=481t4+49t3+2t2+4t,(ℬ𝒜ℬ)1,1=16243t4+1627t3+83t2+163t,(ℬ𝒜ℬ)1,1=481t3+49t2+2t+4,(ℬ𝒜ℬ)1,1=881t3+89t2+4t+8,(ℬ𝒜ℬ)1,1=427t3+43t2+6t+12,(ℬ𝒜ℬ)1,1=1681t3+169t2+8t+16.
So, the noncommutative parentheses are given by
(82)(𝒜,ℬ)03=𝒜3,(𝒜,ℬ)33=ℬ3,(𝒜,ℬ)13=(𝒜2ℬ+ℬ𝒜2+𝒜ℬ𝒜)=(α1,1α1,2α1,3α1,4α2,1α2,2α2,3α2,4α3,1α3,2α3,3α3,4α4,1α4,2α4,3α4,4),
where
(83)α1,1=29t3,α1,2=49t3+5243t4,α1,3=23t3+1243t5+10243t4,α1,4=89t3+46561t6+2243t5+581t4,α2,1=79t2,α2,2=149t2+19t3,α2,3=73t2+281t4+29t3,α2,4=289t2+1243t5+481t4+13t3,α3,1=2t,α3,2=4t+49t2,α3,3=6t+19t3+89t2,α3,4=8t+5243t4+29t3+43t2,α4,1=3,α4,2=6+t,α4,3=9+518t2+2t,α4,4=12+118t3+59t2+3t,(𝒜,ℬ)23=(𝒜ℬ2+ℬ2𝒜+ℬ𝒜ℬ)=(β1,1β1,2β1,3β1,4β2,1β2,2β2,3β2,4β3,1β3,2β3,3β3,4β4,1β4,2β4,3β4,4),
being
(84)β1,1=1726244t6+131458t5+11162t4+2081t3,β1,2=106561t6+291458t5+35243t4+4081t3+178732t7,β1,3=2729t6+471458t5+107486t4+2027t3+1472392t8+578732t7,β1,4=1243t6+651458t5+827t4+8081t3+14251528t9+1118098t8+16561t7,β2,1=132916t5+118t4+718t3+43t2,β2,2=8729t5+754t4+2327t3+83t2+18748t6,β2,3=5243t5+35162t4+7154t3+4t2+152488t7+58748t6,β2,4=23729t5+49162t4+169t3+163t2+1472392t8+113122t7+1729t6,β3,1=11486t4+727t3+53t2+163t,β3,2=14243t4+1727t3+349t2+323t+11458t5,β3,3=19t4+2927t3+539t2+16t+18748t6+51458t5,β3,4=1481t4+4127t3+8t2+643t+178732t7+12187t6+2243t5,β4,1=581t3+23t2+4t+12,β4,2=1381t3+53t2+283t+24+1486t4,β4,3=1754t3+269t2+443t+36+12916t5+5486t4,β4,4=4081t3+379t2+20t+48+126244t6+1729t5+281t4.In Figure1, we show a comparison of the solution computed through the symbolic method presented in this paper, and the numeric solution to the problem calculated by a Runge-Kutta fourth order method with a step of h=0.001. Let us stress that the difference between both solutions at any time is smaller than 10-7.Figure 1
Comparison of the solution computed through the symbolic method presented in this paper, and the numeric solution to the problem calculated by a Runge-Kutta fourth order method with a step ofh=0.001.


### 4.2. Description of the Program

In order to describe the algebraic processor, let us introduce the following test problem:(85)x¨+x˙+x=sint,
with initial conditions
(86)x(0)=0,x˙(0)=1.
This equation describes the (small) angular position in radians x(t) of a forced damped pendulum with a periodic driving force. In this equation, x¨(t) represents the inertia, x˙(t) represents the friction, and sin(t) represents a sinusoidal driving torque applied at the pivot of the pendulum. The exact solution to this problem is given by
(87)x(t)=-cost+e-t/2cos(32t)+3e-t/2sin(32t).
The program proceeds as follows. The first window (Figure 2) allows us introducing the definition of variables, the order of the ODE and its coefficients, and the function of the nonhomogeneous term u(t).Figure 2
The first window shown by the symbolic processor allows us to introduce the parameters which define the differential equation to be solved, including the initial conditions.The next window visualizes the expression of the companion matrix related to the ODE (Figure3).Figure 3
MatrixA.Then, matricesB, exp(Bt/m), C, and exp(Ct/m), which are used in the calculation of the exponential matrix, are computed by using a value of m satisfying (23). These matrices are presented in Figures 4, 5, 6, and 7.Figure 4
MatrixB.Figure 5
Matrixexp(Bt/m).Figure 6
MatrixC.Figure 7
Matrixexp(Ct/m).In Figures8 and 9, we evaluate the exponential matrix eAt in t=1 and for m=1000. The algebraic processor enables us to change these two parameters. We should emphasize that the symbolic solution to the problem is dependent on the value of parameter m. We have computed this solution for different values of this parameter (m=10,100,1000, and 10000). In Figures 10 and 11, we compare the real solution with symbolic solutions x(t) for these values of m. We observe that the symbolic solution gets closer to the real solution as the value of m is increased.Figure 8
Parameters related to the accuracy of the solution to the differential equation.Figure 9
Matrixexp(A).Figure 10
Comparison of real solution (black) and symbolic solutions form=10 (light blue line) and m=100 (yellow line).Figure 11
Comparison of real solution (black) and symbolic solutions form=1000 (dark blue line) and m=10000 (green line).Table1 shows the numerical evolution of the symbolic approximation xS(t), computed for m=1000. These values have been obtained by substituting variable t in the symbolic solution by values 0.1,0.2,…,1.Table 1
Numerical evolution of the symbolic approximationxS(t), for m=1000.
Step
t
x
(
t
)
1
0.1
0.0951665336810
2
0.2
0.1813305224920
3
0.3
0.2594801299955
4
0.4
0.3305856479131
5
0.5
0.3955920028686
6
0.6
0.4554122210115
7
0.7
0.5109218189711
8
0.8
0.5629540836788
9
0.9
0.6122961986812
10
1.0
0.6596861707192In Figure12, we compare the symbolic approximation, computed via the symbolic processor, with the exact solution (87). We also compare these two functions with a numeric solution computed through a Runge-Kutta 4th order method with step h=0.1. We see that the solution symbolically computed fits better to the exact solution to than the described numerical solution. In Table 2, we show the error on the solution computed using both methods. We observe that error obtained in the symbolic approximation is quite smaller that the numeric one. This fact can be better observed in Figure 13, where we show the difference in absolute value between the exact solution and the symbolic and numeric solutions. Nevertheless, it is not the goal of this paper to develop a numeric tool. The symbolic technique we have developed provides an analytical solution which can be used as a kernel to apply perturbation methods to compute the solution to a perturbed differential equation depending on a small parameter.Table 2
Numerical error in the symbolic (xS(t)) and numeric (xR(t)) approximations, with respect to the exact solution x(t).
t
|
x
S
(
t
)
-
x
(
t
)
|
|
x
R
(
t
)
-
x
(
t
)
|
0.1
0.0000047853239
0.0000001672031
0.4
0.0000712453497
0.0000358751267
0.7
0.0002074165163
0.0003544605868
1.0
0.0003710415717
0.0014307873617
1.3
0.0004547161744
0.0038253117836
1.6
0.0003284030320
0.0080436899562
1.9
0.0000622406880
0.0144197157174
2.2
0.0006216534547
0.0230329551851
2.5
0.0011277539714
0.0336675456541
2.8
0.0013572764337
0.0458134712942
3.1
0.0012347982225
0.0587076431208
3.4
0.0008953272672
0.0714078077594
3.7
0.0006002231245
0.0828896894497
4.0
0.0005492334199
0.0921556872924
4.3
0.0007154938930
0.0983431460014
4.6
0.0008300921248
0.1008204330678
4.9
0.0005487250652
0.0992608551151
5.2
0.0002963137375
0.0936868279861
5.5
0.0015298326576
0.0844798593631
5.8
0.0026838757810
0.0723553936143
6.1
0.0032571351779
0.0583046260627
6.4
0.0030289419613
0.0435091879952
6.7
0.0022189180928
0.0292364988851
7.0
0.0013691340100
0.0167257096706
7.3
0.0010043530807
0.0070750555622
7.6
0.0012821439144
0.0011409176473
7.9
0.0018689726696
0.0005417339524
8.2
0.0021436005756
0.0021897330323
8.5
0.0016148602301
0.0091082291225
8.8
0.0002869277587
0.0196147078451
9.1
0.0012796703448
0.0327897203523
9.4
0.0022771057968
0.0474740017777
9.7
0.0021518593816
0.0623710756305
10.0
0.0009803998394
0.0761625070468Figure 12
Exact solution (black), symbolic solution (blue), and numeric solution (green).Figure 13
Error in the symbolic (black) and numeric (blue) approximations.


## 4.1. On the Form of the Exponential Matrix

The aim of this section is to illustrate the form that the approximation of the exponential matrix adopts. For that purpose, let us consider the differential equation given by(65)x(4)+4x⃛+3x¨+2x˙+x=0,
with initial conditions
(66)x(0)=1,x˙(0)=0,x¨(0)=0,x⃛(0)=0.
This equation is transformed into the system of differential equations given by
(67)X˙(t)=AX(t),
where A is the companion matrix
(68)A=(010000100001-1-2-3-4).
According to the technique described in Section 2, we split matrix A into B+C, where
(69)B=(0100001000010000),C=(000000000000-1-2-3-4).
Equation (36),
(70)(eBt/m+H(t,m))m=∑k=0m‍λk(t,m)(eBt/m,J(t,m))km,
allows us to compute a symbolic approximation to the solution to the differential equation. In the problem we are discussing, eBt/m, H(t,m), λ(t,m), and J(t,m) can be written as
(71)eBt/m=(1tmt2(2m)t3(6m)01tmt2(2m)001tm0001),H(t,m)=λ(t,m)J(t,m),
with
(72)λ(t,m)=14(e-4t/m-1),J(t,m)=(t3(6m3)2t3(6m3)3t3(6m3)4t3(6m3)t2(2m2)2t2(2m2)3t2(2m2)4t2(2m2)tm2tm3tm4tm1234).Taking, for instance,m=3, we get that
(73)eAt≃∑k=03‍λk(t,3)(eBt/3,J(t,3))k3=𝒜3+λ(𝒜2ℬ+ℬ𝒜2+𝒜ℬ𝒜)+λ2(𝒜ℬ2+ℬ2𝒜+ℬ𝒜ℬ)+λ3ℬ3.
For the sake of simplicity, we have introduced the following notation:
(74)𝒜=eBt/3,ℬ=J(t,3).Matrices𝒜3, 𝒜2ℬ, ℬ𝒜2, 𝒜ℬ𝒜, 𝒜ℬ2, ℬ2𝒜, ℬ𝒜ℬ, and ℬ3 are computed with the help of the designed symbolic processor, with the result
(75)𝒜3=(1tt22t3601tt22001t0001),𝒜2ℬ=(t36t33t322t33t22t23t222t2t2t3t4t1234),𝒜ℬ2=((𝒜ℬ2)1,1(𝒜ℬ2)1,2(𝒜ℬ2)1,3(𝒜ℬ2)1,4(𝒜ℬ2)2,1(𝒜ℬ2)2,2(𝒜ℬ2)2,3(𝒜ℬ2)2,4(𝒜ℬ2)3,1(𝒜ℬ2)3,2(𝒜ℬ2)3,3(𝒜ℬ2)3,4(𝒜ℬ2)4,1(𝒜ℬ2)4,2(𝒜ℬ2)4,3(𝒜ℬ2)4,4),
where
(76)(𝒜ℬ2)1,1=26561t6+4729t5+481t4+1681t3,(𝒜ℬ2)1,2=46561t6+8729t5+881t4+3281t3,(𝒜ℬ2)1,3=22187t6+4243t5+427t4+1627t3,(𝒜ℬ2)1,4=86561t6+16729t5+1681t4+6481t3,(𝒜ℬ2)2,1=1729t5+281t4+29t3+89t2,(𝒜ℬ2)2,2=2729t5+481t4+49t3+169t2,(𝒜ℬ2)2,3=1243t5+227t4+23t3+83t2,(𝒜ℬ2)2,4=4729t5+881t4+89t3+329t2,(𝒜ℬ2)3,1=1243t4+227t3+23t2+83t,(𝒜ℬ2)3,2=2243t4+427t3+43t2+163t,(𝒜ℬ2)3,3=181t4+29t3+2t2+8t,(𝒜ℬ2)3,4=4243t4+827t3+83t2+323t,(𝒜ℬ2)4,1=1162t3+19t2+t+4,(𝒜ℬ2)4,2=181t3+29t2+2t+8,(𝒜ℬ2)4,3=154t3+13t2+3t+12,(𝒜ℬ2)4,4=281t3+49t2+4t+16,𝒜ℬ𝒜=((𝒜ℬ𝒜)1,1(𝒜ℬ𝒜)1,2(𝒜ℬ𝒜)1,3(𝒜ℬ𝒜)1,4(𝒜ℬ𝒜)2,1(𝒜ℬ𝒜)2,2(𝒜ℬ𝒜)2,3(𝒜ℬ𝒜)2,4(𝒜ℬ𝒜)3,1(𝒜ℬ𝒜)3,2(𝒜ℬ𝒜)3,3(𝒜ℬ𝒜)3,4(𝒜ℬ𝒜)4,1(𝒜ℬ𝒜)4,2(𝒜ℬ𝒜)4,3(𝒜ℬ𝒜)4,4),
where
(77)(𝒜ℬ𝒜)1,1=481t3,(𝒜ℬ𝒜)1,2=4243t4+881t3,(𝒜ℬ𝒜)1,3=2729t5+8243t4+427t3,(𝒜ℬ𝒜)1,4=26561t6+4729t5+481t4+1681t3,(𝒜ℬ𝒜)2,1=29t2,(𝒜ℬ𝒜)2,2=227t3+49t2,(𝒜ℬ𝒜)2,3=181t4+427t3+23t2,(𝒜ℬ𝒜)2,4=1729t5+281t4+29t3+89t2,(𝒜ℬ𝒜)3,1=23t,(𝒜ℬ𝒜)3,2=29t2+43t,(𝒜ℬ𝒜)3,3=127t3+49t2+2t,(𝒜ℬ𝒜)3,4=1243t4+227t3+23t2+83t,(𝒜ℬ𝒜)4,1=1,(𝒜ℬ𝒜)4,2=13t+2,(𝒜ℬ𝒜)4,3=118t2+23t+3,(𝒜ℬ𝒜)4,4=1162t3+19t2+t+4,ℬ3=((ℬ3)1,1(ℬ3)1,2(ℬ3)1,3(ℬ3)1,4(ℬ3)2,1(ℬ3)2,2(ℬ3)2,3(ℬ3)2,4(ℬ3)3,1(ℬ3)3,2(ℬ3)3,3(ℬ3)3,4(ℬ3)4,1(ℬ3)4,2(ℬ3)4,3(ℬ3)4,4),
where
(78)(ℬ3)1,1=14251528t9+1118098t8+16561t7+116561t6+171458t5+481t4+881t3,(ℬ3)1,2=12125764t9+159049t8+26561t7+226561t6+17729t5+881t4+1681t3,(ℬ3)1,3=11417176t9+139366t8+12187t7+112187t6+17486t5+427t4+827t3,(ℬ3)1,4=11062882t9+259049t8+46561t7+446561t6+34729t5+1681t4+3281t3,(ℬ3)2,1=1472392t8+113122t7+1729t6+11729t5+17162t4+49t3+89t2,(ℬ3)2,2=1236196t8+16561t7+2729t6+22729t5+1781t4+89t3+169t2,(ℬ3)2,3=1157464t8+14374t7+1243t6+11243t5+1754t4+43t3+83t2,(ℬ3)2,4=1118098t8+26561t7+4729t6+44729t5+3481t4+169t3+329t2,(ℬ3)3,1=178732t7+12187t6+2243t5+22243t4+1727t3+83t2+163t,(ℬ3)3,2=139366t7+22187t6+4243t5+44243t4+3427t3+163t2+323t,(ℬ3)3,3=126244t7+1729t6+281t5+2281t4+179t3+8t2+16t,(ℬ3)3,4=119683t7+42187t6+8243t5+88243t4+6827t3+323t2+643t,(ℬ3)4,1=126244t6+1729t5+281t4+2281t3+179t2+8t+16,(ℬ3)4,2=113122t6+2729t5+481t4+4481t3+349t2+16t+32,(ℬ3)4,3=18748t6+1243t5+227t4+2227t3+173t2+24t+48,(ℬ3)4,4=16561t6+4729t5+881t4+8881t3+689t2+32t+64,ℬ2𝒜=((ℬ2𝒜)1,1(ℬ2𝒜)1,2(ℬ2𝒜)1,3(ℬ2𝒜)1,4(ℬ2𝒜)2,1(ℬ2𝒜)2,2(ℬ2𝒜)2,3(ℬ2𝒜)2,4(ℬ2𝒜)3,1(ℬ2𝒜)3,2(ℬ2𝒜)3,3(ℬ2𝒜)3,4(ℬ2𝒜)4,1(ℬ2𝒜)4,2(ℬ2𝒜)4,3(ℬ2𝒜)4,4),
where
(79)(ℬ2𝒜)1,1=126244t6+11458t5+1162t4+281t3,(ℬ2𝒜)1,2=178732t7+26561t6+51458t5+5243t4+481t3,(ℬ2𝒜)1,3=1472392t8+578732t7+22187t6+111458t5+17486t4+227t3,(ℬ2𝒜)1,4=14251528t9+1118098t8+16561t7+116561t6+171458t5+481t4+881t3,(ℬ2𝒜)2,1=12916t5+1162t4+118t3+29t2,(ℬ2𝒜)2,2=18748t6+2729t5+5162t4+527t3+49t2,(ℬ2𝒜)2,3=152488t7+58748t6+2243t5+11162t4+1754t3+23t2,(ℬ2𝒜)2,4=1472392t8+113122t7+1729t6+11729t5+17162t4+49t3+89t2,(ℬ2𝒜)3,1=1486t4+127t3+13t2+43t,(ℬ2𝒜)3,2=11458t5+4243t4+527t3+109t2+83t,(ℬ2𝒜)3,3=18748t6+51458t5+481t4+1127t3+179t2+4t,(ℬ2𝒜)3,4=178732t7+12187t6+2243t5+22243t4+1727t3+83t2+163t,(ℬ2𝒜)4,1=1162t3+19t2+t+4,(ℬ2𝒜)4,2=1486t4+481t3+59t2+103t+8,(ℬ2𝒜)4,3=12916t5+5486t4+427t3+119t2+173t+12,(ℬ2𝒜)4,4=126244t6+1729t5+281t4+2281t3+179t2+8t+16,ℬ𝒜2=((ℬ𝒜2)1,1(ℬ𝒜2)1,2(ℬ𝒜2)1,3(ℬ𝒜2)1,4(ℬ𝒜2)2,1(ℬ𝒜2)2,2(ℬ𝒜2)2,3(ℬ𝒜2)2,4(ℬ𝒜2)3,1(ℬ𝒜2)3,2(ℬ𝒜2)3,3(ℬ𝒜2)3,4(ℬ𝒜2)4,1(ℬ𝒜2)4,2(ℬ𝒜2)4,3(ℬ𝒜2)4,4),
where
(80)(ℬ𝒜2)1,1=1162t3,(ℬ𝒜2)1,2=1243t4+181t3,(ℬ𝒜2)1,3=1729t5+2243t4+154t3,(ℬ𝒜2)1,4=26561t6+2729t5+181t4+281t3,(ℬ𝒜2)2,1=118t2,(ℬ𝒜2)2,2=127t3+19t2,(ℬ𝒜2)2,3=181t4+227t3+16t2,(ℬ𝒜2)2,4=2729t5+281t4+19t3+29t2,(ℬ𝒜2)3,1=13t,(ℬ𝒜2)3,2=29t2+23t,(ℬ𝒜2)3,3=227t3+49t2+t,(ℬ𝒜2)3,4=4243t4+427t3+23t2+43t,(ℬ𝒜2)4,1=1,(ℬ𝒜2)4,2=23t+2,(ℬ𝒜2)4,3=29t2+43t+3,(ℬ𝒜2)4,4=481t3+49t2+2t+4,ℬ𝒜ℬ=((ℬ𝒜ℬ)1,1(ℬ𝒜ℬ)1,2(ℬ𝒜ℬ)1,3(ℬ𝒜ℬ)1,4(ℬ𝒜ℬ)2,1(ℬ𝒜ℬ)2,2(ℬ𝒜ℬ)2,3(ℬ𝒜ℬ)2,4(ℬ𝒜ℬ)3,1(ℬ𝒜ℬ)3,2(ℬ𝒜ℬ)3,3(ℬ𝒜ℬ)3,4(ℬ𝒜ℬ)4,1(ℬ𝒜ℬ)4,2(ℬ𝒜ℬ)4,3(ℬ𝒜ℬ)4,4),
where
(81)(ℬ𝒜ℬ)1,1=26561t6+2729t5+181t4+281t3,(ℬ𝒜ℬ)1,1=46561t6+4729t5+281t4+481t3,(ℬ𝒜ℬ)1,1=22187t6+2243t5+127t4+227t3,(ℬ𝒜ℬ)1,1=86561t6+8729t5+481t4+881t3,(ℬ𝒜ℬ)1,1=2729t5+281t4+19t3+29t2,(ℬ𝒜ℬ)1,1=4729t5+481t4+29t3+49t2,(ℬ𝒜ℬ)1,1=2243t5+227t4+13t3+23t2,(ℬ𝒜ℬ)1,1=8729t5+881t4+49t3+89t2,(ℬ𝒜ℬ)1,1=4243t4+427t3+23t2+43t,(ℬ𝒜ℬ)1,1=8243t4+827t3+43t2+83t,(ℬ𝒜ℬ)1,1=481t4+49t3+2t2+4t,(ℬ𝒜ℬ)1,1=16243t4+1627t3+83t2+163t,(ℬ𝒜ℬ)1,1=481t3+49t2+2t+4,(ℬ𝒜ℬ)1,1=881t3+89t2+4t+8,(ℬ𝒜ℬ)1,1=427t3+43t2+6t+12,(ℬ𝒜ℬ)1,1=1681t3+169t2+8t+16.
So, the noncommutative parentheses are given by
(82)(𝒜,ℬ)03=𝒜3,(𝒜,ℬ)33=ℬ3,(𝒜,ℬ)13=(𝒜2ℬ+ℬ𝒜2+𝒜ℬ𝒜)=(α1,1α1,2α1,3α1,4α2,1α2,2α2,3α2,4α3,1α3,2α3,3α3,4α4,1α4,2α4,3α4,4),
where
(83)α1,1=29t3,α1,2=49t3+5243t4,α1,3=23t3+1243t5+10243t4,α1,4=89t3+46561t6+2243t5+581t4,α2,1=79t2,α2,2=149t2+19t3,α2,3=73t2+281t4+29t3,α2,4=289t2+1243t5+481t4+13t3,α3,1=2t,α3,2=4t+49t2,α3,3=6t+19t3+89t2,α3,4=8t+5243t4+29t3+43t2,α4,1=3,α4,2=6+t,α4,3=9+518t2+2t,α4,4=12+118t3+59t2+3t,(𝒜,ℬ)23=(𝒜ℬ2+ℬ2𝒜+ℬ𝒜ℬ)=(β1,1β1,2β1,3β1,4β2,1β2,2β2,3β2,4β3,1β3,2β3,3β3,4β4,1β4,2β4,3β4,4),
being
(84)β1,1=1726244t6+131458t5+11162t4+2081t3,β1,2=106561t6+291458t5+35243t4+4081t3+178732t7,β1,3=2729t6+471458t5+107486t4+2027t3+1472392t8+578732t7,β1,4=1243t6+651458t5+827t4+8081t3+14251528t9+1118098t8+16561t7,β2,1=132916t5+118t4+718t3+43t2,β2,2=8729t5+754t4+2327t3+83t2+18748t6,β2,3=5243t5+35162t4+7154t3+4t2+152488t7+58748t6,β2,4=23729t5+49162t4+169t3+163t2+1472392t8+113122t7+1729t6,β3,1=11486t4+727t3+53t2+163t,β3,2=14243t4+1727t3+349t2+323t+11458t5,β3,3=19t4+2927t3+539t2+16t+18748t6+51458t5,β3,4=1481t4+4127t3+8t2+643t+178732t7+12187t6+2243t5,β4,1=581t3+23t2+4t+12,β4,2=1381t3+53t2+283t+24+1486t4,β4,3=1754t3+269t2+443t+36+12916t5+5486t4,β4,4=4081t3+379t2+20t+48+126244t6+1729t5+281t4.In Figure1, we show a comparison of the solution computed through the symbolic method presented in this paper, and the numeric solution to the problem calculated by a Runge-Kutta fourth order method with a step of h=0.001. Let us stress that the difference between both solutions at any time is smaller than 10-7.Figure 1
Comparison of the solution computed through the symbolic method presented in this paper, and the numeric solution to the problem calculated by a Runge-Kutta fourth order method with a step ofh=0.001.


## 4.2. Description of the Program

In order to describe the algebraic processor, let us introduce the following test problem:(85)x¨+x˙+x=sint,
with initial conditions
(86)x(0)=0,x˙(0)=1.
This equation describes the (small) angular position in radians x(t) of a forced damped pendulum with a periodic driving force. In this equation, x¨(t) represents the inertia, x˙(t) represents the friction, and sin(t) represents a sinusoidal driving torque applied at the pivot of the pendulum. The exact solution to this problem is given by
(87)x(t)=-cost+e-t/2cos(32t)+3e-t/2sin(32t).
The program proceeds as follows. The first window (Figure 2) allows us introducing the definition of variables, the order of the ODE and its coefficients, and the function of the nonhomogeneous term u(t).Figure 2
The first window shown by the symbolic processor allows us to introduce the parameters which define the differential equation to be solved, including the initial conditions.The next window visualizes the expression of the companion matrix related to the ODE (Figure3).Figure 3
MatrixA.Then, matricesB, exp(Bt/m), C, and exp(Ct/m), which are used in the calculation of the exponential matrix, are computed by using a value of m satisfying (23). These matrices are presented in Figures 4, 5, 6, and 7.Figure 4
MatrixB.Figure 5
Matrixexp(Bt/m).Figure 6
MatrixC.Figure 7
Matrixexp(Ct/m).In Figures8 and 9, we evaluate the exponential matrix eAt in t=1 and for m=1000. The algebraic processor enables us to change these two parameters. We should emphasize that the symbolic solution to the problem is dependent on the value of parameter m. We have computed this solution for different values of this parameter (m=10,100,1000, and 10000). In Figures 10 and 11, we compare the real solution with symbolic solutions x(t) for these values of m. We observe that the symbolic solution gets closer to the real solution as the value of m is increased.Figure 8
Parameters related to the accuracy of the solution to the differential equation.Figure 9
Matrixexp(A).Figure 10
Comparison of real solution (black) and symbolic solutions form=10 (light blue line) and m=100 (yellow line).Figure 11
Comparison of real solution (black) and symbolic solutions form=1000 (dark blue line) and m=10000 (green line).Table1 shows the numerical evolution of the symbolic approximation xS(t), computed for m=1000. These values have been obtained by substituting variable t in the symbolic solution by values 0.1,0.2,…,1.Table 1
Numerical evolution of the symbolic approximationxS(t), for m=1000.
Step
t
x
(
t
)
1
0.1
0.0951665336810
2
0.2
0.1813305224920
3
0.3
0.2594801299955
4
0.4
0.3305856479131
5
0.5
0.3955920028686
6
0.6
0.4554122210115
7
0.7
0.5109218189711
8
0.8
0.5629540836788
9
0.9
0.6122961986812
10
1.0
0.6596861707192In Figure12, we compare the symbolic approximation, computed via the symbolic processor, with the exact solution (87). We also compare these two functions with a numeric solution computed through a Runge-Kutta 4th order method with step h=0.1. We see that the solution symbolically computed fits better to the exact solution to than the described numerical solution. In Table 2, we show the error on the solution computed using both methods. We observe that error obtained in the symbolic approximation is quite smaller that the numeric one. This fact can be better observed in Figure 13, where we show the difference in absolute value between the exact solution and the symbolic and numeric solutions. Nevertheless, it is not the goal of this paper to develop a numeric tool. The symbolic technique we have developed provides an analytical solution which can be used as a kernel to apply perturbation methods to compute the solution to a perturbed differential equation depending on a small parameter.Table 2
Numerical error in the symbolic (xS(t)) and numeric (xR(t)) approximations, with respect to the exact solution x(t).
t
|
x
S
(
t
)
-
x
(
t
)
|
|
x
R
(
t
)
-
x
(
t
)
|
0.1
0.0000047853239
0.0000001672031
0.4
0.0000712453497
0.0000358751267
0.7
0.0002074165163
0.0003544605868
1.0
0.0003710415717
0.0014307873617
1.3
0.0004547161744
0.0038253117836
1.6
0.0003284030320
0.0080436899562
1.9
0.0000622406880
0.0144197157174
2.2
0.0006216534547
0.0230329551851
2.5
0.0011277539714
0.0336675456541
2.8
0.0013572764337
0.0458134712942
3.1
0.0012347982225
0.0587076431208
3.4
0.0008953272672
0.0714078077594
3.7
0.0006002231245
0.0828896894497
4.0
0.0005492334199
0.0921556872924
4.3
0.0007154938930
0.0983431460014
4.6
0.0008300921248
0.1008204330678
4.9
0.0005487250652
0.0992608551151
5.2
0.0002963137375
0.0936868279861
5.5
0.0015298326576
0.0844798593631
5.8
0.0026838757810
0.0723553936143
6.1
0.0032571351779
0.0583046260627
6.4
0.0030289419613
0.0435091879952
6.7
0.0022189180928
0.0292364988851
7.0
0.0013691340100
0.0167257096706
7.3
0.0010043530807
0.0070750555622
7.6
0.0012821439144
0.0011409176473
7.9
0.0018689726696
0.0005417339524
8.2
0.0021436005756
0.0021897330323
8.5
0.0016148602301
0.0091082291225
8.8
0.0002869277587
0.0196147078451
9.1
0.0012796703448
0.0327897203523
9.4
0.0022771057968
0.0474740017777
9.7
0.0021518593816
0.0623710756305
10.0
0.0009803998394
0.0761625070468Figure 12
Exact solution (black), symbolic solution (blue), and numeric solution (green).Figure 13
Error in the symbolic (black) and numeric (blue) approximations.


## 5. Conclusions

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



---

*Source: 518194-2013-08-27.xml*