\pagespan

118 \copyrightinfo2021

ERROR ESTIMATES AND BLOW-UP ANALYSIS OF
A FINITE-ELEMENT APPROXIMATION FOR THE
parabolic-elliptic Keller-Segel system

Wenbin Chen Shanghai Key Laboratory for Contemporary Applied Mathematics, School of Mathematical Sciences, Fudan University, Shanghai, 200433, China [email protected] Qianqian Liu School of Mathematical Sciences, Fudan University, Shanghai, 200433, China [email protected] Jie Shen Department of Mathematics, Purdue University, West Lafayette, IN 47907, USA [email protected]

Abstract.

The Keller-Segel equations are widely used for describing chemotaxis in biology. Recently, a new fully discrete scheme for this model was proposed in [46], mass conservation, positivity and energy decay were proved for the proposed scheme, which are important properties of the original system. In this paper, we establish the error estimates of this scheme. Then, based on the error estimates, we derive the finite-time blowup of nonradial numerical solutions under some conditions on the mass and the moment of the initial data.

Key words and phrases:

parabolic-elliptic systems, finite element method, error estimates, finite-time blowup.

1991 Mathematics Subject Classification:

65M12, 35K61, 35K55, 92C17

1. Introduction.

Keller and Segel first proposed a nonlinear model in the 1970s to describe the effect of cell aggregation in [27, 28]. A simplified Keller-Segel model in 2-D is given by

	$\displaystyle\frac{\partial u}{\partial t}=$	$\displaystyle\Delta u-\chi\nabla\cdot(u\nabla v),\quad x\in\Omega,\ t>0,$		(1.1)
	$\displaystyle 0=$	$\displaystyle\Delta v-v+\alpha u,\quad x\in\Omega,\ t>0,$		(1.2)

where $\Omega\subset\mathbb{R}^{2}$ is a bounded domain with smooth boundary $\partial\Omega$ . The unknown $u=u(x,t)$ and $v=v(x,t)$ represent the concentration of the organism and chemoattractant respectively. The parameters $\chi,\alpha$ are positive constants with $\chi$ being the sensitivity of chemotaxis. The model is supplemented with initial conditions

u(x,t=0)=u_{0}(x),\ v(x,t=0)=v_{0}(x),\ x\in\Omega,

and no flux boundary conditions

\frac{\partial u}{\partial\bm{n}}-\chi u\frac{\partial v}{\partial\bm{n}}=0,\ \frac{\partial v}{\partial\bm{n}}=0,\quad x\in\partial\Omega,\ t>0,

where $\bm{n}$ denotes the unit outward normal vector to the boundary $\partial\Omega$ , $\partial/\partial\bm{n}$ represents differentiation along $\bm{n}$ on $\partial\Omega$ .

A different version of the Keller-Segel model consists in replacing (1.2) by

\frac{\partial v}{\partial t}=\Delta v-v+\alpha u,\quad x\in\Omega,\ t>0.

(1.3)

The equation (1.1) describes the motion of the organism $u$ . The term $F=-\nabla u+\chi u\nabla v$ is the flux, and the effect of diffusion $-\Delta u$ and that of chemotaxis $\chi\nabla\cdot(u\nabla v)$ are competing for $u$ to vary. The equation (1.2) describes the change in concentration of the chemoattractant $v$ , it is influenced by the diffusion and the decay of the chemoattractant as well as the growth of the organism. In general, the chemoattractant particles are much smaller than the organism particles, thus it diffuses faster, which means that the diffusion of the chemoattractant will reach the equilibrium state in a relatively short time. The model (1.1)-(1.2) is called parabolic–elliptic system. On the other hand, (1.1) with (1.3) is a parabolic–parabolic system.

The solution of the Keller-Segel model (1.1)-(1.2) has several well-known properties, particularly, it may blow up in finite time. Various aspects and results for the classical Keller-Segel model since 1970, along with some open questions, are summarized in [25]. Positivity, mass conservation and energy dissipation of Keller-Segel equations can be found in [35],[36],[47],[29] and [6], which plays an important role to study the Keller-Segel system. Blanchet, Dolbeault and Perthame presented in [3] a detail proof of the existence of weak solutions when the initial mass is below the critical mass, above which any solution to the parabolic-elliptic systems blows up in finite time in the whole Euclidean space. In [37], Nagai demonstrated the finite-time blowup of nonradial solutions under some assumptions on the mass and the moment of the initial data. As for the parabolic-parabolic systems, Blanchet proved in [2] the optimal critical mass of the solutions in $\mathbb{R}^{d}$ with $d\geq 3$ . Wei proved that for every nonnegative initial data in $L^{1}(R^{2})$ , the 2-D Keller-Segel equation is globally well-posed if and only if the total mass $M\leq 8\pi$ in [49].

Although the large time behavior of the solution of the Keller-Segel model (1.1)-(1.2) has been well studied, there is still much to explore on the numerical side. Since the Keller-Segel equations possess three important properties: positivity, mass conservation and energy dissipation, it is preferable that numerical schemes can preserve these properties at the discrete level. In [26], the existence of weak solutions and upper bounds for the blow-up time for time-discrete (including the implicit Euler, BDF and Runge-Kutta methods) approximations of the parabolic-elliptic Keller-Segel models in the two-dimensional whole space are established. Liu, Li and Zhou proposed a numerical method in [34] which preserves both positivity and asymptotic limit, the proposed numerical method does not generate negative density if initialized properly under a less strict stability condition. Saito and Suzuki presented a finite difference scheme in [42] which satisfies the conservation of a discrete $L^{1}$ norm.

Some finite element methods are proposed in previous works. Saito presented a finite element scheme for parabolic-elliptic systems in [43] that satisfies both positivity and mass conservation properties. Under some assumptions on the regularity of solutions, the error estimates were established. Saito further constructed the finite element methods to the parabolic-parabolic systems in [44] and derived error analysis by using analytical semigroup theory. Gurusamy and Balachandran proposed a finite element method for parabolic-parabolic systems and established the existence of approximate solutions by using Schauder’s fixed point theorem in [23]. Further the error estimates for the approximate solutions in $H^{1}$ -norm were derived.

The discontinuous Galerkin methods can be also used to solve the Keller-Segel equations. Epshteyn and Kurganov developed a family of new interior penalty discontinuous Galerkin methods and proved error estimates for the proposed high-order discontinuous Galerkin methods in [15]. Epshteyn and Izmirlioglu further constructed a discontinuous Galerkin method for Keller-Segel model in [16] and obtained fully discrete error estimates for the proposed scheme. In 2017, Li, Shu and Yang applied the local discontinuous Galerkin (LDG) method to 2D Keller–Segel chemotaxis model in [30], they improved the results upon [15] and gave optimal rate of convergence under special finite element spaces before the blow-up occurs. In 2019, Guo, Li and Yang constructed a consistent numerical energy and prove the energy dissipation with the LDG discretization in [22].

Another important numerical methods for Keller-Segel models are finite volume methods since the positivity property can be naturally preserved. Filbet proposed in [18] a finite volume scheme for the parabolic-elliptic system, and by assuming the CFL condition $\chi\Delta t\mathcal{D}_{\mathcal{T},1}<1$ and the initial datum $n^{0}\geq a^{0}>0$ , he proved existence and uniqueness of the numerical solution by using the Browder fixed point theorem, and showed that the numerical approximation converges to the exact solution under some assumptions. In 2016, Zhou and Saito proposed a finite volume scheme in [52], and established error estimates in $L^{p}$ norm with a suitable $p>2$ for the two dimensional case under some regularity assumptions of solutions and admissible mesh. By focusing on the radially symmetirc solution, they derived some a prior estimates to study the blow-up phenomenon of numerical solution.

There have been growing interests in positivity-preserving analysis for gradient flows with logarithmic energy potential. Some theoretical analysis of the positivity-preserving property and the energy stability have been explored for these numerical schemes for certain systems, such as Cahn-Hilliard systems in [7, 10, 11, 51, 12], the Poisson-Nernst-Planck-Cahn-Hilliard systems in [41], the Poisson-Nernst-Planck systems in [32], the thin film model without slope selection in [31] and a structure-preserving, operator splitting scheme for reaction-diffusion systems in [33]. The techniques of the higher order consistency analysis combined with rough error estimate and refined one have been presented in [32, 13, 14] which will be utilized in the following to obtain the convergence analysis.

Recently, a new approach for constructing positivity preserving schemes was proposed in [46]. The key for this approach is to write $\Delta u$ as $\nabla\cdot(u\nabla\log u)$ in (1.1), and then use a convex splitting idea to construct mass conservative, bound preserving, and uniquely solvable schemes for (1.1)-(1.2) and for (1.1)-(1.3). The main purposes of this paper are to establish the convergence of the fully discrete scheme proposed in [46], and to show the finite-time blowup of numerical solutions under some conditions on the mass and moments of the initial data. More precisely, let $u_{h}^{k}$ be an approximation of $u(\cdot,k\tau)$ , where $\tau>0$ is the time step and $k\in\mathbb{N}$ . Let $\theta=(u_{0},1)$ be the initial mass, and $M^{k}=(u_{h}^{k},\varphi)_{h}$ be the moment of $u_{h}^{k}$ . Our first goal is to establish the error estimates for the fully discrete scheme proposed in [46] (cf. Theorem 6). Another important feature of the Keller-Segel system (1.1)-(1.2) is that the solution may blow up in finite time under certain conditions on the initial data. Our second goal is to show that the numerical solution will also blow up in finite time under similar conditions on the initial data (cf. Theorem 15). Many previous works (see [42, 43, 46]) show that the numerical solution seems to blow up under large initial data by several numerical experiments. However, there is still much to explore on the theoretical proof of blowup phenomenon besides the radial numerical solution in [52] mentioned before.

The rest of the paper is organized as follows. In Section 2, we recall some properties of the classical Keller-Segel equations, including its finite-time blowup behavior. In Section 3, we introduce the fully discrete scheme constructed in [46] and carry out a rigorous error analysis. In Section 4, we show that the numerical solution will blow up in finite time under suitable conditions on the initial data.

2. The Keller-Segel equations

In this section, we recall some properties for the Keller-Segel system (1.1)-(1.2) with no flux boundary conditions. In addition, we assume the initial value $u_{0}\in W^{2,p}(\Omega),\ 1<p<\infty$ , and satisfies

u_{0}\geq 0\quad\text{and}\quad u_{0}\not\equiv 0,\quad\forall x\in\Omega.

It was shown in [35] that there exist some $T>0$ such that (1.1)-(1.2) is well posed in the time interval $[0,T]$ . Moreover, it holds that

Theorem 1.

The Keller-Segel system (1.1)-(1.2) satisfies the following properties:

(i)

Positivity preserving:

$u(x,t)>0,\quad v(x,t)>0\quad\text{on}\quad\Omega\times(0,T].\\$

In fact, it is a consequence of the strong maximum principle [8] .

(ii)

Mass conservation:

\int_{\Omega}u(x,t)\text{d}x=\int_{\Omega}u_{0}(x)\text{d}x,\quad\text{for all}\quad t>0.\\

It is immediately follows from

\dfrac{\text{d}}{\text{d}t}\int_{\Omega}u(x,t)\text{d}x=\int_{\partial\Omega}\left(\dfrac{\partial u}{\partial\bm{n}}-u\dfrac{\partial v}{\partial\bm{n}}\right)\text{d}S=0\ .

As a consequence of (i) and (ii), we obtain the conservation of $L^{1}$ norm, namely

\|u(t)\|_{L^{1}(\Omega)}=\|u_{0}\|_{L^{1}(\Omega)}.

(iii)

Energy decay:

\dfrac{\text{d}F[u(t),v(t)]}{\text{d}t}=-\int_{\Omega}u|\nabla\cdot(\log u-\chi v)|^{2}\text{d}x\leq 0,

where the free energy of (1.1)-(1.2) is defined by

F[u,v]=\int_{\Omega}\left(u(\log u-1)-\chi uv+\dfrac{\chi}{2\alpha}|\nabla v|^{2}+\dfrac{\chi}{2\alpha}v^{2}\right)\text{d}x\ .

The following result is shown in [37].

Lemma 2.

[37] Let $q\in\Omega$ and $0<r_{1}<r_{2}<dist(q,\partial\Omega)$ , where $dist(q,\partial\Omega)$ is the distance between $q$ and $\partial\Omega$ . Then there exist positive constants $C_{1},C_{2}$ depending only on $r_{1},r_{2}$ and $dist(q,\partial\Omega)$ such that for $t\in(0,T]$ ,

	$\displaystyle\dfrac{\text{d}}{\text{d}t}\int_{\Omega}u(x,t)\Phi(x)\text{d}x$	(2.1)
$\displaystyle\leq$	$\displaystyle 4\int_{\Omega}u_{0}(x)\text{d}x-\dfrac{\alpha\chi}{2\pi}\left(\int_{\Omega}u_{0}(x)\text{d}x\right)^{2}+C_{1}\left(\int_{\Omega}u_{0}(x)\text{d}x\right)\left(\int_{\Omega}u(x,t)\Phi(x)\text{d}x\right)$
	$\displaystyle+C_{2}\left(\int_{\Omega}u_{0}(x)\text{d}x\right)^{3/2}\left(\int_{\Omega}u(x,t)\Phi(x)\text{d}x\right)^{1/2},$

where $\Phi(x)=\phi(|x-q|)$ with

\displaystyle\phi(r)=\begin{cases}r^{2}&if\quad 0\leq r\leq r_{1},\\ a_{1}r^{2}+a_{2}r+a_{3}&if\quad r_{1}<r\leq r_{2},\\ r_{1}r_{2}&if\quad r>r_{2},\\ \end{cases}

where $a_{1}=-\dfrac{r_{1}}{r_{2}-r_{1}},\quad a_{2}=\dfrac{2r_{1}r_{2}}{r_{2}-r_{1}},\quad a_{3}=-\dfrac{r_{1}^{2}r_{2}}{r_{2}-r_{1}}$ .

The finite-time blowup behavior is then proved using the above result.

Theorem 3.

[37] Assume that $\int_{\Omega}u_{0}(x)\text{d}x>8\pi/(\alpha\chi)$ , if $\int_{\Omega}u_{0}(x)|x-q|^{2}\text{d}x$ is sufficiently small, then the solution $(u,v)$ to (1.1)-(1.2) blows up in finite time.

Moreover, the following pointwise estimates for $v$ is established in [19]. An application of the Neumann semigroup leads to

	$\displaystyle v(x,t)$	$\displaystyle\geq\left(\int_{0}^{\infty}\dfrac{1}{(4\pi t)^{\frac{n}{2}}}e^{-(t+\frac{(\mathrm{diam}\Omega)^{2}}{4t})}\mathrm{d}t\right)\int_{\Omega}u(x,t)\mathrm{d}x$
		$\displaystyle=\\|u_{0}\\|_{L^{1}(\Omega)}\int_{0}^{\infty}\dfrac{1}{(4\pi t)^{\frac{n}{2}}}e^{-(t+\frac{(\mathrm{diam}\Omega)^{2}}{4t})}\mathrm{d}t,$

for all $x\in\Omega,\ t\in(0,T]$ , whenever $(u,v)$ solves (1.1)-(1.2) in $\Omega\times(0,T]$ for some $T>0$ .

In this paper, we assume that

u(x,t)\geq\epsilon_{0}\quad\text{for some }\epsilon_{0}>0,\quad(x,t)\in\overline{\Omega}\times(0,T].

3. The fully discrete scheme and error estimates

In this section, we describe the fully discrete scheme in [46] for (1.1)-(1.2), construct the error equations and establish the error estimates.

We now give a precise description of our finite element space $X_{h}$ . Given a triangulation $\mathcal{T}$ for $\Omega$ , we let $Z_{h}$ consists of all the vertices excluding those where Dirichlet boundary conditions are prescribed. We define $X_{h}$ to be the finite element space spanned by the piecewise linear continuous functions based on $\mathcal{T}$ . Let $e$ be a triangle of the triangulation $\mathcal{T}$ , and $P_{e,i},i=1,2,3$ be its vertices, we define the quadrature formula

Q_{e}(f)=\dfrac{1}{3}\text{area}(e)\sum_{i=1}^{3}f(P_{e,i})\approx\int_{e}f\mathrm{d}x.

We recall that [48]

\Big{|}Q_{e}(f)-\int_{e}f\mathrm{d}x\Big{|}\leq Ch^{2}\sum_{|\alpha|=2}\|D^{\alpha}f\|_{L^{1}(e)}.

We then define the discrete inner product in $X_{h}$ by

(u,v)_{h}=\sum_{e\in\mathcal{T}}Q_{e}(uv),

the corresponding norm is defined by $\|\cdot\|_{L_{h}^{2}}$ . We have the following estimates in [48] for the quadrature error.

Lemma 4.

Let $\epsilon_{h}(\cdot,\cdot)=(\cdot,\cdot)_{h}-(\cdot,\cdot)$ denote the quadrature error, then we have

|\epsilon_{h}(u,v)|\leq Ch^{2}\|\nabla u\|_{L^{2}(\Omega)}\|\nabla v\|_{L^{2}(\Omega)},\ \forall\ u,v\in X_{h}.

Applying Lemma 4, the norm $\|\cdot\|_{L_{h}^{2}}$ has the following property

\displaystyle\|\nabla u\|_{L_{h}^{2}(\Omega)}=\|\nabla u\|_{L^{2}(\Omega)},\ \forall\ u\in X_{h}.

(3.1)

Let $I_{h}:\;C(\Omega)\rightarrow X_{h}$ be the Lagrange interpolation operator, which has the approximation property [4] that for all $g\in H^{2}(\Omega)\cap H_{0}^{1}(\Omega)$ ,

\|I_{h}g-g\|_{L^{2}(\Omega)}\leq Ch^{2}\|g\|_{H^{2}(\Omega)}\quad\text{and}\quad\|\nabla(I_{h}g-g)\|_{L^{2}(\Omega)}\leq Ch\|g\|_{H^{2}(\Omega)}.

(3.2)

The fully discrete scheme proposed in [46] for (1.1)-(1.2) is to find $(u_{h}^{k+1},v_{h}^{k+1})\in X_{h}\times X_{h}$ such that for all $(\varphi,\psi)\in X_{h}\times X_{h}$ ,

	$\displaystyle\left(\overline{\partial}u_{h}^{k},\varphi\right)_{h}\ +\left(u_{h}^{k}\nabla\left(I_{h}\log u_{h}^{k+1}\right),\nabla\varphi\right)_{h}-\chi\left(u_{h}^{k}\nabla v_{h}^{k},\nabla\varphi\right)_{h}=0,$		(3.3)
	$\displaystyle(\nabla v_{h}^{k+1},\nabla\psi)+(v_{h}^{k+1},\psi)-\alpha(u_{h}^{k+1},\psi)_{h}=0,$		(3.4)

with the initial value $u_{h}^{0}:=I_{h}u_{0}$ . Here, the $(\cdot,\cdot)$ represents the usual $L^{2}$ inner product, the $(\cdot,\cdot)_{h}$ is the discrete inner product defined above, and $\overline{\partial}u_{h}^{k}$ is the forward Euler difference quotient approximating to $\partial_{t}u(t^{k})$ defined by

\overline{\partial}u_{h}^{k}=\frac{u_{h}^{k+1}-u_{h}^{k}}{\tau}.

In this setting, the authors in [46] proved the following.

Lemma 5.

[46] The numerical scheme (3.3)-(3.4) has the following properties:

(i)

Unique solvability:
The scheme (3.3)-(3.4) has a unique solution $(u_{h}^{k+1},v_{h}^{k+1})\in X_{h}\times X_{h}$ .
(ii)

Positivity preserving: If $u_{h}^{k}>0$ , then $u_{h}^{k+1}>0$ .
(iii)

Mass conservation:

$(u_{h}^{k+1},1)_{h}=(u_{h}^{k},1)_{h}=(u_{h}^{0},1)_{h},$

It is immediately derived by taking $\varphi=1$ in (3.3).

(iv)

Energy decay:

\tilde{F}^{k+1}-\tilde{F}^{k}\leq-\tau\left(u_{h}^{k}\nabla(I_{h}(\log u_{h}^{k+1}-\chi v_{h}^{k})),\nabla(I_{h}(\log u_{h}^{k+1}-\chi v_{h}^{k}))\right)_{h}\leq 0,

where the discrete energy of (3.3)-(3.4) is defined by

\tilde{F}^{k}[u,v]=(u_{h}^{k}(\log u_{h}^{k}-1),1)_{h}-\chi(u_{h}^{k},v_{h}^{k})_{h}+\dfrac{\chi}{2\alpha}\|\nabla v_{h}^{k}\|^{2}+\dfrac{\chi}{2\alpha}\|v_{h}^{k}\|^{2}.

We denote by $(u,v)$ the exact solution pair to the original equations (1.1)-(1.3), and all the upper bounds for the exact solution are denoted as $C$ . We set $u^{k}=u(t_{k}),\,v^{k}=v(t_{k})$ , and denote

e_{u}^{k}:=u^{k}-u_{h}^{k},\quad e_{v}^{k}:=v^{k}-v_{h}^{k},\quad\forall k\in\mathbb{N}.

The following theorem is the main result of this section.

Theorem 6.

Assume $u_{0}\in W^{2,p}(\Omega)(1<p<\infty)$ and the exact solution pair $(u,v)$ is smooth enough for a fixed final time $T>0$ . Then, provided $\tau$ and $h$ are sufficiently small and under the mild mesh-sizes requirement $\tau\leq Ch$ , we have the following error estimates

\|e_{u}^{m}\|_{L_{h}^{2}(\Omega)}+\|e_{v}^{m}\|_{H_{h}^{1}(\Omega)}+\Big{(}\tau\sum_{k=0}^{m-1}\|\nabla e_{u}^{k+1}\|_{L_{h}^{2}(\Omega)}^{2}\Big{)}^{1/2}\leq C(\tau+h),\quad\forall m\in\mathbb{N},

where $t_{m}=m\tau\leq T$ , $C>0$ is independent of $\tau$ and $h$ .

The proof for this theorem will be carried out with a sequence of procedures that we describe below.

Remark 7.

The mesh-sizes requirement $\tau\leq Ch$ in Theorem 6 is proposed to obtain a higher order consistency analysis via a perturbation argument, which is needed to get the separation property and the $W^{1,\infty}$ bound for the numerical solution.

3.1. Higher order consistent approximation to (3.3)-(3.4)

In this subsection, we apply the perturbation argument method in [32] to the finite element scheme to construct $f_{1},f_{2},f_{3}$ such that

\hat{u}:=u+hf_{1}+h^{2}f_{2}+h^{3}f_{3},

is consistent with the given numerical scheme (3.3)-(3.4) at the order $O(h^{4})$ . The following lemma is used to construct $f_{1},f_{2},f_{3}$ and the proof is given in Appendix.

By applying a perturbation argument, a higher order $O(h^{4})$ consistency is satisfied for $\hat{u}$ , which is needed to obtain the separation property and a $W^{1,\infty}$ bound for the numerical solution.

Lemma 8.

Suppose that $\tau\leq Ch$ and $u$ is smooth enough, then there exist bounded smooth functions $f_{1},f_{2},f_{3}$ , such that $\hat{u}=u+hf_{1}+h^{2}f_{2}+h^{3}f_{3}$ satisfies

\displaystyle\left(\overline{\partial}\hat{u}^{k},\varphi\right)_{h}+\left(\hat{u}^{k}\nabla I_{h}\log\hat{u}^{k+1},\nabla\varphi\right)_{h}-\chi\left(\hat{u}^{k}\nabla A_{h}\hat{u}^{k},\nabla\varphi\right)_{h}=\langle\mathcal{R}^{k}(\hat{u}),\varphi\rangle,

(3.5)

for all $\varphi\in X_{h}$ , $k\in\mathbb{N}$ , where $A_{h}=\alpha(-\Delta_{h}+I)^{-1}Q_{h}$ and $\langle\cdot,\cdot\rangle$ denotes the duality product satisfying

\displaystyle|\langle\mathcal{R}^{k}(\hat{u}^{k}),\varphi\rangle|\leq Ch^{4}\|\varphi\|_{H^{1}},

(3.6)

where $C$ depends on the regularity of the solution $u$ .

Remark 9.

Under the conditions that the exact solution $u\geq\epsilon_{0}$ for some $\epsilon_{0}>0$ , and $h$ is sufficiently small, we obtain that

\hat{u}\geq\dfrac{\epsilon_{0}}{2}.

(3.7)

Since the correction functions $f_{j},j=1,2,3$ only depend on the exact solution $u$ , they are bounded in $W^{1,\infty}$ norm. Then, we can obtain the following $W^{1,\infty}$ bound for $\hat{u}$ :

\|\hat{u}^{k}\|_{W^{1,\infty}}\leq C,\quad\forall k\geq 0.

(3.8)

3.2. A rough error estimate

In this subsection, we derive the strict separation property and a uniform $W^{1,\infty}$ bound for the numerical solution.

We recall the following inverse estimate in [4, p.111, Lemma 4.5.3].

Lemma 10.

Given a quasi-uniform triangulation $\mathcal{T}$ on domain $\Omega\subset\mathbb{R}^{n}$ , and $X_{h}$ be a finite-dimensional subspace of $W^{l,p}(K)\cap W^{m,q}(K)$ , where $1\leq p\leq\infty,1\leq q\leq\infty$ and $0\leq m\leq l$ . Then there exists a positive constant $C$ such that for all $u\in X_{h}$ , we have

\|u\|_{W^{l,p}(K)}\leq Ch^{m-l+n/p-n/q}\|u\|_{W^{m,q}(K)},

where $C$ is independent of $u$ .

We will also use the following discrete Gronwall inequality in [45, 50].

Lemma 11.

Assume that $\tau>0,B>0,\{a_{k}\},\{b_{k}\},\{\gamma_{k}\}$ are non-negative sequences such that

a_{m}+\tau\sum_{k=1}^{m}b_{k}\leq\tau\sum_{k=1}^{m-1}\gamma_{k}a_{k}+B,\ m\geq 1.

Then

a_{m}+\tau\sum_{k=1}^{m}b_{k}\leq B\mathrm{exp}\left(\tau\sum_{k=1}^{m-1}\gamma_{k}\right),\ m\geq 1.

Define an alternative error function:

\tilde{u}^{k}:=I_{h}\hat{u}^{k}-u_{h}^{k},\quad\forall k\in\mathbb{N}.

Subtracting the numerical scheme (3.3) from the consistency estimate (3.5) implies that

(\overline{\partial}\tilde{u}^{k},\varphi)_{h}=-(\tilde{u}^{k}\nabla\mathcal{V}_{u}^{k}+u_{h}^{k}\nabla\tilde{\mu}_{u}^{k},\nabla\varphi)_{h}+\langle\mathcal{R}^{k},\varphi\rangle,

(3.9)

where

\mathcal{V}_{u}^{k}:=I_{h}\log\hat{u}^{k+1}-\chi A_{h}\hat{u}^{k},

\tilde{\mu}_{u}^{k}:=I_{h}\log\hat{u}^{k+1}-I_{h}\log u_{h}^{k+1}-\chi A_{h}\tilde{u}^{k}.

Since $\mathcal{V}_{u}^{k}$ only depends on the exact solution, we can assume

\|\mathcal{V}_{u}^{k}\|_{W^{2,\infty}}\leq C.

(3.10)

Lemma 12.

The numerical solutions of the scheme (3.3)-(3.4) have the strict separation property and a uniform $W^{1,\infty}$ bound:

\displaystyle u_{h}^{k}\geq\dfrac{\epsilon_{0}}{4},\ \|u_{h}^{k}\|_{W^{1,\infty}}\leq C^{*},

for all $0\leq k\leq T/\tau$ , where $\epsilon_{0}$ and $C^{*}$ are positive constants.

Proof.

We shall first make the following assumption at the previous time step:

\|\tilde{u}^{k}\|_{L_{h}^{2}}\leq Ch^{15/4}.

(3.11)

Then, we will demonstrate that such an assumption will be recovered at the next time step in Section 3.3.

Using the inverse inequality, we obtain a $W^{1,\infty}$ bound for the numerical error function:

\|\tilde{u}^{k}\|_{\infty}\leq\dfrac{C\|\tilde{u}^{k}\|_{L_{h}^{2}}}{h}\leq\dfrac{Ch^{15/4}}{h}\leq Ch^{11/4}\leq 1,

(3.12)

\|\nabla\tilde{u}^{k}\|_{\infty}\leq\dfrac{C\|\tilde{u}^{k}\|_{\infty}}{h}\leq\dfrac{Ch^{11/4}}{h}\leq Ch^{7/4}\leq 1.

A combination of the above with (3.8), we get a $W^{1,\infty}$ bound for $u_{h}^{k}$ at the previous time step:

\|u_{h}^{k}\|_{\infty}\leq\|\hat{u}^{k}\|_{\infty}+\|\tilde{u}^{k}\|_{\infty}\leq C+1\leq C^{*},

\|\nabla u_{h}^{k}\|_{\infty}\leq\|\nabla\hat{u}^{k}\|_{\infty}+\|\nabla\tilde{u}^{k}\|_{\infty}\leq C+1\leq C^{*}.

Because of (3.12), taking $h$ sufficiently small, we have

\|\tilde{u}^{k}\|_{\infty}\leq Ch^{11/4}\leq\dfrac{\epsilon_{0}}{4}.

Then the strict separation property is valid for $u_{h}^{k}$ :

u_{h}^{k}\geq I_{h}\hat{u}^{k}-\|\tilde{u}^{k}\|_{\infty}\geq\dfrac{\epsilon_{0}}{4}.

(3.13)

Taking $\varphi=\tau\tilde{\mu}_{u}^{k}$ in (3.9) leads to

(\tilde{u}^{k+1},\tilde{\mu}_{u}^{k})_{h}+\tau(u_{h}^{k}\nabla\tilde{\mu}_{u}^{k},\nabla\tilde{\mu}_{u}^{k})_{h}=(\tilde{u}^{k},\tilde{\mu}_{u}^{k})_{h}-\tau(\tilde{u}^{k}\nabla\mathcal{V}_{u}^{k},\nabla\tilde{\mu}_{u}^{k})_{h}+\tau\langle\mathcal{R}^{k},\tilde{\mu}_{u}^{k}\rangle.

(3.14)

Now we deal with the left hand side of (3.14). To proceed the first term on the left hand side of (3.14), notice that $((\tilde{u}^{k+1})^{3},1)_{h}=0$ since $(\tilde{u}^{k+1},1)_{h}=(I_{h}\hat{u}^{0}-u_{h}^{0},1)_{h}=0$ stated in Remark 20, using the Hölder inequality, we have

	$\displaystyle(\tilde{u}^{k+1},\tilde{\mu}_{u}^{k})_{h}=(\tilde{u}^{k+1},I_{h}\log\hat{u}^{k+1}-I_{h}\log u_{h}^{k+1})_{h}+(\tilde{u}^{k+1},-\chi A_{h}\tilde{u}^{k})_{h}$	(3.15)
$\displaystyle\geq$	$\displaystyle\Big{(}\tilde{u}^{k+1},\dfrac{1}{\hat{u}^{k+1}}\tilde{u}^{k+1}\Big{)}_{h}+\frac{1}{2\\|\zeta\\|_{\infty}^{2}}\Big{(}(\tilde{u}^{k+1})^{3},1\Big{)}_{h}-C^{\prime}\\|\tilde{u}^{k+1}\\|_{L_{h}^{2}}\\|\tilde{u}^{k}\\|_{L_{h}^{2}}$
$\displaystyle\geq$	$\displaystyle\dfrac{1}{C}\\|\tilde{u}^{k+1}\\|_{L_{h}^{2}}^{2}-\dfrac{1}{C}\\|\tilde{u}^{k+1}\\|_{L_{h}^{2}}^{2}-C^{\prime\prime}\\|\tilde{u}^{k}\\|_{L_{h}^{2}}^{2}$
$\displaystyle=$	$\displaystyle-C^{\prime\prime}\\|\tilde{u}^{k}\\|_{L_{h}^{2}}^{2},$

where $\zeta$ lies between $\hat{u}^{k+1}$ and $u_{h}^{k+1}$ , and (3.8) has been utilized in the second inequality. As for the second term on the left hand side of (3.14), using the strict separation property of the numerical solution (3.13), we have

(u_{h}^{k}\nabla\tilde{\mu}_{u}^{k},\nabla\tilde{\mu}_{u}^{k})_{h}\geq\dfrac{\epsilon_{0}}{4}\|\nabla\tilde{\mu}_{u}^{k}\|_{L_{h}^{2}}^{2}.

(3.16)

Next, we deal with the right hand side of (3.14). We apply the Hölder inequality and the Young inequality:

	$\displaystyle\|(\tilde{u}^{k},\tilde{\mu}_{u}^{k})_{h}\|$	$\displaystyle\leq\\|\tilde{u}^{k}\\|_{-1}\\|\nabla\tilde{\mu}_{u}^{k}\\|_{L_{h}^{2}}\leq C\\|\tilde{u}^{k}\\|_{L_{h}^{2}}\\|\nabla\tilde{\mu}_{u}^{k}\\|_{L_{h}^{2}}$		(3.17)
		$\displaystyle\leq\dfrac{6C}{\epsilon_{0}\tau}\\|\tilde{u}^{k}\\|_{L_{h}^{2}}^{2}+\dfrac{\epsilon_{0}\tau}{24}\\|\nabla\tilde{\mu}_{u}^{k}\\|_{L_{h}^{2}}^{2}.$		(3.17)

An application of the Cauchy-Schwarz inequality and (3.10) leads to

	$\displaystyle\|-(\tilde{u}^{k}\nabla\mathcal{V}_{u}^{k},\nabla\tilde{\mu}_{u}^{k})_{h}\|$	$\displaystyle\leq\\|\tilde{u}^{k}\nabla\mathcal{V}_{u}^{k}\\|_{L_{h}^{2}}\\|\nabla\tilde{\mu}_{u}^{k}\\|_{L_{h}^{2}}\leq\\|\nabla\mathcal{V}_{u}^{k}\\|_{\infty}\\|\tilde{u}^{k}\\|_{L_{h}^{2}}\\|\nabla\tilde{\mu}_{u}^{k}\\|_{L_{h}^{2}}$
		$\displaystyle\leq\dfrac{6C}{\epsilon_{0}}\\|\tilde{u}^{k}\\|_{L_{h}^{2}}^{2}+\dfrac{\epsilon_{0}}{24}\\|\nabla\tilde{\mu}_{u}^{k}\\|_{L_{h}^{2}}^{2}.$

Using the inequality (3.6), we have

|\langle\mathcal{R}^{k},\tilde{\mu}_{u}^{k}\rangle|\leq Ch^{4}\|\nabla\tilde{\mu}_{u}^{k}\|_{2}\leq\dfrac{6C}{\epsilon_{0}}h^{8}+\dfrac{\epsilon_{0}}{24}\|\nabla\tilde{\mu}_{u}^{k}\|_{L_{h}^{2}}^{2}.

(3.18)

Substitution of (3.15)-(3.18) into (3.14) leads to

\dfrac{\epsilon_{0}}{4}\tau\|\nabla\tilde{\mu}_{u}^{k}\|_{L_{h}^{2}}^{2}-C^{\prime\prime}\|\tilde{u}^{k}\|_{L_{h}^{2}}^{2}\leq\dfrac{6C}{\epsilon_{0}\tau}\|\tilde{u}^{k}\|_{L_{h}^{2}}^{2}+\dfrac{6C\tau}{\epsilon_{0}}\|\tilde{u}^{k}\|_{L_{h}^{2}}^{2}+\dfrac{6C}{\epsilon_{0}}h^{8}\tau+\dfrac{\epsilon_{0}\tau}{8}\|\nabla\tilde{\mu}_{u}^{k}\|_{L_{h}^{2}}^{2}.

Then we have the following estimate by applying $\tau\leq Ch$

\tau\|\nabla\tilde{\mu}_{u}^{k}\|_{L_{h}^{2}}\leq Ch^{15/4}.

(3.19)

Again, taking $\varphi=\tilde{u}^{k+1}-\tilde{u}^{k}$ in the error equation (3.9) leads to

\|\tilde{u}^{k+1}-\tilde{u}^{k}\|_{L_{h}^{2}}\leq\tau\left(\|\nabla\cdot(\tilde{u}^{k}\nabla\mathcal{V}_{u}^{k})\|_{L_{h}^{2}}+\|\nabla\cdot(u_{h}^{k}\nabla\tilde{\mu}_{u}^{k})\|_{L_{h}^{2}}+\dfrac{1}{h}\|\mathcal{R}^{k}\|_{H^{-1}}\right).

(3.20)

Now we estimate the first term on the right hand side of (3.20). Using the Young inequality, we have

$\displaystyle\tau\\|\nabla\cdot(\tilde{u}^{k}\nabla\mathcal{V}_{u}^{k})\\|_{L_{h}^{2}}$	$\displaystyle\leq\tau(\\|\tilde{u}^{k}\Delta\mathcal{V}_{u}^{k}\\|_{L_{h}^{2}}+\\|\nabla\tilde{u}^{k}\nabla\mathcal{V}_{u}^{k}\\|_{L_{h}^{2}})$	(3.21)
	$\displaystyle\leq\tau(\\|\Delta\mathcal{V}_{u}^{k}\\|_{\infty}\\|\tilde{u}^{k}\\|_{L_{h}^{2}}+\\|\nabla\mathcal{V}_{u}^{k}\\|_{\infty}\\|\nabla\tilde{u}^{k}\\|_{L_{h}^{2}})$
	$\displaystyle\leq C\tau(\\|\tilde{u}^{k}\\|_{L_{h}^{2}}+\\|\nabla\tilde{u}^{k}\\|_{L_{h}^{2}})$
	$\displaystyle\leq Ch^{15/4},$

where (3.11) and $\tau\leq Ch$ have been used in the last inequality. For the second term on the right hand side of (3.20), we have

$\displaystyle\tau\\|\nabla\cdot(u_{h}^{k}\nabla\tilde{\mu}_{u}^{k})\\|_{L_{h}^{2}}$	$\displaystyle\leq\tau(\\|\nabla u_{h}^{k}\nabla\tilde{\mu}_{u}^{k}\\|_{L_{h}^{2}}+\\|u_{h}^{k}\Delta\tilde{\mu}_{u}^{k}\\|_{L_{h}^{2}})$	(3.22)
	$\displaystyle\leq\tau(\\|\nabla u_{h}^{k}\\|_{\infty}\\|\nabla\tilde{\mu}_{u}^{k}\\|_{L_{h}^{2}}+\\|u_{h}^{k}\\|_{\infty}\\|\Delta\tilde{\mu}_{u}^{k}\\|_{L_{h}^{2}})$
	$\displaystyle\leq C^{*}\tau(\\|\nabla\tilde{\mu}_{u}^{k}\\|_{L_{h}^{2}}+\\|\Delta\tilde{\mu}_{u}^{k}\\|_{L_{h}^{2}})$
	$\displaystyle\leq Ch^{11/4},$

where (3.19) and the inverse inequality have been used in the last inequality. Substitution of (3.21)-(3.22) and (3.6) into (3.20) leads to

\|\tilde{u}^{k+1}-\tilde{u}^{k}\|_{L_{h}^{2}}\leq Ch^{11/4}.

Then, we can obtain a rough estimate for $\tilde{u}^{k+1}$ :

\|\tilde{u}^{k+1}\|_{L_{h}^{2}}\leq\|\tilde{u}^{k+1}-\tilde{u}^{k}\|_{L_{h}^{2}}+\|\tilde{u}^{k}\|_{L_{h}^{2}}\leq Ch^{11/4}.

An application of the inverse inequality implies that

	$\displaystyle\\|\tilde{u}^{k+1}\\|_{\infty}\leq\dfrac{C\\|\tilde{u}^{k+1}\\|_{L_{h}^{2}}}{h}\leq Ch^{7/4}\leq 1,$
	$\displaystyle\\|\nabla\tilde{u}^{k+1}\\|_{\infty}\leq\dfrac{C\\|\tilde{u}^{k+1}\\|_{\infty}}{h}\leq Ch^{3/4}\leq 1.$

We take $h$ sufficiently small such that

\|\tilde{u}^{k+1}\|_{\infty}\leq Ch^{7/4}\leq\dfrac{\epsilon_{0}}{4}.

A combination of above with (3.7) leads to the strict separation property:

u_{h}^{k+1}\geq I_{h}\hat{u}^{k+1}-\|\tilde{u}^{k+1}\|_{\infty}\geq\dfrac{\epsilon_{0}}{4}.

In addition, we can obtain the following $W^{1,\infty}$ bound for the numerical solution at time step $t^{k+1}$ :

	$\displaystyle\\|u_{h}^{k+1}\\|_{\infty}\leq\\|\hat{u}^{k+1}\\|_{\infty}+\\|\tilde{u}^{k+1}\\|_{\infty}\leq C+1\leq C^{*},$
	$\displaystyle\\|\nabla u_{h}^{k+1}\\|_{\infty}\leq\\|\nabla\hat{u}^{k+1}\\|_{\infty}+\\|\nabla\tilde{u}^{k+1}\\|_{\infty}\leq C+1\leq C^{*},$

which completes the proof. ∎

3.3. Recovery of the assumption (3.11)

In this subsection, the assumption will be recovered at the next time step.

Taking $\varphi=\tilde{u}^{k+1}$ in (3.9), we arrive at

\displaystyle\dfrac{1}{2\tau}\left(\|\tilde{u}^{k+1}\|_{L_{h}^{2}}^{2}-\|\tilde{u}^{k}\|_{L_{h}^{2}}^{2}\right)+(u_{h}^{k}\nabla\tilde{\mu}_{u}^{k},\nabla\tilde{u}^{k+1})_{h}\leq-(\tilde{u}^{k}\nabla\mathcal{V}_{u}^{k},\nabla\tilde{u}^{k+1})_{h}+\left\langle R^{k},\tilde{u}^{k+1}\right\rangle.

The second term on the left hand side of above inequality can be rewritten as

		$\displaystyle(u_{h}^{k}\nabla\tilde{\mu}_{u}^{k},\nabla\tilde{u}^{k+1})_{h}=(u_{h}^{k}(\nabla(I_{h}\log\hat{u}^{k+1}-I_{h}\log u_{h}^{k+1})-\chi\nabla A_{h}\tilde{u}^{k}),\nabla\tilde{u}^{k+1})_{h}$
	$\displaystyle=$	$\displaystyle\Big{(}u_{h}^{k}\nabla I_{h}\Big{(}\dfrac{\hat{u}^{k+1}-u_{h}^{k+1}}{\xi}\Big{)},\nabla\tilde{u}^{k+1}\Big{)}_{h}-(\chi u_{h}^{k}\nabla A_{h}\tilde{u}^{k},\nabla\tilde{u}^{k+1})_{h}$
	$\displaystyle=$	$\displaystyle\Big{(}\dfrac{u_{h}^{k}}{I_{h}\xi}\nabla\tilde{u}^{k+1},\nabla\tilde{u}^{k+1}\Big{)}_{h}-\Big{(}\dfrac{u_{h}^{k}\nabla I_{h}\xi}{(I_{h}\xi)^{2}}\tilde{u}^{k+1},\nabla\tilde{u}^{k+1}\Big{)}_{h}-(\chi u_{h}^{k}\nabla A_{h}\tilde{u}^{k},\nabla\tilde{u}^{k+1})_{h},$

where $\xi$ lies between $\hat{u}^{k+1}$ and $u_{h}^{k+1}$ . Utilizing the strict separation property and the $W^{1,\infty}$ bound for the numerical solution, we obtain the following estimate

		$\displaystyle(u_{h}^{k}\nabla\tilde{\mu}_{u}^{k},\nabla\tilde{u}^{k+1})_{h}$
	$\displaystyle\geq$	$\displaystyle\dfrac{\epsilon_{0}}{4C^{}}\\|\nabla\tilde{u}^{k+1}\\|_{2}^{2}-\dfrac{16(C^{})^{2}}{\epsilon_{0}^{2}}\|(\tilde{u}^{k+1},\nabla\tilde{u}^{k+1})_{h}\|-(\chi u_{h}^{k}\nabla A_{h}\tilde{u}^{k},\nabla\tilde{u}^{k+1})_{h}$
	$\displaystyle\geq$	$\displaystyle\dfrac{\epsilon_{0}}{4C^{}}\\|\nabla\tilde{u}^{k+1}\\|_{L_{h}^{2}}^{2}-\dfrac{\epsilon_{0}}{8C^{}}\\|\nabla\tilde{u}^{k+1}\\|_{L_{h}^{2}}^{2}-C\\|\tilde{u}^{k+1}\\|_{L_{h}^{2}}^{2}-(\chi u_{h}^{k}\nabla A_{h}\tilde{u}^{k},\nabla\tilde{u}^{k+1})_{h}$
	$\displaystyle=$	$\displaystyle\dfrac{\epsilon_{0}}{8C^{*}}\\|\nabla\tilde{u}^{k+1}\\|_{L_{h}^{2}}^{2}-C\\|\tilde{u}^{k+1}\\|_{L_{h}^{2}}^{2}-(\chi u_{h}^{k}\nabla A_{h}\tilde{u}^{k},\nabla\tilde{u}^{k+1})_{h},$

where $C$ is a positive constant satisfying $\sqrt{C\epsilon_{0}/8C^{*}}\geq 8(C^{*}/\epsilon_{0})^{2}$ . Combining above inequalities leads to

		$\displaystyle\dfrac{1}{2\tau}\left(\\|\tilde{u}^{k+1}\\|_{L_{h}^{2}}^{2}-\\|\tilde{u}^{k}\\|_{L_{h}^{2}}^{2}\right)+\dfrac{\epsilon_{0}}{8C^{*}}\\|\nabla\tilde{u}^{k+1}\\|_{L_{h}^{2}}^{2}$		(3.23)
	$\displaystyle\leq$	$\displaystyle\ C\\|\tilde{u}^{k+1}\\|_{L_{h}^{2}}^{2}+(\chi u_{h}^{k}\nabla A_{h}\tilde{u}^{k},\nabla\tilde{u}^{k+1})_{h}-(\tilde{u}^{k}\nabla\mathcal{V}_{u}^{k},\nabla\tilde{u}^{k+1})_{h}+\left\langle R^{k},\tilde{u}^{k+1}\right\rangle.$		(3.23)

Using Cauchy-Schwarz inequality and the strict separation property of the numerical solution, we have

		$\displaystyle\|(\chi u_{h}^{k}\nabla A_{h}\tilde{u}^{k},\nabla\tilde{u}^{k+1})_{h}\|\leq\chi\\|u_{h}^{k}\\|_{\infty}\\|\nabla A_{h}\tilde{u}^{k}\\|_{L_{h}^{2}}\\|\nabla\tilde{u}^{k+1}\\|_{L_{h}^{2}}$
	$\displaystyle\leq$	$\displaystyle C\\|\tilde{u}^{k}\\|_{L_{h}^{2}}\\|\nabla\tilde{u}^{k+1}\\|_{L_{h}^{2}}\leq\dfrac{\epsilon_{0}}{32C^{*}}\\|\nabla\tilde{u}^{k+1}\\|_{L_{h}^{2}}^{2}+C\\|\tilde{u}^{k}\\|_{L_{h}^{2}}^{2}.$

An application of (3.10) shows that

\displaystyle|-(\tilde{u}^{k}\nabla\mathcal{V}_{u}^{k},\nabla\tilde{u}^{k+1})_{h}|\leq C\|\tilde{u}^{k}\|_{L_{h}^{2}}\|\nabla\tilde{u}^{k+1}\|_{L_{h}^{2}}\leq\dfrac{\epsilon_{0}}{32C^{*}}\|\nabla\tilde{u}^{k+1}\|_{L_{h}^{2}}^{2}+C\|\tilde{u}^{k}\|_{L_{h}^{2}}^{2}.

The last term on the right hand side of (3.23) can be estimated as follows

\displaystyle\Big{|}\left\langle R^{k},\tilde{u}^{k+1}\right\rangle\Big{|}\leq Ch^{4}\|\nabla\tilde{u}^{k+1}\|_{2}\leq Ch^{8}+\dfrac{\epsilon_{0}}{32C^{*}}\|\nabla\tilde{u}^{k+1}\|_{L_{h}^{2}}^{2}.

Substitution of above into (3.23) leads to

\displaystyle\dfrac{1}{2\tau}\left(\|\tilde{u}^{k+1}\|_{L_{h}^{2}}^{2}-\|\tilde{u}^{k}\|_{L_{h}^{2}}^{2}\right)+\dfrac{\epsilon_{0}}{32C^{*}}\|\nabla\tilde{u}^{k+1}\|_{L_{h}^{2}}^{2}\leq Ch^{8}+C(\|\tilde{u}^{k}\|_{L_{h}^{2}}^{2}+\|\tilde{u}^{k+1}\|_{L_{h}^{2}}^{2}),

multiplying $2\tau$ on both sides and summing this inequality from 0 to $k$ leads to

\displaystyle\|\tilde{u}^{k+1}\|_{L_{h}^{2}}^{2}+\sum_{i=0}^{k}\dfrac{\epsilon_{0}\tau}{16C^{*}}\|\nabla\tilde{u}^{i+1}\|_{L_{h}^{2}}^{2}\leq Ch^{8}+\sum_{i=0}^{k+1}4C\tau\|\tilde{u}^{i}\|_{L_{h}^{2}}^{2}.

Choosing $\tau$ sufficiently small such that $1-4C\tau>1/2$ , using the Gronwall inequality (Lemma 11), we derive that

\displaystyle\|\tilde{u}^{k+1}\|_{L_{h}^{2}}^{2}+\sum_{i=0}^{k}\dfrac{\epsilon_{0}\tau}{8C^{*}}\|\nabla\tilde{u}^{i+1}\|_{L_{h}^{2}}^{2}\leq C\exp\Big{(}\sum_{i=0}^{k}8C\tau\Big{)}h^{8},

then we obtain $\|\tilde{u}^{k+1}\|_{L_{h}^{2}}\leq Ch^{4}\leq Ch^{15/4}$ , the assumption (3.11) is recovered at the next time step.

3.4. Proof of Theorem 6

In this subsection, we shall make use of the strict separation property and the uniform $W^{1,\infty}$ bound for the numerical solution derived in the above to prove Theorem 6.

A weak formulation of (1.1)-(1.2) is

	$\displaystyle(u_{t},\varphi)+(\nabla u,\nabla\varphi)-\chi(u\nabla v,\nabla\varphi)=0,\ \forall\varphi\in H_{0}^{1}(\Omega),$		(3.24)
	$\displaystyle(\nabla v,\nabla\psi)+(v,\psi)=\alpha(u,\psi),\ \forall\psi\in H_{0}^{1}(\Omega).$		(3.25)

Substituting $I_{h}u(t)$ into (3.24) at $t=t_{k+1}$ , we have

	$\displaystyle\left(\overline{\partial}I_{h}u^{k},\ \varphi\right)_{h}+\left(\nabla I_{h}u^{k+1},\ \nabla\varphi\right)_{h}-\chi\left(I_{h}u^{k}\nabla v^{k},\ \nabla\varphi\right)_{h}$	(3.26)
$\displaystyle=$	$\displaystyle\left(\overline{\partial}u^{k}-u_{t}^{k+1},\ \varphi\right)+\left(\overline{\partial}I_{h}u^{k}-\overline{\partial}u^{k},\ \varphi\right)+\left(\nabla(I_{h}u^{k+1}-u^{k+1}),\ \varphi\right)$
	$\displaystyle+\chi\left((u^{k}-I_{h}u^{k})\nabla v^{k},\ \varphi\right)+\chi\left(u^{k+1}\nabla v^{k+1}-u^{k}\nabla v^{k},\ \varphi\right)$
	$\displaystyle+\left(\overline{\partial}I_{h}u^{k},\varphi\right)_{h}-\left(\overline{\partial}I_{h}u^{k},\varphi\right)+\left(\nabla I_{h}u^{k+1},\nabla\varphi\right)_{h}-\left(\nabla I_{h}u^{k+1},\nabla\varphi\right)$
	$\displaystyle+\chi\left(I_{h}u^{k}\nabla v^{k},\nabla\varphi\right)-\chi\left(I_{h}u^{k}\nabla v^{k},\nabla\varphi\right)_{h}$
$\displaystyle:=$	$\displaystyle I_{1}+I_{2}+I_{3}+I_{4}+I_{5}+I_{6}+I_{7}+I_{8}:=R_{1},$

where $R_{1}$ represents the truncation error. Similarly, substituting $I_{h}u(t),I_{h}v(t)$ into (3.25) at $t=t_{k}$ leads to

	$\displaystyle(\nabla I_{h}v^{k},\ \nabla\psi)+(\I_{h}v^{k},\ \psi)=$	$\displaystyle\alpha(I_{h}u^{k},\ \psi)+(\nabla(I_{h}v(t_{k})-v^{k}),\ \psi)$		(3.27)
		$\displaystyle+(I_{h}v^{k}-v^{k},\ \psi)+\alpha(u^{k}-I_{h}u^{k},\ \psi).$

We rewrite the numerical scheme (3.3) as

	$\displaystyle\left(\overline{\partial}u_{h}^{k},\varphi\right)_{h}+\left(\nabla u_{h}^{k+1},\nabla\varphi\right)_{h}-\chi\left(u_{h}^{k}\nabla v_{h}^{k},\nabla\varphi\right)_{h}$	(3.28)
$\displaystyle=$	$\displaystyle\left(u_{h}^{k}\nabla(I-I_{h})\log u_{h}^{k+1},\nabla\varphi\right)_{h}+\left((u_{h}^{k+1}-u_{h}^{k})\nabla\log u_{h}^{k+1},\nabla\varphi\right)_{h}$
$\displaystyle:=$	$\displaystyle I_{9}+I_{10}:=R_{2}.$

We split the error functions as

	$\displaystyle e_{u}^{k}$	$\displaystyle=u^{k}-u_{h}^{k}=(u^{k}-I_{h}u^{k})+(I_{h}u^{k}-u_{h}^{k}):=\rho_{u}^{k}+\sigma_{u}^{k},$
	$\displaystyle e_{v}^{k}$	$\displaystyle=v^{k}-v_{h}^{k}=(v^{k}-I_{h}v^{k})+(I_{h}v^{k}-v_{h}^{k}):=\rho_{v}^{k}+\sigma_{v}^{k}.$

Then using the property of the interpolation (3.2), we have

	$\displaystyle\\|\rho_{u}^{k}\\|_{L^{2}(\Omega)}+h\\|\nabla\rho_{u}^{k}\\|_{L^{2}(\Omega)}\leq Ch^{2}\\|u\\|_{H^{2}(\Omega)},$		(3.29)
	$\displaystyle\\|\rho_{v}^{k}\\|_{L^{2}(\Omega)}+h\\|\nabla\rho_{v}^{k}\\|_{L^{2}(\Omega)}\leq Ch^{2}\\|v\\|_{H^{2}(\Omega)}.$		(3.29)

Subtracting the numerical scheme formulation (3.28) and (3.4) from the weak form (3.26) and (3.27), we obtain the following error equations:

$\displaystyle\left(\overline{\partial}\sigma_{u}^{k},\ \varphi\right)_{h}+\left(\nabla\sigma_{u}^{k+1},\ \nabla\varphi\right)$	${}_{h}=I_{11}+R_{1}-R_{2},\$	(3.30)
$\displaystyle(\nabla\sigma_{v}^{k},\ \nabla\psi)+(\sigma_{v}^{k},\ \psi)=\alpha$	$\displaystyle(\sigma_{u}^{k},\ \psi)-\alpha\epsilon_{h}(u_{h}^{k},\ \psi)+(\nabla(I_{h}v-v),\ \psi)\$	(3.31)
	$\displaystyle+(I_{h}v-v,\ \psi)+\alpha(u-I_{h}u,\ \psi),$

for all $\varphi,\psi\in X_{h},k\geq 1$ , where $R_{1},R_{2}$ are defined before and $I_{11}$ is defined as follows

\displaystyle I_{11}=\chi\left(I_{h}u(t_{k})\nabla v^{k}-u_{h}^{k}\nabla v_{h}^{k},\ \nabla\varphi\right)_{h}.

Taking $\varphi=\sigma_{u}^{k+1}$ in (3.30) leads to

\displaystyle\dfrac{1}{2\tau}\left(\|\sigma_{u}^{k+1}\|_{L_{h}^{2}}^{2}-\|\sigma_{u}^{k}\|_{L_{h}^{2}}^{2}+\|\sigma_{u}^{k+1}-\sigma_{u}^{k}\|_{L_{h}^{2}}^{2}\right)+\|\nabla\sigma_{u}^{k+1}\|_{L_{h}^{2}}^{2}=I_{11}+R_{1}-R_{2}.

(3.32)

Now we estimate the terms on the right-hand side of (3.33). For the first term $I_{11}$ , applying the Cauchy-Schwarz inequality and the Young inequality yields

	$\displaystyle\|I_{11}\|=\|\chi\left(\sigma_{u}^{k}\nabla v^{k}+u_{h}^{k}\nabla e_{v}^{k},\ \nabla\sigma_{u}^{k+1}\right)_{h}\|$	(3.33)
$\displaystyle\leq$	$\displaystyle\ \\|\sigma_{u}^{k}\nabla v^{k}\\|_{L_{h}^{2}(\Omega)}\\|\nabla\sigma_{u}^{k+1}\\|_{L_{h}^{2}(\Omega)}+\\|u_{h}^{k}\nabla e_{v}^{k}\\|_{L_{h}^{2}(\Omega)}\\|\nabla\sigma_{u}^{k+1}\\|_{L_{h}^{2}(\Omega)}$
$\displaystyle\leq$	$\displaystyle\ C\\|\sigma_{u}^{k}\\|_{L_{h}^{2}(\Omega)}^{2}+\dfrac{1}{20}\\|\nabla\sigma_{u}^{k+1}\\|_{L_{h}^{2}(\Omega)}^{2}+C\\|\nabla e_{v}^{k}\\|_{L_{h}^{2}(\Omega)}^{2}.$

In order to estimate $\|\nabla e_{v}^{k}\|_{L_{h}^{2}(\Omega)}^{2}$ above, taking $\psi=\sigma_{v}^{k}$ in (3.31) and applying Lemma 4 leads to

	$\displaystyle\\|\nabla\sigma_{v}^{k}\\|_{L^{2}(\Omega)}^{2}+\\|\sigma_{v}^{k}\\|_{L^{2}(\Omega)}^{2}\leq$	$\displaystyle\ \alpha\\|\sigma_{u}^{k}\\|_{L^{2}(\Omega)}\\|\sigma_{v}^{k}\\|_{L^{2}(\Omega)}+Ch^{2}\\|\nabla u_{h}^{k}\\|_{L^{2}(\Omega)}\\|\nabla\sigma_{v}^{k}\\|_{L^{2}(\Omega)}$
		$\displaystyle+Ch\\|v\\|_{H^{2}(\Omega)}\\|\sigma_{v}^{k}\\|_{L^{2}(\Omega)}+Ch^{2}\\|u\\|_{H^{2}(\Omega)}\\|\sigma_{v}^{k}\\|_{L^{2}(\Omega)}$
	$\displaystyle\leq$	$\displaystyle\ C\\|\sigma_{u}^{k}\\|_{L^{2}(\Omega)}^{2}+\\|\sigma_{v}^{k}\\|_{L^{2}(\Omega)}^{2}+\dfrac{1}{2}\\|\nabla\sigma_{v}^{k}\\|_{L^{2}(\Omega)}^{2}+Ch^{2},$

where the property of the interpolation has been used in the first inequality. Thus we obtain the following estimate for $\|\nabla\sigma_{v}^{k}\|_{L^{2}(\Omega)}$ :

\displaystyle\|\nabla\sigma_{v}^{k}\|_{L^{2}(\Omega)}^{2}\leq C\|\sigma_{u}^{k}\|_{L^{2}(\Omega)}^{2}+Ch^{2}.

(3.34)

Applying Lemma 4 indicates that

\displaystyle\|\nabla\sigma_{v}^{k}\|_{L_{h}^{2}(\Omega)}^{2}\leq C\|\sigma_{u}^{k}\|_{L_{h}^{2}(\Omega)}^{2}+Ch^{2}\|\nabla\sigma_{u}^{k}\|_{L^{2}(\Omega)}^{2}+Ch^{2}.

A combination of the above estimates for $\|\nabla\sigma_{v}^{k}\|_{L_{h}^{2}(\Omega)}$ with (3.29) leads to

\displaystyle\|\nabla e_{v}^{k}\|_{L_{h}^{2}(\Omega)}^{2}\leq C\|\sigma_{u}^{k}\|_{L_{h}^{2}(\Omega)}^{2}+Ch^{2}\|\nabla\sigma_{u}^{k}\|_{L^{2}(\Omega)}^{2}+Ch^{2}.

(3.35)

Substitution of above into (3.33) leads to

\displaystyle|I_{11}|

\displaystyle\leq C\|\sigma_{u}^{k}\|_{L_{h}^{2}(\Omega)}^{2}+\dfrac{1}{20}\|\nabla\sigma_{u}^{k+1}\|_{L_{h}^{2}(\Omega)}^{2}+Ch^{2}\|\nabla\sigma_{u}^{k}\|_{L^{2}(\Omega)}^{2}+Ch^{2}.

(3.36)

Next we estimate the second term $R_{1}=I_{1}+I_{2}+I_{3}+I_{4}+I_{5}+I_{6}+I_{7}+I_{8}$ . For the first term, we can derive [48]

	$\displaystyle\|I_{1}\|$	$\displaystyle=\|\left(\overline{\partial}u^{k}-u_{t}^{k+1},\ \sigma_{u}^{k+1}\right)\|\leq C\tau^{1/2}\Big{(}\int_{t_{k}}^{t_{k+1}}\\|u_{tt}\\|_{L^{2}(\Omega)}^{2}\mathrm{d}s\Big{)}^{1/2}\\|\sigma_{u}^{k+1}\\|_{L^{2}(\Omega)}$		(3.37)
		$\displaystyle\leq C\tau\int_{t_{k}}^{t_{k+1}}\\|u_{tt}\\|_{L^{2}(\Omega)}^{2}\mathrm{d}s+\dfrac{1}{20}\\|\nabla\sigma_{u}^{k+1}\\|_{L^{2}(\Omega)}^{2}.$		(3.37)

An application of the property of the interpolation and the Young inequality leads to

$\displaystyle\|I_{2}\|=$	$\displaystyle\|\left(\overline{\partial}I_{h}u^{k}-\overline{\partial}u^{k},\ \sigma_{u}^{k+1}\right)\|$	(3.38)
$\displaystyle\leq$	$\displaystyle\ \dfrac{1}{\tau}\\|I_{h}(u^{k+1}-u^{k})-(u^{k+1}-u^{k})\\|_{L^{2}(\Omega)}\\|\sigma_{u}^{k+1}\\|_{L^{2}(\Omega)}$
$\displaystyle=$	$\displaystyle\ \dfrac{1}{\tau}\\|(I_{h}u_{t}-u_{t})\tau\\|_{L^{2}(\Omega)}\\|\sigma_{u}^{k+1}\\|_{L^{2}(\Omega)}$
$\displaystyle\leq$	$\displaystyle\ Ch^{2}\\|u_{t}\\|_{H^{2}(\Omega)}\\|\nabla\sigma_{u}^{k+1}\\|_{L^{2}(\Omega)}$
$\displaystyle\leq$	$\displaystyle\ Ch^{4}+\dfrac{1}{20}\\|\nabla\sigma_{u}^{k+1}\\|_{L^{2}(\Omega)}^{2}.$

Using the Cauchy-Schwarz inequality and the property of the interpolation, we have

$\displaystyle\|I_{3}\|=$	$\displaystyle\|\left(\nabla(I_{h}u^{k+1}-u^{k+1}),\ \nabla\sigma_{u}^{k+1}\right)\|$	(3.39)
$\displaystyle\leq$	$\displaystyle\\|\nabla(I_{h}u^{k+1}-u^{k+1})\\|_{L^{2}(\Omega)}\\|\nabla\sigma_{u}^{k+1}\\|_{L^{2}(\Omega)}$
$\displaystyle\leq$	$\displaystyle Ch\\|u^{k+1}\\|_{H^{2}(\Omega)}\\|\nabla\sigma_{u}^{k+1}\\|_{L^{2}(\Omega)}$
$\displaystyle\leq$	$\displaystyle Ch^{2}+\dfrac{1}{20}\\|\nabla\sigma_{u}^{k+1}\\|_{L^{2}(\Omega)}^{2}.$

Similarly, we have

$\displaystyle\|I_{4}\|=$	$\displaystyle\|\chi\left((u^{k}-I_{h}u^{k})\nabla v^{k},\ \nabla\sigma_{u}^{k+1}\right)\|$	(3.40)
$\displaystyle\leq$	$\displaystyle\ \chi\\|u^{k}-I_{h}u^{k}\\|_{L^{2}(\Omega)}\\|\nabla v^{k}\\|_{L^{\infty}(\Omega)}\\|\nabla\sigma_{u}^{k+1}\\|_{L^{2}(\Omega)}$
$\displaystyle\leq$	$\displaystyle\ Ch^{2}\\|u^{k}\\|_{H^{2}(\Omega)}\\|\nabla\sigma_{u}^{k+1}\\|_{L^{2}(\Omega)}$
$\displaystyle\leq$	$\displaystyle\ Ch^{4}+\dfrac{1}{20}\\|\nabla\sigma_{u}^{k+1}\\|_{L^{2}(\Omega)}^{2}.$

Now we apply the Cauchy-Schwarz inequality and the Young inequality:

$\displaystyle\|I_{5}\|=$	$\displaystyle\|\chi\left(\tau u_{t}\nabla v^{k+1}+\tau u^{k}\nabla v_{t},\ \nabla\sigma_{u}^{k+1}\right)\|$	(3.41)
$\displaystyle\leq$	$\displaystyle\ C\left(\\|u_{t}\\|_{L^{\infty}(\Omega)}\\|\nabla v\\|_{L^{\infty}(\Omega)}+\\|u\\|_{L^{\infty}(\Omega)}\\|\nabla v_{t}\\|_{L^{\infty}(\Omega)}\right)\tau\\|\nabla\sigma_{u}^{k+1}\\|_{L^{2}(\Omega)}$
$\displaystyle\leq$	$\displaystyle\ C\tau^{2}+\dfrac{1}{20}\\|\nabla\sigma_{u}^{k+1}\\|_{L^{2}(\Omega)}^{2}.$

Notice that

	$\displaystyle I_{6}$	$\displaystyle=\dfrac{1}{\tau}\left(I_{h}u^{k+1}-I_{h}u^{k},\ \sigma_{u}^{k+1}\right)-\dfrac{1}{\tau}\left(I_{h}u^{k+1}-I_{h}u^{k},\ \sigma_{u}^{k+1}\right)_{h}$
		$\displaystyle=-\epsilon_{h}\left(I_{h}u_{t},\ \sigma_{u}^{k+1}\right),$

therefore, using Lemma 4 leads to the following estimate for $I_{6}$ :

\displaystyle|I_{6}|=|\epsilon_{h}(I_{h}u_{t},\sigma_{u}^{k+1})|\leq Ch^{2}\|\nabla I_{h}u_{t}\|\|\nabla\sigma_{u}^{k+1}\|\leq Ch^{4}+\dfrac{1}{20}\|\nabla\sigma_{u}^{k+1}\|_{L^{2}(\Omega)}^{2}.

(3.42)

We recall the quadrature formula defined before and Lemma 4, and arrive at

\displaystyle I_{7}

\displaystyle=\left(I_{h}u^{k+1},\nabla\sigma_{u}^{k+1}\right)_{h}-\left(I_{h}u^{k+1},\nabla\sigma_{u}^{k+1}\right)=0,

(3.43)

An application of the Cauchy-Schwarz inequality and the property of the interpolation leads to

$\displaystyle\|I_{8}\|$	$\displaystyle=\|\chi\left(I_{h}u^{k}\nabla v^{k},\nabla\sigma_{u}^{k+1}\right)-\chi\left(I_{h}u^{k}\nabla v^{k},\nabla\sigma_{u}^{k+1}\right)_{h}\|$	(3.44)
	$\displaystyle=\|\chi\left(I_{h}u^{k}\nabla(I-I_{h})v^{k},\nabla\sigma_{u}^{k+1}\right)-\chi\left(I_{h}u^{k}\nabla(I-I_{h})v^{k},\nabla\sigma_{u}^{k+1}\right)_{h}\|$
	$\displaystyle\leq Ch\\|I_{h}u^{k}\\|_{L^{\infty}}\\|u^{k}\\|_{L^{2}(\Omega)}\\|\nabla\sigma_{u}^{k+1}\\|_{L^{2}(\Omega)}$
	$\displaystyle\leq Ch^{2}+\dfrac{1}{20}\\|\nabla\sigma_{u}^{k+1}\\|_{L^{2}(\Omega)}^{2}.$

Substituting estimates (3.37)-(3.44) into $R_{1}$ and applying the property $\eqref{norm}$ , we obtain

	$\displaystyle\|R_{1}\|\leq$	$\displaystyle C\tau\int_{t_{k}}^{t_{k+1}}\\|u_{tt}\\|_{L^{2}(\Omega)}^{2}\mathrm{d}s+\dfrac{7}{20}\\|\nabla\sigma_{u}^{k+1}\\|_{L^{2}(\Omega)}^{2}+Ch^{2}+C\tau^{2}$		(3.45)
	$\displaystyle=$	$\displaystyle C\tau\int_{t_{k}}^{t_{k+1}}\\|u_{tt}\\|_{L^{2}(\Omega)}^{2}\mathrm{d}s+\dfrac{7}{20}\\|\nabla\sigma_{u}^{k+1}\\|_{L_{h}^{2}(\Omega)}^{2}+Ch^{2}+C\tau^{2}.$		(3.45)

It remains to bound each term of $R_{2}=I_{9}+I_{10}$ . Now we use the Cauchy-Schwarz inequality and the Young inequality:

$\displaystyle\|I_{9}\|=$	$\displaystyle\|\left(u_{h}^{k}\nabla(I-I_{h})\log u_{h}^{k+1},\nabla\sigma_{u}^{k+1}\right)_{h}\|$	(3.46)
$\displaystyle\leq$	$\displaystyle\ \\|u_{h}^{k}\nabla(I-I_{h})\log u_{h}^{k+1}\\|_{L_{h}^{2}(\Omega)}\\|\nabla\sigma_{u}^{k+1}\\|_{L_{h}^{2}(\Omega)}$
$\displaystyle\leq$	$\displaystyle\ \\|u_{h}^{k}\\|_{L^{\infty}(\Omega)}\\|\nabla(I-I_{h})\log u_{h}^{k+1}\\|_{L_{h}^{2}(\Omega)}\\|\nabla\sigma_{u}^{k+1}\\|_{L_{h}^{2}(\Omega)}$
$\displaystyle\leq$	$\displaystyle\ Ch^{2}+\dfrac{1}{20}\\|\nabla\sigma_{u}^{k+1}\\|_{L_{h}^{2}(\Omega)},$

where we have used the following inequality:

	$\displaystyle\\|\nabla(I-I_{h})\log u_{h}^{k+1}\\|_{L_{h}^{2}(\Omega)}\leq$	$\displaystyle\ Ch\sum_{e}\sum_{\|\alpha\|=2}\Big{(}\int_{e}\|D^{\alpha}\log u_{h}^{k+1}\|^{2}\mathrm{d}x\Big{)}^{1/2}$
	$\displaystyle\leq$	$\displaystyle\ Ch\sum_{e}\Big{(}\int_{e}\dfrac{1}{(u_{h}^{k+1})^{4}}\|\nabla u_{h}^{k+1}\|^{4}\mathrm{d}x\Big{)}^{1/2}$
	$\displaystyle\leq$	$\displaystyle\ Ch\dfrac{1}{(\epsilon_{0}/4)^{4}}\sum_{e}\Big{(}\int_{e}\|\nabla u_{h}^{k+1}\|^{4}\mathrm{d}x\Big{)}^{1/2}$
	$\displaystyle\leq$	$\displaystyle Ch.$

An application of the strict separation property and the $W^{1,\infty}$ bound of the numerical solution leads to

$\displaystyle\|I_{10}\|=$	$\displaystyle\|\left((u_{h}^{k+1}-u_{h}^{k})\nabla\log u_{h}^{k+1},\nabla\sigma_{u}^{k+1}\right)_{h}\|$	(3.47)
$\displaystyle\leq$	$\displaystyle\ \\|\dfrac{u_{h}^{k+1}-u_{h}^{k}}{u_{h}^{k+1}}\nabla u_{h}^{k+1}\\|_{L_{h}^{2}(\Omega)}\\|\nabla\sigma_{u}^{k+1}\\|_{L_{h}^{2}(\Omega)}$
$\displaystyle\leq$	$\displaystyle\ \dfrac{\\|\nabla u_{h}^{k+1}\\|_{L^{\infty}(\Omega)}}{\epsilon_{0}/4}\left(\\|\sigma_{u}^{k+1}-\sigma_{u}^{k}\\|_{L_{h}^{2}(\Omega)}+C\tau\right)\\|\nabla\sigma_{u}^{k+1}\\|_{L_{h}^{2}(\Omega)}$
$\displaystyle\leq$	$\displaystyle\ C(\\|\sigma_{u}^{k+1}-\sigma_{u}^{k}\\|_{L_{h}^{2}(\Omega)}+C\tau)\\|\nabla\sigma_{u}^{k+1}\\|_{L_{h}^{2}(\Omega)}$
$\displaystyle\leq$	$\displaystyle\ C\\|\sigma_{u}^{k+1}-\sigma_{u}^{k}\\|_{L_{h}^{2}(\Omega)}^{2}+C\tau^{2}+\dfrac{1}{20}\\|\nabla\sigma_{u}^{k+1}\\|_{L_{h}^{2}(\Omega)}^{2}.$

Combining the estimates (3.46)-(3.47), we obtain

\displaystyle|R_{2}|\leq C\|\sigma_{u}^{k+1}-\sigma_{u}^{k}\|_{L_{h}^{2}(\Omega)}^{2}+\dfrac{1}{10}\|\nabla\sigma_{u}^{k+1}\|_{L_{h}^{2}(\Omega)}^{2}+Ch^{2}+C\tau^{2}.

(3.48)

Finally, combining (3.36),(3.45) and (3.48) in (3.32), we find

$\displaystyle\dfrac{1}{2\tau}$	$\displaystyle\left(\\|\sigma_{u}^{k+1}\\|_{L_{h}^{2}(\Omega)}^{2}-\\|\sigma_{u}^{k}\\|_{L_{h}^{2}(\Omega)}^{2}+\\|\sigma_{u}^{k+1}-\sigma_{u}^{k}\\|_{L_{h}^{2}(\Omega)}^{2}\right)+\\|\nabla\sigma_{u}^{k+1}\\|_{L_{h}^{2}(\Omega)}^{2}$	(3.49)
$\displaystyle\leq$	$\displaystyle\ C\\|\sigma_{u}^{k}\\|_{L^{2}(\Omega)}^{2}+C\\|\sigma_{u}^{k+1}-\sigma_{u}^{k}\\|_{L_{h}^{2}(\Omega)}^{2}+\dfrac{1}{2}\\|\nabla\sigma_{u}^{k+1}\\|_{L_{h}^{2}(\Omega)}^{2}+Ch^{2}+C\tau^{2}$
	$\displaystyle+Ch^{2}\\|\nabla\sigma_{u}^{k}\\|_{L^{2}(\Omega)}^{2}+C\tau\int_{t_{k}}^{t_{k+1}}\\|u_{tt}\\|_{L^{2}(\Omega)}^{2}\mathrm{d}s.$

Multiplying by $2\tau$ on both sides of (3.49) and summing up from $k=0$ to $m-1$ , we get

	$\displaystyle\\|\sigma_{u}^{m}\\|_{L_{h}^{2}(\Omega)}^{2}-\\|\sigma_{u}^{0}\\|_{L_{h}^{2}(\Omega)}^{2}+(1-2C\tau)\\|\sigma_{u}^{k+1}-\sigma_{u}^{k}\\|_{L_{h}^{2}(\Omega)}^{2}+\tau\sum_{k=0}^{m-1}\\|\nabla\sigma_{u}^{k+1}\\|_{L_{h}^{2}(\Omega)}^{2}$
	$\displaystyle\leq C\tau\sum_{k=0}^{m-1}\\|\sigma_{u}^{k}\\|_{L_{h}^{2}(\Omega)}^{2}+Ch^{2}+C\tau^{2}+C\tau h^{2}\sum_{k=0}^{m-1}\\|\nabla\sigma_{u}^{k}\\|_{L^{2}(\Omega)}^{2}+C\tau^{2}\int_{0}^{T}\\|u_{tt}\\|_{L^{2}}^{2}\mathrm{d}s.$

Assuming $1-2C\tau>0$ , $1-Ch^{2}>1/2$ since $\tau$ and $h$ is small enough, and applying the discrete Gronwall inequality (Lemma 11) to the above leads to

\|\sigma_{u}^{m}\|_{L_{h}^{2}(\Omega)}+\Big{(}\tau\sum_{k=0}^{m-1}\|\nabla\sigma_{u}^{k+1}\|_{L_{h}^{2}(\Omega)}^{2}\Big{)}^{1/2}\leq C(h+\tau),\ \forall m\in\mathbb{N}.

A combination of the above estimates for $\|\sigma_{u}^{m}\|_{L_{h}^{2}(\Omega)}$ and $\|\nabla\sigma_{u}^{k+1}\|_{L_{h}^{2}(\Omega)}$ with (3.29) leads to the desired error estimates for $u$ . Finally, we obtain the error estimates for $v$ from (3.34) and (3.35).

The proof of Theorem 6 is complete.

4. Finite-time blowup

In this section, we discuss whether the solution of the fully discrete scheme (3.3)-(3.4) will blow up in finite time.

We first prove a discrete analog of Lemma 2. Taking $\varphi=I_{h}\Phi$ in (3.3), where $\Phi$ is defined as in Lemma 2, then from (3.2), we have the following error estimate

\|\varphi-\Phi\|_{L^{2}(\Omega)}+h\|\varphi-\Phi\|_{L^{2}(\Omega)}\leq Ch^{2}\|\Phi\|_{H^{2}(\Omega)}.

(4.1)

Lemma 13.

Assume that $u$ is smooth enough for a fix time $T>0$ , let $M^{k}=(u_{h}^{k},\ \varphi)_{h},\ \theta=(u,1)$ . Under the mild mesh-sizes requirement $\tau\leq Ch$ , if $(u_{h}^{0},\varphi)_{h}<\infty$ , then $(u_{h}^{k},\varphi)_{h}<\infty$ and the following inequality holds

\frac{M^{k+1}-M^{k}}{\tau}\leq 4\theta-\dfrac{\alpha\chi}{2\pi}\theta^{2}+C_{1}\theta M^{k}+C_{2}\theta^{3/2}(M^{k}+Ch)^{1/2}+C_{0}h.

Proof.

We rewrite (3.3) as

	$\displaystyle\frac{M^{k+1}-M^{k}}{\tau}$	$\displaystyle=-(\nabla u^{k+1},\nabla\Phi)+\chi(u^{k}\nabla v^{k},\nabla\Phi)+\sum_{i=1}^{8}J_{i}$		(4.2)
		$\displaystyle\leq 4\theta-\dfrac{\alpha\chi}{2\pi}\theta^{2}+C_{1}\theta(u^{k},\Phi)+C_{2}\theta^{3/2}(u^{k},\Phi)^{1/2}+\sum_{i=1}^{8}J_{i},$		(4.2)

where Lemma 2 has been used in the last inequality, and $J_{i}\ (i=1,\cdots,8)$ are defined as follows

	$\displaystyle J_{1}=$	$\displaystyle(\nabla u^{k+1},\ \nabla(\Phi-\varphi)),$
	$\displaystyle J_{2}=$	$\displaystyle(\nabla(u^{k+1}-u_{h}^{k+1}),\ \nabla\varphi),$
	$\displaystyle J_{3}=$	$\displaystyle\left(\nabla u_{h}^{k+1},\nabla\varphi\right)-\left(\nabla u_{h}^{k+1},\nabla\varphi\right)_{h},$
	$\displaystyle J_{4}=$	$\displaystyle\left((u_{h}^{k+1}-u_{h}^{k})\nabla\log u_{h}^{k+1},\nabla\varphi\right)_{h},$
	$\displaystyle J_{5}=$	$\displaystyle\chi\left(u^{k}\nabla v^{k},\nabla(\varphi-\Phi)\right),$
	$\displaystyle J_{6}=$	$\displaystyle\chi(u_{h}^{k}\nabla v_{h}^{k}-u^{k}\nabla v^{k},\nabla\varphi),$
	$\displaystyle J_{7}=$	$\displaystyle\chi(u_{h}^{k}\nabla v_{h}^{k},\nabla\varphi)_{h}-\chi(u_{h}^{k}\nabla v_{h}^{k},\nabla\varphi),$
	$\displaystyle J_{8}=$	$\displaystyle\left(u_{h}^{k}\nabla((I-I_{h})\log u_{h}^{k+1}),\nabla\varphi\right)_{h}.$

Thanks to (4.1), we obtain

		$\displaystyle(u^{k},\Phi)=(u^{k},\Phi-\varphi)+(u^{k}-u_{h}^{k},\varphi)+(u_{h}^{k},\varphi)-(u_{h}^{k},\varphi)_{h}+(u_{h}^{k},\varphi)_{h}$
	$\displaystyle\leq$	$\displaystyle\ \\|u^{k}\\|_{L^{2}(\Omega)}\\|\Phi-\varphi\\|_{L^{2}(\Omega)}+\\|e_{u}^{k}\\|_{L^{2}(\Omega)}\\|\varphi\\|_{L^{2}(\Omega)}+Ch^{2}\\|\nabla u_{h}^{k}\\|_{L^{2}(\Omega)}\\|\nabla\varphi\\|_{L^{2}(\Omega)}+M^{k}$
	$\displaystyle\leq$	$\displaystyle\ Ch+M^{k},$

where the error estimates and $\tau\leq Ch$ have been used in the last inequality. Substitution of above into (4.2) leads to

\frac{M^{k+1}-M^{k}}{\tau}\leq 4\theta-\dfrac{\alpha\chi}{2\pi}\theta^{2}+C_{1}\theta M^{k}+C_{2}\theta^{3/2}(M^{k}+Ch)^{1/2}+\sum_{i=1}^{8}J_{i}.

(4.3)

Next, we estimate $J_{i}\ (i=1,\cdots,8)$ respectively. Utilizing the property of the interpolation operator in (4.1), we have

|J_{1}|\leq\|\nabla u^{k+1}\|_{L^{2}(\Omega)}\|\nabla(\Phi-\varphi)\|_{L^{2}(\Omega)}\leq Ch.

(4.4)

We derive from Theorem 6 that

$\displaystyle\|J_{2}\|$	$\displaystyle=\|(\nabla(u^{k+1}-u_{h}^{k+1}),\ \nabla(\varphi-\Phi))+(\nabla(u^{k+1}-u_{h}^{k+1}),\ \nabla\Phi)\|$	(4.5)
	$\displaystyle=\|(\nabla(u^{k+1}-u_{h}^{k+1}),\ \nabla(\varphi-\Phi))-((u^{k+1}-u_{h}^{k+1}),\Delta\Phi)\|$
	$\displaystyle\leq Ch\\|\nabla e_{u}^{k}\\|_{L^{2}(\Omega)}\\|\Phi\\|_{H^{2}(\Omega)}+\\|e_{u}^{k}\\|_{L^{2}(\Omega)}\\|\Delta\Phi\\|_{L^{2}(\Omega)}$
	$\displaystyle\leq Ch,$

where the property of interpolation and $\tau\leq Ch$ have been used. Noticing the definition of $(\cdot,\cdot)_{h}$ , we have

J_{3}=\left(\nabla u_{h}^{k+1},\nabla\varphi\right)-\left(\nabla u_{h}^{k+1},\nabla\varphi\right)_{h}=0.

(4.6)

An application of the strict separation property, the $W^{1,\infty}$ bound for the numerical solution and $\tau\leq Ch$ indicates that

$\displaystyle\|J_{4}\|$	$\displaystyle=\Big{\|}\big{(}\dfrac{u_{h}^{k+1}-u_{h}^{k}}{u_{h}^{k+1}}\nabla u_{h}^{k+1},\nabla\varphi\Big{)}_{h}\Big{\|}$	(4.7)
	$\displaystyle\leq\Big{\\|}\dfrac{u_{h}^{k+1}-u_{h}^{k}}{u_{h}^{k+1}}\nabla u_{h}^{k+1}\Big{\\|}_{L_{h}^{2}(\Omega)}\\|\nabla\varphi\\|_{L_{h}^{2}(\Omega)}$
	$\displaystyle\leq\dfrac{\\|\nabla u_{h}^{k+1}\\|_{L^{\infty}(\Omega)}}{\epsilon_{0}/4}\\|u_{h}^{k+1}-u_{h}^{k}\\|_{L_{h}^{2}(\Omega)}\\|\nabla\varphi\\|_{L_{h}^{2}(\Omega)}$
	$\displaystyle\leq\dfrac{C^{*}}{\epsilon_{0}/4}(\\|e_{u}^{k+1}-e_{u}^{k}\\|_{L_{h}^{2}(\Omega)}+C\tau)\\|\nabla\varphi\\|_{L_{h}^{2}(\Omega)}$
	$\displaystyle\leq C(\\|e_{u}^{k+1}\\|_{L_{h}^{2}(\Omega)}+\\|e_{u}^{k}\\|_{L_{h}^{2}(\Omega)}+C\tau)\\|\nabla\varphi\\|_{L_{h}^{2}(\Omega)}$
	$\displaystyle\leq Ch.$

Utilizing the property of the interpolation operator in (4.1), we have

$\displaystyle\|J_{5}\|$	$\displaystyle\leq\chi\\|u^{k}\nabla v^{k}\\|_{L^{2}(\Omega)}\\|\nabla(\varphi-\Phi)\\|_{L^{2}(\Omega)}$	(4.8)
	$\displaystyle\leq Ch\\|u^{k}\\|_{L^{\infty}(\Omega)}\\|\nabla v^{k}\\|_{L^{2}(\Omega)}\\|\Phi\\|_{H^{2}(\Omega)}$
	$\displaystyle\leq Ch\\|u^{k}\\|_{L^{\infty}(\Omega)}\\|u^{k}\\|_{L^{2}(\Omega)}\\|\Phi\\|_{H^{2}(\Omega)}$
	$\displaystyle\leq Ch.$

An application of the Cauchy-Schwarz inequality and the error estimates leads to

	$\displaystyle\|J_{6}\|=\|-\chi(u_{h}^{k}\nabla e_{v}^{k}+e_{u}^{k}\nabla v^{k},\nabla\varphi)\|$	(4.9)
$\displaystyle\leq$	$\displaystyle\chi(\\|u_{h}^{k}\\|_{L^{\infty}(\Omega)}\\|\nabla e_{v}^{k}\\|_{L^{2}(\Omega)}\\|\nabla\varphi\\|_{L^{2}(\Omega)}+\\|e_{u}^{k}\\|_{L^{2}(\Omega)}\\|\nabla v^{k}\\|_{L^{2}(\Omega)}\\|\nabla\varphi\\|_{L^{\infty}(\Omega)})$
$\displaystyle\leq$	$\displaystyle Ch,$

where $\tau\leq Ch$ has been used in the last inequality. An application of Lemma 4 to $J_{4}$ leads to

|J_{7}|=|\chi(u_{h}^{k}\nabla v_{h}^{k},\nabla\varphi)_{h}-\chi(u_{h}^{k}\nabla v_{h}^{k},\nabla\varphi)|=0.

(4.10)

Using the property of the interpolation operator $I_{h}$ as in the proof of the Theorem 6, we have

\displaystyle|J_{8}|=

\displaystyle|\left(u_{h}^{k}\nabla((I-I_{h})\log u_{h}^{k+1}),\nabla\varphi\right)_{h}|\leq Ch.

(4.11)

Finally, a substitution of (4.4)-(4.11) into (4.3) implies that

\frac{M^{k+1}-M^{k}}{\tau}\leq 4\theta-\dfrac{\alpha\chi}{2\pi}\theta^{2}+C_{1}\theta M^{k}+C_{2}\theta^{3/2}(M^{k}+Ch)^{1/2}+C_{0}h.

The proof is complete. ∎

Remark 14.

The positive constant $C_{0}$ in Lemma 13 depends on the regularity of the exact solutions.

We can then derive the following discrete analog of Theorem 3.

Theorem 15.

Assume that $\theta=(u,1)>8\pi/(\alpha\chi)$ . If $(u_{h}^{0},\varphi)_{h}$ is sufficiently small, $h$ and $\tau\leq Ch$ are sufficiently small, then the solution $(u_{h}^{k},v_{h}^{k})$ to the fully discrete scheme (3.3)-(3.4) will blow up in finite time, namely the maximal existence time $T_{\max}=k_{\max}\tau$ of the discrete solutions is finite.

Proof.

Obviously, we can derive the following inequality from Lemma 13

\displaystyle\frac{M^{k+1}-M^{k}}{\tau}\leq 4\theta-\dfrac{\alpha\chi}{2\pi}\theta^{2}+C_{1}\theta M^{k}+C_{2}\theta^{3/2}(M^{k})^{1/2}+C_{3}h^{1/2}.

Denote $\beta=-(4\theta-\frac{\alpha\chi}{2\pi}\theta^{2}+C_{3}h^{1/2})$ , we have the following inequality

\displaystyle\frac{M^{k+1}-M^{k}}{\tau}\leq-\beta+C_{1}\theta M^{k}+C_{2}\theta^{3/2}\sqrt{M^{k}}.

(4.12)

Under the condition that $\theta=(u,1)>8\pi/(\alpha\chi)$ , we can choose sufficiently small $h$ such that $\beta>0$ . Since $M^{0}=(u_{h}^{0},\varphi)_{h}$ is sufficiently small, we have

\displaystyle\beta_{0}:=\beta-C_{1}\theta M^{0}-C_{2}\theta^{3/2}\sqrt{M^{0}}>0.

We claim that the following inequality holds for all $k\in\mathbb{N}$

\displaystyle M^{k+1}\leq M^{k}-\tau\beta_{0}.

(4.13)

We prove the above inequality by induction. Using the inequality (4.12) for $k=0$ , we have

\displaystyle M^{1}\leq M^{0}-\tau(\beta-C_{1}\theta M^{0}-C_{2}\theta^{3/2}\sqrt{M^{0}})=M^{0}-\tau\beta_{0}.

Now assume that (4.13) holds for $k-1$ , we have

\displaystyle M^{k}\leq M^{k-1}-\tau\beta_{0}.

Notice that $M^{k}$ is decreasing about $k$ , we have

	$\displaystyle M^{k+1}$	$\displaystyle\leq M^{k}-\tau(\beta-C_{1}\theta M^{k}-C_{2}\theta^{3/2}\sqrt{M^{k}})$
		$\displaystyle\leq M^{k}-\tau(\beta-C_{1}\theta M^{0}-C_{2}\theta^{3/2}\sqrt{M^{0}})$
		$\displaystyle=M^{k}-\tau\beta_{0}.$

Next summing (4.13) over $k$ shows that

\displaystyle M^{k+1}\leq M^{0}-(k+1)\tau\beta_{0}.

Hence, if the solution $(u_{h}^{k},v_{h}^{k})$ exists for all $k\geq 0$ , then $M^{k+1}$ becomes negative provided that $T>M^{0}/\beta_{0}$ . This is a contradiction to the positivity of $M^{k}$ . Thus, the proof is complete. ∎

Remark 16.

Note that in the classical Keller-Segel system, the solution may blow up in finite time. Based on the error estimates, we prove that the numerical solution can also blow up under large initial value. There are several numerical examples in [46] to validate the blowup behavior of the numerical solution to the fully discrete scheme (3.3)-(3.4). The analysis of Theorem 15 depends on the regularity of the solution, it is very interesting whether we can still have similar results under weak regularity, we will continue to conduct on this issue in the future.

5. Conclusion

In this paper, we established error estimates for a fully discrete scheme proposed in [46] for the classical parabolic-elliptic Keller-Segel system, and showed that the numerical solution will blow up in finite time under some assumptions, similar to the situation for the exact solution of the classical parabolic-elliptic Keller-Segel system.

Acknowledgments

W. Chen is supported by the National Natural Science Foundation of China (NSFC) 12071090 and the work of J. Shen is partially supported by NSF Grant DMS-2012585 and AFOSR Grant FA9550-20-1-0309.

Appendix A Appendix

Lemma 17.

Denote $A=\alpha(-\Delta+I)^{-1},A_{h}=\alpha(-\Delta_{h}+I)^{-1}Q_{h}$ , where $Q_{h}$ is defined as $(Q_{h}u,\psi_{h})=(u,\psi_{h})_{h}$ , then the following estimate holds for all $u\in H^{2}(\Omega)$ ,

\displaystyle\|A_{h}u-Au\|_{L^{2}(\Omega)}+h\|A_{h}u-Au\|_{H^{1}(\Omega)}\leq Ch^{2}\|u\|_{H^{2}(\Omega)}.

(A.1)

Proof.

Let $w=Au$ and $w_{h}=A_{h}u$ , we have the following equations:

	$\displaystyle(\nabla w,\nabla\psi)+(w,\psi)=\alpha(u,\psi),\ \forall\psi\in H^{1},$
	$\displaystyle(\nabla w_{h},\nabla\psi)+(w_{h},\psi)=\alpha(u,\psi)_{h},\ \forall\psi\in X_{h}.$

Then we have

\displaystyle(\nabla(w_{h}-w),\nabla\psi)+(w_{h}-w,\psi)=\alpha((u,\psi)_{h}-(u,\psi)).

Taking $\psi=w_{h}-\psi_{h}$ and denoting $w_{h}-w$ by $w_{h}-\psi_{h}-(w-\psi_{h})$ in the above equation leads to

\displaystyle\|w_{h}-\psi_{h}\|_{H^{1}}^{2}\leq\|w-\psi_{h}\|_{H^{1}}\|w_{h}-\psi_{h}\|_{H^{1}}+\alpha((u,w_{h}-\psi_{h})_{h}-(u,w_{h}-\psi_{h})).

From [48] and Lemma 4, we have the following estimate

		$\displaystyle\|(u,w_{h}-\psi_{h})-(u,w_{h}-\psi_{h})_{h}\|$
	$\displaystyle=$	$\displaystyle\|(u-I_{h}u,w_{h}-\psi_{h})-(u-I_{h},w_{h}-\psi_{h})_{h}+(I_{h}u,w_{h}-\psi_{h})-(I_{h}u,w_{h}-\psi_{h})_{h}\|$
	$\displaystyle\leq$	$\displaystyle Ch\|u\|_{H^{1}}\\|w_{h}-\psi_{h}\\|_{L^{2}}+Ch^{2}\\|\nabla u\\|_{L^{2}}\\|\nabla(w_{h}-\psi_{h})\\|_{L^{2}}$
	$\displaystyle\leq$	$\displaystyle Ch\|u\|_{H^{1}}\\|w_{h}-\psi_{h}\\|_{H^{1}}.$

Combing the above estimates with the elliptic regularity estimate leads to

	$\displaystyle\\|w-w_{h}\\|_{H^{1}}$	$\displaystyle\leq C\inf_{\forall\psi_{h}\in X_{h}}\\|w-\psi_{h}\\|_{H^{1}}+Ch\|u\|_{H^{1}}$
		$\displaystyle\leq Ch\|w\|_{H^{2}}+Ch\|u\|_{H^{1}}\leq Ch\\|u\\|_{H^{1}}.$

The $L^{2}$ error estimate can be obtained by using the duality argument

\displaystyle\|w-w_{h}\|_{L^{2}}\leq Ch\|w-w_{h}\|_{H^{1}}+Ch^{2}|u|_{H^{2}}\leq Ch^{2}\|u\|_{H^{2}}.

Combing above estimates with the definitions of $A$ and $A_{h}$ shows that

		$\displaystyle\\|A_{h}u-Au\\|_{L^{2}(\Omega)}+h\\|A_{h}u-Au\\|_{H^{1}(\Omega)}$
	$\displaystyle=$	$\displaystyle\\|w_{h}-w\\|_{L^{2}(\Omega)}+h\\|w_{h}-w\\|_{H^{1}(\Omega)}$
	$\displaystyle\leq$	$\displaystyle Ch^{2}\\|u\\|_{H^{2}(\Omega)},$

which completes the proof. ∎

In order to obtain Lemma 8, we proceed in several steps. Firstly, we deal with $\hat{u}_{1}:=u+hf_{1}$ and construct $f_{1}$ as follows.

Lemma 18.

Assume that $\tau\leq Ch$ and $u$ is smooth enough, then there exists a bounded smooth function $f_{1}$ , such that $\hat{u}_{1}:=u+hf_{1}$ satisfies

\displaystyle\left(\overline{\partial}\hat{u}_{1}^{k},\varphi\right)_{h}+\left(\hat{u}_{1}^{k}\nabla I_{h}\log\hat{u}_{1}^{k+1},\nabla\varphi\right)_{h}-\chi\left(\hat{u}_{1}^{k}\nabla A_{h}\hat{u}_{1}^{k},\nabla\varphi\right)_{h}=\langle\mathcal{R}_{1}^{k},\varphi\rangle,

(A.2)

for all $\varphi\in X_{h}$ , $k\in\mathbb{N}$ , where $A_{h}$ is defined as in Lemma 17 and

\displaystyle|\langle\mathcal{R}^{k}_{1},\varphi\rangle|\leq Ch^{2}\|\varphi\|_{H^{1}}.

Proof.

From (A.2), we have the following equality:

	$\displaystyle(\dfrac{\partial\hat{u}_{1}}{\partial t},\varphi)+(\nabla\hat{u}_{1},\nabla\varphi)-\chi(\hat{u}_{1}\nabla A\hat{u}_{1},\nabla\varphi)$	(A.3)
$\displaystyle=$	$\displaystyle(\dfrac{\partial\hat{u}_{1}}{\partial t},\varphi)+(\nabla\hat{u}_{1},\nabla\varphi)-\chi(\hat{u}_{1}\nabla A\hat{u}_{1},\nabla\varphi)+\langle\mathcal{R}^{k}_{1},\varphi\rangle$
	$\displaystyle-\left(\overline{\partial}\hat{u}_{1}^{k},\varphi\right)_{h}-\left(\hat{u}_{1}^{k}\nabla I_{h}\log\hat{u}_{1}^{k+1},\nabla\varphi\right)_{h}+\chi\left(\hat{u}_{1}^{k}\nabla A_{h}\hat{u}_{1}^{k},\nabla\varphi\right)_{h},$

where the operators $A$ and $A_{h}$ are defined as in Lemma 17. Moreover, the left hand of (A.3) can be rewritten as

		$\displaystyle(\dfrac{\partial\hat{u}_{1}}{\partial t},\varphi)+(\nabla\hat{u}_{1},\nabla\varphi)-\chi(\hat{u}_{1}A\hat{u}_{1},\nabla\varphi)$
	$\displaystyle=$	$\displaystyle h((\dfrac{\partial f_{1}}{\partial t},\varphi)+(\nabla f_{1},\nabla\varphi)-\chi(uAf_{1},\nabla\varphi)-\chi(f_{1}Au,\nabla\varphi))$
		$\displaystyle-h^{2}(f_{1}Af_{1},\nabla\varphi),$

where equation (3.24) has been used.

Step 1: Construction for $f_{1}$ . For any $\varphi\in H^{1}$ , define $\langle\mathcal{R}^{k}_{0}(u),\varphi\rangle$ as

	$\displaystyle\langle\mathcal{R}^{k}_{0}(u),\varphi\rangle:=$	$\displaystyle(\dfrac{\partial u}{\partial t}-\overline{\partial}u^{k},\varphi)+(\nabla(u-u^{k+1}),\nabla\varphi)-(u^{k}\nabla(I_{h}-I)(\log u^{k+1}),\nabla\varphi)_{h}$
		$\displaystyle-((u^{k}-u^{k+1})\nabla\log u^{k+1},\nabla\varphi)_{h}+(\nabla u^{k+1},\nabla\varphi)-(\nabla u^{k+1},\nabla\varphi)_{h}$
		$\displaystyle+\chi(u^{k}\nabla A_{h}u^{k}-u\nabla A_{h}u,\nabla\varphi)_{h}+\chi(u\nabla(A_{h}-A)u,\nabla\varphi)$
	$\displaystyle:=$	$\displaystyle\langle\mathcal{N}(u),\varphi\rangle+\chi(u^{k}\nabla A_{h}u^{k}-u\nabla A_{h}u,\nabla\varphi)_{h}+\chi(u\nabla(A_{h}-A)u,\nabla\varphi),$

we can show that $\frac{1}{h}\mathcal{R}^{k}_{0}(u)\in H^{-1}$ is well defined. Using the Cauchy-Schwarz inequality and the property of the interpolation, we obtain

	$\displaystyle\|(\nabla u^{k+1},\nabla\varphi)-(\nabla u^{k+1},\nabla\varphi)_{h}\|=$	$\displaystyle(\nabla(I-I_{h})u^{k+1},\nabla\varphi)-(\nabla(I-I_{h})u^{k+1},\nabla\varphi)_{h}$
	$\displaystyle\leq$	$\displaystyle Ch\\|u^{k+1}\\|_{H^{2}}\\|\nabla\varphi\\|_{L^{2}},$

then the following estimate holds for $\langle\mathcal{N}(u),\varphi\rangle$ :

\displaystyle\langle\mathcal{N}(u),\varphi\rangle\leq C_{4}h\|\varphi\|_{H^{1}},

(A.4)

and the positive constant

C_{4}=C_{4}(\|u\|_{W^{2,\infty}(0,T;H^{2})}+\|u\|_{W^{1,\infty}(0,T;L^{\infty})}\|\log u\|_{L^{\infty}(0,T;H^{2})}).

An application of Lemma 17 leads to

	$\displaystyle\|(u\nabla(A_{h}-A)u,\nabla\varphi)\|\leq$	$\displaystyle\\|u\\|_{L^{\infty}}\\|\nabla(A_{h}u-Au)\\|_{L^{2}}\\|\nabla\varphi\\|_{L^{2}}$
	$\displaystyle\leq$	$\displaystyle Ch\\|u\\|_{L^{\infty}}\\|u\\|_{H^{2}}\\|\nabla\varphi\\|_{L^{2}}.$

Combining (A.4) with few direct calculations shows that the following estimate holds for $\langle\mathcal{R}^{k}_{0}(u),\varphi\rangle$ :

\displaystyle|\langle\mathcal{R}^{k}_{0}(u),\varphi\rangle|\leq C_{5}h\|\varphi\|_{H^{1}},

where

C_{5}=C_{5}(\|u\|_{W^{2,\infty}(0,T;H^{2})}+\|u\|_{W^{1,\infty}(0,T;L^{\infty})}\|\log u\|_{L^{\infty}(0,T;H^{2})}+\|u\|_{W^{1,\infty}(0,T;L^{\infty})}^{2}),

then $\frac{1}{h}\langle\mathcal{R}^{k}_{0}(u),\varphi\rangle$ is well-defined. Combining $\langle\mathcal{R}^{k}_{0}(u),\varphi\rangle$ with (A.3) leads to the following linear partial differential equation for $f_{1}$ :

\displaystyle(\dfrac{\partial f_{1}}{\partial t},\varphi)+(\nabla f_{1},\nabla\varphi)-\chi(u\nabla Af_{1},\nabla\varphi)-\chi(f_{1}\nabla Au,\nabla\varphi)=\dfrac{1}{h}\langle\mathcal{R}^{k}_{0}(u),\varphi\rangle,

(A.5)

for all $t\in(t_{k},t_{k+1}]$ . From [17, Chapter7.1, Theorem 3], there exists a weak solution $f_{1}$ of (A.5). In addition, from [17, Chapter7.1, Theorem 7], suppose that $u$ is smooth enough such that $\mathcal{R}_{0}^{k}(u)$ is smooth enough in $\overline{\Omega}\times[t_{k},t_{k+1}]$ , and the $m^{th}$ -order compatibility conditions hold for $m=0,1,\cdots$ , then the problem (A.5) has a smooth enough solution $f_{1}$ in $\overline{\Omega}\times[t_{k},t_{k+1}]$ .

Step 2: Construction for $\langle\mathcal{R}^{k}_{1}(\hat{u}_{1}),\varphi\rangle$ . Let $\langle\mathcal{R}^{k}_{1}(\hat{u}_{1}),\varphi\rangle$ be

	$\displaystyle\langle\mathcal{R}^{k}_{1}(\hat{u}_{1}),\varphi\rangle:=$	$\displaystyle h\langle\mathcal{N}(f_{1}),\varphi\rangle+\chi h((u\nabla(A_{h}-A)f_{1}+f_{1}\nabla(A_{h}-A)u,\nabla\varphi))$
	$\displaystyle+$	$\displaystyle\chi h((u^{k}\nabla A_{h}f_{1}^{k}-u\nabla A_{h}f_{1},\nabla\varphi)_{h}+(f_{1}^{k}\nabla A_{h}u^{k}-f_{1}\nabla A_{h}u,\nabla\varphi)_{h})$
	$\displaystyle-$	$\displaystyle(u^{k}\nabla I_{h}(\log(u^{k+1}+hf_{1}^{k+1})-\log u^{k+1}-\dfrac{hf_{1}^{k+1}}{u^{k+1}}),\nabla\varphi)_{h}$
	$\displaystyle-$	$\displaystyle h(u^{k}\nabla(I_{h}-I)\dfrac{f_{1}^{k+1}}{u^{k+1}},\nabla\varphi)_{h}-h((u^{k}-u^{k+1})\nabla\dfrac{f_{1}^{k+1}}{u^{k+1}},\nabla\varphi)_{h}$
	$\displaystyle+$	$\displaystyle\left(\overline{\partial}u^{k},\varphi\right)-\left(\overline{\partial}u^{k},\varphi\right)_{h}+h^{2}(f_{1}^{k}\nabla\dfrac{f_{1}^{k+1}}{u^{k+1}},\nabla\varphi)+\chi h^{2}(f_{1}\nabla Af_{1},\nabla\varphi)$
	$\displaystyle+$	$\displaystyle\chi((u\nabla A_{h}u,\nabla\varphi)_{h}-(u\nabla A_{h}u,\nabla\varphi)).$

Similarly, the following estimate holds for $\langle\mathcal{N}(f_{1}),\varphi\rangle$ as discussed in (A.4):

\displaystyle\langle\mathcal{N}(f_{1}),\varphi\rangle\leq C_{6}h\|\varphi\|_{H^{1}},

(A.6)

where the positive constant

C_{6}=C_{6}(\|f_{1}\|_{W^{2,\infty}(0,T;H^{2})}+\|f_{1}\|_{W^{1,\infty}(0,T;L^{\infty})}\|\log u\|_{L^{\infty}(0,T;H^{2})}).

Combining (A.6) with few calculations yields the following estimate for $\langle\mathcal{R}^{k}_{1}(\hat{u}_{1}),\varphi\rangle$ :

\displaystyle|\langle\mathcal{R}^{k}_{1}(\hat{u}_{1}),\varphi\rangle|\leq C_{7}h^{2}\|\varphi\|_{H^{1}},

where the positive constant

	$\displaystyle C_{7}=C_{7}(\\|u\\|_{W^{1,\infty}(0,T;H^{2})}+(\\|u\\|_{W^{1,\infty}(0,T;L^{\infty})}+\\|f_{1}\\|_{L^{\infty}(0,T;L^{\infty})})\\|\frac{f_{1}}{u}\\|_{L^{\infty}(0,T;H^{2})}$
	$\displaystyle+\\|u\\|_{W^{1,\infty}(0,T;L^{\infty})}\\|f_{1}\\|_{W^{1,\infty}(0,T;L^{\infty})}+\\|u\\|_{L^{\infty}(0,T;H^{2})}^{2}+\\|f_{1}\\|_{L^{\infty}(0,T;L^{\infty})}^{2})+C_{6},$

then the $O(h^{2})$ consistency for $\hat{u}_{1}=u+hf_{1}$ is obtained, which leads to Lemma 18. ∎

Remark 19.

Taking $\varphi=1$ in (A.2) leads to $(\hat{u}_{1}^{k},1)_{h}=(\hat{u}_{1}^{0},1)_{h}$ for $k\in\mathbb{N}$ , i.e., $\hat{u}_{1}$ preserves mass conservation property. Choosing suitable initial condition $f_{1}(x,0)=0$ such that $(f_{1}(x,0),1)_{h}=0$ , we obtain $(\hat{u}_{1}^{k},1)_{h}=(u^{0},1)_{h}$ .

After repeated application of the perturbation argument as illustrated in Lemma 18, Lemma 8 can be proved.

Proof of Lemma 8.

The duality product $\frac{1}{h^{2}}\langle\mathcal{R}^{k}_{1}(\hat{u}_{1}),\varphi\rangle$ is well-defined since the fact that $\frac{1}{h^{2}}\langle\mathcal{R}^{k}_{1}(\hat{u}_{1}),\varphi\rangle$ is uniformly bounded as $h\rightarrow 0,\tau\rightarrow 0$ and $\tau\leq Ch$ . We can construct $f_{2}$ by solving the following linear partial differential equation:

\displaystyle(\dfrac{\partial f_{2}}{\partial t},\varphi)+(\nabla f_{2},\nabla\varphi)-\chi(u\nabla Af_{2},\nabla\varphi)-\chi(f_{2}\nabla Au,\nabla\varphi)=\dfrac{1}{h^{2}}\langle\mathcal{R}^{k}_{1}(\hat{u}_{1}),\varphi\rangle,

(A.7)

for all $t\in(t_{k},t_{k+1}]$ . As discussed in Step 1 above, the problem (A.7) has a smooth enough solution $f_{2}$ in $\overline{\Omega}\times[t_{k},t_{k+1}]$ .

By repeated application of the methods in Step 2 above, we can construct $\langle\mathcal{R}^{k}_{2}(\hat{u}_{2}),\varphi\rangle$ by $\hat{u}_{2}:=\hat{u}_{1}+h^{2}f_{2}$ such that the $O(h^{3})$ consistency for $\hat{u}_{2}$ is arrived:

\displaystyle\left(\overline{\partial}\hat{u}_{2}^{k},\varphi\right)_{h}+\left(\hat{u}_{2}^{k}\nabla I_{h}\log\hat{u}_{2}^{k+1},\nabla\varphi\right)_{h}-\chi\left(\hat{u}_{2}^{k}\nabla A_{h}\hat{u}_{2}^{k},\nabla\varphi\right)_{h}=\langle\mathcal{R}^{k}_{2}(\hat{u}_{2}),\varphi\rangle,

(A.8)

for all $\varphi\in X_{h}$ , $k\in\mathbb{N}$ , where

\displaystyle|\langle\mathcal{R}^{k}_{2}(\hat{u}_{2}),\varphi\rangle|\leq Ch^{3}\|\varphi\|_{H^{1}},

where $C$ is a positive constant depending on the derivatives of $\hat{u}_{2}$ , such that $\frac{1}{h^{3}}\langle\mathcal{R}^{k}_{2}(\hat{u}_{2}),\varphi\rangle$ is well-defined.

Again, by using Step 1 in Lemma 18, the correction function $f_{3}$ can be constructed by the following linear partial differential equation:

\displaystyle(\dfrac{\partial f_{3}}{\partial t},\varphi)+(\nabla f_{3},\nabla\varphi)-\chi(u\nabla Af_{3},\nabla\varphi)-\chi(f_{3}\nabla Au,\nabla\varphi)=\dfrac{1}{h^{3}}\langle\mathcal{R}^{k}_{2}(\hat{u}_{2}),\varphi\rangle,

(A.9)

for all $t\in(t_{k},t_{k+1}]$ , and the problem (A.9) has a smooth enough solution $f_{3}$ in $\overline{\Omega}\times[t_{k},t_{k+1}]$ .

Combing equations (A.5),(A.7) and (A.9) leads to

	$\displaystyle(\dfrac{\partial\hat{u}}{\partial t},\varphi)+(\nabla\hat{u},\nabla\varphi)-\chi(\hat{u}A\hat{u},\nabla\varphi)$	(A.10)
$\displaystyle=$	$\displaystyle\langle\mathcal{R}_{0}^{k}(u),\varphi\rangle+\langle\mathcal{R}_{1}^{k}(\hat{u}_{1}),\varphi\rangle+\langle\mathcal{R}_{2}^{k}(\hat{u}_{2}),\varphi\rangle-h^{4}\chi((f_{2}\nabla Af_{2},\nabla\varphi)$
	$\displaystyle+(f_{1}\nabla Af_{3}+f_{3}\nabla Af_{1},\nabla\varphi))-h^{5}\chi(f_{2}\nabla Af_{3}+f_{3}\nabla Af_{2},\nabla\varphi)$
	$\displaystyle-h^{6}\chi(f_{3}\nabla Af_{3},\nabla\varphi).$

Denote $\langle\mathcal{R}^{k}(\hat{u}),\varphi\rangle$ as follows

	$\displaystyle\langle\mathcal{R}^{k}(\hat{u}),\varphi\rangle$	$\displaystyle:=h^{3}((\overline{\partial}f_{3}^{k}-\dfrac{\partial f_{3}}{\partial t},\varphi)+(\overline{\partial}f_{3}^{k},\varphi)_{h}-(\overline{\partial}f_{3}^{k},\varphi)$
		$\displaystyle+(f_{3}^{k}\nabla(I_{h}-I)\log\hat{u}_{2}^{k+1},\varphi)_{h}-((u^{k+1}-u^{k})\nabla\dfrac{f_{3}^{k+1}}{\hat{u}_{2}^{k+1}},\nabla\varphi)_{h}$
		$\displaystyle+(\nabla f_{3}^{k+1},\nabla\varphi)_{h}-(\nabla f_{3}^{k+1},\nabla\varphi)+((f_{3}^{k}-f_{3}^{k+1})\nabla\log\hat{u}_{2}^{k+1},\nabla\varphi)_{h})$
		$\displaystyle+h^{6}(f_{3}^{k}\nabla\dfrac{f_{3}^{k+1}}{\hat{u}_{2}^{k+1}},\nabla\varphi)_{h}+(\hat{u}^{k}\nabla(\log\hat{u}^{k+1}-\log\hat{u}_{2}^{k+1}-\dfrac{h^{3}f_{3}^{k+1}}{\hat{u}_{2}^{k+1}}),\nabla\varphi)_{h}$
		$\displaystyle+(\hat{u}_{2}^{k}\nabla(\log\hat{u}^{k+1}-\log\hat{u}_{1}^{k+1}-\dfrac{h^{2}f_{2}^{k+1}}{\hat{u}_{1}^{k+1}}),\nabla\varphi)_{h}+h^{4}(f_{2}^{k}\nabla\dfrac{f_{2}^{k+1}}{\hat{u}_{1}^{k+1}},\nabla\varphi)_{h}$
		$\displaystyle+h^{2}((\overline{\partial}f_{2}^{k},\varphi)_{h}-(\overline{\partial}f_{2}^{k},\varphi))+h^{3}\chi((f_{3}^{k}\nabla A_{h}\hat{u}_{2}^{k}+\hat{u}_{2}^{k}\nabla A_{h}f_{3}^{k},\nabla\varphi)_{h}$
		$\displaystyle-(f_{3}\nabla A_{h}\hat{u}_{2}+\hat{u}_{2}\nabla A_{h}f_{3},\nabla\varphi))+h^{6}\chi((f_{3}^{k}\nabla A_{h}f_{3}^{k},\nabla\varphi)_{h}$
		$\displaystyle-(f_{3}\nabla A_{h}f_{3},\nabla\varphi))+h^{4}\chi((f_{2}^{k}\nabla A_{h}f_{2}^{k},\nabla\varphi)_{h}-(f_{2}\nabla A_{h}f_{2},\nabla\varphi))$
		$\displaystyle+h^{2}\chi((f_{2}\nabla A_{h}\hat{u}_{1}+\hat{u}_{1}\nabla A_{h}f_{2},\nabla\varphi)_{h}-(f_{2}\nabla A_{h}\hat{u}_{1}+\hat{u}_{1}\nabla A_{h}f_{2},\nabla\varphi)).$

Combining above with few direct calculations shows the following estimate for $\langle\mathcal{R}^{k}(\hat{u}),\varphi\rangle$

\displaystyle|\langle\mathcal{R}^{k}(\hat{u}),\varphi\rangle|\leq Ch^{4}\|\varphi\|_{H^{1}},

where $C$ depending on the derivatives of $u$ , then the $\frac{1}{h^{4}}\langle\mathcal{R}^{k}(\hat{u}),\varphi\rangle$ is well defined and the $O(h^{4})$ consistency holds for $\hat{u}=u+hf_{1}+h^{2}f_{2}+h^{3}f_{3}$ , which leads to Lemma 8. ∎

Remark 20.

Similarly, taking $\varphi=1$ in (3.5) leads to $(\hat{u}^{k},1)_{h}=(\hat{u}^{0},1)_{h}$ for $k\in\mathbb{N}$ . Choosing the initial condition $f_{i}(x,0)=0$ such that $(f_{i}(x,0),1)_{h}=0$ ( $i=1$ , $2$ , $3$ ), we obtain $(\hat{u}^{k},1)_{h}=(u^{0},1)_{h}$ .

References

[1] R. Adams and J. Fournier, Sobolev Spaces, 2nd Edition, Academic Press, Singapore, 2003.
[2] A. Blanchet and P. Laurençot, The parabolic-parabolic Keller-Segel system with critical diffusion as a gradient flow in $\mathbb{R}^{d},d\geq 3$ , Communications in Partial Differential Equations, 38(4): 658-686, 2012.
[3] A. Blanchet, J. Dolbeault and B. Perthame, Two-dimensional Keller-Segel model: Optimal critical mass and qualitative properties of the solutions, Electronic Journal of Differential Equations, 2006(44): 285–296, 2016.
[4] S. C. Brenner and L. R. Scott, The Mathematical Theory of Finite Element Methods, 3rd Edition, Springer, New York, 2007.
[5] V. Calvez and L. Corrias, The parabolic-parabolic Keller-Segel model in $\mathbb{R}^{2}$ , Communications in Mathematical Sciences, 6(2): 417–447, 2008.
[6] V. Calvez, L. Corrias and M. A. Ebde, Blow-up, concentration phenomenon and global existence for the Keller-Segel model in high dimension, Journal of Differential Equations, 267(11): 561–584, 2019.
[7] W. Chen, C. Wang, X. Wang and S. Wise, Positivity-preserving, energy stable numerical schemes for the Cahn-Hilliard equation with logarithmic potential, Journal of Computational Physics X, 3: 100031, 2019.
[8] J. I. Diaz and T. Nagai, Symmetrization in a parabolic-elliptic system related to chemotaxis, Advances in Mathematical Sciences and Applications, 5(2): 659-680, 1995.
[9] J. Dolbeault and B. Perthame, Optimal critical mass in the two-dimensional Keller-Segel model in $\mathbb{R}^{2}$ . Retour Au Numéro, 339(9): 611–616, 2004.
[10] L. Dong, C. Wang, H. Zhang and Z. Zhang, A positivity-preserving, energy stable and convergent numerical scheme for the Cahn-Hilliard equation with a Flory-Huggins-deGennes energy, Communications in Mathematical Sciences, 17(4): 921-939, 2019.
[11] L. Dong, C. Wang, H. Zhang and Z. Zhang, A positivity-preserving second-order BDF scheme for the Cahn-Hilliard equation with variable interfacial parameters, Communications in Computational Physics, 28(3): 967-998, 2020.
[12] L. Dong, C. Wang, S. Wise and Z. Zhang, A positivity preserving, energy stable scheme for the ternary Cahn-Hilliard system with the singular interfacial parameters, Journal of Computational Physics, 442: 110451, 2021.
[13] C. Duan, C. Liu, C. Wang and X. Yue, Convergence analysis of a numerical Scheme for the porous medium equation by an energetic variational approach, Numerical Mathematics: Theory, Methods and Applications, 13(1): 63-80, 2020.
[14] C. Duan, W. Chen, C. Liu, C. Wang and S. Zhou, Convergence analysis of structure-preserving numerical methods for nonlinear Fokker-Planck equations with nonlocal interactions, Mathematical Methods in the Applied Sciences, accepted and in press, 2021.
[15] Y. Epshteyn and A. Kurganov, New interior penalty discontinuous Galerkin methods for the Keller-Segel Chemotaxis model, SIAM Journal on Numerical Analysis, 47(1): 386–408, 2008.
[16] Y. Epshteyn and A. Izmirlioglu, Fully discrete analysis of a discontinuous finite element method for the Keller-Segel chemotaxis model, Journal of Scientific Computing, 40(1–3): 211–256, 2009.
[17] L. C. Evans, Partial Differential Equations. AMS Graduate studies in Mathematics. AMS: Providence, RI, 1998.
[18] F. Filbet, A finite volume scheme for the Patlak–Keller–Segel chemotaxis model, Numerische Mathematik, 104(4): 457–488, 2006.
[19] K. Fujie, M. Winkler and T. Yokota, Boundedness of solutions to parabolic-elliptic chemotaxis-growth systems with signal-dependent sensitivity, Mathematical Methods in the Applied Sciences, 38(6): 1212–1224, 2015.
[20] J. Fuhrmann, J. Lankeit and M. Winkler, A double critical mass phenomenon in a no-flux-Dirichlet Keller-Segel system, https://arxiv.org/abs/2101.06748, 2021.
[21] H. Gajewski, K. Zacharias, and K. Gröger, Global behavior of a reaction-diffusion system modelling chemotaxis, Mathematische Nachrichten, 195(1): 77–114, 1998.
[22] L. Guo, X. Li, and Y. Yang, Energy dissipative local discontinuous Galerkin methods for keller–segel chemotaxis model, Journal of Scientific Computing, 78: 1387–1404, 2019.
[23] A. Gurusamy, K. Balachandran, Finite element method for solving Keller-Segel chemotaxis system with cross-diffusion, International Journal of Dynamics and Control, 6:539–549,2018.
[24] D. Horstmann, The nonsymmetric case of the Keller-Segel model in chemotaxis: Some recent results, Nonlinear Differential Equations and Applications, 8: 399–423, 2001.
[25] D. Horstmann, From 1970 until present: The Keller-Segel model in chemo-taxis and its consequences, Jahresbericht der Deutschen Mathematiker-Vereinigung, 106(2): 51–69, 2004.
[26] A. Jüengel and O. Leingang, Blow-up of solutions to semi-discrete parabolic-elliptic Keller-Segel models, Discrete and Continuous Dynamical Systems, 24(9): 4755–4782, 2019.
[27] E. F. Keller and L. A. Segel, Initiation on slime mold aggregation viewed as instability. Journal of Theoretical Biology, 26(3): 399–415, 1970.
[28] E. F. Keller and L. A. Segel, Model for chemotaxis, Journal of Theoretical Biology, 30(2): 225–234, 1971.
[29] H. Kozono and Y. Sugiyama, Local existence and finite time blow-up of solutions in the 2-D Keller-Segel system, Journal of Evolution Equations, 8: 353–378, 2008.
[30] X. Li, C. W. Shu, and Y. Yang, Local discontinuous Galerkin method for the Keller-Segel chemotaxis model, Journal of Scientific Computing, 73: 943–967, 2017.
[31] W. Li, W. Chen, C. Wang, Y. Yan, and R. He, A second order energy stable linear scheme for a thin film model without slope selection, Journal of Scientific Computing, 76: 1905–1937, 2018.
[32] C. Liu, C. Wang, S. M. Wise, X. Yue, and S. Zhou, A positivity-preserving, energy stable and convergent numerical scheme for the Poisson-Nernst-Planck system, Mathematics of Computation, 90(331): 2071-2106, 2021.
[33] C. Liu, C. Wang and Y. Wang, A structure-preserving, operator splitting scheme for reaction-diffusion equations with detailed balance, Journal of Computational Physics, 436: 110253, 2021.
[34] J. G. Liu, L. Wang and Z. Zhou, Positivity-preserving and asymptotic preserving method for 2D Keller-Segal equations, Mathematics of Computation, 87(311): 1165–1189, 2018.
[35] T. Nagai, Behavior of solutions to a parabolic-elliptic system modelling chemotaxis, Journal of the Korean Mathematical Society, 37(5): 721–732, 2000.
[36] T. Nagai, T. Senba, and T. Suzuki, Concentration behavior of blow-up solutions for a simplified system of chemotaxis (variational problems and related topics), Rims Kokyuroku, 1181: 140–176, 2001.
[37] T. Nagai, Blowup of nonradial solutions to parabolic–elliptic systems modeling chemotaxis in two-dimensional domains, Journal of Inequalities and Applications, 6: 37-55, 2001.
[38] T. Nagai, Global existence and blowup of solutions to a chemotaxis system, Nonlinear Analysis: Theory, Methods and Applications, 47(2): 777–787, 2001.
[39] T. Nagai, T. Senba, and K. Yoshida, Application of the trudinger-moser inequality to a parabolic system of chemotaxis, Funkc Ekvacioj, 40(3): 411–433, 1997.
[40] L. E. Payne and P. W. Schaefer, Lower bounds for blow-up time in parabolic problems under dirichlet conditions, Journal of Mathematical Analysis and Applications, 328: 1196–1205, 2007.
[41] Y. Qian, C. Wang and S. Zhou, A positive and energy stable numerical scheme for the Poisson-Nernst-Planck-Cahn-Hilliard equations with steric interactions, Journal of Computational Physics, 426: 109908, 2021.
[42] N. Saito and T. Suzuki, Notes on finite difference schemes to a parabolic-elliptic system modelling chemotaxis, Applied Mathematics and Computation, 171(1): 72-90, 2005.
[43] N. Saito, Conservative upwind finite-element method for a simplified Keller-Segel system modelling chemotaxis, IMA Journal of Numerical Analysis, 27(2): 332–365, 2007.
[44] N. Saito, Error analysis of a conservative finite-element approximation for the Keller-Segel system of chemotaxis, Communications on Pure and Applied Analysis, 11(1): 339–364, 2012.
[45] J. Shen, Long time stability and convergence for fully discrete nonlinear Galerkin methods, Applicable Analysis, 38: 201–229, 1990.
[46] J. Shen and J. Xu, Unconditionally bound preserving and energy dissipative schemes for a class of Keller–Segel equations, SIAM Journal on Numerical Analysis, 58(3): 1674–1695, 2020.
[47] T. Suzuki, Free Energy and Self-Interacting Particles, Progr. Nonlinear Differential Equations Appl, 1st Edition, Birkhäuser Boston, 2005.
[48] V. Thomee, Galerkin Finite Element Methods for Parabolic Problems, 2nd Edition, Springer-Verlag, 2006.
[49] D. Wei, Global well-posedness and blow-up for the 2-D Patlak-Keller-Segel equation, Journal of Functional Analysis, 274: 388–401, 2018.
[50] Y. Yan, W. Chen, C. Wang and S. Wise, A second-order energy stable BDF numerical scheme for the Cahn-Hilliard equation, Communications in Computational Physics, 23(2): 572–602, 2018.
[51] M. Yuan, W. Chen, C. Wang, S. Wise and Z. Zhang, An energy stable finite element scheme for the three-component Cahn-Hilliard-type model for macromolecular microsphere composite hydrogels, Journal of Scientific Computing, 87:78, 2021.
[52] G. Zhou and N. Saito, Finite volume methods for a Keller–Segel system: discrete energy, error estimates and numerical blow-up analysis, Numerische Mathematik, 135: 1-47, 2016.

	$\displaystyle\|(\tilde{u}^{k},\tilde{\mu}_{u}^{k})_{h}\|$	$\displaystyle\leq\\|\tilde{u}^{k}\\|_{-1}\\|\nabla\tilde{\mu}_{u}^{k}\\|_{L_{h}^{2}}\leq C\\|\tilde{u}^{k}\\|_{L_{h}^{2}}\\|\nabla\tilde{\mu}_{u}^{k}\\|_{L_{h}^{2}}$		(3.17)
		$\displaystyle\leq\dfrac{6C}{\epsilon_{0}\tau}\\|\tilde{u}^{k}\\|_{L_{h}^{2}}^{2}+\dfrac{\epsilon_{0}\tau}{24}\\|\nabla\tilde{\mu}_{u}^{k}\\|_{L_{h}^{2}}^{2}.$		(3.17)

$\displaystyle\tau\\|\nabla\cdot(\tilde{u}^{k}\nabla\mathcal{V}_{u}^{k})\\|_{L_{h}^{2}}$	$\displaystyle\leq\tau(\\|\tilde{u}^{k}\Delta\mathcal{V}_{u}^{k}\\|_{L_{h}^{2}}+\\|\nabla\tilde{u}^{k}\nabla\mathcal{V}_{u}^{k}\\|_{L_{h}^{2}})$	(3.21)
	$\displaystyle\leq\tau(\\|\Delta\mathcal{V}_{u}^{k}\\|_{\infty}\\|\tilde{u}^{k}\\|_{L_{h}^{2}}+\\|\nabla\mathcal{V}_{u}^{k}\\|_{\infty}\\|\nabla\tilde{u}^{k}\\|_{L_{h}^{2}})$
	$\displaystyle\leq C\tau(\\|\tilde{u}^{k}\\|_{L_{h}^{2}}+\\|\nabla\tilde{u}^{k}\\|_{L_{h}^{2}})$
	$\displaystyle\leq Ch^{15/4},$

$\displaystyle\tau\\|\nabla\cdot(u_{h}^{k}\nabla\tilde{\mu}_{u}^{k})\\|_{L_{h}^{2}}$	$\displaystyle\leq\tau(\\|\nabla u_{h}^{k}\nabla\tilde{\mu}_{u}^{k}\\|_{L_{h}^{2}}+\\|u_{h}^{k}\Delta\tilde{\mu}_{u}^{k}\\|_{L_{h}^{2}})$	(3.22)
	$\displaystyle\leq\tau(\\|\nabla u_{h}^{k}\\|_{\infty}\\|\nabla\tilde{\mu}_{u}^{k}\\|_{L_{h}^{2}}+\\|u_{h}^{k}\\|_{\infty}\\|\Delta\tilde{\mu}_{u}^{k}\\|_{L_{h}^{2}})$
	$\displaystyle\leq C^{*}\tau(\\|\nabla\tilde{\mu}_{u}^{k}\\|_{L_{h}^{2}}+\\|\Delta\tilde{\mu}_{u}^{k}\\|_{L_{h}^{2}})$
	$\displaystyle\leq Ch^{11/4},$

	$\displaystyle\\|\nabla\sigma_{v}^{k}\\|_{L^{2}(\Omega)}^{2}+\\|\sigma_{v}^{k}\\|_{L^{2}(\Omega)}^{2}\leq$	$\displaystyle\ \alpha\\|\sigma_{u}^{k}\\|_{L^{2}(\Omega)}\\|\sigma_{v}^{k}\\|_{L^{2}(\Omega)}+Ch^{2}\\|\nabla u_{h}^{k}\\|_{L^{2}(\Omega)}\\|\nabla\sigma_{v}^{k}\\|_{L^{2}(\Omega)}$
		$\displaystyle+Ch\\|v\\|_{H^{2}(\Omega)}\\|\sigma_{v}^{k}\\|_{L^{2}(\Omega)}+Ch^{2}\\|u\\|_{H^{2}(\Omega)}\\|\sigma_{v}^{k}\\|_{L^{2}(\Omega)}$
	$\displaystyle\leq$	$\displaystyle\ C\\|\sigma_{u}^{k}\\|_{L^{2}(\Omega)}^{2}+\\|\sigma_{v}^{k}\\|_{L^{2}(\Omega)}^{2}+\dfrac{1}{2}\\|\nabla\sigma_{v}^{k}\\|_{L^{2}(\Omega)}^{2}+Ch^{2},$

	$\displaystyle\|I_{1}\|$	$\displaystyle=\|\left(\overline{\partial}u^{k}-u_{t}^{k+1},\ \sigma_{u}^{k+1}\right)\|\leq C\tau^{1/2}\Big{(}\int_{t_{k}}^{t_{k+1}}\\|u_{tt}\\|_{L^{2}(\Omega)}^{2}\mathrm{d}s\Big{)}^{1/2}\\|\sigma_{u}^{k+1}\\|_{L^{2}(\Omega)}$		(3.37)
		$\displaystyle\leq C\tau\int_{t_{k}}^{t_{k+1}}\\|u_{tt}\\|_{L^{2}(\Omega)}^{2}\mathrm{d}s+\dfrac{1}{20}\\|\nabla\sigma_{u}^{k+1}\\|_{L^{2}(\Omega)}^{2}.$		(3.37)

ERROR ESTIMATES AND BLOW-UP ANALYSIS OF A FINITE-ELEMENT APPROXIMATION FOR THE parabolic-elliptic Keller-Segel system

Abstract.

Key words and phrases:

1991 Mathematics Subject Classification:

1. Introduction.

2. The Keller-Segel equations

Theorem 1.

Lemma 2.

Theorem 3.

3. The fully discrete scheme and error estimates

Lemma 4.

Lemma 5.

Theorem 6.

Remark 7.

3.1. Higher order consistent approximation to (3.3)-(3.4)

Lemma 8.

Remark 9.

3.2. A rough error estimate

Lemma 10.

Lemma 11.

Lemma 12.

Proof.

3.3. Recovery of the assumption (3.11)

3.4. Proof of Theorem 6

4. Finite-time blowup

Lemma 13.

Proof.

Remark 14.

Theorem 15.

Proof.

Remark 16.

5. Conclusion

Acknowledgments

Appendix A Appendix

Lemma 17.

Proof.

Lemma 18.

Proof.

Remark 19.

Proof of Lemma 8.

Remark 20.

References

ERROR ESTIMATES AND BLOW-UP ANALYSIS OF
A FINITE-ELEMENT APPROXIMATION FOR THE
parabolic-elliptic Keller-Segel system