Efficient resolution of singularly perturbed coupled systems: Equations of reaction-diffusion

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) Efficient resolution of singularly perturbed coupled systems: Equations of reaction-diffusion C. Clavero1, J.L. Gracia1, F.J. Lisbona1C. Rodrigo1 1Dpto. Matem´atica Aplicada, Universidad de Zaragoza, Pedro Cerbuna 12, 50009 Zaragoza. E-mails: [email protected], [email protected], [email protected] [email protected]. Palabras clave: Singular perturbation, reaction-diffusion problems, uniform convergence, coupled system, Shishkin mesh Resumen In this communication we consider a class of singularly perturbed linear system of reaction-diffusion type coupled in the reaction terms. To approximate its solution, in [3] the backward Euler method and the central difference scheme on a layer–adapted mesh of Shishkin type was used. We propose a new semi-implicit method which decouples the linear system to be solved at each time level and we prove that it is a uniformly convergent scheme (with respect to the diffusion parameters) in the discrete maximum norm. We display some numerical experiments illustrating in practice the theoretical results. From these examples we can see both the uniform convergence of the numerical method and also its efficiency to approximate the solution of the reaction–diffusion system. 1. Introduction We consider reaction–diffusion singularly perturbed problems      L~ε~u ≡∂~u ∂t +Lx,~ε~u =~ f, (x, t)∈Q= Ω ×(0, T] = (0,1) ×(0, T], ~u(0, t) = ~g0(t), ~u(1, t) = ~g1(t),∀t∈[0, T], ~u(x, 0) = ~ 0,∀x∈Ω, (1) where Lx,~ε ≡Ã−ε1∂2 ∂x2 −ε2∂2 ∂x2!+A, A =µa11(x, t)a12(x, t) a21(x, t)a22(x, t)¶. 1 C. Clavero, J.L. Gracia, F.J. Lisbona, C. Rodrigo We denote ~ε = (ε1, ε2)T, with 0 < ε1≤ε2≤1, the vectorial singular perturbation parameter. The coupling matrix Asatisfies aij ≤0 if i6=j, (2) and also aii >0,m´ın (x,t)∈¯ Qaii ≥m´ax (x,t)∈¯ Q|aij|, i 6=j, i, j = 1,2.(3) If condition (3) does not hold we consider the transformation ~v(x, t) = ~u(x, t)e−α0t, with α0>0 sufficiently large, in order to transform diagonal entries. These hypothesis guarantee that the maximum principle holds. Also, we assume that enough regularity and compatibility conditions hold for data of problem (1) in order that ~u ∈C4,2(¯ Q), i.e., the spatial partial derivatives of the solution are continuous up to fourth order and the time partial derivatives are continuous up to second order. For instance, we will suppose the conditions ~g(k) i(0) = ~ 0, i = 0,1, k = 0,1,2, ∂k+k0~ f ∂xk∂tk0(0,0) = ∂k+k0~ f ∂xk∂tk0(1,0) = ~ 0,0≤k+ 2k0≤2,(4) which are an extension of the compatibility conditions for the scalar case (see [4]). These small parameters cause a multiscale character of the solution and, depending on the values of the boundary conditions and the values of the singular perturbation parameters, it can appear two overlapping boundary layers of a width O(ε1/2 i), i = 1,2, in both sides x= 0 and x= 1 of the domain. This behaviour was proven in [5] and more recently in [3]. Examples of this type of problems appear in some areas, by instance, in the study of the flow in fractured porous media (Barenblatt system), in the modelling of diffusion process in bones, considered as a multiple porosity medium, and in diffusion process in electroanalytic chemistry. To approximate efficiently the solution of the system it is convenient to have robust methods in order to achieve a prescribed accuracy using discretization parameters independent of the diffusion coefficients. This type of methods are called uniformly convergent methods. In recent years fitted mesh methods has been used extensively where a classical operator is defined on a layer–adapted mesh. The meshes introduced by Shiskhin have received a special attention and they are piecewise uniform condensing the grid points in the boundary layers. In [3] the authors consider a Shishkin mesh to approximate the solution of problem (1) and on this mesh the backward Euler scheme for the time stepping and the classical central difference scheme for the space discretization are used. In [2] the Crank–Nicolson scheme is used instead of the Euler scheme with the purpose of improving the order of the resulting finite difference scheme. In order to be more efficient in the numerical resolution, in this communication we will use a Jacobi additive scheme for the time discretization which decouples the unknowns and it is well suited to parallel implementation. This is a desirable property, specially for the multidimensional case and also when the system has a larger number of unknowns (for instance in chemical reactions with nreactants). In the section devoted to the numerical experiments we display the corresponding to a system with three equations. Also, in this communication the uniform convergence, with respect to the singular perturbation 2 Efficient resolution of singularly perturbed coupled systems parameters, in the discrete maximum norm of this new scheme is showed. Finally, some numerical experiments that corroborate the order of uniform convergence of this scheme are given. In addition, the computational advantages over the scheme given in [3] are shown. Henceforth, any positive constant is independent of the diffusion parameters ε1,ε2and the discretization parameters Nand ∆t. 2. Numerical scheme To approximate the solution of (1) we consider the backward Euler and the central difference schemes to discretize the time and spatial variables respectively. The numerical solution is defined on the mesh ¯ QN=¯ ΩNׯωN, where, for simplicity, we consider a uniform mesh in the time discretization ¯ωN={k∆t, 0≤k≤M, ∆t=T/M}, and a piecewise uniform mesh ¯ ΩNin the spatial discretization, condensing the mesh points in the boundary layers. The structure of the solution was analyzed in [3] and the authors proved that the solution has two overlapping boundary layers at both sides x= 0 and x= 1 of the domain of a width O(ε1/2 i), i = 1,2. The asymptotic analysis of the error leads to choose the transition points in the mesh as follows τε2= m´ın {1/4,√ε2ln N}, τε1= m´ın {τε2/2,√ε1ln N}.(5) In the subintervals [0, τε1], [τε1, τε2], [τε2,1−τε2], [1−τε2,1−τε1] and [1−τε1,1] we distribute uniformly N/8 + 1, N/8 + 1, N/2 + 1, N/8 + 1 and N/8 + 1 mesh points respectively (see Figure 1). τ ε 1 τε2 1−τ1−τε ε12 N 0 1 N/8N/8 N/2 N/8 N/8 Figura 1: Shishkin mesh and number of intervals. So the mesh points are given by xj=           jhε1, j = 0, . . . , N/8, xN/8+ (j−N/8)hε2, j =N/8 + 1, . . . , N/4, xN/4+ (j−N/4)H, j =N/4 + 1, . . . , 3N/4, x3N/4+ (j−3N/4)hε2, j = 3N/4+1,...,7N/8, x7N/8+ (j−7N/8)hε1, j = 7N/8+1, . . . , N, where hε1=8τε1 N, hε2=8(τε2−τε1) N, H =2(1 −2τε2) N. 3 C. Clavero, J.L. Gracia, F.J. Lisbona, C. Rodrigo On this mesh, we define the following finite difference scheme (I+ ∆tLN x,~ε)~ Un+1 j=~ Un j+ ∆t~ fn+1 j,0< j < N, n = 0, . . . , M −1,(6) where LN x,~ε ≡µ−ε1δ2 −ε2δ2¶+An+1 j, δ2Zj=2 hj+hj+1 µZj+1 −Zj hj+1 −Zj−Zj−1 hj¶, with hj=xj−xj−1, j = 1, . . . , N,An+1 j= (a(xj, tn+1)), ~ fn+1 j=~ f(xj, tn+1) and ~ Un j denotes the approximation of the value ~u(xj, tn). So, ~ U0 j=~ 0,0≤j≤Nand ~ Un+1 0= ~g0(tn+1),~ Un+1 N=~g1(tn+1). In [3] the following result of convergence was proved. Theorem 1 Let ~u(x, t)be the solution of (1) and {~ Un+1 i}the solution of (6). If the coefficients of matrix Asatisfy the positivity conditions (2) and (3), then k~u(xj, tn+1)−~ Un+1 jk¯ QN≤C(N−2+qln2N+ ∆t),0< q < 1,(7) where N, ∆tand qare such that N−q≤C∆t. In scheme (6), a coupled system must be solved at each time level. For a large number of variables the resolution of the linear systems can be very expensive and for this reason we propose a new scheme. We begin written the scheme (6) as follows ~ U0 j=~ 0,0≤j≤N, For n= 0, . . . , M −1, µI+ ∆tT1∆tD1 ∆tD2I+ ∆tT2¶~ Un+1 = ∆t~ fn+1 +~ Un, (8) where Iis the (N−1)–identity matrix, the matrices D1and D2are diagonal D1=diag[a12(xi, tn+1)], D2=diag[a21(xi, tn+1)],1≤i≤N−1, and T1and T2are (N−1×N−1) tridiagonal matrices Tk=·−2εk hj(hj+hj+1),2εk hjhj+1 +akk(xj, tn+1),−2εk hj+1(hj+hj+1)¸,1< j < N −1, k = 1,2. We note that the first and last rows of these matrices are special because the boundary conditions must be considered. An important characteristic, that it is fundamental in the analysis of the convergence, is that T1and T2(and also I+ ∆tT1and I+ ∆tT2) are M–matrices. We propose the following additive scheme to approximate the solution of problem (1) ~ V0 j=~ 0,0≤j≤N, For n= 0, . . . , M −1, µI+ ∆tT1 I+ ∆tT2¶~ Vn+1 = ∆t~ fn+1 +µI−∆tD1 −∆tD2I¶~ Vn. (9) For this scheme we have the following result of convergence. 4 Efficient resolution of singularly perturbed coupled systems Theorem 2 Let ~u(x, t)be the solution of (1) and {~ Vn+1 i}the solution of (9). If the coefficients of matrix Asatisfy the positivity conditions (2) and (3), then k~u(xj, tn+1)−~ Vn+1 jk¯ QN≤C(N−2+qln2N+ ∆t),0< q < 1,(10) where N, ∆tand qare such that N−q≤C∆t. Proof. Firstly we determine an estimation for ~ En+1 =~ Un+1 −~ Vn+1, where ~ Uis the solution of problem (8) and ~ Vis the solution of problem (9). We consider the following decomposition ~ En+1 =~en+1 +~ Fn+1, where ~en+1 =~ Un+1 −~ Vn+1 aux and ~ Fn+1 =~ Vn+1 aux −~ Vn+1, and ~ Vn+1 aux is the solution of the following auxiliary problem µI+ ∆tT1 I+ ∆tT2¶~ Vn+1 aux = ∆t~ fn+1 +µI−∆tD1 −∆tD2I¶~ Un.(11) From (8) and (11), we have that ~en+1 is solution of µI+ ∆tT1 I+ ∆tT2¶~en+1 =µ−∆tD1 −∆tD2¶(~ Un+1 −~ Un). Using that (I+ ∆tT1) and (I+ ∆tT2) are M–matrices, we have that k~en+1k¯ ΩN≤C∆tk~ Un+1 −~ Unk¯ ΩN. Taking into account that k~utk ≤ C(see [3]) and Theorem 1, it is straightforward to deduce that k~en+1k¯ ΩN≤C∆t(N−2+qln2N+ ∆t). On the other hand, ~ Fn+1 is solution of the following problem µI+ ∆tT1 I+ ∆tT2¶~ Fn+1 =µI−∆tD1 −∆tD2I¶~ En. Using that k(I+ ∆tTk)−1k¯ ΩN<1, k∆t(I+ ∆tTk)−1Dkk¯ ΩN<1, k = 1,2 and (3), we deduce that k~ Fn+1k¯ ΩN≤ k~ Enk¯ ΩN. Hence, a recursive argument proves k~ En+1k¯ ΩN≤ n+1 X i=1 k~eik¯ ΩN≤C∆tM(N−2+qln2N+ ∆t) = C(N−2+qln2N+ ∆t).(12) The result follows from Theorem 1, the bound (12) and the triangular inequality k~u(xj, tn+1)−~ Vn+1 jk¯ QN≤ k~u(xj, tn+1)−~ Un+1 jk¯ QN+k~ En+1k¯ QN≤C(N−2+qln2N+ ∆t). A similar result of convergence to Theorem 2 can be established for the additive scheme ~ V0 j=~ 0,0≤j≤N, For n= 0, . . . , M −1, µI+ ∆tT1 ∆tD2I+ ∆tT2¶~ Vn+1 = ∆t~ fn+1 +µI−∆tD1 I¶~ Vn. (13) 5 C. Clavero, J.L. Gracia, F.J. Lisbona, C. Rodrigo 3. Numerical experiments We consider the following test problem ∂u1 ∂t −ε1 ∂2u1 ∂x2+κ(u1−u2) = 1, ∂u2 ∂t −ε2 ∂2u2 ∂x2+κ(u2−u1) = 1,     (x, t)∈(0,1) ×(0,1], ~u(0, t) = ~u(1, t) = 0, t ∈[0, T], ~u(x, 0) = 0, x ∈[0,1]. (14) This coupled system is used to model the flow in fractured porous media. The first equation of the system models the flow in the fracture system and the second equation models the flow in the porous matrix structure (see [1]). In these equations uiare the fluid pressures, εithe permeabilities and κis the coefficient that control the exchange of fluid between the pores and the fractures. In the numerical experiments we take κ= 1. We use a variant of the double mesh principle to estimate the pointwise errors |~ Un i− ~u(xi, tn)|in the mesh points {(xi, tn)}. We calculate a new approximation {ˆ ~ Un i}on the mesh {(ˆxi,ˆ tn)}that contains the mesh points of the original mesh and their midpoints ˆx2i=xi, i = 0, . . . , N, ˆx2i+1 = (xi+xi+1)/2, i = 0, . . . , N −1, ˆ t2n=tn, n = 0, . . . , M, ˆ t2n+1 = (tn+tn+1)/2, n = 0, . . . , M −1. At the mesh points of the coarse mesh we calculate the maximum errors and the uniform errors by ~ d~ε,N,∆t= m´ax 0≤n≤Mm´ax 0≤i≤N|~ Un i−ˆ ~ U2n 2i|,~ dN,∆t= m´ax Sd~ε,N,∆t,(15) where the singular perturbation parameters take values on the set S={(ε1, ε2)|ε2= 20,2−2, . . . , 2−30, ε1=ε2,2−2ε2, . . . , 2−58,2−60}, to permit that the maximum errors stabilize. From values (15) we determinate the corresponding orders of convergence in a standard way ~p =log(~ d~ε,N,∆t/~ d~ε,2N,∆t/2) log 2 , ~puni =log(~ dN,∆t/~ d2N,∆t/2) log 2 . In Tables 1 and 2 we display the numerical results for the schemes (8) and (9) respectively. The first row corresponds to the first variable and the second one to the second variable. The spatial discretization parameter takes the values N= 64,128,256,512,1024 and the time discretization parameter ∆t= 0,1,0,1/2,0,1/22,0,1/23,0,1/24. From these tables we observe that for this problem the scheme (8) gives better errors than the scheme (9) but the errors are reduced in the same proportion as far as Nincreases. To show the advantages of the additive scheme we compare the CPU time of both schemes (8) and (9). It is clear that scheme (9) become more efficient when the number of equations increases. For this reason, we introduce the following test problem ∂u1 ∂t −ε1 ∂2u1 ∂x2+κ(2u1−u2−u3) = 1, ∂u2 ∂t −ε2 ∂2u2 ∂x2+κ(2u2−u1−u3) = 1, ∂u3 ∂t −ε3 ∂2u3 ∂x2+κ(2u3−u1−u2) = 1,              (x, t)∈(0,1) ×(0,1], ~u(0, t) = ~u(1, t) = 0, t ∈[0, T], ~u(x, 0) = 0, x ∈[0,1], (16) 6 Efficient resolution of singularly perturbed coupled systems Tabla 1: Scheme (8): Uniform errors ~ dN,∆tand uniform orders of convergence ~puni for problem (14). N=64 N=128 N=256 N=512 N=1024 ε1, ε2∈S∆t= 0,1 ∆t= 0,1/2 ∆t= 0,1/22∆t= 0,1/23∆t= 0,1/24 [~ dN,∆t]10.123E-1 0.638E-2 0.326E-2 0.165E-2 0.827E-3 [~puni]10.950 0.969 0.985 0.993 [~ dN,∆t]20.101E-1 0.542E-2 0.281E-2 0.143E-2 0.724E-3 [~puni]20.893 0.945 0.972 0.986 Tabla 2: Scheme (9): Uniform errors ~ dN,∆tand uniform orders of convergence ~puni for problem (14). N=64 N=128 N=256 N=512 N=1024 ε1, ε2∈S∆t= 0,1 ∆t= 0,1/2 ∆t= 0,1/22∆t= 0,1/23∆t= 0,1/24 [~ dN,∆t]10.438E-1 0.224E-1 0.119E-1 0.621E-2 0.313E-2 [~puni]10.968 0.912 0.939 0.986 [~ dN,∆t]20.456E-1 0.230E-1 0.121E-1 0.627E-2 0.315E-2 [~puni]20.985 0.929 0.948 0.992 and we use similar schemes to (8) and (9) to approximate its solution. Now three overlapping boundary layers can appear at each end of the spatial domain and then we consider six transition points τεiand 1 −τεiwith i= 1,2,3,where τε3= m´ın {1/4,√ε3ln N}, τε2= m´ın {τε3/2,√ε2ln N}, τε1= m´ın {τε2/2,√ε1ln N}, (17) and now we distribute uniformly N/12 + 1 mesh in the fine meshes and N/2 + 1 points in the coarse mesh. In the numerical experiments the singular perturbation parameters take values on the set ˆ S={(ε1, ε2, ε3)|ε3= 20,2−2, . . . , 2−30, ε2=ε3,2−2ε2, . . . , 2−38,2−40, ε1=ε2,2−2ε2, . . . , 2−58,2−60}. In Tables 3 and 4 we show the numerical results for problem (16) observing the same behaviour than in the previous test problem (14). Tabla 3: Scheme (8): Uniform errors ~ dN,∆tand uniform orders of convergence ~puni for problem (16). N=64 N=128 N=256 N=512 N=1024 ε1, ε2, ε3∈ˆ S∆t= 0,1 ∆t= 0,1/2 ∆t= 0,1/22∆t= 0,1/23∆t= 0,1/24 [~ dN,∆t]10.135E-01 0.658E-02 0.338E-02 0.171E-02 0.861E-03 [~puni]11.036 0.962 0.981 0.990 [~ dN,∆t]20.153E-01 0.686E-02 0.311E-02 0.157E-02 0.790E-03 [~puni]21.161 1.141 0.984 0.994 [~ dN,∆t]30.182E-01 0.695E-02 0.281E-02 0.143E-02 0.724E-03 [~puni]31.391 1.304 0.973 0.986 Finally, in Table 5 we show the CPU time (in seconds) to complete the previous Tables 1–4. The PC used for the timing results is a Pentium IV with 2,6Mhz. From Table 7 C. Clavero, J.L. Gracia, F.J. Lisbona, C. Rodrigo Tabla 4: Scheme (9): Uniform errors ~ dN,∆tand uniform orders of convergence ~puni for problem (16). N=64 N=128 N=256 N=512 N=1024 ε1, ε2, ε3∈ˆ S∆t= 0,1 ∆t= 0,1/2 ∆t= 0,1/22∆t= 0,1/23∆t= 0,1/24 [~ dN,∆t]10.762E-01 0.413E-01 0.227E-01 0.121E-01 0.614E-02 [~puni]10.884 0.864 0.913 0.973 [~ dN,∆t]20.762E-01 0.413E-01 0.227E-01 0.121E-01 0.614E-02 [~puni]20.884 0.864 0.912 0.973 [~ dN,∆t]30.784E-01 0.421E-01 0.229E-01 0.121E-01 0.615E-02 [~puni]30.898 0.878 0.919 0.977 5 we can deduce the computational advantages of the additive scheme proposed in this communication. Tabla 5: Schemes (8) and (9): CPU times for problems (14) and (16). Scheme (8) Scheme (9) Test problem (14) 348.219” 228.125” (ε1, ε2∈S) Test problem (16) 34749.2340” 4601.6880” (ε1, ε2, ε3∈ˆ S) Agradecimientos This research was partially supported by the project MEC/FEDER MTM2004-01905 and the Diputaci´on General de Arag´on. Referencias [1] G.I. Barenblatt, I.P. Zheltov, I.N. Kochina, Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks, J. Appl. Math. and Mech. 24 (1960) 1286–1303. [2] C. Clavero, J.L. Gracia, F. Lisbona, Time dependent singularly perturbed coupled reaction-diffusion systems: a high accurate uniformly convergent method, Submitted [3] J.L. Gracia, F. Lisbona, A uniformly convergent scheme for a system of reaction–diffusion equations, To appear in J. Comp. Appl. Math. [4] O.A. Ladyzhenskaya, V.A. Solonikov, N. N. Ural’tseva, Linear and quasilinear equations of parabolic type. Translations of Mathematical Monographs, 23, American Mathematical Society, USA, 1968. [5] G.I. Shishkin, Mesh approximation of singularly perturbed boundary-value problems for systems of elliptic and parabolic equations, Comp. Maths. Math. Phys. 35 (1995) 429–446. 8