Sinc-galerkin solution for nonlinear two-point boundary value problems with applications to chemical reactor theory

Available online at www.sciencedirect.com SCIENCe

ELSEVIER

#DIRECT e

MATHEMATICAL AND COMPUTER MODELLING

Mathematical and Computer Modelling 42 (2005) 1237-1244 www.elsevier.com/locate/mcm

S i n c - G a l e r k i n S o l u t i o n for N o n l i n e a r T w o - P o i n t B o u n d a r y Value P r o b l e m s w i t h A p p l i c a t i o n s to C h e m i c a l R e a c t o r T h e o r y A. SAADATMANDI Department of Applied Mathematics Amirkabir University of Technology, Tehran, Iran M . RAZZAGHI* Department of Mathematics and Statistics Mississippi State University, Mississippi State, MS 39762, U.S.A. and Department of Applied Mathematics Amirkabir University of Technology, Tehran, Iran razzaghi~ math.msstate, edu. M. DEHGHAN Department of Applied Mathematics Amirkabir University of Technology, Tehran, Iran

(Received October 2003; revised and accepted April 2005) Abstract--The Sinc-Galerkin method is presented for solving nonlinear two-point boundary value problems for second order differential equations. A problem arising from chemical reactor theory is then considered. Properties of the Sinc-Galerkin method are utilized to reduce the computation of nonlinear two-point boundary value problems to some algebraic equations. The method is computationally attractive and applications are demonstrated through an illustrative example. @ 2005 Elsevier Ltd. All rights reserved. Keywords--Sinc,

Galerkin, Two-point, Chemical, Reactor.

1. I N T R O D U C T I O N In this paper, we first consider the nonlinear differential equations

au" + bu' + a (x, u (x)) = 0,

(1)

~o~ (0) + c1~' (0) = o,

(2)

dou (1) + d]u' (1) = O,

(3)

with boundary conditions

*Author to whom all correspondence should be addressed. The work of the first author was supported by Iran Telecommunication Research Center. 0895-7177/05/$ - see front matter @ 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.mcm.2005.04.008

Typeset by fit.&~-TEX

A. SAADATMANDI et

1238

al.

where b, co, cl, do, dz, a ~ 0 are given constants and G is an analytic function. We then solve a problem which falls into this category from chemical reactor theory. The mathematical model for an adiabatic tubular chemical reactor which processes an irreversible exothermic chemical reaction. For steady state solutions, this model can be reduced to the ordinary differential equation given in [1] by u" - At' + A# (/3 - u) exp (u) = O,

(4)

with boundary conditions u' (0) = Au (0),

u' (1) = O.

(5)

The unknown u represents the steady state temperature of the reaction, and the parameters A, #, and/3 represent the Peclet number, the Damkohler number and the dimensionless adiabatic temperature rise, respectively. This problem has been studied by several authors [2-4] who have demonstrated numerically the existence of solutions (sometimes multiple solutions), for particular ranges. A numerical solution of (4) with boundary conditions in (5) is reported in [1]. In [1], by first using Green's function technique, the problem is converted into a Hammerstein integral equation and then the solution is obtained by using Adomian's method. Sinc methods have increasingly been recognized as powerful tools for attacking problems in applied physics and engineering [5]. The books [5,6] provide excellent overviews of methods based on Sinc functions for solving ordinary and partial differential equations and integral equations. Sine methods have also been employed as forward solvers in the solution of inverse problems [7,8]. In [9], a Sinc collocation method for solving linear two-point boundary value problems for secondorder differential equations with mixed boundary conditions is presented and in [10], a SineGalerkin technique is applied to calculate the wind-driven current in a sea. There are several reasons to approximate by Sine functions. First, they are easily implemented and give good accuracy for problems with singularities; approximation by Sinc function are typified by errors of the form O ( e x p ( - c / h ) ) , where c > 0 is a constant and h is a step size. Second, approximation by Sine functions handles singularities in the problem. The effect of any such singularities will appear in some form in any scheme of numerical solution, and it is well known that polynomial methods do not perform well near singularities. Finally , these kinds of approximation yield both an effective and rapidly convergent scheme for solving the problem and circumvents the instability problems that one typically encounters in some difference methods [11]. In the present paper, we first apply the Sinc-Galerkin method for solving (1)-(3) and then solve (4),(5). Our method consists of reducing the solution of (1)-(3) to a set of algebraic equations by expanding the u(z) as Sine function with unknown coefficients. The properties of Sine function are then utilized to evaluate the unknown coefficients. The paper is organized as follows. Section 2 is devoted to the basic formulation of the Sine function required for our subsequent development. In Section 3 the Sinc-Galerkin method is used to approximate the solution for (1)-(3). In Section 4, we report our numerical finding for (4),(5) and demonstrate the accuracy of the proposed numerical scheme by considering a numerical example.

2. S I N C

FUNCTION

PROPERTIES

Sinc function properties are discussed thoroughly in [5,6]. In this section an overview of the basic formulation of the Sinc function required for our subsequent development is presented. 2.1. T h e Sinc F u n c t i o n The Sine function is defined on the whole real line, - o c < z < cx~, by Sine (z) =

sin (rcz) - - , rrz 1,

z¢O,

z=0.

Sinc-Galerkin Solution

1239

For h > 0, and k = 0, ±1, + 2 , . . . , the translated Sinc functions with evenly spaced nodes are given by sin [Tr/h (z -- kh)]

S ( k , h ) ( z ) = S i n c ( Z - ~ kh)- -

7r/h(z-kh)

' z#kh, z kh.

1,

=

The Sine functions form an interpolatory set of functions, i.e.,

S(k,h)(jh)=dkj =

1, k = j , O, k7 ~j.

If a function f(z) is defined on the real axis, then for h > 0 the series

C(f,h)(z)=

~

f(kh) Sine

k=--oo

is called the Whittaker cardinal expansion of f whenever this series converges. The properties of Whittaker cardinal expansion have been extensively studied in [61. This properties are derived in the infinite strip Ds of the complex w-plane, where for d > 0,

D s = { w = t +is : [s] < d <_ 2 } . Approximations can be constructed for infinite, semi-infinite and finite intervals. To construct approximations on the interval (0, 1), which is used in this paper, the eye-shaped domain in the z-plane

DE=

z=x+iy:


,

is mapped conformally onto the infinite strip Ds via

This is shown in Figure 1. iS W-plane

T III1~11111

/lll,~lllll , I I I l ll l l l l T,~

x

I

I I I I I II

DE

I I I

I

Figure 1. Relationship between the domains

DE

and

~"

Ds

Ds.

The basis functions on (0, 1) are taken to be the composite translated Sinc functions,

Sk(z)- S(k,h)o¢(z) = Sinc(¢(Z)hkh ) ,

z E DE,

(6)

where S(k, h) e ¢(z) is defined by S(k, h)(¢(z)). The inverse map of w = ¢(z) is z=q5 -l(w)=~(w)-

exp (w) l+exp(w)"

Thus, we may define the inverse images of the real line and of the evenly spaced nodes {kh}k~__~ as

F = { ¢ ( t ) E DE : --oc < t < oc} and

ekh

zk = ~ (kh) - 1 + ekh' respectively.

k = O, -t-1, +2,...,

1240

A. SAADATMANDI et al.

2.2. Sinc F u n c t i o n Interpolation and Quadrature The class of functions known as exponential error estimates exist for infinite Sinc interpolation and quadrature are denoted by B(DE) and is now defined. DEFINITION. Let B(DE) be the class of functions F which are analytic in

DE ,

satisfy

IF (z) dz I ~ O, t :i:oc, (t+L) where L = {iv : ]vI < d ~_ 7r/2}, and on the boundary of DE, 'denoted ODE), satisfy ¢

N (F) = /

IF (z) dzl < 0 0 .

Jo D~ Interpolation and quadrature rules for function in B(DE) are defined in the following theorems whose proof are found in [5,6]. THEOREM 1. If O'F C B(DE) then, for all z C F, F(z)=

where

~

F(zj)S(j,h)o¢(z)+EF,

sin (~¢ (z)/h) ~

¢' (~) F (~) d~

THEOREM 2. IY F E B(DE) and there are positive constants a and M so that

F(z)

¢'(z)

-< M e x p ( - a I¢(z)]) '

for all z G F then, N

fF F ( z ) d z = h

~

F (zj) ¢'(zj----~+ O (exp ( - a N h ) ) + 0 (exp (-2~rd/h))

(7)

j=-N

Hence, if h = V / ~ / a N , O(exp(-(27vdaN)l/2) ).

the exponential order of the Sinc trapezoidal quadrature rule is

COROLLARY. An important special case in (7) occurs when the integrand has the form H(z)S (/, h) o¢(z). Due to the interpolation

S (l,h) o ¢(zj) = S (l,h) (jh) = 0,

if 1 ¢ j,

the Sinc quadrature rule is a weighted point evaluation of the order of the method ~v H (zl) H (z) S (l, h) o ¢ (z) dz = h ~ + 0 (exp (-27rd/h)).

(8)

As seen in Theorem 1, the Sinc approximations are derived using the residue theorem with contour integrals where the integrands are independent. Given these error terms, the above expressions show Sinc interpolation and quadrature on B(DE) converge exponentially [10]. The Sinc-Galerkin method requires derivatives of composite Sinc functions evaluated at the nodes. The expressions required for the present discussion are [5] 5 (°) [S(j,h) oO(z)][z=zk { 1, j = k , J'k= = 0, j ~ k , (~(1) d [S (j, h) o ¢ (z)] 1 / 0, j,k = d-¢ ~=~k = ~ [. (--1)k-J

~(2) "j,k

d~ (z)] ~ zk 1 db 2 [S(j,h) o ¢ = =~

/-

j = k, j ¢ k,

3 ' - 2 ( - 1 ) k-j

(k - y)~

j=k, j~k,

Sinc-Galerkin Solution

3.

THE

1241

SINC-GALERKIN

METHOD

For the boundary conditions in (2),(3), the Sinc basis functions in (6) do not have a derivative when z tends to 0 and 1. Thus, we modify the Sine basis functions as

S(k,h) o¢(z) ¢,(z) Now, the derivatives of the modified Sine basis functions are defined as z approaches 0 or 1. We also add boundary basis functions that are cubic polynomials. These polynomials are given by

uo(z) = (3e, - c0)(1 - z) 2 + (co - 2cl)(1 - z) 3, ul (z) = ( d 0 + 3 d l ) z 2 - ( d o + 2 d l ) z 3. Note that the derivatives at z = 0 and 1 for these basis functions are defined, and also satisfy boundary conditions in (2),(3). For the purpose of illustrating the exposition of our method we define

Lu (x) - au" (x) + bu' (x),

(9)

then, (1) is now given by Lu (x) = - a

(x, u (x)).

In the Sinc-Galerkin technique, the approximate solution for u(x) conditions in (2),(3) is represented by

(z) =

+

in (1) subject to boundary

(z),

(10)

where U b (Z) =

(11)

C__N__lU 0 (Z) -Jr-C N + l U 1 (Z)

and

N

Zg

k=-N

S(k,h) o¢(z) =z (l-z)

¢'(z)

E

ckS(k,h) o ¢ ( z ) .

(12)

k=-N

The approximate solution in (10) satisfy the boundary conditions in (2),(3). The 2 N + 3 coefficients {ek}~=_N_l N+ 1 are determined by orthogonalizing the residual Lua(z) + G(z, Ua(Z)) with respect to the Sine basis functions {sj}N+2N_ 1 in (6). The inner product is defined by

(f, g} =

j•01

f (z) g (z) w (z) dz,

where w(z) is a weight function given by [10] w (z) = 1 / v / ~ ( z ) = v/z (1 - z). So, we have the discrete Sinc-Galerkin system as

(Lua +G(z, ua),Sj) = 0 ,

j=-N-1,...,N

+ I.

This leads to

c N-1 + (Luh, Sj} +CN+l = - < G ( z , Ua),Sj)

j=-N-1,...,N

+ I.

(13)

1242

A. SAADATMANDI

et al.

With the Sinc quadrature rule in (8) we evaluate the inner products in (13). First, we evaluate (Luo, Sj) and {Lul, $5} as 1

(Luo, Sj} =

~0

Luo (z) Sj (z) dz = h Luo (zj) ~/¢, (z) (4, (zj))~/~

(14)

and

(nul, Sj) =

j•O 1

Sj (z) dz

Lul (z) ~

Lul (zj)

(15)

= h (~)t (zj))3/2"

The nodal points are zj = eJh/(1 + edh). The numerators of the right-hand side of (14),(15) are

Luo = 3b (2cl - co) z 2 + (-6aco + 12acl - 6bcl + 4bco) z - 6acl + 4aco -bco, Lul = -3b (2dl + do) z 2 + (-6ado - 12adl + 2bdo + 6bdl) z ÷ 2a (do ÷ 3dl). By using (9), the inner product {Luh, Sj} in (13) is given by

(Luh, Sy) =

Lun (z) Sj (z) w (z) dz

= a

/01

]l

~,~ (z) Sj (z) w (z) dz + b

/01

(16) I

~ (z) Sj (z) ~ (z) dz.

Integrating by parts from the first integral in the right-hand side of (16) we get

/o I

/o 1

l, lzl lzll z

Since w(O) = w(1) = O, the first term in the right-hand side of (17) is zero, also by integrating by parts from the second term we obtain

1 ~'h (z) ~zd (sj (z) ~ (z)) dz = ~ (z) Gd (Sj (z) w (z)) o -

uh (z) ~

(Sj (z) w (z)) dz.

(18)

Using (12) we get uh(O) = uh(1) = O, we also have

± (sj (4 w (z)) = d s ~ (~) ¢' (z) ~ (z) + s5 (z) ~' (z) dz = (dsj(z)+~ d~2 (Sj (z) ~ (z)) = ~

- d g S j (z) (¢' (z)) ~ ~ (z) +

Sj (z)

1 (1 -

(19)

2z) Sj (z)) ~ ,

(z) ~ (~) + Sj (z) ~' (z)

Sj (~) (¢" (z) ~ (z) + 2¢' (z) ~' (z)) + S, (z) ~" (z) (2O)

= (~¢-S¢2Sj(z)-- 1S j

.

From (17), (18), and (20) we get f01 u~ (z) Sj (z) w (z) dz = 9~ol uh (z) ( - d2 ~Sj

(z) - ~1S j (z) ) (¢' (z)) 3/2 dz

(21)

Integrating by parts from the second integral in the right-hand side of (16) we get

~01 ~',~ (z) sj (z) ~ (z) d~ = [~,~(~) Sj (z) ~ (z)] 1 - fO 1 ~,~ (~) d ( s j ( z ) w ( z ) ) dz.

(22)

Sinc-Galerkin Solution

1243

Using (19),(22) we obtain

L 1~'. (~) sj (z) ~(~) dz = - L 1~h (~) (d

1(l_2z)~j (z))~7(z)dz"

Sj (z) + 7

(23)

From (16), (21), and (23), we get

(Luh, S;) = a L i~h (z)

(d~ ~-g~s; (z) - is 4 ; (~)) (¢' (z))~l~ dz -b

/o I ~ ( ~ )

S ; ( z ) + 1 (i - 2z) Sj (z)) X/~(z) dz. 7

Applying, the Sinc quadrature rule, we have

N
/ 5j,k (2) - (1/4) 5j,(0)k

N bh E

ck I

k=-N

a5(2)

h,g 1)

¢' (zk) ~t~

¢' (zk) 312

j,k

Q,~(1) ek

(1/2) (1 _ 2zk) 5~02)

~j,k + . . . . .

~9' (Zk) 3/2

(24)

''j,k

4¢ t (Zk) 1/2 Jr- 2(y

(~k)3/--''~

Furthermore, the inner product (G(z, Ua), Sj) is given by

(~ (Z, Ua),Sj) = h ~(zj'ua (zj)) (zj)) ~/~ '

(25)

(¢'

where Ua (zj)

cj = I ub (zj) + ¢' (zj~'

j

- N , . . . , N,

j = - N - 1, N + 1.

t Ub (Zj), Substituting (14)-(25) in (13) we get a(~(2) Luo(zj) + ~N ck j,k C-N-1 (¢t(Zj))3/2 (~t k=-N (Zk) 112

(1) bSj,k ¢, (Zk) 312

Lul (zj)

-[-CN+I (¢,

(Zj)) 3/2 --

j =-N-1,...,N

( 4¢' a b(1 - 2zk)~ 5(o)1 + --J.~] (Zk) 1/2 2~9! (Zk) 312 )

(26)

a (zj, Ua (zj)) ¢, (Zj) 3/2

'

+ I.

(26) generates a set of 2N + 3 nonlinear algebraic equations which can be solved for the unknown coefLients cj in (11),(12) by using Newton's method. Consequently ua(z) given in (10) can be calculated. 4. N U M E R I C A L

EXAMPLE

To validate the application of Sinc-Galerkin method to (4),(5), we use particular values of the parameters, )~ = 10, /3 = 3, and # = 0.02. For such values for the parameters, a unique solution is guaranteed by the contraction mapping principle [1]. For this problem, by using the contraction mapping principle, the required integrations cannot be done analytically and are evaluated numerically using the trapezoidal rule. For Sinc-Galerkin method in the absence of additional information it is usually safe to take a = 1 and d _< Tr/4. These are conservative choices which work reasonable well in most cases [121. Here, we choose c~ = 1, d = 7c/4, N = 10, and N = 20 and also truncate the numerical results after the sixth decimal points. Table 1 gives a comparison of the results from the contraction

1244

A. SAADATMANDI

et al.

Table 1. Results for Example 1. N = Sinc-Galerkin 10 N=20

x

Contraction Principle

Shooting Method

Method in [1]

0.0

0.006079

0.006048

0.006048

0.006049

0.006048

0.2

0.018224

0.018192

0.018192

0.018197

0.018192

0.4

0.030456

0.030424

0.030424

0.030437

0.030424

0.6

0.042701

0.042669

0.042669

0.042649

0.042669

0.8

0.054401

0.054371

0.054371

0.054383

0.054371

1.0

0.061459

0.061458

0.061458

0.061459

0.061458

mapping principle, the shooting method, the Adomian's method and present method with N = 10 and N = 20. As shown in Table 1, the results using present method with N = 20 agree with those of the shooting method and Adomian's method up to the sixth decimal place. The results from the contraction mapping principle agree to at least three decimal place, but, these will involve errors due to the numerical integration procedure used. 5.

CONCLUSION

The Sinc-Galerkin method is used to solve nonlinear two-point boundary value problems with applications to chemical reactor theory. Properties of the Sinc-Galerkin method are utilized to reduce the computation of this problem to some algebraic equation. The method is computationally attractive and applications are demonstrated through an illustrative example.

REFERENCES 1. N.M. Madbouly, D.F. McGhee, and G.F. Roach, Adomian's method for Hammerstein integral equations arising from chemical reactor theory, Applied Mathematics and Computation. 117, 341-249, (2001). 2. A. Poore, A tubular chemical reactor model, A Collection of Nonlinear Model Problems Contributed to the Proceeding of the A M S - S I A M , 28-31, (1989). 3. R. Heinemann and A. Poore, Multiplicity stability and oscillatory dynamics of the tubular reactor, Chemical Engineering Science 36, 1411-1419, (1981). 4. R. Heinemann and A. Poore, The effect of activiation energy on tubular reactor multiplicity, Chemical Engineering Science 37, 128 131, (1982). 5. F. Stenger, Numerical Methods Based on Sine and Analytic Functions, Springer-Verlag, New York, (1993). 6. J. Lund, and K. Bowers, Sine Methods for Quadrature and Differential Equations, SIAM, Philadelphia, PA, (1992). 7. J. Lund and C. Vogel., A Fully-Galerkin method for the solution of an inverse problem in a parabolic partial differential equation, numerical solution of an inverse, Inverse Problems 6~ 205-217, (1990). 8. R. Smith and K. Bowers, A Sinc-Galerkin estimation of diffusivity in parabolic problems, Inverse Problems 9, 113 135, (1993). 9. B. Bialecki, Sinc-collocation methods for two-point boundary value problems, I M A J. Numer. Anal 11, 357-375, (1991). 10. D.F. Winter, K. Bowers, and J. Lund, Wind-driven currents in a sea with a variable Eddy viscosity calculated via a Sinc-Galerkin technique, Internt. J. Numer. Methods Fluids 33, 1041-1073, (2000). 11. M. S. Sababheh, and A.M.N. Al-khaled, Some convergence results on Sinc interpolation, J. Inequal. Pure and Appl. Math 4, 32 48, (2003). 12. J.L. Mueller and T.S. Shores, A new Sinc-Calcrkin method for convection-diffusion equations with mixed boundary conditions, Computers Math. Applic. 47 (4/5), 803-822, (2004).