J Sci Comput (2018) 77:971–1000 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 / Published online: 1 June 2018 © The Author(s) 2018 Abstract We propose and analyse a ﬁnite 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 ﬁeld. 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 sufﬁcient 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 speciﬁc 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 speciﬁc 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 972 J Sci Comput (2018) 77:971–1000 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 ﬂat 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 ﬁeld 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 ﬁnite 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, sufﬁcient 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 ﬁeld. 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 ﬁnite 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 sufﬁcient conditions under which a region is invariant at the semi- and fully-discrete levels when a lumped evolving surface ﬁnite 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 ﬁnd 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 sufﬁcient 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 (2018) 77:971–1000 973 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 ﬁrst 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, sufﬁcient 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 ﬁndings 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 ﬁeld 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 deﬁned by G := Γ(t ) ×{t }. The material velocity T T t ∈[0,T ] v : G → R of Γ(t ) is deﬁned 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 974 J Sci Comput (2018) 77:971–1000 ∂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 ) deﬁned 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 fulﬁls 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 sufﬁciently 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 deﬁned by ∂g ˜ ∂ g := + v ·∇˜ g, (6) ∂t where ∇ is the standard gradient in R and g ˜ is any differentiable extension of g deﬁned on a neighborhood of G . Deﬁnition (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 sufﬁciently 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 ). Speciﬁcally, 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 sufﬁciently 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 sufﬁciently 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 (2018) 77:971–1000 975 Remark 1 (Surfaces without boundary) Lemmas 2–4 hold on surfaces with or without bound- ary, i.e. ∂Γ (t ) =∅ or ∂Γ (t ) =∅, respectively. Speciﬁcally, 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 ﬂuxes 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 ﬂux 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 ﬂuxes q , k = 1,..., r, are tangent to Γ(t ), we can apply the integration-by-parts formula (8) to the ﬁrst 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 deﬁned 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 ﬂux 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 sufﬁcient condition under which system (16) possesses an invariant region. To this end, we give the following deﬁnitions. Deﬁnition 1 (Dilation rates)The minimum and maximum instantaneous dilation rates are deﬁned 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 976 J Sci Comput (2018) 77:971–1000 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 deﬁned 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 Deﬁnition 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 deﬁne the constants ∗ ∗ μ if σ ≥ 0, μ if σ ≥ 0, ∗ k ∗ min max k μ := μ := (20) ∗ ∗ μ if σ < 0, μ if σ < 0, max k min ∗ ∗ where μ and μ are deﬁned 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 ﬂat 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 fulﬁlled, respectively. ∗ ∗ Notice that on stationary surfaces, since μ = μ = 0, then conditions (21)–(22) reduce to the inward ﬂux 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 (2018) 77:971–1000 977 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 ﬂux 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. Deﬁnition 3 (Maximum principle)For k = 1, the scalar equation (16) fulﬁls the maximum principle if, for any initial condition, i.e. u(·, 0), the solution u fulﬁls 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 fulﬁls 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 deﬁnition 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 ﬁrst 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 978 J Sci Comput (2018) 77:971–1000 for all k = 1,..., r. Therefore, the variational formulation seeks to ﬁnd 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 ﬁnite element method (ESFEM) studied in [1] for the approx- imation of the variational problem (29), we present its lumped counterpart, the lumped evolving surface ﬁnite element method (LESFEM). 3.1 Surface Triangulation and Some Deﬁnitions Following [9], given h > 0, called meshsize, a triangulation Γ (t ) of the evolving surface Γ(t ) is deﬁned 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 deﬁned in Lemma 1; h δ δ –Forall t ∈[0, T ], the normal projection a(·, t ) : U (t ) → Γ(t ) deﬁned 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, deﬁned 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 deﬁned 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 ﬁxed time t ∈[0, T ],let V (t ) be the space of piecewise linear functions on Γ (t ) h h deﬁned by V (t ) := {ϕ ∈ C (Γ (t )) ϕ is linear afﬁne for each Z ∈ Z }. (31) h h | h 123 J Sci Comput (2018) 77:971–1000 979 Let V be the space of time-dependent, spatially piecewise linear functions deﬁned 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 sufﬁciently smooth function U ∈ V is deﬁned by ∂U ∂ U := + I (v) ·∇U,(x, t ) ∈ G , h h,T ∂t where v is the material velocity. For our purposes, we deﬁne the discrete counterpart of the dilation rates introduced in Deﬁnition 1. Deﬁnition 4 (Discrete dilation rates)The minimum and maximum discrete instantaneous dilation rates are deﬁned 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 deﬁned 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, deﬁned in (35) fulﬁl ∂ χ = 0, i = 1,..., N . (37) Hence, for the functions U ,k = 1,..., r deﬁned 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 980 J Sci Comput (2018) 77:971–1000 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 ﬁnite element method (LESFEM), applied to the variational formulation (29), seeks to ﬁnd 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: ﬁnd 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 ﬁnite element method (ESFEM) in [1] due to the presence of the interpolant operator on the ﬁrst 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 deﬁned 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 (2018) 77:971–1000 981 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 coefﬁcients 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 coefﬁcients 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 speciﬁc 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 μ deﬁned 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 sufﬁcient 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 ﬂows ∇ ·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 deﬁned 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 fulﬁl the following property I (χ χ ) = I (χ χ )∇ · I (v), ∀i = 1,..., N . (49) h i j h i j Γ h dt Γ (t ) Γ (t ) h h 123 982 J Sci Comput (2018) 77:971–1000 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 fulﬁls 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 deﬁned 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 ﬂows ∇ · 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 (2018) 77:971–1000 983 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 ﬁeld 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 ﬂow 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 deﬁned in (18)and (34), respectively, for arbitrary surfaces that evolve with the material velocity (48). The velocity ﬁeld (48) corresponds to the speciﬁc 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 speciﬁc 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 ﬁeld (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 984 J Sci Comput (2018) 77:971–1000 ∗ ∗ 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 ﬂow on isotropically growing smooth or triangulated surfaces) Let Γ(t ) be a smooth surface that evolves with the velocity ﬁeld (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 deﬁned in Sect. 4.2.From(60)and (61), J (θ , t ) fulﬁls 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 ﬁrst 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 speciﬁc 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 speciﬁc 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 (2018) 77:971–1000 985 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 proﬁles Theorem 2 (Semi-discrete maximum principle for the linear heat equation (70)) If the veloc- ity ﬁeld v fulﬁls μ + β ≥ 0, (71) min with μ as deﬁned in (34), and the triangulation Γ (t ) meets the Delaunay condition for min h all t ∈[0, T ], then the LESFEM solution of (70) fulﬁls 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 986 J Sci Comput (2018) 77:971–1000 −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 fulﬁls −1 (43) from the Delaunay condition. Then, it follows that −dM A|ξ|≤ 0. Hence, it sufﬁces −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 ﬁeld v fulﬁls μ + β ≥ 0, ∀t ≥ 0, (78) min with μ as deﬁned in (34), and the triangulation Γ meets the Delaunay condition for min h all t > 0, then the LESFEM–IMEX Euler solution of (70) fulﬁls 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 satisﬁes τβ ≤ 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) fulﬁls 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 (2018) 77:971–1000 987 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 ﬁrst 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 deﬁne f (τ ) := 1−τβ and g(τ ) = e . These functions fulﬁl 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 fulﬁlled for any sufﬁciently small τ . –if f (0) ≤ g (0), then condition (89) is fulﬁlled 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 justiﬁed 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 988 J Sci Comput (2018) 77:971–1000 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 fulﬁls 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 deﬁne the constants μ if σ ≥ 0, μ if σ ≥ 0, min k max μ := μ := (96) μ if σ < 0, μ if σ < 0, max k min where μ and μ are the dilation rates deﬁned 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 ) satisﬁes 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 sufﬁces to verify that, for all (U ,..., U ) ∈ , k = 1,..., r,and i = 1,..., N, 1 r 123 J Sci Comput (2018) 77:971–1000 989 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 sufﬁcient 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 deﬁned 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 τ fulﬁls τ 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 sufﬁces 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 990 J Sci Comput (2018) 77:971–1000 with μ as deﬁned 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 speciﬁc 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 fulﬁlled 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 μ deﬁned in (18). The results in this section min max are conﬁned 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 (2018) 77:971–1000 991 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 speciﬁc 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 μ deﬁned 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 992 J Sci Comput (2018) 77:971–1000 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 μ deﬁned 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 (2018) 77:971–1000 993 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 ﬁnd 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 ﬁnal time. The initial domain Γ(0) is the unit sphere S ,that evolves under the velocity ﬁeld 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 fulﬁls μ = = 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 μ = 1fulﬁll(78)and (80) for each τ> 0. Hence, the min LESFEM–IMEX Euler solution to (124) fulﬁls 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 = 0fulﬁl h (0) i −1 h (0) = 0.4013 and h (0) ≈ , i = 2,..., 8. The corresponding timesteps τ fulﬁl 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 ﬁnest 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 994 J Sci Comput (2018) 77:971–1000 Fig. 2 Numerical Example 1 on the linear heat equation (124): numerical solution at different times obtained on the ﬁnest 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 ﬁgure 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 coefﬁcient 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 ﬁeld (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 ﬁnal time we choose T = 100. From 123 J Sci Comput (2018) 77:971–1000 995 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 modiﬁed kinetics deﬁned in (129)–(130) Theorem 1, it follows that the corresponding dilation rates fulﬁl ∗ ∗ μ = μ ≈ 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 modiﬁed inward ﬂux 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 modiﬁed 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)fulﬁl 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|, )) fulﬁls (u (x), u (x)) ∈ for all x ∈ Γ(0) 1,0 2,0 and is shown in the ﬁrst 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 996 J Sci Comput (2018) 77:971–1000 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 ﬁnest mesh Γ at h,7 different times are shown in Fig. 5. In particular, at the ﬁnal 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 (2018) 77:971–1000 997 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 ﬁnite 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 ﬂow 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 sufﬁcient conditions for the linear heat equation on evolving surfaces to fulﬁl 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 sufﬁcient 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 classiﬁed 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 ﬁrst 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 998 J Sci Comput (2018) 77:971–1000 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 Scientiﬁc 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 ﬁeld (48)and let φ(t ) be the respective growth function, as deﬁned 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 coefﬁcient. 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 proﬁle η(t ) deﬁned in (134) thus becomes λt η(t ) = exp t − 2ln(t + 1) + , t ∈[0, 1]. (135) 12(t + 1) 123 J Sci Comput (2018) 77:971–1000 999 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 ﬁnite 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 ﬁnite 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 ﬁnite elements. J. R. Soc. Interface 9(76), 3027–3044 (2012) 13. Elliott, C.M., Venkataraman, C.: Error analysis for an ALE evolving surface ﬁnite 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 ﬁnite 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 ﬁnite 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 1000 J Sci Comput (2018) 77:971–1000 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 ﬁnite 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 inﬂux. 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 ﬁnite 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 ﬁnite 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 efﬁcient adaptive multigrid solver for the optimal control of phase ﬁeld formulations of geometric evolution laws. Commun. Comput. Phys. 21(1), 65–92 (2017)
Journal of Scientific Computing – Springer Journals
Published: Jun 1, 2018
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
Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly
Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.
Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.
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.
“Hi guys, I cannot tell you how much I love this resource. Incredible. I really believe you've hit the nail on the head with this site in regards to solving the research-purchase issue.”
Daniel C.
“Whoa! It’s like Spotify but for academic articles.”
@Phil_Robichaud
“I must say, @deepdyve is a fabulous solution to the independent researcher's problem of #access to #information.”
@deepthiw
“My last article couldn't be possible without the platform @deepdyve that makes journal papers cheaper.”
@JoseServera
DeepDyve Freelancer | DeepDyve Pro | |
---|---|---|
Price | FREE | $49/month |
Save searches from | ||
Create lists to | ||
Export lists, citations | ||
Read DeepDyve articles | Abstract access only | Unlimited access to over |
20 pages / month | ||
PDF Discount | 20% off | |
Read and print from thousands of top scholarly journals.
Already have an account? Log in
Bookmark this article. You can see your Bookmarks on your DeepDyve Library.
To save an article, log in first, or sign up for a DeepDyve account if you don’t already have one.
All DeepDyve websites use cookies to improve your online experience. They were placed on your computer when you launched this website. You can change your cookie settings through your browser.
ok to continue