Finite Element numerical schemes for a chemo-attraction and consumption model^†^†thanks: This research work has been supported by Proyecto PGC2018-098308-B-I00, financed by FEDER / Ministerio de Ciencia e Innovación - Agencia Estatal de Investigación, Spain.

F. Guillén-González Dpto. Ecuaciones Diferenciales y Análisis Numérico and IMUS, Universidad de Sevilla, Facultad de Matemáticas, C/ Tarfia, S/N, 41012 Sevilla (SPAIN). Email: [email protected] G. Tierra Department of Mathematics, University of North Texas. Email: [email protected]

Abstract

This work is devoted to design and study efficient and accurate numerical schemes to approximate a chemo-attraction model with consumption effects, which is a nonlinear parabolic system for two variables; the cell density and the concentration of the chemical signal that the cell feel attracted to.

We present several finite element schemes to approximate the system, detailing the main properties of each of them, such as conservation of cells, energy-stability and approximated positivity. Moreover, we carry out several numerical simulations to study the efficiency of each of the schemes and to compare them with others classical schemes.

2010 Mathematics Subject Classification. 35K51, 35Q92, 65M12, 65M60, 92C17.

Keywords: Chemo-Attraction and consumption model, finite elements, energy-stability, approximated positivity.

1 Introduction

Chemotaxis can be defined as the orientation or movement of an organism or cell in relation to chemical agents. This movement can be towards a higher (attractive) or lower (repulsive) concentration of the chemical substance. At the same time, the presence of living organisms can produce or consume the chemical substance, producing nontrivial dynamics of the living organisms and chemical substances. In 1971, Keller and Segel [16] introduced the model (1) that can be considered the first realistic attempt to capture the chemotactic response of bacteria towards chemical agents in a bounded spatial domain $\Omega\subset\mathbb{R}^{d}$ ( $d=1,2,3$ ) during a time interval $[0,T]$ , where the cell population density $u(\boldsymbol{x},t)$ moves towards the concentration of the chemical substance $v(\boldsymbol{x},t)$ , which is produced by the cell population with a rate $\mu>0$ :

\left\{\begin{array}[]{rcl}u_{t}-\Delta u+\chi\nabla\cdot(u\nabla v)&=&0\,,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\alpha\,v_{t}-\Delta v+v-\mu\,u&=&0\,,\end{array}\right.

(1)

where $\chi>0$ denotes the chemo-sensitivity parameter, $\alpha\in\{0,1\}$ determines the character of the chemical equation, being parabolic when $\alpha=1$ and elliptic when $\alpha=0$ . Since then, many models following the same spirit have been proposed and studied mathematically (check [3, 14, 15] for reviews on the development of this topic during the last years).

In this work we focus on developing numerical schemes for a system where the cell population is attracted by the chemical substance, which is consumed by the cell population with a rate proportional to the amount of living organisms:

\left\{\begin{array}[]{rcl}u_{t}-\Delta u+\chi\,\nabla\cdot(u\nabla v)&=&0\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr v_{t}-\Delta v+\mu\,u\,v&=&0\,.\end{array}\right.

(2)

There are several works that have focused on studying the analytical properties of model (2). In [20] the corresponding $d$ -dimensional problem was studied for $d\geq 2$ , proving global (in time) regularity and uniqueness whether the following criterium holds

0<\chi\leq\frac{1}{6(d+1)\|v_{0}\|_{L^{\infty}(\Omega)}}.

For $3$ -dimensional domains and arbitrarily large initial data, in [21] it is showed that this type of system admits at least one global weak solution. Moreover, it is asserted that such solutions at least eventually (i.e., for large enough times) become bounded and smooth, and that they approach the unique relevant constant steady state in the large time limit.

In recent times several groups have focused on numerical analysis for this type of attractive chemotaxis models with consumption/production of chemical substance and many relevant works have been produced. It is not our intention to provide a detailed review of all the works that have been produced in recent years, but rather to provide the reader with some interesting works that are related with the work that we present in this text. In [6] the authors investigate nonnegativity of exact and Finite Element (FE) numerical solutions to a generalized Keller-Segel model, under certain standard assumptions. Marrocco presented in [17] a new formulation of the Keller-Segel system, based on the introduction of a new variable and he approximated this new system via a mixed FE technique. Saito in [19] focused on presenting an error analysis of an approximation for the Keller-Segel system using a semi-implicit time discretization with a time-increment control and Baba-Tabata’s conservative upwind FE approximation [1], that allows to show the positivity and mass conservation properties of the scheme.

Several works have also focused on numerical schemes to approximate the simplified Keller-Segel system, where the parabolic equation for the concentration of the chemical substance is replaced by an elliptic one (taking $\alpha=0$ in (1)), arriving at a parabolic-elliptic system. Filbet presented in [9] a well-posed Finite Volume scheme to discretize the simplified Keller-Segel, providing a priori estimates and convergence. In [25] a Finite Volume approximation of the simplified Keller-Segel system is also considered, presenting a linear scheme that satisfies both positivity and mass conservation, deriving some inequalities on the discrete free energy and under some assumptions they establish error estimates in $L^{p}$ norm with a suitable $p>2$ for the $2$ -dimensional case. A simplified Keller-Segel system with additional cross diffusion is presented in [4]. The main feature of this model is that there exists a new entropy functional yielding gradient estimates for the cell density and chemical concentration. The authors also present in [4] a Finite Volume scheme that satisfy positivity preservation, mass conservation, entropy stability, and (under additional assumptions) entropy dissipation. Moreover, the existence of a discrete solution and its numerical convergence to the continuous solution is proved.

There are more related models that have been studied recently. For instance, one related type of models are the ones that focus on repulsive chemotaxis systems with the cell population producing chemical substance, that is, a system like (1) with $\chi<0$ . We refer the reader to [11, 12] (and the references therein) where the authors focused on studying unconditionally energy stable and mass-conservative FE numerical schemes, by introducing the gradient of the chemical concentration variable, for chemo-repulsive systems with quadratic ( $-\mu u^{2}$ ) and linear production terms ( $-\mu u$ ) in $\eqref{eq:KellerSegel}_{2}$ , respectively.

On the other hand, it has been experimentally observed [22] that the chemotactic motion in liquid environments affects substantially the migration of cells and this fact has increased the interest of studying the coupling of chemotaxis systems with the Navier-Stokes equations. For instance, Winkler in [23] has studied the $d$ -dimensional problem ( $d=2,3$ ) of attractive chemotaxis models with consumption of chemical substance, showing that under suitable regularity assumptions on the initial data, the chemotaxis-Navier-Stokes system admits a unique global classical solution ( $d=2$ ) and the simplified chemotaxis-Stokes system possesses at least one global weak solution ( $d=3$ ). Moreover, in [7] the authors construct numerical approximations for the same type of system. The presented approximations are based on using the Finite Element method, obtaining optimal error estimates and convergence towards regular solutions. Finally, we would like to mention another related work [13], where Finite Elements together with singular potentials have been used in the context of the Cahn-Hilliard equation to achieve energy-stable numerical schemes that satisfy approximated positivity properties.

This work is organized as follows. In Section 2 we present the attractive chemotaxis with consumption model that we have considered, its main properties and a reformulation that will allow us to design numerical schemes satisfying some energy laws, getting in particular an energy stable scheme in one-dimensional domains. The numerical schemes are developed and studied in Section 3. In Section 4 we report some numerical experiments that we have performed to study the efficiency and the accuracy of the schemes and to compare them with others classical schemes. Finally, the conclusions of our work are presented in Section 5.

2 The model

In this work, we consider an attractive-consumption chemotaxis model in a bounded domain $\Omega\subset\mathbb{R}^{d}$ $(d=1,2,3)$ given by the following parabolic system of PDEs:

\left\{\begin{array}[]{rll}u_{t}-\Delta u+\chi\,\nabla\cdot(u\nabla v)=0&\quad\mbox{in}\ \Omega\,,&t>0\,,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr v_{t}-\Delta v+\mu\,u\,v=0&\quad\mbox{in}\ \Omega\,,&t>0\,,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle\nabla u\cdot{\boldsymbol{n}}=\nabla v\cdot{\boldsymbol{n}}=0&\quad\mbox{on}\ \partial\Omega\,,&t>0\,,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr u(0,\boldsymbol{x})=u_{0}(\boldsymbol{x})\geq 0,\ v(0,\boldsymbol{x})=v_{0}(\boldsymbol{x})>0&\quad\mbox{in}\ \Omega\,,&\end{array}\right.

(3)

where $u(t,\boldsymbol{x})\geq 0$ denotes the cell population density and $v(t,\boldsymbol{x})>0$ denotes the concentration of the chemical substance that the cell feel attracted to and ${\boldsymbol{n}}$ denotes the outward normal vector to $\partial\Omega$ . Moreover, $(u_{0},v_{0})$ represents the initial density and concentration, $\chi>0$ is the chemo-sensitivity coefficient and $\mu>0$ is the consumption one.

2.1 Properties of the model

It is known that any regular enough solution $(u,v)$ of problem (3) satisfies the following properties:

1.

Positivity of $u$ and strictly positivity of $v$ ([20])
If $u_{0}\geq 0$ in $\Omega$ then $u(t,\cdot)\geq 0$ in $\Omega$ for any $t>0$ . Assuming $u\in L^{\infty}((0,\infty)\times\Omega)$ and $v_{0}\geq v_{0}^{min}>0$ in $\Omega$ then $v(t,\cdot)>0$ in $\Omega$ for any $t>0$ . In fact, one has the lower bound

$v(t,\cdot)\geq v_{0}^{min}\exp(-\mu\,t\,\|u\|_{L^{\infty}((0,t)\times\Omega)}).$
2.

Maximum principle for $v$ ([8])
One has $0\leq v(t,\cdot)\leq\|v_{0}\|_{L^{\infty}(\Omega)}$ in $\Omega$ for any $t>0$ . In fact, $\|v(t,\cdot)\|_{L^{\infty}(\Omega)}$ is a non-increasing function.

Cell density conservation. Integrating equation (3)₁,

\frac{d}{dt}\left(\int_{\Omega}u(t,\cdot)\right)=0\,,\quad\mbox{that is,}\quad\int_{\Omega}u(t,\cdot)=\int_{\Omega}u_{0}\,,\quad\forall\,t>0\,.

4.

Weak regularity for $v$ . Testing equation (3)₂ by $v$ , one has

$\frac{d}{dt}\|v(t,\cdot)\|_{L^{2}(\Omega)}^{2}+\|\nabla v(t,\cdot)\|_{L^{2}(\Omega)}^{2}\leq 0\,.$ (4)

Estimate of a singular functional. Assuming $u(t,x)>0$ for all $(t,x)$ , one has the time differential inequality

\frac{d}{dt}\left(\int_{\Omega}G(u)d\boldsymbol{x}\right)+\frac{1}{2}\int_{\Omega}\frac{1}{u^{2}}|\nabla u|^{2}d\boldsymbol{x}\,\leq\,\frac{\chi^{2}}{2}\int_{\Omega}|\nabla v|^{2}d\boldsymbol{x}\,,

(5)

where

G(u)=-\log(u)+C

is a convex function and the right hand side of (5) belongs to $L^{1}(0,+\infty)$ owing to (4). Relation (5) is derived testing equation (3)₁ by $G^{\prime}(u)=-1/u$ ,

\frac{d}{dt}\left(\int_{\Omega}G(u)d\boldsymbol{x}\right)+\int_{\Omega}\frac{1}{u^{2}}|\nabla u|^{2}d\boldsymbol{x}\,=\,-\chi\int_{\Omega}\frac{1}{u}\nabla u\cdot\nabla v\,d\boldsymbol{x}\,,

hence (5) holds by using Hölder inequality. This differential inequality (5) is analogous to the one derived in [24] for a Keller-Segel system with singular sensitivity.

Energy law [23]. Assuming $u(t,x)>0$ for all $(t,x)$ , one has the energy law:

\displaystyle\frac{d}{dt}E(u,v)\,+\mu\,D_{1}(u)\,+\chi\,D_{2}(v)+\mu\,\chi\,D_{3}(u,v)\,=\,\chi R(v)\,,

(6)

where

\begin{array}[]{rcl}E(u,v)&:=&\displaystyle\frac{\mu}{4}\int_{\Omega}F(u)d\boldsymbol{x}+\frac{\chi}{2}\int_{\Omega}|\nabla(\sqrt{v})|^{2}d\boldsymbol{x}\,,\quad\hbox{with $F^{\prime}(u)=\log(u)$,}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr D_{1}(u)&:=&\displaystyle\frac{1}{4}\int_{\Omega}\frac{1}{u}|\nabla u|^{2}d\boldsymbol{x}\,=\,\int_{\Omega}|\nabla(\sqrt{u})|^{2}d\boldsymbol{x}\geq 0\,,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr D_{2}(v)&:=&\displaystyle\int_{\Omega}(\Delta(\sqrt{v}))^{2}d\boldsymbol{x}+\ \frac{1}{3}\int_{\Omega}\frac{1}{v}|\nabla(\sqrt{v})|^{4}d\boldsymbol{x}\geq 0\,,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr D_{3}(u,v)&:=&\displaystyle\frac{1}{2}\int_{\Omega}u\,|\nabla(\sqrt{v})|^{2}d\boldsymbol{x}\geq 0\,,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr R(v)&:=&-\displaystyle\frac{2\chi}{3}\int_{\Omega}\frac{1}{\sqrt{v}}\Big{[}\nabla\cdot\Big{(}\big{|}\nabla(\sqrt{v})\big{|}^{2}\nabla(\sqrt{v})\Big{)}-3\big{(}\nabla(\sqrt{v})\big{)}^{t}\Big{(}\nabla\big{(}\nabla(\sqrt{v})\big{)}\Big{)}^{t}\nabla(\sqrt{v})\Big{]}d\boldsymbol{x}\,,\end{array}

Energy law (6) can be proved following the same ideas presented in Theorem 2.2 below. Moreover, in the particular case of $1$ -dimensional domains, the energy $E(u,v)$ is dissipative due to the right-hand term of (6) vanishes. In fact,

R(v)\,=\,-\frac{2\chi}{3}\int_{\Omega}\frac{1}{\sqrt{v}}\Big{[}\Big{(}\big{(}(\sqrt{v})_{x}\big{)}^{3}\Big{)}_{x}-3\big{(}(\sqrt{v})_{x}\big{)}^{2}(\sqrt{v})_{xx}\Big{]}dx\,=\,0\,.

On the contrary, in higher dimensions it is not clear the sign of $R(v)$ , preventing the possibility of obtaining a dissipative energy law without introducing constraints on the physical parameters.

Remark 2.1.

Notice that functional $G(u)=-\log(u)+C$ in the inequality (5) is more singular (for $u=0$ ) than the energy potential $F(u)=u\log(u)-u+1$ in (6). These singularities will be crucial to prove the approximate positivity of some of the numerical schemes presented in this work. A similar idea has been considered in [13] in the context of the Cahn-Hilliard equation.

2.2 Reformulation of the problem. The $(u,v,{\boldsymbol{\sigma}})$ problem

In order to develop a numerical scheme satisfying a discrete version of the energy estimate (6), we need to reformulate problem (3). The idea is to rewrite the $v$ -equation (3)₂ multiplying it by $w^{\prime}(v)=1/(2\sqrt{v})$ (hence $w(v)=\sqrt{v})$ , as follows

(\sqrt{v})_{t}\,-\,\frac{\Delta((\sqrt{v})^{2})}{2\sqrt{v}}\,+\,\frac{\mu}{2}u\sqrt{v}\,=\,0\,.

Then, by introducing the notation $w:=\sqrt{v}>0$ (due to the positivity of $v$ ) we obtain the PDE

w_{t}\,-\,\frac{1}{w}|\nabla w|^{2}\,-\,\Delta w\,+\,\frac{\mu}{2}uw\,=\,0\,.

(7)

Now we can take the gradient with respect to $\boldsymbol{x}$ of (7) to obtain:

(\nabla w)_{t}\,+\,\frac{1}{w^{2}}|\nabla w|^{2}\nabla w\,-\,\frac{1}{w}\nabla\big{(}|\nabla w|^{2}\big{)}\,-\,\nabla\Delta w\,+\,\frac{\mu}{2}u\nabla w\,+\,\frac{\mu}{2}w\,\nabla u\,=\,0\,.

(8)

Introducing a new unknown ${\boldsymbol{\sigma}}:=\nabla w=\nabla(\sqrt{v})$ we can reformulate three of the terms in (8):

•

Use the definition of ${\boldsymbol{\sigma}}$ to rewrite the term $\displaystyle\frac{1}{w^{2}}|\nabla w|^{2}\nabla w$ as

\frac{1}{w^{2}}|\nabla w|^{2}\nabla w\,=\,\frac{1}{3v}|{\boldsymbol{\sigma}}|^{2}{\boldsymbol{\sigma}}+\frac{2}{3v}|{\boldsymbol{\sigma}}|^{2}\nabla(\sqrt{v})\,.

•

Use $\sigma$ and a truncation of $u$ to rewrite term $\displaystyle\frac{\mu}{2}\nabla w\,u$ as $\displaystyle\frac{\mu}{2}{\boldsymbol{\sigma}}\,u_{+}$ , where $u_{+}$ denotes the positive part of $u$ ( $u_{+}(\boldsymbol{x})=\max\{u(\boldsymbol{x}),0\}$ ).
•

Rewrite the term $\displaystyle\frac{\mu}{2}w\,\nabla u$ as $\displaystyle\frac{\mu}{2}\sqrt{v}\,\,\nabla u$ .

Using these considerations, we formally arrive at the following $(u,v,{\boldsymbol{\sigma}})$ reformulation of the problem (3):

\left\{\begin{array}[]{rll}u_{t}-\Delta u+2\chi\nabla\cdot(u\,\sqrt{v}\,{\boldsymbol{\sigma}})=0&\quad\mbox{in}\ \Omega\,,&t>0\,,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle v_{t}-\Delta v+\mu\,u\,v=0&\quad\mbox{in}\ \Omega\,,&t>0\,,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle{\boldsymbol{\sigma}}_{t}-\frac{2}{\sqrt{v}}(\nabla{\boldsymbol{\sigma}})^{t}{\boldsymbol{\sigma}}+\frac{1}{3v}|{\boldsymbol{\sigma}}|^{2}{\boldsymbol{\sigma}}+\frac{2}{3v}|{\boldsymbol{\sigma}}|^{2}\nabla(\sqrt{v})&&\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle-\nabla(\nabla\cdot{\boldsymbol{\sigma}})+\mbox{rot}(\mbox{rot}\,{\boldsymbol{\sigma}})+\frac{\mu}{2}{\boldsymbol{\sigma}}\,u_{+}+\frac{\mu}{2}\sqrt{v}\,\,\nabla u=0&\quad\mbox{in}\ \Omega\,,&t>0\,,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle\nabla u\cdot{\boldsymbol{n}}=\nabla v\cdot{\boldsymbol{n}}={\boldsymbol{\sigma}}\cdot{\boldsymbol{n}}=(\mbox{rot}\,{\boldsymbol{\sigma}}\times{\boldsymbol{n}})\big{|}_{tang}=0&\quad\mbox{on}\ \partial\Omega\,,&t>0\,,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle u(\boldsymbol{x},0)=u_{0}(\boldsymbol{x})\geq 0\,,v(0,\boldsymbol{x})=v_{0}(\boldsymbol{x})>0\,,{\boldsymbol{\sigma}}(\boldsymbol{x},0)={\boldsymbol{\sigma}}_{0}(\boldsymbol{x})&\quad\mbox{in}\ \Omega\,,&\end{array}\right.

(9)

where $\mbox{rot}\,{\boldsymbol{\sigma}}$ denotes the well-known rotational operator (also called curl) which is scalar for 2D domains and vectorial for 3D ones. We have introduced the term $\mbox{rot}(\mbox{rot}\,{\boldsymbol{\sigma}})(=\mbox{rot}(\mbox{rot}\,\nabla v)=0)$ as in [12] to complete the ${\boldsymbol{H}}^{1}_{{\boldsymbol{\sigma}}}$ -norm of ${\boldsymbol{\sigma}}$ , where:

{\boldsymbol{H}}^{1}_{{\boldsymbol{\sigma}}}(\Omega)\,:=\,\{{\boldsymbol{\sigma}}\in{\boldsymbol{H}}^{1}(\Omega):{\boldsymbol{\sigma}}\cdot{\boldsymbol{n}}=0\,\mbox{ on }\partial\Omega\}\,,

and

\|{\boldsymbol{\sigma}}\|_{{\boldsymbol{H}}_{\boldsymbol{\sigma}}^{1}}^{2}\,:=\,\|{\boldsymbol{\sigma}}\|_{{\boldsymbol{L}}^{2}}^{2}+\|\mbox{rot}\,{\boldsymbol{\sigma}}\|_{{\boldsymbol{L}}^{2}}^{2}+\|\nabla\cdot{\boldsymbol{\sigma}}\|_{L^{2}}^{2}\,,

with $\|\cdot\|_{{\boldsymbol{H}}_{\boldsymbol{\sigma}}^{1}}$ being equivalent to $\|\cdot\|_{{\boldsymbol{H}}^{1}}$ .

Theorem 2.2.

Any regular enough solution $(u,v,{\boldsymbol{\sigma}})$ of (9) satisfy the following energy law:

\displaystyle\frac{d}{dt}{E}(u,{\boldsymbol{\sigma}})\,+\mu\,{D}_{1}(u)\,+\chi\,{D}_{2}(v,{\boldsymbol{\sigma}})\,+\chi\mu\,{D}_{3}(u,{\boldsymbol{\sigma}})\,=\,R(v,{\boldsymbol{\sigma}})\,,

(10)

where

\begin{array}[]{rclrcl}{E}(u,{\boldsymbol{\sigma}})&:=&\displaystyle\frac{\mu}{4}\int_{\Omega}F(u)d\boldsymbol{x}\,+\,\frac{\chi}{2}\int_{\Omega}|{\boldsymbol{\sigma}}|^{2}d\boldsymbol{x}\,,&{D}_{1}(u)&:=&\displaystyle\frac{1}{4}\int_{\Omega}\frac{1}{u}|\nabla u|^{2}\,d\boldsymbol{x}\,,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr{D}_{2}(v,{\boldsymbol{\sigma}})&:=&\displaystyle\int_{\Omega}\Big{(}|\nabla\cdot{\boldsymbol{\sigma}}|^{2}+|\mbox{rot}\,{\boldsymbol{\sigma}}|^{2}\Big{)}d\boldsymbol{x}+\frac{1}{3}\int_{\Omega}\frac{1}{v}|{\boldsymbol{\sigma}}|^{4}d\boldsymbol{x}\,,&{D}_{3}(u,{\boldsymbol{\sigma}})&:=&\displaystyle\frac{1}{2}\int_{\Omega}u_{+}\,|{\boldsymbol{\sigma}}|^{2}d\boldsymbol{x}\,.\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr R(v,{\boldsymbol{\sigma}})&:=&-\displaystyle\frac{2\chi}{3}\int_{\Omega}\frac{1}{\sqrt{v}}\Big{(}\nabla\cdot(|{\boldsymbol{\sigma}}|^{2}{\boldsymbol{\sigma}})-3{\boldsymbol{\sigma}}^{t}(\nabla{\boldsymbol{\sigma}})^{t}{\boldsymbol{\sigma}}\Big{)}d\boldsymbol{x}\,,&&\end{array}

(11)

with ${D}_{1}(u),\,{D}_{2}(v,\sigma),\,{D}_{3}(u,\sigma)\geq 0$ .

Proof.

The key argument of this proof is to test by functions that allow us to relate the chemotaxis term in (9)₁ with one of the consumption terms in (9)₂. With this idea in mind, we test (9)₁ by $\displaystyle\frac{\mu}{4}F^{\prime}(u)$ to obtain

\frac{d}{dt}\left(\frac{\mu}{4}\int_{\Omega}F(u)d\boldsymbol{x}\right)\,+\,\mu\,{D}_{1}(u)\,-\,\frac{\chi\mu}{2}\int_{\Omega}u\,\sqrt{v}\,{\boldsymbol{\sigma}}\cdot\nabla(F^{\prime}(u))\,d\boldsymbol{x}\,=\,0\,.

(12)

Testing (9)₃ by $\chi{\boldsymbol{\sigma}}$ we obtain

\begin{array}[]{c}\displaystyle\frac{d}{dt}\left(\frac{\chi}{2}\int_{\Omega}|{\boldsymbol{\sigma}}|^{2}d\boldsymbol{x}\right)\,-\,2\chi\int_{\Omega}\frac{1}{\sqrt{v}}{\boldsymbol{\sigma}}^{t}(\nabla{\boldsymbol{\sigma}})^{t}{\boldsymbol{\sigma}}\,d\boldsymbol{x}\,+\,\frac{\chi}{3}\int_{\Omega}\frac{1}{v}|{\boldsymbol{\sigma}}|^{4}d\boldsymbol{x}\,+\,\frac{2\chi}{3}\int_{\Omega}\frac{1}{v}|{\boldsymbol{\sigma}}|^{2}\nabla(\sqrt{v})\cdot{\boldsymbol{\sigma}}d\boldsymbol{x}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle\,+\,\chi\int_{\Omega}(\nabla\cdot{\boldsymbol{\sigma}})^{2}d\boldsymbol{x}\,+\,\chi\int_{\Omega}|\mbox{rot}\,{\boldsymbol{\sigma}}|^{2}d\boldsymbol{x}\,+\,\frac{\chi\mu}{2}\int_{\Omega}u_{+}\,|{\boldsymbol{\sigma}}|^{2}d\boldsymbol{x}\,+\,\frac{\chi\mu}{2}\int_{\Omega}\sqrt{v}\,\nabla u\,\cdot\,{\boldsymbol{\sigma}}\,d\boldsymbol{x}\,=\,0\,.\end{array}

(13)

Using integration by parts we rewrite the fourth term of (13) as

\frac{2\chi}{3}\int_{\Omega}\frac{1}{v}|{\boldsymbol{\sigma}}|^{2}\nabla(\sqrt{v})\cdot{\boldsymbol{\sigma}}d\boldsymbol{x}\,=\,-\frac{2\chi}{3}\int_{\Omega}|{\boldsymbol{\sigma}}|^{2}\nabla\left(\frac{1}{\sqrt{v}}\right)\cdot{\boldsymbol{\sigma}}d\boldsymbol{x}\,=\,\frac{2\chi}{3}\int_{\Omega}\frac{1}{\sqrt{v}}\Big{(}\nabla\cdot(|{\boldsymbol{\sigma}}|^{2}{\boldsymbol{\sigma}})\Big{)}d\boldsymbol{x}\,,

Finally, using $\nabla u=u\nabla(\ln(u))=u\nabla(F^{\prime}(u))$ and adding equations (12) and (13), the terms $\int_{\Omega}u\,\sqrt{v}\,{\boldsymbol{\sigma}}\cdot\nabla(F^{\prime}(u))\,d\boldsymbol{x}$ and $\int_{\Omega}\sqrt{v}\,\nabla u\,\cdot\,{\boldsymbol{\sigma}}\,d\boldsymbol{x}$ cancel, and the desired relation (10) holds. ∎

Corollary 2.3.

In the particular case of considering one-dimensional domains ( $1D$ ), the right side of relation (10) vanishes, implying the following energy dissipative law of the system:

\displaystyle\frac{d}{dt}{E}(u,\sigma)\,+\mu\,{D}_{1}(u)\,+\chi\,{D}_{2}(v,\sigma)\,+\chi\mu\,{D}_{3}(u,\sigma)\,=\,0\,\,.

(14)

Proof.

Since, in one-dimensional domains variable $\sigma$ is a scalar quantity, then the term $R(v,\sigma)$ reads:

R(v,\sigma)\,=\,\displaystyle-\frac{2\chi}{3}\int_{\Omega}\frac{1}{\sqrt{v}}\Big{(}\partial_{x}(\sigma^{3})-3\sigma(\partial_{x}\sigma)\sigma\Big{)}dx\,=\,0\,.

(15)

∎

3 Numerical Schemes

We discretize the time interval $[0,T]$ using Finite Differences and the spatial domain $\Omega\subset\mathbb{R}^{d}$ $d=1,2,3$ using Finite Elements with a shape-regular and quasi-uniform family of triangulations of $\Omega$ , denoted by $\{\mathcal{T}_{h}\}_{h>0}$ . For the sake of simplicity we consider a constant time step $\Delta t:=T/N$ , where $N$ represents the total number of time intervals considered and we denote by $\delta_{t}$ the (backward) discrete time derivative

\delta_{t}u^{n+1}\,:=\,\frac{u^{n+1}-u^{n}}{\Delta t}\,.

In our analysis we will consider $C^{0}$ -FE spaces of order $1$ (denoted by $\mathbb{P}_{1}$ ) for the approximation of $(u,v,{\boldsymbol{\sigma}})$ via the discrete spaces $U_{h}$ , $V_{h}$ and $\Sigma_{h}$ . Additionally, we consider in our schemes mass-lumping ideas [5] to help us achieve positivity of the unknowns in some of the proposed schemes. In order to do that, we introduce the discrete semi-inner product on $C^{0}(\overline{\Omega})$ and its induced discrete seminorm:

(\phi,\psi)_{h}\,:=\,\int_{\Omega}I_{h}(\phi\psi)\,\qquad\mbox{ and }\qquad|\phi|_{h}\,:=\,\sqrt{(\phi,\phi)_{h}}\,.

(16)

with $I_{h}(f(x))$ denoting the nodal $\mathbb{P}_{1}$ -interpolation of the function $f(x)$ .

3.1 UV-schemes

In this section we present three different schemes to approximate the weak formulation that appears when we test (3)₁ and (3)₂ by regular enough test functions $\bar{u}$ and $\bar{v}$ , that is, find $(u,v):[0,T]\times\overline{\Omega}\mapsto\mathbb{R}^{2}$ such that for all $(\bar{u},\bar{v}):\overline{\Omega}\mapsto\mathbb{R}^{2}$ :

\left\{\begin{array}[]{rcl}(u_{t},\bar{u})+(\nabla u,\nabla\bar{u})-\chi\,(u\nabla v,\nabla\bar{u})&=&0\,,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr(v_{t},\bar{v})+(\nabla v,\nabla\bar{v})+\mu\,(uv,\bar{v})&=&0\,.\end{array}\right.

(17)

3.1.1 Scheme UV

•

[Step 1] Given $(u^{n},v^{n})\in U_{h}\times V_{h}$ , find $u^{n+1}\in U_{h}=\mathbb{P}_{1}$ solving the linear problem:

\displaystyle\left(\delta_{t}u^{n+1},\bar{u}\right)_{h}+\Big{(}\nabla u^{n+1},\nabla\bar{u}\Big{)}-\chi\Big{(}u^{n+1}\nabla v^{n},\nabla\bar{u}\Big{)}=0\quad\forall\,\bar{u}\in U_{h}\,.

(18)

•

[Step 2] Given $(u^{n+1},v^{n})\in U_{h}\times V_{h}$ , find $v^{n+1}\in V_{h}=\mathbb{P}_{1}$ solving the linear problem:

\displaystyle\left(\delta_{t}v^{n+1},\bar{v}\right)_{h}+\Big{(}\nabla v^{n+1},\nabla\bar{v}\Big{)}+\mu\Big{(}(u^{n+1})_{+}v^{n+1},\bar{v}\Big{)}_{h}=0\quad\forall\,\bar{v}\in V_{h}\,.

(19)

Note that scheme UV is linear, decoupled and conservative, because taking $\bar{u}=1$ in (18),

\int_{\Omega}u^{n+1}d\boldsymbol{x}\,=\,\int_{\Omega}u^{n}d\boldsymbol{x}\,=\,\cdots\,=\,\int_{\Omega}u^{0}d\boldsymbol{x}\,=:\,m_{0}\,.

Theorem 3.1 ( $u^{n+1}$ -problem).

Given $v^{n}\in V_{h}$ , assuming $\Delta t$ small enough satisfying

\Delta t<\frac{2}{\chi^{2}\|\nabla v^{n}\|_{\infty}},

(20)

then there exist a unique solution $u^{n+1}\in U_{h}$ solving (18).

Proof.

Since problem (18) is a squared algebraic linear system, it suffices to prove uniqueness (that implies existence). In fact, it suffices to prove that if $u\in U_{h}$ is a solution of the corresponding homogeneous problem

\frac{1}{\Delta t}(u,\bar{u})_{h}+(\nabla u,\nabla\bar{u})-\chi\Big{(}u\nabla v^{n},\nabla\bar{u}\Big{)}=0,

then $u=0$ . In fact, testing by $\bar{u}=u$ and applying Holder and Young inequalities, one has

\displaystyle\frac{1}{\Delta t}\int_{\Omega}I_{h}\Big{(}u^{2}\Big{)}d\boldsymbol{x}+\frac{1}{2}\int_{\Omega}|\nabla u|^{2}d\boldsymbol{x}\leq\displaystyle\frac{\chi^{2}}{2}\int_{\Omega}u^{2}|\nabla v^{n}|^{2}d\boldsymbol{x}\leq\displaystyle\frac{\chi^{2}}{2}\|\nabla v^{n}\|_{\infty}\int_{\Omega}u^{2}d\boldsymbol{x}\,.

Using the relation $\int_{\Omega}u^{2}d\boldsymbol{x}\leq\int_{\Omega}I_{h}\Big{(}u^{2}\Big{)}d\boldsymbol{x}$ (see [10]), we obtain

\displaystyle\left(\frac{1}{\Delta t}-\frac{\chi^{2}}{2}\|\nabla v^{n}\|_{\infty}\right)\int_{\Omega}u^{2}d\boldsymbol{x}+\frac{1}{2}\int_{\Omega}|\nabla u|^{2}d\boldsymbol{x}\leq 0\,.

Therefore, assuming hypothesis (20), we have $u=0$ . ∎

Theorem 3.2 ( $v^{n+1}$ -problem).

There exist a unique solution $v^{n+1}\in V_{h}$ solving (19). Moreover, if the triangulation $\{\mathcal{T}_{h}\}$ is acute, that is, all angles of the simplices are less or equal than ${\pi}/2$ , then the discrete maximum principle holds, that is,

\mbox{if }v^{n}>0\quad\mbox{then}\quad v^{n+1}>0,\quad\mbox{ and if }v^{n}\leq M\quad\mbox{then}\quad v^{n+1}\leq M\,.

(21)

Finally, the following weak estimate holds:

\Delta t\sum_{n=0}^{N-1}\int_{\Omega}|\nabla v^{n+1}|^{2}dx\,\leq\,\int_{\Omega}I_{h}\big{(}(v^{0})^{2}\big{)}dx\,.

(22)

Proof.

Existence and uniqueness. Since problem (19) is linear, it suffices to prove that if $v\in V_{h}$ is a solution of the homogeneous problem

\frac{1}{\Delta t}(v,\bar{v})_{h}+(\nabla v,\nabla\bar{v})+\mu\Big{(}(u^{n+1})_{+}v,\bar{v}\Big{)}_{h}=0\,,

then $v=0$ . Indeed, testing by $\bar{v}=v\in V_{h}$ we obtain

\frac{1}{\Delta t}\int_{\Omega}I_{h}\Big{(}v^{2}\Big{)}d\boldsymbol{x}+\int_{\Omega}|\nabla v|^{2}d\boldsymbol{x}+\mu\int_{\Omega}I_{h}\Big{(}(u^{n+1})_{+}v^{2}\Big{)}d\boldsymbol{x}=0\,.

Thus previous relation implies that $v=0$ .

Discrete Maximum Principle. Assume that $v^{n}>0$ . We can define the following problem: Taking $z^{n}:=\displaystyle\min_{x\in\Omega}v^{n}>0$ , find $z^{n+1}\in\mathbb{R}$ such that

\frac{z^{n+1}-z^{n}}{\Delta t}+\mu\|u^{n+1}_{+}\|_{\infty}z^{n+1}=0\,,

that is,

z^{n+1}=\frac{z^{n}}{1+\mu\Delta t\|u^{n+1}_{+}\|_{\infty}}>0\,.

Looking at $z^{n+1}$ as a constant function, $\nabla z^{n+1}=0$ , then

\left(\frac{z^{n+1}-z^{n}}{\Delta t},\bar{v}\right)_{h}+(\nabla z^{n+1},\nabla\bar{v})+\mu\big{(}\|u^{n+1}_{+}\|_{\infty}z^{n+1},\bar{v}\big{)}_{h}=0\,,\quad\forall\,\bar{v}\in V_{h}.

(23)

Let see that $z^{n+1}$ solving (23) satisfy $z^{n+1}\leq v^{n+1}$ . For this, we define $w^{n+1}=v^{n+1}-z^{n+1}\in V_{h}$ ( $w^{n}=v^{n}-z^{n}\geq 0$ ). Subtracting (19) with (23) we obtain:

\left(\frac{w^{n+1}-w^{n}}{\Delta t},\bar{v}\right)_{h}+(\nabla w^{n+1},\nabla\bar{v})+\mu\big{(}u^{n+1}_{+}w^{n+1}+(u^{n+1}_{+}-\|u^{n+1}_{+}\|_{\infty})z^{n+1},\bar{v}\big{)}_{h}=0\,,\quad\forall\,\bar{v}\in V_{h}.

Testing by $\bar{v}=I_{h}(w^{n+1}_{-})\in V_{h}$ (with $w_{-}:=\min\{w,0\}$ ) and using the relation $I_{h}(v\,v_{-})\,=\,I_{h}((v_{-})^{2})$ :

\begin{array}[]{c}\displaystyle\frac{1}{\Delta t}\int_{\Omega}I_{h}\left((w^{n+1}_{-})^{2}\right)d\boldsymbol{x}-\frac{1}{\Delta t}\int_{\Omega}I_{h}(w^{n}w^{n+1}_{-})d\boldsymbol{x}+(\nabla I_{h}(w^{n+1}_{+})+\nabla I_{h}(w^{n+1}_{-}),\nabla I_{h}(w_{-}^{n+1}))\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle+\mu\int_{\Omega}I_{h}\Big{(}z^{n+1}(u_{+}^{n+1}-\|u^{n+1}_{+}\|_{\infty})w^{n+1}_{-}\Big{)}d\boldsymbol{x}+\mu\int_{\Omega}I_{h}\Big{(}u^{n+1}_{+}(w_{-}^{n+1})^{2}\Big{)}d\boldsymbol{x}\leq 0\,.\end{array}

Using that we are considering an acute triangulation, that is, the interior angles of the triangles or tetrahedra are less or equal than ${\pi}/2$ , we can deduce [5]:

(\nabla I_{h}(w^{n+1}_{+}),\nabla I_{h}(w_{-}^{n+1}))\geq 0\,.

Hence, due to the positivity of all the integrands, we can deduce that $I_{h}(w^{n+1}_{-})=0$ so $w^{n+1}\geq 0$ and therefore $v^{n+1}=w^{n+1}+z^{n+1}\geq z^{n+1}>0$ .

On the other hand, we can rewrite (19) as

\frac{1}{\Delta t}\displaystyle(v^{n+1}-M,\bar{v})_{h}+\Big{(}\nabla(v^{n+1}-M),\nabla\bar{v}\Big{)}+\mu\Big{(}(u^{n+1})_{+}(v^{n+1}-M+M),\bar{v}\Big{)}_{h}=\frac{1}{\Delta t}\displaystyle(v^{n}-M,\bar{v})_{h}

(24)

for all $\bar{v}\in V_{h}$ . Then, testing (24) by $\bar{v}=I_{h}((v^{n+1}-M)_{+})$ (and using that $v^{n}\leq M$ ) we obtain

\begin{array}[]{c}\displaystyle\frac{1}{\Delta t}\int_{\Omega}I_{h}\Big{[}\big{(}(v^{n+1}-M)_{+}\big{)}^{2}\Big{]}d\boldsymbol{x}+\displaystyle\int_{\Omega}\Big{[}\nabla\big{(}I_{h}(v^{n+1}-M)_{+}\big{)}\Big{]}^{2}d\boldsymbol{x}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle+\displaystyle\int_{\Omega}\nabla\big{(}I_{h}(v^{n+1}-M)_{-}\big{)}\cdot\nabla\big{(}I_{h}(v^{n+1}-M)_{+}\big{)}d\boldsymbol{x}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle+\mu\int_{\Omega}I_{h}\Big{(}(u^{n+1})_{+}\big{(}[(v^{n+1}-M)_{+}]^{2}+M(v^{n+1}-M)_{+}\big{)}\Big{)}d\boldsymbol{x}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle\leq\frac{1}{\Delta t}\int_{\Omega}I_{h}((v^{n}-M)(v^{n+1}-M)_{+})d\boldsymbol{x}\leq 0\,.\end{array}

Using that the interior angles of the triangles or tetrahedra are less or equal than $\pi/2$ we can deduce

\displaystyle\int_{\Omega}\nabla\big{(}I_{h}(v^{n+1}-M)_{-}\big{)}\cdot\nabla\big{(}I_{h}(v^{n+1}-M)_{+}\big{)}d\boldsymbol{x}\geq 0\,.

that implies $I_{h}\big{[}\big{(}(v^{n+1}-M)_{+}\big{)}^{2}\big{]}=0$ and therefore $v^{n+1}\leq M$ .

Weak Estimate. Testing (19) by $v^{n+1}$ we obtain

\begin{array}[]{rcl}\displaystyle\frac{1}{2\Delta t}\int_{\Omega}I_{h}\big{(}(v^{n+1})^{2}\big{)}d\boldsymbol{x}-\frac{1}{2\Delta t}\int_{\Omega}I_{h}\big{(}(v^{n})^{2}\big{)}d\boldsymbol{x}+\frac{1}{2\Delta t}\int_{\Omega}I_{h}\big{(}(v^{n+1}-v^{n})^{2}\big{)}d\boldsymbol{x}&&\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle+\int_{\Omega}|\nabla v^{n+1}|^{2}d\boldsymbol{x}+\int_{\Omega}I_{h}\big{(}(u^{n+1})_{+}(v^{n+1})^{2}\big{)}d\boldsymbol{x}&=&0\,.\end{array}

Adding the previous relation for $n$ from $0$ to $N-1$ , and multiplying by $2\Delta t$

\int_{\Omega}I_{h}\big{(}|v^{N}|^{2}\big{)}d\boldsymbol{x}\displaystyle+2\,\Delta t\sum_{n=0}^{N-1}\int_{\Omega}|\nabla v^{n+1}|^{2}d\boldsymbol{x}\leq\int_{\Omega}I_{h}\big{(}|v^{0}|^{2}\big{)}d\boldsymbol{x}\,,

and estimate (22) is derived.

∎

3.1.2 Scheme UV-ND (Nonlinear Diffusion)

The main idea of this scheme is to rewrite the diffusion term in a way that we can use it for obtaining a discrete version of inequality (5). In order to do that, for any $\varepsilon>0$ small enough, we introduce a new functional $G_{\varepsilon}(u)$ such that is a $C^{2}$ -approximation of $G(u)$ . This can be achieved by defining

G_{\varepsilon}^{\prime\prime}(u)\,:=\,\left\{\begin{array}[]{ccrcl}\frac{1}{\varepsilon^{2}}&\mbox{ if }&&u&<\varepsilon\,,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\frac{1}{u^{2}}&\mbox{ if }&\displaystyle\varepsilon\leq&u&\leq\frac{1}{\varepsilon}\,,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\varepsilon^{2}&\mbox{ if }&\frac{1}{\varepsilon}\leq&u&\,,\end{array}\right.

(25)

with $G_{\varepsilon}^{\prime}(u)$ , $G_{\varepsilon}(u)$ being the corresponding integral functions such that $G_{\varepsilon}^{\prime}(u)=-1/u$ and $G_{\varepsilon}(u)=-\log(u)+\frac{1}{\varepsilon}$ when $\varepsilon\leq u\leq 1/\varepsilon$ , assuring that $G_{\varepsilon}(u)\geq 0$ for all $u\in\mathbb{R}$ .

Then the scheme UV-ND reads:

•

[Step 1] Find $u^{n+1}\in U_{h}=\mathbb{P}_{1}$ solving the nonlinear problem:

\displaystyle\left(\delta_{t}u^{n+1},\bar{u}\right)_{h}+\Big{(}(u^{n+1})^{2}\nabla(I_{h}G_{\varepsilon}^{\prime}(u^{n+1})),\nabla\bar{u}\Big{)}-\chi\Big{(}u^{n+1}\nabla v^{n},\nabla\bar{u}\Big{)}=0\quad\forall\,\bar{u}\in U_{h}\,.

(26)

Note that (26) is nonlinear and conservative (taking $\bar{u}=1$ one has $\int_{\Omega}u^{n+1}d\boldsymbol{x}=\int_{\Omega}u^{n}d\boldsymbol{x}$ ).

•

[Step 2] Find $v^{n+1}\in V_{h}=\mathbb{P}_{1}$ solving the same linear problem (19) presented in scheme UV. Consequently, Theorem 3.2 also holds for this scheme.

Notice that parameter $\varepsilon>0$ has been introduced for the spatial approximation, and then it will be taken as $\varepsilon=\varepsilon(h)\to 0$ as $h\to 0$ .

Theorem 3.3.

If $(u^{n+1},v^{n+1})$ solves scheme UV-ND, then the following discrete inequality holds:

\displaystyle\delta_{t}\left(\int_{\Omega}I_{h}G_{\varepsilon}(u^{n+1})d\boldsymbol{x}\right)\,+\,\frac{1}{2}\int_{\Omega}(u^{n+1})^{2}\big{|}\nabla\big{(}I_{h}G_{\varepsilon}^{\prime}(u^{n+1})\big{)}\big{|}^{2}d\boldsymbol{x}\,\leq\,\frac{\chi^{2}}{2}\int_{\Omega}|\nabla v^{n+1}|^{2}d\boldsymbol{x}\,.

(27)

Proof.

By using Taylor expansion and the convexity of $G_{\varepsilon}$ , due to $G_{\varepsilon}^{\prime\prime}(u)\geq 0$ for all $u$ , we have

\displaystyle\int_{\Omega}I_{h}\left(({u^{n+1}-u^{n}})G_{\varepsilon}^{\prime}(u^{n+1})\right)d\boldsymbol{x}\geq\displaystyle(G_{\varepsilon}(u^{n+1}),1)_{h}-(G_{\varepsilon}(u^{n}),1)_{h}\,,

(28)

which is nonnegative because $G_{\varepsilon}^{\prime}(u)$ is an increasing function ( $G_{\varepsilon}^{\prime}(u)\sim-1/u$ ). Therefore, testing (26) by $I_{h}G^{\prime}_{\varepsilon}(u^{n+1})$ we have

\delta_{t}\left(\int_{\Omega}I_{h}G_{\varepsilon}(u^{n+1})d\boldsymbol{x}\right)\,+\,\int_{\Omega}(u^{n+1})^{2}\big{|}\nabla\big{(}I_{h}G_{\varepsilon}^{\prime}(u^{n+1})\big{)}\big{|}^{2}d\boldsymbol{x}\,\leq\,\chi\int_{\Omega}u^{n+1}\nabla\big{(}I_{h}G_{\varepsilon}^{\prime}(u^{n+1})\big{)}\cdot\nabla v^{n+1}\,d\boldsymbol{x}\,.

Finally, by using Hölder inequality we arrive at expression (27). ∎

Corollary 3.4.

If $(u^{n+1},v^{n+1})$ solves scheme UV-ND, the following estimates hold

\int_{\Omega}I_{h}G_{\varepsilon}(u^{n+1})d\boldsymbol{x}\leq C\,,\quad\forall\,n\,,

(29)

\Delta t\sum_{n=0}^{N-1}\int_{\Omega}(u^{n+1})^{2}\big{|}\nabla\big{(}I_{h}G_{\varepsilon}^{\prime}(u^{n+1})\big{)}\big{|}^{2}d\boldsymbol{x}\,\leq\,C\,,

(30)

where $C$ is a constant that bounds both $\int_{\Omega}I_{h}G_{\varepsilon}(u^{0})dx$ and $\int_{\Omega}I_{h}\big{(}(v^{0})^{2}\big{)}dx$ . Moreover, the following estimates hold

\int_{\Omega}\big{(}I_{h}(u_{-}^{n})\big{)}^{2}dx\leq C\,\varepsilon^{2}\quad\mbox{ and }\quad\|u^{n}\|_{L^{1}}\leq m_{0}+C\varepsilon\quad\forall\,n\geq 1\,.

(31)

Proof.

Adding up (27), using estimate (22) and $G_{\varepsilon}(u)\geq 0$ for all $u\in\mathbb{R}^{d}$ , we can derive the estimates (29) and (30). On the other hand, estimate (31) is deduced from the inequality

\frac{1}{\varepsilon^{2}}(u_{-})^{2}\,\leq\,G_{\varepsilon}(u)\,,\quad\forall\,u\in\mathbb{R}\,,

following the same arguments presented in [12]. ∎

Remark 3.5.

From estimate (31), we can say that scheme UV-ND satisfies an approximate positivity property for $u^{n}$ , because $I_{h}(u_{-}^{n})\rightarrow 0$ in $L^{2}(\Omega)$ as $\varepsilon\rightarrow 0$ , with $\mathcal{O}(\varepsilon)$ accuracy rate.

3.1.3 Scheme UV-NS (Nonlinear Sensitivity)

In this section we present another scheme developed to satisfy a discrete version of inequality (5). The key point now is to rewrite the sensitivity term in a nonlinear way introducing $d$ new functionals $\Lambda^{i}_{\varepsilon}(u):U_{h}\rightarrow\mathbb{P}_{0}$ ( $i=1,\dots,d$ ) such that they satisfy

\Big{(}\Lambda^{i}_{\varepsilon}(u)\partial_{x_{i}}\big{(}I_{h}G_{\varepsilon}^{\prime}(u)\big{)}\Big{)}^{2}\,=\,\partial_{x_{i}}u\,\partial_{x_{i}}\big{(}I_{h}G_{\varepsilon}^{\prime}(u)\big{)}\qquad\forall\,i=1,\dots,d\,,

(32)

that is, $\Lambda^{i}_{\varepsilon}(u)$ are constant by elements functions such that (32) holds in each element of the triangulation. In fact, (32) holds in any dimension by imposing the constraint of considering a structured mesh and the choice of $U_{h}=\mathbb{P}_{1}$ , as it has been done with related expressions in [2, 12, 13].

Remark 3.6 (How to construct $\Lambda^{i}_{\varepsilon}(u)$ ).

In the one-dimensional case, the domain $\Omega=[a,b]$ can be splitted into $N$ subintervals named $I_{j}$ with $I_{j}=[x_{j},x_{j+1}]$ ( $1\leq j\leq N$ ), being $x_{j}$ the nodes of the partition. Moreover, the discrete derivative with respect to $x$ can be defined as the vector of length $N$ with components:

\delta_{x}u\big{|}_{I_{j}}=\frac{u_{j+1}-u_{j}}{|I_{j}|}\,,

with $u_{j}\sim u(x_{j})$ and $|I_{j}|$ denoting the length of the interval $I_{j}$ . Then, in the one-dimensional case, $\Lambda^{i}_{\varepsilon}(u)$ can be constructed in the following way:

\Lambda^{1}_{\varepsilon}(u)\big{|}_{I_{j}}\,=\,\left\{\begin{array}[]{cl}\displaystyle\frac{\pm\sqrt{\delta_{x}(u)\big{|}_{I_{j}}\delta_{x}\big{(}I_{h}G_{\varepsilon}^{\prime}(u)\big{)}\big{|}_{I_{j}}}}{\delta_{x}\big{(}I_{h}G_{\varepsilon}^{\prime}(u)\big{)}\big{|}_{I_{j}}}&\quad\mbox{ if }u_{j}\neq u_{j+1}\,,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr u_{j}&\quad\mbox{ if }u_{j}=u_{j+1}\,,\end{array}\right.

(33)

and choosing the sign $\pm$ such that $sign(\Lambda^{1}_{\varepsilon}(u)\big{|}_{I_{j}})=sign((u_{j+1}+u_{j})/2)$ .
The definition (33) can be extended to higher dimensional domains just by using the same construction for each functional $\Lambda^{i}_{\varepsilon}(u)$ , where $I_{j}$ represents now the intervals in the corresponding $i$ -direction.

Now, we state the new scheme UV-NS as:

•

[Step 1] Find $u^{n+1}\in U_{h}=\mathbb{P}_{1}$ solving the nonlinear problem:

\displaystyle\left(\delta_{t}u^{n+1},\bar{u}\right)_{h}+\Big{(}\nabla u^{n+1},\nabla\bar{u}\Big{)}-\chi\Big{(}\Lambda_{\varepsilon}(u^{n+1})\nabla v^{n},\nabla\bar{u}\Big{)}=0\quad\forall\,\bar{u}\in U_{h}\,,

(34)

with $\Lambda_{\varepsilon}(u^{n+1})=\mbox{diag}(\Lambda^{i}_{\varepsilon}(u^{n+1}))_{i=1,\dots,d}$ . Note that (34) is nonlinear and conservative.

•

[Step 2] Find $v^{n+1}\in V_{h}=\mathbb{P}_{1}$ solving the same linear problem (19) presented in scheme UV.

Theorem 3.7.

If $(u^{n+1},v^{n+1})$ solves scheme UV-NS, then the following discrete inequality holds:

\displaystyle\delta_{t}\left(\int_{\Omega}I_{h}G_{\varepsilon}(u^{n+1})d\boldsymbol{x}\right)\,+\,\frac{1}{2}\int_{\Omega}\left|\Lambda_{\varepsilon}(u^{n+1})\nabla\big{(}I_{h}G_{\varepsilon}^{\prime}(u^{n+1})\big{)}\right|^{2}d\boldsymbol{x}\,\leq\,\frac{\chi^{2}}{2}\int_{\Omega}|\nabla v^{n+1}|^{2}d\boldsymbol{x}\,.

(35)

Proof.

Testing (34) by $I_{h}G^{\prime}_{\varepsilon}(u^{n+1})$ we have

\begin{array}[]{c}\displaystyle\delta_{t}\left(\int_{\Omega}I_{h}G_{\varepsilon}(u^{n+1})d\boldsymbol{x}\right)\,+\,\int_{\Omega}\nabla u^{n+1}\cdot\nabla\big{(}I_{h}G^{\prime}_{\varepsilon}(u^{n+1})\big{)}d\boldsymbol{x}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle\,\leq\,\chi\int_{\Omega}\Lambda_{\varepsilon}(u^{n+1})\nabla\big{(}I_{h}G_{\varepsilon}^{\prime}(u^{n+1})\big{)}\cdot\nabla v^{n+1}d\boldsymbol{x}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle\leq\frac{1}{2}\int_{\Omega}\left|\Lambda_{\varepsilon}(u^{n+1})\nabla\big{(}I_{h}G_{\varepsilon}^{\prime}(u^{n+1})\big{)}\right|^{2}d\boldsymbol{x}\,+\,\frac{\chi^{2}}{2}\int_{\Omega}|\nabla v^{n+1}|^{2}d\boldsymbol{x}\,,\end{array}

and taking into account (32) we can deduce expression (35). ∎

Remark 3.8.

The estimates presented in Corollary 3.4 also holds in this scheme, substituting estimate (30) by the corresponding one

\Delta t\sum_{n=0}^{N-1}\int_{\Omega}\left|\Lambda_{\varepsilon}(u^{n+1})\,\nabla\big{(}I_{h}G_{\varepsilon}^{\prime}(u^{n+1})\big{)}\right|^{2}d\boldsymbol{x}\,\leq\,C\,.

(36)

In particular, estimate (31) holds, giving approximate positivity for scheme UV-NS in the sense of Remark 3.8.

3.2 Scheme UVS

In this section we present a scheme that approximates the weak formulation obtained when we test (9) by regular enough functions $(\bar{u},\bar{v},\bar{{\boldsymbol{\sigma}}}):\overline{\Omega}\mapsto\mathbb{R}^{1\times 1\times d}$ by using a Galerkin approximation and mass lumping ideas. In order to do that, we first introduce a function $a_{\varepsilon}$ to define a truncation function of $F(u)$ (and its derivatives), namely $F_{\varepsilon}$ , by considering $\displaystyle F_{\varepsilon}^{\prime\prime}(u):=1/{a_{\varepsilon}(u)}$ as an approximation of $1/u$ , where $a_{\varepsilon}(u)$ is a $C^{1}$ -truncation function of $a(u)=u$ defined as

a_{\varepsilon}(u)\,:=\,\left\{\begin{array}[]{cl}\varepsilon&\quad\mbox{ if }\quad u\leq\varepsilon\,,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr u&\quad\mbox{ if }\quad u\in\left(\varepsilon,\frac{1}{\varepsilon}\right)\,,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\frac{1}{\varepsilon}&\quad\mbox{ if }\quad u\geq\frac{1}{\varepsilon}\,.\end{array}\right.

(37)

In fact, the corresponding integration constants that arise in computations from $F_{\varepsilon}^{\prime\prime}(u)$ to $F_{\varepsilon}(u)$ are fixed considering that $F_{\varepsilon}^{\prime}(u)=\ln(u)$ and $F_{\varepsilon}(u)=u\ln(u)-u+1$ when $u\in(\varepsilon,+\infty)$ .

The proposed numerical scheme UVS is:

•

[Step 1] Find $(u^{n+1},{\boldsymbol{\sigma}}^{n+1})\in U_{h}\times\boldsymbol{\Sigma}_{h}=\mathbb{P}_{1}\times\mathbb{P}_{1}$ with ${\boldsymbol{\sigma}}^{n+1}|_{\partial\Omega}=0$ and solving the coupled and nonlinear system

\left\{\begin{array}[]{rcl}\displaystyle\left(\delta_{t}u^{n+1},\bar{u}\right)_{h}+(\nabla u^{n+1},\nabla\bar{u})-2\chi\big{(}u^{n+1}\,\sqrt{v^{n}}{\boldsymbol{\sigma}}^{n+1},\nabla\bar{u}\big{)}&=&0\,,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle\left(\delta_{t}{\boldsymbol{\sigma}}^{n+1},\bar{{\boldsymbol{\sigma}}}\right)-2\left(\frac{1}{\sqrt{v^{n}}}(\nabla{\boldsymbol{\sigma}}^{n+1})^{t}{\boldsymbol{\sigma}}^{n+1},\bar{{\boldsymbol{\sigma}}}\right)+\frac{1}{3}\left(\frac{1}{v^{n}}|{\boldsymbol{\sigma}}^{n+1}|^{2}{\boldsymbol{\sigma}}^{n+1},\bar{{\boldsymbol{\sigma}}}\right)&&\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle+\frac{2}{3}\left(\frac{1}{v^{n}}\nabla(\sqrt{v^{n}})|{\boldsymbol{\sigma}}^{n+1}|^{2},\bar{{\boldsymbol{\sigma}}}\right)+(\nabla\cdot{\boldsymbol{\sigma}}^{n+1},\nabla\cdot\bar{{\boldsymbol{\sigma}}})+(\mbox{rot}\,{\boldsymbol{\sigma}},^{n+1}\mbox{rot}\,\bar{{\boldsymbol{\sigma}}})&&\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle+\frac{\mu}{2}({\boldsymbol{\sigma}}^{n+1}\,u^{n+1}_{+},\bar{{\boldsymbol{\sigma}}})+\frac{\mu}{2}\Big{(}\sqrt{v^{n}}\,u^{n+1}\,\nabla(I_{h}F_{\varepsilon}^{\prime}(u^{n+1})),\bar{{\boldsymbol{\sigma}}}\Big{)}&=&0\,,\end{array}\right.

(38)

for all $(\bar{u},\bar{{\boldsymbol{\sigma}}})\in U_{h}\times\boldsymbol{\Sigma}_{h}$ with $\bar{{\boldsymbol{\sigma}}}|_{\partial\Omega}=0$ . In the last term of the ${\boldsymbol{\sigma}}$ -system, $\nabla u$ has been approximated by $u\nabla I_{h}F_{\varepsilon}(u)$ . Note that (38) is a nonlinear and conservative problem.

•

[Step 2] Find $v^{n+1}\in V_{h}=\mathbb{P}_{1}\subset H^{1}(\Omega)$ solving the same linear problem (19) presented in scheme UV.

Now we present a result that provides a discrete version of the estimate (10) derived in Theorem 2.2.

Theorem 3.9.

If $(u^{n+1},{\boldsymbol{\sigma}}^{n+1})$ solves Step 1 (38), then the following discrete version of (10) holds:

\displaystyle\delta_{t}{E}_{h}^{\varepsilon}(u^{n+1},{\boldsymbol{\sigma}}^{n+1})\,+\mu\,{D}_{1,h}^{\varepsilon}(u^{n+1})\,+\chi\,{D}_{2}(v^{n},{\boldsymbol{\sigma}}^{n+1})\,+\chi\mu\,{D}_{3}(u^{n+1},{\boldsymbol{\sigma}}^{n+1})\,\leq\,R(v^{n},{\boldsymbol{\sigma}}^{n+1})\,,

(39)

where ${D}_{2}(v,{\boldsymbol{\sigma}})$ , ${D}_{3}(u,{\boldsymbol{\sigma}})$ and $R(v,{\boldsymbol{\sigma}})$ are defined in (11) and

{E}_{h}^{\varepsilon}(u,{\boldsymbol{\sigma}})\,:=\,\displaystyle\frac{\mu}{4}\int_{\Omega}I_{h}F_{\varepsilon}(u)d\boldsymbol{x}\,+\,\frac{\chi}{2}\int_{\Omega}|{\boldsymbol{\sigma}}|^{2}d\boldsymbol{x}\quad\mbox{ and }\quad{D}_{1,h}^{\varepsilon}(u)\,:=\,\frac{1}{4}\int_{\Omega}\nabla u\cdot\nabla\big{(}I_{h}F_{\epsilon}^{\prime}(u)\big{)}\,d\boldsymbol{x}\,.

Proof.

Testing (38)₁ by $\bar{u}=\displaystyle\frac{\mu}{4}I_{h}F^{\prime}_{\varepsilon}(u^{n+1})$ and using Taylor expansion and convexity of $F_{\varepsilon}$ (as we have done for $G_{\varepsilon}$ in (28)) we have

\delta_{t}\left(\frac{\mu}{4}\int_{\Omega}I_{h}F_{\varepsilon}(u^{n+1})d\boldsymbol{x}\right)+\mu{D}_{1,h}^{\varepsilon}(u^{n+1})-\frac{\chi\mu}{2}\int_{\Omega}u^{n+1}\,\sqrt{v^{n}}{\boldsymbol{\sigma}}^{n+1}\cdot\nabla(I_{h}F^{\prime}_{\varepsilon}(u^{n+1}))\,d\boldsymbol{x}\,\leq\,0\,.

(40)

Testing (38)₂ by $\bar{{\boldsymbol{\sigma}}}=\chi{\boldsymbol{\sigma}}^{n+1}$ we have

\begin{array}[]{l}\displaystyle\delta_{t}\left(\frac{\chi}{2}\int_{\Omega}({\boldsymbol{\sigma}}^{n+1})^{2}dx\right)-\displaystyle 2\chi\int_{\Omega}\frac{1}{\sqrt{v^{n}}}({\boldsymbol{\sigma}}^{n+1})^{t}(\nabla({\boldsymbol{\sigma}}^{n+1}))^{t}{\boldsymbol{\sigma}}^{n+1}\,d\boldsymbol{x}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle+\frac{2\chi}{3}\int_{\Omega}\frac{1}{v^{n}}|{\boldsymbol{\sigma}}^{n+1}|^{2}\nabla(\sqrt{v^{n}})\cdot{\boldsymbol{\sigma}}^{n+1}d\boldsymbol{x}+\displaystyle\frac{\chi}{3}\int_{\Omega}\frac{1}{v^{n}}|{\boldsymbol{\sigma}}^{n+1}|^{4}d\boldsymbol{x}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle+\chi\int_{\Omega}\Big{(}(\nabla\cdot{\boldsymbol{\sigma}}^{n+1})^{2}+|\mbox{rot}\,{\boldsymbol{\sigma}}^{n+1}|^{2}\Big{)}d\boldsymbol{x}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle+\displaystyle\frac{\chi\mu}{2}\int_{\Omega}|{\boldsymbol{\sigma}}^{n+1}|^{2}\,u^{n+1}_{+}d\boldsymbol{x}+\frac{\chi\mu}{2}\int_{\Omega}\sqrt{v^{n}}\,u^{n+1}\,\nabla(I_{h}F_{\varepsilon}^{\prime}(u^{n+1}))\cdot{\boldsymbol{\sigma}}^{n+1}d\boldsymbol{x}=0\,.\end{array}

(41)

Using integration by parts:

\int_{\Omega}\frac{1}{v^{n}}|{\boldsymbol{\sigma}}^{n+1}|^{2}\nabla(\sqrt{v^{n}})\cdot{\boldsymbol{\sigma}}^{n+1}d\boldsymbol{x}=-\int_{\Omega}|{\boldsymbol{\sigma}}^{n+1}|^{2}\nabla\left(\frac{1}{\sqrt{v^{n}}}\right)\cdot{\boldsymbol{\sigma}}^{n+1}d\boldsymbol{x}=\int_{\Omega}\frac{1}{\sqrt{v^{n}}}\nabla\cdot(|{\boldsymbol{\sigma}}^{n+1}|^{2}{\boldsymbol{\sigma}}^{n+1})d\boldsymbol{x},

then the desired relation (39) holds by adding equations (40) and (41) .

∎

Corollary 3.10.

In the particular case of considering one-dimensional domains ( $1D$ ), relation (39) implies:

\displaystyle\delta_{t}{E}_{h}^{\varepsilon}(u^{n+1},\sigma^{n+1})\,+\mu\,{D}_{1,h}^{\varepsilon}(u^{n+1})\,+\chi\,{D}_{2}(v^{n},\sigma^{n+1})\,+\mu\chi\,{D}_{3}(u^{n+1},\sigma^{n+1})\,\leq\,0\,,

(42)

with ${D}_{1,h}^{\varepsilon}(u^{n+1}),\,{D}_{2}(v^{n},\sigma^{n+1}),\,{D}_{3}(u^{n+1},\sigma^{n+1})\geq 0$ . In particular, scheme UVS is unconditional energy-stable with respect to the energy ${E}_{h}^{\varepsilon}(u,\sigma)$ defined in (11)₁.

Proof.

From Theorem 3.9, we only need to prove that ${D}_{1,h}^{\varepsilon}(u^{n+1})\geq 0$ and $R(v^{n},\sigma^{n+1})=0$ . Since $u^{n+1}\in U_{h}=\mathbb{P}_{1}$ and we are working in $1D$ , we can write:

{D}_{1,h}^{\varepsilon}(u^{n+1})\,=\,\frac{1}{4}\int_{\Omega}(u_{x}^{n+1})(I_{h}\,F_{\varepsilon}^{\prime}(u^{n+1}))_{x}\,dx\,=\,\frac{1}{4}\sum_{j=1}^{J}\left(\frac{u^{n+1}_{j+1}-u^{n+1}_{j}}{h}\right)\left(\frac{F_{\varepsilon}^{\prime}(u_{j+1}^{n+1})-F_{\varepsilon}^{\prime}(u_{j}^{n+1})}{h}\right)\,,

which is nonnegative because $F_{\varepsilon}^{\prime}(u)$ is an increasing function ( $F_{\varepsilon}^{\prime}(u)\sim\ln u$ ). Finally, in one-dimensional domains variable $\sigma$ is a scalar quantity, so the term $R(v^{n},\sigma^{n+1})$ reads:

R(v^{n},\sigma^{n+1})\,=\,-\displaystyle\frac{2\chi}{3}\int_{\Omega}\frac{1}{\sqrt{v^{n}}}\Big{(}\partial_{x}((\sigma^{n+1})^{3})-3(\partial_{x}\sigma^{n+1})(\sigma^{n+1})^{2}\Big{)}dx\,=\,0\,.

∎

Corollary 3.11.

In the particular case of considering one-dimensional domains ( $1D$ ), the following estimates hold

\int_{\Omega}I_{h}F_{\varepsilon}(u^{n+1})dx\,\leq\,C\quad\forall\,n\geq 1\,,

(43)

\Delta t\sum_{n=0}^{N-1}\Big{(}\mu\,{D}_{1,h}^{\varepsilon}(u^{n+1})\,+\chi\,{D}_{2}^{\varepsilon}(v^{n},\sigma^{n+1})\,+\mu\chi\,{D}_{3}(u^{n+1},\sigma^{n+1})\Big{)}\,\leq\,C\,,

(44)

where $C>0$ is a constant that bounds $\int_{\Omega}I_{h}F_{\varepsilon}(u^{0})dx$ . Moreover, the following estimates also hold

\int_{\Omega}\big{(}I_{h}(u_{-}^{n})\big{)}^{2}dx\,\leq\,C\,\varepsilon\quad\mbox{ and }\quad\|u^{n}\|_{L^{1}}\leq m_{0}+C\sqrt{\varepsilon}\quad\forall\,n\geq 1\,.

(45)

Proof.

The proof follows the same arguments presented in Corollary 3.4 considering now the estimate

\frac{1}{\varepsilon}(u_{-})^{2}\,\leq\,F_{\varepsilon}(u),\quad\forall\,u\in\mathbb{R}.

∎

Remark 3.12.

From estimate (45), we can say that the one-dimensional version of scheme UVS satisfies an approximate positivity property for $u^{n}$ , because $I_{h}(u_{-}^{n})\rightarrow 0$ in $L^{2}(\Omega)$ as $\varepsilon\rightarrow 0$ , with $\mathcal{O}(\sqrt{\varepsilon})$ order.

4 Numerical simulations in $1D$ domains

The aim of this section is to report the numerical results obtained carrying out simulations using the schemes presented through the paper in one-dimensional domains. The idea is to illustrate the type of dynamics exhibit by chemo-attraction and consumption models and to compare the effectiveness of each of the presented schemes. All the simulations have been performed using MATLAB software [18].

We consider a regular partition of the spatial domain $\overline{\Omega}=[a,b]$ denoted by $T_{h}:=\bigcup_{i=1}^{J-1}I_{i}$ , where $J$ denotes the number of nodes, $\{x_{j}\}_{j=1}^{J}$ the coordinates of these nodes and $h$ the size of the mesh, that for simplicity we assume that is constant in the whole domain. We will compare the schemes presented in Section 3 together with a conservative and positive Finite Volume scheme that can be viewed as a Finite Element scheme with artificial diffusion, and it has been derived following the ideas presented in [25].

The physical and discrete parameters for each example will be detailed in each subsection except the value of the truncation parameter and the tolerance used for the iterative methods for the nonlinear schemes, that are going to be always chosen as:

\varepsilon\,=\,h^{2}\quad\mbox{ and }\quad C_{tol}\,=\,10^{-8}\,.

(46)

The section is organized as follows: Firstly we introduce the iterative algorithms to approximate the nonlinear problems and the Finite Element scheme with artificial diffusion equivalent to a conservative and positive Finite Volume scheme. Then we present one example with the complete dynamics, that is, until it reaches the equilibrium configuration. The purpose of this example is to show that the system tends to a flat/constant equilibrium configuration of $u$ (with $u_{eq}=m_{0}=\frac{1}{|\Omega|}\int_{\Omega}u^{0}dx$ ) and $v=0$ , while the energy $E(u,v)$ decreases and the volume of $u$ remains constant. After that, we compare the ability of each scheme to maintain the positivity of the $u$ unknown using a choice of the initial condition and physical parameters designed in such a way that in some regions $u$ tends to be very close to zero while in other parts the value of $u$ is far away form zero. In the third part we focus on studying how each scheme capture the evolution of the energy in time using different values of the spatial and time discretization parameters. Finally, we perform a numerical study of the experimental order of convergence in space for each scheme.

4.1 Iterative Methods

Some of the presented schemes are nonlinear. In the following we detail the iterative algorithms that we have considered for approximating each of the nonlinear systems.

4.1.1 Scheme UV-ND (Nonlinear Diffusion)

Iterative algorithm to approximate the nonlinear problem (26): Find $u_{\ell+1}\in U_{h}$ such that

\begin{array}[]{rcl}\displaystyle\frac{1}{\Delta t}(u_{\ell+1},\bar{u})_{h}+(\nabla u_{\ell+1},\nabla\bar{u})&=&\displaystyle\frac{1}{\Delta t}(u^{n},\bar{u})_{h}+(\nabla(u_{\ell}),\nabla\bar{u})+\chi\Big{(}u_{\ell}\nabla v^{n},\nabla\bar{u}\Big{)}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle\hfil&&-\Big{(}(u_{\ell})^{2}\nabla(I_{h}G_{\varepsilon}^{\prime}(u^{\ell})),\nabla\bar{u}\Big{)}\end{array}

(47)

Stopping criterium: Iterates until $\frac{\|u_{\ell+1}-u_{\ell}\|_{H^{1}(\Omega)}}{\|u_{\ell}\|_{H^{1}(\Omega)}}\,\leq\,C_{tol}\,,$ with $C_{tol}$ a given tolerance.

4.1.2 Scheme UV-NS (Nonlinear Sensitivity)

Iterative algorithm to approximate the nonlinear problem (34): Find $u_{\ell+1}\in U_{h}$ such that

\displaystyle\frac{1}{\Delta t}(u_{\ell+1},\bar{u})_{h}+(\nabla u_{\ell+1},\nabla\bar{u})=\displaystyle\frac{1}{\Delta t}(u^{n},\bar{u})_{h}+\chi\Big{(}\Lambda_{\varepsilon}(u_{\ell})\nabla v^{n},\nabla\bar{u}\Big{)}

(48)

Stopping criterium: Iterate until $\frac{\|u_{\ell+1}-u_{\ell}\|_{H^{1}(\Omega)}}{\|u_{\ell}\|_{H^{1}(\Omega)}}\,\leq\,C_{tol}\,$ .

4.1.3 Scheme UVS

We consider the following iterative algorithm to approximate this nonlinear problem (38):

Substep.1

Find $u_{\ell+1}\in U_{h}$ such that

\displaystyle\frac{1}{\Delta t}(u_{\ell+1},\bar{u})_{h}+(\nabla(u_{\ell+1}),\nabla\bar{u})\,=\,\frac{1}{\Delta t}(u^{n},\bar{u})_{h}+2\chi\big{(}u_{\ell}\,\sqrt{v^{n}}{\boldsymbol{\sigma}}_{\ell},\nabla\bar{u}\big{)}\,,\quad\forall\,\bar{u}\in U_{h}.

(49)

Substep.2

Find ${\boldsymbol{\sigma}}_{\ell+1}\in\boldsymbol{\Sigma}_{h}$ such that, for all $\bar{{\boldsymbol{\sigma}}}\in\boldsymbol{\Sigma}_{h}$ ,

\begin{array}[]{rcl}\displaystyle\frac{1}{\Delta t}({\boldsymbol{\sigma}}_{\ell+1},\bar{{\boldsymbol{\sigma}}})&+&(\nabla\cdot{\boldsymbol{\sigma}}^{\ell+1},\nabla\cdot\bar{{\boldsymbol{\sigma}}})+(\mbox{rot}\,{\boldsymbol{\sigma}},^{\ell+1}\mbox{rot}\,\bar{{\boldsymbol{\sigma}}})+\frac{\mu}{2}({\boldsymbol{\sigma}}_{\ell+1}\,(u_{\ell+1})_{+},\bar{{\boldsymbol{\sigma}}})\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle\hfil&=&\displaystyle\frac{1}{\Delta t}({\boldsymbol{\sigma}}^{n},\bar{{\boldsymbol{\sigma}}})-\frac{1}{3}\left(\frac{1}{v^{n}}|{\boldsymbol{\sigma}}^{\ell}|^{2}{\boldsymbol{\sigma}}^{\ell},\bar{{\boldsymbol{\sigma}}}\right)\displaystyle+2\left(\frac{1}{\sqrt{v^{n}}}(\nabla{\boldsymbol{\sigma}}^{\ell})^{t}{\boldsymbol{\sigma}}^{\ell},\bar{{\boldsymbol{\sigma}}}\right)\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle\hfil&-&\displaystyle\frac{2}{3}\left(\frac{1}{v^{n}}\nabla(\sqrt{v^{n}})|{\boldsymbol{\sigma}}^{\ell}|^{2},\bar{{\boldsymbol{\sigma}}}\right)-\displaystyle\frac{\mu}{2}\Big{(}\sqrt{(}v^{n})\,u_{\ell}\,\nabla(I_{h}F_{\varepsilon}^{\prime}(u_{\ell})),\bar{{\boldsymbol{\sigma}}}\Big{)}\,.\end{array}

(50)

Substep 3:

Stopping criterium. Iterate until $\frac{\|(u_{\ell+1},{\boldsymbol{\sigma}}_{\ell+1})-(u_{\ell},{\boldsymbol{\sigma}}_{\ell})\|_{H^{1}(\Omega)\times\textbf{H}^{1}(\Omega)}}{\|(u_{\ell},{\boldsymbol{\sigma}}_{\ell})\|_{H^{1}(\Omega)\times\textbf{H}^{1}(\Omega)}}\,\leq\,C_{tol}\,.$

4.2 Scheme UV-AD (Artificial Diffusion) in $1D$

We introduce the upwind Finite Volume scheme presented in [25] that can be viewed as a Finite Element scheme with artificial diffusion. For this, we define the interior control volume as $K_{j}=[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}]$ for $j=2,\dots,J-1$ , while the boundary control volumes are defined as $K_{1}=[x_{1},x_{1+\frac{1}{2}}]$ and $K_{J}=[x_{J-\frac{1}{2}},x_{J}]$ , with $x_{j+\frac{1}{2}}=(x_{j}+x_{j+1})/2$ . Moreover we consider the notation for discrete spatial derivatives $\delta_{x}u_{j+\frac{1}{2}}:=(u_{j+1}-u_{j})/h$ for $j=1,\dots,J-1$ . Using this notation the boundary conditions are discretized as

\delta_{x}u_{1-\frac{1}{2}}\,=\,\delta_{x}u_{J+\frac{1}{2}}\,=\,0\quad\quad\mbox{ and }\quad\quad\delta_{x}v_{1-\frac{1}{2}}\,=\,\delta_{x}v_{J+\frac{1}{2}}\,=\,0\,.

The proposed upwind finite volume scheme is the following:

•

[Step 1] For all $j=1,\dots,J$ , find $u_{j}^{n+1}$ solving the linear problem:

\begin{array}[]{c}\displaystyle|K_{j}|\delta_{t}u^{n+1}_{j}\,-\,\Big{(}\delta_{x}u^{n+1}_{j+\frac{1}{2}}-\delta_{x}u^{n+1}_{j-\frac{1}{2}}\Big{)}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\,+\,\chi\Big{(}(\delta_{x}v^{n}_{j+\frac{1}{2}})_{+}u_{j}^{n+1}\,+\,(\delta_{x}v^{n}_{j+\frac{1}{2}})_{-}u_{j+1}^{n+1}\,-\,(\delta_{x}v^{n}_{j-\frac{1}{2}})_{+}u_{j-1}^{n+1}\,-\,(\delta_{x}v^{n}_{j-\frac{1}{2}})_{-}u_{j}^{n+1}\Big{)}=0.\end{array}

(51)

•

[Step 2] For all $j=1,\dots,J$ , find $v_{j}^{n+1}$ solving the linear problem:

|K_{j}|\delta_{t}v^{n+1}_{j}\,-\,\Big{(}\delta_{x}v^{n+1}_{j+\frac{1}{2}}-\delta_{x}v^{n+1}_{j-\frac{1}{2}}\Big{)}+\mu\,u_{j}^{n+1}v_{j}^{n+1}=0\,.

(52)

In fact, it is easy to check that scheme (51)-(52) is equivalent to the following linear Finite Element scheme with artificial diffusion, denoted by UV-AD:

•

[Step 1] Find $u^{n+1}\in U_{h}=\mathbb{P}_{1}\subset H^{1}(\Omega)$ solving the linear problem:

\displaystyle\left(\delta_{t}u^{n+1},\bar{u}\right)_{h}+\Big{(}u^{n+1}_{x},\bar{u}_{x}\Big{)}+h\frac{\chi}{2}\Big{(}|v_{x}^{n}|u^{n+1}_{x},\bar{u}_{x}\Big{)}-\chi\Big{(}u^{n+1}(v^{n})_{x},\bar{u}_{x}\Big{)}=0\quad\forall\,\bar{u}\in U_{h}\,.

(53)

Note that (53) is conservative, because taking $\bar{u}=1$ one has $\int_{\Omega}u^{n+1}dx=\int_{\Omega}u^{n}dx\,.$

•

[Step 2] Find $v^{n+1}\in V_{h}=\mathbb{P}_{1}\subset H^{1}(\Omega)$ solving the same linear problem (19) presented in scheme UV.

4.3 Example I: Dynamic towards constants

In this case we consider as initial configuration two smooth functions with the same amplitude and the physical parameters shown below:

\left\{\begin{array}[]{rcl}u_{0}&=&1.0001+\cos(5\pi x)\,,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr v_{0}&=&1.0001+\cos(2\pi x)\,,\end{array}\right.\quad[0,T]\,=\,[0,0.3]\,,\quad\chi\,=\,100\quad\mbox{ and }\quad\mu=1000\,.

(54)

In Figure 1 we present the dynamics of the system using scheme UV with discretization parameters $h=10^{-3}$ and $\Delta t=10^{-6}$ , where we can observe that initially the system tends to accumulate the cell population density close to the boundaries (due to the attraction of the cell population towards the chemical substance), then the consumption effects dominates the dynamics producing that the amount of $v$ decreases and finally the system start to minimize the gradients of $u$ and the system reaches the expected equilibrium configuration, that is, a flat configuration of $u$ and $v=0$ . Moreover, we can see in Figure 2 that the energy is decreasing until the system reaches the equilibrium configuration and the amount of each unknown is presented in Figure 3, showing that the amount of chemical substance decreases to zero (due to the consumption effects) and the volume cell population density remains constant in time.

Refer to caption — Figure 1: Example I: Dynamic of scheme UV using data in (54). The evolution in time of $u$ and $v$ is presented from Left to Right and from Top to Bottom.

4.4 Example II: Positivity test

The initial configuration and the physical parameters for this example are:

\left\{\begin{array}[]{rcl}u_{0}&=&1.1-e^{-(\frac{x-0.5}{0.1})^{2}}\,,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr v_{0}&=&2-e^{-(\frac{x-0.5}{0.01})^{2}}\,,\end{array}\right.\quad[0,T]\,=\,[0,10^{-4}]\,,\quad\chi\,=\,100\quad\mbox{ and }\quad\mu=1\,.

(55)

The initial configuration corresponds with two symmetric smooth functions with the physical parameters chosen such that the attraction effects dominates over the consumption ones, in order to produce a dynamic with the variable $u$ tending to zero in the central region, and with high gradients on it, in order to test how well the schemes maintain at a discrete level the positivity of the unknown $u$ . The exact dynamics of this example is presented in Figure 4 and it has been computed using scheme UV with parameters $h=10^{-5}$ and $\Delta t=10^{-9}$ .

In Table 1 we present the minimum values achieved by $u$ in the whole domain $\Omega$ during the complete time interval $[0,T]$ (i.e. $\min_{t\in[0,T]}\min_{x\in\Omega}u$ ), for each of the schemes using different values of the discretization parameters $h$ and $\Delta t$ . Scheme UV is able to maintain $u$ positive using the two time steps considered ( $\Delta t=10^{-7},10^{-8}$ ) but only for $h$ small enough, in fact when $h\leq 1/5000$ . This fact can be explained by taking into account that this discrete formulation when $h$ is small enough is the closest one to the continuous in space problem, which satisfies a maximum principle. Schemes UV-ND and UV-NS do not work well when $\Delta t=10^{-7}$ , because the iterative methods (49) and (50) are not convergent for that choice of the time step. The reason why we obtain some values for scheme UV-NS is because we have added an additional condition to the iterative loop, which states that if after $100$ iterations the stopping criterium is not satisfied, the algorithm has to take the last obtained result and move forward (hoping that this idea might help the scheme in case that the desired tolerance is not achieved but the error is getting close to it). This idea does not prevent scheme UV-ND to crash, but it helps somehow scheme UV-AD to obtain some intermediate results, not satisfying a good approximation of scheme UV-AD. On the other hand, if we consider $\Delta t=10^{-8}$ , then the iterative schemes associated with both schemes, UV-ND and UV-NS, are convergent and both schemes are able to achieve the expected approximated positivity (the obtained values are not strictly positive but they are very small in absolute value). The obtained results are in agreement with the results in Remarks 3.5 and 3.8, which state that schemes UV-ND and UV-NS will satisfy the positivity constraint approximatively. For scheme UV-AD we can observe that it is always capture the positivity of the unknown $u$ as expected, given that upwind schemes were designed with that purpose in mind. Finally, scheme UVS is not reliable for $\Delta t=10^{-7}$ due to the fact that for small values of $h$ the simulations crashes. But it performs very well when $\Delta t=10^{-8}$ achieving approximated positivity when $h=1/1000$ and strictly positivity when $h\leq 1/5000$ . The obtained results are in agreement with Remark 3.12, which states that scheme UVS will satisfy the positivity constraint approximatively.

$\displaystyle\min_{t\in[0,T]}\min_{x\in\Omega}u$	$h=1/100$	$h=1/500$	$h=1/1000$	$h=1/5000$	$h=1/10000$
Scheme UV
$\Delta t=10^{-7}$	$-0.3414$	$-0.0674$	$-3.24\times 10^{-4}$	$6.45\times 10^{-22}$	$7.96\times 10^{-22}$
$\Delta t=10^{-8}$	$-0.3450$	$-0.0838$	$-4.50\times 10^{-4}$	$2.61\times 10^{-22}$	$3.28\times 10^{-22}$
Scheme UV-ND
$\Delta t=10^{-7}$	$\times$	$\times$	$\times$	$\times$	$\times$
$\Delta t=10^{-8}$	$-3.13\times 10^{-4}$	$-6.90\times 10^{-6}$	$-9.20\times 10^{-7}$	$-1.43\times 10^{-8}$	$-2.57\times 10^{-9}$
Scheme UV-NS
$\Delta t=10^{-7}$	$-0.0011$	$-1.44\times 10^{-5}$	$-2.25\times 10^{-6}$	$-0.1514$	$-0.1622$
$\Delta t=10^{-8}$	$-0.0011$	$-1.40\times 10^{-5}$	$-1.86\times 10^{-6}$	$-2.68\times 10^{-8}$	$-4.98\times 10^{-9}$
Scheme UV-AD
$\Delta t=10^{-7}$	$4.84\times 10^{-5}$	$1.55\times 10^{-10}$	$3.15\times 10^{-14}$	$5.98\times 10^{-20}$	$7.33\times 10^{-21}$
$\Delta t=10^{-8}$	$4.84\times 10^{-5}$	$1.54\times 10^{-10}$	$2.82\times 10^{-14}$	$3.16\times 10^{-20}$	$3.44\times 10^{-21}$
Scheme UVS
$\Delta t=10^{-7}$	$-0.1425$	$-0.0065$	$-1.50\times 10^{-4}$	$\times$	$\times$
$\Delta t=10^{-8}$	$-0.1454$	$-0.0118$	$-1.63\times 10^{-4}$	$2.32\times 10^{-22}$	$2.90\times 10^{-22}$

Table 1: Example II: Evolution in time of the minimum values achieved by

u

in the domain

\Omega

in the whole interval

[0,T]

for different values of the discretization parameters

h

and

\Delta t

4.5 Example III: Energy stability test

The aim of this test is to compare the approximation of the evolution of the energy for the different schemes. The initial configuration are two smooth functions with different amplitude and the physical parameters are:

\left\{\begin{array}[]{rcl}u_{0}&=&4(2.0001+\cos(7\pi x))\,,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr v_{0}&=&3(2.0001+\cos(12\pi x))\,,\end{array}\right.\quad[0,T]\,=\,[0,10^{-4}]\,,\quad\chi\,=\,30\quad\mbox{ and }\quad\mu=10000\,.

(56)

The parameters are chosen in order that consumption effects are stronger than the attraction ones. We present the exact dynamics of the system in Figure 5, which has been computed using scheme UV with discretization parameters $h=10^{-5}$ and $\Delta t=10^{-8}$ . The evolution of the energy for schemes UV, UV-ND, UV-NS and UV-AD is presented in Figure 6, and for scheme UVS in Figure 7. We can observe how schemes UV, UV-ND, UV-NS and UV-AD approximates well the evolution of the energy once the discretization parameters are small enough (for this example the requirements are $h\leq 10^{-3}$ and $\Delta t\leq 10^{-7}$ ). In fact, it is interesting to mention that scheme UV is able to produce a reasonable approximation even when $\Delta t=10^{-6}$ while schemes UV-ND and UV-NS crashes. Moreover, it is interesting to check Figure 8 to understand why scheme UV-AD does not perform well with the choice $\Delta t=10^{-6}$ . The reason is that this scheme has been designed as an upwind scheme, which are a type of schemes famous for behaving well when approximating transport effects due to some (incompresible) velocity, while maintaining the positivity of the solution variable. But, in the chemotaxis model considered in this work, the transport of $u$ is due to the term $\nabla\cdot(u\nabla v)$ so that the transport velocity is represented by $\nabla v$ which in general is not incompresible (because $\nabla\cdot(\nabla v)\not=0$ in general). Then, it produces nonphysical oscillations in the solution if the time step is not small enough (although the scheme always preserve the positivity of the solution as expected). On the other hand, scheme UVS does not seem to produce as good approximations of the evolution of the exact energy $E(u,v)$ as the other schemes are producing with the considered time steps. In fact, what it is known by Corollary 3.10 is that this scheme satisfy the energy law (39) for $E(u,{\boldsymbol{\sigma}})$ , but it is not clear that this property implies that the evolution of energy $E(u,{\boldsymbol{\sigma}})$ has to exactly match the evolution of energy $E(u,v)$ . As a matter of fact, we can observe in Figure 7, that those energies do not match in this example for the considered time steps.

4.6 Example IV: Approximation and error estimates test

In this example we perform a numerical error estimate study in space for each of the schemes presented in the manuscript. The initial configuration and the physical parameters considered for this test are:

\left\{\begin{array}[]{rcl}u_{0}&=&3(1.0001+\cos(8\pi x))\,,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr v_{0}&=&5(1.0001+\cos(7\pi x))\,,\end{array}\right.\quad[0,T]\,=\,[0,10^{-4}]\,,\quad\chi\,=\,10\quad\mbox{ and }\quad\mu=1500\,.

(57)

We will compute the EOC (Experimental Order of Convergence) taking $\Delta t=10^{-9}$ and using as reference (or exact) solution the one obtained by solving the system using scheme UV with spatial discretization parameter $h=10^{-5}$ at final time $T=10^{-4}$ . The very non-trivial dynamics of this system are presented in Figure 9. We observe in the dynamics that the cell population density $u$ moves towards the regions of high concentration of chemical substance, which is consumed at a very high rate, producing that the location of high concentration regions of chemical substance $v$ changes fast in time, so the cell population density needs to rearrange itself also in a very fast way in order to move towards the new locations with high concentration of substance $v$ .

We now introduce some additional notation. The individual errors using discrete norms and the convergence rate between two consecutive meshes of size $h$ and $\tilde{h}$ are defined as:

e(\phi):=\|\phi_{exact}-\phi_{h}\|_{L^{2}(\Omega)}\quad\mbox{ and }\quad r(\phi):=\left[\log\left(\frac{e(\phi)}{\tilde{e}(\phi)}\right)\right]\left[\log\left(\frac{h}{\tilde{h}}\right)\right]^{-1}\,.

The convergence history for a sequence of triangulations for all the schemes is presented in Table 2. We can observe that schemes UV, UV-ND, UV-NS and UV-AD show order of convergence, and between these four schemes it is clear that UV, UV-ND and UV-NS performs better than UV-AD, due to the fact that this last scheme introduces extra numerical dissipation to help the scheme to achieve the positivity of $u$ , but at the same time this numerical dissipation prevents the scheme to achieve a higher order of convergence. On the other hand, scheme UVS does not seem to achieve order of convergence when compared with a reference solution computed using scheme UV (although the errors are small, for instance smaller than the ones obtained with scheme UV-AD). But it achieves optimal orders when the study is performed using a reference solution computed using scheme UVS itself with $h=10^{-5}$ . These results make evident that the scheme UVS is well behaved although it needs very small discretization parameters to achieve exactly the same solution than the one obtained by the other schemes.

Scheme UV
h	$e(u)$	$r(u)$	$e(v)$	$r(v)$	$e(v_{x})$	$r(v_{x})$
1/200	$0.0344$	$-$	$0.0046$	$-$	$0.7893$	$-$
1/400	$0.0079$	$2.1171$	$0.0012$	$1.9859$	$0.2038$	$1.9536$
1/600	$0.0035$	$2.0093$	$0.0005$	$1.9943$	$0.0904$	$2.0052$
1/800	$0.0020$	$1.9916$	$0.0003$	$1.9968$	$0.0508$	$2.0037$
1/1000	$0.0013$	$2.0000$	$0.0002$	$1.9980$	$0.0325$	$2.0025$
Scheme UV-ND
1/200	$0.0355$	$-$	$0.0032$	$-$	$0.6240$	$-$
1/400	$0.0119$	$1.5785$	$0.0009$	$1.8814$	$0.1772$	$1.8162$
1/600	$0.0060$	$1.6656$	$0.0004$	$1.8811$	$0.0852$	$1.8049$
1/800	$0.0037$	$1.7198$	$0.0002$	$1.9069$	$0.0500$	$1.8564$
1/1000	$0.0025$	$1.8046$	$0.0001$	$1.9353$	$0.0327$	$1.8958$
Scheme UV-NS
1/200	$0.0298$	$-$	$0.00400$	$-$	$0.7671$	$-$
1/400	$0.0074$	$2.0132$	$0.00099$	$2.0069$	$0.1946$	$1.9787$
1/600	$0.0032$	$2.0473$	$0.00044$	$1.9965$	$0.0859$	$2.0163$
1/800	$0.0018$	$1.9928$	$0.00024$	$1.9981$	$0.0482$	$2.0071$
1/1000	$0.0012$	$1.9987$	$0.00015$	$1.9991$	$0.0308$	$2.0044$
Scheme UV-AD
1/200	$0.1032$	$-$	$0.0232$	$-$	$1.4694$	$-$
1/400	$0.0558$	$0.8865$	$0.0124$	$0.9029$	$0.8023$	$0.8730$
1/600	$0.0384$	$0.9197$	$0.0085$	$0.9405$	$0.5568$	$0.9010$
1/800	$0.0294$	$0.9321$	$0.0064$	$0.9569$	$0.4275$	$0.9185$
1/1000	$0.0238$	$0.9388$	$0.0052$	$0.9662$	$0.3472$	$0.9328$
Scheme UVS
Reference solution computed with scheme UV
1/200	$0.0782$	$-$	$0.0071$	$-$	$0.9256$	$-$
1/400	$0.0311$	$1.3284$	$0.0037$	$0.9605$	$0.2546$	$1.8621$
1/600	$0.0258$	$0.4587$	$0.0032$	$0.3135$	$0.1353$	$1.5588$
1/800	$0.0247$	$0.1545$	$0.0031$	$0.1274$	$0.0998$	$1.0568$
1/1000	$0.0244$	$0.0587$	$0.0031$	$0.0580$	$0.0868$	$0.6286$
Scheme UVS
Reference solution computed with scheme UVS
1/200	$0.0757$	$-$	$0.0063$	$-$	$0.9193$	$-$
1/400	$0.0195$	$1.9564$	$0.0018$	$1.8221$	$0.2368$	$1.9565$
1/600	$0.0087$	$1.9990$	$0.0008$	$1.9479$	$0.1059$	$1.9859$
1/800	$0.0048$	$2.0317$	$0.0005$	$2.0183$	$0.0599$	$1.9800$
1/1000	$0.0031$	$1.9533$	$0.0003$	$1.9576$	$0.0384$	$1.9959$

Table 2: Example IV: Experimental absolute errors and order of convergences for schemes UV, UV-ND, UV-NS, UV-AD and UVS

5 Conclusions

In this work we have proposed and studied numerical schemes to approximate a chemo-attraction and consumption model, and we have compared them numerically in $1D$ domains. We have focused on designing schemes in such a way that they maintain the main properties of the continuous problem at the discrete level. In fact, the three more challenging properties that we have identified and focused on are: $(a)$ positivity, $(b)$ dissipative energy law (6) and $(c)$ an estimate of a singular functional (5). We have developed several schemes:

•

Schemes UV-ND and UV-NS. Satisfying discrete versions of $(a)$ and $(c)$ .
•

Scheme UVS. Satisfying a discrete version of $(a)$ and $(b)$ .

We have compared these schemes with a straightforward Scheme UV and with Scheme UV-AD (an upwind Finite Volume scheme known to satisfy $(a)$ and that it can be reinterpreted as a Finite Element scheme with artificial numerical dissipation).

The numerical tests reported in this work have been designed to check how well the schemes behave with respect to positivity approximation, the approximation of the evolution in time of the energy and the behavior of the error as the size of the spatial mesh is reduced (error estimates).

Scheme UVS has been designed to satisfy the energy property $(b)$ , having to rewrite the original system introducing regularization terms that to be properly approximated they might need small choices of the discrete parameters. The numerical results illustrate that the scheme perform well in all the tests if the discretization parameters are small enough (as small as the one needed for the rest of the schemes), although it is not able to exactly capture the decreasing evolution of the energy $E(u,v)$ . This fact it is not a surprise because the scheme UVS is designed to satisfy a relation for a modified energy $E(u,{\boldsymbol{\sigma}})$ , which will only approximates well to $E(u,v)$ for very small discretization parameters. The other four schemes (UV, UV-ND, UV-NS and UV-AD) performs well in the three proposed tests but we need to remark that schemes UV, UV-ND and UV-NS performs much better than UV-AD in the approximation test, because schemes UV, UV-ND and UV-NS show second order convergence for all the unknowns, while UV-AD achieves only first order (due to the artificial dissipation introduced to maintain the positivity). Moreover, we have detected that some spurious oscillations might appear when using scheme UV-AD with a not very small time step, due to the fact that upwind schemes can have problems when the transport velocity is not incompressible. Although both schemes UV-ND and UV-NS perform equally well in the numerical tests, it is worth to mention that scheme UV-ND is much more flexible than scheme UV-NS, because the latter needs the requirement of considering structured meshes in order to satisfy the derived properties.

Finally, this work emphasizes two interesting points that can be extended to other problems where there are dissipative energy laws and maximum/minimum principles involved:
(I) For numerical schemes that satisfy a energy stability property with respect to a modified energy, it is important to keep in mind that even if in the continuous case the two energies are equivalent, the discrete versions will coincide only if the discretization parameters are really small.
(II) Finite Volume schemes achieve positivity by introducing artificial numerical dissipation in the system, but this dissipation can interfere with the accuracy of the scheme. Moreover, when a transport term appears with no incompressible velocity and the time step is not carefully chosen, the solutions can exhibit nonphysical oscillations.

References

[1] K. Baba and T. Tabata. On a conservative upwind finite-element scheme for convective diffusion equations. RAIRO Anal. Numér. 15 (1981) 3-25.
[2] J. W. Barrett and J. F. Blowey. Finite element approximation of a nonlinear cross-diffusion population model. Numer. Math. 98 (2004), no. 2, 195-221.
[3] N. Bellomo, A. Bellouquid, Y. Tao and M. Winkler. Toward a mathematical theory of Keller-Segel models of pattern formation in biological tissues. Mathematical Models and Methods in Applied Sciences, 25 (2015) 1663-1763
[4] M. Bessemoulin-Chatard and A. Jüngel. A finite volume scheme for a Keller-Segel model with additional cross-diffusion. IMA J. Numer. Anal. 34 (2014) 96-122.
[5] P.G. Ciarlet and P.A. Raviart. Maximum principle and uniform convergence for the finite element method. Comput. Methods Appl. Mech. Engrg. 2 (1973) 17-31.
[6] P. De Leenheer, J. Gopalakrishnan and E. Zuhr. Nonnegativity of exact and numerical solutions of some chemotactic models. Comput. Math. Appl. 66 (2013) 356-375
[7] A. Duarte-Rodríguez, M.A. Rodríguez-Bellido, D.A. Rueda-Gómez and E.J. Villamizar-Roa. Numerical analysis for a Chemotaxis-Navier-Stokes system. ESAIM: M2AN (2020) doi: 10.1051/m2an/2020039
[8] L.C. Evans. Partial Differential Equations. American Mathematical Society. (2010)
[9] F. Filbet. A finite volume scheme for the Patlak-Keller-Segel chemotaxis model. Numer. Math. 104 (2006) 457-488
[10] F. Guillén-González and J.V. Gutiérrez-Santacreu. From a cell model with active motion to a Hele-Shaw-like system: a numerical approach. Numerische Mathematik, 143 (2019) 107-137.
[11] F. Guillén-González, M.A. Rodríguez-Bellido and D.A. Rueda-Gómez. Study of a chemo-repulsion model with quadratic production. Part II: Analysis of an unconditional energy-stable fully discrete scheme. Computers and Mathematics with Applications, 80 (2020) 636-652.
[12] F. Guillén-González, M.A. Rodríguez-Bellido and D.A. Rueda-Gómez. Unconditionally energy stable fully discrete schemes for a chemo-repulsion model. Mathematics of Computation 88 (2019) 2069-2099
[13] F. Guillén-González and G. Tierra. Energy-stable and boundedness preserving numerical schemes for the Cahn-Hilliard equation with degenerate mobility. Submitted arXiv:2301.04913.
[14] D. Horstmann. From 1970 until present: the Keller-Segel model in chemotaxis and its consequences I. Jahresber. Deutsch. Math.-Verein. 105 (2003) 103-165
[15] D. Horstmann. From 1970 until present: the Keller-Segel model in chemotaxis and its consequences. II. Jahresber. Deutsch. Math.-Verein. 106 (2004) 51-69
[16] E.F. Keller and L.A. Segel. Traveling bands of chemotactic bacteria: a theoretical analysis. J. Theoret. Biol. 30 (1971) 377-380
[17] A. Marrocco. Numerical simulation of chemotactic bacteria aggregation via mixed finite elements. M2AN Math. Model. Numer. Anal. 37 (2003) 617-630
[18] MATLAB R2019b version 9.7.0. Natick, Massachusetts: The MathWorks Inc. (2019)
[19] N. Saito. Error analysis of a conservative finite-element approximation for the Keller-Segel system of chemotaxis. Commun. Pure Appl. Anal. 11 (2012) 339-364
[20] Y. Tao. Boundedness in a chemotaxis model with oxygen consumption by bacteria. J. Math. Anal. Appl. 381 (2011) 521-529
[21] Y. Tao and M. Winkler. Eventual smoothness and stabilization of large-data solutions in a three-dimensional chemotaxis system with consumption of chemoattractant. J. Differential Equations 252 (2012) 2520-2543
[22] I. Tuval, L. Cisneros, C. Dombrowski, C.W. Wolgemuth, J.O. Kessler and R.E. Goldstein. Bacterial swimming and oxygen transport near contact lines. P.N.A.S. 102 (2005) 2277-2282
[23] M. Winkler. Global large-data solutions in a chemotaxis Navier-Stokes system modeling cellular swimming in fluid drops. Communications in Partial Differential Equations 37 (2012) 319-351
[24] M. Winkler. The two-dimensional Keller-Segel system with singular sensitivity and signal absorption: Global large-data solutions and their relaxation properties. Mathematical Models and Methods in Applied Sciences 26 (2016) 987-1024
[25] G. Zhou and N. Saito. Finite volume methods for a Keller-Segel system: discrete energy, error estimates and numerical blow-up analysis. Numer. Math. 135 (2017) 265-311

Finite Element numerical schemes for a chemo-attraction and consumption model††thanks: This research work has been supported by Proyecto PGC2018-098308-B-I00, financed by FEDER / Ministerio de Ciencia e Innovación - Agencia Estatal de Investigación, Spain.

Abstract

1 Introduction

2 The model

2.1 Properties of the model

Remark 2.1.

2.2 Reformulation of the problem. The (u,v,𝝈)(u,v,{\boldsymbol{\sigma}}) problem

Theorem 2.2.

Proof.

Corollary 2.3.

Proof.

3 Numerical Schemes

3.1 UV-schemes

3.1.1 Scheme UV

Theorem 3.1 (un+1u^{n+1}-problem).

Proof.

Theorem 3.2 (vn+1v^{n+1}-problem).

Proof.

3.1.2 Scheme UV-ND (Nonlinear Diffusion)

Theorem 3.3.

Proof.

Corollary 3.4.

Proof.

Remark 3.5.

3.1.3 Scheme UV-NS (Nonlinear Sensitivity)

Remark 3.6 (How to construct Λεi​(u)\Lambda^{i}_{\varepsilon}(u)).

Theorem 3.7.

Proof.

Remark 3.8.

3.2 Scheme UVS

Theorem 3.9.

Proof.

Corollary 3.10.

Proof.

Corollary 3.11.

Proof.

Remark 3.12.

4 Numerical simulations in 1​D1D domains

4.1 Iterative Methods

4.1.1 Scheme UV-ND (Nonlinear Diffusion)

4.1.2 Scheme UV-NS (Nonlinear Sensitivity)

4.1.3 Scheme UVS

4.2 Scheme UV-AD (Artificial Diffusion) in 1​D1D

4.3 Example I: Dynamic towards constants

4.4 Example II: Positivity test

4.5 Example III: Energy stability test

4.6 Example IV: Approximation and error estimates test

5 Conclusions

References

Finite Element numerical schemes for a chemo-attraction and consumption model^†^†thanks: This research work has been supported by Proyecto PGC2018-098308-B-I00, financed by FEDER / Ministerio de Ciencia e Innovación - Agencia Estatal de Investigación, Spain.

2.2 Reformulation of the problem. The $(u,v,{\boldsymbol{\sigma}})$ problem

Theorem 3.1 ( $u^{n+1}$ -problem).

Theorem 3.2 ( $v^{n+1}$ -problem).

Remark 3.6 (How to construct $\Lambda^{i}_{\varepsilon}(u)$ ).

4 Numerical simulations in $1D$ domains

4.2 Scheme UV-AD (Artificial Diffusion) in $1D$