- Email: [email protected]
A Wavelet-Galerkin method for a singularly perturbed convection-dominated diffusion equation
Applied Mathematics and Computation 181 (2006) 1635–1644 www.elsevier.com/locate/amc
A Wavelet-Galerkin method for a singularly perturbed convection-dominated diffusion equation Mohamed El-Gamel Department of Mathematical Sciences, Faculty of Engineering, Mansoura University, Egypt
Abstract There are few techniques available to numerically solve singularly perturbed parabolic problems. In this paper we show that the Wavelet-Galerkin method is a very effective tool in numerically solving such problems. The method is then tested on several examples and a comparison with the method of reduction order is made. It is shown that the Wavelet-Galerkin method yields better results. 2006 Elsevier Inc. All rights reserved. Keywords: Wavelet-Galerkin method; Singularly perturbed; Convection-dominated diffusion equation
1. Introduction Wavelet analysis is a new numerical concept which allows one to represent a function in terms of a set of basis functions, called wavelets, which are localized both in location and scale. As we noted earlier, spectral bases are infinitely differentiable, but have global support. On the other hand, basis functions used in finite difference or finite element methods have small compact support, but have poor continuity properties. As a result, spectral methods have good spectral localization, but poor spatial localization, while finite difference and finite element methods have good localization, but poor spectral localization. Wavelet bases seem to combine the advantages of both spectral and finite difference (or finite element) bases. One can expect that numerical methods based on wavelet bases are able to attain good spatial and spectral resolution. In wavelet applications to the solution of partial differential equations the most frequently used wavelets are those with compact support introduced by Daubechies [4]. Exploration of usage of Daubechies wavelets to solve partial differential equations has been undertaken by a number of investigators such as [1,2,9,11,13,18–20,25]. Additional references can be found in the recent review by El-Gamel [7,8]. The main objective of this paper is to present a new numerical approach for solving the singularly perturbed convection-dominated diffusion equation
E-mail address: [email protected] 0096-3003/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.03.017
1636
M. El-Gamel / Applied Mathematics and Computation 181 (2006) 1635–1644
ou o2 u ou ¼ 2 þ gðx; tÞ; ot ox ox
ð1:1Þ
0 < x < 1; t > 0;
subject to the boundary conditions uð0; tÞ ¼ 0;
uð1; tÞ ¼ 0;
ð1:2Þ
and the initial condition uðx; 0Þ ¼ u0 ðxÞ;
Pm
ð1:3Þ i
where (0 < 1) and g(x, t) = h(t)f(x), where f ðxÞ ¼ i¼0 ai x , is a polynomial of degree m in x, otherwise, we approximate f by such a polynomial if necessary. Singularly perturbed problems appear in many branches of engineering, such as fluid mechanics, heat transfer, and problems in structural mechanics posed over thin domains. For more details see [15,22]. Theorems that list conditions for the existence and uniqueness of solutions of such problems are thoroughly discussed in [22]. This kind of problem has been investigated by many researchers [3,6,10,12,14,21,23]. The organization of the paper is as follows. In Section 2, as a background, the fundamental properties of wavelet functions and connection coefficients are described. The Wavelet-Galerkin method is developed for second-order singularly perturbed boundary value problems with homogeneous boundary condition in Section 3. In Section 4, we show how wavelets are used to solve singularly perturbed convection-dominated diffusion equation. Some numerical examples are presented in Section 5. Finally, Section 6 provides conclusions of the study. 2. Orthonormal wavelet bases 2.1. Multiresolution in L2(R) A multiresolution analysis of L2(R) is defined [4] as a sequence of closed subspaces Vj with the following properties: (1) (2) (3) (4)
Vj Vj+1 "j 2 Z, f ðxÞ 2 V j () f ð2xÞ 2 V jþ1 , f ðxÞ 2 V 0 () f ðx þ 1Þ 2 V 0 , ¨j2ZVj is dense in L2(R); ˙jVj2Z = {0}.
The existence of a scaling function /(x) is required of which the translates generate a basis in each Vj via V j ¼ spanf/ji gi2Z with /ji ¼ 2j=2 /ð2j x iÞ;
j; i 2 Z:
In the classical case this basis is orthonormal, so that h/ji ; /jk iR ¼ dik ; with hf ; giR ¼
Z
ð2:1Þ
1
f ðxÞ gðxÞdx; 1
being the usual inner product. Let the Wj denote a subspace complementing the subspace Vj in Vj+1, i.e., V jþ1 ¼ V j W j : Each element of Vj+1 can be uniquely written as the sum of an element in Vj and an element in Wj, which contains the details needed to pass from an approximation at level j to an approximation at level j + 1.
M. El-Gamel / Applied Mathematics and Computation 181 (2006) 1635–1644
1637
Based on the function /(x) one can find w(x), the so-called mother wavelet, of which the translates and dilates constitute orthonormal bases of the spaces Wj W j ¼ spanfwji gi2Z ; generated by the wavelets wji ¼ 2j=2 wð2j x iÞ;
j; i 2 Z:
Each function f 2 L2(R), can now be expressed as 1 X X X cj0 i /j0 i ðxÞ þ d ji wji ðxÞ; f ðxÞ ¼ j¼j0
i2Z
ð2:2Þ
i2Z
where cji ¼ hf ; /ji iR
d ji ¼ hf ; wji iR :
Of course, in numerical applications the sum in (2.2) are truncated which corresponds to the projection of f into a subspace of VJ L2(R). The decomposition (2.2) is orthogonal, as, by construction, hwji ; wlk iR ¼djl dik ; hwji ; /lk iR ¼0;
j P l:
in addition to (2.1). 2.2. Regularity and local decay of wavelet coefficients It is well known that the local or global regularity of a function is closely related to the decay of its wavelet coefficients. The latter directly determines the error being made when truncating a wavelet sum at some scale. Depending on the type of norm and whether global or local characterization is considered, various relations of the this kind have been developed. See [4,16] for an overview. Let us recall here only briefly as an example the case of an a-Lipschitz function taken from [5]. Suppose f 2 L2(R), then for [a, b] R the function f is a-Lipschitz for any x0 2 [a, b]. i.e., a
jf ðx0 þ hÞ f ðx0 Þj 6 Cjhj ; if and only if there exists a constant A such that jhf ; wji ij 6 A2jað1=2Þ ; for any (j, i) with i/2j 2 [a, b]. 2.3. Connection coefficients The compactly supported wavelets of Daubechies [4] have proven to be a useful tool in the numerical solutions of partial differential equations [1,7–9,11,13,18–20,25] for recent numerical solutions of partial differential equations using Wavelet-Galerkin techniques. In each of the above papers the notion of connection coefficients has played a key role in the approximations of partial derivatives as well as nonlinear terms in the wavelet discretization of the differential equations. Connection coefficients and some of their properties and algorithms for computing them are described in [17]. For solving Eq. (1.1) in the bounded interval [0, 1], we need to calculate the following two connection coefficients: Z x M mk ðxÞ ¼ y m /ðy kÞdy; ð2:3Þ 0 Z x Cnk ðxÞ ¼ /ðnÞ ðy kÞ/ðyÞdy: ð2:4Þ 0
1638
M. El-Gamel / Applied Mathematics and Computation 181 (2006) 1635–1644
The computation of the connection coefficient (2.3) and (2.4) is reported by Chen et al. [17] with the exception that the normalizing conditions are replaced by " # m1 X X n K n Cnk ðmÞ ¼ n!H1 ðmÞ ð1Þ K n Cnk ðL 1Þ : 2L6k6mLþ1
k¼mLþ2
3. Singularly perturbed second-order boundary value problem In this section, we shall study Wavelet-Galerkin scheme for the singularly perturbed boundary value problem u00 þ au0 þ bu ¼ f ðxÞ
ð3:1Þ
for 0 < x < 1;
subject to boundary conditions uð0Þ ¼ 0;
uð1Þ ¼ 0;
ð3:2Þ
where (0 < 1), and a, b are constants and f(x) is a polynomial of degree any order in x, in [0, 1] otherwise, we approximate them by such a polynomial if necessary. Let the solution uJ(x) of the problem be approximated by its Jth level wavelet series on the interval (0, 1), i.e., J 2X 1
uJ ðxÞ ¼
uk /Jk ðxÞ;
k 2 Z;
ð3:3Þ
k¼2L
where /Jk(x) = 2J/2/(2Jx k), J > 0, and uk, k = 2 L, 3 L, . . . , 2J 1 are 2J + L 2 unknown coefficients to be determined. Here, the integer J is used to control the smoothness of the solution. The larger integer J is used, the more accurate a solution can be obtained. However, the number of equations required to solve the unknown coefficients is increased. In (3.3), the parameter L implies that the wavelet associated with the set of L Daubechies filter coefficients is used as the solution bases. Substituting the wavelet series approximation uJ(x) in Eq. (3.3) for u(x) in Eq. (3.1), yields
J 2X 1
J J 2X 1 2X 1 d2 d / / ðxÞ þ a u ðxÞ þ b uk /Jk ðxÞ ¼ f ðxÞ: k Jk Jk 2 dx dx k¼2L k¼2L
uk
k¼2L
ð3:4Þ
To determine the coefficient uk, we take the inner product of both sides of Eq. (3.4) with /Jl, J J J Z 1 Z 1 Z 1 2X 1 2X 1 2X 1 uk /00Jk ðxÞ/Jl ðxÞdx þ a uk /0J ;k ðxÞ/Jl ðxÞdx þ b uk /Jk ðxÞ/Jl ðxÞdx 0
k¼2L
¼
Z
k¼2L
f ðxÞ/Jl ðxÞdx;
0
We assume that f ðxÞ ¼ J 2X 1
0
k¼2L
0
1
cJkl uk
k¼2L
þa
l ¼ 2 L; 3 L; . . . ; 2J 1:
Pm
J 2X 1
l¼0 ai x
bJkl uk
k¼2L
i
ð3:5Þ
, is a polynomial of degree m in x. We write the equation in (3.5) as
þb
J 2X 1
aJkl uk ¼ d Jml ;
l ¼ 2 L; 3 L; . . . ; 2J 1;
ð3:6Þ
k¼2L
where cJkl ¼
Z
1
/00Jk ðxÞ/Jl ðxÞdx; 0
aJkl ¼
Z
0
bJkl ¼
Z
1
/0J ;k ðxÞ/Jl ðxÞdx;
ð3:7Þ
0 1
/Jk ðxÞ/Jl ðxÞdx
ð3:8Þ
M. El-Gamel / Applied Mathematics and Computation 181 (2006) 1635–1644
1639
and d Jml ¼
Z
1
0
m X
ai xi /Jl ðxÞdx ¼
i¼0
m X i¼0
ai 2
ðiþ12ÞJ
M il ð2J Þ:
ð3:9Þ
Eq. (3.6) can be further put into the matrix-vector form A1 U ¼ D;
ð3:10Þ
where A1 ¼ C þ aB þ bA; C ¼ ½cJkl 2L6k;l62J 1 ; A¼
½aJkl 2L6k;l62J 1 ;
B ¼ ½bJkl 2L6k;l62J 1 ; D¼
½d Jml 2L6l62J 1
ð3:11Þ ð3:12Þ
and U ¼ ½u2L ; u3L ; . . . ; u2J 1 s ; where s denotes the matrix transpose. Now we have a linear system of 2J + L 2 equations of the 2J + L 2 unknown coefficients. We can obtain the coefficients of the approximate solution by solving this linear system. The solution U gives the coefficients in the Wavelet-Galerkin approximation uJ(x) of u(x). 4. Singularly perturbed convection-dominated diffusion equation In this section, we will solve the problem given by Eqs. (1.1)–(1.3). Using the finite difference method and setting ut ¼
ui ui1 ; Dt
then, Eq. (1.1) becomes ui ui1 d2 ui dui ¼ 2 þ gðx; ti Þ: Dt dx dx Rewrite this equation as d2 ui dui 1 i 1 i1 u ¼ u gðx; ti Þ: 2 Dt Dt dx dx
ð4:1Þ
Let the solution uJ(x, ti) of the problem be approximated by its J-th level wavelet series on the interval (0, 1), i.e., uJ ðx; ti Þ ¼
J 2X 1
uik /Jk ðxÞ;
k 2 Z;
ð4:2Þ
k¼2L
Substituting the wavelet series approximation uJ(x, ti) for ui in Eq. (4.1), yields 2X J J J 2X 1 2X 1 1 d2 i d i 1 1 i1 i u u / ðxÞ / ðxÞ u / ðxÞ ¼ u gðx; ti Þ: k Jk k Jk k Jk 2 dx Dt k¼2L Dt dx k¼2L k¼2L To determine the coefficient uik , we take the inner product of both sides of Eq. (4.3) with /Jl, 2X J J J 2X 1 2X 1 1 1 00 0 i i uk h/Jk ðxÞ; /Jl ðxÞi uk h/Jk ðxÞ; /Jl ðxÞi uik h/Jk ðxÞ; /Jl ðxÞi Dt k¼2L k¼2L k¼2L i1 u ¼ h1; /Jl ðxÞi hðti Þhf ðxÞ; /Jl ðxÞi; l ¼ 2 L; 3 L; . . . ; 2J 1: Dt
ð4:3Þ
ð4:4Þ
1640
M. El-Gamel / Applied Mathematics and Computation 181 (2006) 1635–1644
Using the notations defined in (3.7)–(3.9), we write (4.4) as 2X J J J i1 2X 1 2X 1 1 1 u J i J i J i ckl uk bkl uk akl uk ¼ d J0l þ hðti Þd Jml : Dt Dt k¼2L k¼2L k¼2L
ð4:5Þ
The above equations can be further put into the matrix-vector form HU ¼ G;
ð4:6Þ
where
1 A H¼CB Dt
and i1 u G¼ Djm¼0 þ hðti ÞD; Dt and A, B, C and D are defined by Eqs. (3.11) and (3.12). Eq. (4.6) is a linear system of (2J + L 2) equations in (2J + L 2) unknown coefficients. By solving this system, we obtain an approximate solution at resolution level J at time t = ti. Repeating these steps until we get desired time. If we substitute uik ¼ uk ðti Þ in Eq. (4.2) we obtain an approximate solution to the problem (1.1)–(1.3) at time t = ti. 5. Numerical examples The four examples included in this section were selected in order to illustrate the performance of the Wavelet-Galerkin method in solving singularly perturbed problems. For the sake of comparison only, we consider the second and third examples discussed by Reddy [21], who used the method of reduction of order to obtain his numerical solution, and the last example discussed by Kang [24], who used the finite-difference streamlinediffusion method to obtain his numerical solution. In the Wavelet-Galerkin solutions, Daubechies 6 wavelets are used because they give sufficiently accurate results and each example was run for a sequence of J = 5, 6, 7, 9. We use absolute error which is defined as absolute error ¼ juexact solution U Wavelet-Galerkin j:
Example 5.1. Consider the boundary value problem
d2 u þ u ¼ x x2 2; dx2
0
and uð0Þ ¼ 0;
uð1Þ ¼ 0:
The exact solution is given by uðxÞ ¼ xð1 xÞ: By taking different values of at level J = 9 the maximum absolute error is displayed in Table 1. Example 5.2 [21]. Consider the boundary value problem
d2 u du ¼ 1 þ 2x; þ dx2 dx
and uð0Þ ¼ 0;
uð1Þ ¼ 1:
0
M. El-Gamel / Applied Mathematics and Computation 181 (2006) 1635–1644
1641
Table 1 Maximum absolute error for Example 5.1
Maximum absolute error
102 104 105 107 109 1011
1.8966E007 1.7915E009 2.0134E010 1.7957E012 3.6970E014 2.9587E014
Table 2 Maximum absolute error for Example 5.2
Method of reduction of order [21]
Wavelet-Galerkin method
103 104
0.171E04 0.136E04
0.452E07 0.105E07
The exact solution is given by uðxÞ ¼ xðx þ 1 2Þ þ
ð2 1Þð1 ex= Þ : 1 e1=
Table 2 exhibits a comparison between the maximum absolute errors obtained by using the Wavelet-Galerkin method with the analogous results of Reddy [21], who used the method of reduction of order for solving this problem. Table 3 exhibits the exact solution, the Wavelet-Galerkin solution and the absolute errors at level J = 9. Example 5.3 [21]. Consider the boundary value problem
d2 u du u ¼ 0; þ dx2 dx
0
and uð0Þ ¼ 1;
uð1Þ ¼ 1:
The exact solution is given by uðxÞ ¼
ðem2 1Þem1 x þ ð1 em1 Þem2 x ; em2 em1
where m1 ¼
1 þ
pffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ 4 ; 2
and
m2 ¼
1
pffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ 4 : 2
Table 3 The exact solution, the Wavelet-Galerkin solutions, and the absolute errors for Example 5.2 at level J = 9 x
Exact solution
Wavelet-Galerkin solution
Absolute error 1.0E07
0.0 1/8 2/8 3/8 4/8 5/8 6/8 7/8 1.0
0.0 0.85920 0.68735 0.48425 0.2499 0.0157 0.31255 0.64065 1.0
0.0 0.8592000 0.6873500 0.4842500 0.2499000 0.01569999 0.31254999 0.64064999 1.0
0.0 0.1 0.0 0.0 0.0 0.1 0.1 0.1 0.0
1642
M. El-Gamel / Applied Mathematics and Computation 181 (2006) 1635–1644
Table 4 exhibits a comparison between the maximum absolute errors obtained by using the Wavelet-Galerkin method with the analogous results of Reddy [21], who used the method of reduction of order for solving this problem. Table 5 exhibits the exact solution, the Wavelet-Galerkin solution and the absolute errors at level J = 9. Example 5.4 [24]. Consider the boundary value problem ou o2 u ou ¼ 2 þ gðx; tÞ; ot ox ox
0 < x < 1; t > 0
and uð0; tÞ ¼ 1 þ
1 et= ; 1 e1=
uð1; tÞ ¼ 1
and uðx; 0Þ ¼ 1; where gðx; tÞ ¼
½t2 t þ ð1 xÞ ð1xÞt= e : ð1 e1= Þ
The exact solution is given by uðx; tÞ ¼ 1 þ
1 eð1xÞt= : 1 e1=
Table 6 exhibits a comparison between the maximum errors obtained by using the Wavelet-Galerkin method with the analogous results of Kang [24], who used the finite-difference streamline-diffusion method at time t = 0.995. All computations associated with the above examples were performed by using MATLAB. Table 4 Maximum absolute error for Example 5.3
Method of reduction of order [21]
Wavelet-Galerkin method
103 104
0.1763E03 0.1920E04
0.1789E07 0.4083E08
Table 5 The exact solution, the Wavelet-Galerkin solutions, and the absolute errors for Example 5.3 at level J = 9 x
Exact solution
Wavelet-Galerkin solution
Absolute error 1.0E08
0.0 1/8 2/8 3/8 4/8 5/8 6/8 7/8 1.0
1.0 0.416898489 0.472401974 0.535294876 0.606560980 0.687315047 0.778820249 0.882507932 1.0
1.0 0.416898485 0.472401971 0.535294874 0.606560979 0.687315046 0.778820249 0.882507931 1.0
0.0 0.4 0.3 0.2 0.1 0.1 0.0 0.1 0.0
Table 6 Maximum absolute error for Example 5.4 at time t = 0.995
Finite-difference method [24]
Wavelet-Galerkin method
102
0.0124
0.543E04
M. El-Gamel / Applied Mathematics and Computation 181 (2006) 1635–1644
1643
6. Conclusions We have presented and illustrated the Wavelet-Galerkin method for singularly perturbed boundary value problems. An obvious advantage of this procedure is that it loads Daubechies’s constants and calculates the scaling functions, the connection coefficients as well as the rest of the constants only once. This leads to a considerable saving of the computational time and improves numerical results through the reduction of round-off errors. We have implemented the present method on four examples by taking different values of . Computational results are presented in tables. Here we have given results for only few values, although the solutions are computed at all points. The approximate solution is compared with the solution of the method of reduction of order. It can be observed from the results that the Wavelet-Galerkin method yields better results, which shows the efficiency of the method. Acknowledgement The author would like to express sincere thanks to Prof. A. Mohsen at Cairo University, for helpful discussions over the manuscript. References [1] A. Avudainayagam, Wavelet-Galerkin method for integro-differential equations, Appl. Numer. Math. 32 (2000) 247–254. [2] A. Avudainayagam, C. Vani, Wavelet-Galerkin solutions of quasilinear hyperbolic conservation equations, Commun. Numer. Meth. Eng. 15 (1999) 589–601. [3] T. Aziz, A. Khan, A spline method for second-order singularly perturbed boundary-value problems, J. Comput. Appl. Math. 147 (2002) 445–452. [4] I. Daubechies, Ten Lectures on Wavelets, Captial City Press, Vermont, 1992. [5] J. Frohich, K. Schneider, An adaptive wavelet-Galerkin algorithm for one-and two-dimensional flam computations, Eur. J. B/Fluids 13 (1994) 439–471. [6] M. El-Gamel, J.R. Cannon, On the solution of second order singularly-perturbed boundary value problem by the Sinc-Galerkin method, Z. Angew. Math. Phys. 56 (2005) 45–58. [7] M. El-Gamel, A. Zayed, A comparison between the Wavelet-Galerkin and the Sinc-Galerkin methods in solving nonhomogeneous heat equations, in: Inverse Problem, Image Analysis, and Medical Imaging, in: Zuhair Nashed, Otmar Scherzer (Eds.), Contemporary Mathematics of the American Mathematical Society, vol. 313, AMS, Providence, 2002, pp. 97–116. [8] M. El-Gamel, Wavelet algorithm for the numerical solution of nonhomogeneous time-dependent problems, Int. J. Differ. Eqn. Appl. 9 (2004) 169–185. [9] R. Glowiniski, W.M. Lawton, M. Ravachol, E. Tenenbaum, Wavelet solutions of linear and nonlinear elliptic, parabolic and hyperbolic problems in one space dimension, Comput. Meth. Appl. Sci. Eng. (1990) 55–120 (chapter 4). [10] S. Gonzalez-Pinto, L. Casasu´s, P. Gonza´lez-Vera, A numerical scheme to approximate the solution of a singularly perturbed nonlinear differential equation, J. Comput. Appl. Math. 35 (1991) 217–225. [11] S. Jaffard, Wavelet methods for fast resolution of elliptic problems, SIAM J. Numer. Anal. 29 (1992) 965–986. [12] L. Jiang, X. Yue, Local exponentially fitted finite element schemes for singularly perturbed convection–diffusion problems, J. Comput. Appl. Math. 132 (2001) 277–293. [13] F. Jin, T.Q. Ye, Instability analysis of prismatic members by wavelet-Galerkin method, Adv. Eng. Softw. 30 (1999) 361–367. [14] M. Kadalbajoo, A. Rao, The alternating group explicit (AGE) method for singularly perturbed boundary value problems, Appl. Math. Comput. 68 (1995) 125–142. [15] J. Kevorkian, J. Cole, Perturbation Methods in Applied Mathematics, Springer-Verlag, 1981. [16] S. Mallat, A Wavelet Tour of Signal Processing, Academic Press, New York, 1999. [17] M.Q. Chen, C.H. Wang, Y.P. Shin, The computation of wavelet-Galerkin approximation on a bounded interval, Int. J. Numer. Meth. Eng. 39 (1996) 2921–2944. [18] M.Q. Chen, C.H. Wang, Y.P. Shin, A Wavelet-Galerkin method for solving population balanced equations, Comput. Chem. Eng. 20 (1996) 131–145. [19] S. Qian, J. Weiss, Wavelets and the numerical solution of boundary value problems, Appl. Math. Lett. 6 (1993) 47–52. [20] S. Qian, J. Weiss, Wavelets and the numerical solution of partial differential equations, J. Comput. Phys. 106 (1993) 155–175. [21] Y. Reddy, P. Chakravathy, Method of reduction of order for solving singularly perturbed two point boundary value problem, Appl. Math. Comput. 136 (2003). [22] H.G. Roos, M. Stynes, L. Tobiska, Numerical Methods for Singularly Perturbed Differential Equations, Springer-Verlag, 1996. [23] A. Singh, A numerical method for singularly perturbed systems of linear two point boundary value problems using partial decoupling, J. Comput. Appl. Math. 40 (1992) 55–62.
1644
M. El-Gamel / Applied Mathematics and Computation 181 (2006) 1635–1644
[24] T. Kang, D. Yu, Some a posteriori error estimates of the finite-difference streamline-diffusion method for convection-dominated diffusion equations, Adv. Comput. Math. 15 (2001) 193–218. [25] Jin-Chao Xu, Wei-Chang Shann, Galerkin-Wavelet methods for two-point boundary value problems, Numer. Math. 63 (1992) 123– 144.
A Wavelet-Galerkin method for a singularly perturbed convection-dominated diffusion equation Mohamed El-Gamel Department of Mathematical Sciences, Faculty of Engineering, Mansoura University, Egypt
Abstract There are few techniques available to numerically solve singularly perturbed parabolic problems. In this paper we show that the Wavelet-Galerkin method is a very effective tool in numerically solving such problems. The method is then tested on several examples and a comparison with the method of reduction order is made. It is shown that the Wavelet-Galerkin method yields better results. 2006 Elsevier Inc. All rights reserved. Keywords: Wavelet-Galerkin method; Singularly perturbed; Convection-dominated diffusion equation
1. Introduction Wavelet analysis is a new numerical concept which allows one to represent a function in terms of a set of basis functions, called wavelets, which are localized both in location and scale. As we noted earlier, spectral bases are infinitely differentiable, but have global support. On the other hand, basis functions used in finite difference or finite element methods have small compact support, but have poor continuity properties. As a result, spectral methods have good spectral localization, but poor spatial localization, while finite difference and finite element methods have good localization, but poor spectral localization. Wavelet bases seem to combine the advantages of both spectral and finite difference (or finite element) bases. One can expect that numerical methods based on wavelet bases are able to attain good spatial and spectral resolution. In wavelet applications to the solution of partial differential equations the most frequently used wavelets are those with compact support introduced by Daubechies [4]. Exploration of usage of Daubechies wavelets to solve partial differential equations has been undertaken by a number of investigators such as [1,2,9,11,13,18–20,25]. Additional references can be found in the recent review by El-Gamel [7,8]. The main objective of this paper is to present a new numerical approach for solving the singularly perturbed convection-dominated diffusion equation
E-mail address: [email protected] 0096-3003/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.03.017
1636
M. El-Gamel / Applied Mathematics and Computation 181 (2006) 1635–1644
ou o2 u ou ¼ 2 þ gðx; tÞ; ot ox ox
ð1:1Þ
0 < x < 1; t > 0;
subject to the boundary conditions uð0; tÞ ¼ 0;
uð1; tÞ ¼ 0;
ð1:2Þ
and the initial condition uðx; 0Þ ¼ u0 ðxÞ;
Pm
ð1:3Þ i
where (0 < 1) and g(x, t) = h(t)f(x), where f ðxÞ ¼ i¼0 ai x , is a polynomial of degree m in x, otherwise, we approximate f by such a polynomial if necessary. Singularly perturbed problems appear in many branches of engineering, such as fluid mechanics, heat transfer, and problems in structural mechanics posed over thin domains. For more details see [15,22]. Theorems that list conditions for the existence and uniqueness of solutions of such problems are thoroughly discussed in [22]. This kind of problem has been investigated by many researchers [3,6,10,12,14,21,23]. The organization of the paper is as follows. In Section 2, as a background, the fundamental properties of wavelet functions and connection coefficients are described. The Wavelet-Galerkin method is developed for second-order singularly perturbed boundary value problems with homogeneous boundary condition in Section 3. In Section 4, we show how wavelets are used to solve singularly perturbed convection-dominated diffusion equation. Some numerical examples are presented in Section 5. Finally, Section 6 provides conclusions of the study. 2. Orthonormal wavelet bases 2.1. Multiresolution in L2(R) A multiresolution analysis of L2(R) is defined [4] as a sequence of closed subspaces Vj with the following properties: (1) (2) (3) (4)
Vj Vj+1 "j 2 Z, f ðxÞ 2 V j () f ð2xÞ 2 V jþ1 , f ðxÞ 2 V 0 () f ðx þ 1Þ 2 V 0 , ¨j2ZVj is dense in L2(R); ˙jVj2Z = {0}.
The existence of a scaling function /(x) is required of which the translates generate a basis in each Vj via V j ¼ spanf/ji gi2Z with /ji ¼ 2j=2 /ð2j x iÞ;
j; i 2 Z:
In the classical case this basis is orthonormal, so that h/ji ; /jk iR ¼ dik ; with hf ; giR ¼
Z
ð2:1Þ
1
f ðxÞ gðxÞdx; 1
being the usual inner product. Let the Wj denote a subspace complementing the subspace Vj in Vj+1, i.e., V jþ1 ¼ V j W j : Each element of Vj+1 can be uniquely written as the sum of an element in Vj and an element in Wj, which contains the details needed to pass from an approximation at level j to an approximation at level j + 1.
M. El-Gamel / Applied Mathematics and Computation 181 (2006) 1635–1644
1637
Based on the function /(x) one can find w(x), the so-called mother wavelet, of which the translates and dilates constitute orthonormal bases of the spaces Wj W j ¼ spanfwji gi2Z ; generated by the wavelets wji ¼ 2j=2 wð2j x iÞ;
j; i 2 Z:
Each function f 2 L2(R), can now be expressed as 1 X X X cj0 i /j0 i ðxÞ þ d ji wji ðxÞ; f ðxÞ ¼ j¼j0
i2Z
ð2:2Þ
i2Z
where cji ¼ hf ; /ji iR
d ji ¼ hf ; wji iR :
Of course, in numerical applications the sum in (2.2) are truncated which corresponds to the projection of f into a subspace of VJ L2(R). The decomposition (2.2) is orthogonal, as, by construction, hwji ; wlk iR ¼djl dik ; hwji ; /lk iR ¼0;
j P l:
in addition to (2.1). 2.2. Regularity and local decay of wavelet coefficients It is well known that the local or global regularity of a function is closely related to the decay of its wavelet coefficients. The latter directly determines the error being made when truncating a wavelet sum at some scale. Depending on the type of norm and whether global or local characterization is considered, various relations of the this kind have been developed. See [4,16] for an overview. Let us recall here only briefly as an example the case of an a-Lipschitz function taken from [5]. Suppose f 2 L2(R), then for [a, b] R the function f is a-Lipschitz for any x0 2 [a, b]. i.e., a
jf ðx0 þ hÞ f ðx0 Þj 6 Cjhj ; if and only if there exists a constant A such that jhf ; wji ij 6 A2jað1=2Þ ; for any (j, i) with i/2j 2 [a, b]. 2.3. Connection coefficients The compactly supported wavelets of Daubechies [4] have proven to be a useful tool in the numerical solutions of partial differential equations [1,7–9,11,13,18–20,25] for recent numerical solutions of partial differential equations using Wavelet-Galerkin techniques. In each of the above papers the notion of connection coefficients has played a key role in the approximations of partial derivatives as well as nonlinear terms in the wavelet discretization of the differential equations. Connection coefficients and some of their properties and algorithms for computing them are described in [17]. For solving Eq. (1.1) in the bounded interval [0, 1], we need to calculate the following two connection coefficients: Z x M mk ðxÞ ¼ y m /ðy kÞdy; ð2:3Þ 0 Z x Cnk ðxÞ ¼ /ðnÞ ðy kÞ/ðyÞdy: ð2:4Þ 0
1638
M. El-Gamel / Applied Mathematics and Computation 181 (2006) 1635–1644
The computation of the connection coefficient (2.3) and (2.4) is reported by Chen et al. [17] with the exception that the normalizing conditions are replaced by " # m1 X X n K n Cnk ðmÞ ¼ n!H1 ðmÞ ð1Þ K n Cnk ðL 1Þ : 2L6k6mLþ1
k¼mLþ2
3. Singularly perturbed second-order boundary value problem In this section, we shall study Wavelet-Galerkin scheme for the singularly perturbed boundary value problem u00 þ au0 þ bu ¼ f ðxÞ
ð3:1Þ
for 0 < x < 1;
subject to boundary conditions uð0Þ ¼ 0;
uð1Þ ¼ 0;
ð3:2Þ
where (0 < 1), and a, b are constants and f(x) is a polynomial of degree any order in x, in [0, 1] otherwise, we approximate them by such a polynomial if necessary. Let the solution uJ(x) of the problem be approximated by its Jth level wavelet series on the interval (0, 1), i.e., J 2X 1
uJ ðxÞ ¼
uk /Jk ðxÞ;
k 2 Z;
ð3:3Þ
k¼2L
where /Jk(x) = 2J/2/(2Jx k), J > 0, and uk, k = 2 L, 3 L, . . . , 2J 1 are 2J + L 2 unknown coefficients to be determined. Here, the integer J is used to control the smoothness of the solution. The larger integer J is used, the more accurate a solution can be obtained. However, the number of equations required to solve the unknown coefficients is increased. In (3.3), the parameter L implies that the wavelet associated with the set of L Daubechies filter coefficients is used as the solution bases. Substituting the wavelet series approximation uJ(x) in Eq. (3.3) for u(x) in Eq. (3.1), yields
J 2X 1
J J 2X 1 2X 1 d2 d / / ðxÞ þ a u ðxÞ þ b uk /Jk ðxÞ ¼ f ðxÞ: k Jk Jk 2 dx dx k¼2L k¼2L
uk
k¼2L
ð3:4Þ
To determine the coefficient uk, we take the inner product of both sides of Eq. (3.4) with /Jl, J J J Z 1 Z 1 Z 1 2X 1 2X 1 2X 1 uk /00Jk ðxÞ/Jl ðxÞdx þ a uk /0J ;k ðxÞ/Jl ðxÞdx þ b uk /Jk ðxÞ/Jl ðxÞdx 0
k¼2L
¼
Z
k¼2L
f ðxÞ/Jl ðxÞdx;
0
We assume that f ðxÞ ¼ J 2X 1
0
k¼2L
0
1
cJkl uk
k¼2L
þa
l ¼ 2 L; 3 L; . . . ; 2J 1:
Pm
J 2X 1
l¼0 ai x
bJkl uk
k¼2L
i
ð3:5Þ
, is a polynomial of degree m in x. We write the equation in (3.5) as
þb
J 2X 1
aJkl uk ¼ d Jml ;
l ¼ 2 L; 3 L; . . . ; 2J 1;
ð3:6Þ
k¼2L
where cJkl ¼
Z
1
/00Jk ðxÞ/Jl ðxÞdx; 0
aJkl ¼
Z
0
bJkl ¼
Z
1
/0J ;k ðxÞ/Jl ðxÞdx;
ð3:7Þ
0 1
/Jk ðxÞ/Jl ðxÞdx
ð3:8Þ
M. El-Gamel / Applied Mathematics and Computation 181 (2006) 1635–1644
1639
and d Jml ¼
Z
1
0
m X
ai xi /Jl ðxÞdx ¼
i¼0
m X i¼0
ai 2
ðiþ12ÞJ
M il ð2J Þ:
ð3:9Þ
Eq. (3.6) can be further put into the matrix-vector form A1 U ¼ D;
ð3:10Þ
where A1 ¼ C þ aB þ bA; C ¼ ½cJkl 2L6k;l62J 1 ; A¼
½aJkl 2L6k;l62J 1 ;
B ¼ ½bJkl 2L6k;l62J 1 ; D¼
½d Jml 2L6l62J 1
ð3:11Þ ð3:12Þ
and U ¼ ½u2L ; u3L ; . . . ; u2J 1 s ; where s denotes the matrix transpose. Now we have a linear system of 2J + L 2 equations of the 2J + L 2 unknown coefficients. We can obtain the coefficients of the approximate solution by solving this linear system. The solution U gives the coefficients in the Wavelet-Galerkin approximation uJ(x) of u(x). 4. Singularly perturbed convection-dominated diffusion equation In this section, we will solve the problem given by Eqs. (1.1)–(1.3). Using the finite difference method and setting ut ¼
ui ui1 ; Dt
then, Eq. (1.1) becomes ui ui1 d2 ui dui ¼ 2 þ gðx; ti Þ: Dt dx dx Rewrite this equation as d2 ui dui 1 i 1 i1 u ¼ u gðx; ti Þ: 2 Dt Dt dx dx
ð4:1Þ
Let the solution uJ(x, ti) of the problem be approximated by its J-th level wavelet series on the interval (0, 1), i.e., uJ ðx; ti Þ ¼
J 2X 1
uik /Jk ðxÞ;
k 2 Z;
ð4:2Þ
k¼2L
Substituting the wavelet series approximation uJ(x, ti) for ui in Eq. (4.1), yields 2X J J J 2X 1 2X 1 1 d2 i d i 1 1 i1 i u u / ðxÞ / ðxÞ u / ðxÞ ¼ u gðx; ti Þ: k Jk k Jk k Jk 2 dx Dt k¼2L Dt dx k¼2L k¼2L To determine the coefficient uik , we take the inner product of both sides of Eq. (4.3) with /Jl, 2X J J J 2X 1 2X 1 1 1 00 0 i i uk h/Jk ðxÞ; /Jl ðxÞi uk h/Jk ðxÞ; /Jl ðxÞi uik h/Jk ðxÞ; /Jl ðxÞi Dt k¼2L k¼2L k¼2L i1 u ¼ h1; /Jl ðxÞi hðti Þhf ðxÞ; /Jl ðxÞi; l ¼ 2 L; 3 L; . . . ; 2J 1: Dt
ð4:3Þ
ð4:4Þ
1640
M. El-Gamel / Applied Mathematics and Computation 181 (2006) 1635–1644
Using the notations defined in (3.7)–(3.9), we write (4.4) as 2X J J J i1 2X 1 2X 1 1 1 u J i J i J i ckl uk bkl uk akl uk ¼ d J0l þ hðti Þd Jml : Dt Dt k¼2L k¼2L k¼2L
ð4:5Þ
The above equations can be further put into the matrix-vector form HU ¼ G;
ð4:6Þ
where
1 A H¼CB Dt
and i1 u G¼ Djm¼0 þ hðti ÞD; Dt and A, B, C and D are defined by Eqs. (3.11) and (3.12). Eq. (4.6) is a linear system of (2J + L 2) equations in (2J + L 2) unknown coefficients. By solving this system, we obtain an approximate solution at resolution level J at time t = ti. Repeating these steps until we get desired time. If we substitute uik ¼ uk ðti Þ in Eq. (4.2) we obtain an approximate solution to the problem (1.1)–(1.3) at time t = ti. 5. Numerical examples The four examples included in this section were selected in order to illustrate the performance of the Wavelet-Galerkin method in solving singularly perturbed problems. For the sake of comparison only, we consider the second and third examples discussed by Reddy [21], who used the method of reduction of order to obtain his numerical solution, and the last example discussed by Kang [24], who used the finite-difference streamlinediffusion method to obtain his numerical solution. In the Wavelet-Galerkin solutions, Daubechies 6 wavelets are used because they give sufficiently accurate results and each example was run for a sequence of J = 5, 6, 7, 9. We use absolute error which is defined as absolute error ¼ juexact solution U Wavelet-Galerkin j:
Example 5.1. Consider the boundary value problem
d2 u þ u ¼ x x2 2; dx2
0
and uð0Þ ¼ 0;
uð1Þ ¼ 0:
The exact solution is given by uðxÞ ¼ xð1 xÞ: By taking different values of at level J = 9 the maximum absolute error is displayed in Table 1. Example 5.2 [21]. Consider the boundary value problem
d2 u du ¼ 1 þ 2x; þ dx2 dx
and uð0Þ ¼ 0;
uð1Þ ¼ 1:
0
M. El-Gamel / Applied Mathematics and Computation 181 (2006) 1635–1644
1641
Table 1 Maximum absolute error for Example 5.1
Maximum absolute error
102 104 105 107 109 1011
1.8966E007 1.7915E009 2.0134E010 1.7957E012 3.6970E014 2.9587E014
Table 2 Maximum absolute error for Example 5.2
Method of reduction of order [21]
Wavelet-Galerkin method
103 104
0.171E04 0.136E04
0.452E07 0.105E07
The exact solution is given by uðxÞ ¼ xðx þ 1 2Þ þ
ð2 1Þð1 ex= Þ : 1 e1=
Table 2 exhibits a comparison between the maximum absolute errors obtained by using the Wavelet-Galerkin method with the analogous results of Reddy [21], who used the method of reduction of order for solving this problem. Table 3 exhibits the exact solution, the Wavelet-Galerkin solution and the absolute errors at level J = 9. Example 5.3 [21]. Consider the boundary value problem
d2 u du u ¼ 0; þ dx2 dx
0
and uð0Þ ¼ 1;
uð1Þ ¼ 1:
The exact solution is given by uðxÞ ¼
ðem2 1Þem1 x þ ð1 em1 Þem2 x ; em2 em1
where m1 ¼
1 þ
pffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ 4 ; 2
and
m2 ¼
1
pffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ 4 : 2
Table 3 The exact solution, the Wavelet-Galerkin solutions, and the absolute errors for Example 5.2 at level J = 9 x
Exact solution
Wavelet-Galerkin solution
Absolute error 1.0E07
0.0 1/8 2/8 3/8 4/8 5/8 6/8 7/8 1.0
0.0 0.85920 0.68735 0.48425 0.2499 0.0157 0.31255 0.64065 1.0
0.0 0.8592000 0.6873500 0.4842500 0.2499000 0.01569999 0.31254999 0.64064999 1.0
0.0 0.1 0.0 0.0 0.0 0.1 0.1 0.1 0.0
1642
M. El-Gamel / Applied Mathematics and Computation 181 (2006) 1635–1644
Table 4 exhibits a comparison between the maximum absolute errors obtained by using the Wavelet-Galerkin method with the analogous results of Reddy [21], who used the method of reduction of order for solving this problem. Table 5 exhibits the exact solution, the Wavelet-Galerkin solution and the absolute errors at level J = 9. Example 5.4 [24]. Consider the boundary value problem ou o2 u ou ¼ 2 þ gðx; tÞ; ot ox ox
0 < x < 1; t > 0
and uð0; tÞ ¼ 1 þ
1 et= ; 1 e1=
uð1; tÞ ¼ 1
and uðx; 0Þ ¼ 1; where gðx; tÞ ¼
½t2 t þ ð1 xÞ ð1xÞt= e : ð1 e1= Þ
The exact solution is given by uðx; tÞ ¼ 1 þ
1 eð1xÞt= : 1 e1=
Table 6 exhibits a comparison between the maximum errors obtained by using the Wavelet-Galerkin method with the analogous results of Kang [24], who used the finite-difference streamline-diffusion method at time t = 0.995. All computations associated with the above examples were performed by using MATLAB. Table 4 Maximum absolute error for Example 5.3
Method of reduction of order [21]
Wavelet-Galerkin method
103 104
0.1763E03 0.1920E04
0.1789E07 0.4083E08
Table 5 The exact solution, the Wavelet-Galerkin solutions, and the absolute errors for Example 5.3 at level J = 9 x
Exact solution
Wavelet-Galerkin solution
Absolute error 1.0E08
0.0 1/8 2/8 3/8 4/8 5/8 6/8 7/8 1.0
1.0 0.416898489 0.472401974 0.535294876 0.606560980 0.687315047 0.778820249 0.882507932 1.0
1.0 0.416898485 0.472401971 0.535294874 0.606560979 0.687315046 0.778820249 0.882507931 1.0
0.0 0.4 0.3 0.2 0.1 0.1 0.0 0.1 0.0
Table 6 Maximum absolute error for Example 5.4 at time t = 0.995
Finite-difference method [24]
Wavelet-Galerkin method
102
0.0124
0.543E04
M. El-Gamel / Applied Mathematics and Computation 181 (2006) 1635–1644
1643
6. Conclusions We have presented and illustrated the Wavelet-Galerkin method for singularly perturbed boundary value problems. An obvious advantage of this procedure is that it loads Daubechies’s constants and calculates the scaling functions, the connection coefficients as well as the rest of the constants only once. This leads to a considerable saving of the computational time and improves numerical results through the reduction of round-off errors. We have implemented the present method on four examples by taking different values of . Computational results are presented in tables. Here we have given results for only few values, although the solutions are computed at all points. The approximate solution is compared with the solution of the method of reduction of order. It can be observed from the results that the Wavelet-Galerkin method yields better results, which shows the efficiency of the method. Acknowledgement The author would like to express sincere thanks to Prof. A. Mohsen at Cairo University, for helpful discussions over the manuscript. References [1] A. Avudainayagam, Wavelet-Galerkin method for integro-differential equations, Appl. Numer. Math. 32 (2000) 247–254. [2] A. Avudainayagam, C. Vani, Wavelet-Galerkin solutions of quasilinear hyperbolic conservation equations, Commun. Numer. Meth. Eng. 15 (1999) 589–601. [3] T. Aziz, A. Khan, A spline method for second-order singularly perturbed boundary-value problems, J. Comput. Appl. Math. 147 (2002) 445–452. [4] I. Daubechies, Ten Lectures on Wavelets, Captial City Press, Vermont, 1992. [5] J. Frohich, K. Schneider, An adaptive wavelet-Galerkin algorithm for one-and two-dimensional flam computations, Eur. J. B/Fluids 13 (1994) 439–471. [6] M. El-Gamel, J.R. Cannon, On the solution of second order singularly-perturbed boundary value problem by the Sinc-Galerkin method, Z. Angew. Math. Phys. 56 (2005) 45–58. [7] M. El-Gamel, A. Zayed, A comparison between the Wavelet-Galerkin and the Sinc-Galerkin methods in solving nonhomogeneous heat equations, in: Inverse Problem, Image Analysis, and Medical Imaging, in: Zuhair Nashed, Otmar Scherzer (Eds.), Contemporary Mathematics of the American Mathematical Society, vol. 313, AMS, Providence, 2002, pp. 97–116. [8] M. El-Gamel, Wavelet algorithm for the numerical solution of nonhomogeneous time-dependent problems, Int. J. Differ. Eqn. Appl. 9 (2004) 169–185. [9] R. Glowiniski, W.M. Lawton, M. Ravachol, E. Tenenbaum, Wavelet solutions of linear and nonlinear elliptic, parabolic and hyperbolic problems in one space dimension, Comput. Meth. Appl. Sci. Eng. (1990) 55–120 (chapter 4). [10] S. Gonzalez-Pinto, L. Casasu´s, P. Gonza´lez-Vera, A numerical scheme to approximate the solution of a singularly perturbed nonlinear differential equation, J. Comput. Appl. Math. 35 (1991) 217–225. [11] S. Jaffard, Wavelet methods for fast resolution of elliptic problems, SIAM J. Numer. Anal. 29 (1992) 965–986. [12] L. Jiang, X. Yue, Local exponentially fitted finite element schemes for singularly perturbed convection–diffusion problems, J. Comput. Appl. Math. 132 (2001) 277–293. [13] F. Jin, T.Q. Ye, Instability analysis of prismatic members by wavelet-Galerkin method, Adv. Eng. Softw. 30 (1999) 361–367. [14] M. Kadalbajoo, A. Rao, The alternating group explicit (AGE) method for singularly perturbed boundary value problems, Appl. Math. Comput. 68 (1995) 125–142. [15] J. Kevorkian, J. Cole, Perturbation Methods in Applied Mathematics, Springer-Verlag, 1981. [16] S. Mallat, A Wavelet Tour of Signal Processing, Academic Press, New York, 1999. [17] M.Q. Chen, C.H. Wang, Y.P. Shin, The computation of wavelet-Galerkin approximation on a bounded interval, Int. J. Numer. Meth. Eng. 39 (1996) 2921–2944. [18] M.Q. Chen, C.H. Wang, Y.P. Shin, A Wavelet-Galerkin method for solving population balanced equations, Comput. Chem. Eng. 20 (1996) 131–145. [19] S. Qian, J. Weiss, Wavelets and the numerical solution of boundary value problems, Appl. Math. Lett. 6 (1993) 47–52. [20] S. Qian, J. Weiss, Wavelets and the numerical solution of partial differential equations, J. Comput. Phys. 106 (1993) 155–175. [21] Y. Reddy, P. Chakravathy, Method of reduction of order for solving singularly perturbed two point boundary value problem, Appl. Math. Comput. 136 (2003). [22] H.G. Roos, M. Stynes, L. Tobiska, Numerical Methods for Singularly Perturbed Differential Equations, Springer-Verlag, 1996. [23] A. Singh, A numerical method for singularly perturbed systems of linear two point boundary value problems using partial decoupling, J. Comput. Appl. Math. 40 (1992) 55–62.
1644
M. El-Gamel / Applied Mathematics and Computation 181 (2006) 1635–1644
[24] T. Kang, D. Yu, Some a posteriori error estimates of the finite-difference streamline-diffusion method for convection-dominated diffusion equations, Adv. Comput. Math. 15 (2001) 193–218. [25] Jin-Chao Xu, Wei-Chang Shann, Galerkin-Wavelet methods for two-point boundary value problems, Numer. Math. 63 (1992) 123– 144.