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] J.L. Gracia, F. Lisbona, A uniformly convergent scheme for a system of reaction–diffusion equations,
To appear in J. Comp. Appl. Math. 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

Clavero García, Carmelo

Gracia Lozano, José Luis

LIsbona Cortés, Francisco

Rodrigo Cardiel, Carmen

idUS. Depósito de Investigación Universidad de Sevilla

XX Congreso de Ecuaciones Diferenciales y AplicacionesX Congreso de Matema´tica AplicadaSevilla, 24-28 septiembre 2007(pp. 1–8)Efficient resolution of singularly perturbed coupled systems:Equations of reaction-diffusionC. Clavero1, J.L. Gracia1, F.J. Lisbona1 C. Rodrigo11 Dpto. Matema´tica Aplicada, Universidad de Zaragoza, Pedro Cerbuna 12, 50009 Zaragoza. E-mails:clavero@unizar.es, jlgracia@unizar.es, lisbona@unizar.es carmenr@unizar.es.Palabras clave: Singular perturbation, reaction-diffusion problems, uniform convergence, coupledsystem, Shishkin meshResumenIn this communication we consider a class of singularly perturbed linear system ofreaction-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 meshof Shishkin type was used. We propose a new semi-implicit method which decouplesthe linear system to be solved at each time level and we prove that it is a uniformlyconvergent scheme (with respect to the diffusion parameters) in the discrete maximumnorm. We display some numerical experiments illustrating in practice the theoreticalresults. From these examples we can see both the uniform convergence of the numericalmethod and also its efficiency to approximate the solution of the reaction–diffusionsystem.1. IntroductionWe 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)whereLx,~ε ≡(−ε1 ∂2∂x2−ε2 ∂2∂x2)+A, A =(a11(x, t) a12(x, t)a21(x, t) a22(x, t)).1C. Clavero, J.L. Gracia, F.J. Lisbona, C. RodrigoWe denote ~ε = (ε1, ε2)T , with 0 < ε1 ≤ ε2 ≤ 1, the vectorial singular perturbationparameter. The coupling matrix A satisfiesaij ≤ 0 if i 6= j, (2)and alsoaii > 0, mı´n(x,t)∈Q¯aii ≥ ma´x(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 guaranteethat the maximum principle holds.Also, we assume that enough regularity and compatibility conditions hold for data ofproblem (1) in order that ~u ∈ C4,2(Q¯), i.e., the spatial partial derivatives of the solutionare continuous up to fourth order and the time partial derivatives are continuous up tosecond 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, dependingon the values of the boundary conditions and the values of the singular perturbationparameters, it can appear two overlapping boundary layers of a width O(ε1/2i ), i = 1, 2,in both sides x = 0 and x = 1 of the domain. This behaviour was proven in [5] and morerecently in [3]. Examples of this type of problems appear in some areas, by instance, inthe study of the flow in fractured porous media (Barenblatt system), in the modellingof diffusion process in bones, considered as a multiple porosity medium, and in diffusionprocess in electroanalytic chemistry.To approximate efficiently the solution of the system it is convenient to have robustmethods in order to achieve a prescribed accuracy using discretization parameters inde-pendent of the diffusion coefficients. This type of methods are called uniformly convergentmethods. In recent years fitted mesh methods has been used extensively where a classicaloperator is defined on a layer–adapted mesh. The meshes introduced by Shiskhin havereceived a special attention and they are piecewise uniform condensing the grid pointsin the boundary layers. In [3] the authors consider a Shishkin mesh to approximate thesolution of problem (1) and on this mesh the backward Euler scheme for the time steppingand the classical central difference scheme for the space discretization are used. In [2] theCrank–Nicolson scheme is used instead of the Euler scheme with the purpose of improvingthe order of the resulting finite difference scheme.In order to be more efficient in the numerical resolution, in this communication we willuse a Jacobi additive scheme for the time discretization which decouples the unknownsand it is well suited to parallel implementation. This is a desirable property, specially forthe multidimensional case and also when the system has a larger number of unknowns(for instance in chemical reactions with n reactants). In the section devoted to the nume-rical experiments we display the corresponding to a system with three equations. Also, inthis communication the uniform convergence, with respect to the singular perturbation2Efficient resolution of singularly perturbed coupled systemsparameters, in the discrete maximum norm of this new scheme is showed. Finally, somenumerical experiments that corroborate the order of uniform convergence of this schemeare given. In addition, the computational advantages over the scheme given in [3] areshown.Henceforth, any positive constant is independent of the diffusion parameters ε1, ε2 andthe discretization parameters N and ∆t.2. Numerical schemeTo approximate the solution of (1) we consider the backward Euler and the centraldifference schemes to discretize the time and spatial variables respectively. The numericalsolution is defined on the meshQ¯N = Ω¯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 Ω¯N in the spatial discretization, condensing the mesh pointsin the boundary layers. The structure of the solution was analyzed in [3] and the authorsproved that the solution has two overlapping boundary layers at both sides x = 0 andx = 1 of the domain of a width O(ε1/2i ), i = 1, 2. The asymptotic analysis of the errorleads to choose the transition points in the mesh as followsτε2 = mı´n {1/4,√ε2 lnN}, τε1 = mı´n {τε2/2,√ε1 lnN}. (5)In the subintervals [0, τε1 ], [τε1 , τε2 ], [τε2 , 1−τε2 ], [1−τε2 , 1−τε1 ] and [1−τε1 , 1] we distributeuniformly N/8+ 1, N/8+ 1, N/2+ 1, N/8+ 1 and N/8+ 1 mesh points respectively (seeFigure 1).τε1τε2 1−τ1−τε ε12N0 1N/8N/8 /2 N/8 N/8Figura 1: Shishkin mesh and number of intervals.So the mesh points are given byxj =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,wherehε1 =8τε1N, hε2 =8(τε2 − τε1)N, H =2(1− 2τε2)N.3C. Clavero, J.L. Gracia, F.J. Lisbona, C. RodrigoOn this mesh, we define the following finite difference scheme(I +∆tLNx,~ε)~Un+1j = ~Unj +∆t ~fn+1j , 0 < j < N, n = 0, . . . ,M − 1, (6)whereLNx,~ε ≡( −ε1δ2−ε2δ2)+An+1j , δ2Zj =2hj + hj+1(Zj+1 − Zjhj+1− Zj − Zj−1hj),with hj = xj − xj−1, j = 1, . . . , N , An+1j = (a(xj , tn+1)), ~fn+1j = ~f(xj , tn+1) and ~Unjdenotes the approximation of the value ~u(xj , tn). So, ~U0j = ~0, 0 ≤ j ≤ N and ~Un+10 =~g0(tn+1), ~Un+1N = ~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+1i } the solution of (6). If thecoefficients of matrix A satisfy the positivity conditions (2) and (3), then‖~u(xj , tn+1)− ~Un+1j ‖Q¯N ≤ C(N−2+q ln2N +∆t), 0 < q < 1, (7)where N,∆t and q are such that N−q ≤ C∆t.In scheme (6), a coupled system must be solved at each time level. For a large numberof variables the resolution of the linear systems can be very expensive and for this reasonwe propose a new scheme. We begin written the scheme (6) as follows~U0j = ~0, 0 ≤ j ≤ N,For n = 0, . . . ,M − 1,(I +∆tT1 ∆tD1∆tD2 I +∆tT2)~Un+1 = ∆t ~fn+1 + ~Un,(8)where I is the (N − 1)–identity matrix, the matrices D1 and D2 are diagonalD1 = diag[a12(xi, tn+1)], D2 = diag[a21(xi, tn+1)], 1 ≤ i ≤ N − 1,and T1 and T2 are (N − 1×N − 1) tridiagonal matricesTk =[ −2εkhj(hj + hj+1),2εkhjhj+1+ akk(xj , tn+1),−2εkhj+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 boundaryconditions must be considered. An important characteristic, that it is fundamental in theanalysis of the convergence, is that T1 and T2 (and also I + ∆tT1 and I + ∆tT2) areM–matrices.We propose the following additive scheme to approximate the solution of problem (1)~V 0j = ~0, 0 ≤ j ≤ N,For n = 0, . . . ,M − 1,(I +∆tT1I +∆tT2)~V n+1 = ∆t ~fn+1 +(I −∆tD1−∆tD2 I)~V n.(9)For this scheme we have the following result of convergence.4Efficient resolution of singularly perturbed coupled systemsTheorem 2 Let ~u(x, t) be the solution of (1) and {~V n+1i } the solution of (9). If thecoefficients of matrix A satisfy the positivity conditions (2) and (3), then‖~u(xj , tn+1)− ~V n+1j ‖Q¯N ≤ C(N−2+q ln2N +∆t), 0 < q < 1, (10)where N,∆t and q are such that N−q ≤ C∆t.Proof. Firstly we determine an estimation for ~En+1 = ~Un+1 − ~V n+1, where ~U is thesolution of problem (8) and ~V is the solution of problem (9). We consider the followingdecomposition~En+1 = ~en+1 + ~Fn+1,where ~en+1 = ~Un+1 − ~V n+1aux and ~Fn+1 = ~V n+1aux − ~V n+1, and ~V n+1aux is the solution of thefollowing auxiliary problem(I +∆tT1I +∆tT2)~V n+1aux = ∆t ~fn+1 +(I −∆tD1−∆tD2 I)~Un. (11)From (8) and (11), we have that ~en+1 is solution of(I +∆tT1I +∆tT2)~en+1 =( −∆tD1−∆tD2)(~Un+1 − ~Un).Using that (I +∆tT1) and (I +∆tT2) are M–matrices, we have that‖~en+1‖Ω¯N ≤ C∆t‖~Un+1 − ~Un‖Ω¯N .Taking into account that ‖~ut‖ ≤ C (see [3]) and Theorem 1, it is straightforward to deducethat‖~en+1‖Ω¯N ≤ C∆t(N−2+q ln2N +∆t).On the other hand, ~Fn+1 is solution of the following problem(I +∆tT1I +∆tT2)~Fn+1 =(I −∆tD1−∆tD2 I)~En.Using that ‖(I + ∆tTk)−1‖Ω¯N < 1, ‖∆t(I + ∆tTk)−1Dk‖Ω¯N < 1, k = 1, 2 and (3), wededuce that‖~Fn+1‖Ω¯N ≤ ‖ ~En‖Ω¯N .Hence, a recursive argument proves‖ ~En+1‖Ω¯N ≤n+1∑i=1‖~ei‖Ω¯N ≤ C∆tM(N−2+q ln2N +∆t) = C(N−2+q ln2N +∆t). (12)The result follows from Theorem 1, the bound (12) and the triangular inequality‖~u(xj , tn+1)− ~V n+1j ‖Q¯N ≤ ‖~u(xj , tn+1)− ~Un+1j ‖Q¯N + ‖ ~En+1‖Q¯N ≤ C(N−2+q ln2N +∆t).A similar result of convergence to Theorem 2 can be established for the additive scheme~V 0j = ~0, 0 ≤ j ≤ N,For n = 0, . . . ,M − 1,(I +∆tT1∆tD2 I +∆tT2)~V n+1 = ∆t ~fn+1 +(I −∆tD1I)~V n.(13)5C. Clavero, J.L. Gracia, F.J. Lisbona, C. Rodrigo3. Numerical experimentsWe 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 equationof the system models the flow in the fracture system and the second equation models theflow in the porous matrix structure (see [1]). In these equations ui are the fluid pressures,εi the permeabilities and κ is the coefficient that control the exchange of fluid between thepores and the fractures. In the numerical experiments we take κ = 1.We use a variant of the double mesh principle to estimate the pointwise errors |~Uni −~u(xi, tn)| in the mesh points {(xi, tn)}. We calculate a new approximation { ~ˆUni } on themesh {(xˆi, tˆn)} that contains the mesh points of the original mesh and their midpointsxˆ2i = xi, i = 0, . . . , N, xˆ2i+1 = (xi + xi+1)/2, i = 0, . . . , N − 1,tˆ2n = tn, n = 0, . . . ,M, tˆ2n+1 = (tn + tn+1)/2, n = 0, . . . ,M − 1.At the mesh points of the coarse mesh we calculate the maximum errors and the uniformerrors by~d~ε,N,∆t = ma´x0≤n≤Mma´x0≤i≤N|~Uni − ~ˆU2n2i |, ~dN,∆t = ma´xSd~ε,N,∆t, (15)where the singular perturbation parameters take values on the setS = {(ε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 corres-ponding 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 thetime discretization parameter ∆t = 0,1, 0,1/2, 0,1/22, 0,1/23, 0,1/24. From these tables weobserve that for this problem the scheme (8) gives better errors than the scheme (9) butthe errors are reduced in the same proportion as far as N increases.To show the advantages of the additive scheme we compare the CPU time of bothschemes (8) and (9). It is clear that scheme (9) become more efficient when the numberof 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)6Efficient resolution of singularly perturbed coupled systemsTabla 1: Scheme (8): Uniform errors ~dN,∆t and uniform orders of convergence ~puni forproblem (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]1 0.123E-1 0.638E-2 0.326E-2 0.165E-2 0.827E-3[~puni]1 0.950 0.969 0.985 0.993[~dN,∆t]2 0.101E-1 0.542E-2 0.281E-2 0.143E-2 0.724E-3[~puni]2 0.893 0.945 0.972 0.986Tabla 2: Scheme (9): Uniform errors ~dN,∆t and uniform orders of convergence ~puni forproblem (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]1 0.438E-1 0.224E-1 0.119E-1 0.621E-2 0.313E-2[~puni]1 0.968 0.912 0.939 0.986[~dN,∆t]2 0.456E-1 0.230E-1 0.121E-1 0.627E-2 0.315E-2[~puni]2 0.985 0.929 0.948 0.992and we use similar schemes to (8) and (9) to approximate its solution. Now three overlap-ping boundary layers can appear at each end of the spatial domain and then we considersix transition points τεi and 1− τεi with i = 1, 2, 3, whereτε3 = mı´n {1/4,√ε3 lnN}, τε2 = mı´n {τε3/2,√ε2 lnN}, τε1 = mı´n {τε2/2,√ε1 lnN},(17)and now we distribute uniformly N/12+ 1 mesh in the fine meshes and N/2+ 1 points inthe coarse mesh. In the numerical experiments the singular perturbation parameters takevalues on the setSˆ = {(ε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 samebehaviour than in the previous test problem (14).Tabla 3: Scheme (8): Uniform errors ~dN,∆t and uniform orders of convergence ~puni forproblem (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]1 0.135E-01 0.658E-02 0.338E-02 0.171E-02 0.861E-03[~puni]1 1.036 0.962 0.981 0.990[~dN,∆t]2 0.153E-01 0.686E-02 0.311E-02 0.157E-02 0.790E-03[~puni]2 1.161 1.141 0.984 0.994[~dN,∆t]3 0.182E-01 0.695E-02 0.281E-02 0.143E-02 0.724E-03[~puni]3 1.391 1.304 0.973 0.986Finally, in Table 5 we show the CPU time (in seconds) to complete the previousTables 1–4. The PC used for the timing results is a Pentium IV with 2,6Mhz. From Table7C. Clavero, J.L. Gracia, F.J. Lisbona, C. RodrigoTabla 4: Scheme (9): Uniform errors ~dN,∆t and uniform orders of convergence ~puni forproblem (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]1 0.762E-01 0.413E-01 0.227E-01 0.121E-01 0.614E-02[~puni]1 0.884 0.864 0.913 0.973[~dN,∆t]2 0.762E-01 0.413E-01 0.227E-01 0.121E-01 0.614E-02[~puni]2 0.884 0.864 0.912 0.973[~dN,∆t]3 0.784E-01 0.421E-01 0.229E-01 0.121E-01 0.615E-02[~puni]3 0.898 0.878 0.919 0.9775 we can deduce the computational advantages of the additive scheme proposed in thiscommunication.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ˆ)AgradecimientosThis research was partially supported by the project MEC/FEDER MTM2004-01905and the Diputacio´n General de Arago´n.Referencias[1] G.I. Barenblatt, I.P. Zheltov, I.N. Kochina, Basic concepts in the theory of seepage of homogeneousliquids 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-diffusionsystems: 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 parabolictype. Translations of Mathematical Monographs, 23, American Mathematical Society, USA, 1968.[5] G.I. Shishkin, Mesh approximation of singularly perturbed boundary-value problems for systems ofelliptic and parabolic equations, Comp. Maths. Math. Phys. 35 (1995) 429–446.8

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

https://idus.us.es/xmlui/bitstream/handle/11441/34929/Efficient%20resolution%20of%20singularly%20perturbed%20coupled%20systems.pdf?sequence=1&isAllowed=y

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

Abstract

Similar works

Full text

Available Versions

idUS. Depósito de Investigación Universidad de Sevilla