Numerical Preservation of Velocity Induced Invariant Regions for Reaction–Diffusion Systems on Evolving Surfaces

Numerical Preservation of Velocity Induced Invariant Regions for Reaction–Diffusion Systems on... J Sci Comput https://doi.org/10.1007/s10915-018-0741-7 Numerical Preservation of Velocity Induced Invariant Regions for Reaction–Diffusion Systems on Evolving Surfaces 1 2 Massimo Frittelli · Anotida Madzvamuse · 1 3 Ivonne Sgura · Chandrasekhar Venkataraman Received: 13 December 2017 / Revised: 4 May 2018 / Accepted: 5 May 2018 © The Author(s) 2018 Abstract We propose and analyse a finite element method with mass lumping (LESFEM) for the numerical approximation of reaction–diffusion systems (RDSs) on surfaces in R that evolve under a given velocity field. A fully-discrete method based on the implicit–explicit (IMEX) Euler time-discretisation is formulated and dilation rates which act as indicators of the surface evolution are introduced. Under the assumption that the mesh preserves the Delaunay regularity under evolution, we prove a sufficient condition, that depends on the dilation rates, for the existence of invariant regions (i) at the spatially discrete level with no restriction on the mesh size and (ii) at the fully-discrete level under a timestep restriction that depends on the kinetics, only. In the specific case of the linear heat equation, we prove a semi- and a fully-discrete maximum principle. For the well-known activator-depleted and Thomas reaction–diffusion models we prove the existence of a family of rectangles in the phase space that are invariant only under specific growth laws. Two numerical examples are provided to computationally demonstrate (i) the discrete maximum principle and optimal convergence for the heat equation on a linearly growing sphere and (ii) the existence of an invariant region for the LESFEM–IMEX Euler discretisation of a RDS on a logistically growing surface. Anotida Madzvamuse A.Madzvamuse@sussex.ac.uk Massimo Frittelli massimo.frittelli@unisalento.it Ivonne Sgura ivonne.sgura@unisalento.it Chandrasekhar Venkataraman cv28@st-andrews.ac.uk Dipartimento di Matematica e Fisica “E. De Giorgi”, Università del Salento, via per Arnesano, 73100 Lecce, Italy Department of Mathematics, School of Mathematical and Physical Sciences, University of Sussex, Pevensey III, 5C15, Brighton BN1 9QH, UK School of Mathematics and Statistics, University of St Andrews, Fife KY16 9SS, UK 123 J Sci Comput Keywords Evolving surface · Dilation rate · Heat equation · Maximum principle · Reaction–diffusion · Invariant region 1 Introduction Reaction–diffusion systems (RDSs) are a well-known class of partial differential equations that arise from the mathematical modelling of numerous phenomena taking place on a space- time domain, see for instance [31]. In many applications the spatial domain is a curved surface rather than a flat domain, which may be time-dependent. Among the applications of RDSs on surfaces we mention brain growth [26], cell migration [3], chemotaxis [12], developmental biology [28], electrodeposition [24] and phase field modeling [42]. The growing interest toward PDEs on evolving surfaces has stimulated the development of several numerical methods for such problems, among which we mention (but not limited to) embedding methods [2], kernel methods [18], implicit boundary integral methods [5,35], surface finite element methods (SFEM) [10] and some of their recent variations and extensions [13,16,17,20,23, 40]. An interesting property of some RDSs is the existence of invariant regions.Aninvariant region is a subset  of the phase-space such that, if the initial condition takes values in ,then the solution of the RDS takes values in  at all times. We recall that, for scalar equations, the well-known notion of maximum principle is equivalent to the invariance of all the regions of the form [0, M ], M > 0. For RDSs on a stationary surface, sufficient conditions have been found for a region to be invariant at the continuous level, see [39]. To the best of the authors’ knowledge, the extension of these results to RDSs on evolving surfaces has not been considered in the literature. In this paper we will focus on surfaces that evolve according to a prescribed material velocity field. More complicated forms of the evolution law are beyond the scope of this manuscript and are a subject of our current studies. For a numerical method it is interesting to understand if the invariant regions of the continuous problem are preserved under discretisation. In [15,16]weprovedthat, forRDSs on stationary surfaces, the lumped surface finite element method (LSFEM), combined with an implicit–explicit (IMEX) Euler method, preserves the invariant regions of the continuous problem. The purpose of the present paper is to extend these results to the case of evolving surfaces, in particular (i) we prove a semi- and a fully-discrete maximum principle for the heat equation with a linear source term and (ii) we provide sufficient conditions under which a region is invariant at the semi- and fully-discrete levels when a lumped evolving surface finite element method (LESFEM) and an IMEX Euler timestepping are considered. In particular we quantify the impact of surface evolution (measured through the dilation rate)onthe existence of invariant regions and we find that surface growth or contraction respectively fosters or inhibits the invariance of a given region in the phase-space. Crucial in our analysis is the assumption that the mesh preserves the Delaunay property under evolution, which is not true for an arbitrary surface evolution law. A class of surface evolution laws for which the Delaunay property is automatically preserved is that of isotropic growth [29], which has biological applications [7,27,33,36]. To the best of the authors’ knowledge, an adaptive strategy for the preservation of the Delaunay property under a generic evolution law is still an open problem. An attempt in this direction is the work in [22]. For the special case of isotropic growth, we provide fully practical sufficient conditions for the existence of invariant regions at the semi- and fully-discrete levels. As an application of our general theory, we classify some classes of invariant regions, depending on the growth rate of 123 J Sci Comput the evolving surface, for two well-known RD models in the literature: the activator-depleted (or Schnakenberg, also known as the Brusselator model) and the Thomas models. Finally, we provide two numerical examples. In the first example we experimentally show that the LESFEM–IMEX Euler method, applied to the heat equation with a linear source term on a linearly growing sphere, exhibits optimal convergence rates in space and time. In the second example we consider the Thomas RDS on an exponentially growing Dupin ring cyclide thereby showing (i) the existence of invariant regions for the fully discretised model and (ii) the violation of this region in the absence of mass lumping. The present paper is structured as follows. First, in Sect. 2 we recall (i) the derivation of RDSs on evolving surfaces and (ii) some basic notions about invariant regions and we introduce the notion of dilation rate, which is crucial in our analysis. In Sect. 3 we introduce the LESFEM for the space discretisation of RDSs on evolving surfaces and we carry out a fully- discrete scheme using the IMEX Euler timestepping. Section 4 deals with the characterisation of the dilation rate in terms of the material velocity. In particular, we compute exactly the dilation rate for the class of isotropic growth laws. In Sect. 5 we prove a semi- and a fully- discrete maximum principle for the linear heat equation on evolving surfaces with a linear source term. We prove, in Sect. 6, sufficient conditions for the existence of invariant regions for RDSs of arbitrarily many equations on evolving surfaces at the semi- and fully-discrete levels. Section 7 presents some classes of invariant regions for the activator-depleted and the Thomas RD models on evolving surfaces, respectively. Numerical examples are presented in Sect. 8. Finally, in Sect. 9 we conclude and discuss our findings with an eye for future extensions of the present work. 2 Reaction–Diffusion Equations on an Evolving Surface 2.1 Preliminaries and Basic Results 2 3 For t ∈[0, T ],let Γ(t ) be a C orientable surface in R , represented as the zero-level set 1 2 3 3 of a signed distance function d ∈ C ([0, T ], C (R )), i.e. Γ(t ) ={x ∈ R | d(x, t ) = 0}, with ∇d(x, t ) = 0for t ∈[0, T ] and x ∈ Γ(t ). Hence, the outward unit normal vector field on Γ(t ) is given by ∇d(x, t ) n(x, t ) := , t ∈[0, T ], x ∈ Γ(t ), (1) ∇d(x, t ) where · denotes the Euclidean norm in R . Following [10, Section 5], we assume that 3 1 2 there exists a mapping G : Γ(0) ×[0, T]→ R , G ∈ C ([0, T ], C (Γ (0))), such that for all t ∈[0, T ], G(Γ (0), t ) = Γ(t ) and G(·, t ) is a diffeomorphism between Γ(0) and Γ(t ). The space-time surface G is defined by G := Γ(t ) ×{t }. The material velocity T T t ∈[0,T ] v : G → R of Γ(t ) is defined by ∂G v(G(x , t ), t ) = (x , t ), ∀x ∈ Γ(0), t ∈[0, T ]. (2) 0 0 0 ∂t 1 2 3 Vice-versa, if v ˜ ∈ C ([0, T ], C (R )) is an extension of v, i.e. v ˜ (G(x , t ), t ) = v(G(x , t ), t ) for x ∈ Γ(0) and t ∈[0, T ], the mapping G (and thus the time-dependent 0 0 surface Γ(t )) is recovered by solving, for each x ∈ Γ(0), the Cauchy problem 123 J Sci Comput ∂G (x , t ) = v ˜ (G(x , t ), t ), t ∈[0, T ], 0 0 (3) ∂t G(x , 0) = x . 0 0 For t ∈[0, T ] and δ> 0, let U (t ) be the open neighbourhood of Γ(t ) defined by U (t ) := {(x, t ) ∈ R ×[0, T]:|d(x, t )| <δ}. (4) We recall from [10] the following property. Lemma 1 (Fermi coordinates) For any t ∈[0, T ] there exists δ> 0 such that for any x ∈ Γ(t ) there exists a unique a(x) ∈ Γ(t ) that fulfils x = a(x) + d(x)n(a(x)), (5) Hence, for t ∈[0, T ], every point x ∈ U (t ) can be described by d(x) and a(x), which are called normal coordinates or Fermi coordinates of x. The function a(·, t ) : U (t ) → Γ(t ) is called the normal projection onto Γ(t ). For any sufficiently smooth g : G → R and g : G → R ,let ∇ g, Δ g and ∇ · g T T Γ(t ) Γ(t ) Γ(t ) denote the tangential gradient, the Laplace-Beltrami operator and the tangential divergence of g on Γ(t ), respectively (see [10] for the details). In the following, we will write Γ instead of Γ(t ) to simplify the notation. Furthermore, let ∂ g denote the material derivative of g defined by ∂g ˜ ∂ g := + v ·∇˜ g, (6) ∂t where ∇ is the standard gradient in R and g ˜ is any differentiable extension of g defined on a neighborhood of G . Definition (6)is intrinsic, i.e. it does not depend on the choice of the extension g ˜ (see [11] for further details). Let us recall some basic results from [10]. Lemma 2 (Integration by parts) If g : G → R is sufficiently smooth and t ∈[0, T ],it holds that g · μ = ∇ · g − (g · n)(∇ · n), t ∈[0, T ], (7) Γ Γ ∂Γ (t ) Γ(t ) Γ(t ) where μ : ∂Γ (t ) → R is the outward conormal unit vector on ∂Γ (t ),i.e.normalto ∂Γ (t ) and tangent to Γ(t ). Specifically, if g is tangent to Γ ,i.e. g · n = 0, it holds that g · μ = ∇ · g, t ∈[0, T ]. (8) ∂Γ (t ) Γ(t ) Lemma 3 (Transport formula) For sufficiently smooth g : G → R, it holds that g = ∂ g + g∇ · v , t ∈[0, T ]. (9) dt Γ(t ) Γ(t ) Lemma 4 (Green’s formula on surfaces) For sufficiently smooth f, g : G → R, it holds that ∇ f ·∇ g =− f Δ g + f ∇ g · μ, t ∈[0, T ]. (10) Γ Γ Γ Γ Γ(t ) Γ(t ) ∂Γ (t ) 123 J Sci Comput Remark 1 (Surfaces without boundary) Lemmas 2–4 hold on surfaces with or without bound- ary, i.e. ∂Γ (t ) =∅ or ∂Γ (t ) =∅, respectively. Specifically, if ∂Γ (t ) =∅, then the boundary integrals in (7), (8)and (10)vanish. 2.2 Derivation of the Reaction–Diffusion Model in Strong Form Suppose we are given r ∈ N species u : Γ(t ) → R, k = 1,..., r,and let q : Γ(t ) → R , k k k = 1,..., r, be their fluxes tangent to Γ(t ). We recall from [1] the derivation of a system of r equations for u := (u ,..., u ) that accounts for (i) the diffusion on the surface, (ii) the 1 r flux across the boundary (if non-empty) and (iii) the production rates f (u), k = 1,..., r, of the given species. To this end, let R(0) be a portion of Γ(0) and let R(t ) = G(R(0), t ) be the portion of Γ(t ) corresponding to the initial portion R(0). Notice that R(t ) is itself an evolving surface, hence formulae (7) through (10) still hold if Γ(t ) is replaced by R(t ).We consider a mass balance on R(t ) of the form u =− q · μ + f (u), k = 1,..., r, t ∈[0, T ]. (11) k k k dt R (t ) ∂R (t ) R (t ) Since the fluxes q , k = 1,..., r, are tangent to Γ(t ), we can apply the integration-by-parts formula (8) to the first term on the right hand side of (11). Then (11) becomes u =− ∇ · q + f (u), k = 1,..., r, t ∈[0, T ]. (12) k Γ k k dt R (t ) R (t ) R (t ) By applying the transport formula (9) to the left hand side of (12), we obtain ∂ u + u ∇ · v +∇ · q = f (u), k = 1,..., r, t ∈[0, T ], (13) k k Γ Γ k k R (t ) R (t ) where v is the material velocity defined in Sect. 2.1.Since R(t ) is an arbitrary portion, we conclude that ∂ u + u ∇ · v +∇ · q = f (u), k = 1,..., r, t ∈[0, T ]. (14) k k Γ Γ k k We assume q corresponds to a diffusive flux according to Fick’s law as follows: q =−d ∇ u , ∀k = 1,..., r, (15) k k Γ k where d , k = 1,..., r, are positive diffusivity constants. Inserting (15)into(14), we end up with the reaction–diffusion system of the form ∂ u + u ∇ · v = d Δ u + f (u), k = 1,..., r, t ∈[0, T ]. (16) k k Γ k Γ k k 2.3 Invariant Regions and Maximum Principle In this section we recall basic notions concerning invariant regions for systems of the form (16) and conjecture a sufficient condition under which system (16) possesses an invariant region. To this end, we give the following definitions. Definition 1 (Dilation rates)The minimum and maximum instantaneous dilation rates are defined by ∗ ∗ H (t ) := min ∇ · v(x, t ) and H (t ) := max ∇ · v(x, t ), t ∈[0, T ], (17) Γ Γ min max x∈Γ(t ) x∈Γ(t ) 123 J Sci Comput respectively. When the minimum and maximum instantaneous dilation rates coincide, we call ∗ ∗ ∗ H (t ) := H (t ) = H (t ) the instantaneous dilation rate.The minimum and maximum min max global dilation rates are defined by ∗ ∗ ∗ ∗ μ := min H (t ) and μ := max H (t ), (18) min min max max t ∈[0,T ] t ∈[0,T ] respectively. When the minimum and maximum global dilation rates coincide, we call μ := ∗ ∗ μ = μ the global dilation rate. max min Definition 2 (Invariant regions) For the system (16),aregion  in the phase-space R is said to be an invariant region if, whenever the initial condition u is in , the solution u stays in  as long as it exists and is unique. We focus our attention on regions  ⊂ R of hyper rectangular shape, that is to say of the form := [σ , σ ], (19) k=1 where σ ∈ R ∪{−∞} and σ ∈ R ∪{+∞} for all k = 1,..., r. For instance, if σ = 0and k k σ =+∞ for all k = 1,..., r,then  is the positive orthant in R , which means that the solution of the RDS stays positive at all times. Consider the (r − 1)-dimensional hyperfaces :=  ∩{u = σ },  :=  ∩{u = σ }, k = 1,..., r. k k k k k k For k = 1,..., r, we define the constants ∗ ∗ μ if σ ≥ 0, μ if σ ≥ 0, ∗ k ∗ min max k μ := μ := (20) ∗ ∗ μ if σ < 0, μ if σ < 0, max k min ∗ ∗ where μ and μ are defined in (18). min max Next, we conjecture a criterion under which a hyper-rectangle is invariant for system (16). This criterion holds true in the stationary cases (when μ = 0): (i) when Γ is a stationary monodimensional domain in R (see [38]), (ii) when Γ is a stationary k-dimensional domain in R , k ∈ N (see [6]) and (iii) when Γ is a stationary Riemannian manifold without boundary (see [39]). In the case of isotropically evolving flat domains, the invariance of the positive orthant was studied in [41]. To the best of the author’s knowledge, the case of evolving surfaces has not been studied at the continuous level. Hence, we introduce at the continuous level the following conjecture. Conjecture 1 Let  be a hyper-rectangle as in (19) in the phase space of (16) and let f be Lipschitz on .If ∗ r f (u)< μ σ , ∀u ∈  ∩ R , ∀k = 1,..., r, (21) k k k ∗ r f (u)>μ σ , ∀u ∈  ∩ R , ∀k = 1,..., r, (22) k k then  is an invariant region for (16). In particular, when σ =+∞ and σ =−∞,then r r ∩ R and  ∩ R are respectively empty, and so (21) and (22) are automatically fulfilled, respectively. ∗ ∗ Notice that on stationary surfaces, since μ = μ = 0, then conditions (21)–(22) reduce to the inward flux condition considered in [39, Chapter 4]. In Sect. 4 next, we will prove the discrete counterpart of Conjecture 1 obtained by discretizing the RDS in space with the lumped version of the ESFEM method [10] and IMEX Euler in time. 123 J Sci Comput We remark that, on stationary domains, some systems are known to possess an invariant region which do not meet the strict inequalities (21)–(22). For instance, for many mass-action laws, the positive orthant is invariant [4,16] even though the flux of f is tangent to this region, instead of strictly inward. For the scalar case k = 1, the notion of invariant region reduces to the well-known concept of maximum principle. Definition 3 (Maximum principle)For k = 1, the scalar equation (16) fulfils the maximum principle if, for any initial condition, i.e. u(·, 0), the solution u fulfils min 0, min u(y, 0) ≤ u(x, t ) ≤ max 0, max u(y, 0) ,(x, t ) ∈ G . (23) y∈Γ(0) y∈Γ(0) In particular, if the initial condition is nonnegative, u(y, 0) ≥ 0for y ∈ Γ(0), the solution fulfils 0 ≤ u(x, t ) ≤ max u(y, 0), (x, t ) ∈ G . (24) y∈Γ(0) Notice that the maximum principle corresponds to the fact that every (monodimensional) region of the form  =[σ , σ ], with σ ≤ 0, σ ≥ 0, is invariant. 2.4 Derivation of the Variational Formulation Following [1], we derive the variational formulation of system (16). To this end, for each t ∈[0, T ] we multiply equations (16) by the respective test functions ϕ ,...,ϕ ∈ 1 r 2 1 • • 2 −1 L ([0, T ]; H (Γ (t ))) with ∂ ϕ ,...,∂ ϕ ∈ L ([0, T ]; H (Γ (t ))) and integrate over 1 r Γ(t ): (ϕ ∂ u + ϕ u ∇ · v − ϕ f (u)) = d ϕ Δ u , (25) k k k k Γ k k k k Γ k Γ(t ) Γ(t ) 1 −1 for all k = 1,..., r. For the rigorous definition of the Sobolev spaces H and H on a manifold see [21], while for the Bochner spaces L ([0, T ]; B), with B being any Banach space, see [34]. By applying the Green formula (10) to the right hand side of (25) we obtain (ϕ ∂ u + ϕ u ∇ · v − ϕ f (u)) + d ∇ ϕ ·∇ u = d ϕ ∇ u · μ, k k k k Γ k k k Γ k Γ k k k Γ k Γ(t ) Γ(t ) ∂Γ (t ) (26) for all k = 1,..., r. We assume that either Γ(t ) has no boundary, i.e. ∂Γ (t ) =∅,or homogeneous Neumann boundary condition are enforced, i.e. ∇ u · μ = 0on ∂Γ (t ),so Γ k that the last term in (26) vanishes (see Remark 1). Furthermore, by observing that ∂ (ϕ u ) = k k • • ϕ ∂ u + u ∂ ϕ ,(26) becomes k k k k • • ∂ (ϕ u ) = u ∂ ϕ − ϕ u ∇ · v + ϕ f (u) − d ∇ ϕ ·∇ u , k k k k k k Γ k k k Γ k Γ k Γ(t ) Γ(t ) Γ(t ) (27) for all k = 1,..., r. By applying the transport property to the first term on the left hand side of (27), we have ϕ u − u ∂ ϕ = ϕ f (u) − d ∇ ϕ ·∇ u , (28) k k k k k k k Γ k Γ k dt Γ(t ) Γ(t ) Γ(t ) Γ(t ) 123 J Sci Comput for all k = 1,..., r. Therefore, the variational formulation seeks to find u ,..., u ∈ 1 r 2 1 • • 2 −1 L ([0, T ]; H (Γ (t ))) with ∂ u ,...,∂ u ∈ L ([0, T ]; H (Γ (t ))) such that, for each 1 r t ∈[0, T ], u ϕ − u ∂ ϕ + d ∇ u ·∇ ϕ = f (u)ϕ , (29) k k k k k Γ k Γ k k k dt Γ(t ) Γ(t ) Γ(t ) Γ(t ) 2 1 • • 2 −1 for all ϕ ,...,ϕ ∈ L ([0, T ]; H (Γ (t ))) with ∂ ϕ ,...,∂ ϕ ∈ L ([0, T ]; H (Γ (t ))). 1 r 1 r 3 Lumped Evolving Surface Finite Element Method Following the evolving surface finite element method (ESFEM) studied in [1] for the approx- imation of the variational problem (29), we present its lumped counterpart, the lumped evolving surface finite element method (LESFEM). 3.1 Surface Triangulation and Some Definitions Following [9], given h > 0, called meshsize, a triangulation Γ (t ) of the evolving surface Γ(t ) is defined by Γ (t ) = Z (t ), Z (t )∈Z (t ) where Z (t ) is a set of evolving triangles, with x (t ), i = 1,..., N ∈ N, being the overall h i evolving nodes, such that – The nodes evolve with the exact material velocity, i.e. x ˙ (t ) = v(x (t ), t ) for t ∈[0, T ] i i and i = 1,..., N ; –For all t ∈[0, T ] and for any two distinct triangles Z (t ) and Z (t ) in Z (t ),the 1 2 h intersection Z (t ) ∩ Z (t ) is either empty, or a node, or a complete edge; 1 2 –For all t ∈[0, T ] and Z (t ) ∈ Z (t ), Z (t ) ⊂ U (t ),where U is as defined in Lemma 1; h δ δ –Forall t ∈[0, T ], the normal projection a(·, t ) : U (t ) → Γ(t ) defined in Lemma 1 is a one-to-one mapping between Γ (t ) and Γ(t ), i.e. a(Γ (t ), t ) = Γ(t ). h h We assume that, for each t ∈[0, T ], Γ (t ) meets the Delaunay condition, defined as follows. Let e be an edge of Γ (t ) and let Z and Z be the faces of Γ (t ) sharing the edge e.Let α h 1 2 h 1 and α be the angles in Z and Z opposite to e, respectively. The triangulation Γ (t ) is said 2 1 2 h to meet the Delaunay condition if, for all t ∈[0, T ] and for any edge e of Γ (t ), it holds that α + α ≤ π, (30) 1 2 see for instance [16]. The space-time triangulated surface G is defined by h,T G := Γ (t ) ×{t }. h,T h t ∈[0,T ] Since Γ (t ) is piecewise planar, there exists a time-differentiable mapping G : Γ (0) × h h h [0, T]→ R such that, for all t ∈[0, T ], G (Γ (0), t ) = Γ (t ) and, for every facet Z ∈ Z , h h h h G (·, t ) is a diffeomorphism between Z (0) and Z (t ). For a fixed time t ∈[0, T ],let V (t ) be the space of piecewise linear functions on Γ (t ) h h defined by V (t ) := {ϕ ∈ C (Γ (t )) ϕ is linear affine for each Z ∈ Z }. (31) h h | h 123 J Sci Comput Let V be the space of time-dependent, spatially piecewise linear functions defined by V := {ϕ : G → R | ϕ(·, t ) ∈ V (t ) for each t ∈[0, T ]}. (32) h h,T h Given t ∈[0, T ] and a function η ∈ C (Γ (t )), its linear interpolant I η is the unique h h function in V (t ) such that I η(x (t )) = η(x (t )), i = 1,... N . h i i The discrete material derivative of a sufficiently smooth function U ∈ V is defined by ∂U ∂ U := + I (v) ·∇U,(x, t ) ∈ G , h h,T ∂t where v is the material velocity. For our purposes, we define the discrete counterpart of the dilation rates introduced in Definition 1. Definition 4 (Discrete dilation rates)The minimum and maximum discrete instantaneous dilation rates are defined by H (t ) := ess inf ∇ · I (v)(x, t ), min x∈Γ (t ) Γ h h h H (t ) := ess sup ∇ · I (v)(x, t ), t ∈[0, T ], (33) max Γ h x∈Γ (t ) h respectively. When the minimum and maximum discrete instantaneous dilation rates coincide, we call H (t ) := H (t ) = H (t ) the discrete instantaneous dilation rate.The minimum min max and maximum discrete global dilation rates are defined by μ := min H (t ), μ := max H (t ), (34) min min max max t ∈[0,T ] t ∈[0,T ] respectively. When the minimum and maximum discrete global dilation rates coincide, we call μ := μ = μ the discrete global dilation rate. min max For every i = 1,..., N,the i-th Lagrange basis function χ is the unique V function such i h that χ (x (t ), t ) = δ , t ∈[0, T ], i, j = 1,... N , (35) i j ij where δ is the usual Kronecker symbol. The components U ,..., U ∈ V of the spatially ij 1 r h discrete solution may be expressed in the Lagrange basis as U (x, t ) = ξ (t )χ (x, t ), (x, t ) ∈ G , k = 1,..., r. (36) k k,i i h,T i =1 3.2 Preliminary Results on Triangulated Surfaces We recall from [9] the following property of the basis functions. Lemma 5 (Transport property of the basis functions) The basis functions χ ,i = 1,..., N, defined in (35) fulfil ∂ χ = 0, i = 1,..., N . (37) Hence, for the functions U ,k = 1,..., r defined in (36) it holds that ∂ U (x, t ) = ξ (t )χ (x, t ), (x, t ) ∈ G , k = 1,..., r. (38) k k,i i h,T i =1 123 J Sci Comput We recall from [10, Lemma 5.6 and Remark 5.7] the following preliminary result. Lemma 6 (Leibniz formula on triangulated surfaces) For any time-differentiable U, V ∈ V , it holds that • • UV = ∂ UV + U ∂ V + UV ∇ · I (v). (39) Γ h dt Γ (t ) Γ (t ) Γ (t ) Γ (t ) h h h h 3.3 Lumped Evolving Surface Finite Element Method The lumped evolving surface finite element method (LESFEM), applied to the variational formulation (29), seeks to find U ,..., U ∈ V such that 1 r h I (U χ ) − U ∂ χ + d ∇ U ·∇ χ = I ( f (U)χ ), (40) h k i k i k Γ k Γ i h k i h h h dt Γ (t ) Γ (t ) Γ (t ) Γ (t ) h h h h for all k = 1,..., r and i = 1,..., N . Thanks to the transport property (37) of the basis functions, formulation (40) is equivalent to: find U ,..., U ∈ V such that 1 r h I (U χ ) + d ∇ U ·∇ χ = I ( f (U)χ ), (41) h k i k Γ k Γ i h k i h h dt Γ (t ) Γ (t ) Γ (t ) h h h for all k = 1,..., r and i = 1,..., N . The LESFEM method (41) differs from the evolving surface finite element method (ESFEM) in [1] due to the presence of the interpolant operator on the first and the last terms in (41). By expressing U ,..., U according to (36), the matrix 1 r form of (41)is (M ξ ) + d Aξ = Mf (ξ ,...,ξ ), k = 1,..., r, (42) k k k k 1 r dt where A and M are the (time-dependent) stiffness and lumped mass matrices defined by A (t ) = ∇ χ ·∇ χ , i, j = 1,..., N , ij Γ i Γ j h h Γ (t ) χ if i = j, Γ (t ) M (t ) = I (χ χ ) = i, j = 1,..., N , ij h i j 0if i = j, Γ (t ) for t ∈[0, T ], respectively. Notice that the lumped mass matrix M is diagonal. The Delaunay condition (30) holds if and only if A (t ) ≤ 0, i = j, t ∈[0, T ]. (43) ij The above fact, together with the structure of the lumped mass matrix implies that, for any s ≥ 0, −1 (M + sA) M (t ) ≥ 0, t ∈[0, T ]; (44) −1 (M + sA) M (t )1 = 1, t ∈[0, T ]. (45) See [16] for further details. 123 J Sci Comput 3.4 Time Discretisation We are now concerned with the time discretisation of the spatially discrete system (42)arising from the LESFEM. We discretise system (42) by means of the IMEX (IMplicit–EXplicit) Euler method, i.e by treating diffusion implicitly and reaction terms explicitly. To this end, let τ> 0 be a timestep, let t := nτ for all n = 0,..., N with N = ,let A and n T T n n n M be the stiffness and lumped mass matrices at time t , respectively, let (ξ ,...,ξ ) be 1 r n n n the coefficients of the numerical solution at time t ,and let f := f (ξ ,...,ξ ) for each n k k 1 0 0 k = 1,..., r.If (ξ ,...,ξ ) are the coefficients of the spatially discrete initial datum, the 1 r IMEX Euler time discretisation of (42)is n+1 n+1 n n M ξ − M ξ k k n+1 n+1 n n + d A ξ = M f , k = 1,..., r, n = 0,..., N , (46) k T k k or equivalently n+1 n+1 n+1 −1 n n n ξ = (M + τ d A ) M (ξ + τ f ), k = 1,..., r, n = 0,..., N . (47) k T k k 4 Characterisation of Surface Growth The purpose of this section is to characterise surface growth in terms of the material velocity v, with specific regard to isotropic growth. In fact, the lumped mass M (t ) and stiffness A(t ) matrices depend on v. In particular: 1. For an arbitrary triangulated surface that evolves with an arbitrary material velocity, we dM bound the time derivative of the lumped mass matrix in terms of the constants μ min dt and μ defined in (34), i.e. in terms of the divergence ∇ · I (v) of the discrete max Γ h material velocity. We will need this result to prove a sufficient condition for the existence of invariant regions for the semi- and fully-discrete schemes; 2. For an arbitrary smooth or triangulated surface that evolves with an arbitrary material velocity, we characterise the velocity flows ∇ ·v and ∇ · I (v) in terms of the mappings Γ Γ h G and G introduced in Sects. 2.1 and 3.1, respectively; 3. For an arbitrary smooth or triangulated surface that evolves isotropically in space, that is v(x, t ) = S(t )x,(x, t ) ∈ R ×[0, T ], (48) where S :[0, T]→ R is an arbitrary smooth function, we compute exactly ∇ · v and ∇ · I (v) in terms of v. This result will yield a fully practical criterion to detect the Γ h invariant regions of a given RDS in the case of isotropic surface evolution. 4.1 Bounding the Rate of Change of the Mass Matrix in Terms of the Dilation Rates dM In this section we bound the time derivative of M in terms of the discrete dilation rates dt dM defined in (34). To this end, we prove the following characterisation of . dt Lemma 7 (Transport formula for the lumped mass matrix M) The entries of the lumped mass matrix M fulfil the following property I (χ χ ) = I (χ χ )∇ · I (v), ∀i = 1,..., N . (49) h i j h i j Γ h dt Γ (t ) Γ (t ) h h 123 J Sci Comput Proof By choosing U = I (χ χ ) for any i, j = 1,..., N and V = 1 in the Leibniz formula h i j (39)wehave I (χ χ ) = ∂ I (χ χ ) + I (χ χ )∇ · I (v). (50) h i j h i j h i j Γ h dt Γ (t ) Γ (t ) Γ (t ) h h h Now, if i = j,then I (χ χ ) = χ , otherwise, if i = j, I (χ χ ) = 0. Then, from the h i j i h i j transport property (37)wehave ∂ I (χ χ ) = 0for all i, j = 1,..., N . Equation (50) thus h i j implies the desired result (49). In some proofs we will need the following corollary of Lemma 7. Corollary 1 (Consequence of the transport formula for the lumped mass matrix M) The dM diagonal matrix fulfils the estimates dt dm ii μ m (t ) ≤ (t ) ≤ μ m (t ), i = 1,..., N , t ∈[0, T ], (51) min ii max ii dt where μ and μ are defined in (34). min max 4.2 Characterising Velocity Flows in Terms of the Mappings G and G We wish to characterise the continuous and discrete velocity flows ∇ · v and ∇ · I (v) in Γ Γ h terms of the mappings G and G introduced in Sects. 2.1 and 3.1, respectively, for arbitrary smooth or triangulated surfaces that evolve under an arbitrary material velocity. To this end, let Γ(t ) be an arbitrary evolving smooth surface and let (A, X ) be any local parametrisation of Γ(0),where A ⊂ R is an open connected set and X : A → Γ(0) is a differentiable map such that its Jacobian J is full-rank on A.Let B be a measurable subset of A.For all t,the portion X (B) of Γ(0) evolves into G(X (B), t ) ⊂ Γ(t ). By choosing f (x, t ) = 1 for each (x, t ) ∈ G in the transport formula (9), we have 1 = ∇ · v. (52) dt G(X (B),t ) G(X (B),t ) Let e ˆ , i = 1, 2, 3, be the standard basis vectors in R .For (θ , t ) ∈ B ×[0, T ],let J (x, t ) ∈ 2,3 R be the (spatial) Jacobian of the function G(X (θ ), t ) and let ⎡ ⎤ e ˆ e ˆ e ˆ 1 2 3 ˜ ⎣ ⎦ J (θ , t ) := . (53) J (θ , t ) The surface integrals in (52) can be written as integrals on the planar domain B by using the parametrisation G(X (·), t ) : B → G(X (B), t ). Hence, (52) becomes ˜ ˜ det(J (θ , t ))dθ = ∇ · v(G(X (θ ), t )) det(J (θ , t ))dθ, (54) dt B B where · denotes the Euclidean norm in R , or equivalently, ˜ ˜ det(J (θ , t ))dθ = ∇ · v(G(X (θ ), t )) det(J (θ , t ))dθ. (55) dt B B Since (55) holds for any measurable subset B of A, then it holds that ˜ ˜ det(J (θ , t ))=∇ · v(G(X (θ ), t )) det(J (θ , t )). (56) dt 123 J Sci Comput By applying the chain rule, (56) is equivalent to ∇ · v(G(X (θ ), t )) = ln  det(J (θ , t )). (57) dt Given any triangulated surface Γ (t ) (which evolves under the discrete velocity field I (v)), h h we notice that every facet Z (t ) ∈ Z (t ) is smooth and thus parametrisable. For any facet of the initial triangulated surface, Z (0) ∈ Z (0),let (A, X ) be a parametrisation of Z (0),as 2,3 described above. For (θ , t ) ∈ A ×[0, T ],let J (θ , t ) ∈ R be the Jacobian of the function G (X (θ ), t ) and let ⎡ ⎤ e ˆ e ˆ e ˆ 1 2 3 ˜ ⎣ ⎦ J (θ , t ) := . (58) J (θ , t ) A similar reason as above yields the following discrete counterpart of (57). ∇ · I v(G (X (θ ), t )) = ln  det(J (θ , t )). (59) Γ h h h dt Relations (57)and (59) are useful in that they (i) express the velocity flow without tangential derivatives and (ii) can be computed exactly when G and G are known explicitly, e.g. for the isotropic growth as discussed in the next subsection. 4.3 Computing the Dilation Rates for the Isotropic Growth In this section we compute the dilation rates defined in (18)and (34), respectively, for arbitrary surfaces that evolve with the material velocity (48). The velocity field (48) corresponds to the specific case of isotropic evolution, see for instance [8,29] for the case of evolving planar domains and [1,32] for the general case. In particular, for suitable choices of the function S(t ), the growth law (48) admits some specific cases such as uniform, exponential, logistic and periodic growth, see [1,8,29,32]. From (3), it is easy to show that, with the velocity field (48), each x ∈ Γ evolves to the point 0 0 G(x , t ) = exp S(τ )dτ x , t ∈[0, T ], (60) 0 0 therefore the evolution induced by an isotropic growth is a time-dependent dilation of the initial surface. The function φ(t ) := exp S(τ )dτ , t ∈[0, T ], (61) that appears in (60) is known as the growth function, see for instance [29]. Remark 2 (Properties of isotropic growth) Isotropic growth preserves the angles of triangu- lated surfaces. This has two consequences: 1. if Γ (0) meets the Delaunay condition, then Γ (t ) retains the Delaunay condition for all h h t ∈[0, T ]; 2. if A(0) and M (0) are the stiffness and the mass matrices at t = 0, then A(t ) = A(0), M (t ) = φ (t )M (0), t ∈[0, T ]. (62) Hence, in implementations, A(t )and M (t ) need not be computed at each time step. 123 J Sci Comput ∗ ∗ In the following result we compute the dilation rates μ , μ , μ and μ on an min max max min arbitrary smooth or triangulated surface that evolves with the material velocity (48), in terms of S(t ). Theorem 1 (Velocity flow on isotropically growing smooth or triangulated surfaces) Let Γ(t ) be a smooth surface that evolves with the velocity field (48) and let Γ (t ) be the corresponding triangulated surface. Then, the instantaneous dilation rates satisfy H (t ) = 2S(t ) = H (t ), t ∈[0, T ]. (63) Hence, it follows that ∗ ∗ μ = μ = 2min S(t ), and μ = μ = 2max S(t ). (64) min max min max t ∈[0,T ] t ∈[0,T ] Proof Let Γ (t ) be an evolving smooth surface and let (A, X ), J (θ , t ) and J (θ , t ) as defined in Sect. 4.2.From(60)and (61), J (θ , t ) fulfils J (θ , t ) = φ(t )J (θ ), (θ , t ) ∈ A ×[0, T ], (65) 2,3 where J : A → R is the Jacobian of X. It follows that ˜ ˜ det J (θ , t )= φ (t ) det J (θ , 0),(θ, t ) ∈ A ×[0, T ], (66) which implies that ˜ ˜ ln  det J (θ , t )−ln  det J (θ , 0)=2ln φ(t ) =2 S(τ )dτ, (θ, t ) ∈ A×[0, T ]. (67) By differentiating (67), we have ln  det J (θ , t )= 2S(t ), (θ , t ) ∈ A ×[0, T ]. (68) dt By combining (57), (68) and dropping the parameterised coordinates θ,wehave ∇ · v(x, t ) = 2S(t ), (x, t ) ∈ G , (69) Γ T which proves the first equality in (63). Analogously, we prove the second equality in (63)by using (59). This completes the proof. For the uniform, exponential, logistic and periodic growths, the functions φ(t ), S(t ) and the ∗ ∗ dilation rates μ , μ , μ and μ are detailed in Table 1, see also [29,Table 1] for min max min max the case of evolving planar domains, while the corresponding plots are depicted in Fig. 1. 5 Linear Heat Equation and Discrete Maximum Principles We consider, for k = 1, the specific case of linear heat equation on an evolving surface Γ(t ): ∂ u + u∇ · v = dΔ u − βu,β ∈ R, (70) Γ Γ and we prove the semi- and fully-discrete maximum principles for the case when μ + β ≥ min 0. Equation (70) is a special case of the general system (16) that we are interested in. However, we start with this specific case as (i) it provides more insights on the effect of growth on stability, (ii) we are able to prove a better timestep stability condition and (iii) to make the reader familiar with the demonstrative techniques. The following result, that we proved in [16, Theorem 2.1] for the special case of stationary surfaces, addresses the maximum principle at the semi-discrete level. 123 J Sci Comput Table 1 Particular types of growth with their respective growth functions φ(t ), S(t ) functions and constants ∗ ∗ μ , μ , μ ,and μ max min max min Type of growth Growth function φ(t ) S(t)μ = μ μ = μ min max min max r 2r Linear rt + 1 2r rt + 1 rT + 1 Exponential exp(rt ) r 2r 2r K exp(Krt ) rK (K − 1) 2rK (K − 1) Logistic 2r (K − 1) K − 1 + exp(Krt ) K − 1 + exp(Krt ) K − 1 + exp(KrT ) √ √ r sin(rt ) 2r 3 2r 3 Periodic 2 − cos(rt ) − 2 − cos(rt ) 3 3 The constant r > 0 is the growth rate. For the logistic growth, K > 1isthe carrying capacity, i.e. the square root of the ratio between the asymptotical and the initial area of the surface (see [29]) φ (t) 10 1.5 10 1 r=0.5 0 5 10 0 5 10 0 5 10 0 5 10 r=1 r=1.5 1.5 1.5 1.5 1 1 1 S(t) 0.5 0.5 0.5 −1 0 0 0 0 5 10 0 5 10 0 5 10 0 5 10 t t t t Fig. 1 Plots of the growth functions φ(t ) (top row) and the corresponding S(t ) functions (bottom row) listed in Table 1 for K = 2and r = 0.5, 1, 1.5. From left to right: linear, exponential, logistic and periodic growth profiles Theorem 2 (Semi-discrete maximum principle for the linear heat equation (70)) If the veloc- ity field v fulfils μ + β ≥ 0, (71) min with μ as defined in (34), and the triangulation Γ (t ) meets the Delaunay condition for min h all t ∈[0, T ], then the LESFEM solution of (70) fulfils the following discrete maximum principle min 0, min ξ (0) ≤ ξ (t ) ≤ max 0, max ξ (0) , i = 1,..., N , t > 0. j i j j =1,...,N j =1,...,N (72) Proof From (42), the LESFEM spatial discretisation of (70)is (Mξ) + dAξ =−β Mξ. (73) dt By applying the chain rule, (73) becomes dM M ξ + ξ + dAξ =−β Mξ. (74) dt 123 J Sci Comput −1 By multiplying (74) on the left by M ,wehave dM −1 −1 ξ =−dM Aξ − M ξ − βξ. (75) dt −1 −1 dM All we have to prove is that the ODE (75) is dissipative, i.e. −dM A|ξ|− M |ξ|− dt β|ξ|≤ 0.For every t ∈[0, T ], M is diagonal with positive diagonal entries and A fulfils −1 (43) from the Delaunay condition. Then, it follows that −dM A|ξ|≤ 0. Hence, it suffices −1 dM to prove that −(M + β I )|ξ|≤ 0, that is true provided dt dM −1 M + β I ≥ 0. (76) dt By using (51), condition (76)istrueif (μ + β)I ≥ 0, (77) min which holds true from assumption (71). This completes the proof. The next theorem shows the same result for the LESFEM–IMEX Euler full-discretisation of (70), under a timestep restriction. This result holds true for the special case of stationary surfaces (see [16, Theorem 2.2]). Theorem 3 (Fully-discrete maximum principle for the linear heat equation (70)) If the veloc- ity field v fulfils μ + β ≥ 0, ∀t ≥ 0, (78) min with μ as defined in (34), and the triangulation Γ meets the Delaunay condition for min h all t > 0, then the LESFEM–IMEX Euler solution of (70) fulfils the following minimum- maximum principle 0 n 0 min 0, min ξ ≤ ξ ≤ max 0, max ξ , i = 1,..., N , n = 0,..., N , j j j =1,...,N j =1,...,N (79) if the timestep satisfies τβ ≤ 1. (80) In particular, there is no timestep restriction if β ≤ 0. Proof The full-discretisation (47) of the heat equation (70) can be written as n+1 n+1 n+1 −1 n+1 n+1 −1 n n ξ = (M + τdA ) M (M ) M (1 − τβ)ξ , n = 0,..., N . (81) From (44)–(45)wehave n+1 n+1 −1 n+1 (M + τdA ) M ≥ 0, n = 0,..., N , (82) n+1 n+1 −1 n+1 (M + τdA ) M 1 = 1, n = 0,..., N . (83) Then scheme (81) fulfils the discrete maximum principle if n+1 −1 n n n (M ) M (1 − τβ)ξ ≥ min{0,ξ }, n = 0,..., N ; (84) n+1 −1 n n n (M ) M (1 − τβ)ξ ≤ max{0,ξ }, n = 0,..., N . (85) 123 J Sci Comput n+1 −1 n Since (M ) M is diagonal with strictly positive diagonal entries, conditions (84)–(85) aretrueprovided 1 − τβ ≥ 0; (86) n −1 n+1 (1 − τβ)I ≤ (M ) M , n = 0,..., N . (87) Condition (86) is true under assumption (80). In order to prove (87), we need to estimate n+1 n M as a function of M . To this end, by applying Gronwall’s lemma to the first inequality in (51), we have n+1 n τμ min M ≥ M e , n = 0,..., N . (88) By using (88), condition (87)istrueif τμ min 1 − τβ ≤ e . (89) τμ min Let us now define f (τ ) := 1−τβ and g(τ ) = e . These functions fulfil f (0) = g(0) = 1, f is linear and g is non-concave for all μ ∈ R.Then min –if f (0)> g (0), then condition (89) is not fulfilled for any sufficiently small τ . –if f (0) ≤ g (0), then condition (89) is fulfilled for every τ> 0. Now, condition f (0) ≤ g (0) means −β ≤ μ , which is true from assumption (78). This min completes the proof. Remark 3 (Interplay between material velocity and source term) Relation (71) implies that – domain growth (μ > 0) can enable the discrete maximum principle even for β< 0; min – local domain contraction (μ < 0) can prevent the discrete maximum principle even min for β ≥ 0. This interplay is justified by observing that domain evolution implies a dilution effect, explained as follows. By choosing ϕ = 1 in the variational formulation (29) with k = 1 and f (u) =−βu, we obtain u =−β u, t ∈[0, T ]. (90) dt Γ(t ) Γ(t ) If |Γ(t )| denotes the surface area of Γ(t ) and u(t ):= u denotes the mean value |Γ(t )| Γ(t ) of u,(90) becomes (|Γ(t )|u(t )) =−β|Γ(t )|u(t ), t ∈[0, T ]. (91) dt By solving (91)for u(t ),weobtain dt β|Γ(t )|+ |Γ(t )| dt u(t )=− u(t ), t ∈[0, T ]. (92) dt |Γ(t )| By choosing g = 1 in the transport formula (9), we have |Γ(t )|= ∇ · v ≥|Γ(t )|μ , t ∈[0, T ]. (93) min dt Γ(t ) By combining (92)and (93)wehave u(t )≤−(β + μ )u(t ), t ∈[0, T ]. (94) min dt 123 J Sci Comput From (94), the dilution effect arising from surface growth can be interpreted as the dampening or uplifting effect of μ on u(t ). The estimate (94) implies that u(t ) is non-increasing min if β + μ ≥ 0, (95) min which is the continuous counterpart of (71). We conclude that condition (71) is consistent with the interpretation of surface growth in terms of dilution effect. Remark 4 (Interplay between timestep restriction and source term) Relation (80) implies that the timestep restriction needed for guaranteeing the discrete maximum principle is indepen- dent of the material velocity and it only depends on the stiffness parameter β of the source term. In particular, when the source term is nonnegative (i.e. when β ≤ 0), the LESFEM– IMEX Euler fully-discrete scheme unconditionally fulfils the discrete maximum principle. 6 Reaction–Diffusion Systems and Invariant Regions In this section we prove, for the semi- and full-discretisations of RDSs of the form (16), a criterion to test if a hyper-rectangle in the phase-space is invariant. In the case k = 1 of scalar equations, the notion of invariant region collapses to that of minimum-maximum principle, considered in the previous section for the special case of the linear heat equation. We assume that the Delaunay regularity of the mesh is preserved under evolution. For k = 1,..., r,we define the constants μ if σ ≥ 0, μ if σ ≥ 0, min k max μ := μ := (96) μ if σ < 0, μ if σ < 0, max k min where μ and μ are the dilation rates defined in (34). In the following theorem we min max prove that, under similar assumptions of Conjecture 1,  is an invariant region for the solution obtained from the semi-discrete scheme (42). Hence, the following theorem extends [16, Theorem 3.3] to the case of evolving surfaces. Theorem 4 (Invariant rectangles for (42)) Let  be a hyper-rectangle as in (19) in the phase space of (42),let f be Lipschitz on . If the triangulation Γ (t ) satisfies the Delaunay condition for all t ≥ 0 and f (U)< μ σ , ∀U ∈  ∩ R , ∀k = 1,..., r, (97) k k k k f (U)>μ σ , ∀U ∈  ∩ R , ∀k = 1,..., r, (98) k k then  is an invariant region for (42). Proof The semi-discrete method (42) can be written, after applying the chain rule to the term d −1 (M ξ ) and multiplying on the left by M as dt dM −1 −1 ˙ ¯ ¯ ξ =−d M Aξ + f (ξ ,...,ξ ) − M ξ , k = 1,..., r. (99) k k k k 1 r k dt Since M is diagonal with positive diagonal entries and A ≤ 0for i = j from the assumption ij of Delaunay regularity, proceeding as in the proof of [16, Theorem 3.3], it suffices to verify that, for all (U ,..., U ) ∈ , k = 1,..., r,and i = 1,..., N , 1 r 123 J Sci Comput dm ii −1 f (U ,..., σ ,..., U ) − m σ < 0, (100) k 1 k r k ii dt dm ii −1 f (U ,...,σ ,..., U ) − m σ > 0, (101) k 1 r k ii k dt where σ and σ are as in (19). Using relation (51), conditions (100)–(101) hold if, for all k = 1,..., r, f (U ,..., σ ,..., U ) − μ σ < 0, (102) k 1 k r k f (U ,...,σ ,..., U ) − μ σ > 0, (103) k 1 r k k with μ and μ as in (96), that is true from assumptions (97)–(98). This completes the proof. The following theorem provides a sufficient condition for regions to be invariant for the LESFEM–IMEX Euler scheme (47) and extends [16, Theorem 3.4]. In contrast to the semi- discrete case, we relax the strict inequalities (21)–(22) with conditions (105)–(106), in which we use the perturbed dilation rates μ ˜ and μ given by k k τμ μ if σ ≥ 0, max ⎪ min k ⎨ e − 1 if σ > 0, μ ˜ := μ := (104) τμ k max k τ e − 1 ˜ μ if σ ≤ 0, if σ < 0, k min respectively, and μ and μ are defined in (34). Observe that μ ˜ → μ and μ → μ min max k k k as τ → 0. Theorem 5 (Invariant rectangles for (47)) Let  be a hyper-rectangle as in (19) in the phase space of (42),let f be Lipschitz on . If the triangulation Γ (t ) meets the Delaunay condition for all t ≥ 0 and f (U) ≤ σ μ ˜ , ∀U ∈  ∩ R , ∀k = 1,..., r, (105) k k k k f (U) ≥ σ μ , ∀U ∈  ∩ R , ∀k = 1,..., r, (106) k k k k then  is an invariant region for (47) if the timestep τ fulfils τ L ≤ 1, k = 1,..., r, (107) where, for all k = 1,..., r, L is the Lipschitz constant of f . k k Proof The fully-discrete scheme (47) can be written as n+1 n+1 n+1 −1 n+1 n+1 −1 n n n ξ = (M + τdA ) M (M ) M (ξ + τ f ), (108) k k n ∈ N ∪{0}, k = 1,..., r. Since the mesh meets the Delaunay assumption at all times, the matrix properties (82)–(83) hold. Then, it suffices to prove that n+1 −1 n n n σ 1 ≤ (M ) M (ξ + τ f ) ≤ σ 1, k = 1,..., r, (109) k k k where 1 is the column vector of ones. We will prove the two inequalities in (109) in turn. From (51), the inequality on the right side of (109) holds true if n n τ μ ξ + τ f ≤ σ e 1, k = 1,..., r, (110) k k 123 J Sci Comput with μ as defined in (96). Suppose σ ≥ 0. From assumption (105) we can estimate f as follows n n f ≤ σ μ + L (σ 1 − ξ ), k = 1,..., r. (111) k k k k k k From (111), condition (110) holds true provided n τ μ ξ (1 − τ L ) + τ μ σ + τ L σ 1 ≤ σ e 1, k = 1,..., r. (112) k k k k k k k From assumption (107), since ξ ≤ σ 1,then(112) holds true if τ μ ˜ k σ (1 − τ L ) + τ μ σ + τ L σ ≤ σ e , k = 1,..., r, (113) k k k k k k k that is to say τ μ 1 + τ μ ≤ e , k = 1,..., r, (114) which holds true for each τ ∈ R. Suppose, instead, σ < 0. From assumption (105) we can estimate f as follows τ μ e − 1 n n f ≤ σ + L (σ 1 − ξ ), k = 1,..., r. (115) k k k k k From (115), condition (110) holds true provided n τ μ τ μ k k ξ (1 − τ L ) + σ (e − 1 + τ L )1 ≤ σ e 1, k = 1,..., r. (116) k k k k From assumption (107), since ξ ≤ σ 1,then(116) holds true if τ μ τ μ k k σ (1 − τ L )1 + σ (e − 1 + τ L )1 ≤ σ e 1, k = 1,..., r. (117) k k k k k As (117) always holds with the equality, we conclude that the second inequality in (109)is true under assumptions (105)and (107). Similarly, the inequality on the left side of (109) holds under assumptions (106)and (107). This completes the proof. Remark 5 (Sharper timestep restriction) In the specific case of the linear heat equation (70), estimate (80) in Theorem 3 is sharper than estimate (107) in Theorem 5. In fact, since the Lipschitz constant of the source term is L =|β|, the timestep restriction (107) is fulfilled for τ |β|≤ 1, that is more restrictive than condition (80). 7 Velocity-Induced Invariant Regions for RD Models Now, we consider two different RDSs that are well-known in the literature and prove, at the discrete level, the existence of discrete invariant hyper-rectangles for these RDSs, depending on the global discrete dilation rates μ and μ defined in (18). The results in this section min max are confined to the spatially discrete level, but from Conjecture 1, we claim that the same results holds at the continuous level. In the special case of stationary surfaces (i.e. when μ = μ = 0), we obtain invariant hyper-rectangles that have been studied in the min max literature (see [4,6,19]). It is worth remarking that, even though we consider two RD models for illustrative purposes, the following analysis can be easily extended to other types of RDSs. To mention two examples, for the well-known Hodgkin-Huxley model and for the DIB model for electrodeposition considered in [24,25], the invariant region study is carried out in [14]. 123 J Sci Comput 7.1 RDS with Activator-Depleted Kinetics Let us consider an RDS with the well-known non-dimensional activator-depleted kinetics, also known as Schnakenberg or Brusselator kinetics (see for instance [1,37]), on evolving surfaces • 2 ∂ u + u ∇ · v − Δ u = f (u , u ) := γ(a − u + u u ), 1 1 Γ Γ 1 1 1 2 1 2 (118) • 2 ∂ u + u ∇ · v − dΔ u = f (u , u ) := γ(b − u u ), 2 2 Γ Γ 2 2 1 2 2 where a, b and γ are positive parameters and d is a positive diffusion rate. The model describes a system of two interacting chemicals, in which u ≥ 0and u ≥ 0 are the respective 1 2 concentrations. For this reason, we focus our attention on invariant regions contained in the positive ortant. In the following theorem we prove that: (i) the positive orthant is invariant for (118) regardless of μ and μ . At the continuous level, the result holds in the specific min max case of stationary planar domains, see [4]. (ii) when μ > 0, the model possesses invariant min stripes (depending on μ ) in the positive orthant. min Theorem 6 (Velocity-induced invariant regions for the activator-depleted model (118)) For the LESFEM spatial discretisation of (118), the following statements hold: 1. For any value of the constants μ , and μ defined in (34), the positive orthant min max + 2 := [0, +∞[ is invariant. 2. If μ > 0 and σ is a constant such that min 2 γ b σ ≥ , (119) min then the stripe  =[0, +∞[×[0, σ ] is invariant. Proof In order to prove Statements (1) and (2) we have to verify conditions (21)–(22). For Statement (1), we observe that + + –  := {0}×[0, σ ]⊂  := {0}×[0, +∞[ and, for (u , u ) ∈  ,wehave 2 1 2 1 1 f (u , u ) = f (u , u ) = γ a > 0; 1 2 1 1 2 + + –  := [0, σ ]×{0}⊂  := [0, +∞[×{0} and, for (u , u ) ∈  ,wehave 1 1 2 2 2 f (u , u ) = f (u , u ) = γ b > 0. 1 2 2 1 2 This proves Statement (1). For Statement (2), let μ > 0 and we assume for the moment min that the strict inequality holds in (119). Then the set  := [0, +∞[×{σ } is contained in 1 2 the region γ a − (γ + μ )u min (u,v) ∈ R |u > 0,v > , γ u in which f (u , u ) := f (u , u ) − μ u < 0. This proves Statement (2) when the strict 1 1 2 1 1 2 min 1 inequality holds in (119). Otherwise, observe that =[0, +∞[×[0, σ ]= [0, +∞[×[0, σ + ε], (120) 2 2 ε>0 i.e.  is the intersection of invariant regions and is thus invariant. This completes the proof of Statement (2). 123 J Sci Comput 7.2 RDS with Thomas kinetics Let us consider an RDS with the non-dimensional Thomas kinetics (see for instance [31,p. 78]), on evolving surfaces u u • 1 2 ∂ u + u ∇ · v − Δ u = f (u , u ) := γ a − u − ρ , 1 1 Γ Γ 1 1 1 2 1 2 1+u +Ku (121) • u u ⎪ 1 2 ⎩ ∂ u + u ∇ · v − dΔ u = f (u , u ) := γ α(b − u ) − ρ , 2 2 Γ Γ 2 2 1 2 2 1+u +Ku where α, a, b, γ , K and ρ are positive constants and d is a positive diffusion rate. The model describes a system of two interacting chemicals, in which u ≥ 0and u ≥ 0are 1 2 the respective concentrations. For this reason, we focus our attention on invariant regions contained in the positive orthant. Theorem 7 (Velocity-induced invariant regions for the Thomas model (121)) For the LESFEM spatial discretisation of (121), the following statements hold: 1. For any value of the constants μ , and μ defined in (34), the positive orthant min max + 2 := [0, +∞[ is invariant. 2. If μ > −γ min(1,α) and σ and σ are two constants such that min 1 2 γ a γαb σ ≥ , σ ≥ , (122) 1 2 γ + μ γα + μ min min then the region  =[0, σ ]×[0, σ ] is invariant. 1 2 Proof To prove Statements (1) and (2), we have to verify conditions (21)–(22). For Statement (1), observe that –for (u , u ) ∈  := {0}×[0, σ ],wehave f (u , u ) = f (u , u ) = a > 0; 1 2 2 1 2 1 1 2 –for (u , u ) ∈  := [0, σ ]×{0},wehave f (u , u ) = f (u , u ) = αb > 0. 1 2 1 1 2 2 1 2 This proves Statement (1). For Statement (2), let μ > −γ min(1,α) and we assume for min the moment that the strict inequalities hold in (122). Then, observe that • the set  := {σ }×[0, σ ] is contained in the region 1 1 2 1 + u + Ku (u,v) ∈ R |u > 0,v > (γ a − (μ + γ)u) , min γρu in which f (u , u ) := f (u , u ) − μ u < 0; 1 2 1 1 2 min 1 • the set  := [0, σ ]×{σ } is contained in the region 2 1 2 γαb(1 + u + Ku ) (u,v) ∈ R |u > 0,v > , γρu + (γ α + μ )(1 + u + Ku ) min in which f (u , u ) := f (u , u ) − μ u < 0. 1 2 2 1 2 min 2 This proves Statement (2) when the strict inequalities hold in (122). Otherwise, we have that =[0, σ ]×[0, σ ]= [0, σ + ε]×[0, σ + ε], (123) 1 2 1 2 ε>0 i.e.  is the intersection of invariant regions and is thus invariant. 123 J Sci Comput 8 Numerical Examples The purpose of this section is to provide two numerical examples in which we (i) estimate the experimental order of convergence of the LESFEM and (ii) experimentally show the ability of Theorem 5 to find invariant regions of RDSs on evolving surfaces at the discrete level. 8.1 Numerical Example 1: Linear Heat Equation on an Evolving Sphere In this example, we wish to estimate the experimental order of convergence of the LESFEM. As a test problem, we consider the linear heat equation given by ∂ u + u∇ · v − Δ u = u, (124) Γ Γ We choose T = 1 to be the final time. The initial domain Γ(0) is the unit sphere S ,that evolves under the velocity field v(x, t ) := ,(x, t ) ∈ R ×[0, 1], (125) t + 1 and undergoes linear growth for r = 1, see Table 1. In particular, the domain Γ(t ) at time t ∈[0, T ] is a sphere whose radius is given by the growth function φ(t ) = t + 1and the 2r minimum dilation rate fulfils μ = = 1(seeTable 1). min rT +1 In order to determine the experimental order of convergence, we consider the analytical solution to (124)given by xyz t u(x , y, z, t ) = exp t − 2log(t + 1) − ,(x , y, z, t ) ∈ R ×[0, 1]. (t + 1) t + 1 (126) See “Appendix” for the derivation of the solution (126). The constants β =−1and μ = 1fulfill(78)and (80) for each τ> 0. Hence, the min LESFEM–IMEX Euler solution to (124) fulfils a discrete maximum principle unconditionally on τ . In order to appreciate the quadratic convergence in space we solve the problem on a sequence of eight Delaunay meshes Γ , i = 1,..., 8, whose mesh sizes at t = 0fulfil h (0) i −1 h (0) = 0.4013 and h (0) ≈ , i = 2,..., 8. The corresponding timesteps τ fulfil 1 i i τ = 4e-2 and τ = τ , i = 2,..., 8. In Fig. 2 we show a sequence of snap shots 1 i i −1 i −1 of the evolution of the numerical solution on the finest mesh Γ . The experimental order of ∞ 2 convergence is computed by measuring the error, in L ([0, T ], L (Γ (t ))) norm, between the numerical solution U and the piecewise linear interpolant I (u) of the exact solution. The result is shown in Fig. 3: the convergence is experimentally optimal in that it is quadratic in the meshsize and linear in the timestep. 8.2 Numerical Example 2: RDS with Thomas Kinetics and Invariant Regions In this example, we show that the LESFEM–IMEX Euler preserves the invariant regions of RDSs on evolving surfaces. Let us consider the following RDS with Thomas kinetics u u • 1 2 ⎨ ∂ u + u ∇ · v − d Δ u = f (u , u ) := γ a − u − ρ ; 1 1 Γ 1 Γ 1 1 1 2 1 1+u +Ku (127) • u u ⎪ 1 2 ⎩ ∂ u + u ∇ · v − Δ u = f (u , u ) := γ α(b − u ) − ρ , 2 2 Γ Γ 2 2 1 2 2 1+u +Ku 123 J Sci Comput Fig. 2 Numerical Example 1 on the linear heat equation (124): numerical solution at different times obtained on the finest mesh Γ with N = 16962 gridpoints, initial meshsize h = 3.542e−2 and timestep τ = 3.116e−4. 8 8 Plotted values range from −0.1924 (blue) to 0.1924 (red) (Color figure online) −2 2.2 EOC −3 E ∞ 2 L (L ) 1.8 −4 1.6 −1 −1 ∞ 2 Fig. 3 Numerical Example 1 on the linear heat equation (124): Error in L (0, T , L (Γ (t ))) norm (left) and experimental rate of convergence (right). The quadratic convergence in space is optimal considered in Sect. 7.2, with a = 150, b = 100, ρ = 13, K = 0.05 as in [30]and γ = 1for illustrative purposes. With these parameters, system (127) admits, in the absence of domain growth, the homogeneous steady state P ≈ (37.7382, 25.1588) calculated using Newton- Raphson method on a stationary domain. Notice that the diffusion coefficient for u has been normalised to 1 for convenience, and we choose d = 0.01. The initial domain Γ(0) is the Dupin ring cyclide considered in [18, Appendix B], rescaled for convenience, given by ⎧ ⎫ ⎨ 2 ⎬ 261 39 3249 3 2 2 2 2 D := (x , y, z) ∈ R : 9(x + y + z ) + − 4 6x − − y = 0 , ⎩ 100 10 25 ⎭ an orientable surface without boundary that is topologically equivalent to a torus. The surface rK (K −1) evolves under the velocity field (48) with S(t ) = > 0, t ∈[0, T ], with r = 0.2, K −1+exp(Krt ) K = 3 and undergoes a logistic growth, see Table 1. As final time we choose T = 100. From 123 J Sci Comput Fig. 4 Numerical Example 2 on the Thomas RDS (127): invariant region for the LESFEM–IMEX Euler full discretisation of (127) under the timestep restriction τ ≤ 1.676e−3. This region is obtained by considering the nullclines of the modified kinetics defined in (129)–(130) Theorem 1, it follows that the corresponding dilation rates fulfil ∗ ∗ μ = μ ≈ 0, and μ = μ = 0.8. (128) min max min max From Theorem (7), model (127) possesses arbitrarily large bounded invariant regions con- tained in the positive orthant. Let =[σ , σ ]×[σ , σ ]≈[0.3366, 126.4194]×[13.2938, 45.8182] 1 2 1 2 be the smallest bounded region that (i) contains P and (ii) meets the modified inward flux conditions (105)–(106) with τ =¯ τ := 2e-3 and thus (105)–(106) hold for all 0 <τ ≤¯ τ.As illustrated in Fig. 4,  is obtained by considering the nullclines of the modified kinetics f (u , u ) := f (u , u ) − μ u , k = 1, 2; (129) 1 2 k 1 2 min k τμ ¯ max e − 1 f (u , u ) := f (u , u ) − u , k = 1, 2. (130) 1 2 k 1 2 k τ¯ It is easy to see that, on a region of the form  =[σ , σ ]×[σ , σ ], the Lipschitz 1 2 1 2 constants L and L of the kinetics f and f of (127)fulfil 1 2 1 2 " " 2 2 2 2 ˆ ˆ L ≤ L := γ (1 + ρσ ) + ρ , and L ≤ L := γ (ρσ ) + (α + ρ) . 1 1 2 2 2 2 Hence, the timestep restriction (107) becomes τ ≤ ≈ 1.676e-3. It follows that ˆ ˆ max(L ,L ) 1 2 is invariant for the LESFEM–IMEX Euler full discretisation under the timestep restriction 0 <τ ≤ min τ, ¯ ≈ 1.676e-3. We choose τ = 1e-3. The region  is smaller ˆ ˆ max(L ,L ) 1 2 than the invariant region provided in Theorem 7, which has a simple analytical expression but is not optimal. The following initial condition u (x , y, z) = σ + (σ − σ )ψ (x , y, z); 1,0 1 1 1 u (x , y, z) = σ + (σ − σ )ψ (x , y, z), 2,0 2 2 2 where ψ(x , y, z) := 1 − 25(min(|y|, )) fulfils (u (x), u (x)) ∈  for all x ∈ Γ(0) 1,0 2,0 and is shown in the first snap shot of Fig. 5. We solve the problem with τ = 1e-3 on a sequence of seven meshes Γ (t ), i = 1,..., 7, with decreasing initial meshsizes h (0), h,i i 123 J Sci Comput Fig. 5 Numerical Example 2 on the Thomas RDS (127): snap shots of the U component of the LESFEM– IMEX Euler numerical solution at different times Table 2 Numerical Example 2 on the Thomas RDS (127): for ESFEM,  is not invariant for any i = 1,..., 7, as the method violates the minimum of U for i = 5, 6, 7 and the minimum of U for each i = 1,..., 7 1 2 ih (0) min η min η min η min η i t ∈[τ,T ] t ∈[τ,T ] 1 t ∈[τ,T ] t ∈[τ,T ] 2 1 2 1 1.190e+00 8.482e−02 1.841e−01 − 8.001e−01 5.214e−01 2 8.537e−01 8.005e−02 5.281e−01 − 5.318e−01 6.397e−01 3 5.898e−01 5.860e−02 1.759e−01 − 5.083e−01 7.626e−01 4 4.273e−01 2.259e−02 1.793e−01 − 6.559e−01 5.147e−01 5 3.011e−01 − 7.288e−02 1.791e−01 − 6.338e−01 8.181e−01 6 2.114e−01 − 3.555e−01 1.241e−01 − 5.017e−01 4.297e−01 7 1.531e−01 − 5.376e−01 1.321e−01 − 3.816e−01 6.152e−01 i = 1,..., 7, with both the LESFEM–IMEX Euler and ESFEM-IMEX Euler methods. Snap shots of the LESFEM–IMEX Euler numerical solution obtained on the finest mesh Γ at h,7 different times are shown in Fig. 5. In particular, at the final time T = 100 (see the last snap shot of Fig. 5), the surface is stationary up to machine precision and the numerical solution has reached a stationary pattern. For a given numerical solution (U , U ) on the mesh Γ , 1 2 h,i i = 1,..., 7, consider the following functions η (t ) := min (σ − U (x, t )), η (t ) := min (U (x, t ) − σ ), k = 1, 2. k k k x∈Γ (t ) x∈Γ (t ) h,i h,i These functions are the oriented distances of the numerical solution (U , U ) from the edges 1 2 of . If the oriented distances η and η , k = 1, 2, stay positive at all times, it means that (U , U ) is in  at all times. For all i = 1,..., 7, we show the minima over the time interval 1 2 [τ, T ] (i.e. excluding the initial data) of η and η , k = 1, 2, for both the ESFEM-IMEX Euler and LESFEM–IMEX Euler methods in Tables 2 and 3, respectively. We observe that ESFEM-IMEX Euler violates  for all i = 1,..., 7, while LESFEM–IMEX Euler preserves for all i = 1,..., 7. 123 J Sci Comput Table 3 Numerical Example 2 on the Thomas RDS (127): for LESFEM, the region  is invariant for all i = 1,..., 7 ih (0) min η min η min η min η i t ∈[τ,T ] t ∈[τ,T ] 1 t ∈[τ,T ] t ∈[τ,T ] 2 1 2 1 1.190e+00 1.060e−01 1.738e−01 7.6017752487788e−02 2.665e−01 2 8.537e−01 1.060e−01 4.796e−01 7.6017752487775e−02 3.676e−01 3 5.898e−01 1.060e−01 1.958e−01 7.6017752487774e−02 5.336e−01 4 4.273e−01 1.060e−01 2.163e−01 7.6017752487768e−02 6.679e−01 5 3.011e−01 1.060e−01 1.934e−01 7.6017752487765e−02 9.708e−01 6 2.114e−01 1.060e−01 1.909e−01 7.6017752487761e−02 9.401e−01 7 1.531e−01 1.060e−01 1.899e−01 7.6017752487756e−02 9.367e−01 The minima of η (and thus the minima of U ) coincide up to machine precision 9 Conclusions In this paper we have considered a lumped evolving surface finite element method for the spatial discretisation of reaction–diffusion systems on evolving surfaces, by extending substantially the counterpart on stationary surfaces studied in [16]. We have obtained a fully- discrete scheme by applying the IMEX Euler timestepping to the spatially discrete method. In Sect. 4 we have presented a characterisation of the tangential flow of the continuous and discrete material velocities on evolving smooth and triangulated surfaces, respectively. As a consequence we have obtained, in Theorem 1, an explicit expression for the dilation rates for the special case of isotropic growth. In Theorems 2 and 3 of Sect. 5 we have presented sufficient conditions for the linear heat equation on evolving surfaces to fulfil the maximum principle at the semi- and fully-discrete levels, respectively. In particular, at the fully-discrete level, no timestep restriction is needed if the source parameter β in (70) is nonpositive. In Theorems 4 and 5 of Sect. 6 we provided sufficient conditions under which a hyper-rectangle in the phase space of an RDS on evolving surfaces is invariant at the semi- and fully-discrete levels, respectively. In particular, at the fully-discrete level, a timestep restriction depending on the Lipschitz constants of the kinetics is needed. In Sect. 7 we classified some families of invariant regions for two well-known RD models on evolving surfaces: the activator-depleted and the Thomas models. Two numerical examples are presented in Sect. 8 that experimentally show (i) the optimal convergence (i.e. quadratic in space and linear in time) of the proposed method for the linear heat equation on a linearly evolving sphere and (ii) the existence of a bounded invariant region for an RDS with Thomas kinetics on a logistically growing Dupin ring cyclide at the discrete level. The mathematical and numerical analysis of more complicated settings such as non- isotropic chemically driven is a non-trivial exercise that requires new mathematical tools and techniques. Such a study is beyond the scope of this paper. In order to provide a first step in this direction, we have considered the simplest form of growth, uniform isotropic, which has allowed us to state precisely the conditions necessary for the preservation of invariant regions as well as stating the conditions for positivity of solutions. In future studies we wish to extend the method by considering an ALE approach as in [13] in order to preserve good mesh properties under evolution. This would involve the study of convection RDSs on evolving surfaces, where the convection term arises from the ALE mapping. 123 J Sci Comput Acknowledgements The authors (MF, AM, IS CV) would like to thank the Isaac Newton Institute for Math- ematical Sciences for its hospitality during the programme [Coupling Geometric PDEs with Physics for Cell Morphology, Motility and Pattern Formation] supported by EPSRC Grant Number EP/K032208/1. MF’s and IS’s research work has been performed under the auspices of the Italian National Group for Scientific Calculus (GNCS-INdAM). This work (AM, CV) is partly supported by the EPSRC Grant Number EP/J016780/1 and the Leverhulme Trust Research Project Grant (RPG-2014-149). AM acknowledges funding from the European Union Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant Agree- ment No. 642866, the Commission for Developing Countries, and was partially supported by a grant from the Simons Foundation. AM is a Royal Society Wolfson Research Merit Award Holder funded generously by the Wolfson Foundation. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 Interna- tional License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Appendix: Derivation of the Analytical Solution to Experiment 8.1 2 2 Let ψ : S → R be an eigenfunction of the Laplace-Beltrami operator on S with eigen- value λ ∈ R, i.e. Δ 2 ψ = λψ. Let the surface evolve with an arbitrary isotropic velocity field (48)and let φ(t ) be the respective growth function, as defined in (61). For each t ∈[0, T ], Γ(t ) is a sphere of radius φ(t ). By using the chain rule, it is easy to show that x λ x Δ ψ = ψ ,(x, t ) ∈ G , (131) Γ(t ) T φ(t ) φ (t ) φ(t ) i.e. ψ is an eigenfunction of Δ and its eigenvalue decays quadratically with the Γ(t ) φ(t ) growth function. We look for solutions to (124)ofthe form u(x, t ) = η(t )ψ ,(x, t ) ∈ G , (132) φ(t ) where η(t ) is an unknown time-dependent coefficient. By using (63), (131)and (132)in (124), and cancelling ψ on both sides, we have φ(t ) $ % λd η( ˙ t ) = η(t ) −β − 2S(t ) + , t ∈[0, T ], (133) φ (t ) which leads to $ % λd η(t ) = η(0) exp −β − 2S(τ ) + dτ, t ∈[0, T ]. (134) φ (τ ) In (134) we choose • a linear growth with r = 1 (hence S(t ) = and φ(t ) = t + 1for t ∈[0, 1],see t +1 Table 1); • η(0) = 1, β =−1and d = for illustrative purposes. The profile η(t ) defined in (134) thus becomes λt η(t ) = exp t − 2ln(t + 1) + , t ∈[0, 1]. (135) 12(t + 1) 123 J Sci Comput By combining (132) with (135), we have x λt u(x , y, z, t ) = ψ exp t − 2ln(t + 1) + ,(x , y, z, t ) ∈ G . (136) t + 1 12(t + 1) Finally, we choose the eigenfunction ψ(x , y, z) := xyz, (x , y, z) ∈ S , whose eigenvalue is λ =−12. Hence, (136) becomes xyz t u(x , y, z, t ) = exp t − 2ln(t + 1) − ,(x , y, z, t ) ∈ G . (137) (t + 1) t + 1 References 1. Barreira, R., Elliott, C.M., Madzvamuse, A.: The surface finite element method for pattern formation on evolving biological surfaces. J. Math. Biol. 63(6), 1095–1119 (2011) 2. Bertalmıo, M., Cheng, L.T., Osher, S., Sapiro, G.: Variational problems and partial differential equations on implicit surfaces. J. Comput. Phys. 174(2), 759–780 (2001) 3. Blazakis, K.N., Madzvamuse, A., Reyes-Aldasoro, C.C., Styles, V., Venkataraman, C.: Whole cell tracking through the optimal control of geometric evolution laws. J. Comput. Phys. 297, 495–514 (2015) 4. Chellaboina, V., Bhat, S.P., Haddad, W.M., Bernstein, D.S.: Modeling and analysis of mass-action kinetics. IEEE Control Syst. 29, 60–78 (2009) 5. Chen, C, Kublik, C, Tsai, R.: An implicit boundary integral method for interfaces evolving by Mullins– Sekerka dynamics. In: Maekawa, Y., Jimbo, S. (eds.) Mathematics for Nonlinear Phenomena—Analysis and Computation, pp. 1–21. Springer (2017) 6. Chueh, K.N., Conley, C.C., Smoller, J.A.: Positively invariant regions for systems of nonlinear diffusion equations. Indiana Univ. Math. J. 26, 373–392 (1977) 7. Corson, F., Hamant, O., Bohn, S., Traas, J., Boudaoud, A., Couder, Y.: Turning a plant tissue into a living cell froth through isotropic growth. Proc. Nat. Acad. Sci. 106(21), 8453–8458 (2009) 8. Crampin, E.J., Gaffney, E.A., Maini, P.K.: Reaction and diffusion on growing domains: scenarios for robust pattern formation. Bull. Math. Biol. 61(6), 1093–1120 (1999) 9. Dziuk, G., Elliott, C.M.: Finite elements on evolving surfaces. IMA J. Numer. Anal. 27(2), 262–292 (2007) 10. Dziuk, G., Elliott, C.M.: Finite element methods for surface PDEs. Acta Numer. 22, 289–396 (2013) 11. Dziuk, G., Elliott, C.M.: L -estimates for the evolving surface finite element method. Math. Comput. 82(281), 1–24 (2013) 12. Elliott, C.M., Stinner, B., Venkataraman, C.: Modelling cell motility and chemotaxis with evolving surface finite elements. J. R. Soc. Interface 9(76), 3027–3044 (2012) 13. Elliott, C.M., Venkataraman, C.: Error analysis for an ALE evolving surface finite element method. Numer. Methods Partial Differ. Equ. 31(2), 459–499 (2015) 14. Frittelli, M: Numerical Methods for Partial Differential Equations on Stationary and Evolving Surfaces. Ph.D. thesis, Dipartimento di Matematica e Fisica “E. De Giorgi”, Università del Salento (2018) (to appear) 15. Frittelli, M., Madzvamuse, A., Sgura, I., Venkataraman, C.: Lumped finite elements for reaction-cross- diffusion systems on stationary surfaces. Comput. Math. Appl. 74(12), 3008–3023 (2017). https://doi. org/10.1016/j.camwa.2017.07.044 16. Frittelli, M., Madzvamuse, A., Sgura, I., Venkataraman, C.: Preserving invariance properties of reaction– diffusion systems on stationary surfaces. IMA J. Numer. Anal. (2017). https://doi.org/10.1093/imanum/ drx058 17. Frittelli, M., Sgura, I.: Virtual element method for the Laplace–Beltrami equation on surfaces. ESAIM Math. Model. Numer. Anal. (2017). https://doi.org/10.1051/m2an/2017040 18. Fuselier, E.J., Wright, G.B.: A high-order kernel method for diffusion and reaction–diffusion equations on surfaces. J. Sci. Comput. 56(3), 535–565 (2013) 19. González-Olivares, E., Ramos-Jiliberto, R.: Dynamic consequences of prey refuges in a simple model system: more prey, fewer predators and enhanced stability. Ecol. Model. 166, 135–146 (2003) 20. Grande, J., Reusken, A.: A higher order finite element method for partial differential equations on surfaces. SIAM J. Numer. Anal. 54(1), 388–414 (2016) 21. Hebey, E., Robert, F.: Sobolev spaces on manifolds. In: Krupka, D., Saunders, D. (eds.) Handbook of Global Analysis, pp. 375–415 (2008) 123 J Sci Comput 22. Kovács, B.: Computing arbitrary Lagrangian Eulerian maps for evolving surfaces. arXiv preprint arXiv:1612.01701 (2016) 23. Kovács, B., Li, B., Lubich, C., Guerra, C.A.P.: Convergence of finite elements on an evolving surface driven by diffusion on the surface. Numer. Math. 137(3), 643–689 (2017) 24. Lacitignola, D., Bozzini, B., Frittelli, M., Sgura, I.: Turing pattern formation on the sphere for a mor- phochemical reaction–diffusion model for electrodeposition. Commun. Nonlinear Sci. Numer. Simul. 48, 484–508 (2017). https://doi.org/10.1016/j.cnsns.2017.01.008 25. Lacitignola, D., Bozzini, B., Sgura, I.: Spatio-temporal organization in a morphochemical electrode- position model: Hopf and Turing instabilities and their interplay. Eur. J. Appl. Math. 26(2), 143–173 (2015) 26. Lefèvre, J., Mangin, J.F.: A reaction–diffusion model of human brain development. PLoS Comput. Biol. 6(4), e1000,749 (2010) 27. Li, H., Lin, Y., Heath, R.M., Zhu, M.X., Yang, Z.: Control of pollen tube tip growth by a rop gtpase- dependent pathway that leads to tip-localized calcium influx. Plant Cell 11(9), 1731–1742 (1999) 28. Madzvamuse, A., Barreira, R.: Exhibiting cross-diffusion-induced patterns for reaction-diffusion systems on evolving domains and surfaces. Phys. Rev. E 90(4), 043,307 (2014) 29. Madzvamuse, A., Gaffney, E.A., Maini, P.K.: Stability analysis of non-autonomous reaction-diffusion systems: the effects of growing domains. J. Math. Biol. 61(1), 133–164 (2010) 30. Madzvamuse, A., Wathen, A.J., Maini, P.K.: A moving grid finite element method applied to a model biological pattern generator. J. Comput. Phys. 190(2), 478–500 (2003) 31. Murray, J.D.: Mathematical Biology. II Spatial Models and Biomedical Applications (Interdisciplinary Applied Mathematics), vol. 18. Springer, New York (2001) 32. Plaza, R.G., Sanchez-Garduno, F., Padilla, P., Barrio, R.A., Maini, P.K.: The effect of growth and curvature on pattern formation. J. Dyn. Differ. Equ. 16(4), 1093–1121 (2004) 33. Pruyne, D., Bretscher, A.: Polarization of cell growth in yeast. I. Establishment and maintenance of polarity states. J. Cell Sci. 113(3), 365–375 (2000) 34. Quarteroni, A., Valli, A.: Numerical Approximation of Partial Differential Equations, vol. 23. Springer, New York (2008) 35. Ren, K., Tsai, R., Zhong, Y.: An implicit boundary integral method for computing electric potential of macromolecules in solvent. arXiv:1709.08070v4 (2018) 36. Sahlin, P., Jönsson, H.: A modeling study on how cell division affects properties of epithelial tissues under isotropic growth. PLoS ONE 5(7), e11,750 (2010) 37. Schnakenberg, J.: Simple chemical reaction systems with limit cycle behaviour. J. Theor. Biol. 81(3), 389–400 (1979) 38. Smoller, J.: Shock Waves and Reaction Diffusion Equations. Springer, New York (1994) 39. Taylor, M.E.: Partial Differential Equations. III. Springer, New York (2011) 40. Tuncer, N., Madzvamuse, A.: Projected finite elements for systems of reaction–diffusion equations on closed evolving spheroidal surfaces. Commun. Comput. Phys. 21(3), 718–747 (2017) 41. Venkataraman, C., Lakkis, O., Madzvamuse, A.: Global existence for semilinear reaction–diffusion sys- tems on evolving domains. J. Math. Biol. 64(1–2), 41–67 (2012) 42. Yang, F.W., Venkataraman, C., Styles, V., Madzvamuse, A.: A robust and efficient adaptive multigrid solver for the optimal control of phase field formulations of geometric evolution laws. Commun. Comput. Phys. 21(1), 65–92 (2017) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Journal of Scientific Computing Springer Journals

Numerical Preservation of Velocity Induced Invariant Regions for Reaction–Diffusion Systems on Evolving Surfaces

Free
30 pages

Loading next page...
 
/lp/springer_journal/numerical-preservation-of-velocity-induced-invariant-regions-for-h0bp7vnlWG
Publisher
Springer US
Copyright
Copyright © 2018 by The Author(s)
Subject
Mathematics; Algorithms; Computational Mathematics and Numerical Analysis; Mathematical and Computational Engineering; Theoretical, Mathematical and Computational Physics
ISSN
0885-7474
eISSN
1573-7691
D.O.I.
10.1007/s10915-018-0741-7
Publisher site
See Article on Publisher Site

Abstract

J Sci Comput https://doi.org/10.1007/s10915-018-0741-7 Numerical Preservation of Velocity Induced Invariant Regions for Reaction–Diffusion Systems on Evolving Surfaces 1 2 Massimo Frittelli · Anotida Madzvamuse · 1 3 Ivonne Sgura · Chandrasekhar Venkataraman Received: 13 December 2017 / Revised: 4 May 2018 / Accepted: 5 May 2018 © The Author(s) 2018 Abstract We propose and analyse a finite element method with mass lumping (LESFEM) for the numerical approximation of reaction–diffusion systems (RDSs) on surfaces in R that evolve under a given velocity field. A fully-discrete method based on the implicit–explicit (IMEX) Euler time-discretisation is formulated and dilation rates which act as indicators of the surface evolution are introduced. Under the assumption that the mesh preserves the Delaunay regularity under evolution, we prove a sufficient condition, that depends on the dilation rates, for the existence of invariant regions (i) at the spatially discrete level with no restriction on the mesh size and (ii) at the fully-discrete level under a timestep restriction that depends on the kinetics, only. In the specific case of the linear heat equation, we prove a semi- and a fully-discrete maximum principle. For the well-known activator-depleted and Thomas reaction–diffusion models we prove the existence of a family of rectangles in the phase space that are invariant only under specific growth laws. Two numerical examples are provided to computationally demonstrate (i) the discrete maximum principle and optimal convergence for the heat equation on a linearly growing sphere and (ii) the existence of an invariant region for the LESFEM–IMEX Euler discretisation of a RDS on a logistically growing surface. Anotida Madzvamuse A.Madzvamuse@sussex.ac.uk Massimo Frittelli massimo.frittelli@unisalento.it Ivonne Sgura ivonne.sgura@unisalento.it Chandrasekhar Venkataraman cv28@st-andrews.ac.uk Dipartimento di Matematica e Fisica “E. De Giorgi”, Università del Salento, via per Arnesano, 73100 Lecce, Italy Department of Mathematics, School of Mathematical and Physical Sciences, University of Sussex, Pevensey III, 5C15, Brighton BN1 9QH, UK School of Mathematics and Statistics, University of St Andrews, Fife KY16 9SS, UK 123 J Sci Comput Keywords Evolving surface · Dilation rate · Heat equation · Maximum principle · Reaction–diffusion · Invariant region 1 Introduction Reaction–diffusion systems (RDSs) are a well-known class of partial differential equations that arise from the mathematical modelling of numerous phenomena taking place on a space- time domain, see for instance [31]. In many applications the spatial domain is a curved surface rather than a flat domain, which may be time-dependent. Among the applications of RDSs on surfaces we mention brain growth [26], cell migration [3], chemotaxis [12], developmental biology [28], electrodeposition [24] and phase field modeling [42]. The growing interest toward PDEs on evolving surfaces has stimulated the development of several numerical methods for such problems, among which we mention (but not limited to) embedding methods [2], kernel methods [18], implicit boundary integral methods [5,35], surface finite element methods (SFEM) [10] and some of their recent variations and extensions [13,16,17,20,23, 40]. An interesting property of some RDSs is the existence of invariant regions.Aninvariant region is a subset  of the phase-space such that, if the initial condition takes values in ,then the solution of the RDS takes values in  at all times. We recall that, for scalar equations, the well-known notion of maximum principle is equivalent to the invariance of all the regions of the form [0, M ], M > 0. For RDSs on a stationary surface, sufficient conditions have been found for a region to be invariant at the continuous level, see [39]. To the best of the authors’ knowledge, the extension of these results to RDSs on evolving surfaces has not been considered in the literature. In this paper we will focus on surfaces that evolve according to a prescribed material velocity field. More complicated forms of the evolution law are beyond the scope of this manuscript and are a subject of our current studies. For a numerical method it is interesting to understand if the invariant regions of the continuous problem are preserved under discretisation. In [15,16]weprovedthat, forRDSs on stationary surfaces, the lumped surface finite element method (LSFEM), combined with an implicit–explicit (IMEX) Euler method, preserves the invariant regions of the continuous problem. The purpose of the present paper is to extend these results to the case of evolving surfaces, in particular (i) we prove a semi- and a fully-discrete maximum principle for the heat equation with a linear source term and (ii) we provide sufficient conditions under which a region is invariant at the semi- and fully-discrete levels when a lumped evolving surface finite element method (LESFEM) and an IMEX Euler timestepping are considered. In particular we quantify the impact of surface evolution (measured through the dilation rate)onthe existence of invariant regions and we find that surface growth or contraction respectively fosters or inhibits the invariance of a given region in the phase-space. Crucial in our analysis is the assumption that the mesh preserves the Delaunay property under evolution, which is not true for an arbitrary surface evolution law. A class of surface evolution laws for which the Delaunay property is automatically preserved is that of isotropic growth [29], which has biological applications [7,27,33,36]. To the best of the authors’ knowledge, an adaptive strategy for the preservation of the Delaunay property under a generic evolution law is still an open problem. An attempt in this direction is the work in [22]. For the special case of isotropic growth, we provide fully practical sufficient conditions for the existence of invariant regions at the semi- and fully-discrete levels. As an application of our general theory, we classify some classes of invariant regions, depending on the growth rate of 123 J Sci Comput the evolving surface, for two well-known RD models in the literature: the activator-depleted (or Schnakenberg, also known as the Brusselator model) and the Thomas models. Finally, we provide two numerical examples. In the first example we experimentally show that the LESFEM–IMEX Euler method, applied to the heat equation with a linear source term on a linearly growing sphere, exhibits optimal convergence rates in space and time. In the second example we consider the Thomas RDS on an exponentially growing Dupin ring cyclide thereby showing (i) the existence of invariant regions for the fully discretised model and (ii) the violation of this region in the absence of mass lumping. The present paper is structured as follows. First, in Sect. 2 we recall (i) the derivation of RDSs on evolving surfaces and (ii) some basic notions about invariant regions and we introduce the notion of dilation rate, which is crucial in our analysis. In Sect. 3 we introduce the LESFEM for the space discretisation of RDSs on evolving surfaces and we carry out a fully- discrete scheme using the IMEX Euler timestepping. Section 4 deals with the characterisation of the dilation rate in terms of the material velocity. In particular, we compute exactly the dilation rate for the class of isotropic growth laws. In Sect. 5 we prove a semi- and a fully- discrete maximum principle for the linear heat equation on evolving surfaces with a linear source term. We prove, in Sect. 6, sufficient conditions for the existence of invariant regions for RDSs of arbitrarily many equations on evolving surfaces at the semi- and fully-discrete levels. Section 7 presents some classes of invariant regions for the activator-depleted and the Thomas RD models on evolving surfaces, respectively. Numerical examples are presented in Sect. 8. Finally, in Sect. 9 we conclude and discuss our findings with an eye for future extensions of the present work. 2 Reaction–Diffusion Equations on an Evolving Surface 2.1 Preliminaries and Basic Results 2 3 For t ∈[0, T ],let Γ(t ) be a C orientable surface in R , represented as the zero-level set 1 2 3 3 of a signed distance function d ∈ C ([0, T ], C (R )), i.e. Γ(t ) ={x ∈ R | d(x, t ) = 0}, with ∇d(x, t ) = 0for t ∈[0, T ] and x ∈ Γ(t ). Hence, the outward unit normal vector field on Γ(t ) is given by ∇d(x, t ) n(x, t ) := , t ∈[0, T ], x ∈ Γ(t ), (1) ∇d(x, t ) where · denotes the Euclidean norm in R . Following [10, Section 5], we assume that 3 1 2 there exists a mapping G : Γ(0) ×[0, T]→ R , G ∈ C ([0, T ], C (Γ (0))), such that for all t ∈[0, T ], G(Γ (0), t ) = Γ(t ) and G(·, t ) is a diffeomorphism between Γ(0) and Γ(t ). The space-time surface G is defined by G := Γ(t ) ×{t }. The material velocity T T t ∈[0,T ] v : G → R of Γ(t ) is defined by ∂G v(G(x , t ), t ) = (x , t ), ∀x ∈ Γ(0), t ∈[0, T ]. (2) 0 0 0 ∂t 1 2 3 Vice-versa, if v ˜ ∈ C ([0, T ], C (R )) is an extension of v, i.e. v ˜ (G(x , t ), t ) = v(G(x , t ), t ) for x ∈ Γ(0) and t ∈[0, T ], the mapping G (and thus the time-dependent 0 0 surface Γ(t )) is recovered by solving, for each x ∈ Γ(0), the Cauchy problem 123 J Sci Comput ∂G (x , t ) = v ˜ (G(x , t ), t ), t ∈[0, T ], 0 0 (3) ∂t G(x , 0) = x . 0 0 For t ∈[0, T ] and δ> 0, let U (t ) be the open neighbourhood of Γ(t ) defined by U (t ) := {(x, t ) ∈ R ×[0, T]:|d(x, t )| <δ}. (4) We recall from [10] the following property. Lemma 1 (Fermi coordinates) For any t ∈[0, T ] there exists δ> 0 such that for any x ∈ Γ(t ) there exists a unique a(x) ∈ Γ(t ) that fulfils x = a(x) + d(x)n(a(x)), (5) Hence, for t ∈[0, T ], every point x ∈ U (t ) can be described by d(x) and a(x), which are called normal coordinates or Fermi coordinates of x. The function a(·, t ) : U (t ) → Γ(t ) is called the normal projection onto Γ(t ). For any sufficiently smooth g : G → R and g : G → R ,let ∇ g, Δ g and ∇ · g T T Γ(t ) Γ(t ) Γ(t ) denote the tangential gradient, the Laplace-Beltrami operator and the tangential divergence of g on Γ(t ), respectively (see [10] for the details). In the following, we will write Γ instead of Γ(t ) to simplify the notation. Furthermore, let ∂ g denote the material derivative of g defined by ∂g ˜ ∂ g := + v ·∇˜ g, (6) ∂t where ∇ is the standard gradient in R and g ˜ is any differentiable extension of g defined on a neighborhood of G . Definition (6)is intrinsic, i.e. it does not depend on the choice of the extension g ˜ (see [11] for further details). Let us recall some basic results from [10]. Lemma 2 (Integration by parts) If g : G → R is sufficiently smooth and t ∈[0, T ],it holds that g · μ = ∇ · g − (g · n)(∇ · n), t ∈[0, T ], (7) Γ Γ ∂Γ (t ) Γ(t ) Γ(t ) where μ : ∂Γ (t ) → R is the outward conormal unit vector on ∂Γ (t ),i.e.normalto ∂Γ (t ) and tangent to Γ(t ). Specifically, if g is tangent to Γ ,i.e. g · n = 0, it holds that g · μ = ∇ · g, t ∈[0, T ]. (8) ∂Γ (t ) Γ(t ) Lemma 3 (Transport formula) For sufficiently smooth g : G → R, it holds that g = ∂ g + g∇ · v , t ∈[0, T ]. (9) dt Γ(t ) Γ(t ) Lemma 4 (Green’s formula on surfaces) For sufficiently smooth f, g : G → R, it holds that ∇ f ·∇ g =− f Δ g + f ∇ g · μ, t ∈[0, T ]. (10) Γ Γ Γ Γ Γ(t ) Γ(t ) ∂Γ (t ) 123 J Sci Comput Remark 1 (Surfaces without boundary) Lemmas 2–4 hold on surfaces with or without bound- ary, i.e. ∂Γ (t ) =∅ or ∂Γ (t ) =∅, respectively. Specifically, if ∂Γ (t ) =∅, then the boundary integrals in (7), (8)and (10)vanish. 2.2 Derivation of the Reaction–Diffusion Model in Strong Form Suppose we are given r ∈ N species u : Γ(t ) → R, k = 1,..., r,and let q : Γ(t ) → R , k k k = 1,..., r, be their fluxes tangent to Γ(t ). We recall from [1] the derivation of a system of r equations for u := (u ,..., u ) that accounts for (i) the diffusion on the surface, (ii) the 1 r flux across the boundary (if non-empty) and (iii) the production rates f (u), k = 1,..., r, of the given species. To this end, let R(0) be a portion of Γ(0) and let R(t ) = G(R(0), t ) be the portion of Γ(t ) corresponding to the initial portion R(0). Notice that R(t ) is itself an evolving surface, hence formulae (7) through (10) still hold if Γ(t ) is replaced by R(t ).We consider a mass balance on R(t ) of the form u =− q · μ + f (u), k = 1,..., r, t ∈[0, T ]. (11) k k k dt R (t ) ∂R (t ) R (t ) Since the fluxes q , k = 1,..., r, are tangent to Γ(t ), we can apply the integration-by-parts formula (8) to the first term on the right hand side of (11). Then (11) becomes u =− ∇ · q + f (u), k = 1,..., r, t ∈[0, T ]. (12) k Γ k k dt R (t ) R (t ) R (t ) By applying the transport formula (9) to the left hand side of (12), we obtain ∂ u + u ∇ · v +∇ · q = f (u), k = 1,..., r, t ∈[0, T ], (13) k k Γ Γ k k R (t ) R (t ) where v is the material velocity defined in Sect. 2.1.Since R(t ) is an arbitrary portion, we conclude that ∂ u + u ∇ · v +∇ · q = f (u), k = 1,..., r, t ∈[0, T ]. (14) k k Γ Γ k k We assume q corresponds to a diffusive flux according to Fick’s law as follows: q =−d ∇ u , ∀k = 1,..., r, (15) k k Γ k where d , k = 1,..., r, are positive diffusivity constants. Inserting (15)into(14), we end up with the reaction–diffusion system of the form ∂ u + u ∇ · v = d Δ u + f (u), k = 1,..., r, t ∈[0, T ]. (16) k k Γ k Γ k k 2.3 Invariant Regions and Maximum Principle In this section we recall basic notions concerning invariant regions for systems of the form (16) and conjecture a sufficient condition under which system (16) possesses an invariant region. To this end, we give the following definitions. Definition 1 (Dilation rates)The minimum and maximum instantaneous dilation rates are defined by ∗ ∗ H (t ) := min ∇ · v(x, t ) and H (t ) := max ∇ · v(x, t ), t ∈[0, T ], (17) Γ Γ min max x∈Γ(t ) x∈Γ(t ) 123 J Sci Comput respectively. When the minimum and maximum instantaneous dilation rates coincide, we call ∗ ∗ ∗ H (t ) := H (t ) = H (t ) the instantaneous dilation rate.The minimum and maximum min max global dilation rates are defined by ∗ ∗ ∗ ∗ μ := min H (t ) and μ := max H (t ), (18) min min max max t ∈[0,T ] t ∈[0,T ] respectively. When the minimum and maximum global dilation rates coincide, we call μ := ∗ ∗ μ = μ the global dilation rate. max min Definition 2 (Invariant regions) For the system (16),aregion  in the phase-space R is said to be an invariant region if, whenever the initial condition u is in , the solution u stays in  as long as it exists and is unique. We focus our attention on regions  ⊂ R of hyper rectangular shape, that is to say of the form := [σ , σ ], (19) k=1 where σ ∈ R ∪{−∞} and σ ∈ R ∪{+∞} for all k = 1,..., r. For instance, if σ = 0and k k σ =+∞ for all k = 1,..., r,then  is the positive orthant in R , which means that the solution of the RDS stays positive at all times. Consider the (r − 1)-dimensional hyperfaces :=  ∩{u = σ },  :=  ∩{u = σ }, k = 1,..., r. k k k k k k For k = 1,..., r, we define the constants ∗ ∗ μ if σ ≥ 0, μ if σ ≥ 0, ∗ k ∗ min max k μ := μ := (20) ∗ ∗ μ if σ < 0, μ if σ < 0, max k min ∗ ∗ where μ and μ are defined in (18). min max Next, we conjecture a criterion under which a hyper-rectangle is invariant for system (16). This criterion holds true in the stationary cases (when μ = 0): (i) when Γ is a stationary monodimensional domain in R (see [38]), (ii) when Γ is a stationary k-dimensional domain in R , k ∈ N (see [6]) and (iii) when Γ is a stationary Riemannian manifold without boundary (see [39]). In the case of isotropically evolving flat domains, the invariance of the positive orthant was studied in [41]. To the best of the author’s knowledge, the case of evolving surfaces has not been studied at the continuous level. Hence, we introduce at the continuous level the following conjecture. Conjecture 1 Let  be a hyper-rectangle as in (19) in the phase space of (16) and let f be Lipschitz on .If ∗ r f (u)< μ σ , ∀u ∈  ∩ R , ∀k = 1,..., r, (21) k k k ∗ r f (u)>μ σ , ∀u ∈  ∩ R , ∀k = 1,..., r, (22) k k then  is an invariant region for (16). In particular, when σ =+∞ and σ =−∞,then r r ∩ R and  ∩ R are respectively empty, and so (21) and (22) are automatically fulfilled, respectively. ∗ ∗ Notice that on stationary surfaces, since μ = μ = 0, then conditions (21)–(22) reduce to the inward flux condition considered in [39, Chapter 4]. In Sect. 4 next, we will prove the discrete counterpart of Conjecture 1 obtained by discretizing the RDS in space with the lumped version of the ESFEM method [10] and IMEX Euler in time. 123 J Sci Comput We remark that, on stationary domains, some systems are known to possess an invariant region which do not meet the strict inequalities (21)–(22). For instance, for many mass-action laws, the positive orthant is invariant [4,16] even though the flux of f is tangent to this region, instead of strictly inward. For the scalar case k = 1, the notion of invariant region reduces to the well-known concept of maximum principle. Definition 3 (Maximum principle)For k = 1, the scalar equation (16) fulfils the maximum principle if, for any initial condition, i.e. u(·, 0), the solution u fulfils min 0, min u(y, 0) ≤ u(x, t ) ≤ max 0, max u(y, 0) ,(x, t ) ∈ G . (23) y∈Γ(0) y∈Γ(0) In particular, if the initial condition is nonnegative, u(y, 0) ≥ 0for y ∈ Γ(0), the solution fulfils 0 ≤ u(x, t ) ≤ max u(y, 0), (x, t ) ∈ G . (24) y∈Γ(0) Notice that the maximum principle corresponds to the fact that every (monodimensional) region of the form  =[σ , σ ], with σ ≤ 0, σ ≥ 0, is invariant. 2.4 Derivation of the Variational Formulation Following [1], we derive the variational formulation of system (16). To this end, for each t ∈[0, T ] we multiply equations (16) by the respective test functions ϕ ,...,ϕ ∈ 1 r 2 1 • • 2 −1 L ([0, T ]; H (Γ (t ))) with ∂ ϕ ,...,∂ ϕ ∈ L ([0, T ]; H (Γ (t ))) and integrate over 1 r Γ(t ): (ϕ ∂ u + ϕ u ∇ · v − ϕ f (u)) = d ϕ Δ u , (25) k k k k Γ k k k k Γ k Γ(t ) Γ(t ) 1 −1 for all k = 1,..., r. For the rigorous definition of the Sobolev spaces H and H on a manifold see [21], while for the Bochner spaces L ([0, T ]; B), with B being any Banach space, see [34]. By applying the Green formula (10) to the right hand side of (25) we obtain (ϕ ∂ u + ϕ u ∇ · v − ϕ f (u)) + d ∇ ϕ ·∇ u = d ϕ ∇ u · μ, k k k k Γ k k k Γ k Γ k k k Γ k Γ(t ) Γ(t ) ∂Γ (t ) (26) for all k = 1,..., r. We assume that either Γ(t ) has no boundary, i.e. ∂Γ (t ) =∅,or homogeneous Neumann boundary condition are enforced, i.e. ∇ u · μ = 0on ∂Γ (t ),so Γ k that the last term in (26) vanishes (see Remark 1). Furthermore, by observing that ∂ (ϕ u ) = k k • • ϕ ∂ u + u ∂ ϕ ,(26) becomes k k k k • • ∂ (ϕ u ) = u ∂ ϕ − ϕ u ∇ · v + ϕ f (u) − d ∇ ϕ ·∇ u , k k k k k k Γ k k k Γ k Γ k Γ(t ) Γ(t ) Γ(t ) (27) for all k = 1,..., r. By applying the transport property to the first term on the left hand side of (27), we have ϕ u − u ∂ ϕ = ϕ f (u) − d ∇ ϕ ·∇ u , (28) k k k k k k k Γ k Γ k dt Γ(t ) Γ(t ) Γ(t ) Γ(t ) 123 J Sci Comput for all k = 1,..., r. Therefore, the variational formulation seeks to find u ,..., u ∈ 1 r 2 1 • • 2 −1 L ([0, T ]; H (Γ (t ))) with ∂ u ,...,∂ u ∈ L ([0, T ]; H (Γ (t ))) such that, for each 1 r t ∈[0, T ], u ϕ − u ∂ ϕ + d ∇ u ·∇ ϕ = f (u)ϕ , (29) k k k k k Γ k Γ k k k dt Γ(t ) Γ(t ) Γ(t ) Γ(t ) 2 1 • • 2 −1 for all ϕ ,...,ϕ ∈ L ([0, T ]; H (Γ (t ))) with ∂ ϕ ,...,∂ ϕ ∈ L ([0, T ]; H (Γ (t ))). 1 r 1 r 3 Lumped Evolving Surface Finite Element Method Following the evolving surface finite element method (ESFEM) studied in [1] for the approx- imation of the variational problem (29), we present its lumped counterpart, the lumped evolving surface finite element method (LESFEM). 3.1 Surface Triangulation and Some Definitions Following [9], given h > 0, called meshsize, a triangulation Γ (t ) of the evolving surface Γ(t ) is defined by Γ (t ) = Z (t ), Z (t )∈Z (t ) where Z (t ) is a set of evolving triangles, with x (t ), i = 1,..., N ∈ N, being the overall h i evolving nodes, such that – The nodes evolve with the exact material velocity, i.e. x ˙ (t ) = v(x (t ), t ) for t ∈[0, T ] i i and i = 1,..., N ; –For all t ∈[0, T ] and for any two distinct triangles Z (t ) and Z (t ) in Z (t ),the 1 2 h intersection Z (t ) ∩ Z (t ) is either empty, or a node, or a complete edge; 1 2 –For all t ∈[0, T ] and Z (t ) ∈ Z (t ), Z (t ) ⊂ U (t ),where U is as defined in Lemma 1; h δ δ –Forall t ∈[0, T ], the normal projection a(·, t ) : U (t ) → Γ(t ) defined in Lemma 1 is a one-to-one mapping between Γ (t ) and Γ(t ), i.e. a(Γ (t ), t ) = Γ(t ). h h We assume that, for each t ∈[0, T ], Γ (t ) meets the Delaunay condition, defined as follows. Let e be an edge of Γ (t ) and let Z and Z be the faces of Γ (t ) sharing the edge e.Let α h 1 2 h 1 and α be the angles in Z and Z opposite to e, respectively. The triangulation Γ (t ) is said 2 1 2 h to meet the Delaunay condition if, for all t ∈[0, T ] and for any edge e of Γ (t ), it holds that α + α ≤ π, (30) 1 2 see for instance [16]. The space-time triangulated surface G is defined by h,T G := Γ (t ) ×{t }. h,T h t ∈[0,T ] Since Γ (t ) is piecewise planar, there exists a time-differentiable mapping G : Γ (0) × h h h [0, T]→ R such that, for all t ∈[0, T ], G (Γ (0), t ) = Γ (t ) and, for every facet Z ∈ Z , h h h h G (·, t ) is a diffeomorphism between Z (0) and Z (t ). For a fixed time t ∈[0, T ],let V (t ) be the space of piecewise linear functions on Γ (t ) h h defined by V (t ) := {ϕ ∈ C (Γ (t )) ϕ is linear affine for each Z ∈ Z }. (31) h h | h 123 J Sci Comput Let V be the space of time-dependent, spatially piecewise linear functions defined by V := {ϕ : G → R | ϕ(·, t ) ∈ V (t ) for each t ∈[0, T ]}. (32) h h,T h Given t ∈[0, T ] and a function η ∈ C (Γ (t )), its linear interpolant I η is the unique h h function in V (t ) such that I η(x (t )) = η(x (t )), i = 1,... N . h i i The discrete material derivative of a sufficiently smooth function U ∈ V is defined by ∂U ∂ U := + I (v) ·∇U,(x, t ) ∈ G , h h,T ∂t where v is the material velocity. For our purposes, we define the discrete counterpart of the dilation rates introduced in Definition 1. Definition 4 (Discrete dilation rates)The minimum and maximum discrete instantaneous dilation rates are defined by H (t ) := ess inf ∇ · I (v)(x, t ), min x∈Γ (t ) Γ h h h H (t ) := ess sup ∇ · I (v)(x, t ), t ∈[0, T ], (33) max Γ h x∈Γ (t ) h respectively. When the minimum and maximum discrete instantaneous dilation rates coincide, we call H (t ) := H (t ) = H (t ) the discrete instantaneous dilation rate.The minimum min max and maximum discrete global dilation rates are defined by μ := min H (t ), μ := max H (t ), (34) min min max max t ∈[0,T ] t ∈[0,T ] respectively. When the minimum and maximum discrete global dilation rates coincide, we call μ := μ = μ the discrete global dilation rate. min max For every i = 1,..., N,the i-th Lagrange basis function χ is the unique V function such i h that χ (x (t ), t ) = δ , t ∈[0, T ], i, j = 1,... N , (35) i j ij where δ is the usual Kronecker symbol. The components U ,..., U ∈ V of the spatially ij 1 r h discrete solution may be expressed in the Lagrange basis as U (x, t ) = ξ (t )χ (x, t ), (x, t ) ∈ G , k = 1,..., r. (36) k k,i i h,T i =1 3.2 Preliminary Results on Triangulated Surfaces We recall from [9] the following property of the basis functions. Lemma 5 (Transport property of the basis functions) The basis functions χ ,i = 1,..., N, defined in (35) fulfil ∂ χ = 0, i = 1,..., N . (37) Hence, for the functions U ,k = 1,..., r defined in (36) it holds that ∂ U (x, t ) = ξ (t )χ (x, t ), (x, t ) ∈ G , k = 1,..., r. (38) k k,i i h,T i =1 123 J Sci Comput We recall from [10, Lemma 5.6 and Remark 5.7] the following preliminary result. Lemma 6 (Leibniz formula on triangulated surfaces) For any time-differentiable U, V ∈ V , it holds that • • UV = ∂ UV + U ∂ V + UV ∇ · I (v). (39) Γ h dt Γ (t ) Γ (t ) Γ (t ) Γ (t ) h h h h 3.3 Lumped Evolving Surface Finite Element Method The lumped evolving surface finite element method (LESFEM), applied to the variational formulation (29), seeks to find U ,..., U ∈ V such that 1 r h I (U χ ) − U ∂ χ + d ∇ U ·∇ χ = I ( f (U)χ ), (40) h k i k i k Γ k Γ i h k i h h h dt Γ (t ) Γ (t ) Γ (t ) Γ (t ) h h h h for all k = 1,..., r and i = 1,..., N . Thanks to the transport property (37) of the basis functions, formulation (40) is equivalent to: find U ,..., U ∈ V such that 1 r h I (U χ ) + d ∇ U ·∇ χ = I ( f (U)χ ), (41) h k i k Γ k Γ i h k i h h dt Γ (t ) Γ (t ) Γ (t ) h h h for all k = 1,..., r and i = 1,..., N . The LESFEM method (41) differs from the evolving surface finite element method (ESFEM) in [1] due to the presence of the interpolant operator on the first and the last terms in (41). By expressing U ,..., U according to (36), the matrix 1 r form of (41)is (M ξ ) + d Aξ = Mf (ξ ,...,ξ ), k = 1,..., r, (42) k k k k 1 r dt where A and M are the (time-dependent) stiffness and lumped mass matrices defined by A (t ) = ∇ χ ·∇ χ , i, j = 1,..., N , ij Γ i Γ j h h Γ (t ) χ if i = j, Γ (t ) M (t ) = I (χ χ ) = i, j = 1,..., N , ij h i j 0if i = j, Γ (t ) for t ∈[0, T ], respectively. Notice that the lumped mass matrix M is diagonal. The Delaunay condition (30) holds if and only if A (t ) ≤ 0, i = j, t ∈[0, T ]. (43) ij The above fact, together with the structure of the lumped mass matrix implies that, for any s ≥ 0, −1 (M + sA) M (t ) ≥ 0, t ∈[0, T ]; (44) −1 (M + sA) M (t )1 = 1, t ∈[0, T ]. (45) See [16] for further details. 123 J Sci Comput 3.4 Time Discretisation We are now concerned with the time discretisation of the spatially discrete system (42)arising from the LESFEM. We discretise system (42) by means of the IMEX (IMplicit–EXplicit) Euler method, i.e by treating diffusion implicitly and reaction terms explicitly. To this end, let τ> 0 be a timestep, let t := nτ for all n = 0,..., N with N = ,let A and n T T n n n M be the stiffness and lumped mass matrices at time t , respectively, let (ξ ,...,ξ ) be 1 r n n n the coefficients of the numerical solution at time t ,and let f := f (ξ ,...,ξ ) for each n k k 1 0 0 k = 1,..., r.If (ξ ,...,ξ ) are the coefficients of the spatially discrete initial datum, the 1 r IMEX Euler time discretisation of (42)is n+1 n+1 n n M ξ − M ξ k k n+1 n+1 n n + d A ξ = M f , k = 1,..., r, n = 0,..., N , (46) k T k k or equivalently n+1 n+1 n+1 −1 n n n ξ = (M + τ d A ) M (ξ + τ f ), k = 1,..., r, n = 0,..., N . (47) k T k k 4 Characterisation of Surface Growth The purpose of this section is to characterise surface growth in terms of the material velocity v, with specific regard to isotropic growth. In fact, the lumped mass M (t ) and stiffness A(t ) matrices depend on v. In particular: 1. For an arbitrary triangulated surface that evolves with an arbitrary material velocity, we dM bound the time derivative of the lumped mass matrix in terms of the constants μ min dt and μ defined in (34), i.e. in terms of the divergence ∇ · I (v) of the discrete max Γ h material velocity. We will need this result to prove a sufficient condition for the existence of invariant regions for the semi- and fully-discrete schemes; 2. For an arbitrary smooth or triangulated surface that evolves with an arbitrary material velocity, we characterise the velocity flows ∇ ·v and ∇ · I (v) in terms of the mappings Γ Γ h G and G introduced in Sects. 2.1 and 3.1, respectively; 3. For an arbitrary smooth or triangulated surface that evolves isotropically in space, that is v(x, t ) = S(t )x,(x, t ) ∈ R ×[0, T ], (48) where S :[0, T]→ R is an arbitrary smooth function, we compute exactly ∇ · v and ∇ · I (v) in terms of v. This result will yield a fully practical criterion to detect the Γ h invariant regions of a given RDS in the case of isotropic surface evolution. 4.1 Bounding the Rate of Change of the Mass Matrix in Terms of the Dilation Rates dM In this section we bound the time derivative of M in terms of the discrete dilation rates dt dM defined in (34). To this end, we prove the following characterisation of . dt Lemma 7 (Transport formula for the lumped mass matrix M) The entries of the lumped mass matrix M fulfil the following property I (χ χ ) = I (χ χ )∇ · I (v), ∀i = 1,..., N . (49) h i j h i j Γ h dt Γ (t ) Γ (t ) h h 123 J Sci Comput Proof By choosing U = I (χ χ ) for any i, j = 1,..., N and V = 1 in the Leibniz formula h i j (39)wehave I (χ χ ) = ∂ I (χ χ ) + I (χ χ )∇ · I (v). (50) h i j h i j h i j Γ h dt Γ (t ) Γ (t ) Γ (t ) h h h Now, if i = j,then I (χ χ ) = χ , otherwise, if i = j, I (χ χ ) = 0. Then, from the h i j i h i j transport property (37)wehave ∂ I (χ χ ) = 0for all i, j = 1,..., N . Equation (50) thus h i j implies the desired result (49). In some proofs we will need the following corollary of Lemma 7. Corollary 1 (Consequence of the transport formula for the lumped mass matrix M) The dM diagonal matrix fulfils the estimates dt dm ii μ m (t ) ≤ (t ) ≤ μ m (t ), i = 1,..., N , t ∈[0, T ], (51) min ii max ii dt where μ and μ are defined in (34). min max 4.2 Characterising Velocity Flows in Terms of the Mappings G and G We wish to characterise the continuous and discrete velocity flows ∇ · v and ∇ · I (v) in Γ Γ h terms of the mappings G and G introduced in Sects. 2.1 and 3.1, respectively, for arbitrary smooth or triangulated surfaces that evolve under an arbitrary material velocity. To this end, let Γ(t ) be an arbitrary evolving smooth surface and let (A, X ) be any local parametrisation of Γ(0),where A ⊂ R is an open connected set and X : A → Γ(0) is a differentiable map such that its Jacobian J is full-rank on A.Let B be a measurable subset of A.For all t,the portion X (B) of Γ(0) evolves into G(X (B), t ) ⊂ Γ(t ). By choosing f (x, t ) = 1 for each (x, t ) ∈ G in the transport formula (9), we have 1 = ∇ · v. (52) dt G(X (B),t ) G(X (B),t ) Let e ˆ , i = 1, 2, 3, be the standard basis vectors in R .For (θ , t ) ∈ B ×[0, T ],let J (x, t ) ∈ 2,3 R be the (spatial) Jacobian of the function G(X (θ ), t ) and let ⎡ ⎤ e ˆ e ˆ e ˆ 1 2 3 ˜ ⎣ ⎦ J (θ , t ) := . (53) J (θ , t ) The surface integrals in (52) can be written as integrals on the planar domain B by using the parametrisation G(X (·), t ) : B → G(X (B), t ). Hence, (52) becomes ˜ ˜ det(J (θ , t ))dθ = ∇ · v(G(X (θ ), t )) det(J (θ , t ))dθ, (54) dt B B where · denotes the Euclidean norm in R , or equivalently, ˜ ˜ det(J (θ , t ))dθ = ∇ · v(G(X (θ ), t )) det(J (θ , t ))dθ. (55) dt B B Since (55) holds for any measurable subset B of A, then it holds that ˜ ˜ det(J (θ , t ))=∇ · v(G(X (θ ), t )) det(J (θ , t )). (56) dt 123 J Sci Comput By applying the chain rule, (56) is equivalent to ∇ · v(G(X (θ ), t )) = ln  det(J (θ , t )). (57) dt Given any triangulated surface Γ (t ) (which evolves under the discrete velocity field I (v)), h h we notice that every facet Z (t ) ∈ Z (t ) is smooth and thus parametrisable. For any facet of the initial triangulated surface, Z (0) ∈ Z (0),let (A, X ) be a parametrisation of Z (0),as 2,3 described above. For (θ , t ) ∈ A ×[0, T ],let J (θ , t ) ∈ R be the Jacobian of the function G (X (θ ), t ) and let ⎡ ⎤ e ˆ e ˆ e ˆ 1 2 3 ˜ ⎣ ⎦ J (θ , t ) := . (58) J (θ , t ) A similar reason as above yields the following discrete counterpart of (57). ∇ · I v(G (X (θ ), t )) = ln  det(J (θ , t )). (59) Γ h h h dt Relations (57)and (59) are useful in that they (i) express the velocity flow without tangential derivatives and (ii) can be computed exactly when G and G are known explicitly, e.g. for the isotropic growth as discussed in the next subsection. 4.3 Computing the Dilation Rates for the Isotropic Growth In this section we compute the dilation rates defined in (18)and (34), respectively, for arbitrary surfaces that evolve with the material velocity (48). The velocity field (48) corresponds to the specific case of isotropic evolution, see for instance [8,29] for the case of evolving planar domains and [1,32] for the general case. In particular, for suitable choices of the function S(t ), the growth law (48) admits some specific cases such as uniform, exponential, logistic and periodic growth, see [1,8,29,32]. From (3), it is easy to show that, with the velocity field (48), each x ∈ Γ evolves to the point 0 0 G(x , t ) = exp S(τ )dτ x , t ∈[0, T ], (60) 0 0 therefore the evolution induced by an isotropic growth is a time-dependent dilation of the initial surface. The function φ(t ) := exp S(τ )dτ , t ∈[0, T ], (61) that appears in (60) is known as the growth function, see for instance [29]. Remark 2 (Properties of isotropic growth) Isotropic growth preserves the angles of triangu- lated surfaces. This has two consequences: 1. if Γ (0) meets the Delaunay condition, then Γ (t ) retains the Delaunay condition for all h h t ∈[0, T ]; 2. if A(0) and M (0) are the stiffness and the mass matrices at t = 0, then A(t ) = A(0), M (t ) = φ (t )M (0), t ∈[0, T ]. (62) Hence, in implementations, A(t )and M (t ) need not be computed at each time step. 123 J Sci Comput ∗ ∗ In the following result we compute the dilation rates μ , μ , μ and μ on an min max max min arbitrary smooth or triangulated surface that evolves with the material velocity (48), in terms of S(t ). Theorem 1 (Velocity flow on isotropically growing smooth or triangulated surfaces) Let Γ(t ) be a smooth surface that evolves with the velocity field (48) and let Γ (t ) be the corresponding triangulated surface. Then, the instantaneous dilation rates satisfy H (t ) = 2S(t ) = H (t ), t ∈[0, T ]. (63) Hence, it follows that ∗ ∗ μ = μ = 2min S(t ), and μ = μ = 2max S(t ). (64) min max min max t ∈[0,T ] t ∈[0,T ] Proof Let Γ (t ) be an evolving smooth surface and let (A, X ), J (θ , t ) and J (θ , t ) as defined in Sect. 4.2.From(60)and (61), J (θ , t ) fulfils J (θ , t ) = φ(t )J (θ ), (θ , t ) ∈ A ×[0, T ], (65) 2,3 where J : A → R is the Jacobian of X. It follows that ˜ ˜ det J (θ , t )= φ (t ) det J (θ , 0),(θ, t ) ∈ A ×[0, T ], (66) which implies that ˜ ˜ ln  det J (θ , t )−ln  det J (θ , 0)=2ln φ(t ) =2 S(τ )dτ, (θ, t ) ∈ A×[0, T ]. (67) By differentiating (67), we have ln  det J (θ , t )= 2S(t ), (θ , t ) ∈ A ×[0, T ]. (68) dt By combining (57), (68) and dropping the parameterised coordinates θ,wehave ∇ · v(x, t ) = 2S(t ), (x, t ) ∈ G , (69) Γ T which proves the first equality in (63). Analogously, we prove the second equality in (63)by using (59). This completes the proof. For the uniform, exponential, logistic and periodic growths, the functions φ(t ), S(t ) and the ∗ ∗ dilation rates μ , μ , μ and μ are detailed in Table 1, see also [29,Table 1] for min max min max the case of evolving planar domains, while the corresponding plots are depicted in Fig. 1. 5 Linear Heat Equation and Discrete Maximum Principles We consider, for k = 1, the specific case of linear heat equation on an evolving surface Γ(t ): ∂ u + u∇ · v = dΔ u − βu,β ∈ R, (70) Γ Γ and we prove the semi- and fully-discrete maximum principles for the case when μ + β ≥ min 0. Equation (70) is a special case of the general system (16) that we are interested in. However, we start with this specific case as (i) it provides more insights on the effect of growth on stability, (ii) we are able to prove a better timestep stability condition and (iii) to make the reader familiar with the demonstrative techniques. The following result, that we proved in [16, Theorem 2.1] for the special case of stationary surfaces, addresses the maximum principle at the semi-discrete level. 123 J Sci Comput Table 1 Particular types of growth with their respective growth functions φ(t ), S(t ) functions and constants ∗ ∗ μ , μ , μ ,and μ max min max min Type of growth Growth function φ(t ) S(t)μ = μ μ = μ min max min max r 2r Linear rt + 1 2r rt + 1 rT + 1 Exponential exp(rt ) r 2r 2r K exp(Krt ) rK (K − 1) 2rK (K − 1) Logistic 2r (K − 1) K − 1 + exp(Krt ) K − 1 + exp(Krt ) K − 1 + exp(KrT ) √ √ r sin(rt ) 2r 3 2r 3 Periodic 2 − cos(rt ) − 2 − cos(rt ) 3 3 The constant r > 0 is the growth rate. For the logistic growth, K > 1isthe carrying capacity, i.e. the square root of the ratio between the asymptotical and the initial area of the surface (see [29]) φ (t) 10 1.5 10 1 r=0.5 0 5 10 0 5 10 0 5 10 0 5 10 r=1 r=1.5 1.5 1.5 1.5 1 1 1 S(t) 0.5 0.5 0.5 −1 0 0 0 0 5 10 0 5 10 0 5 10 0 5 10 t t t t Fig. 1 Plots of the growth functions φ(t ) (top row) and the corresponding S(t ) functions (bottom row) listed in Table 1 for K = 2and r = 0.5, 1, 1.5. From left to right: linear, exponential, logistic and periodic growth profiles Theorem 2 (Semi-discrete maximum principle for the linear heat equation (70)) If the veloc- ity field v fulfils μ + β ≥ 0, (71) min with μ as defined in (34), and the triangulation Γ (t ) meets the Delaunay condition for min h all t ∈[0, T ], then the LESFEM solution of (70) fulfils the following discrete maximum principle min 0, min ξ (0) ≤ ξ (t ) ≤ max 0, max ξ (0) , i = 1,..., N , t > 0. j i j j =1,...,N j =1,...,N (72) Proof From (42), the LESFEM spatial discretisation of (70)is (Mξ) + dAξ =−β Mξ. (73) dt By applying the chain rule, (73) becomes dM M ξ + ξ + dAξ =−β Mξ. (74) dt 123 J Sci Comput −1 By multiplying (74) on the left by M ,wehave dM −1 −1 ξ =−dM Aξ − M ξ − βξ. (75) dt −1 −1 dM All we have to prove is that the ODE (75) is dissipative, i.e. −dM A|ξ|− M |ξ|− dt β|ξ|≤ 0.For every t ∈[0, T ], M is diagonal with positive diagonal entries and A fulfils −1 (43) from the Delaunay condition. Then, it follows that −dM A|ξ|≤ 0. Hence, it suffices −1 dM to prove that −(M + β I )|ξ|≤ 0, that is true provided dt dM −1 M + β I ≥ 0. (76) dt By using (51), condition (76)istrueif (μ + β)I ≥ 0, (77) min which holds true from assumption (71). This completes the proof. The next theorem shows the same result for the LESFEM–IMEX Euler full-discretisation of (70), under a timestep restriction. This result holds true for the special case of stationary surfaces (see [16, Theorem 2.2]). Theorem 3 (Fully-discrete maximum principle for the linear heat equation (70)) If the veloc- ity field v fulfils μ + β ≥ 0, ∀t ≥ 0, (78) min with μ as defined in (34), and the triangulation Γ meets the Delaunay condition for min h all t > 0, then the LESFEM–IMEX Euler solution of (70) fulfils the following minimum- maximum principle 0 n 0 min 0, min ξ ≤ ξ ≤ max 0, max ξ , i = 1,..., N , n = 0,..., N , j j j =1,...,N j =1,...,N (79) if the timestep satisfies τβ ≤ 1. (80) In particular, there is no timestep restriction if β ≤ 0. Proof The full-discretisation (47) of the heat equation (70) can be written as n+1 n+1 n+1 −1 n+1 n+1 −1 n n ξ = (M + τdA ) M (M ) M (1 − τβ)ξ , n = 0,..., N . (81) From (44)–(45)wehave n+1 n+1 −1 n+1 (M + τdA ) M ≥ 0, n = 0,..., N , (82) n+1 n+1 −1 n+1 (M + τdA ) M 1 = 1, n = 0,..., N . (83) Then scheme (81) fulfils the discrete maximum principle if n+1 −1 n n n (M ) M (1 − τβ)ξ ≥ min{0,ξ }, n = 0,..., N ; (84) n+1 −1 n n n (M ) M (1 − τβ)ξ ≤ max{0,ξ }, n = 0,..., N . (85) 123 J Sci Comput n+1 −1 n Since (M ) M is diagonal with strictly positive diagonal entries, conditions (84)–(85) aretrueprovided 1 − τβ ≥ 0; (86) n −1 n+1 (1 − τβ)I ≤ (M ) M , n = 0,..., N . (87) Condition (86) is true under assumption (80). In order to prove (87), we need to estimate n+1 n M as a function of M . To this end, by applying Gronwall’s lemma to the first inequality in (51), we have n+1 n τμ min M ≥ M e , n = 0,..., N . (88) By using (88), condition (87)istrueif τμ min 1 − τβ ≤ e . (89) τμ min Let us now define f (τ ) := 1−τβ and g(τ ) = e . These functions fulfil f (0) = g(0) = 1, f is linear and g is non-concave for all μ ∈ R.Then min –if f (0)> g (0), then condition (89) is not fulfilled for any sufficiently small τ . –if f (0) ≤ g (0), then condition (89) is fulfilled for every τ> 0. Now, condition f (0) ≤ g (0) means −β ≤ μ , which is true from assumption (78). This min completes the proof. Remark 3 (Interplay between material velocity and source term) Relation (71) implies that – domain growth (μ > 0) can enable the discrete maximum principle even for β< 0; min – local domain contraction (μ < 0) can prevent the discrete maximum principle even min for β ≥ 0. This interplay is justified by observing that domain evolution implies a dilution effect, explained as follows. By choosing ϕ = 1 in the variational formulation (29) with k = 1 and f (u) =−βu, we obtain u =−β u, t ∈[0, T ]. (90) dt Γ(t ) Γ(t ) If |Γ(t )| denotes the surface area of Γ(t ) and u(t ):= u denotes the mean value |Γ(t )| Γ(t ) of u,(90) becomes (|Γ(t )|u(t )) =−β|Γ(t )|u(t ), t ∈[0, T ]. (91) dt By solving (91)for u(t ),weobtain dt β|Γ(t )|+ |Γ(t )| dt u(t )=− u(t ), t ∈[0, T ]. (92) dt |Γ(t )| By choosing g = 1 in the transport formula (9), we have |Γ(t )|= ∇ · v ≥|Γ(t )|μ , t ∈[0, T ]. (93) min dt Γ(t ) By combining (92)and (93)wehave u(t )≤−(β + μ )u(t ), t ∈[0, T ]. (94) min dt 123 J Sci Comput From (94), the dilution effect arising from surface growth can be interpreted as the dampening or uplifting effect of μ on u(t ). The estimate (94) implies that u(t ) is non-increasing min if β + μ ≥ 0, (95) min which is the continuous counterpart of (71). We conclude that condition (71) is consistent with the interpretation of surface growth in terms of dilution effect. Remark 4 (Interplay between timestep restriction and source term) Relation (80) implies that the timestep restriction needed for guaranteeing the discrete maximum principle is indepen- dent of the material velocity and it only depends on the stiffness parameter β of the source term. In particular, when the source term is nonnegative (i.e. when β ≤ 0), the LESFEM– IMEX Euler fully-discrete scheme unconditionally fulfils the discrete maximum principle. 6 Reaction–Diffusion Systems and Invariant Regions In this section we prove, for the semi- and full-discretisations of RDSs of the form (16), a criterion to test if a hyper-rectangle in the phase-space is invariant. In the case k = 1 of scalar equations, the notion of invariant region collapses to that of minimum-maximum principle, considered in the previous section for the special case of the linear heat equation. We assume that the Delaunay regularity of the mesh is preserved under evolution. For k = 1,..., r,we define the constants μ if σ ≥ 0, μ if σ ≥ 0, min k max μ := μ := (96) μ if σ < 0, μ if σ < 0, max k min where μ and μ are the dilation rates defined in (34). In the following theorem we min max prove that, under similar assumptions of Conjecture 1,  is an invariant region for the solution obtained from the semi-discrete scheme (42). Hence, the following theorem extends [16, Theorem 3.3] to the case of evolving surfaces. Theorem 4 (Invariant rectangles for (42)) Let  be a hyper-rectangle as in (19) in the phase space of (42),let f be Lipschitz on . If the triangulation Γ (t ) satisfies the Delaunay condition for all t ≥ 0 and f (U)< μ σ , ∀U ∈  ∩ R , ∀k = 1,..., r, (97) k k k k f (U)>μ σ , ∀U ∈  ∩ R , ∀k = 1,..., r, (98) k k then  is an invariant region for (42). Proof The semi-discrete method (42) can be written, after applying the chain rule to the term d −1 (M ξ ) and multiplying on the left by M as dt dM −1 −1 ˙ ¯ ¯ ξ =−d M Aξ + f (ξ ,...,ξ ) − M ξ , k = 1,..., r. (99) k k k k 1 r k dt Since M is diagonal with positive diagonal entries and A ≤ 0for i = j from the assumption ij of Delaunay regularity, proceeding as in the proof of [16, Theorem 3.3], it suffices to verify that, for all (U ,..., U ) ∈ , k = 1,..., r,and i = 1,..., N , 1 r 123 J Sci Comput dm ii −1 f (U ,..., σ ,..., U ) − m σ < 0, (100) k 1 k r k ii dt dm ii −1 f (U ,...,σ ,..., U ) − m σ > 0, (101) k 1 r k ii k dt where σ and σ are as in (19). Using relation (51), conditions (100)–(101) hold if, for all k = 1,..., r, f (U ,..., σ ,..., U ) − μ σ < 0, (102) k 1 k r k f (U ,...,σ ,..., U ) − μ σ > 0, (103) k 1 r k k with μ and μ as in (96), that is true from assumptions (97)–(98). This completes the proof. The following theorem provides a sufficient condition for regions to be invariant for the LESFEM–IMEX Euler scheme (47) and extends [16, Theorem 3.4]. In contrast to the semi- discrete case, we relax the strict inequalities (21)–(22) with conditions (105)–(106), in which we use the perturbed dilation rates μ ˜ and μ given by k k τμ μ if σ ≥ 0, max ⎪ min k ⎨ e − 1 if σ > 0, μ ˜ := μ := (104) τμ k max k τ e − 1 ˜ μ if σ ≤ 0, if σ < 0, k min respectively, and μ and μ are defined in (34). Observe that μ ˜ → μ and μ → μ min max k k k as τ → 0. Theorem 5 (Invariant rectangles for (47)) Let  be a hyper-rectangle as in (19) in the phase space of (42),let f be Lipschitz on . If the triangulation Γ (t ) meets the Delaunay condition for all t ≥ 0 and f (U) ≤ σ μ ˜ , ∀U ∈  ∩ R , ∀k = 1,..., r, (105) k k k k f (U) ≥ σ μ , ∀U ∈  ∩ R , ∀k = 1,..., r, (106) k k k k then  is an invariant region for (47) if the timestep τ fulfils τ L ≤ 1, k = 1,..., r, (107) where, for all k = 1,..., r, L is the Lipschitz constant of f . k k Proof The fully-discrete scheme (47) can be written as n+1 n+1 n+1 −1 n+1 n+1 −1 n n n ξ = (M + τdA ) M (M ) M (ξ + τ f ), (108) k k n ∈ N ∪{0}, k = 1,..., r. Since the mesh meets the Delaunay assumption at all times, the matrix properties (82)–(83) hold. Then, it suffices to prove that n+1 −1 n n n σ 1 ≤ (M ) M (ξ + τ f ) ≤ σ 1, k = 1,..., r, (109) k k k where 1 is the column vector of ones. We will prove the two inequalities in (109) in turn. From (51), the inequality on the right side of (109) holds true if n n τ μ ξ + τ f ≤ σ e 1, k = 1,..., r, (110) k k 123 J Sci Comput with μ as defined in (96). Suppose σ ≥ 0. From assumption (105) we can estimate f as follows n n f ≤ σ μ + L (σ 1 − ξ ), k = 1,..., r. (111) k k k k k k From (111), condition (110) holds true provided n τ μ ξ (1 − τ L ) + τ μ σ + τ L σ 1 ≤ σ e 1, k = 1,..., r. (112) k k k k k k k From assumption (107), since ξ ≤ σ 1,then(112) holds true if τ μ ˜ k σ (1 − τ L ) + τ μ σ + τ L σ ≤ σ e , k = 1,..., r, (113) k k k k k k k that is to say τ μ 1 + τ μ ≤ e , k = 1,..., r, (114) which holds true for each τ ∈ R. Suppose, instead, σ < 0. From assumption (105) we can estimate f as follows τ μ e − 1 n n f ≤ σ + L (σ 1 − ξ ), k = 1,..., r. (115) k k k k k From (115), condition (110) holds true provided n τ μ τ μ k k ξ (1 − τ L ) + σ (e − 1 + τ L )1 ≤ σ e 1, k = 1,..., r. (116) k k k k From assumption (107), since ξ ≤ σ 1,then(116) holds true if τ μ τ μ k k σ (1 − τ L )1 + σ (e − 1 + τ L )1 ≤ σ e 1, k = 1,..., r. (117) k k k k k As (117) always holds with the equality, we conclude that the second inequality in (109)is true under assumptions (105)and (107). Similarly, the inequality on the left side of (109) holds under assumptions (106)and (107). This completes the proof. Remark 5 (Sharper timestep restriction) In the specific case of the linear heat equation (70), estimate (80) in Theorem 3 is sharper than estimate (107) in Theorem 5. In fact, since the Lipschitz constant of the source term is L =|β|, the timestep restriction (107) is fulfilled for τ |β|≤ 1, that is more restrictive than condition (80). 7 Velocity-Induced Invariant Regions for RD Models Now, we consider two different RDSs that are well-known in the literature and prove, at the discrete level, the existence of discrete invariant hyper-rectangles for these RDSs, depending on the global discrete dilation rates μ and μ defined in (18). The results in this section min max are confined to the spatially discrete level, but from Conjecture 1, we claim that the same results holds at the continuous level. In the special case of stationary surfaces (i.e. when μ = μ = 0), we obtain invariant hyper-rectangles that have been studied in the min max literature (see [4,6,19]). It is worth remarking that, even though we consider two RD models for illustrative purposes, the following analysis can be easily extended to other types of RDSs. To mention two examples, for the well-known Hodgkin-Huxley model and for the DIB model for electrodeposition considered in [24,25], the invariant region study is carried out in [14]. 123 J Sci Comput 7.1 RDS with Activator-Depleted Kinetics Let us consider an RDS with the well-known non-dimensional activator-depleted kinetics, also known as Schnakenberg or Brusselator kinetics (see for instance [1,37]), on evolving surfaces • 2 ∂ u + u ∇ · v − Δ u = f (u , u ) := γ(a − u + u u ), 1 1 Γ Γ 1 1 1 2 1 2 (118) • 2 ∂ u + u ∇ · v − dΔ u = f (u , u ) := γ(b − u u ), 2 2 Γ Γ 2 2 1 2 2 where a, b and γ are positive parameters and d is a positive diffusion rate. The model describes a system of two interacting chemicals, in which u ≥ 0and u ≥ 0 are the respective 1 2 concentrations. For this reason, we focus our attention on invariant regions contained in the positive ortant. In the following theorem we prove that: (i) the positive orthant is invariant for (118) regardless of μ and μ . At the continuous level, the result holds in the specific min max case of stationary planar domains, see [4]. (ii) when μ > 0, the model possesses invariant min stripes (depending on μ ) in the positive orthant. min Theorem 6 (Velocity-induced invariant regions for the activator-depleted model (118)) For the LESFEM spatial discretisation of (118), the following statements hold: 1. For any value of the constants μ , and μ defined in (34), the positive orthant min max + 2 := [0, +∞[ is invariant. 2. If μ > 0 and σ is a constant such that min 2 γ b σ ≥ , (119) min then the stripe  =[0, +∞[×[0, σ ] is invariant. Proof In order to prove Statements (1) and (2) we have to verify conditions (21)–(22). For Statement (1), we observe that + + –  := {0}×[0, σ ]⊂  := {0}×[0, +∞[ and, for (u , u ) ∈  ,wehave 2 1 2 1 1 f (u , u ) = f (u , u ) = γ a > 0; 1 2 1 1 2 + + –  := [0, σ ]×{0}⊂  := [0, +∞[×{0} and, for (u , u ) ∈  ,wehave 1 1 2 2 2 f (u , u ) = f (u , u ) = γ b > 0. 1 2 2 1 2 This proves Statement (1). For Statement (2), let μ > 0 and we assume for the moment min that the strict inequality holds in (119). Then the set  := [0, +∞[×{σ } is contained in 1 2 the region γ a − (γ + μ )u min (u,v) ∈ R |u > 0,v > , γ u in which f (u , u ) := f (u , u ) − μ u < 0. This proves Statement (2) when the strict 1 1 2 1 1 2 min 1 inequality holds in (119). Otherwise, observe that =[0, +∞[×[0, σ ]= [0, +∞[×[0, σ + ε], (120) 2 2 ε>0 i.e.  is the intersection of invariant regions and is thus invariant. This completes the proof of Statement (2). 123 J Sci Comput 7.2 RDS with Thomas kinetics Let us consider an RDS with the non-dimensional Thomas kinetics (see for instance [31,p. 78]), on evolving surfaces u u • 1 2 ∂ u + u ∇ · v − Δ u = f (u , u ) := γ a − u − ρ , 1 1 Γ Γ 1 1 1 2 1 2 1+u +Ku (121) • u u ⎪ 1 2 ⎩ ∂ u + u ∇ · v − dΔ u = f (u , u ) := γ α(b − u ) − ρ , 2 2 Γ Γ 2 2 1 2 2 1+u +Ku where α, a, b, γ , K and ρ are positive constants and d is a positive diffusion rate. The model describes a system of two interacting chemicals, in which u ≥ 0and u ≥ 0are 1 2 the respective concentrations. For this reason, we focus our attention on invariant regions contained in the positive orthant. Theorem 7 (Velocity-induced invariant regions for the Thomas model (121)) For the LESFEM spatial discretisation of (121), the following statements hold: 1. For any value of the constants μ , and μ defined in (34), the positive orthant min max + 2 := [0, +∞[ is invariant. 2. If μ > −γ min(1,α) and σ and σ are two constants such that min 1 2 γ a γαb σ ≥ , σ ≥ , (122) 1 2 γ + μ γα + μ min min then the region  =[0, σ ]×[0, σ ] is invariant. 1 2 Proof To prove Statements (1) and (2), we have to verify conditions (21)–(22). For Statement (1), observe that –for (u , u ) ∈  := {0}×[0, σ ],wehave f (u , u ) = f (u , u ) = a > 0; 1 2 2 1 2 1 1 2 –for (u , u ) ∈  := [0, σ ]×{0},wehave f (u , u ) = f (u , u ) = αb > 0. 1 2 1 1 2 2 1 2 This proves Statement (1). For Statement (2), let μ > −γ min(1,α) and we assume for min the moment that the strict inequalities hold in (122). Then, observe that • the set  := {σ }×[0, σ ] is contained in the region 1 1 2 1 + u + Ku (u,v) ∈ R |u > 0,v > (γ a − (μ + γ)u) , min γρu in which f (u , u ) := f (u , u ) − μ u < 0; 1 2 1 1 2 min 1 • the set  := [0, σ ]×{σ } is contained in the region 2 1 2 γαb(1 + u + Ku ) (u,v) ∈ R |u > 0,v > , γρu + (γ α + μ )(1 + u + Ku ) min in which f (u , u ) := f (u , u ) − μ u < 0. 1 2 2 1 2 min 2 This proves Statement (2) when the strict inequalities hold in (122). Otherwise, we have that =[0, σ ]×[0, σ ]= [0, σ + ε]×[0, σ + ε], (123) 1 2 1 2 ε>0 i.e.  is the intersection of invariant regions and is thus invariant. 123 J Sci Comput 8 Numerical Examples The purpose of this section is to provide two numerical examples in which we (i) estimate the experimental order of convergence of the LESFEM and (ii) experimentally show the ability of Theorem 5 to find invariant regions of RDSs on evolving surfaces at the discrete level. 8.1 Numerical Example 1: Linear Heat Equation on an Evolving Sphere In this example, we wish to estimate the experimental order of convergence of the LESFEM. As a test problem, we consider the linear heat equation given by ∂ u + u∇ · v − Δ u = u, (124) Γ Γ We choose T = 1 to be the final time. The initial domain Γ(0) is the unit sphere S ,that evolves under the velocity field v(x, t ) := ,(x, t ) ∈ R ×[0, 1], (125) t + 1 and undergoes linear growth for r = 1, see Table 1. In particular, the domain Γ(t ) at time t ∈[0, T ] is a sphere whose radius is given by the growth function φ(t ) = t + 1and the 2r minimum dilation rate fulfils μ = = 1(seeTable 1). min rT +1 In order to determine the experimental order of convergence, we consider the analytical solution to (124)given by xyz t u(x , y, z, t ) = exp t − 2log(t + 1) − ,(x , y, z, t ) ∈ R ×[0, 1]. (t + 1) t + 1 (126) See “Appendix” for the derivation of the solution (126). The constants β =−1and μ = 1fulfill(78)and (80) for each τ> 0. Hence, the min LESFEM–IMEX Euler solution to (124) fulfils a discrete maximum principle unconditionally on τ . In order to appreciate the quadratic convergence in space we solve the problem on a sequence of eight Delaunay meshes Γ , i = 1,..., 8, whose mesh sizes at t = 0fulfil h (0) i −1 h (0) = 0.4013 and h (0) ≈ , i = 2,..., 8. The corresponding timesteps τ fulfil 1 i i τ = 4e-2 and τ = τ , i = 2,..., 8. In Fig. 2 we show a sequence of snap shots 1 i i −1 i −1 of the evolution of the numerical solution on the finest mesh Γ . The experimental order of ∞ 2 convergence is computed by measuring the error, in L ([0, T ], L (Γ (t ))) norm, between the numerical solution U and the piecewise linear interpolant I (u) of the exact solution. The result is shown in Fig. 3: the convergence is experimentally optimal in that it is quadratic in the meshsize and linear in the timestep. 8.2 Numerical Example 2: RDS with Thomas Kinetics and Invariant Regions In this example, we show that the LESFEM–IMEX Euler preserves the invariant regions of RDSs on evolving surfaces. Let us consider the following RDS with Thomas kinetics u u • 1 2 ⎨ ∂ u + u ∇ · v − d Δ u = f (u , u ) := γ a − u − ρ ; 1 1 Γ 1 Γ 1 1 1 2 1 1+u +Ku (127) • u u ⎪ 1 2 ⎩ ∂ u + u ∇ · v − Δ u = f (u , u ) := γ α(b − u ) − ρ , 2 2 Γ Γ 2 2 1 2 2 1+u +Ku 123 J Sci Comput Fig. 2 Numerical Example 1 on the linear heat equation (124): numerical solution at different times obtained on the finest mesh Γ with N = 16962 gridpoints, initial meshsize h = 3.542e−2 and timestep τ = 3.116e−4. 8 8 Plotted values range from −0.1924 (blue) to 0.1924 (red) (Color figure online) −2 2.2 EOC −3 E ∞ 2 L (L ) 1.8 −4 1.6 −1 −1 ∞ 2 Fig. 3 Numerical Example 1 on the linear heat equation (124): Error in L (0, T , L (Γ (t ))) norm (left) and experimental rate of convergence (right). The quadratic convergence in space is optimal considered in Sect. 7.2, with a = 150, b = 100, ρ = 13, K = 0.05 as in [30]and γ = 1for illustrative purposes. With these parameters, system (127) admits, in the absence of domain growth, the homogeneous steady state P ≈ (37.7382, 25.1588) calculated using Newton- Raphson method on a stationary domain. Notice that the diffusion coefficient for u has been normalised to 1 for convenience, and we choose d = 0.01. The initial domain Γ(0) is the Dupin ring cyclide considered in [18, Appendix B], rescaled for convenience, given by ⎧ ⎫ ⎨ 2 ⎬ 261 39 3249 3 2 2 2 2 D := (x , y, z) ∈ R : 9(x + y + z ) + − 4 6x − − y = 0 , ⎩ 100 10 25 ⎭ an orientable surface without boundary that is topologically equivalent to a torus. The surface rK (K −1) evolves under the velocity field (48) with S(t ) = > 0, t ∈[0, T ], with r = 0.2, K −1+exp(Krt ) K = 3 and undergoes a logistic growth, see Table 1. As final time we choose T = 100. From 123 J Sci Comput Fig. 4 Numerical Example 2 on the Thomas RDS (127): invariant region for the LESFEM–IMEX Euler full discretisation of (127) under the timestep restriction τ ≤ 1.676e−3. This region is obtained by considering the nullclines of the modified kinetics defined in (129)–(130) Theorem 1, it follows that the corresponding dilation rates fulfil ∗ ∗ μ = μ ≈ 0, and μ = μ = 0.8. (128) min max min max From Theorem (7), model (127) possesses arbitrarily large bounded invariant regions con- tained in the positive orthant. Let =[σ , σ ]×[σ , σ ]≈[0.3366, 126.4194]×[13.2938, 45.8182] 1 2 1 2 be the smallest bounded region that (i) contains P and (ii) meets the modified inward flux conditions (105)–(106) with τ =¯ τ := 2e-3 and thus (105)–(106) hold for all 0 <τ ≤¯ τ.As illustrated in Fig. 4,  is obtained by considering the nullclines of the modified kinetics f (u , u ) := f (u , u ) − μ u , k = 1, 2; (129) 1 2 k 1 2 min k τμ ¯ max e − 1 f (u , u ) := f (u , u ) − u , k = 1, 2. (130) 1 2 k 1 2 k τ¯ It is easy to see that, on a region of the form  =[σ , σ ]×[σ , σ ], the Lipschitz 1 2 1 2 constants L and L of the kinetics f and f of (127)fulfil 1 2 1 2 " " 2 2 2 2 ˆ ˆ L ≤ L := γ (1 + ρσ ) + ρ , and L ≤ L := γ (ρσ ) + (α + ρ) . 1 1 2 2 2 2 Hence, the timestep restriction (107) becomes τ ≤ ≈ 1.676e-3. It follows that ˆ ˆ max(L ,L ) 1 2 is invariant for the LESFEM–IMEX Euler full discretisation under the timestep restriction 0 <τ ≤ min τ, ¯ ≈ 1.676e-3. We choose τ = 1e-3. The region  is smaller ˆ ˆ max(L ,L ) 1 2 than the invariant region provided in Theorem 7, which has a simple analytical expression but is not optimal. The following initial condition u (x , y, z) = σ + (σ − σ )ψ (x , y, z); 1,0 1 1 1 u (x , y, z) = σ + (σ − σ )ψ (x , y, z), 2,0 2 2 2 where ψ(x , y, z) := 1 − 25(min(|y|, )) fulfils (u (x), u (x)) ∈  for all x ∈ Γ(0) 1,0 2,0 and is shown in the first snap shot of Fig. 5. We solve the problem with τ = 1e-3 on a sequence of seven meshes Γ (t ), i = 1,..., 7, with decreasing initial meshsizes h (0), h,i i 123 J Sci Comput Fig. 5 Numerical Example 2 on the Thomas RDS (127): snap shots of the U component of the LESFEM– IMEX Euler numerical solution at different times Table 2 Numerical Example 2 on the Thomas RDS (127): for ESFEM,  is not invariant for any i = 1,..., 7, as the method violates the minimum of U for i = 5, 6, 7 and the minimum of U for each i = 1,..., 7 1 2 ih (0) min η min η min η min η i t ∈[τ,T ] t ∈[τ,T ] 1 t ∈[τ,T ] t ∈[τ,T ] 2 1 2 1 1.190e+00 8.482e−02 1.841e−01 − 8.001e−01 5.214e−01 2 8.537e−01 8.005e−02 5.281e−01 − 5.318e−01 6.397e−01 3 5.898e−01 5.860e−02 1.759e−01 − 5.083e−01 7.626e−01 4 4.273e−01 2.259e−02 1.793e−01 − 6.559e−01 5.147e−01 5 3.011e−01 − 7.288e−02 1.791e−01 − 6.338e−01 8.181e−01 6 2.114e−01 − 3.555e−01 1.241e−01 − 5.017e−01 4.297e−01 7 1.531e−01 − 5.376e−01 1.321e−01 − 3.816e−01 6.152e−01 i = 1,..., 7, with both the LESFEM–IMEX Euler and ESFEM-IMEX Euler methods. Snap shots of the LESFEM–IMEX Euler numerical solution obtained on the finest mesh Γ at h,7 different times are shown in Fig. 5. In particular, at the final time T = 100 (see the last snap shot of Fig. 5), the surface is stationary up to machine precision and the numerical solution has reached a stationary pattern. For a given numerical solution (U , U ) on the mesh Γ , 1 2 h,i i = 1,..., 7, consider the following functions η (t ) := min (σ − U (x, t )), η (t ) := min (U (x, t ) − σ ), k = 1, 2. k k k x∈Γ (t ) x∈Γ (t ) h,i h,i These functions are the oriented distances of the numerical solution (U , U ) from the edges 1 2 of . If the oriented distances η and η , k = 1, 2, stay positive at all times, it means that (U , U ) is in  at all times. For all i = 1,..., 7, we show the minima over the time interval 1 2 [τ, T ] (i.e. excluding the initial data) of η and η , k = 1, 2, for both the ESFEM-IMEX Euler and LESFEM–IMEX Euler methods in Tables 2 and 3, respectively. We observe that ESFEM-IMEX Euler violates  for all i = 1,..., 7, while LESFEM–IMEX Euler preserves for all i = 1,..., 7. 123 J Sci Comput Table 3 Numerical Example 2 on the Thomas RDS (127): for LESFEM, the region  is invariant for all i = 1,..., 7 ih (0) min η min η min η min η i t ∈[τ,T ] t ∈[τ,T ] 1 t ∈[τ,T ] t ∈[τ,T ] 2 1 2 1 1.190e+00 1.060e−01 1.738e−01 7.6017752487788e−02 2.665e−01 2 8.537e−01 1.060e−01 4.796e−01 7.6017752487775e−02 3.676e−01 3 5.898e−01 1.060e−01 1.958e−01 7.6017752487774e−02 5.336e−01 4 4.273e−01 1.060e−01 2.163e−01 7.6017752487768e−02 6.679e−01 5 3.011e−01 1.060e−01 1.934e−01 7.6017752487765e−02 9.708e−01 6 2.114e−01 1.060e−01 1.909e−01 7.6017752487761e−02 9.401e−01 7 1.531e−01 1.060e−01 1.899e−01 7.6017752487756e−02 9.367e−01 The minima of η (and thus the minima of U ) coincide up to machine precision 9 Conclusions In this paper we have considered a lumped evolving surface finite element method for the spatial discretisation of reaction–diffusion systems on evolving surfaces, by extending substantially the counterpart on stationary surfaces studied in [16]. We have obtained a fully- discrete scheme by applying the IMEX Euler timestepping to the spatially discrete method. In Sect. 4 we have presented a characterisation of the tangential flow of the continuous and discrete material velocities on evolving smooth and triangulated surfaces, respectively. As a consequence we have obtained, in Theorem 1, an explicit expression for the dilation rates for the special case of isotropic growth. In Theorems 2 and 3 of Sect. 5 we have presented sufficient conditions for the linear heat equation on evolving surfaces to fulfil the maximum principle at the semi- and fully-discrete levels, respectively. In particular, at the fully-discrete level, no timestep restriction is needed if the source parameter β in (70) is nonpositive. In Theorems 4 and 5 of Sect. 6 we provided sufficient conditions under which a hyper-rectangle in the phase space of an RDS on evolving surfaces is invariant at the semi- and fully-discrete levels, respectively. In particular, at the fully-discrete level, a timestep restriction depending on the Lipschitz constants of the kinetics is needed. In Sect. 7 we classified some families of invariant regions for two well-known RD models on evolving surfaces: the activator-depleted and the Thomas models. Two numerical examples are presented in Sect. 8 that experimentally show (i) the optimal convergence (i.e. quadratic in space and linear in time) of the proposed method for the linear heat equation on a linearly evolving sphere and (ii) the existence of a bounded invariant region for an RDS with Thomas kinetics on a logistically growing Dupin ring cyclide at the discrete level. The mathematical and numerical analysis of more complicated settings such as non- isotropic chemically driven is a non-trivial exercise that requires new mathematical tools and techniques. Such a study is beyond the scope of this paper. In order to provide a first step in this direction, we have considered the simplest form of growth, uniform isotropic, which has allowed us to state precisely the conditions necessary for the preservation of invariant regions as well as stating the conditions for positivity of solutions. In future studies we wish to extend the method by considering an ALE approach as in [13] in order to preserve good mesh properties under evolution. This would involve the study of convection RDSs on evolving surfaces, where the convection term arises from the ALE mapping. 123 J Sci Comput Acknowledgements The authors (MF, AM, IS CV) would like to thank the Isaac Newton Institute for Math- ematical Sciences for its hospitality during the programme [Coupling Geometric PDEs with Physics for Cell Morphology, Motility and Pattern Formation] supported by EPSRC Grant Number EP/K032208/1. MF’s and IS’s research work has been performed under the auspices of the Italian National Group for Scientific Calculus (GNCS-INdAM). This work (AM, CV) is partly supported by the EPSRC Grant Number EP/J016780/1 and the Leverhulme Trust Research Project Grant (RPG-2014-149). AM acknowledges funding from the European Union Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant Agree- ment No. 642866, the Commission for Developing Countries, and was partially supported by a grant from the Simons Foundation. AM is a Royal Society Wolfson Research Merit Award Holder funded generously by the Wolfson Foundation. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 Interna- tional License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Appendix: Derivation of the Analytical Solution to Experiment 8.1 2 2 Let ψ : S → R be an eigenfunction of the Laplace-Beltrami operator on S with eigen- value λ ∈ R, i.e. Δ 2 ψ = λψ. Let the surface evolve with an arbitrary isotropic velocity field (48)and let φ(t ) be the respective growth function, as defined in (61). For each t ∈[0, T ], Γ(t ) is a sphere of radius φ(t ). By using the chain rule, it is easy to show that x λ x Δ ψ = ψ ,(x, t ) ∈ G , (131) Γ(t ) T φ(t ) φ (t ) φ(t ) i.e. ψ is an eigenfunction of Δ and its eigenvalue decays quadratically with the Γ(t ) φ(t ) growth function. We look for solutions to (124)ofthe form u(x, t ) = η(t )ψ ,(x, t ) ∈ G , (132) φ(t ) where η(t ) is an unknown time-dependent coefficient. By using (63), (131)and (132)in (124), and cancelling ψ on both sides, we have φ(t ) $ % λd η( ˙ t ) = η(t ) −β − 2S(t ) + , t ∈[0, T ], (133) φ (t ) which leads to $ % λd η(t ) = η(0) exp −β − 2S(τ ) + dτ, t ∈[0, T ]. (134) φ (τ ) In (134) we choose • a linear growth with r = 1 (hence S(t ) = and φ(t ) = t + 1for t ∈[0, 1],see t +1 Table 1); • η(0) = 1, β =−1and d = for illustrative purposes. The profile η(t ) defined in (134) thus becomes λt η(t ) = exp t − 2ln(t + 1) + , t ∈[0, 1]. (135) 12(t + 1) 123 J Sci Comput By combining (132) with (135), we have x λt u(x , y, z, t ) = ψ exp t − 2ln(t + 1) + ,(x , y, z, t ) ∈ G . (136) t + 1 12(t + 1) Finally, we choose the eigenfunction ψ(x , y, z) := xyz, (x , y, z) ∈ S , whose eigenvalue is λ =−12. Hence, (136) becomes xyz t u(x , y, z, t ) = exp t − 2ln(t + 1) − ,(x , y, z, t ) ∈ G . (137) (t + 1) t + 1 References 1. Barreira, R., Elliott, C.M., Madzvamuse, A.: The surface finite element method for pattern formation on evolving biological surfaces. J. Math. Biol. 63(6), 1095–1119 (2011) 2. Bertalmıo, M., Cheng, L.T., Osher, S., Sapiro, G.: Variational problems and partial differential equations on implicit surfaces. J. Comput. Phys. 174(2), 759–780 (2001) 3. Blazakis, K.N., Madzvamuse, A., Reyes-Aldasoro, C.C., Styles, V., Venkataraman, C.: Whole cell tracking through the optimal control of geometric evolution laws. J. Comput. Phys. 297, 495–514 (2015) 4. Chellaboina, V., Bhat, S.P., Haddad, W.M., Bernstein, D.S.: Modeling and analysis of mass-action kinetics. IEEE Control Syst. 29, 60–78 (2009) 5. Chen, C, Kublik, C, Tsai, R.: An implicit boundary integral method for interfaces evolving by Mullins– Sekerka dynamics. In: Maekawa, Y., Jimbo, S. (eds.) Mathematics for Nonlinear Phenomena—Analysis and Computation, pp. 1–21. Springer (2017) 6. Chueh, K.N., Conley, C.C., Smoller, J.A.: Positively invariant regions for systems of nonlinear diffusion equations. Indiana Univ. Math. J. 26, 373–392 (1977) 7. Corson, F., Hamant, O., Bohn, S., Traas, J., Boudaoud, A., Couder, Y.: Turning a plant tissue into a living cell froth through isotropic growth. Proc. Nat. Acad. Sci. 106(21), 8453–8458 (2009) 8. Crampin, E.J., Gaffney, E.A., Maini, P.K.: Reaction and diffusion on growing domains: scenarios for robust pattern formation. Bull. Math. Biol. 61(6), 1093–1120 (1999) 9. Dziuk, G., Elliott, C.M.: Finite elements on evolving surfaces. IMA J. Numer. Anal. 27(2), 262–292 (2007) 10. Dziuk, G., Elliott, C.M.: Finite element methods for surface PDEs. Acta Numer. 22, 289–396 (2013) 11. Dziuk, G., Elliott, C.M.: L -estimates for the evolving surface finite element method. Math. Comput. 82(281), 1–24 (2013) 12. Elliott, C.M., Stinner, B., Venkataraman, C.: Modelling cell motility and chemotaxis with evolving surface finite elements. J. R. Soc. Interface 9(76), 3027–3044 (2012) 13. Elliott, C.M., Venkataraman, C.: Error analysis for an ALE evolving surface finite element method. Numer. Methods Partial Differ. Equ. 31(2), 459–499 (2015) 14. Frittelli, M: Numerical Methods for Partial Differential Equations on Stationary and Evolving Surfaces. Ph.D. thesis, Dipartimento di Matematica e Fisica “E. De Giorgi”, Università del Salento (2018) (to appear) 15. Frittelli, M., Madzvamuse, A., Sgura, I., Venkataraman, C.: Lumped finite elements for reaction-cross- diffusion systems on stationary surfaces. Comput. Math. Appl. 74(12), 3008–3023 (2017). https://doi. org/10.1016/j.camwa.2017.07.044 16. Frittelli, M., Madzvamuse, A., Sgura, I., Venkataraman, C.: Preserving invariance properties of reaction– diffusion systems on stationary surfaces. IMA J. Numer. Anal. (2017). https://doi.org/10.1093/imanum/ drx058 17. Frittelli, M., Sgura, I.: Virtual element method for the Laplace–Beltrami equation on surfaces. ESAIM Math. Model. Numer. Anal. (2017). https://doi.org/10.1051/m2an/2017040 18. Fuselier, E.J., Wright, G.B.: A high-order kernel method for diffusion and reaction–diffusion equations on surfaces. J. Sci. Comput. 56(3), 535–565 (2013) 19. González-Olivares, E., Ramos-Jiliberto, R.: Dynamic consequences of prey refuges in a simple model system: more prey, fewer predators and enhanced stability. Ecol. Model. 166, 135–146 (2003) 20. Grande, J., Reusken, A.: A higher order finite element method for partial differential equations on surfaces. SIAM J. Numer. Anal. 54(1), 388–414 (2016) 21. Hebey, E., Robert, F.: Sobolev spaces on manifolds. In: Krupka, D., Saunders, D. (eds.) Handbook of Global Analysis, pp. 375–415 (2008) 123 J Sci Comput 22. Kovács, B.: Computing arbitrary Lagrangian Eulerian maps for evolving surfaces. arXiv preprint arXiv:1612.01701 (2016) 23. Kovács, B., Li, B., Lubich, C., Guerra, C.A.P.: Convergence of finite elements on an evolving surface driven by diffusion on the surface. Numer. Math. 137(3), 643–689 (2017) 24. Lacitignola, D., Bozzini, B., Frittelli, M., Sgura, I.: Turing pattern formation on the sphere for a mor- phochemical reaction–diffusion model for electrodeposition. Commun. Nonlinear Sci. Numer. Simul. 48, 484–508 (2017). https://doi.org/10.1016/j.cnsns.2017.01.008 25. Lacitignola, D., Bozzini, B., Sgura, I.: Spatio-temporal organization in a morphochemical electrode- position model: Hopf and Turing instabilities and their interplay. Eur. J. Appl. Math. 26(2), 143–173 (2015) 26. Lefèvre, J., Mangin, J.F.: A reaction–diffusion model of human brain development. PLoS Comput. Biol. 6(4), e1000,749 (2010) 27. Li, H., Lin, Y., Heath, R.M., Zhu, M.X., Yang, Z.: Control of pollen tube tip growth by a rop gtpase- dependent pathway that leads to tip-localized calcium influx. Plant Cell 11(9), 1731–1742 (1999) 28. Madzvamuse, A., Barreira, R.: Exhibiting cross-diffusion-induced patterns for reaction-diffusion systems on evolving domains and surfaces. Phys. Rev. E 90(4), 043,307 (2014) 29. Madzvamuse, A., Gaffney, E.A., Maini, P.K.: Stability analysis of non-autonomous reaction-diffusion systems: the effects of growing domains. J. Math. Biol. 61(1), 133–164 (2010) 30. Madzvamuse, A., Wathen, A.J., Maini, P.K.: A moving grid finite element method applied to a model biological pattern generator. J. Comput. Phys. 190(2), 478–500 (2003) 31. Murray, J.D.: Mathematical Biology. II Spatial Models and Biomedical Applications (Interdisciplinary Applied Mathematics), vol. 18. Springer, New York (2001) 32. Plaza, R.G., Sanchez-Garduno, F., Padilla, P., Barrio, R.A., Maini, P.K.: The effect of growth and curvature on pattern formation. J. Dyn. Differ. Equ. 16(4), 1093–1121 (2004) 33. Pruyne, D., Bretscher, A.: Polarization of cell growth in yeast. I. Establishment and maintenance of polarity states. J. Cell Sci. 113(3), 365–375 (2000) 34. Quarteroni, A., Valli, A.: Numerical Approximation of Partial Differential Equations, vol. 23. Springer, New York (2008) 35. Ren, K., Tsai, R., Zhong, Y.: An implicit boundary integral method for computing electric potential of macromolecules in solvent. arXiv:1709.08070v4 (2018) 36. Sahlin, P., Jönsson, H.: A modeling study on how cell division affects properties of epithelial tissues under isotropic growth. PLoS ONE 5(7), e11,750 (2010) 37. Schnakenberg, J.: Simple chemical reaction systems with limit cycle behaviour. J. Theor. Biol. 81(3), 389–400 (1979) 38. Smoller, J.: Shock Waves and Reaction Diffusion Equations. Springer, New York (1994) 39. Taylor, M.E.: Partial Differential Equations. III. Springer, New York (2011) 40. Tuncer, N., Madzvamuse, A.: Projected finite elements for systems of reaction–diffusion equations on closed evolving spheroidal surfaces. Commun. Comput. Phys. 21(3), 718–747 (2017) 41. Venkataraman, C., Lakkis, O., Madzvamuse, A.: Global existence for semilinear reaction–diffusion sys- tems on evolving domains. J. Math. Biol. 64(1–2), 41–67 (2012) 42. Yang, F.W., Venkataraman, C., Styles, V., Madzvamuse, A.: A robust and efficient adaptive multigrid solver for the optimal control of phase field formulations of geometric evolution laws. Commun. Comput. Phys. 21(1), 65–92 (2017)

Journal

Journal of Scientific ComputingSpringer Journals

Published: Jun 1, 2018

References

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off