Journal of Computational and Applied Mathematics 184 (2005) 464 – 474 www.elsevier.com/locate/cam A numerical algorithm for some singularly perturbed boundary value problems Basem S. Attili∗ UAEU, Department of Mathematics and Computer Science, College of Science, P.O. Box 17551, Al-Ain, United Arab Emirates Received 4 December 2003; received in revised form 29 August 2004 Abstract A numerical algorithm is proposed to solve singularly perturbed linear two-point value problems. The method starts with a partial decoupling of the system to obtain two independent subsystems, fast and slow components. Each subsystem is then solved separately. A second-order finite difference scheme is used for this purpose. Numerical examples will be presented to show the efficiency of the method. © 2005 Elsevier B.V. All rights reserved. Keywords: Singularly perturbed systems; Two point B.V.P’s; Finite differences 1. Introduction Singularly perturbed boundary value problems are common in applied sciences. They often occur in, for example, optimal control, [24] and convection–diffusion equations. The presence of the perturbation parameter leads to difficulties when numerical techniques are used to solve such problems and convergence will not be uniform. This is since boundary layers appear in these problems, see [17,23–25]. One of the simplest examples of problems of such a type in a one-dimensional case is the problem of stationary diffusion process with a reacting substance that has the form d2 u(t) − c(t)u(t) = f (t), dt 2 Subject to u(t) = (t),
2 t ∈ (0, 1) ∗ Tel.: +9713 7134286; fax: +97 3767 1291. E-mail address: b.attili@uaeu.ac.ae. 0377-0427/$ - see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2005.01.021 (1.1) B.S. Attili / Journal of Computational and Applied Mathematics 184 (2005) 464 – 474 465 where c(t)
c0 > 0 and
is a small parameter that characterizes the diffusion coefficient of the substance matter. When
tends to zero, boundary layers appear in a neighborhood of the boundary. If we consider a special case of (1.1) where c(t) = 1 and f (t) = −1 with some boundary conditions; that is, (1.1) becomes d2 u(t) − u(t) = −1, t ∈ (0, 1) dt 2 Subject to u(0) = u(1) = 0.
2 (1.2) The exact solution to (2.1) has the form u(t) = u(t;
) = 1 − e−t/
+ e−(1−t)/
. 1 + e−1/
(1.3) Note that 0  u(t)  1 and lim
−→0 u(1/2) = 1 and the maximum value of u(t) occurs at t = 21 . Using the second-order central difference approximation to the second derivative of the form t¯t (z(t)) = [(t − t )z(t)] , h (1.4) where t z(t) = (z(t + h) − z(t)) h t z(t) = (z(t) − z(t − h)) . h and It is known, see [14,15], that the error in this formula (in the l − ∞ norm) depends on the parameter
and the step size h; that is |u(t) − z(t)|  K(
).h2 . (1.5) Moreover, for sufficiently small values of
; that is, for
=
(h) = 1/ h, the maximum error becomes larger than some positive constant K̄, see [28,29] max |u(t) − z(t)|
K > 0 D (1.6) as h −→ 0. The estimate (1.5) shows that the error tends to zero and hence the scheme is given by
2 t¯t z(t) − z(t) = −1; z(0) = z(1) = 0, where t¯t is as given by (1.4) converges as h −→ 0, while (1.6) shows that this difference scheme does not converge uniformly with respect to the small parameter
. Even in the case when only the approximate solution is required, finite difference schemes and finite element methods produced unsatisfactory results, see [21,26,27]. It was shown in [29] that the results of using classical methods are also unsatisfactory even when a very fine grid is used. This suggests having numerical methods where the error in the approximate solution tends to zero as h −→ 0 (N −→ ∞) independently of the parameter
; that is, uniform convergence is desired. Some special schemes can be found in for example Il’in [15] where 466 B.S. Attili / Journal of Computational and Applied Mathematics 184 (2005) 464 – 474 special finite difference methods are used. Fitted schemes that allow the use of the meshes with an arbitrary distribution of nodes can be found in [1,2,6,11,15]. Problems that arise in viscoelastic flow were considered by Allen and Southwell [2], Attili [3], Davies et al. [5] and Karageorghis et al. [16]. Specially fitted schemes were used by Allen and Southwell [2] and some kind of regularization of the singularity was suggested by others before applying special schemes of boundary value solvers. In this work, we will consider a more general problem rather than (1.1), which was used as an illustration. The approach we are going to follow starts with decoupling of the system (1.1) followed by applying a difference scheme of order 2. This will be done in Section 3. Some preliminaries will be given in Section 2. Also, in Section 3, we will describe the algorithm used to solve the resulting system. Finally, numerical examples in comparison with their exact solutions will be presented in Section 4. 2. Preliminaries We will consider a general framework which is a system of first-order differential equations that has the form dx1 = A11 (t)x1 + A12 (t)x2 + F1 (t), dt dx2 = A21 (t)x1 + A22 (t)x2 + F2 (t),
dt (2.1) where 0 < t < 1, and subject to the linear boundary conditions of the form L(x1 (0), x1 (1), x2 (0), x2 (1)) = 0. Here xi ∈ R ni , i = 1, 2,
> 0 is a small parameter, Aij (t) are matrices with appropriate dimension and F1 (t), F2 (t) are given real vector functions of appropriate dimension. It is assumed that the functions Aij (t), F1 (t) and F2 (t) together with L are such that the system given by (2.1) has a unique bounded solution. This system can be considered as a control stochastic linear system where F1 and F2 represent the control. The analytical aspects of systems like (2.1) through asymptotic analysis are well established, see [17,23–25]. Studies conclude that the solution has two states; namely, one varies slowly and the other rapidly with respect to t. Also, the presence of the small parameter
leads to what is called boundary layer phenomenon where the solution does not converge uniformly at end points. Similar behavior occurs when such nonuniformity occurs at an interior point called interior layer phenomenon, see [9,10,13,18,19]. As a result of such boundary layer behavior and the fast–slow states of the solutions, the two point boundary value problem given by (2.1) cannot be solved accurately using typical and usual numerical techniques. Different suggestions and approaches were used to overcome this problem. For example, Mattheji [22] considered partial decoupling to (2.1) before applying a difference scheme. Fitted methods on arbitrary distribution of nodes were considered by Il’in [15] while adaptive methods were considered in a way such that the nodes are distributed to achieve parameter uniform convergence. B.S. Attili / Journal of Computational and Applied Mathematics 184 (2005) 464 – 474 467 3. Decoupling of the system As mentioned before, the solution of the system given by (2.1) comprises of two states one varying slowly and the other rapidly with respect to the independent variable t. Due to the nonuniform behavior and the slow–fast property of the solution, partial decoupling is needed. The purpose of such a decoupling is to transform the system in order to obtain two independent subsystems, one for the slow and the other for the fast state. Before presenting the partial decoupling and for what follows, |v| denotes the supremum norm for vector v, while for a matrix |G| denotes the vector induced supremum norm and G = supt∈[0,1] |G(t)|, where G(t) is a vector or a matrix function. For our purpose and to obtain the decoupling, we will introduce the variable x3 (new slow variable) given by x3 = x1 −
Nx2 . (3.1) If A22 is an invertible matrix, we may associate a system of lower dimension of the form x3′ = (A11 − NA21 )x3 + (F1 − NF2 ), (3.2) N(t) = N0 (t) +
N1 (t) + O(
2 ) (3.3) where −1 with N0 = A12 A−1 22 and N1 = ((A11 − N0 A21 )N0 − Ṅ0 )A22 continuously differentiable matrix functions with the dot represents differentiation, see [9]. As a result, we have the following theorem. Theorem 1. Assume that the matrix functions Aij ; i, j = 1, 2 are continuously differentiable such that Aij (t)  bij and Ȧij (t)  cij where bij and cij are positive constants. If any eigenvalue (t) of A22 (t) satisfies the conditional stability conditions Re((t))  − c3 or Re((t))
c4 with c3 and c4 positive constants. Then there exists an
0 > 0 and conditionally differentiable matrix functions N(t), N0 (t) and N1 (t) given by (3.3) such that for 0 <
<
0 , the system (2.1) is transformed to x3′ = (A11 − NA21 )x3 + (F1 − NF2 ),
x2′ = A21 x3 + (
A21 N + A22 )x2 + F2 . (3.4) Similar results can be found in [30] and a more advanced and detailed one can be found in [9]. For the proof and details, the reader can refer to these references. Now the system given by (3.4) is essential for the numerical treatment of the problem given by (2.1). The first part of (3.4) can be solved for x3 using any two point boundary value problem solver provided that n1 appropriate linearly independent boundary conditions are available. Once x3 is computed, the second part of (3.4) can be solved for x2 . Note that this part can be written in the form
x2′ = F (t)x2 + g(t), 0 < t < 1, (3.5) where F (t) = (
A21 N + A22 ) and g(t) = A21 x3 + F2 . For the decoupling of the boundary conditions, note that (3.5) can be written in the form
d(N −1 x2 ) = (D −
N −1 Ṅ )(N −1 x2 ) + N −1 g dt (3.6) 468 B.S. Attili / Journal of Computational and Applied Mathematics 184 (2005) 464 – 474 with D = N −1 ((
A21 N + A22 )N −
Ṅ), where D = diag(D1 , D2 ). With G = N −1 g = N −1 (A21 x3 + F2 ), we can split (3.6) into
u̇1 = D1 u1 + G1 ,
u̇2 = D2 u2 + G2 , (3.7) where x2 = (u1 , u2 ). Leading to the decoupling of the boundary conditions into L(z(0), z(1), u1 (0), u1 (1), u2 (0), u2 (1)) = 0. It should be noted that one can do without separating the boundary conditions. Instead of such separation, we can use the uniform asymptotic expansion of the system given in (2.1). Then numerically compute the dominant part of the fundamental solution matrix, which in turn can be introduced in the boundary conditions too obtain the initial value of the solution of the given boundary value problem. Details of such uniform expansions of systems of the form given by (2.1) can be found in [7,8]. We will treat (3.5) numerically using an appropriate difference scheme that converges uniformly to the solution x2 assuming there are n2 linearly independent boundary conditions. For that purpose, let hj be a grid defined on [0, 1] with h = maxj hj . Here t0 = 0, tN = 1 and tk = kh; k = 0, 1, . . . , N. Let zk = z(tk ) be the numerical approximation to the solution at the grid points and similarly Fk = F (tk ), gk = g(tk ) and x2k = x2 (tk ). The error in the approximation will be given by ek = zk − x2k . With this notation, we will have Euler’s forward implicit scheme given by
 zk+1 − zk
= Fk+1 zk+1 + gk+1 (3.8) hk and Euler’s backward explicit scheme is given by 


zk+1 I zk = − gk , Fk + hk hk (3.9) where I is the identity matrix with appropriate order. These schemes given by (3.8) and (3.9) can be shown to have a quadratic rate of convergence as stated in the following result. Theorem 2. Assume F (t) and g(t) are twice continuously differentiable. Assume also that F (t) , F −1 (t) and g(t) are bounded and the eigenvalues of F (t) are such that Re((t))  − c3 or Re((t))
c4 with c3 and c4 positive constants. Then if hk is chosen such that 

, 2
F (t) , (3.10) min hk > max k F (t) then the schemes given by (3.8) and (3.9) solves (3.5) with zk bounded and the error ek is O(h2 ) if h F (t) ?1. (3.11) The proof of the theorem can be obtained from a result given by Kreiss et al. [20] and so it is omitted. Note that there is an assumption on the minimum of the grid length given by (3.10) although in general one B.S. Attili / Journal of Computational and Applied Mathematics 184 (2005) 464 – 474 469 requires h <
as a precondition on the grid to make any standard method accurate. This latter assumption will produce large round-off errors and that is why (3.10) is needed. Also that for the proof, condition (3.10) is a strong one and a weaker assumption can do the job; that is, instead of (3.10), the proof follows with the assumption 2
hk −1 < Fk+1 < . hk
(3.12) Remark 1. If (3.11) is not satisfied, the error in the methods will be O(h). In this case the differentiability conditions on F (t) and g(t) can be relaxed. It is well known that when the eigenvalues of F (t) are all negative (positive) then the forward (backward) integration is stable, see [4]. This suggests the use of some transformation on (3.5) in order for F (t) to be block-diagonalized; that is, we require the system (3.5) to be transformed through the use of N(t) (whose existence and uniqueness is guaranteed, see [12]) to
d(N −1 x2 ) = (D −
N −1 Ṅ )(N −1 x2 ) + N −1 g, dt (3.13) where D(t) = N −1 F (t)N (t) = diag(D1 , D2 ) with Re((D1 (t)) < 0 and Re((D2 )) > 0. Note that the stability of the system given by (3.13) will not change since the conditional stability of F (t) implies that of D and (D −
N −1 Ṅ), see [18, Chapter 5]. The same can be said about A22 and
A21 N + A22 . The system matrix F (t) has a special block structure. It is important for it to remain so throughout the whole interval. If such a structure changes, one must find subintervals (ai , ai+1 ); i = 0, 1, . . . , n − 1 on which the block structure does not alter. This means different equations of the form given by (3.13) needs to be solved on each subinterval and ai and ai+1 ; i = 0, 1, . . . , k has to be among the grid points. These details will be combined in the following algorithm: Algorithm 3. Step I: Compute N0 (t) = A12 A−1 22 , N1 (t) = ((A11 − N0 A21 )N0 − Ṅ0 )A−1 22 . Step II: Set N(t) = N0 (t) +
N1 (t). Step III: Form the transformed system (3.4). Step IV: For the first part of (3.4) (i) Determine the partition points a0 , a1 , . . . , an−1 . (ii) On each subinterval, use the QR-algorithm to compute the block-diagonalizing transformation N0 , N1 , . . . , Nn−1 . Step V: Compute (i) Di = Ni−1 ((A22 +
A21 N)Ni −
Ṅi ), where Di = diag(D1i , D2i ) with Re((D1i (t))) < 0 and Re((D2i (t))) > 0. (ii) Gi = (GT1i , GT2i )T = Ni−1 (A2i z + g2 ). (iii) ui = (uT1i , uT2i )T = Ni−1 u. 470 B.S. Attili / Journal of Computational and Applied Mathematics 184 (2005) 464 – 474 Step VI: On each subinterval split the second part of (3.4) into
u̇1i = D1i u1i + G1i ,
u̇2i = D2i u2i + G2i . (3.14) Step VII: Find the grid points on each subinterval by monitoring (3.11) and (3.12) hij 2
< |Di−1 (ti,j +1 )| < hij
and hij |Di (ti,j +1 )| > 1, where ti,j +1 refers to the (j + 1)st grid point on the ith subinterval. Step VIII: Combine the discritized form of both (3.13) and the first part of (3.4) and solve using 
wi,j +1 − wi,j +1 = K1 K2 [wi,j +1 − wi,j ] + K1 [Hi,j +1 − Hi,j ],
hi,j where w = [z, y1 , y2 ] the jth numerical approximation at the ith step and  1  − 2 Im 0 0 K1 = with n1 + n2 = n, 0 0n1 0 0 0 In 2   0 0 A11 − NA21 K2 = 0 D1i 0 0 0 D2i subject to the boundary conditions L(z(0), z(1), y1i (0), y1i (1), y2i (0), y2i (1)) = 0. 4. Numerical results In this section, we will illustrate the techniques developed in the previous sections by some numerical examples. Two examples from O’Mally [24] will be considered. One of the examples exhibits boundary layers at the end points and the other at an interior point. Example 1. Consider the second-order linear equation
ẍ − 2t ẋ = 0, x(−1) = −1, −1 < t < 1, x(1) = 2. (4.1) The solution for −1 < t < 1 is given by 1 t 2 −1 e
+ O(1) t and graphically for
= 0.01 by Fig. 1. Transform (4.1) into a first-order system by letting ẋ = x1 and ẍ = ẋ1 to obtain x(t) = 2 + ẋ = x1 ,
ẋ1 = 2tx 1 . (4.2) B.S. Attili / Journal of Computational and Applied Mathematics 184 (2005) 464 – 474 471 Fig. 1. Example 1, exact solution with
= 0.01. Fig. 2. Example 1, numerical solution with
= 0.01. Table 1
= 0.01 Example 1 Example 2 a1 −0.98 −0.05
= 0.0001 a2 0.98 0.05 a1 −0.995 −0.003
= 0.000001 a2 0.998 0.0037 a1 −0.999 −0.0017 a2 0.999 0.0076 subject to same boundary conditions. According to the notation of (2.1) we have A11 = 0, A12 = 1, f1 (t) = 0, A21 = 0, A22 = 2t, f2 (t) = 0. As a result the matrices used in (3.4) will be N0 = A12 A−1 22 = 1/t, −1 3 N1 = ((A11 − N0 A21 )N0 − Ṅ0 )A22 = −1/4t . The choice of stretching points for Example 1 as well as Example 2 are given in Table 1. 472 B.S. Attili / Journal of Computational and Applied Mathematics 184 (2005) 464 – 474 Fig. 3. Example 1, numerical solution with
= 0.0001. Fig. 4. Example 1, stability condition violated. Fig. 5. Example 2, numerical solution with
= 0.01. The numerical solution for
= 0.01 and
= 0.0001 are given in Figs. 2 and 3, respectively. When violating condition 3.10 resulted in the solution given in Fig. 4. B.S. Attili / Journal of Computational and Applied Mathematics 184 (2005) 464 – 474 473 Fig. 6. Example 2, numerical solution with
= 0.0001. Example 2. The second example has the form
ẍ + 2t ẋ − 2x = 0, x(−1) = −1, −1 < t < 1, x(1) = 2. (4.3) Rewriting (4.3) as a first-order system leads to ẋ = x1 ;
ẋ1 = −2tx 1 + 2x. (4.4) subject to same boundary conditions as in Example 1. Again the stretching points are given in Table 1 and the solution for
= 0.01 and
= 0.0001 are given in Figs. 5 and 6, respectively. References [1] M.V. Alexeevski, On difference scheme for a singularly perturbed parabolic equations, Uchen. Zap. (Scientific Papers), Natural Sci. 1 (1980) 3–13. [2] D.N. Allen, R.V. Southwell, Relaxation method applied to determine the motion, in two dimensions, of viscous fluid past a fixed cylinder, Quart. J. Mech. Appl. Math. 8 (1955) 129–145. [3] B.S. Attili, Initial value methods for the primary two point boundary value problems in modeling viscoelastic flows, Internat. J. Comput. Math. 74 (2000) 379–390. [4] W.A. Coppel, Stability and Asymptotic Behavior of Differential Equations, Heath, Boston, 1965. [5] A.R. Davies, A. Karageorghis, T.N. Phillips, Spectral Galerkin methods for the primary two point boundary value problem in modeling viscoelastic flows, Internat. J. Numer. Methods Engrg. 26 (1988) 647–662. [6] E.P. Doolan, J.H. Miller, W.A. Schilders, Uniform Numerical Methods for Problems With Initial and Boundary Layers, Boole Press, Dublin, 1980. [7] V. Dragan, A. Halanay, Uniform asymptotic expansions for the fundamental matrix singularly perturbed linear systems and applications, in: F. Verhulst (Ed.), Asymptotic Analysis II, Lecture Notes in Mathematics, vol. 985, Springer, Berlin, 1983, pp. 215–245. [8] V. Dragan, A. Halanay, Stabilization of Linear Systems, Birkhauser, Basel, 1999. [9] V. Dragan, T. Morozan, Exponential stability for a class of singularly perturbed Ito differential equations, EJQTDE, Proceedings of the Sixth Colloquium QTDE, vol. 7, 2000, pp. 1–13. [10] V. Dragan, T. Morozan, Exponential stability for a class of linear time-varying singularly perturbed stochastic systems, Dynam. Contin. Discrete Impulsive Systems Series B 9 (2002) 213–233. 474 B.S. Attili / Journal of Computational and Applied Mathematics 184 (2005) 464 – 474 [11] K.V. Emel’yanov, Difference scheme for a three-dimensional elliptic equation with small parameters at highest-order derivatives, Kraevye Zadachi Dlya Uravn. Mat. Fiz. (Boundary Value Problems for Equations of Mathematical Physics) (1973) 30–42. [12] H. Gingold, A method of global block-diagonalization for matrix valued functions, SIAM J. Math. Anal. 9 (1978) 1076. [13] E.Ya. Gorelova, Stability of singularity perturbed stochastic systems, Automat. Remote Control 58 (1997) 112–121. [14] P.W. Hemker, G.I. Shishkin, Discrete approximation of singularly perturbed parabolic PDE’s with a discontinuous initial condition, in: J.J. Miller (Ed.), Bail VI Proceedings, Boole Press, 1992, pp. 3–4. [15] M.A. Il’in, Difference scheme for a differential equation with a small parameter at the highest-order, Mat. Zametki 6 (1969) 237–248. [16] A. Karageorghis, T.N. Phillips, A.R. Davies, Spectral collocation methods for the primary two point boundary value problem in modeling viscoelastic flows, Internat. J. Numer. Methods Engrg. 26 (1988) 805–813. [17] J. Kevorkian, J.D. Cole, Perturbation Methods in Applied Mathematics, Springer, Berlin, 1981. [18] P.V. Kokotovic, H.K. Khalil, J. O’Reilly, Singular Perturbation Methods in Control Analysis and Design, Academic Press, New York, 1986. [19] P. Kokotovic, H.K. Khalil, J. O’Reilly, Singular Perturbations Methods in Control Analysis and Design, SIAM Appl. Math. (1999). [20] H.O. Kreiss, N.K. Nichols, D.I. Brown, Numerical methods for stiff two-point boundary value problems, SIAM J. Numer. Anal. 23 (1986) 325. [21] G.I. Marchuk, Methods of Numerical Mathematics, Springer, New York, 1982. [22] R.M. Mattheji, On approximating smooth solution of linear singularly perturbed ODE, in: P.W. Hemker, J.J.H. Miller (Eds.), Numerical Analysis of Singular Perturbation Problems, Academic Press, New York, 1979, p. 456. [23] A.H. Nayfeh, Introduction to Perturbation Techniques, Wiley, New York, 1981. [24] R.E. O’Mally Jr., Introduction to Singular Perturbation, Academic Press, New York, 1974. [25] R.E. O’Mally Jr., Singular perturbation and optimal control, in: Mathematical Control Theory, Lecture Notes in Mathematics, vol. 680, Springer, New York, 1978, p. 170. [26] A.A. Samarski, E.S. Nikolaev, Methods of Solving a Grids Equations, Nauka, Moscow, 1978. [27] A.A. Samarski, Theory of Difference Schemes, Nauka, Moscow, 1980. [28] G.I. Shishkin, Approximation of solutions to singularly perturbed boundary value problems with corner boundary layer, Doklady Akad. Nauk SSSR 296 (1987) 39–43. [29] G.I. Shishkin, Grid approximation of singularly perturbed boundary value problems or quasi-linear parabolic equations in the case of complete degeneracy in spatial variables, Sov. J. Numer. Anal. Math. Modelling 6 (1991) 243–261. [30] A. Singh, M.K. Kadalbajoo, Partial decoupling of slow and fast variables, J. Math. Anal. Appl. 155 (1991) 46.