A Petrov–Galerkin finite element method for simulating chemotaxis models on stationary surfaces

Computers and Mathematics with Applications xxx (xxxx) xxx

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

A Petrov–Galerkin finite element method for simulating chemotaxis models on stationary surfaces✩ Shubo Zhao a , Xufeng Xiao a , Jianping Zhao a,b , Xinlong Feng a,b , a b

∗

College of Mathematics and System Sciences, Xinjiang University, Urumqi 830046, PR China Institute of Mathematics and Physics, Xinjiang University, Urumqi 830046, PR China

article

info

Article history: Received 5 June 2019 Received in revised form 4 December 2019 Accepted 20 January 2020 Available online xxxx Keywords: Surface chemotaxis models Surface finite element method Petrov–Galerkin method Positivity preservation Mass conservation

a b s t r a c t In this paper, we present a Petrov–Galerkin finite element method for a class of chemotaxis models defined on surfaces, which describe the movement by one community in reaction to one chemical or biological signal on manifolds. It is desired for numerical methods to satisfy discrete maximum principle and discrete mass conservation property, which is a challenge due to the singular behavior of numerical solution. Thus a Petrov– Galerkin method is combined with an effective mass conservation factor to overcome the challenge. Furthermore, we prove two facts, this method maintains positivity and discrete mass conservation property. In addition, decoupled approach is applied based on the gradient and Laplacian recoveries to solve the coupling system. The relevant stability analyses is provided. Finally, numerical simulations of blowing-up problems and pattern formulations demonstrate the effectiveness of the proposed method. © 2020 Elsevier Ltd. All rights reserved.

1. Introduction Chemotaxis models were firstly introduced by Keller and Segel [1,2] for mathematical description of chemotactic processes. The processes describe species/cells sense and migrate towards/backwards the domains of higher concentrations of chemical/biological signals such as toxins, partners, preys and predators. The models play a vitally important role in simulating and predicting species migration and microbial activity. The density of species/cells rapidly grows in a small neighborhood of certain points or curves, since they continually agglomerate in the high or low concentration region. So the density blows up at accumulation region with the effect of chemotactic concentration [3–6]. In the case of numerical simulation, the blow-up phenomenon or steep gradient behavior of exact solution may lead to nonphysical oscillations or even negative densities/concentrations, if the employed numerical algorithm does not satisfy the discrete maximum principle (DMP). Thus, capturing the blow-up behavior of numerical solution is a challenging problem. For the existing numerical methods for chemotaxis partial differential equations (PDEs), most of them provide high resolution and positivity preserving schemes. The prevalent numerical techniques include finite element methods [7,8], flux-corrected finite element methods [9–12], finite volume methods [13–16] and discontinuous Galerkin methods [17–21]. ✩ This work is in parts supported by the NSF of China (No. U19A2079, No. 11671345, No. 11971377) and the Xinjiang Provincial University Research Foundation of China. ∗ Corresponding author at: Institute of Mathematics and Physics, Xinjiang University, Urumqi 830046, PR China. E-mail addresses: [email protected] (S. Zhao), [email protected] (X. Xiao), [email protected] (J. Zhao), [email protected] (X. Feng). https://doi.org/10.1016/j.camwa.2020.01.019 0898-1221/© 2020 Elsevier Ltd. All rights reserved.

Please cite this article as: S. Zhao, X. Xiao, J. Zhao et al., A Petrov–Galerkin finite element method for simulating chemotaxis models on stationary surfaces, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.019.

2

S. Zhao, X. Xiao, J. Zhao et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

For simulating complex models for biological and medical applications, we need to consider the migration of species/cells on manifolds. There are a lot of discussions on numerical methods for chemotaxis models in two and three dimensions. Nevertheless, it is challenging to develop efficient numerical methods for such equations on manifolds. Some simulating results can be found in [22–24], and these results are obtained by finite element methods and flux-corrected finite element methods. In this paper, we present a Petrov–Galerkin finite element method for the chemotaxis models defined on surface. The method, called Mizukami–Hughes method, was introduced by Mizukami and Hughes in [25] for conforming linear triangular finite elements in planar case. They showed that the method satisfies DMP and produces nonphysical oscillation solution for the steady state convection-diffusion equation. This method is based on the observation that the velocity field can be changed in the direction perpendicular to the gradient of solution without affecting the solution of original equation. Knobloch extended the method to steady state convection-diffusion-reaction equations and three-dimensional problems [26]. The method depends on the unknown discrete solution, so the resulting method is nonlinear. In this paper, original question is rewritten as a convection-diffusion-reaction type PDEs for the convenience of using Mizukami–Hughes method. The spatial discretization is the surface finite element method [27–31]. Mizukami–Hughes method is applied to construct the non-negative type matrix. The mass matrix of standard surface finite element is modified as the lumped mass form [32,33]. These operations can ensure the positivity of numerical solution. Although the rewritten formula is very convenient for Mizukami–Hughes method, it does not satisfy the mass conservation property. The reason to the problem is that the co-normal vector τh to an element edge are not necessarily collinear. As described above, Mizukami–Hughes method is an effective method for singular solution problem, and we expect that the method can be an available tool for chemotaxis models. Thus we combine mass correction method and Mizukami–Hughes method to remedy the flaws [34]. Inspired by conservative Allen-Cahn equation [35] and mass conservation for level set method [36], we construct a suitable correction factor that can produce non-oscillation solution and maintain the discrete mass conservation property. The singular behavior of numerical solution makes the construction of appropriate mass correction factor challenging. Therefore, the characteristics of numerical solution should be considered when constructing the mass correction equation. This paper is organized as follows. In Section 2, we define the chemotaxis models on general surfaces. Then, we briefly introduce the surface finite element method for the chemotaxis models in Section 3. In Section 4, the Mizukami–Hughes method is combined with a mass correction method to solve the chemotaxis problems. The corresponding positivity preservation analysis results are shown in Section 5. In Section 6, the demonstrated numerical examples are performed to test the proposed method. And some concluding remarks are given in the last section. 2. Problem formulation The generic form of chemotaxis model defined on surfaces can be mathematically described using the following system

⎧ ut = αu ∆Γ u − ∇Γ · (χ (c)u∇Γ c) + f (u)u, ⎪ ⎪ ⎨ c = α ∆ c − β c + g(u)u, t c Γ ⎪ u|t =0 = u0 , c |t =0 = c0 , ⎪ ⎩ nΓ · ∇Γ u = 0, nΓ · ∇Γ c = 0,

x ∈ Γ , t ∈ (0, T ],

(a)

x ∈ Γ , t ∈ (0, T ],

(b)

x ∈ Γ,

(c)

x ∈ ∂ Γ or ∂ Γ = ∅, t ∈ (0, T ],

(d)

(2.1)

where the surface concerned Γ can be expressed as a zero-level set function form

Γ = {x = (x, y, z) ∈ U |φ (x) = 0}

(2.2)

U ⊂ R is an open tube domain and the function φ ∈ C (U ) satisfies ∇φ (x) ̸ = 0. u(x, t) denotes the species densities and c(x, t) stands for the chemoattractant concentration, both of which are defined on the general surface Γ . αu and αc are the non-negative diffusion coefficients. β ≥ 0 describes the depletion of chemoattractant. χ (c) ≥ 0 is the chemoattractant sensitivity. f (u) denotes the rate of growth of species. g(u) ≥ 0 denotes the rate of production of chemoattractant. nΓ is the co-normal vector which is normal to ∂ Γ and tangent to Γ . For the surface gradient operator, we have the following definition 3

2

∇Γ η = (I − n(x)n(x)T )∇η = P (x)∇η,

(2.3)

3

here, n(x) ∈ R is a unit normal vector and satisfies n(x) =

∇φ (x) . |∇φ (x)|

Thus, the surface divergence operator is given by ∇Γ · η = trace(∇Γ η), and the Laplace-Beltrami operator ∆Γ on surface can be defined as ∇Γ · ∇Γ . Eqs. (2.1)(c) and (2.1)(d) are the initial condition and the homogeneous Neumann boundary condition, respectively. Notice that the boundary condition vanishes when Γ is a closed surface, ∂ Γ = ∅. In this work, we consider the system (2.1) defined on a closed surface. Please cite this article as: S. Zhao, X. Xiao, J. Zhao et al., A Petrov–Galerkin finite element method for simulating chemotaxis models on stationary surfaces, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.019.

S. Zhao, X. Xiao, J. Zhao et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

3

In particular, we can get the following chemotaxis models defined on surfaces by choosing the special coefficients.

{ Keller–Segel model:

ut = ∆Γ u − ∇Γ · (χ (c)u∇Γ c),

(2.4)

ct = ∆Γ c − c + u.

{ Pattern formation model:

ut = αu ∆Γ u − ∇Γ · (χ (c)u∇Γ c) + u2 (1 − u), ct = ∆Γ c − β c + u.

(2.5)

The solution of the system (2.4) may appear singular, spiky and even blow-up phenomenon in finite time when the initial condition satisfies specific condition [6,37]. For the system (2.5), the blow-up phenomenon is prevented by the reactive term u2 (1 − u), but in general, χ (c) ≫ αu causes strong convection-dominated problem. The theory of blow-up solution has not yet been studied on surface, but the related results can be observed by numerical simulation in this or existing studies [23,24]. Steep gradients, spiky and convection-dominated may rise to nonphysical solution (the species density may become negative) if the numerical method does not satisfy DMP. Thus, it is a challenging problem to construct effective numerical methods for the system (2.1). 3. Surface finite element method In the section, we briefly introduce the surface finite element method for the system (2.1). 3.1. Triangulated surfaces and basic notation We use a piecewise polynomial surface Γh to approximate Γ . The discrete { }N surface Γh is constructed by a set of non-overlapped triangles {ek }M which are generated by a set of vertices xj j=1 on the smooth surface Γ , i.e., k=1

Γh =

M ⋃

ek .

k=1

We assume that the triangulated surface Γh is of Delaunay type and has following important properties: • There is a constant ρ for all ek ∈ Γh , such that maxek ∈Γh hk ≤ ρ minek ∈Γh hk , where hk denotes the length of the longest edge of the triangles ek , and the mesh size h = maxek ∈Γh hk . π • The triangulation Γh is of weakly acute type, that is the magnitude of all the angles in triangles ek not more than . 2 The first assumption amount to a quasi-uniformity property of the triangulation Γh , and the second implies that the stiffness matrix of the Poisson equation, which is discretized using piecewise linear function on Γh , is an M-matrix [38]. There is a piecewise bijection between each point x ∈ Γh and a(x) ∈ Γ , such that x = a(x) + d(x)n(a(x)), where d(x) is the oriented distance function for Γ with |d(x)| = dist(x, Γ ), ∇ d(x) = n(a(x)) and |∇ d(x)| = 1. The theoretical justification can be found in [27,28]. Then, any function η(x) defined on the discrete surface Γh can be lifted to Γ by

ηl (a(x)) = η(x),

a(x) ∈ Γ .

(3.1)

We also can define the inverse lifting

ξ −l (x) = ξ (a(x)),

x ∈ Γh .

l −l Thus, we can get the initial value u− 0 and c0 and other given data defined on the discrete surface Γh by the inverse lifting.

3.2. Surface FEM for the surface chemotaxis problems In practice, Eq. (2.1)(a) may also be rewritten as the following convection-diffusion-reaction model which is convenient for using Mizukami–Hughes method, ut = αu ∆Γ u − χ (c)∇Γ c · ∇Γ u − ∇Γ · (χ (c)∇Γ c)u + f (u)u.

(3.2)

TPhe finite element space based on the triangulated surface Γh is chosen as Sh = Φ ∈ C 0 (Γh ) | Φ is a linear affine for each ek ∈ Γh

{

}

{ } = span ψj | j = 1, . . . , N ⊂ H 1 (Γh ), Please cite this article as: S. Zhao, X. Xiao, J. Zhao et al., A Petrov–Galerkin finite element method for simulating chemotaxis models on stationary surfaces, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.019.

4

S. Zhao, X. Xiao, J. Zhao et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

where N is the number of the vertices of Γh , and ψj denotes the nodal basis function that is ψj (xi ) = δij , (i, j = 1, . . . , N). Thus, the functions Uh ∈ Sh and Ch ∈ Sh can be denoted as Uh =

N ∑

Uh,j ψj (x) and Ch =

j=1

N ∑

Ch,j ψj (x),

x ∈ Γh ,

j=1

real constants Uh,j and Ch,j (j = 1, . . . , N) are the unknown nodal values of the solution Uh and Ch at the vertice xj . n−1 Let 0 = t 0 < t 1 < · · · < t n−1 < t n < · · · < t Nt = T be a ∑ partition with δ t = t n − ∑tN nfor (0, T ]. The finite N n U ψ (x) and C = element approximations of u−l (t n ) and c −l (t n ) are defined as Uh = j h j=1 Ch,j ψj (x), respectively. j=1 h,j Then, the full-discrete variational form of (2.1) can be defined as following by employing the semi-implicit schemes: for ∀Φ ∈ Sh (Γh ), find Uhn , Chn ∈ Sh (Γh ) such that

⎧ n−1 n n n n n n ⎪ ⎨M(Uh , Φ ) + D(Uh , Φ ) + C (Uh , Ch , Φ ) + Ru (Uh , Ch , Φ ) = F (Uh , Φ ), (a) n−1 M(Chn , Φ ) + D(Chn , Φ ) + Rc (Chn , Φ ) = G (Uh , Φ ), (b) ⎪ ⎩ −l −l Uh |t =0 = u0 (x), Ch |t =0 = c0 (x). (c)

(3.3)

Here,

⎧ Ψhn − Ψhn−1 ⎪ n ⎪ , Φ )Γh , Ψhk ∈ {Uhk , Chk }, M ( Ψ , Φ ) = ( ⎪ h ⎪ ⎪ δt n ⎪ { } ⎪ ⎪ ⎪ D(Θhn , Φ ) = (α∇Γh Θhn , ∇Γh Φ )Γh , (α, Θhk ) ∈ (αu , Uhk ), (αc , Chk ) , ⎪ ⎪ ⎨ C (Uhn , Chn , Φ ) = (χ (Chn )∇Γh Chn · ∇Γh Uhn , Φ )Γh , (3.4) ⎪ n n n n n ⎪ ⎪ , C )U , Φ ) · ( χ (C ) ∇ R (U , C , Φ ) = ( ∇ Γh Γh h u Γh ⎪ h h h h ⎪ ⎪ ⎪ ⎪Rc (Chn , Φ ) = (β Chn , Φ )Γh , ⎪ ⎪ ⎪ ⎩ n−1 n−1 n−1 n−1 n−1 n−1 F (Uh , Φ ) = (f (Uh )Uh , Φ )Γh , G (Uh , Φ ) = (g(Uh )Uh , Φ )Γh , ∫ with the inner products (W , V )Γh = Γ WVdσ for any functions W , V ∈ H 1 (Γh ) and ∇Γh is determined by nh that is h the unit normal vector of Γh . Since the surface considered is closed, the boundary terms vanish in above variation form. Finally, we lift the solution Uhn and Chn to elements of the smooth surface Γ as uh = (Uhn )l and ch = (Chn )l , and regard them as numerical solutions of (2.1).

4. The Mizukami–Hughes method and mass correction It is known that standard finite element method is invalid for solving convection-diffusion-reaction equations involving steep gradient solution and convection-dominated in planar case. The similar problems also exist for solving the surface PDEs [39]. As mentioned above, convection-dominated problem and steep gradient solution are the difficulties in solving chemotaxis models. Thus, we introduce a DMP-satisfying Petrov–Galerkin finite element method — Mizukami–Hughes method to eliminate nonphysical oscillation. In addition, a mass correction method is proposed to guarantee the discrete mass conservation. We define the local basic function ψlk (l = 1, 2, 3) on each element ek satisfying ψlk (vm ) = δlm , here, vm (m = 1, 2, 3) are the vertices of ek . The Mizukami–Hughes method is employed to solve the chemotaxis problem (2.1) with the weighting local functions defined on each element

˜k = ψ k + γ , ψ l l l

l = 1, 2, 3,

˜j to denote the weighting global basis function where γl are constants that are determined in later paragraph, and we use ψ ˜ ˜j . When xi and vm represent the same defined on the whole computed domain. Assume that Λj is the support set of ψ

˜k (v ) = ψ ˜j , thus ψ ˜j (xi )|ek holds. point, we have ek ⊂ Λ m l { } ˜ denote the weighting test function satisfying Φ ˜ ∈ ˜ ˜j | j = 1, . . . , N . Thus the numerical Let Φ Sh (Γh ) = span ψ ˜ ∈˜ solution obtained by Mizukami–Hughes method for the system (2.1) can be defined by: for ∀Φ ∈ Sh (Γh ), Φ Sh (Γh ), find Uhn , Chn ∈ Sh (Γh ) such that ⎧ n−1 ˜ n ˜ n ˜ n n ˜ n n ˜ ⎪ ⎨M(Uh , Φ ) + D(Uh , Φ ) + C (Uh , Ch , Φ ) + Ru (Uh , Ch , Φ ) = F (Uh , Φ ), n−1 M(Chn , Φ ) + D(Chn , Φ ) + Rc (Chn , Φ ) = G (Uh , Φ ), ⎪ ⎩ 0 l −l 0 Uh = u− 0 (x), Ch = c0 (x).

(a) (b)

(4.1)

(c)

Please cite this article as: S. Zhao, X. Xiao, J. Zhao et al., A Petrov–Galerkin finite element method for simulating chemotaxis models on stationary surfaces, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.019.

S. Zhao, X. Xiao, J. Zhao et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

5

Let τh be the unit co-normal of element ek . We find that, in fact, the scheme (4.1) does not possess the discrete conservation property in the case f (u) = 0 i.e.,

∫ Γh

Uhn dσh −

∫ = −δ t Γ

∫

Uhn−1 dσh

Γh

(χ∇Γh Chn · ∇Γh Uhn + χ Uhn ∆Γh Chn )dσh

⎞ ⎛h ∫ ∫ ∑∫ χ Uhn ∆Γh Chn dσh ⎠ χ Uhn ∆Γh Chn dσh + = −δ t ⎝ χ∇Γh · (Uhn ∇Γh Chn )dσh − = −δ t

Γh

χ Uhn ∇Γh Chn · τh dlh ̸= 0,

∂ ek

ek ∈Γh

Γh

ek

ek ∈Γh

∑∫

(4.2)

here, the Green formula is applied in the above equation. The inequality comes from the fact that the co-normal vector

τh to an element edge are not necessarily collinear. For simplicity, we let Mn = −δ t

∑∫ ∂ ek

ek ∈Γh

χ Uhn ∇Γh Chn · τh dlh .

(4.3)

In other words, the solution Uhn of the scheme (4.1) breaks mass exchange in the case f (u) ̸ = 0, i.e.,

∫ Γh

Uhn dσh −

= Mn + δ t

∫ Γh

∫ Γh

Uhn−1 dσh

f (Uhn−1 )Uhn−1 dσh ̸ = δ t

(4.4)

∫ Γh

f (Uhn−1 )Uhn−1 dσh .

Thus, we present the following scheme satisfying the discrete conservation property

⎧ ∗ n−1 ˜ ) + D(Uh∗ , Φ ˜ ) + C (Uh∗ , Chn , Φ ˜ ) + Ru (Uh∗ , Chn , Φ ˜ ) = F (Uhn−1, Φ ˜ ), M (Uh , Φ ⎪ ⎪ ⎪ ⎪ ⎨−M∗ (U n , Φ ) = R(θ n , Φ ), h h n−1 ⎪ M(Chn , Φ ) + D(Chn , Φ ) + Rc (Chn , Φ ) = G (Uh , Φ ), ⎪ ⎪ ⎪ ⎩ 0 l −l 0 Uh = u− 0 (x), Ch = c0 (x),

(a) (b)

(4.5)

(c) (d)

here, M∗ (Ψh , Φ ) = (

Uh∗ − Ψh

δt

, Φ )Γh and R(θhn , Φ ) = (θhn , Φ )Γh ,

∑N

(4.6)

n n Uh∗ denotes the intermediate variable and θhn = j=1 θh,j ψj is the mass correction factor. We require that Γh θh dσh = n M − . Actually, (4.5)(b) can be treated as the process of mass correction. We solve Eq. (4.5)(c) first to obtain Chn , and then δt solve Eq. (4.5)(a) to get Uh∗ , finally the numerical solution Uhn can be obtained by solving the correction equation (4.5)(b).

∫

Next, we give the appropriate definition of θhn that allows the scheme (4.5) to satisfy the discrete conservation law,

while keeping the positivity of numerical solution. In general, θhn can be defined in the average formula

θhn = −

Mn

δt

∫

Γh

dσh

.

(4.7)

However, the formula (4.7) may break the positivity of numerical solution in the process of mass correction due to the special shape of the numerical solution. In this paper, we define θhn as

θhn,j = θhn (xj ) = −Ra(xj )

Mn

δt

= −Raj

Mn

δt

,

(4.8)

Please cite this article as: S. Zhao, X. Xiao, J. Zhao et al., A Petrov–Galerkin finite element method for simulating chemotaxis models on stationary surfaces, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.019.

6

S. Zhao, X. Xiao, J. Zhao et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

where the weighting function Ra(x) needs to meet the following requirement

∫

Ra(x)dσh

1= Γh

=

∫ ∑ N

Raj ψj (x)dσh

Γh j=1

=

N ∑

∫ Raj

Γh

j=1

(4.9) N 1∑

ψj (x)dσh =

3

Raj |Ωj |,

j=1

here, |Ωj | is the area of patch of point xj , that is the total area of the elements containing point xj . We note that node correction factor θhn,j = 0 if Mn = 0, which means the corrected mass is zero. And node correction factor θhn,j > 0 if Mn < 0, which does not break the positivity of numerical solutions. However, an inappropriate θhn,j may break the positivity of numerical solution when θhn,j < 0. We define the weighting function Ra(x) as the following Ra(xj ) =

3Rb(xj )

|Ωj |

N ∑

,

Rb(xj ) = 1.

(4.10)

i=1

The key is to give a suitable definition of the parameter Rb(x). The main idea of the mass correction is to allocate the corrected mass Mn to each vertex reasonably, an appropriate way is to allocate more corrected mass to some vertices with more mass. From the point of view, the parameter Rb (xj ) related to the correction factor can be defined as Rb (xj ) =

vertex-mass total mass

Uh∗,j ψj dσh

∫ = ∫

Γh

Γh

∑N

i=1

ψj dσh = ∫ ∗ i=1 Uh,i Γh ψi dσh Uh∗,j |Ωj | = ∑N . ∗ i=1 Uh,i |Ωi | Uh∗,j

∫

= ∑N

Uh∗,i ψi dσh 1 ∗ U 3 h,j

Γh

∑N 1

i=1

3

|Ωj |

(4.11)

Uh∗,i |Ωi |

Thus, the mass correction factor θhn can be denoted as

θhn (xj ) = −

3Uh∗,j Mn

δt

∑N

i=1

Uh∗,i |Ωi |

.

(4.12)

Obviously, the magnitude of corrected mass at each vertex is related to the shape of numerical solution. We can check that the scheme (4.5) does preserve the discrete conservation property. By Eq. (4.5)(b) and the correction factor (4.12), it is easy to show

∫ Γh

σ −

Uhn d h

∫

Uh dσh = δ t ∗

Γh

= δt

∫ Γh N ∑

θhn (x)dσh θhn,j

∫

j=1

= δt

N ∑

−

δt ∑N

j=1

Γh

ψj (x)dσh

3Uh∗,j Mn

∑N

∗ i=1 Uh,i |Ωi |

j=1

Uh∗,j |Ωj |

i=1

Uh∗,i |Ωi |

= −Mn ∑N

·

|Ωj |

(4.13)

3

= −Mn , and according to Eq. (4.5)(a), we get

∫ Γh

Uh∗ dσh −

∫ Γh

Uhn−1 dσh = Mn .

(4.14)

Please cite this article as: S. Zhao, X. Xiao, J. Zhao et al., A Petrov–Galerkin finite element method for simulating chemotaxis models on stationary surfaces, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.019.

S. Zhao, X. Xiao, J. Zhao et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

7

Thus, we have the discrete conservation property

∫ Γh

∫

Uhn dσh −

Γh

∫

Uhn−1 dσh =

Γh

Uh∗ dσh −

∫ Γh

Uhn−1 dσh − Mn

(4.15)

= M − M = 0. n

n

The results show that the key mass correction factor we proposed is valid. The proof of the positivity of the scheme (4.5) is shown in later paragraphs. 4.1. Decoupled approach and lumped mass In order to avoid solving nonlinear system, decoupled approach and operator recoveries techniques are considered. We give priority to solve (4.5)(c) M(Chn , Φ ) + D(Chn , Φ ) + Rc (Chn , Φ ) = G (Uh

n−1

, Φ)

to get in each moment. Then we obtain the value of ∇Γh Chn and ∆Γh Chn at each surface node by the recoveries techniques. Obviously, Eq. (4.5)(c) is equivalent to the following algebraic equation Chn

A + δ t(D + β A) Cn = δ tAUng −1 + ACn−1 ,

(

)

where, A = (aij ) is the mass matrix with the entries aij = (ψi , ψj )Γh , and D = (dij ) with dij = (∇Γh ψi , ∇Γh ψj )Γh is the stiffness matrix of the system (4.1), the components of Un , Cn and Ung −1 are the nodal values of Uhn , Chn and g(Uhn−1 )Uhn−1 , respectively. For the DMP, the lumped mass technique is employed to ensure the non-negative solution. Thus, the mass matrix A is approximated by the lumped mass matrix A with the entries

{∑ N aij =

0,

j=1

aij ,

i = j, i ̸ = j.

Thus, we have the non-negative solution Chn by solving the following algebraic equation A + δ t(D + β A) Cn = δ tAUng −1 + ACn−1 .

)

(

4.2. Operator recoveries and solving Uhn by Mizukami–Hughes method To solve (4.5)(a) and (4.5)(b), we require that the velocity field bk = ∇Γh Chn |ek of each element ek is defined as the barycentre (average) velocity. Hence, the vector of ∇Γh Chn should be well-defined at surface nodes by using the gradient recover technique. In this paper, the superconvergent patch recovery technique is employed on the surface mesh to obtain the accurate surface gradient information. This method was first proposed by Zienkiewicz et al. on planar mesh, and modified by Wei et al. for surface mesh, more details see [40,41]. ˜ ), we also need to recover ∆Γh Chn,j which is defined on surface nodes xj by For computing term Ru (Uh∗ , Chn , Φ

∆Γh Chn,j

=

1

αc

(

Chn,j − Chn,−j 1

δt

) +β

Chn,j

−

g(Uhn,−j 1 )Uhn,−j 1

.

(4.16)

By a piecewise linear interpolation, we have

∆Γh Chn =

N ∑

∆Γh Chn,j ψj .

j=1

Thus, we can substitute these recoveries into (4.5)(a) for solving Uh∗ . The algebraic equation is equivalent to (4.5)(a) is Au + δ t D + B(Cn ) + Au (Cn )

(

(

))

U∗ = δ tAu Unf −1 + Au Un−1 ,

(4.17)

˜j )Γh , diffusion matrix D = (∇Γh ψi , ∇Γh ψ ˜j )Γh where, Au = (˜ aij ) is the mass matrix with entries ˜ aij = (ψi , ψ (∇Γh ψi , ∇Γh ψj )Γh , convection matrix B(Cn ) with the entries

=

˜j )Γh . bij = (χ (Chn )∇Γh Chn · ∇Γh ψi , ψ The matrix Au (Cn ) is composed of the entries

˜j )Γh a∗ij = (∇Γh · (χ (Chn,i )∇Γh Chn,i )ψi , ψ ˜j )Γh = (χ (Chn,i )∆Γh Chn,i + χ ′ (Chn,i )|∇Γh Chn,i |2 )(ψi , ψ ˜j )Γh , = Li (ψi , ψ Please cite this article as: S. Zhao, X. Xiao, J. Zhao et al., A Petrov–Galerkin finite element method for simulating chemotaxis models on stationary surfaces, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.019.

8

S. Zhao, X. Xiao, J. Zhao et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 1. Definition of vertex zones and edge zones.

and the components of Unf −1 are the nodal values of f (Uhn−1 )Uhn−1 . The point of the Mizukami − Hughes method is the way of defining the constants γl on every element. Mizukami and Hughes require the constants γl to satisfy the following conditions in every element:

γl ≥ −

1 3

and

3 ∑

γl = 0,

l = 1, 2, 3,

(4.18)

i=1

and that the local convection matrix Bek (Cn ) with entries e

k ˜ blmk = (χ (Chn )bk · ∇Γh ψm , ψlk )ek

= (bk · ∇Γh ψmk )

∫

˜k dσ χ (Chn )ψ l ∫ 1 n k dσ , = χ (Ch )(bk · ∇Γh ψm )( + γl ) ek

3

l, m = 1, 2, 3.

ek e

k Obviously, the sign of the entries blmk is determined by bk · ∇Γh ψm . Thus Mizukami and Hughes construct the non-negative ek n type local convection matrix B (C ) by changing the vector bk in a direction that perpendiculars to ∇Γh Uhn |ek on each element without affecting the solution. Next, we introduce the way of selecting the velocity bk . We use v1 , v2 and v3 to denote the vertices of each element ek . Each element ek is divided into six zones by three lines crossing the barycentre and parallel to the edges of the element. Then, define the vertex zones VZl and the edge zones EZl as shown in Fig. 1. The boundary of two adjacent zones is classified into the respective vertex zone VZl . Let us first assume that bk · ∇Γh Uhn |ek ̸ = 0. If the velocity bk points into the vertex zone VZ1 as shown in Fig. 1, the following results can be obtained

bk ∈ VZ1 ⇔ bk · ∇Γh ψ1k > 0, bk · ∇Γh ψ2k ≤ 0, bk · ∇Γh ψ3k ≤ 0, then the non-negative type matrix Bek (Cn ) may be constructed by defining the constants

γ1 =

2 3

,

1

γ2 = γ3 = − . 3

However, it is impossible to obtain any matrix of non-negative type when the velocity bk points into the edge zone EZ1 , since bk ∈ EZ1 ⇔ bk · ∇Γh ψ1k < 0, bk · ∇Γh ψ2k > 0, bk · ∇Γh ψ3k > 0. Fortunately, the velocity bk can be replaced by ˜ bk = bk + ε w, where w is perpendicular to ∇Γh Uhn |ek and ε is a constant, then, we have

˜ bk · ∇Γh Uhn |ek = bk · ∇Γh Uhn |ek . Please cite this article as: S. Zhao, X. Xiao, J. Zhao et al., A Petrov–Galerkin finite element method for simulating chemotaxis models on stationary surfaces, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.019.

S. Zhao, X. Xiao, J. Zhao et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

9

Thus, the entries of the local convection matrix Bek (Cn ) can be replaced by following formulation when the velocity bk points into EZ1 e

k blmk = χ (Chn )(˜ bk · ∇Γh ψm )(

1

+ γl )

3

∫

dσ ,

l, m = 1, 2, 3.

ek

According to the above results, we can construct the non-negative type matrix Bek (Cn ) by changing bk such that the velocity points into other vertex zones VZl (l = 2, 3). Assuming that bk · ∇Γh Uhn |ek ̸ = 0, then Mizukami and Hughes show that the constants γl should be defined as follows:

˜ bk ∈ VZ2 and ˜ bk ∈ / VZ3

⇒

γ2 =

˜ bk ∈ / VZ2 and ˜ bk ∈ VZ3

⇒

γ3 =

˜ bk ∈ VZ2 and ˜ bk ∈ VZ3

2 3 2 3

1

, γ1 = γ3 = − , 3 1

, γ1 = γ2 = − , 1

bk · ∇ ψ

3

3|bk · ∇ ψ |

γ1 = − , γ2 =

⇒

(4.19)

3

k Γh 2 k Γh 1

, γ3 =

1 3

− γ2 .

For the case of bk · ∇Γh Uhn |ek = 0, it is no longer a matter of importance how to define the value of γl . For the sake of simplicity, we set

γ1 = −

1

and γ2 = γ3 =

3

1 6

.

In order to implement Mizukami–Hughes method, we need to determine the direction of the velocity ˜ bk . That is, we should know which zones (VZ2 or VZ3 ) the velocity ˜ bk points into. We denote d12 =

v2 − v1 , |v2 − v1 | bk

db =

|bk |

,

dw =

d13 = w

|w|

v3 − v1 , |v3 − v1 |

v3 − v2 , |v3 − v2 |

d23 =

.

The direction of the velocity ˜ bk can be determined by the following conditions

˜ bk ∈ VZ2 and ˜ bk ∈ / VZ3 ⇔ (d13 − dw ) · (db − dw ) ≤ 0 or (d13 + dw ) · (db + dw ) ≤ 0, ˜ bk ∈ / VZ2 and ˜ bk ∈ VZ3 ⇔ (d12 − dw ) · (db − dw ) ≤ 0 or (d12 + dw ) · (db + dw ) ≤ 0, ˜ bk ∈ VZ2 and ˜ bk ∈ VZ3

(4.20)

⇔ (d12 − dw ) · (d13 + dw ) ≥ 0 or (d13 − dw ) · (d12 + dw ) ≥ 0. It is not necessary to consider ε if the velocity flow bk ∈ VZ1 . When bk points into EZ1 , we set the ε such that if ˜ bk ∈ VZ2 and ˜ bk ∈ / VZ3 , then ε = ε0 ∈ R | max{˜ bk · ∇Γh ψ2k = (bk + ε0 w) · ∇Γh ψ2k } ,

{

}

if ˜ bk ∈ / VZ2 and ˜ bk ∈ VZ3 , then ε = ε0 ∈ R | max{˜ bk · ∇Γh ψ3k = (bk + ε0 w) · ∇Γh ψ3k } .

{

}

In the case of ˜ bk ∈ VZ2 and ˜ bk ∈ VZ3 , the second and third rows of the local convection matrix Bek (Cn ) are denoted as ek b2m

=χ

e

(Chn )

(

(bk + ε2 w) · ∇Γh ψ

k m

) 1 (

) 1

k b3m = χ (Chn ) (bk + ε3 w) · ∇Γh ψmk (

(

3 3

+ γ2 ) + γ3 )

∫ ∫ek

dσ , dσ ,

ek

respectively, and ε2 , ε3 are the undermined coefficients. We consider the following two different cases: (d12 − dw ) · (d23 + dw ) ≥ 0 or (d23 − dw ) · (d12 + dw ) ≥ 0

⇒

⎧ ⎨

bk · d12 ε2 = − , w·d } ⎩ε = {ε ∈ R12| max{˜ bk · ∇Γh ψ3k = (bk + ε0 w) · ∇Γh ψ3k } , 3 0

Please cite this article as: S. Zhao, X. Xiao, J. Zhao et al., A Petrov–Galerkin finite element method for simulating chemotaxis models on stationary surfaces, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.019.

10

S. Zhao, X. Xiao, J. Zhao et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

(d13 − dw ) · (d23 − dw ) ≤ 0 or (d23 + dw ) · (d13 + dw ) ≤ 0

⎧ { } ⎨ε2 = ε0 ∈ R | max{˜ bk · ∇Γh ψ2k = (bk + ε0 w) · ∇Γh ψ2k } , bk · d13 ⇒ ⎩ ε3 = − . w · d13 The above choices of the constant γl and the velocity flow bk ensure that the local convection matrix Bek (Cn ) is of

non-negative type, so we have the non-negative type matrix B(Cn ). But non-negative type matrix B(Cn ) is not enough for non-negative numerical solution Uhn , thus we approximate the matrix Au and Au (Cn ) with the lumped mass matrix Au and Au (Cn ), respectively. We can get the non-negative solution Uh∗ by solving the following algebraic equation Au + δ t D + B(Cn ) + Au (Cn )

(

(

))

U∗ = δ tAu Unf −1 + Au Un−1 .

Finally, the non-negative solution Uhn can be obtained by solving the following algebraic equation AUn = A(U∗ + δ t Θ n−1 ),

(4.21)

which corresponds to (4.5)(b). Remark 4.1. Mizukami–Hughes method is nonlinear method even for linear problem, since the constants γl and the velocity flow ˜ bk depend on ∇Γh Uhn . Thus, damped Newton iterative method is applied to solve nonlinear problem in each time step, and set the nonlinear iteration initial value Un0 = Un−1 . For improving the convergence of the nonlinear iterative, we can also consider an improved Mizukami–Hughes method that is proposed by Knobloch. The improved Mizukami– Hughes method choose the constants γl which continuously depend on the orientation of the flow velocity bk and the orientation of ∇Γh Uhn , more details refer to [26]. 5. The analysis of positivity preservation In this section, we discuss about positivity preservation of the proposed method. We only consider the minimum norm of Uh and Ch , since the final solutions are the liftings of Uh and Ch . Let u0 ≥ 0 and c0 ≥ 0, which leads to Uh0 ≥ 0 and Ch0 ≥ 0. It is noteworthy that the triangular mesh Γh is required to be Delaunay-type in the previous section, which ensures that the entries of the stiffness matrix D satisfy dij = (∇Γh ψi , ∇Γh ψj )Γh =

{

> 0, ≤ 0,

i = j, i ̸ = j.

In order to facilitate the analysis of the positive property, we introduce the following lemma: Lemma 5.1 ([32]). If S = (sij ) be a positive definite symmetric matrix with sii > 0 for all i = 1, 2, . . . , N and sij ≤ 0 for all i ̸ = j, then S −1 ≥ 0. The system (4.5) can be rewritten as the following matrix form

⎧ ∗ −1 ⎪ ⎨U = Su Fu , Un = (Su∗ )−1 Fu∗ , ⎪ ⎩ n C = Sc−1 Fc ,

(a) (5.1)

(b) (c)

where,

( ) ⎧ n−1 n n n−1 ⎪ ⎨Su = Au + δ t D + B(C ) + Au (C ) , Fu = δ tAu Uf + Au U , Su∗ = A, Fu∗ = A(U∗ + δ t Θ n ), ⎪ ⎩ Sc = A + δ t(D + β A), Fc = δ tAUng −1 + ACn−1 n−1

and we assume that U ≥ 0 and C are given by the following theorem.

n−1

(a) (b)

(5.2)

(c),

≥ 0 are given to the system. Then, sufficient conditions of positivity preservation

Theorem 5.1. Let the initial condition of the system (2.1) satisfies u0 ≥ 0 and c0 ≥ 0, and the triangular mesh Γh is of Delaunay type. Then the numerical solution Un and Cn of (5.1) satisfy that Un ≥ 0, Cn ≥ 0 for all n = 1, 2, . . . , M with time step size δ t ≤ min(T1 , T2 ), where

⎧ n−1 n−1 ⎪ ⎨ + ∞, (⏐ ⏐) f (Uh )Uh ≥ 0, n−1 ⏐ ⏐ Uh,i T1 = ⏐ ⏐ ⎪ , 0)⏐ , else, ⏐min( ⎩ i=min n−1 n−1 1,...,N ⏐ ⏐ f (Uh,i )Uh,i

(5.3)

Please cite this article as: S. Zhao, X. Xiao, J. Zhao et al., A Petrov–Galerkin finite element method for simulating chemotaxis models on stationary surfaces, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.019.

S. Zhao, X. Xiao, J. Zhao et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

11

and

⏐) (⏐ ⎧ ⎨ min ⏐⏐min( 1 , 0)⏐⏐ , ⏐ i=1,...,N ⏐ T2 = Li ⎩ + ∞,

Li ≥ 0, ∀i ∈ {1, 2, . . . , N },

(5.4)

else.

Proof. Firstly, we consider the positivity preservation of Cn . The coefficients αc , β are non-negative in the system (2.1), and matrices A, D are of non-negative type, which can ensure that Sc−1 ≥ 0

(5.5)

by Lemma 5.1. Since g(u) is non-negative and A is diagonal, we have Fc ≥ 0,

(5.6) n

which implies that the solution C ≥ 0. Next, we discuss the positivity preservation of U∗ . It is easy to get that Fu ≥ 0, when f (Uhn−1 )Uhn−1 ≥ 0. If f (Uhn−1 )Uhn−1 < 0, we must require the time step δ t to satisfy

⏐) (⏐ ⏐ ⏐ Uhn,−i 1 ⏐ ⏐ δ t ≤ min ⏐min( , 0)⏐ . n − 1 n − 1 i=1,...,N ⏐ ⏐ f (Uh,i )Uh,i

(5.7)

To ensure Su−1 ≥ 0,

(5.8)

the matrix A + δ tAu (Cn ) is required to be non-negative type. Consider matrix Au (Cn ) with entries ∗

{∑ N

aij =

j=1

0,

⏐∑ ⏐ N

˜j )Γh = sign(Li ) ⏐ L i (ψi , ψ

j=1

i ̸ = j,

⏐ ⏐

˜j )Γh ⏐ , L i (ψi , ψ

i=j,

we assume that Li ̸ = 0, thus a sufficient condition for (5.8) is

⏐) (⏐ ⏐ ⏐ 1 ⏐ δ t ≤ min ⏐min( , 0)⏐⏐ . i=1,...,N Li

(5.9)

On the other hand, it is obvious that the matrix Su satisfies that Su−1 ⩾ 0, if Li ⩾ 0 for all i = 1, 2, . . . , N. We can get the non-negative solution Uh∗ with time step δ t satisfying (5.7) and (5.9). Finally, we consider the positivity preservation of Un , that is the mass correction can ensure the positivity of numerical solution. For Eq. (5.2)(b), we need to check the condition Fu∗ ≥ 0, that is U∗ + δ t Θ n ≥ 0. Combining Eqs. (4.8), (4.10) and (4.11), the elements in vector U∗ + δ t Θ n can be denoted as

( 3Uh∗,i Mn = Uh∗,i 1 − Uh∗,i − ∑N ∗ i=1 Uh,i |Ωi |

Mn

)

∑N 1

∗ i=1 Uh,i |Ωi |

3

= Uh∗,i 1 − ∫

(

Γh

Mn Uh∗ dσh

(5.10)

) .

With (4.3) and (4.5)(b), we have

∫ Γh

Uh∗ dσh − Mn =

∫ Γh

Uhn−1 (1 + δ tf (Uhn−1 ))dσh .

(5.11)

Obviously, when time step δ t satisfies the condition (5.7) we can get

∫ Γh

Uh∗ dσh − Mn ≥ 0

(5.12)

which implies Fu∗ = A(U∗ + δ t Θ n ) ≥ 0. Then, according to (Su∗ )−1 = (A)−1 ≥ 0, we can derive Un ≥ 0. In particular, time step δ t does not need to satisfy the condition (5.7) when f (u) = 0. It is noteworthy that mass correction does not impose stricter restrictions on time steps δ t. Remark 5.1. Let all assumptions in Theorem 5.1 hold. Then condition δ t ≤ T2 is required for Keller–Segel model to obtain non-negative numerical solution, and δ t ≤ min(T1 , T2 ) is needed for pattern formation model. In the practical calculation, the condition (5.4) can be touched by trial calculation of gradually decreasing time step-size sequences.

Please cite this article as: S. Zhao, X. Xiao, J. Zhao et al., A Petrov–Galerkin finite element method for simulating chemotaxis models on stationary surfaces, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.019.

12

S. Zhao, X. Xiao, J. Zhao et al. / Computers and Mathematics with Applications xxx (xxxx) xxx Table 1 Example 6.1: Accuracy test at T = 0.1 for the proposed method. h

H1 − u

order

L2 − u

order

H1 − c

order

L2 − c

order

0.1236 0.0996 0.0809 0.0685

9.72e−3 7.62e−3 5.99e−3 4.98e−3

– 1.12 1.15 1.11

7.86e−3 5.10e−3 3.45e−3 2.50e−3

– 2.00 1.90 1.93

1.04e−2 7.97e−3 6.30e−3 5.23e−3

– 1.23 1.13 1.11

8.91e−4 5.88e−4 3.86e−4 2.79e−4

– 1.93 2.02 1.95

6. Numerical simulation In the section, we demonstrate the proposed method make it possible to perform numerical simulation of chemotaxis ( ) max |Unk − Unk−1 | ( n ) is less processes on surface. Convergence of nonlinear iteration is achieved when the tolerance ε = max Uk−1 than 1 × 10−5 . Example 6.1 (Accuracy Test). We solve the following problem

⎧ 2 ⎨ut = ∆Γ u − ∇Γ · (u∇Γ c) + u (1 − u) + f , c = ∆Γ c − c + u + g , ⎩t u|t =0 = u0 , c |t =0 = c0 , on a closed sphere

φ (x) = x2 + y2 + z 2 − 0.25 = 0. We choose suitable f and g such that the exact solutions are u(x, t) = et sin(π x), c(x, t) = et (x21 + x22 x3 ), and let δ t = O(h2 ) and terminal time T = 0.1. In Table 1, we present the convergence results for the proposed method. From the table, we can observe optimal rate of convergence. 6.1. Chemotaxis on closed surfaces We consider the Keller–Segel chemotaxis model posed on a closed sphere and a peanut-like surface. The purpose of the following numerical experiments is to show the effectiveness of the proposed method for simulating the blow-up behavior of the Chemotaxis system. Example 6.2 (Blow-up on a Sphere). Let us consider the Keller–Segel chemotaxis model (2.4) with χ (c) = 1

{

ut = ∆Γ u − ∇Γ · (u∇Γ c),

(6.1)

ct = ∆Γ c − c + u, on a closed sphere

φ (x) = x2 + y2 + z 2 − 0.25 = 0, and let the initial conditions are

{

u0 (x) = 600e−60(x

2 +y2 +(z −0.5)2 )

−30(x2 +y2 +(z −0.5)2 )

c0 (x) = 300e

, .

The mesh size h = 0.0323 and time step size δ t = 1 × 10−5 are employed. In Fig. 2, we plot the evolution of the peak at some time instance to more clearly show the shape of numerical solution, the numerical results are non-oscillatory and everywhere positive. Since the peaks of u0 (x) and c0 (x) are located at the same position, the peak position of Uhn does not change throughout the simulation. In Figs. 3(a)–3(c), we show the evolution of the maximum of Uhn with different mesh size. Obviously, the maximum amplitude of Uhn increases by a factor of 4 as mesh size decreases by 1/2, which indicates a blow up of Uhn in O(h2 ) fashion. The results are consistent with those in the literature [42]. Fig. 3(d) show the total mass for several mesh sizes h, we observe very good mass conservation.

Please cite this article as: S. Zhao, X. Xiao, J. Zhao et al., A Petrov–Galerkin finite element method for simulating chemotaxis models on stationary surfaces, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.019.

S. Zhao, X. Xiao, J. Zhao et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

13

Fig. 2. Example 6.2: Numerical results of the model (6.1) on the sphere.

Example 6.3 (Blow-up on a Peanut-like Surface). In this example, we consider the system (6.1) defined on a peanut-like surface

φ (x) = ((2x − 1)2 + 4y2 + 4z 2 )((2x + 1)2 + 4y2 + 4z 2 ) − 1.5 = 0 and the initial conditions are changed to

{ u0 (x) =

200e−20((x−0.4)

2 +y2 +(z −0.3)2 )

200e−20((x+0.4)

2 +y2 +(z −0.3)2 )

, ,

x ≥ 0, x < 0,

and c0 (x) = 100e−10(x

2 +y2 +(z −0.25)2 )

.

The function u0 (x) has two symmetric peaks on the surface, and the peak of c0 (x) is located at the symmetric point. We set the mesh size h = 0.0617 and time step size δ t = 1 × 10−5 . The numerical simulations are illustrated in Fig. 4, we notice that the species migrate towards the peak of the chemoattractant concentrations, and the numerical solution obtained is everywhere nonnegative. The evolution of maximum of Uhn and the total mass are shown in Fig. 5. 6.2. Pattern formations In this part, we use the proposed method to simulate the Pattern formation model defined on surface to investigate the applicability of the methods. The model describes motile cells of Escherichia coli forms the complex space–time pattern. Please cite this article as: S. Zhao, X. Xiao, J. Zhao et al., A Petrov–Galerkin finite element method for simulating chemotaxis models on stationary surfaces, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.019.

14

S. Zhao, X. Xiao, J. Zhao et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 3. Example 6.2: The evolution of max(Uhn ) and total mass.

Example 6.4. In order to observe the pattern formed, we use the larger sphere

φ (x) = x2 + y2 + z 2 − 4 = 0 as the computed domain. The parameters of the system (2.5) are set as αu = 0.0625, χ = 8.5 and β = 32, we have ut = 0.0625∆Γ u − 8.5∇Γ · (u∇Γ c) + u2 (1 − u),

{

(6.2)

ct = ∆Γ c − 32c + u,

and the initial conditions are given by

{ u0 (x) =

( √

1 + 0.5 cos2 π 1,

)

x2 + y2 + (z − 2)2 ,

√

x2 + y2 + (z − 2)2 < 1.5,

else,

1 . The mesh size h = 0.1225 and the time step δ t = 1 × 10−2 are employed. We show the temporal 32 evolution of the cell distribution in Fig. 6. The numerical simulation appears that the initial density of species placed in the local domain propagates into the whole domain in a perforated stripe type pattern. The results are similar to that obtained by other numerical methods [9,15,23,24]. and c0 (x) =

7. Conclusions We studied a modified Petrov–Galerkin finite element method for chemotaxis models defined on general static surfaces. The details about selecting velocity fields and weighting test functions were introduced for algorithm implementation. Firstly, the fully discretized form was obtained by a semi-implicit scheme that based on the gradient and Laplacian recoveries. Secondly, the lumped mass method was employed to guarantee the non-oscillatory numerical solution on the Delaunay-type triangulated surfaces. Furthermore, mass correction method was combined with Mizukami–Hughes method to satisfy the discrete mass conservation property. The selection criteria of time step size was given by stability analyses. We also find that the proposed mass correction equation does not pose any restrictions on time steps. Finally, Please cite this article as: S. Zhao, X. Xiao, J. Zhao et al., A Petrov–Galerkin finite element method for simulating chemotaxis models on stationary surfaces, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.019.

S. Zhao, X. Xiao, J. Zhao et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

15

Fig. 4. Example 6.3: Numerical results of (6.1) on the peanut-like surface.

Fig. 5. Example 6.3: The evolution of max(Uhn ) and total mass.

the numerical examples were performed to verify the effectiveness and stability of the proposed method. Thus the correction factor we introduced can also be used in other method. We will pay more attention to extend this method to the chemotaxis equations defined on evolving surface in future. CRediT authorship contribution statement Shubo Zhao: Methodology, Validation, Writing - original draft. Xufeng Xiao: Software, Visualization. Jianping Zhao: Visualization, Writing - review & editing. Xinlong Feng: Methodology, Funding acquisition, Supervision, Validation. Please cite this article as: S. Zhao, X. Xiao, J. Zhao et al., A Petrov–Galerkin finite element method for simulating chemotaxis models on stationary surfaces, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.019.

16

S. Zhao, X. Xiao, J. Zhao et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 6. Example 6.4: Numerical results of (6.2) on the sphere.

Acknowledgments The authors would like to thank the editor and referees for their valuable comments and suggestions which helped us to improve the results of this article. References [1] E. Keller, L. Segel, Initiation of slime mold aggregation viewed as an instability, J. Theoret. Biol. 26 (1970) 399–415. [2] E. Keller, L. Segel, Model for chemotaxis, J. Theor. Biol. 30 (1971) 225–234. [3] T. Nagai, Blowup of nonradial solutions to parabolic–elliptic systems modeling chemotaxis in two-dimensional domains, J. Inequal. Appl. 6 (2001) 37–55. [4] D. Horstmann, G. Wang, Blow-up in a chemotaxis model without symmetry assumptions, European J. Appl. Math. 12 (2001) 159–177. [5] L. Corriasa, B. Perthame, Asymptotic decay for the solutions of the parabolic Keller–Segel chemotaxis system in critical spaces, Math. Comput. Modelling 47 (2008) 755–764. [6] D. Horstmann, M. Winkler, Boundedness vs. blow-up in a chemotaxis system, J. Differential Equations 215 (2005) 52–107. [7] N. Saito, Conservative upwind finite-element method for a simplified Keller–Segel system modelling chemotaxis, IMA J. Numer. Anal. 27 (2) (2007) 332–365. [8] J. Zhang, J. Zhu, R. Zhang, Characteristic splitting mixed finite element analysis of Keller–Segel chemotaxis models, Appl. Math. Comput. 278 (2016) 33–44. [9] R. Strehl, A. Sokolov, D. Kuzmin, et al., A flux-corrected finite element method for chemotaxis problems, Comput. Methods Appl. Math. 10 (2) (2010) 219–232. [10] R. Strehl, A. Sokolov, D. Kuzmin, et al., A positivity-preserving finite element method for chemotaxis problems in 3D, J. Comput. Appl. Math. 239 (2013) 290–303. [11] R. Strehl, A. Sokolov, S. Turek, Efficient, accurate and flexible finite element solvers for chemotaxis problems, Comput. Math. Appl. 64 (3) (2012) 175–189. [12] X. Huang, X. Xiao, J. Zhao, X. Feng, An efficient operator-splitting FEM-FCT algorithm for 3D chemotaxis models, Eng. Comput. (2019) http://dx.doi.org/10.1007/s00366-019-00771-8. [13] M. Akhmouch, M.B. Amine, A corrected decoupled scheme for chemotaxis models, J. Comput. Appl. Math. 323 (2017) 36–52.

Please cite this article as: S. Zhao, X. Xiao, J. Zhao et al., A Petrov–Galerkin finite element method for simulating chemotaxis models on stationary surfaces, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.019.

S. Zhao, X. Xiao, J. Zhao et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

17

[14] B. Andreianov, M. Bendahmane, M. Saad, Finite volume methods for degenerate chemotaxis model, J. Comput. Appl. Math. 235 (14) (2011) 4015–4031. [15] A. Chertock, A. Kurganov, A second-order positivity preserving central-upwind scheme for chemotaxis and haptotaxis models, Numer. Math. 111 (2) (2008) 169–205. [16] F. Filbet, A finite volume scheme for the Patlak–Keller–Segel chemotaxis model, Numer. Math. 104 (4) (2006) 457–488. [17] Y. Epshteyn, Discontinuous Galerkin methods for the chemotaxis and haptotaxis models, J. Comput. Appl. Math. 224 (1) (2009) 168–181. [18] Y. Epshteyn, A. Izmirlioglu, Fully discrete analysis of a discontinuous finite element method for the Keller–Segel chemotaxis model, J. Sci. Comput. 40 (1–3) (2009) 211–256. [19] Y. Epshteyn, A. Kurganov, New interior penalty discontinuous Galerkin methods for the Keller–Segel chemotaxis model, SIAM J. Numer. Anal. 47 (1) (2008) 386–408. [20] X.H. Li, C.W. Shu, Y. Yang, Local discontinuous Galerkin method for the Keller–Segel chemotaxis model, J. Sci. Comput. 73 (2–3) (2017) 943–967. [21] R. Zhang, J. Zhu, A.F.D. Loula, et al., Operator splitting combined with positivitypreserving discontinuous Galerkin method for the chemotaxis model, J. Comput. Appl. Math. 302 (2016) 312–326. [22] C.M. Elliott, B. Stinner, C. Venkataraman, Modelling cell motility and chemotaxis with evolving surface finite elements, J. R. Soc. Interface 9 (2012) 3027–3044. [23] A. Sokolov, R. Strehl, S. Turek, Numerical simulation of chemotaxis models on stationary surfaces, Discrete Contin. Dyn. Syst. Ser. B 18 (10) (2013) 2689–2704. [24] X. Xiao, X. Feng, Y. He, Numerical simulations for the chemotaxis models on surfaces via a novel characteristic finite element method, Comput. Math. Appl. 78 (1) (2019) 20–34. [25] A. Mizukami, T.J.R. Hughes, A Petrov–Galerkin finite element method for convectiondominated flows: an accurate upwinding technique for satisfying the maximum principle, Comput. Methods Appl. Mech. Engrg. 50 (1985) 181–193. [26] P. Knobloch, Improvements of the Mizukami–Hughes method for convection–diffusion equations, Comput. Methods Appl. Mech. Engrg. 196 (2006) 579–594. [27] G. Dziuk, C.M. Elliott, Finite element methods for surface PDEs, Acta Numer. 22 (2013) 289–396. [28] G. Dziuk, C.M. Elliott, Surface finite elements for parabolic equations, J. Comput. Math. 25 (2007) 385–407. [29] X. Xiao, K. Wang, X. Feng, A lifted local Galerkin method for solving the reaction–diffusion equations on implicit surfaces, Comput. Phys. Comm. 231 (2018) 107–113. [30] X. Xiao, X. Feng, Z. Li, A gradient recovery-based adaptive finite element method for convection–diffusion-reaction equations on surfaces, Internat. J. Numer. Methods Engrg. 120 (7) (2019) 901–917. [31] X. Xiao, J. Zhao, X. Feng, A layers capturing type H-adaptive finite element method for convection–diffusion-reaction equations on surfaces, Comput. Methods Appl. Mech. Engrg. 361 (2020) 112792. [32] V. Thomée, Galerkin Finite Element Methods for Parabolic Problems, Springer-Verlag, Berlin, 1984. [33] X. Xiao, X. Feng, J. Yuan, The lumped mass finite element method for surface parabolic problems: error estimates and maximum principle, Comput. Math. Appl. 76 (3) (2018) 488–507. [34] S. Zhao, X. Xiao, Z. Tan, X. Feng, Two types of spurious oscillations at layers diminishing methods for convection–diffusion-reaction equations on surface, Numer. Heat Transfer A 74 (7) (2018) 1387–1404. [35] J. Rubinstein, P. Sternberg, Nonlocal reaction–diffusion equations and nucleation, IMA J. Appl. Math. 48 (3) (1992) 249–264. [36] A. Salih, S.G. Moulic, A mass conservation scheme for level set method applied to multiphase incompressible flows, Int. J. Comput. Methods Eng. Sci. Mech. 14 (4) (2013) 271–289. [37] B. Perthame, Transport Equations in Biology, in: Frontiers in Mathematics, 2007. [38] E. Burman, A. Ern, Stabilization Galerkin approximatiion of convection–diffusion-reaction equation: discrete maximum principle and convergence, Math. Comp. 74 (2005) 1637–1652. [39] M.A. Olshanskii, A. Reusken, X. Xu, A stabilized finite element method for advection-diffusion equations on surfaces, IMA J. Numer. Anal. 34 (2014) 732–758. [40] O.C. Zienkiewicz, J.Z. Zhu, J. Wu, Superconvergent patch recovery techniques-some further tests, Commun. Numer. Methods Eng. 9 (3) (1993) 251–258. [41] H. Wei, L. Chen, Y. Huang, Superconvergence and gradient recovery of linear finite elements for the Laplace–Beltrami operator on general surfaces, SIAM J. Numer. Anal. 48 (5) (2010) 1920–1943. [42] J. Liu, L. Wang, Z. Zhou, Positivity-preserving and asymptotic preserving method for 2D Keller–Segal equations, Math. Comp. 87 (311) (2017) 1165–1189.

Please cite this article as: S. Zhao, X. Xiao, J. Zhao et al., A Petrov–Galerkin finite element method for simulating chemotaxis models on stationary surfaces, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.019.