On the numerical solution of stiff systems

Applied Mathematics and Computation 170 (2005) 230–236 www.elsevier.com/locate/amc

On the numerical solution of stiff systems Nuran Guzel, Mustafa Bayram

*

Yildiz Teknik Universitesi, Fen-Edebiyat Fakultesi, Matematik Bolumu, 34210 Davutpasa-Esenler-Istanbul, Turkey

Abstract In this paper, we use power series method to solve stiff ordinary differential equations of the first order and an ordinary differential equation of any order by converting it into a system of differential of the order one. Theoretical considerations has been discussed and some examples were presented to show the ability of the method for linear and nonlinear systems of differential equations. We use MAPLE computer algebra systems for numerical calculations [G. Frank, MAPLE V, CRC Press Inc., Florida, 1996]. Ó 2005 Elsevier Inc. All rights reserved. Keywords: Stiff system; Power series; Numerical solution of the ordinary differential equations; MAPLE

1. Introduction A system 8 0 of first order differential equations can be considered as: y 1 ¼ f1 ðx; y 1 ; . . . ; y n Þ > > > < y 02 ¼ f2 ðx; y 1 ; . . . ; y n Þ .. > . > > : 0 y n ¼ fn ðx; y 1 ; . . . ; y n Þ *

Corresponding author. E-mail address: [email protected] (M. Bayram).

0096-3003/$ - see front matter Ó 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2004.11.035

ð1:1Þ

N. Guzel, M. Bayram / Appl. Math. Comput. 170 (2005) 230–236

231

with initial condition y i ðx0 Þ ¼ y 0;i ;

i ¼ 1; 2; . . . ; n;

where each equation represents the first derivative of one of the unknown functions as a mapping depending on the independent variable x, and n unknown functions f1, . . . , fn. Since every ordinary differential equation of order n can be written as a system consisting of n ordinary differential equation of order one, we restrict our study to a system of differential equations of the first order. That is f and y are vector functions for which we assumed sufficient differentiability [3,4,6,7]. The solutions of (1.1) can be assumed that y ¼ y 0 þ ex;

ð1:2Þ

where e is a vector function which is the same size as y0. Substitute (1.2) into (1.1) and neglect higher order term. We have the linear equation of e in the form Ae ¼ B; ð1:3Þ where A and B is a constant matrixes. Solve this equation of (1.3), the coefficients of x in (1.2) can be determined. Repeating above procedure for higher order terms, we can get the arbitrary order power series of the solutions for (1.1).

2. The method We define another type power series in the form f ðxÞ ¼ f0 þ f1 x þ f2 x2 þ    þ ðfn þ p1 e1 þ    þ pm em Þxn ;

ð2:1Þ

where p1, p2, . . . ,pm are constant. e1, e2, . . . ,em are basis of vector e, m is size of vector e. y is a vector with m elements in (1.2). Every element can be represented by the Power series in (2.1). y i ¼ y i;0 þ y i;1 x þ y i;2 x2 þ    þ ei xn ;

ð2:2Þ

where yi is ith element of y. Substitute (2.2) into (1.1), we can get the following fi ¼ ðfi;n þ pi;1 e1 þ    þ pi;m em Þxnj þ Qðxnjþ1 Þ;

ð2:3Þ

where fi is ith element of f(y, y 0 , x) in (1.1) and j is 0 if f(y, y 0 , x) have y 0 , 1 if do not. From (2.3) and (1.3), we can determine the linear equation in (1.3) as follow Ai;j ¼ P i;j ;

ð2:4Þ

Bi ¼ fi;n :

ð2:5Þ

232

N. Guzel, M. Bayram / Appl. Math. Comput. 170 (2005) 230–236

Solve this linear equation, we have ei (i = 1, . . . , m). Substitute ei into (2.2), we have yi (i = 1, . . . ,m) which are polynomials of degree n. Repeating this procedure from (2.2) to (2.5), we can get the arbitrary order Power series of the solution for differential-algebraic equations in (1.1). Let step size of x to be h and substitute it into the power series of y and y 0 , we have y and y 0 at x = x0 + h. If we repeat above procedure, we have numerical solution of differential-algebraic equations in (1.1) [1,2,5,8]. This method has been used to obtained approximate numerical and theoretical solutions of a large class of differential-algebraic equation [9–11] references therein.

3. Applications Example 1. The test problem consider the following differential equation system, y 01 ¼ 1002y 1 þ 1000y 22 ð3:1Þ y 02 ¼ y 1  y 2 ð1 þ y 2 Þ with initial condition y 1 ð0Þ ¼ 1 and y 2 ð0Þ ¼ 1: The exact solution is y 1 ðxÞ ¼ e2x and y 2 ðxÞ ¼ ex : From boundary condition, the solutions of (3.1) can be supposed as y 1 ðxÞ ¼ y 0;1 þ e1 x ) y 1 ðxÞ ¼ 1 þ e1 x; y 2 ðxÞ ¼ y 0;2 þ e2 x ) y 2 ðxÞ ¼ 1 þ e2 x: Substitute (3.2) into (3.1) and neglect higher order terms, we have e1 þ 2 þ QðxÞ ¼ 0; 1 þ e2 þ QðxÞ ¼ 0:

ð3:2Þ

ð3:3Þ

These formulas correspond to (2.3). The linear equation that corresponds to (2.4) can be given in the following Ae ¼ B;

ð3:4Þ

where 

1 A¼ 0

 0 ; 1



 2 B¼ ; 1



 e1 : e¼ e2

From Eq. (3.4), we have linear equation      2 1 0 e1 ¼ : 1 0 1 e2

N. Guzel, M. Bayram / Appl. Math. Comput. 170 (2005) 230–236

233

Solve this linear equation. We have   2 e¼ 1 and y 1 ðxÞ ¼ 1  2x; y 2 ðxÞ ¼ 1  x:

ð3:5Þ

From (3.5) the solutions of (3.1) can be supposed as y 1 ðxÞ ¼ 1  2x þ e1 x2 ;

ð3:6Þ

y 2 ðxÞ ¼ 1  x þ e2 x2 :

In this manner, substitute (3.6) into (3.1) and neglect higher order terms, we have ð2e1  4Þx þ Qðx2 Þ ¼ 0;

ð3:7Þ

ð2e2  1Þx þ Qðx2 Þ ¼ 0; where 

2 A¼ 0

 0 ; 2

  4 B¼ ; 1



 e1 e¼ : e2

From (3.7) we have linear equation      2 0 e1 4 ¼ : 0 2 e2 1 Solve this linear equation, we have   2 e¼ 1 : 2

Therefore y 1 ðxÞ ¼ 1  2x þ 2x2 ; 1 y 2 ðxÞ ¼ 1  x þ x2 : 2 Repeating above procedure we have 4 2 y 1 ¼ 1  2x þ 2x2  x3 þ x4    ; 3 3 1 2 1 3 1 4 y2 ¼ 1  x þ x  x þ x    : 2 6 24

ð3:8Þ

234

N. Guzel, M. Bayram / Appl. Math. Comput. 170 (2005) 230–236

The power series solution of the given Eq. (3.1) are coinciding with the exact solutions of (3.1). Example 2. Let us consider following system of differential equation y 01 ¼ 20y 1  0:25y 2  19:75y 3 ; y 02 y 03

y 1 ð0Þ ¼ 1;

¼ 20y 1  20:25y 2 þ 0:25y 3 ;

y 2 ð0Þ ¼ 0;

¼ 20y 1  19:75y 2  0:25y 3 ;

y 3 ð0Þ ¼ 1:

The analytic solution of the problem is y 1 ¼ ½e1=2t þ e20t ðcosð20tÞ þ sinð20tÞÞ =2; y 2 ¼ ½e1=2t  e20t ðcosð20tÞ  sinð20tÞÞ =2; y 3 ¼ ½e1=2t þ e20t ðcosð20tÞ  sinð20tÞÞ =2: If we apply the power series method to the given equation system, following solution is obtained, 1 3199 2 85333 3 3413333 4 33 x þ x  x  x5 y1 ¼ 1  x  4 16 32 256 250000 12799999 6 x ; þ 36 y2 ¼

79 3199 2 1 3386667 4 126221226 6 x x  x3 þ x  106666x5 þ x ; 4 16 96 256 355 81 3201 2 1 3413333 4 213333333 5 x x þ x3 þ x  x 4 16 96 256 2000 21688755 6 x : þ 61

y 3 ¼ 1 þ

The numerical results are illustrated in Tables 1–3. Example 3. As a second example, we consider a system representing a nonlinear reaction was taken from Hull [8]. Table 1 Comparison of theoretical and numerical values of the y1 in Example 1 t

Theoretical (y1)

Numerical (y1)

Error

0.00 0.02 0.04 0.06 0.08 0.1

1.0000000000 0.9342476725 0.8080890256 0.685011668 0.6128014753 0.6645036013

1.0000000000 0.9342476725 0.8080890254 0.6850115655 0.6128014680 0.6645035734

0.000000 0.000000 0.200E9 0.130E8 0.730E8 0.279E7

N. Guzel, M. Bayram / Appl. Math. Comput. 170 (2005) 230–236

235

Table 2 Comparison of theoretical and numerical values of the y2 in Example 1 t

Theoretical (y2)

Numerical (y2)

Error

0.00 0.02 0.04 0.06 0.08 0.1

0.0000000000 0.3168396725 0.4947663589 0.5716675667 0.5902094749 0.5978369336

0.0000000000 0.3168396744 0.4947664157 0.5716679544 0.5902109247 0.5978407963

0.00000000 0.19000E8 0.56800E7 0.38770E6 0.14498E5 0.38627E5

Table 3 Comparison of theoretical and numerical values of the y3 in Example 1 t

Theoretical (y3)

Numerical (y3)

Error

0.00 0.02 0.04 0.06 0.08 0.1

1.0000000000 0.6732101613 0.4854323144 0.3987779668 0.3705799643 0.3533924907

1.0000000000 0.6732101615 0.4854323234 0.3987780692 0.3705805394 0.3533946846

0.0000000 0.20000E9 0.90000E8 0.10240E6 0.57510E6 0.21939E5

dy 1 ¼ y 1 ; dt dy 2 ¼ y 1  y 22 ; dt dy 3 ¼ y 22 : dt The initial conditions are given by y1(0) = 1, y2(0) = 0 and y3(0) = 0. We have calculated first six terms of the power series as follows. 1 1 1 1 y 1 ¼ 1  x þ x2  x3 þ x4  x5    ; 2 3! 4! 5! 1 1 5 4 y 2 ¼ x  x2  x3 þ x4 þ x5    ; 2 3! 4! 5! 2 6 2 y 3 ¼ x3  x4  x5    : 3! 4! 5!

4. Conclusion The method has been applied to the solution of stiff systems of differential equations (nonlinear stiff and nonstiff equations system). A technique has also been considered to avoid the computation of consistent initial conditions. The

236

N. Guzel, M. Bayram / Appl. Math. Comput. 170 (2005) 230–236

numerical example has been presented to show that the approach is promising and the research is worth to continue in this direction. Power series method has proposed for solving differential equation systems in this study. This method is very simple and effective for most of differential equations system.

References [1] P. HenriciApplied Computational Complex Analysis, vol. 1, John Wiely & Sons, New York, 1974, Chapter 1. [2] G. Corliss, Y.F. Chang, Solving ordinary differential equations using Taylor series, ACM Trans. Math. Soft. 8 (1982) 114–144. [3] K.E. Brenan, S.L. Campbell, L.R. Petzold, Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations, North-Holland, Amsterdam, 1989. [4] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II: Stiff and DifferentialAlgebreis Problems, Springer-Verlag, 1991. [5] W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Numerical Recipes, Cambridge University Press, Cambridge, 1988. [6] U.M. Ascher, L.R. Petzold, Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2003. [7] P. Amodio, F. Mazzia, Numerical solution of differential-algebraic equations and computation of consistent initial/boundary conditions, J. Comput. Appl. Math. 87 (1997) 135–146. [8] T.E. Hull, W.H. Enright, B.M. Fellen, A.E. Sedgwick, Comparing numerical methods for ordinary differential equations, SIAM J. Numer. Anal. 9 (1972) 603. [9] G. Frank, MAPLE V, CRC Press Inc., Florida, 1996. [10] E. C ¸ elik, E. Karaduman, M. Bayram, Numerical method to solve chemical differentialalgebraic equations, Int. J. Quant. Chem. 89 (5) (2002) 447–451. [11] E. C ¸ elik, M. Bayram, Arbitrary order numerical method for solving differential-algebraic equation by Pade´ series, Appl. Math. Comput. 137 (2003) 57–65.