ISSN 03617688, Programming and Computer Software, 2013, Vol. 39, No. 3, pp. 150–157. © Pleiades Publishing, Ltd., 2013. Original Russian Text © N.A. Malaschonok, M.A. Rybakov, 2013, published in Programmirovanie, 2013, Vol. 39, No. 3.

Symbolic–Numerical Solution of Systems of Linear Ordinary Differential Equations with Required Accuracy N. A. Malaschonok and M. A. Rybakov Institute of Mathematics, Physics and Informatics, Tambov State University, ul. Internatsional’naya 33, Tambov, 392021 Russia email: [email protected], [email protected] Received July 2, 2012

Abstract—In the paper, a symbolic–numerical algorithm for solving systems of ordinary linear differential equations with constant coefficients and compound righthand sides. The algorithm is based on the Laplace transform. A part of the algorithm determines the error of calculation that is sufficient for the required accu racy of the solution of the system. The algorithm is efficient in solving systems of differential equations of large size and is capable of choosing methods for solving the algebraic system (the image of the Laplace trans form) depending on its type; by doing so the efficiency of the solution of the original system is optimized. The algorithm is a part of the library of algorithms of the Mathpar system [15]. DOI: 10.1134/S0361768813030043

1. INTRODUCTION Modern computer algebra systems provide some methods (both numerical and symbolic–numerical) for solving systems of linear differential equations. It is no mere chance that the Laplace transform is fre quently used in recent years for solving differential equations (see, e.g., [2, 3, 13, 14]). This gives a solu tion in analytic form, which is very important in appli cations. We present a symbolic–numerical algorithm based on the Laplace transform. Below we will be concerned with systems of linear differential equations with constant coefficients and with compound righthand sides, whose components are polynomials of exponentials, which ensures the applicability of the symbolic Laplace transform. As a result, the problem amounts to solving algebraic sys tems of linear equations with polynomial coefficients. The solution of the original systems is obtained by applying the inverse Laplace transform. It is this step that involves numerical methods for calculating roots of polynomials. We propose an algorithm for deter mining the error of calculations that is sufficient for ensuring the desired calculation accuracy of the solu tion of the original system. The advantages of the proposed algorithm are as follows: • the size of the original system and the order of the derivatives do not play a crucial role; • the most efficient way of solving the algebraic system is chosen depending on the form (size, density, etc.); by doing so, we control the efficiency of solving the original system of differential equations by this method;

• the required accuracy of solution is guaranteed by the choice of the corresponding error of the numer ical portion of calculation. 2. ALGORITHM FOR SOLVING SYSTEMS OF DIFFERENTIAL EQUATIONS BY THE LAPLACE TRANSFORM METHOD Consider the system N

n

∑∑a

l (k) kj x j

l = 1, …, n,

= fl ,

j = 1k = 0

(1)

l

a kj ∈ R, of n differential equations of order N (N is the leading coefficient over all unknown functions) in n unknown (k) k functions with the initial conditions x j (0) = x 0j , k = 0, …, N – 1. The righthand side of the system con tains the functions that can be reduced to the form i

i = 1, …, I l ,

1 tl

i+1

i

f l ( t ) = f l ( t ),

tl < t < tl = 0,

I +1 tl l

, = ∞,

where i

Sl

i f l(t)

=

∑P

i

b ls t i ls ( t )e ,

s=1

i = 1, …, I l , i

and P ls (t) =

150

∑

i

li m c t . m = 0 sm M ls

l = 1, …, n,

(2)

SYMBOLIC–NUMERICAL SOLUTION

Here, xj, j = 1, …, n are the unknown functions of

Il – 1

Fl ( p ) =

(k) xj

are the derivatives of xj, k = argument t, t ≥ 0, and 0, …, N, of order k. The proposed algorithm consists of three main steps. Step 1. The Laplace transform of a given system of differential equations We denote the inverse images of the Laplace trans form of functions xj(t) and fl(t) as Xj(p) and Fl(p), respectively. Step 1.1. The Laplace transform of the lefthand side of system (1) with account for the initial condi tions is taken formally in accordance with properties of the Laplace transform. This yields n

l k a kj p X j ( p )

–

j = 1k = 0

∑∑

l k d jk ( p )x 0j ,

n

∑a

l i–k i + 1, j p

–

=

∑

i

k

ψ ls ( t – t l )e

i k

i

k

b ls t l b ls ( t – t l )

e

.

s=1 i

k

i

i

k

Here, ψ ls (t – t l ) = P ls (t) and ψ ls (t – t l ) = i

ik γ (t j = 0 lsj

k

J ls

mk

– t l ) j, and the coefficients γ lsj are calcu lated by the formula k

M ls – j ik γ lsj

∑c

=

⎛ lk s, i + r ⎜

i + r ⎞⎟ ( t k ) r . l ⎝ r ⎠

r=0

Finally, function fl(t) is represented as follows: Il – 1

fl ( t ) =

∑ [ φ ( t – t )η ( t – t ) – φ ( t – t i l

i l

i l

i l

i+1 )η ( t l

i+1

– tl

)]

i=1

I

+ φ l l ( t – t Il )η ( t – t Il ). i, i

Step 1.3. We denote by Φ l (p) the image of the i

i

+ F l ( p ),

n N–1

(3)

i

∑

∑∑

=

l = 1, …, n. Step 2.1. Let us denote

Sl

=

(5)

l k d jk ( p )x 0j

j = 1k = 0

and l is the equation number. Step 1.2. Functions fl(t) on the righthand sides of system (1) are prepared to the Laplace transform by applying the Heaviside function η(t), which equals 1 for negative values of the argument and 0 for the remaining ones. We write fl(t) as the sum of functions k tl )

l k kj p X j ( p )

j = 1k = 0

∑∑d

i=k

i φl ( t

(4)

N

∑∑a

l = 1, …, n,

N–1

i f l(t)

I,I

( p ) ] + Φ l l l ( p ).

i=1

j = 1k = 0

l

i, i + 1

– Φl

n N–1

where d jk ( p ) =

i, i l (p)

Clearing denominators in expression (4) for each l = 1, …, n, we obtain sums of exponentials with poly nomial coefficients in the numerator. Step 2. Solving a linear algebraic system with poly nomial coefficients Consider the algebraic system of n linear equations with polynomial coefficients in unknown functions Xj, j = 1, …, n:

n N–1

N

∑∑

∑ [Φ

151

i

function φ l (t – t l )η(t – t l ). Since the image of the α ( t – t* )

η(t – Laplace transform of the function (t – t*)n e – t*p n! t*) is  e , the application of the Laplace n+1 (p – α) transform of fl(t) yields PROGRAMMING AND COMPUTER SOFTWARE

l k jk ( p )x 0j

= S l ( p ).

j = 1k = 0

In this notation, system (5) is reduced to the form n

N

∑∑a

j = 1k = 0

l k kj p X j ( p )

= S l ( p ) + F l ( p ),

(6)

l = 1, …, n. For any l = 1, …, n, we reduce the expressions in the righthand sides of system (6) to a common denomi nator. Step 2.2. Recall that, in system (6), polynomials of the main matrix have real coefficients (as is seen from (1)). The method for constructing the solutions is cho sen depending on the form of the system (for example, padic method [9], modular methods, and methods based on determinant identities [7, 8] in the case where coefficients are rational, etc.). Step 3. The inverse Laplace transform Step 3.1. We denote the determinant of the main matrix of system (6) as D(p). The solution of system (6); i.e., functions Xj(p), j = 1, …, n, are represented as the 1  and the sum of fractions for which product of  D(p) the numerators are sums of exponentials with polyno mial coefficients, and the denominators, are natural powers of binomials of the form (p – α). In order to symbolically perform the inverse Laplace transform, we represent Xj(p) as sums of exponentials whose coef ficients are sums of partial fractions. Step 3.2. To prepare to perform the inverse Laplace transform, we take the partial fraction expansion of the rational coefficients of exponentials: A/(p – p*)v, p* ∈ C. To this end, it is required to find roots of the polyno mial D(p). This is precisely the step of the algorithm that involves employing numerical method. At this

Vol. 39

No. 3

2013

152

MALASCHONOK, RYBAKOV

step, we evaluate accuracy of these calculations (see Section 2). As a result, each function Xj(p) is repre sented in the form –α p A mk Xj ( p ) =  e m . (7) β mk m k ( p – p jk ) Step 3.4. The inverse Laplace transforms for func tions (7) are formally written as follows: A mk xj ( t ) =   ( β mk – 1 )!

∑∑

∑∑ m

k

β mk – 1 p ik ( t – α m )

× ( t – αm ) e η ( t – α m ), j = 1, …, n. Note that, in general, the functions x˜ l (t) obtained as the inverse Laplace transforms Xj(p) from the approximate roots of D(p) are complexvalued. It is possible either to turn to realvalued functions by transforming exponentials or take their real parts as solutions of the given system. The resulting error can be included in the required accuracy estimate. 3. ACCURACY ESTIMATE At Step 3.2 of the algorithm, we numerically find roots of the polynomial D(p). The problem is to deter mine what accuracy in the calculation of roots is suffi cient to obtain the required accuracy of the solution of a system of differential equations. We consider all functions and perform all calcula I

tions and estimates in the interval [0, T], where T > t l l for all l = 1, …, n, and is sufficiently large for the prob lem in question. Let x˜ l , l = 1, …, n, be a solution of system (1) obtained with the use of the approximate roots of D(p). We say that the solution of the initial sys tem in [0, T] has accuracy ε if max t ∈ [ 0, T ] x l ( t ) – x˜ l ( t ) < ε, l = 1, …, n. (8) It is required to find the error Δ of calculation of roots of D(p) that is sufficient to guarantee the accuracy ε. l

Let T i (p) be the (i, l)th minor of the matrix of sys tem (6). The solution Xl(p) of system (6) can be written in the form Dl ( p ) X l ( p ) =  , D(P) where

Theorem. For any ε, there exists a Δ such that (8) holds provided that |pr – p *r | < Δ. Proof. Consider the polynomial D(p + Δeiα), α ∈ [0, 2π]. Let us denote Dl ( p ) (9) X˜ l ( p ) =  . iα D ( p + Δe ) It is required to find Δ for which (8) holds. Let us estimate the inverse Laplace transform of the difference Dl ( p ) Dl ( p )   –  . D ( p ) D ( p + Δe iα ) The linearity of the Laplace transform makes it possible to independently estimate the inverse Laplace transforms n

∑

i=1

l

Ti ( p ) F l ( p )  – D(p)

∑ [ F ( p ) + S ( p ) ]T ( p ). l

l

n

∑

i=1

l

Ti ( p ) – S l ( p )  D(p)

We denote by pr the exact roots, and by p *r , r = 1, …, n, the approximate roots of the polynomial D(p). The accuracy of calculation is smaller than Δ provided that |pr – p *r | < Δ.

∑

(10)

n

l

Ti ( p ) S l ( p )  . iα D ( p + Δe ) i=1

∑

(11)

Let us write (10) as l ⎛ T li ( p ) Ti ( p ) ⎞ F l ( p ) ⎜  ⎟ .  –  D ( p ) D ( p + Δe iα )⎠ ⎝ i=1 n

∑

(12)

Both the function Fl(p) and the bracketed differ ence have inverse Laplace transforms. The inverse Laplace transform Ψl(t) of expression (10) can be obtained as the sum of convolutions of fl(t) and the i

inverse Laplace transform Ω l (t) of the function l

l

Ti ( p ) Ti ( p ) .  –  D ( p ) D ( p + Δe iα )

(13)

Let us represent (13) in the form l

l

l

iα

iα

l

T i ( p ) T i ( p + Δe ) T i ( p + Δe ) – T i ( p )  +  .  –  iα D ( p ) D ( p + Δe iα ) D ( p + Δe ) i

(14)

i

Let Q l (t) and Q l ( t ) be the inverse Laplace transforms l

l

l

iα

Ti ( p ) Ti ( p ) T i ( p + Δe ) ,  and   –  of the functions  iα D(p) D(p) D ( p + Δe ) respectively. From the properties of the Laplace trans form, it follows that ΔT

i

i

Q l ( t ) = ( 1 – e )Q l ( t ).

l i

i=1

l

Ti ( p ) F l ( p )   iα D ( p + Δe ) i=1

and

m

Dl ( p ) =

n

(15)

i

i

Now, let us estimate Q l (t). The function Q l (t) can be written in the form μ

μ –μ ⎛ r i ⎞ pt t r = Ꮾ rμ ⎟ e r , ⎜ ( μ r – μ )! )⎠ ⎝ r=1 μ=1 R

i Ql ( t )

∑ ∑

PROGRAMMING AND COMPUTER SOFTWARE

Vol. 39

No. 3

(16)

2013

SYMBOLIC–NUMERICAL SOLUTION

where pr are roots of D(p) of multiplicity μr and R is the number of different roots. The coefficients Ꮾ rμ can be calculated and esti mated as Taylor coefficients of the function (p – l r μr Ti ( p ) p )   in the δneighborhood of the point pr. Let D(p) δ < 1/2min|pr – ps |, r, s = 1, …, R, δ < 1. We should have Δ ≤ δ; otherwise, we take Δ = δ. From the Cauchy inequalities for the Taylor coeffi cients, it follows that

i mn M e ρT . Q l ( t ) ≤  il mn c mn δ Consider the last term in (14). By Lagrange’s formula,

i

Ꮾ rμ ≤ ⎛⎝ max p – pr i

l

μ r T i ( p )⎞ μ = δ ( p – p r ) ⎠ /δ . D(p)

Consider the polynomial D(p) = l μr Ti ( p )

∑

ν mn c p . ν=0 ν

Since

l

iα 1 = Δe   2πi

i

l

i Δ  Qi ( t ) . Q l ( t ) ≤  l δ–Δ

l

⎧ mn ⎫ cν ⎪ ⎪ ⎪ ⎪ =0  , ρ = max ⎨ 1; ν ⎬ c mn ⎪ ⎪ ⎪ ⎪ ⎩ ⎭ which may be taken so as to have

∑

l

Ti ( p ) .

The element of the matrix whose determinant is l denoted by T i (p) is the polynomial γ η, λ (p), where λ is the row number, λ ≠ l, and η is the column number, η ≠ i. Let us denote mη, λ = max p = ρ γ η, λ ( p ) . From the Hadamard inequalities, we find that n

l max p = ρ T i ( p )

≤

∏

∏

Ωl ( t ) ≤ ⎛ e ⎝ i

ΔT

ρT Δ mn – 1 + ⎞  M il e . δ – Δ⎠ c δ mn mn

The inverse Laplace transform fi(t) of the function Fi(p) is a compound function. Each component fi(t) is a sum of exponentials with polynomial coefficients. We estimate it by finding maxima of the coefficients over the corresponding intervals and values of the exponentials on their righthand ends. Let us denote max t ∈ [ 0, T ] f i ( t ) = Mf i . (23) Estimating from above the convolution of the inverse Laplace transforms of functions (12), we obtain the fol lowing estimate for Ψl(t): max t ∈ [ 0, T ] Ψ l ( t ) ≤ ( e

ΔT

n

mnT  Mf M e ρT . – 2 )  i il mn c mn δ i = 1

∑

Let us denote

Laplace transform. Denoting

n

∑

(22)

∑

η = 1, η ≠ i λ = 1, λ ≠ l

M il =

i

Using (15), (18), (20), and (21), we estimate Ω l (t):

mnT  Mf M e ρT Λ l =  (24) i il mn c mn δ i = 1 and consider (11). The polynomial Si(p) has no inverse

2 m η, λ .

Let us denote n

(21)

n

n

∑

l

iα

T i ( p + Δe ) – T i ( p )  : iα D ( p + Δe )

we need to estimate T i ( p ) in the δneighborhood of pr . We denote by ρ the radius of the disk that contains all zeros of D(p). For such ρ, one may consider the number

=ρ

(19)

the inverse Laplace transform Q l ( t ) of the function

∏

l

∫

ξ–p

l

Ti ( ξ )    dξ. iα ξ – p Δe ( – ) ( ξ – p ) =δ

l l iα Δ T i ( p ) – T i ( p + Δe ) ≤  M il . (20) δ–Δ Like in the previous case, we obtain the inequality for

l

T i ( p ) ≤ max p

iα

Using (19), we have:

s = 1, s ≠ r

=δ

l

(18)

T i ( p ) – T i ( p + Δe )

Ti ( p ) ( p – p r )  =  , R D(p) μ c mn p – pr s

max p – pr

153

2 m η, λ .

(17)

η = 1, η ≠ i λ = 1, λ ≠ l

Taking into account that |p – ps | ≥ δ and μ ≤ |mur, we obtain: M il i Ꮾ rμ ≤  . mn c mn δ i

As a result, we obtain estimate Q l (t) in the interval [0, T]: PROGRAMMING AND COMPUTER SOFTWARE

Vol. 39

∑

l n S (p) T i (p) i=1 i

= Kl(p),

Kl ( p )  = Nl(p), we write (11) as follows: D(p) iα

K l ( p ) K l ( p + Δe ) Kl ( p ) Kl ( p )  =   –   –  iα D ( p ) D ( p + Δe ) D ( p ) D ( p + Δe iα ) 1 + Δe  2πi iα

No. 3

2013

∫

ξ–p

N (p) l   dξ. iα ξ – p Δe ( – ) ( ξ – p ) =δ

(25)

154

MALASCHONOK, RYBAKOV

The inverse Laplace transforms Rl(t) of the func tions Nl(p) can be represented in the form (16). Similar estimates result in the inequalities: max t ∈ [ 0, T ] R l ( t ) n

mnT ≤   max p mn c mn δ i = 1 From (5), it is clear that

∑

≤ρ

ρT

S i ( p ) M il e .

m–1

max p

≤ ρ Si ( p ) ≤

∑a

i ν + 1, j

ρ

ν–k

k

x oj .

(26)

(27)

ν=k

We set n

ρT mnT (28) L l =   max p ≤ ρ S i ( p ) M il e . mn c mn δ i = 1 Proceeding in the same way, we obtain estimates for Θl(t), which is the inverse Laplace transform of func tion (11): ΔT Δ ⎞ L . Θ l ( t ) ≤ ⎛ e – 1 +  (29) l ⎝ δ – Δ⎠ Finally, using (22) and (29), we obtain the estimates for xl(t): max t ∈ [ 0, T ] x l ( t ) – x˜ l ( t ) (30) ΔT Δ ≤ ⎛ e – 1 + ⎞ ( Λ l + L l ). ⎝ δ – Δ⎠

∑

k Let us assume that  δ < Δ < δ, k ≥ 0. Then, ine k+1 quality (30) takes the form max t ∈ [ 0, T ] x l ( t ) – x˜ l ( t ) (31) ΔT ≤ ( e – 1 + k ) ( Λ l + L l ). Based on (9) and (31), we require that ΔT

(32) ( e – 1 + k ) ( Λ l + L l ) < ε. Using (32), we obtain an estimate for Δ for any l (we denote it by Δl): –1 1 (33) Δ l <  ln ( ε ( Λ l + L l ) + 1 – k ). T Let us take –1

(34) k < min l ( Λ l + L l ) . Finally, for the desired Δ, we take Δ = minlΔl. The proof of the theorem is complete. Thus, the algorithm for finding an error that is suf ficient to ensure the desired accuracy of the solution is as follows. Step 4. Algorithm for error calculation Step 4.1. Given the matrix of system (6), calculate the numbers m η, λ = max p = ρ γ η, λ ( p ) and Mil (see (17)). Step 4.2. Calculate Mfi (see (23)). Step 4.3. Calculate Λl (see (24)) from the values of δ and ρ.

Step 4.4. Calculate Ll (see (28)), estimate the poly nomials Sl(p) as in (27). Step 4.5. The final result for Δl and Δ = minlΔl is by formulas (33) and (34). Thus, in the course of the proof of the theorem, we have obtained formulas for evaluating the error Δ by the required accuracy ε of the solution of system (1). These formulas are given in Step 4. 4. ON THE COMPLEXITY OF THE ALGORITHM The complexity of the algorithm is governed by the complexities of the three main operations: solution of the linear algebraic system with polynomial coeffi cients at Step (2.2), solution of the linear algebraic sys tem with constant coefficients at Step (3.3), and cal culation of roots of the polynomial at Step 3.2. At Step (2.2), the complexity essentially depends on the calculation of the inverse matrix and multipli cation of the polynomial matrix by the polynomial vector. The complexity of various algorithms for the above calculations, as well as that of solving system at Step 3.3, is calculated accordingly (see, e.g., [6]). Here, one of the most efficient methods is the padic method ([9]), while modular method is useful in par allel computations. The complexity of Step 3.2 is cal culated, for example, in [1]. 5. EXAMPLE Consider the following system of differential equa tions: ⎧ x ''' 1 – x '1 – 2x 1 – x ''' 2 + x2 = f1 ⎨ ⎩ 3x ''' 1 + x '' 1 – 2x '1 + x ''' 2 + x2 = f2 . i

The functions f1 and f2 and the numbers t l are as follows: 1

t

f1 = e , 1

2

2 2t

t 1 = 0,

1

t 1 = 1;

2t

t 2 = 0,

1

t 2 = 1.

f1 = t e , 2

t

f 2 = te , f 2 = e , The initial conditions: 0

x 01 = 5,

1

x 01 = 30,

1

x 02 = 20.

2

2

2

x 01 = 10,

0

x 02 = 4,

2

x 02 = 14, Step 1.1. 3

3

– 2X 1 – pX 2 + p X 1 + X 2 – p X 2 2

2

– ( 10 – 4p – 4p + 5 ( – 1 + p ) ); 2

3

3

– 2pX 1 + p X 1 + 3p X 1 + X 2 + p X 2 2

– ( 110 + 14p + 4p + 10 ( 1 + 3p ) 2

+ 5 ( – 2 + p + 3p ) ). Step 1.2.

PROGRAMMING AND COMPUTER SOFTWARE

Vol. 39

No. 3

2013

SYMBOLIC–NUMERICAL SOLUTION 2

1

2

1

2

1

2

1

f 1 := ( f 1 – f 1 )UnitStep [ ( t – t 1 ) ] + f 1 UnitStep [ t ]; f 2 := ( f 2 – f 2 )UnitStep [ ( t – t 2 ) ] + f 2 UnitStep [ t ]. Step 1.3. 1–p

2–p

The final result for an arbitrary T is given below: Step 4.3. Λ1 = (6240T3 + 20160T2)e5T, Λ2 = (6240T3 + 6720T2)e5T. Step 4.4. 4T

2

1 e e ( p – 2p + 2) ; F 1 =  –  +  –1+p –1+p (– 2 + p) 2–p

L 1 = 9438240e ;

e e p . 1  –  F 2 =  +  2 – 2 + p ( – 1 + p ) ( – 1 + p )2 Step 2.1. 3

– 2X 1 – pX 1 + p X 1 + X 2 – p X 2 1–p

2 2 1  e = 10 – 4p – 4p + 5 ( – 1 + p ) +  –  –1+p –1+p 2–p

p *2 = –0.594937842169665 – i0.830713582043548, p *3 = –0.594937842169665 + i0.830713582043548,

2

( p – 2p + 2) ; + e 3 (– 2 + p) 2

p *4 = 0.355937297682321 – i0.513128324882554,

3

p 5* = 0.355937297682321 + i0.513128324882554,

3

– 2pX 1 + p X 1 + 3p X 1 + X 2 + p X 2 2

2

= 110 + 14p + 4p + 10 ( 1 + 3p ) + 5 ( – 2 + p + 3p ) 2–p

1–p

e p 1 e +  + 2 – 2 . – 2 + p (– 1 + p) (– 1 + p) Step 2.2. –p

2

2

X 1 ( p ) = e ( 8e + 2e – 4ep – 4e p + 2ep 2 2

3

2 3

+ 2e p + 9ep – 6e p – 19ep 2 4

5

2

4

2 5

+ 12e p + 11ep – 8e p 6

2 6

3

– 2ep + 2e p )/ ( ( – 2 + p ) ( – 1 + p ) ( – 2 + p 2

3

4

5

6

– p – 4p – 3p + p + 4p ) ) 2

3

+ ( – 856 + 1692p – 982p + 1061p – 1991p 5

6

+ 1398p – 412p + 160p 8

9

4

3

3

4

5

6

– p – 4p – 3p + p + 4p ) ); –p

2

2

X 2 ( p ) = e ( – 8e – 32ep + 24e p + 64ep 2 2

3

2 3

4

2

2 4

– 28e p – 32ep + 17e p – 24ep – 5e p 5

2 5

6

2 6

+ 34ep – 3e p – 14ep + 5e p + 2ep 2 7

3

7

2

– 2e p )/ ( ( – 2 + p ) ( – 1 + p ) ( – 2 + p 2

3

4

5

6

– p – 4p – 3p + p + 4p ) 2

3

+ ( 1776 – 4576p + 3568p – 1404p + 2465p 5

6

– 2751p + 841p + 133p 8

9

10

3

3

4

4

7

+ 2p – 68p + 16p )/ ( ( – 2 + p ) ( – 1 + p ) 2

5

2

6

( – 2 + p – p – 4p – 3p + p + 4p ) ). Step 3.1. 2

3

4

p *6 = 1.22800108897469. Finally, the solution of this system of differential equations is as follows: x1(t) = 10.031249e–t – 1.25et + 5.538602e1.228001t + 2e0.355937t(–3.735568cos(0.513128t) + 15.529795sin(0.513128t)) + 2e–0.594937t(–0.924357cos(0.830713t) + 0.061193sin(0.830713t)); x2(t) = 10.03125e–t + 0.5et – 8.948223e1.228001t + 0.5ett + 2e0.355937t(–0.493116cos(0.513128t) + 33.959275sin(0.513128t)) + 2e–0.594938t(1.701602cos(0.830713t) + 0.929609sin(0.830713t)).

7

– 95p + 20p )/ ( ( – 2 + p ) ( – 1 + p ) ( – 2 + p 2

4T

L 2 = 3874080e .

Let us take, for example, k = 1 minl(Λl + Ll)–1 and 2 T = 2. For ε = 0.01, we have Δ = 8.06 × 10–14, and for ε = 0.001, Δ = 7.99 × 10–15. We take Δ = 10–15 and find approximate roots of D(p). p *1 = –1.00000000000000,

1–p

3

155

5

6

D ( p ) = – 2 + p – p – 4p – 3p + p + 4p . PROGRAMMING AND COMPUTER SOFTWARE

6. SOFTWARE PACKAGES FOR SOLVING SYSTEMS OF DIFFERENTIAL EQUATIONS The algorithm presented here was implemented in the Mathpar computer algebra system. For the system itself and the user guide, we refer the reader to the website http://mathpar.com. The developed software package LaplaceTransform consists of three main classes. 1. Class SystemLDE includes procedures aimed at solving systems of differential equations. Procedure systDifEquationToMatrix calcu lates the polynomial matrix A(p) by the Laplace trans form of the lefthand side of system (1). Procedure initCondLaplaceTransform calcu lates the vector function ᏿(p) with elements Sl(p), l = 1, …, n, which is obtained as the result of the Laplace transform of the initial conditions. primeFractionForSingle is Procedure designed for expanding the fraction 1/D(p) with a fac tored denominator D(p) into a sum of partial fractions. Vol. 39

No. 3

2013

156

MALASCHONOK, RYBAKOV

The numerators of these fractions are found by the method of indeterminate coefficients. Procedure pimeFraction is designed for expand V ing sums of fractions 1/w v ( p ) with factored v=1 denominators wv(p) into sums of fractions (for Step 3.1). 2. Class LaplaceTransform involves procedures for calculation of the direct Laplace transform of func tions. Procedure vectorLaplaceTransform calcu lates the Laplace transform (Step 1.2, 1.3) of the right hand side of system (1). 3. Class InverseLaplaceTransform contains proce dures for calculation of the inverse Laplace transform. Procedure vectorLaplaceTransform calculates the inverse Laplace transform of vector functions X(p) with elements Xj(p). As a result, we have a solution x(t) of system (1). The software package laplaceTransform uses the following software packages of the Mathpar system: number, matrix, polynom, func. The software package number is aimed at specifying numerical sets and operations on them. Besides, the package number contains two additional classes: Prod uct and SumOfProducts. 1. Class Product contains procedures designed for operations on objects of the form:

∑

∏α

ni i ;

i

here n, i ∈ ⺞, and αi are scalar objects of the Mathpar system (numbers, polynomials, functions). All ele ments in the product are ordered. 2. Class SumOfProducts involves procedures for operations on objects of the form

∑∏ j

n α ijij .

7. CONCLUSIONS In conclusion, we restate advantages of the algo rithm proposed in this paper: 1. The Laplace transform paves the way for symbol ical solution of systems of differential equations, since it leads to solution of an algebraic system of linear equations. 2. Depending on properties of the algebraic system, the most efficient method for solving this system is chosen. 3. The size of the system and the order of the deriv atives may be arbitrary. 4. The numerical stage of the solution involves esti mation of the error sufficient for obtaining the required accuracy of the solution of the system of dif ferential equations. 5. The algorithm is implemented in the Mathpar computer algebra system. 8. REMARKS Solution of the above example in Mathpar on a 1 Gb PC took about one second. The same example was implemented in Mapple on a PC with 4 Gb operating memory: ⎧d d d d [> dsys := ⎨  ⎛  ⎛  x ( t )⎞ ⎞ – ⎛  x ( t )⎞ – 2x ( t ) ⎝ dt ⎝ dt ⎠ ⎠ ⎝ dt ⎠ dt ⎩ d d d –  ⎛  ⎛  y ( t )⎞ ⎞ + y ( t ) ⎠⎠ dt ⎝ dt ⎝ dt 2 2t

t

d d d 3  ⎛  ⎛  y ( t )⎞ ⎞ + y ( t ) ⎠⎠ dt ⎝ dt ⎝ dt

i

Procedure primeFractions of class SumOfProd ucts performs expansion of objects into sums of partial fractions. Procedure reduction collects similar terms in the object SumOfProducts. Package martix is aimed at operations on numerical and functional matrices. Procedure matrixAdjointDet of package matrix enables one to calculate the determinant and adjoint matrix for the matrix of polynomials at Step 2. Procedure solveLAE calculates the solution of system at Step 2. Package polynom is designed for operations on polynomials. One of the procedures of the package, namely procedure factorOfPol, factors a polyno mial of one variable in the field of complex numbers. Package func is designed for operations on super positions of transcendental and rational functions (see [15]). The reader is referred to [15] for a description of these software packages.

t

= ( t e – e )Heaviside ( t – 1 ) + e Heaviside ( t ),

2t

t

t

= ( e – te )Heaviside ( t – 1 ) + te Heaviside ( t ), D ( D ( x ) ) ( 0 ) = 30, D ( x ) ( 0 ) = 10,

D ( D ( y ) ) ( 0 ) = 20,

D ( y ) ( 0 ) = 14,

x ( 0 ) = 5,

y(0) = 4 } [> dsolve(dsys). The memory overfull message was flagged after sev eral minutes of computation. REFERENCES 1. Akritas, A.G., Elements of Computer Algebra with Appli cations, New York: Wiley, 1989. 2. Burghelea, D. and Haller, S., Laplace transform, dynam ics and spectral geometry. arXiv:math.DG/0405037v2. Cited January 17, 2005. 3. Dahiya, R.S. and SaberiNadjafi, J., Theorems on n dimensional Laplace transforms and their applications, Electr. J. Differential Equations, Conf. 02, 1999, pp. 61–74.

PROGRAMMING AND COMPUTER SOFTWARE

Vol. 39

No. 3

2013

SYMBOLIC–NUMERICAL SOLUTION 4. Davenport, J.H., Siret, Y., and Tournier, E., Computer Algebra: Systems and Algorithms for Algebraic Computa tion, London: Academic, 1988. 5. Gathen von zur, J. and Gerhard, J., Modern Computer Algebra, Cambridge: Cambridge Univ. Press, 1999. 6. Malaschonok, G.I., Algorithm for the solution of sys tems of linear equations in commutative rings, Effective Methods in Algebraic Geometry, Mora, T. and Traverso, C., Eds., Progress in Mathematics, vol. 94, Birkhauser, 1991, pp. 289–298. 7. Malaschonok, G.I., Recursive method for the solution of systems of linear equations, in: Computational Math ematics (15th IMACS World Congress, vol. I, Berlin, August 1997), Berlin: Wissenschaft Technik, 1997, pp. 475–480. 8. Malaschonok, G.I., Effective matrix methods in com mutative domains, in Formal Power Series and Algebraic Combinatorics, Berlin: Springer, 2000, pp. 506–517. 9. Malaschonok, G.I., Solution of systems of linear equa tions by the padic method, Program. Comp. Software, 2003, vol. 29, no. 2, pp. 59–71. 10. Malaschonok, N., An algorithm to settle the necessary exactness in Laplace transform method, in Computer

PROGRAMMING AND COMPUTER SOFTWARE

11.

12.

13.

14.

15.

Vol. 39

157

Science and Information Technologies, Yerevan, 2005, pp. 453–456. Malaschonok, N., Estimations in Laplace transform method for solutions of differential equations in sym bolic computations, in Differential Equations and Com puter Algebra Systems, Minsk, 2005, pp. 195–199. Malaschonok, N., Parallel Laplace method with assured accuracy for solutions of differential equations by symbolic computations, in Computer Algebra and Scientific Computing, CASC 2006, Berlin: Springer, 2006, pp. 251–261. Mizutani, N. and Yamada, H., Modified Laplace trans formation method and its applications to anharmonic oscillator, Int. J. Mod. Phys. A, 1997, vol. 12, no. 31, pp. 5687–5709. Podlubny, I., The Laplace transform method for linear differential equations of the fractional order. arXiv:funct an/9710005v1. Cited October 30, 1997. Malaschonok, G.I., Project of parallel computer alge bra, Tambov Univ. Reports. Series: Natural Techn. Sci., 2010, vol.15, no. 6, pp. 1724–1729.

Translated by A. Alimov

No. 3

2013