Numerical solution of differential–algebraic equation systems and applications

Applied Mathematics and Computation 154 (2004) 405–413 www.elsevier.com/locate/amc

Numerical solution of differential–algebraic equation systems and applications Ercan C ß el_ık *, Mustafa Bayram € niversitesi, Fen-Edebiyat Fak€ultesi, Matematik B€ol€um€u, 25240 Erzurum, Turkey Atat€ urk U

Abstract In this paper, the method is developed to differential–algebraic equations systems. More effective method is presented and illustrated by numerical example. The solution for a differential–algebraic equation can be expanded up to arbitrary order using MAPLE computer algebra systems. First we calculate power series of the given equations system then transform it into Pade series form, which give an arbitrary order for solving differential–algebraic equation numerically. Ó 2003 Elsevier Inc. All rights reserved. Keywords: Differential–algebraic equation; Arbitrary order; Power series; Pade series

1. Introduction Differential–algebraic equations are normally obtained when modeling chemical engineering systems. Chemical processes are modeled dynamically using differential–algebraic equations. Chemical processes are inherently nonlinear and multivariable and are typically modeled by coupled differential and algebraic equation. A system of DAEs is characterized by its index, which is the number of differentiations required to convert it into a system of ODEs. DAEs with index > 1 are generally hard to solve and are still under active research. DAEs with index-1 are usually encountered. Several popular numerical solvers

*

Corresponding author. E-mail address: [email protected] (E. C ß el_ık).

0096-3003/$ - see front matter Ó 2003 Elsevier Inc. All rights reserved. doi:10.1016/S0096-3003(03)00719-7

406

E. C ß el_ık, M. Bayram / Appl. Math. Comput. 154 (2004) 405–413

for index-1 DAEs are available at present, e.g. DASSL [1], RADAU5 [2] and _ LIMEX [3]. A system of initial value differential–algebraic equations (DAEs) can be written in its general nonlinear implicit or fully implicit form as in Eq. (1.1a), subject to initial conditions in Eq. (1.1b), where F 2 Rn , y 2 Rn and x 2 R: F ðx; yðxÞ; y 0 ðxÞÞ ¼ 0; yðx0 Þ ¼ y0 ;

y 0 ðx0 Þ ¼ y00

ð1:1aÞ ð1:1bÞ

the existence of algebraic constraints on the variables are expressed by the singularity of the Jacobian matrix oF =oy 0 . The initial conditions for y and y 0 defined by Eq. (1.1b) must be consistent in order for Eq. (1.1a) to admit solution [4].

2. The method A differential–algebraic equation has the form F ðy; y 0 ; xÞ ¼ 0;

ð2:1Þ

with initial values yðx0 Þ ¼ y0 ;

y 0 ðx0 Þ ¼ y1 ;

where F and y are both vector functions for which we assumed sufficient differentiability and the initial values to be consistent, i.e. F ðy0 ; y00 ; x0 Þ ¼ 0:

ð2:2Þ

The solutions of (2.1) can be assumed that y ¼ y0 þ y1 x þ ex2 ;

ð2:3Þ

where e is a vector function which is the same size as y0 and y00 . Substitute (2.3) into (2.1) and neglect higher order term, we have the linear equation of e in the form Ae ¼ B;

ð2:4Þ

where A and B are constant matrices. Solving Eq. (2.4); the coefficients of x2 in (2.3) can be determined. Repeating above procedure for higher order terms, we can get the arbitrary order power series of the solutions for (2.1). The Power series given by above procedure can be transformed into Pade series and we have numerical solution of differential–algebraic equation in (2.1) [5,6].

E. C ß el_ık, M. Bayram / Appl. Math. Comput. 154 (2004) 405–413

407

3. Power series of solution for differential–algebraic equations We define another type power series in the form f ðxÞ ¼ f0 þ f1 x þ f2 x2 þ    þ ðfn þ p1 e1 þ    þ pm em Þxn ;

ð3:1Þ

where p1 ; p2 ; . . . ; pm are constants. e1 ; e2 ; . . . ; em are basis of vector e, m is size of vector e. y is a vector in (2.3) with m elements. Every element can be represented by the power series in (3.1). Therefore we can write yi ¼ yi;0 þ yi;1 x þ yi;2 x2 þ    þ ei xn ;

ð3:2Þ

where yi is ith element of y. Substitute (3.2) into (2.1), we can get fi ¼ ðfi;n þ pi;1 e1 þ    þ pi;m em Þxnj þ Oðxnjþ1 Þ;

ð3:3Þ

where fi is ith element of f ðy; y 0 ; xÞ in (2.1) and if f ðy; y 0 ; xÞ has y 0 then j is 1, otherwise 0. From (3.3) and (3.4), we can determine the linear equation in (2.4) as follow Ai;j ¼ Pi;j ;

ð3:4Þ

Bi ¼ fi;n :

ð3:5Þ

Solve this linear equation, we have ei ði ¼ 1; . . . ; mÞ. Substitute ei into (3.2); we have yi ði ¼ 1; . . . ; mÞ which are polynomials of degree n. Repeating this procedure from (3.2) to (3.5), we can get the arbitrary order power series solution for differential–algebraic equations in (2.1). Let step size of x to be h and substitute it into the power series of y and derivative of y, we have y and y 0 at x ¼ x0 þ h. If we repeat above procedure, we can have numerical solution of differential–algebraic equations in (2.1). But we do not want to obtain numerical solution of differential–algebraic equation from power series. We would like to put (2.1) into Pade series form then obtain numerical solution of it [5,6].

4. Pade´ approximation (series) Suppose that we are given a power series f ðxÞ, so that f ðxÞ ¼

1 X i¼0

ai x i :

P1

i¼0

ai xi , representing a function

ð4:1Þ

408

E. C ß el_ık, M. Bayram / Appl. Math. Comput. 154 (2004) 405–413

A Pade approximant is a rational fraction ½L=M ¼

p0 þ p1 x þ    þ pL xL ; q0 þ q1 x þ    þ qM x M

ð4:2Þ

which has a Maclaurin expansion which agrees with (4.1) as far as possible. Notice that in (4.2) there are L þ 1 numerator coefficients and M þ 1 denominator coefficients. There is a more or less irrelevant common factor between them, and for definiteness we take q0 ¼ 1. This choice turns out to be an essential part of the precise definition and (4.2) is our conventional notation with this choice for q0 . So there are L þ 1 independent numerator coefficients and M independent denominator coefficients, making L þ M þ 1 unknown coefficients in all. This number suggests that normally the ½L=M ought to fit the power series (4.1) through the orders 1; x; x2 ; . . . ; xLþM in the notation of formal power series, 1 X

ai x i ¼

i¼0

p0 þ p1 x þ    þ pL xL þ OðxLþMþ1 Þ: q0 þ q1 x þ    þ qM x M

ð4:3Þ

Multiply the both side of (4.3) by the denominator of right side in (4.3) and compare the coefficients of both sides in (4.3). We have al þ

m X

alk qk ¼ pl

ðl ¼ 0; . . . ; MÞ;

ð4:4Þ

alk qk ¼ 0

ðl ¼ M þ 1; . . . ; M þ LÞ:

ð4:5Þ

k¼1

al þ

L X k¼1

Solve the linear equation in (4.5), we have qk ðk ¼ 1; . . . ; LÞ. And substitute qk into (4.4), we have pl ðl ¼ 0; . . . ; MÞ. Therefore, ½L=M P1 wei have constructed aLþM Pade approximant, which agrees with a x through order x . If i¼0 i M 6 L 6 M þ 2, where M and L are the degree of numerator and denominator in Pade series, respectively, then Pade series gives an A-stable formula for an ordinary differential equation [5,6].

5. The numerical example We consider the following differential–algebraic equation y10 þ y3 y20  ðy2 þ 1Þy30 ¼ y1 þ 1 þ sinðxÞ; ðy3 þ 1Þy10 þ y1 y20 ¼ ex ; y1 y2 y3  e

x

sinðxÞ cosðxÞ ¼ 0

ð5:1Þ

E. C ß el_ık, M. Bayram / Appl. Math. Comput. 154 (2004) 405–413

and initial values 0 1 1 yð0Þ ¼ @ 0 A; 1

1 1 y 0 ð0Þ ¼ @ 1 A: 0

409

0

ð5:2Þ

The exact solution is y1 ðxÞ ¼ ex ; y2 ðxÞ ¼ sinðxÞ;

ð5:3Þ

y3 ðxÞ ¼ cosðxÞ: From initial values, the solutions of (5.1) can be supposed as y1 ðxÞ ¼ 1  x þ e1 x2 ; y2 ðxÞ ¼ x þ e2 x2 ;

ð5:4Þ

2

y3 ðxÞ ¼ 1 þ e3 x : Substitute (5.4) into (5.1) and neglect higher order terms, we have e2 x2 þ Qðx3 Þ ¼ 0; ð2 þ 4e1 þ 2e2 Þx þ Qðx2 Þ ¼ 0;

ð5:5Þ

2

ð2 þ 2e1 þ 2e2  2e3 Þx þ Qðx Þ ¼ 0: These formula corresponds to (5.3). The linear equation that corresponds to (5.4) can be given in the following Ae ¼ B;

ð5:6Þ

where 0

0 A ¼ @4 2

1 2 2

1 0 0 A; 2

0 1 0 B ¼ @2A 2 from Eq. (5.6), we have 0 1 0:5 e¼@ 0 A 0:5

ð5:7Þ

ð5:8Þ

ð5:9Þ

410

E. C ß el_ık, M. Bayram / Appl. Math. Comput. 154 (2004) 405–413

and y1 ðxÞ ¼ 1  x þ 0:5x2 ; y2 ðxÞ ¼ x;

ð5:10Þ

y3 ðxÞ ¼ 1  0:5x2 : From (5.10) the solutions of (5.1) can be supposed as y1 ðxÞ ¼ 1  x þ 0:5x2 þ e1 x3 ; y2 ðxÞ ¼ x þ e2 x3 ; 2

ð5:11Þ

y3 ðxÞ ¼ 1  0:5x þ e3 x

3

in like manner. Substitute (5.1) into (5.11) and neglect higher order terms, we have ð0:166667 þ e2 Þx3 þ Qðx4 Þ ¼ 0; ð1:5 þ 6e1 þ 3e2 Þx2 þ Qðx3 Þ ¼ 0; 2

ð5:12Þ

3

ð1 þ 3e1 þ 3e2  3e3 Þx þ Qðx Þ ¼ 0: From (5.12) we have linear equation 0 10 1 0 1 0 1 0 e1 0:166667 @ 6 3 0 A@ e2 A ¼ @ 1:5 A 3 3 3 1 e3 solve this linear equation, we have 0 1 0:166667 e ¼ @ 0:166667 A 0

ð5:13Þ

ð5:14Þ

therefore y1 ðxÞ ¼ 1  x þ 0:5x2  0:166667x3 ; y2 ðxÞ ¼ x  0:166667x3 ;

ð5:15Þ

2

y3 ðxÞ ¼ 1  0:5x : Repeating above procedure we have y1 ðxÞ ¼ 1  x þ 0:5x2  0:166666667x3 þ 0:04166666667x4  0:008333333333x5 þ 0:001388888889x6  0:0001984126984x7 þ 0:00002480158730x8  0:2755731922x9 ; y2 ðxÞ ¼ x  0:1666666667x3 þ 0:008333333333x5  0:0001984126984x7 þ 0:2755731922105 x9 ; y3 ðxÞ ¼ 1  0:5x2 þ 0:04166666667x4  0:001388888889x6 þ 0:00002480158730x8 :

E. C ß el_ık, M. Bayram / Appl. Math. Comput. 154 (2004) 405–413 120

y1(x)

100

Pade approx.

y1(x)

80

60

40

20

0 -6

-4

-2

0

2

4

6

x

Fig. 1. Values of y1 ðxÞ, its [5/4] Pade approximant.

y2(x)

1.0

Pade approx.

y2(x)

0.5

0.0

-0.5

-1.0 -4

-2

0

2

4

x

Fig. 2. Values of y2 ðxÞ, its [5/4] Pade approximant.

y3(x)

1.0

Pade approx.

y3(x)

0.5

0.0

-0.5

-1.0 -4

-2

0

x

2

4

Fig. 3. Values of y3 ðxÞ, its [5/4] Pade approximant.

411

412

E. C ß el_ık, M. Bayram / Appl. Math. Comput. 154 (2004) 405–413

Power series y1 ðxÞ can be transformed into Pade series (Fig. 1). P ¼ ½5=4 ¼ 1  0:5555555556x þ 0:1388888889x2  0:01984126984x3 þ 0:001653439153x4  0:00006613756614x5 =1  0:4444444444x þ 0:08333333333x2 þ 0:007936507937x3 þ 0:0003306878307x4 : Table 1 Numerical solution of y1 ðxÞ (Panel A), y2 ðxÞ (Panel B) and y3 ðxÞ (Panel C) in (5.1) x

ex

P (Pade)

jex  P (Pade)j

Panel A 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

0.9048374180 0.8187307531 0.7408182207 0.6703200460 0.6065306597 0.5488116311 0.4965853038 0.4493289641 0.4065696597 0.3678794412

0.9048374182 0.8187307533 0.7408182204 0.6703200456 0.6065306602 0.5488116358 0.4965853040 0.4493289638 0.4065696595 0.3678794404

0.2109 0.2109 0.3109 0.4109 0.5109 0.3109 0.2109 0.3109 0.2109 0.8109

sinðxÞ

q (Pade)

j sinðxÞ  q (Pade)j

0.09983341665 0.1986693308 0.2955202067 0.3894183423 0.4794255386 0.5646424734 0.6442176872 0.7173560909 0.7833269096 0.8414709848

0.09983341664 0.1986693308 0.2955202067 0.3894183424 0.4794255383 0.5646424732 0.6442176876 0.7173560931 0.7833269172 0.8414710074

0.1109 0 0 0.1109 0.3109 0.2109 0.4109 0.22108 0.76108 0.226107

cosðxÞ

r (Pade)

j cosðxÞ  r (Pade)j

0.9950041653 0.9800665778 0.9553364891 0.9210609940 0.8775825619 0.8253356149 0.7648421873 0.6957067093 0.6216099683 0.5403023059

0.09983341664 0.1986693308 0.2955202067 0.3894183424 0.4794255383 0.5646424732 0.6442176876 0.7173560931 0.7833269172 0.8414710074

0.1109 0.5109 0.7109 0.3109 0.5109 0.17108 0.106107 0.399107 0.1273106 0.3599106

Panel B 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

Panel C 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

E. C ß el_ık, M. Bayram / Appl. Math. Comput. 154 (2004) 405–413

413

Power series y2 ðxÞ can be transformed into Pade series (Fig. 2). q ¼ ½5=4 ¼ x  0:1338383838x3 þ 0:003312890813x5 =1 þ 0:03282828283x2 þ 0:0004509379509x4 : Power series y3 ðxÞ can be transformed into Pade series (Fig. 3). r ¼ ½5=4 ¼ 1  0:4563492063x2 þ 0:02070105820x4 =1 þ 0:04365079365x2 þ 0:0008597883598x4 : Substitute step size h into y1 , y2 , y3 . We have y1 ðhÞ, y2 ðhÞ and y3 ðhÞ. We also have y10 ðhÞ, y20 ðhÞ and y30 ðhÞ. We show Table 1 for the solution of (5.1) by above numerical method. The numerical values on tables are coinciding with the exact solutions of (5.1).

6. Conclusion A Pade approximation method has proposed for solving differential–algebraic equations in this study. This method is very simple an effective for most of differential–algebraic equations. The computations associated with the example discussed above were performed by using Maple V. The proposed algorithm produced a rapidly convergent series. The Pade approximants that often show superior performance over series approximants provide a promising tool for using in applied fields (scientific applications).

References [1] K.E. Brenan, S.L. Campbell, L.R. Petzold, Numerical Solution of Initial-Value Problems in Differential–Algebraic Equations, North-Holland, New York, 1989. [2] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential– Algebraic Problems, Springer-Verlag, New York, 1991. [3] P. Deuflhard, E. Hairer, J. Zugck, One step and extrapolation methods for differential–algebraic systems, Numerical Mathematics 51 (1987) 501–516. [4] V.V. Murata, E.C. Biscaia, Structural and symbolic techniques for automatic characterization of differential–algebraic equations, Computers and Chemical Engineering 21 (1997) 829–834. [5] E. C ß elik, M. Bayram, Arbitrary order numerical method for solving differential–algebraic equation by Pade series, Applied Mathematics and Computation 137 (2003) 57–65. [6] E. C ß elik, E. Karaduman, M. Bayram, Numerical solutions of chemical differential–algebraic equations, Applied Mathematics and Computation 139 (2003) 259–264.