Lumped finite element method for reaction-diffusion systems on compact surfaces

We propose and analyse a novel surface finite element method that preserves the invariant regions of systems of semilinear parabolic equations on closed compact surfaces in $\mathbb{R}^3$ under discretisation. We also provide a fully-discrete scheme by applying the implicit-explicit (IMEX) Euler method in time. We prove the preservation of the invariant rectangles of the continuous problem under spatial and full discretizations. For scalar equations, these results reduce to the well-known discrete maximum principle. Furthermore, we prove optimal error bounds for the semi- and fully-discrete methods, that is the convergence rates are quadratic in the meshsize and linear in the timestep. Numerical experiments are provided to support the theoretical findings. In particular we provide examples in which, in the absence of lumping, the numerical solution violates the invariant region leading to blow-up due to the nature of the kinetics.


Introduction
Partial differential equations (PDEs) of the form of reaction-diffusion systems (RDSs) have been extensively employed to model many different processes in a wide range of fields such as biology [45,33,47,24], chemistry [8,61], electrochemistry [6,36] and finance [44,4]. In many applications the domain of integration is a stationary or an evolving curved surface, rather than a planar region. For instance, surface RDSs have been applied to the study of biological patterning [3], tumour growth [9], superconductivity [14], metal dealloying [18], biomembrane modeling [21], cell motility [22] and phase separation [57], just to mention a few examples. In this paper we consider RDSs of arbitrarily many equations on a stationary surface of the form: for i = 1, . . . , r, where Γ is a smooth stationary orientable surface of codimension one in R 3 without boundary, ∆ Γ is the Laplace-Beltrami operator on Γ (which is defined as the tangential divergence of the tangential gradient, see [16] for the definitions), d i are strictly positive diffusion coefficients and u 0i are smooth, bounded functions.
A key feature of many RDSs is the existence of invariant regions. A region Σ in the phase space R r is said to be invariant for (1.1) if, whenever the initial condition has values in Σ, the solution of (1.1) stays in Σ as long as it is defined. Knowing that a given model possesses an invariant region is useful in a couple of ways. First, when solving RDSs arising from applications, solutions are usually meaningful as long as they range within a limited set of values (an example is given by mass-action laws [11], in which solutions are required to be componentwise nonnegative). Second, an invariant region provides an a-priori bound on the analytical solution which can be helpful, for instance, when studying the convergence of numerical methods. Sufficient conditions for a region to be invariant for a given RDS were given in [56,12] on planar domains and in [58, p. 335-353] on stationary surfaces. In both cases, for distinct d i 's, the only possible invariant regions for (1.1) are (bounded or unbounded) hyper-rectangles in R r , that is to say regions in the form with m i , M i ∈ R ∪ {−∞, +∞} for all i = 1, . . . , r, whereas if some d i coincide, more general regions are allowed to be invariant. Since we are addressing general diffusion coefficients, we will consider invariant hyper-rectangles (1.2) only. Among the literature RD models having invariant hyper-rectangles we recall the Gierer-Meinhardt [35], Hodgkin-Huxley [25], FitzHugh-Nagumo [50], Oregonator [63], Rosenzweig-Macarthur [55,29], and the spatially extended Lotka-Volterra [1] models. For r = 1 in (1.1), i.e. scalar semilinear equations, we remark that the min-max condition and the maximum principle, given by respectively, correspond to particular invariant regions, given by Σ = [min Γ u 0 , max Γ u 0 ] and Σ = [0, max Γ u 0 ], respectively. The increasing interest from applications in RDSs on manifolds has stimulated the development of numerical methods for such systems. Among the methods for PDEs on stationary surfaces we recall: finite differences [62], the spectral method of lines [9], closest point methods (see [42] and references therein), kernel methods (see [54] and references therein), embedding methods (see [5] and references therein), and surface finite element methods (SFEM) (see [15,16,60] and references therein). In this paper we consider a lumped mass surface finite element method (LSFEM) for the spatial discretization of Eqs. (1.1). We recall that finite elements with mass lumping have already been applied to reaction-diffusion systems on planar domains, see for example [46,26].
To carry out a fully-discrete scheme we will follow an implicit-explicit (IMEX) approach, i.e. by treating diffusion implicitly and reactions explicitly. Among the class of IMEX methods, we will consider the simplest one, the IMEX Euler scheme considered for example in [43,37]. IMEX methods have been widely applied in fluid dynamics, combined with spectral methods on planar domains [7,32], in reaction-diffusion problems, in combination with finite differences in space on planar domains [52], with finite elements on stationary planar domains [20], on evolving planar domains [43], and with the closest point method on stationary surfaces [41]. An error analysis of finite element approximations with IMEX timestepping for semilinear systems on evolving domains is carried out in [37].
When numerically approximating RDSs, it is desirable for the considered numerical method to preserve invariant rectangles of the continuous problem. For the scalar case, that is the maximum principle, on planar domains, works in this direction cover the homogeneous heat equation (see [10] and references therein), reaction-diffusion equations [46,23,20], anisotropic reaction-diffusion [38] and reaction convection-diffusion equations [39]. For reaction-diffusion systems of many equations on planar domains, the problem is addressed in [31]. The aforementioned works consider different spatial methods. Most of them require the disctretisation to be sufficiently refined, in order to preserve invariant rectangles and maximum principles. A notable exception is the lumped finite element method (LFEM) [20,46,10,38,39]. To the best of our knowledge, numerical methods for surface RDSs preserving the invariant rectangles of the continuous problem have not yet been presented. As far as we know, only a time dependent discrete maximum principle for a scalar diffusion problem on evolving surfaces is given in a recent work [34], in which the evolving surface finite element method (ESFEM) is applied. This motivates the present study in which we introduce the LSFEM, which not only preserves the invariant rectangles at the discrete level, but also requires no restriction on the mesh size.
The main contributions of this paper are twofold. First, we prove discrete maximum principles for the LSFEM semi-discretisation and IMEX Euler-LSFEM full discretisation for a class of semilinear parabolic equations, and the preservation of invariant rectangles under discretisation for weakly coupled (i.e. coupled only through the reaction kinetics) semilinear (i.e. in which only the kinetics are nonlinear) RDSs (1.1). Second, we prove optimal error bounds for the semi-discrete and fully discrete schemes. Among the numerical tests, we provide an example of RDS possessing an invariant region, in which the SFEM blows-up, while the LSFEM preserves the region.
The present article is structured as follows. In Section 2 we consider a semilinear scalar parabolic equation on a closed orientable surface in strong and weak formulation, we present its LSFEM space discretization, its Euler IMEX/LSFEM full discretization and prove the preservation of the maximum principle under spatial and full discretization in Theorems 1 and 2, respectively. In Section 3 we consider a general RDS of arbitrarily many equations on closed orientable surfaces, we derive its LSFEM space discretization, its Euler IMEX/LSFEM time discretization and prove the preservation of the invariant rectangles under spatial and full discretizations in Theorems 4 and 5, respectively. In Section 4, optimal error estimates for both the semi-and fully-discrete methods introduced in the previous sections are proven in Theorems 7 and 8, respectively. Numerical experiments are shown in Section 5.

The continuous problem
We start by considering scalar parabolic PDEs in order to illustrate the main ideas behind the approach described in this work and to introduce the analysis in a less technical setting.
Let Γ be a compact, orientable, smooth surface of codimension one in R 3 without boundary. We assume that Γ is represented as the zero level set of a sufficiently smooth signed distance function d, defined in an open neighbourhood W of Γ such that ∇d( The normal unit vector on Γ is then defined by We assume that every point x ∈ W may be uniquely represented as with a(x) ∈ Γ. A sufficient condition on the thickness of W such that this property holds is given in [16]. We briefly recall the definitions of Sobolev and Bochner spaces on surfaces. For q ∈ N ∪ {0}, the Sobolev space H q (Γ) is the space of functions u : Γ → R such that, for all i = 0, . . . , q, the i-th order tangential derivatives, meant in a distributional sense, are L 2 (Γ), . For further details on Sobolev and Bochner spaces on surfaces we refer the interested reader to [28], [30] or [58].
In this section we consider the following semilinear parabolic equation posed on Γ: where the dot denotes the time derivative, d > 0, α ≥ 1, β ≥ 0, endowed with the nonnegative The requirement that α ≥ 1 is needed to make the source term u α be Lipschitz in a neighbourhood of u = 0, which is a necessary condition for the existence and uniqueness of a solution at all positive times. The conditions β ≥ 0 and u 0 ≥ 0 together are needed to guarantee the maximum principle (1.3). The homogeneous heat equation is obtained as a special case for β = 0. The weak formulation of the problem seeks to find a u ∈ L 2 ([0, T ];

Space discretization
As mentioned previously, in the present work our focus is on finite element discretisations. We now present the necessary notation and concepts needed to describe the numerical method. Given h > 0, a triangulated surface Γ h ⊂ W is defined by where K h is a set of finitely many non degenerate triangles, whose diameters do not exceed h and whose vertices {x i } N i=1 lie on Γ, such that, for a(x) as defined in (2.1), a |Γ h (x) is a one-to-one map between Γ and Γ h ⊂ W .
Following [16], we define lifts and unlifts. Given a function V : Next, let S h be the space of piecewise linear functions on Γ h defined by and S ℓ h be its lifted counterpart: be the nodal basis of S h defined by For v ∈ C 0 (Γ h ), the piecewise linear interpolant I h (v) of v is the function in S h given by We define the following semi-discrete problem: We express the semi-discrete solution U in terms of the nodal basis (2.4) as We then define the lumped mass matrixM = (m ij ) and the stiffness matrix A = (a ij ), respectively, bym We recall that the mass matrix used in the standard SFEM [15,16] is defined by Hence, the semi-discrete problem (2.6) can be expressed as the following ODE system: where ξ := (ξ 1 , . . . , ξ N ) T . In the following we will show that, under suitable assumptions on the triangulation K h , this method fulfills a discrete maximum principle, that is the discrete version of (1.3). To this end we introduce a regularity assumption for the triangulation on the mesh K h which mimicks the standard Delaunay condition on planar domains and then we show how it affects the properties of the stiffness matrix A in (2.9). Let e be an edge of the triangulation K h and let K 1 and K 2 be the triangles sharing the edge e. Let α 1 and α 2 be the angles in K 1 and K 2 opposite to e, respectively. For every edge e in K h we require that α 1 + α 2 ≤ π. (2.11) This condition is represented in Fig. 1. The following result extends to triangulated surfaces a characterization of (2.11) given in [59] for the planar case. Figure 1: Schematic representation of condition (2.11) for triangles K 1 and K 2 .
Proof. Let x i and x j be two distinct nodes of K h . If x i and x j are not neighbours, then Otherwise, let e be the edge connecting x i and x j . Since the intersection of the support of the pyramidal functions χ i and χ j is K 1 ∪ K 2 (see Fig. 1) then we can write Let T 1 and T 2 be two direct isometries (that is with det(J 1 ) = det(J 2 ) = 1) that map K 1 and K 2 into two triangles K 0 1 and K 0 2 contained in the xy plane, respectively, and let J 1 and J 2 be the Jacobians of T 1 and T 2 , respectively. Then, expression (2.13) can be written equivalently as (2.14) Since det(J 1 ) = det(J 2 ) = 1 and ∇ K 0 1 and ∇ K 0 2 both collapse to the standard gradient ∇ in R 2 , expression (2.14) then becomes It is known that (see [59]) this expression only depends on the transformed angles α 0 1 = α 1 , α 0 2 = α 2 and is given by which is nonpositive if and only if α 1 + α 2 ≤ π. This completes the proof. Now, let 1 and 0 be the vector of ones and the null vector in R N , respectively. As shown in [59] (pages 272-273), the structure (2.12) of the stiffness matrix, together with the diagonal structure (2.8) of the lumped mass matrix, imply that, for every s > 0,M + sA is an M-matrix. It then follows that

Time discretization
By applying the IMEX Euler scheme (i.e. treating diffusion implicitly and the reactions explicitly), with time step τ > 0, to (2.10) we obtain the fully-discrete schemē , or equivalently, We remark that, for β = 0 (for the case of the homogeneous heat equation), the timestepping scheme collapses to the standard implicit Euler method.

Semi-and fully-discrete maximum principles
It is known that the lumped FEM fulfills a discrete maximum principle for the homogeneous heat equation on planar domains, see [59]. This result has been generalized to general diffusion problems in divergence form in [46]. The purpose of this section is to extend this result to equation (2.2), which includes, as a special case, the homogeneous heat equation on Γ.
Theorem 1 (Maximum principle for (2.6)). The semi-discrete solution ξ(t) of (2.6) fulfills the following maximum principle Proof. We rewrite (2.10) asξ Consider the auxiliary equationξ where |ξ| and sign(ξ) are the componentwise absolute value and the componentwise sign function of ξ, respectively. If µ = max x∈Γ h v h (x), it is sufficient to prove that the solution of the ODE system (2.21) does not escape the region Σ = [0, µ] N , i.e. we have to prove that, for every ε > 0, the solution of (2.22) does not leave the regionΣ := [−ε, µ] N . To this end, we have to prove that the vector field on the right-hand-side of (2.22), computed on every (N − 1)-dimensional face ofΣ, points toward the interior ofΣ. To this end, let ξ be a point on ∂Σ. This means that there exists i = 1, . . . , N such that ξ i ∈ {−ε, µ}. Suppose ξ i = µ; in the case ξ i = −ε the proof is analogous. Then All we have to prove is thatξ i is negative. To this end, we prove that: 2. The i th component of the vector −dM −1 Aξ is nonpositive. In fact, sinceM is a diagonal matrix, this component is given by We can split the sum on the right-hand-side by isolating the a ii ξ i term: Since a ij ≤ 0 for i = j from Lemma 1 and ξ j ≤ ξ i for j = i from (2.23), expression (2.25) is less than or equal From the definition of A, the right hand side of (2.26) is equal to Since Γ h has no boundary, N j=1 χ i ≡ 1 and thus By combining (2.24)-(2.28), we finally have The above points 1. and 2. imply the desired result thatξ i is negative. This completes the proof.
Proof. From the matrix properties (2.16) and (2.17) we have that, for every τ > 0, In order for the scheme (2.19) to fulfill the maximum principle (2.30), it remains to determine a condition on τ such that i.e. the spatial maximum of the fully-discrete solution is not increasing in time. The inequality (2.34) may be rewritten elementwise as which completes the proof.

Reaction-diffusion systems on surfaces
In this section we consider a more general class of surface PDEs that are reaction-diffusion systems of arbitrarily many equations. Analogously to the semilinear parabolic equation (2.2), we apply a lumped finite element space discretization and an IMEX Euler time discretization. We prove that the LSFEM preserves the invariant hyper-rectangles for the semi-discrete and the fully-discrete problems. For the latter case a time step restriction is required.

The continuous problem
If Γ is a compact orientable surface in R 3 without boundary, as in the previous section, and r ∈ N, let us consider the following reaction-diffusion system of r equations on Γ: where f 1 , . . . , f r are C 2 (Γ r ; R) reaction kinetics and a C 2 (Γ) initial condition is given. As remarked in the Introduction, the following arguments still hold for systems on surfaces with boundary and homogeneous Neumann boundary conditions, i.e. zero conormal derivative on ∂Γ [16]. Then, as a special case, planar bounded domains in R 2 with zero-flux boundary conditions could be included in our study. We will confine the present analysis to the case of compact surfaces without boundary to simplify the presentation. In vector form, system (3.1) is given by In order to write the corresponding vector formulation we extend all the spatial norms considered throughout the paper to vector-valued functions v : Γ → R r or V : Γ h → R r as follows. Given a function space S, we consider the tensor product norm on S r defined by , respectively. Following [2], we introduce the following vector notation: We can now write the sum of the equations (3.3) as

Space discretisation
Analogous to the spatially discretized semilinear parabolic equation (2.6), we define the following space discretization for the reaction-diffusion system (3. By expressing each component u i according to (2.7), we have the following matrix form whereM and A are the lumped mass matrix and the stiffness matrix defined in (2.8) and (2.9), respectively.

Time discretization
By applying the IMEX Euler method to (3.6) we obtain the following fully-discrete method for all n ∈ N ∪ {0} and φ n 1 , . . . , φ n r ∈ S h . We can write the sum of equations (3.8) as for all n ∈ N ∪ {0} and φ n ∈ (S h ) r . System (3.8) can be written in matrix form as that can be obtained equivalently by applying the IMEX Euler method directly to the ODE system (3.7).

Preservation of the invariant rectangles
In this section we investigate an interesting property of the lumped finite element discretization of reaction-diffusion systems, which does not hold in the absence of lumping: the preservation of invariant hyper-rectangles. A numerical counterexample will be given in Section 5. This preservation property is crucial when the continuous system is known to have an invariant rectangle for two reasons: (i) the solution might be physically meaningless outside a certain range of feasible values, containing the rectangle and (ii) it is a tool to prove stability estimates and error bounds for the semi-and fully-discrete solutions. We recall the following definition given in [56,58].
Definition 1. For the system (3.1), a region Σ in the phase-space R r is said to be a positively invariant region if, whenever the initial condition u 0 is in Σ, u stays in Σ as long as it exists and is unique.
The following theorem has been proven in [56] when Γ is a monodimensional domain in R, in [12] when Γ is a k-dimensional domain in R k , k ∈ N (zero-flux boundary conditions are enforced if the domain is not the whole space) and in [58] in the case in which Γ is a Riemannian manifold without boundary. This result provides a sufficient condition for Σ to be a positively invariant region in the phase space.
Theorem 3 (Invariant rectangles for the continuous system (3.1) [58]). Let Σ = r k=1 [m k , M k ] be a hyper-rectangle in the phase space of (3.1), let f be Lipschitz on Σ and let n be the unit outward vector defined piecewise on ∂Σ. If
Further assumptions on f such that the global existence and uniqueness of the solution are ensured can be found in [58]. We remark that some systems are known to possess an invariant region which do not meet the strict inequality (3.11). For instance, for many mass-action laws, the positive orthant is invariant [11] even though the flow of f is tangent to this region, instead of strictly inward.
In the following theorems we prove that, under the same assumptions, Σ is an invariant region for the semi-discrete (3.6) and for the fully-discrete solution (3.10) conditionally on τ , as well. Furthermore, in the fully-discrete case, we will relax the strict inequality (3.11) by requiring non-outward flows, only.
Theorem 4 (Invariant rectangles for (3.6)). Let Σ = r k=1 [m k , M k ] be a hyper-rectangle in the phase space, let f be Lipschitz on Σ and let n be the outward unit normal defined piecewise on ∂Σ. If f (u) · n(u) < 0, ∀u ∈ ∂Σ, (3.12) then Σ is an invariant region for the semi-discrete problem (3.6).
Proof. We rewrite the semi-discrete problem (3.7) as (3.13) It suffices to prove that the rN -dimensional rectangleΣ = r k=1 [m k , M k ] N is an invariant region for the ODE system (3.13), i.e. we have to prove that the vector field . . . All we have to prove is thatξ k,i is negative. To see this, consider that 2. the i th component of the vector −d 1M −1 Aξ k is nonpositive. In fact, sinceM is a diagonal matrix, this component is given by We can split the sum on the right-hand-side by isolating the a ii ξ k,i term: Since a ij ≤ 0 for i = j from Lemma 1 and ξ k,j ≤ ξ k,i for j = i from (3.14), expression (3.16) is less than or equal to From the definition of A (2.9), expression (3.17) can be rewritten as Since Γ h has no boundary, N j=1 χ i ≡ 1 and thus These two claims imply the desired fact, i.e. thatξ k,i is negative. This completes the proof.
The following theorem is a fully-discrete counterpart of the previous one. Observe that the strictly inward flux condition (3.12) is replaced by a weaker requirement. This makes the fully-discrete scheme (3.10) somehow more stable than the spatially discrete one (3.6). The reason for this is that, given a trajectory u(x, t) whose time derivative vanishes at (x,t), the function t → u(x, t) might still be strictly monotonic, this means that a trajectory may escape a region Σ even though the flux of the kinetic is tangent to ∂Σ. In order for the fully-discrete scheme (3.10) to fulfill the theorem, it remains to ensure that m k ≤ ξ n k,i + τ f k (ξ n 1,i , . . . , ξ n r,i ) ≤ M k , ∀ k = 1, . . . , r, ∀ i = 1, . . . , N, ∀ n ∈ N ∪ {0}. (3.23) for all k = 1, . . . , r and n ∈ N. If f k (ξ n 1,i , . . . , ξ n k,i ) > 0, then

Stability and error analysis
In this section we will prove stability estimates and optimal L ∞ ([0, T ], L 2 (Γ)) error bounds for the semi-discrete (3.6) and the fully-discrete (3.10) solutions of the reaction-diffusion system To this end, let us introduce some preliminaries and some basic notations. The lumped L 2 product (see for instance [59,46,48,26]) defined by where I h is given in (2.5), induces the following norm on S h which is equivalento to · L 2 (Γ h ) , uniformly with respect to h (see [51] for the proof): Let us define the "broken" Sobolev space endowed with the norm defined by For the error in the lumped quadrature rule (4.1), if U ∈ H 2 h (Γ h ) and V ∈ S h , then the following estimate holds (see [46]): Inequalities (4.2) and (4.3) have been proven on planar triangulations in [48] and [46], respectively. However, since the respective proofs are done piecewise on each triangle, they can be trivially extended to triangulated surfaces with an affine map argument.
The following equivalences between the norms of a function U defined on Γ h and its lifted counterpart U ℓ can be found in [16].
Lemma 2. Let, K ∈ K h ,K := a(T ) ⊂ Γ, where the map a(x) is given in (2.1), and U : K → R. If the norms exist, then the following inequalities hold where ∇ 2 K and ∇ 2K are the tangential Hessian on K andK, respectively.
From the previous Lemma we derive the following estimate for the broken H 2 norm of U .
Proof. Let K ∈ K h . Then, from (4.4)-(4.6), we have (4.8) Now, from (4.8), we have (4.9) We remark that, in the last inequality of (4.9), the exact equality might not hold, since, being u −ℓ only H 2 h (Γ h ), its gradient ∇ Γ h u −ℓ might have finite jumps across the edges of the triangulation K h . This completes the proof.
When lifting integrals, a geometric error must be taken into account. The following equalities hold (see [16, p.317 where δ ℓ h : Γ → R and R T h : Γ → R 3,3 are functions such that (see [16, p.310 For the following proofs we need to define the seminorm | · | D on (H 1 (Γ)) r and (H 1 (Γ h )) r by 14) respectively. Since the diffusion matrix D is diagonal with positive entries, it holds that i.e. the norms (4.14) and (4.15) are equivalent to | · | H 1 (Γ) and | · | H 1 (Γ h ) , respectively.
The following stability estimates are carried out with the usual energy argument. However, thanks to the existence of an invariant region, the estimates will not depend exponentially on time, as the proofs will not rely on Grönwall's lemma. Moreover, we require that the reaction kinetics f in (3.2) are Lipschitz only in the invariant region, instead of being globally Lipschitz.
for all T > 0, where C is a constant independent of T and u 0 .
Proof. By setting ϕ = u in (3.5) we have Combining (4.16) and (4.20) we have Since u ∈ Σ at all times and f is bounded on Σ, we obtain To prove the second estimate, we set ϕ =u in (3.5) and obtain: Combining (4.22) and (4.23) we have Combining (4.24) with (4.16), we have from which we obtain estimate (4.19).
The following lemmas show analogous estimates for the semi-and fully-discrete problems.
Lemma 5 (Stability estimates for the semi-discrete system (3.6)). If U is the solution of for all T > 0, where C is a constant independent of T and U 0 .
Proof. We proceed exactly as in the previous Lemma in order to obtain analogous estimates in the norm · h and then we use the equivalence (4.2) between the norms · h and · L 2 (Γ h ) on S h , uniformly in h.

28)
for all n = 1, . . . , T τ and T > 0, where C is a constant independent of T and U 0 .
Proof. By testing (3.9) with φ i = U i+1 we have After multiplying by τ , Cauchy-Schwarz inequality yields Since U i and U i+1 ∈ Σ and f is Lipschitz on Σ, the last term on the right hand side is bounded by some constant C > 0: Young's inequality yields We sum for i = 0, . . . , n to obtain By using (4.2), the equivalence between | · | D and | · | H 1 (Γ h ) and n = 1, . . . , T τ , (4.27) follows immediately.

Cauchy-Schwarz inequality yields
Since f is Lipschitz -and thus bounded-on Σ, say max Σ f = C, we can bound the last term in the right hand side as follows: Young's inequality yields Rearranging terms and multiplying by 2 we have By summing (4.29) for i = 0, . . . , n we have By using (4.2), the equivalence between | · | D and | · | H 1 (Γ h ) and n = 1, . . . , T τ , (4.28) finally follows.
To prove the convergence of the semi-and fully-discrete methods, we will adopt the surface Ritz projection considered in [13,19,40].

Definition 2. Given u : [0, T ] → H 1 (Γ), the Ritz projection of u is the unique function
(4.30) We remark that this definition is different from the one considered in [17]. The following error estimates for the Ritz projection can be found in [13,19].
If u is a vector function, we will denote withŪ its componentwise Ritz projection and the estimates (4.31)-(4.32) still hold in the tensor product norms (3.4). An L ∞ ([0, T ], L 2 (Γ)) error bound for the semi-discrete solution has been proven in [46] on planar domains. Here we extend this result to triangulated surfaces. where C(u, T ) is a constant depending on u and T .
Proof. Let us write the error as It remains to show the convergence for θ ℓ in (4.34). For the sake of simplicity, we derive an estimate for θ in the norm · h and then we will use (4.2) and (4.4) to estimate θ ℓ L 2 (Γ) . The continuous problem (3.3), the semi-discrete formulation (3.6), the definition of Ritz projection (4.30) and the relations (4.10), (4.11), imply that (4.37) In (4.37) we choose φ = θ. For the first term of (4.37) we observe that We estimate the single terms on the right hand side of (4.37) in turn. By using the Cauchy-Schwarz inequality, the Lipschitz continuity of f , the definition of θ, (4.2), (4.4) and (4.35), we have that By using the estimate (4.3) for ε h , (4.7), the regularity assumptions f ∈ C 2 (Σ) and u ∈ L ∞ ([0, T ], H 2 (Γ)), and by applying the chain rule to the composite function f (u) it follows that
The previous theorems imply that our semi-and fully-discrete methods exhibit optimal convergence rates, that is to say quadratic in the mesh size and linear in the time step.

Numerical tests
In this section we solve some test problems numerically to show that the LSFEM combined with the IMEX Euler in time: • exhibits the optimal convergence rate predicted in Theorem 8 (Experiments 5.1, 5.4); • fulfills the discrete maximum principle for the homogeneous heat equation, whilst the standard SFEM does not (Experiment 5.2); • preserves the invariant rectangles of reaction-diffusion systems, whilst the standard SFEM does not (Experiment 5.3).
The simulations have been carried out using MATLAB. The linear system arising at each timestep is solved with MATLAB's "backslash" command. The code is available on request.

Experiment 1: The linear heat equation and its convergence
In this experiment we solve the parabolic equation (2.2) in the linear case α = 1 on the unit sphere Γ = {(x, y, z) ∈ R 3 |x 2 + y 2 + z 2 = 1}: with d = 1 24 and β = 1 2 , to test the convergence rate of the LSFEM method. The exact solution of (5.1) is In this experiment, as well as in the following ones, the problem is solved on a sequence of eight meshes Γ i , i = 0, . . . , 7 with decreasing meshsizes h i ≈ √ 2 −i h 0 and corresponding time steps τ i = 2 −i τ 0 (see parameter values in Tab. 1), so that τ i is approximately proportional to h 2 i in order to reveal the quadratic convergence, with respect to the mesh size, of the method. All of the τ i fulfill the stability condition given in Theorem 2. For every i = 0, . . . , 7 the L ∞ ([0, T ], L 2 (Γ h )) error between the numerical solution U and the interpolant I h (u) of the exact solution is measured and the numerical results are reported in Table 1. The lumped solution at the final time T = 1 obtained on the finest mesh is shown in Figure 2 (left), as well as its planar projection through spherical coordinates x = cos φ cos ψ, y = cos φ sin ψ, z = sin ψ, In this test example we observe that the lumped SFEM is more accurate than the standard SFEM and the predicted second order convergence in space is attained.

Experiment 2: The homogeneous heat equation and the maximum principle
We solve the parabolic equation (2.2) for the homogeneous case β = 0 on the unit sphere Γ with d = 0.1 and the nonnegative compactly supported H 1 (Γ) initial datum 2) The minima of the computed numerical solution, obtained for every choice of (h, τ ), are reported in Table 2. In Figure 3 we show the LSFEM solution obtained on the finest mesh

Experiment 3: Reaction-diffusion system and the preservation of the invariant rectangle
In this experiment we consider the reaction-diffusion system with Rosenzweig-MacArthur kinetics [29,26] on the unit sphere Γ, where α, a, b, c, d are positive constants. This system has been numerically solved in [26] on a planar domain with LFEM in combination with an implicit Euler time discretization. However, since the theory developed in [26] addresses a problem on domains of more general dimension (n ≤ 3) there is no discrete maximum principle and the authors consider modified kinetics to ensure the positiveness of the numerical solution. The present example shows that, on two dimensional manifolds, lumping guarantees the preservation of the invariant region without needing modified kinetics.  When c = d and 0 < α < 1 √ 2 for every 0 < ε < 1 − 2 √ a, the rectangle is an invariant region for (5.3), see for instance the analysis in [29]. An easy way to see this is to observe that, for every ε, ε ′ > 0, the rectangle fulfills condition (3.11). Then, since the intersection of invariant regions is still invariant, also Σ is invariant for (5.3). The H 1 (Γ) initial datum with 0 < r < 1, is contained in the invariant region Σ. Furthermore, for 0 < α < 1, it is easy to verify that, on Σ, the Lipschitz constants L 1 and L 2 of the kinetics in (5.3) fulfill In the following we choose d 1 = d 2 = 1e-2, α = 1e-3, a = 10, b = 1e-2, c = d = 1, r = 0.2, and ε = 1e-7. With these settings the invariant region (5.4) becomes

Experiment 4: Reaction-diffusion system with activator-depleted kinetics and its convergence
In this example, we test the convergence rate of the method on a reaction-diffusion system on the unit sphere Γ with well-studied activator-depleted substrate kinetics [27,49,53,45] with an additional forcing term on the right hand side:   the forcing terms f i , i = 1, 2 are approximated with the standard and the lumped quadrature rule given by respectively. We observe that the standard quadrature rule is exact for piecewise linear functions, whilst the lumped one is only exact when the product of the functions is piecewise linear. For this reason, the LSFEM is expected to produce larger errors than the SFEM. The L 2 errors and experimental convergence rates are plotted in Fig. 5 together with the LSFEM solution obtained on the finest mesh at the final time T = 1. As expected, the LSFEM exhibits slightly larger errors than the SFEM. Nonetheless, they have the same convergence rate, in agreement with our theoretical findings.

Conclusions
In Section 2 we have considered a lumped surface finite element method (LSFEM), for a class of semilinear parabolic problems on surfaces, by extending its planar counterpart [46] inspired by the ideas in [16]. Time discretisation is carried out by applying the IMEX Euler method in time. We have shown in Theorem 1 that the spatially discrete problem fulfills a discrete maximum principle. In particular, we have proved that: no restriction on the timestep is required in the homogeneous case (thus extending the result of [59] to surfaces); the time step restriction (2.31) is required in the presence of the nonlinear reaction terms in (2.2). In Section 3 we have applied the LSFEM to general systems of arbitrarily many reaction-diffusion equations. In analogy to the continuous setting (see [12]), in Theorem 4 we have shown that, under the sole assumption of Delaunay regularity for the mesh, the strictly inward flux condition (3.11) is sufficient for a rectangle in the phase space to be invariant for the spatially discrete scheme. For the fully-discrete problem arising from IMEX Euler we have shown in Theorem 5 that, under the time step restriction (3.22) involving the Lipschitz constants of the reaction kinetics, condition (3.11) is not only still sufficient to ensure a hyper-rectangle to be invariant, but can be even weakened by requiring non-outward fluxes (3.21). To the best of the authors' knowledge, Theorems 4 and 5 are a novelty even on planar domains. For both the semi-and fully-discrete formulations of the reaction-diffusion systems considered in Section 3, including the parabolic problem of Section 2 as a special case, an optimal L 2 (Γ) error bound has been proven in Section 4. The numerical examples in Section 5 confirm our theoretical findings. The usefulness of LSFEM is illustrated in Experiments 5.2 and 5.3. In particular, we have shown that in the absence of lumping the numerical solutions of the homogeneous heat equation violates the maximum principle (Exp 5.2) and the numerical solution of a classical predator-prey model blows-up instead of being bounded by the invariant rectangle. Emerging applications encourage the extension of the present study to the case of evolving surfaces, which is beyond the scope of this work and will be addressed in future studies.