High-accuracy positivity-preserving numerical method for Keller-Segel model

1.   Introduction

In this paper, the finite-difference method (FDM) is applied to solve the following two-dimensional (2D) Keller-Segel model ^[1,2]:

The model (1.1) describes the evolutionary process of cell density $u(x, y, t)$ and a chemical stimulus (chemoattractant) concentration $v(x, y, t)$ over a time $t$ and location $(x, y)$, where $\Omega = \{(x, y)|a\leqslant x, y \leqslant b\} \subset \mathcal{R}^2$ is a convex bounded domain and $a$ and $b$ are constants. $\nabla$ and $\Delta$ are gradient and Laplacian operators, respectively. $\chi > 0$ represents the chemotactic sensitivity constant; $d > 0$ denotes the diffusion rate of cells. In addition, the initial conditions related with (1.1) are given as

and the boundary conditions are assumed to be homogeneous Neumann boundary conditions, that is,

where $\partial \Omega$ is the boundary of $\Omega$ and $\textbf{n}$ is the outward normal of $\partial \Omega$. With this condition (1.3), the total mass

is conserved as temporal evolution. Besides the above, the model (1.1) has the following form of free energy:

From mathematical analysis, the following equation can be verified by direct calculation of Eq. (1.4),

Many early works show that the free energy (1.4) is decreasing over time and is mainly employed to demonstrate the existence of solutions for the chemotaxis system; see ^[3,4] and the references therein.

Chemotaxis refers to the directional movement, which includes toward or away from the higher concentrations of cells or microorganisms that are stimulated by chemical substances in the external environment along the gradient directions of the concentration for stimuli. Chemotaxis phenomena play a crucial role in numerous intricate biological evolutions, such as bacterial aggregation, angiogenesis, pattern formation, embryonic development and so on^[5]. From as early the 1950 to 1980, the chemotaxis phenomena have been extensively studied by many applied mathematicians and biologists, and a class of models of partial differential systems closely associated with taxis have been proposed; see ^{[1,2,6,7,8,9,10]}. Among them the most classical is the above-named system (1.1) (first proposed by Keller and Segel ^[1,2] in 1970), which simulated the aggregation phenomenon for Amoebae and Dictyostelium, as well as the traveling wave migration phenomenon for Escherichia coli in a capillary filled with nutrients.

Since the model was proposed, many researchers have systematically analyzed the properties of its solution, including its global existence, asymptotic profile ^[11], global boundedness ^[12,13], finite-time blow-up ^[14,15,16], etc. Particularly, if the initial mass $\iint_\Omega u(x, y)dxdy$ of the cells in the 2D case satisfies a critical threshold value, its solution will blow up in finite time ^{[7,8,9,14,15,16]}. This blow-up denotes a mathematical concept of the bacterial aggregation arising in real biotic environments ^[7,17,18]. However, it is arduous if we want to obtain the analytical solutions of the model due to the strong nonlinear characteristics of the model itself. Meanwhile, it is still difficult to better numerically capture the blow-up or spike solutions. Therefore, it is desperately needed to establish a high accuracy and more stable numerical method to investigate the properties of the solutions for the chemotaxis system given by (1.1)–(1.3). At the same time, the numerical investigation of the chemotaxis system also facilitates the theoretical exploration of its dynamic behavior.

In recent years, some numerical methods have been involved via the investigation of the chemotaxis systems ^{[3,4,16,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35]}. For example, Saito and Suzuki ^[19] proposed a conservative numerical scheme by using the FDM to solve a parabolic-elliptic coupling chemotaxis system. Meanwhile, in order to obtain a positivity-preserving scheme under total mass that is conservative, Saito ^[20] proposed a conservative scheme by using the FDM for the system given by (1.1)–(1.3) with nonlinear diffusion. Xiao et al. ^[21] derived a semi-implicit scheme by using a characteristic finite element method (FEM) to simulate the blow-up solutions of the chemotaxis system on surfaces, as well as pattern formulation and aggregation phenomenon of bacteria. And, their method has second-order accuracy in $L_2$-norm and $H_1$-norm errors. Epshteyn and Kurganov ^[22] introduced chemotaxis concentration gradient variables to rewrite the original Keller-Segel model given by (1.1)–(1.3), and they designed an internal penalty discontinuous Galerkin (DG) method for the rewritten system to border on the solution of the original system. Li et al. ^[23] used a local DG (i.e., LDG) method to modify that in ^[22], and they obtained the optimal convergence rates based on a specific finite element space before blow-up for the chemotaxis system. They also deduced a positivity-preserving $P^1$ LDG scheme and proved that it is $L^1$-stable. The LDG method was also used to solve the system given by (1.1)–(1.3) in ^[3], and the energy dissipation with the LDG discretization was proved. Sulman and Nguyen ^[24] proposed an adaptive moving mesh method by applying an implicit-explicit FEM to solve the system given by (1.1)–(1.3). The method has second-order accuracy in the spatiotemporal directions, and it is noteworthy that the obtained solutions are positive at all time steps if the initial values of the system (1.1) are positive. Qiu et al. ^[25] proposed a new scheme by using a interface-corrected direct DG method to solve this model (1.1), and their method satisfies the positivity-preserving requiremet without losing third-order accuracy. Based on the gradient flow structure, Shen and Xu ^[4] proposed a class of numerical schemes for solving the Keller-Segel model given by (1.1)–(1.3). Among them, the first-order accuracy scheme satisfies the mass conservation, bounded positivity, unique solvability and energy dissipation of the original differential equation, while the second-order accuracy scheme satisfies the first three properties. Chen et al. ^[16] analyzed the error of the numerical scheme proposed in ^[4], and they deduced the finite-time blow-up of non-radial numerical solutions under certain assumptions. Based on the generalized smoothed particle hydrodynamics meshless method, Dehghan and Abbaszadeh ^[26] established a second-order-accuracy numerical scheme for some chemotaxis models with the blow-up phenomena during tumor growth, and they numerically simulated the blow-up problems of the original chemotaxis model. Filbet ^[27] employed the finite volume method (FVM) to approximate the solution of the system given by (1.1)–(1.3), and they simulated the blow-up problems. Chertock and Kurganov ^[28] developed a center-upward scheme with second-order accuracy by using the FVM based on their previous work in ^[29] and the blow-up problems of the system (1.1)–(1.3), the pattern formation of bacteria were also simulated. Epshteyn ^[30] proposed an upwind-difference potentials scheme to solve the problems in complex geometries for the system given by (1.1)–(1.3). In addition, some other numerical methods ^{[31,32,33,34,35]}, such as the fractional step (or operator splitting) method ^[31,32], hybrid finite-volume-finite-difference methods ^[33], generalized FDM ^[35], etc., were also used to solve the chemotaxis system given by (1.1)–(1.3).

Although some numerical methods mentioned above for solving the system given by (1.1)–(1.3) can achieve high accuracy in the spatial direction, such as those in ^[25,30] (third-order accuracy) and ^[3,33] (fourth-order accuracy), most numerical methods have low accuracy, especially in the temporal direction; see ^{[4,16,19,20,21,22,23,24,26,27,28,29,31,32,34,35]}. Meanwhile, many high-accuracy numerical methods are non-compact in space, and the stability conditions are relatively harsh. In other words, if we want to employ these methods to solve real problems, we must take a small time step length to satisfy their stability conditions, which will expend expensive computational time. However, the high-order compact (HOC) difference scheme has attracted many researchers because of its strong advantages, such as fewer computational nodes, small computational errors, better numerical stability and non-complicated boundaries. Meanwhile, the backward differentiation formula with fourth-order accuracy (BDF-4), which is an $A$-stabilized method and appears first in ^[36], has been verified to be a feasible method to obtain high accuracy, and it has a relatively large stability condition range ^[37,38]. In addition, the major difficulty in solving the system comes from the nonlinearities, such as the chemotaxis term $\nabla\cdot(\chi u \nabla v)$. One issue is that the coupling form will increase the discrete difficulty when we want to obtain the high accuracy and satisfy mass conservative schemes. Another difficulty is positivity preservation since $u$ and $v$ in the system (1.1)-(1.3) have real biological significance in complex chemotaxis phenomenon and cannot have negative values. To acquire the high-accuracy, positivity-preserving numerical solutions for the system given by (1.1)–(1.3), in this work, our purpose was to derive a compact difference scheme to approximate the solutions of the original chemotaxis system, and the scheme has space-time fourth-order accuracy and is stable, as well as positivity-preserving.

The remainder of this paper is organized as follows. Some preliminary preparations with basic symbols, definitions and theorems are provided in Section 2. In Section 3, we deduce an HOC scheme for the system given by (1.1)–(1.3), and give the computational strategies for the initial time steps and the nonlinear terms. In Section 4, a time advancement algorithm combined with a multigrid method and a positivity-preserving algorithm are proposed. In Section 5, some numerical examples are employed to verify the accuracy, stability, positivity-preserving property, mass conservation and energy dissipation. The finitetime blow-up problems for the chemotaxis system given by (1.1)–(1.3) are simulated by using the proposed method. Finally the conclusion is provided in the end.

2.   Preliminaries

First, the domain $\{(x, y, t)|a\leqslant x, y\leqslant b, 0\leqslant t\leqslant T\}$ is divided by uniform meshes $N^2\times M$, $M, N\in \textbf{Z}^{+}$. Denote $h = (b-a)/N$ to represent the spatial step size, and $\tau = T/M$ stands for the temporal step length. Let $x_i = a+ih, y_j = a+jh$ and $t_n = n\tau$, $0\leqslant i, j\leqslant N$, $0\leqslant n \leqslant M$, and mark the mesh points $(x_i, y_j, t_n)$. Figure 1 shows the 2D spatial mesh point stencil.

Second, we define $\Omega_h = \{(x_i, y_j)|0\leqslant i, j\leqslant N\}$; its discrete boundary $\partial\Omega_h = \{(0, j), (N, j)|0\leqslant j\leqslant N\}\cup \{(i, 0), (i, N)|1\leqslant i\leqslant N-1\}$; let $\Omega_{\tau} = \{t_n|0\leqslant n\leqslant M\}$ and $\Omega_{h\tau} = \Omega_h\times \Omega_{\tau}$. For any mesh function $w\in\mathcal{W}_{h\tau} = \{w_{{i, j}}^n|0\leqslant i, j\leqslant N, 0\leqslant n\leqslant M\}$ defined on $\Omega_{h\tau}$, we have

Next, for simplicity, we let $\frac{\partial^pw}{\partial \varsigma^p}(x, y, t):\triangleq w_{\varsigma^p}(x, y, t)$, $p\in \textbf{Z}^{+}$, where $\varsigma$ represents $x$, $y$ or $t$. And, denote $\textbf{A} = [A_1, A_2, A_3, A_4, A_5] = [-25,48,-36,16,-3]$, $\textbf{B} = [B_1, B_2, B_3] = [1,10,1]$, $\textbf{S} = [S_1, S_2, S_3, S_4] = [-1, \frac{3}{4}, -\frac{1}{3}, \frac{1}{16}]$. To obtain a higher accuracy on the boundary described by Eq.(1.3), we employ the Taylor expansion with the Peano remainder to cope with Eq. (1.3). Thus, the following theorem holds.

Theorem 2.1. Denote $\mathcal{W}_h = \{w(x_i, y_j)|0\leqslant i, j \leqslant N\}$; for mapping $w:\Omega_h\rightarrow \mathcal{W}_h$, we have the following:

(1) If $w(x, y)\in C^{5, 0}([x_0, x_4]\times[y_0, y_N])$, then $w_x(x_0, y_j) = \frac{1}{12h}\sum\limits_{k = 1}^5A_kw(x_{k-1}, y_j)+O_{0, j}(h^4);$

(2) If $w(x, y)\in C^{5, 0}([x_{N-4}, x_N]\times[y_0, y_N])$, then $w_x(x_N, y_j) = -\frac{1}{12h}\sum\limits_{k = N-5}^{N-1}A_{N-k}w(x_{k+1}, y_j)+O_{N, j}(h^4);$

(3) If $w(x, y)\in C^{0, 5}([x_0, x_N]\times[y_0, y_4])$, then $w_y(x_i, y_0) = \frac{1}{12h}\sum\limits_{k = 1}^{5}A_{k}w(x_i, y_{k-1})+O_{i, 0}(h^4);$

(4) If $w(x, y)\in C^{0, 5}([x_0, x_N]\times[y_{N-4}, y_N])$, then $w_y(x_i, y_N) = -\frac{1}{12h}\sum\limits_{k = N-5}^{N-1}A_{N-k}w(x_i, y_{k+1})+O_{i, N}(h^4),$ where the local truncation errors are

Proof. First, we prove (1). We assume that $w(x, y)\in C^{k+1, 0}([x_{i-1}, x_{i+1}]\times[y_{0}, y_N])$, and, according to the Taylor expansion with the Peano remainder, we have

where $1\leqslant i\leqslant N-1$ and $0\leqslant j\leqslant N$. We suppose that $w(x, y)\in C^{5, 0}([x_{0}, x_{4}]\times[y_{0}, y_N])$ and expand $w(x_1, y_j)$, $w(x_2, y_j)$, $w(x_3, y_j)$ and $w(x_4, y_j)$ at $(x_0, y_j)$ by using Eq. (2.1) above, that is, we respectively take $k = 1, 2, 3, 4$, and have

Then, we perform the following the operation: Eq. (2.2)$_{k = 1}$+$\alpha\; \times$ Eq. (2.2)$_{k = 2}$+$\beta\; \times$ Eq. (2.2)$_{k = 3}$+$\gamma\; \times$ Eq. (2.2)$_{k = 4}$ for Eq. (2.2), and denote $\textbf{E} = [E_1, E_2, E_3, E_4] = [1, \alpha, \beta, \gamma]$, where $\alpha, \beta$ and $\gamma,$ are constants; then, we can obtain

Suppose that the coefficients of $w_{x^2}(x_0, y_j)$, $w_{x^3}(x_0, y_j)$ and $w_{x^4}(x_0, y_j)$ in Eq. (2.3) equal to $0$; we obtain

Substituting them into Eq. (2.3), denote $\textbf{A} = [-25,48,-36,16,-3]$ and $\textbf{S} = [-1, \frac{3}{4}, -\frac{1}{3}, \frac{1}{16}]$, and we have

where $O_{0, j}(h^4) = \frac{h^4}{6}\sum\limits_{k = 1}^{4}\int_0^1k^5S_kw_{x^5}(x_0+ksh, y_j)(1-s)^4ds$. Similarly, (2) and (3) are also easily obtained. The proof is complete.

According to Theorem 2.1, and by considering $w_x(x_0, y_j) = 0$, $w_x(x_N, y_j) = 0$, $w_y(x_i, y_0) = 0$ and $w_y(x_i, y_N) = 0$, we can easily derive the following fourth-order-accuracy boundary approximation formulas:

In addition, we suppose that $w(x, y, t)\in C^{6, 6, 5} ([x_{i-1}, x_{i+1}]\times[y_{j-1}, y_{j+1}]\times[t_n, t_{n-1}])$ for $1\leqslant i, j\leqslant N-1$, $1\leqslant n \leqslant M$, and we denote $\textbf{I} = [1,-24,81,-64]$ to derive the following truncation errors based on the Taylor expansion with the Peano remainder, that is,

Similarly, we can easily derive the expressions of $(O_{yy})_{i, j}^n(h^4)$, $(O_{y})_{i, j}^n(h^4)$, $(O_{x})_{i, j}^{n-1}(h^4)$, $(O_{y})_{i, j}^{n-1}(h^4)$, $(O_{xx})_{i, j}^{n-\frac{1}{2}}(h^4)$ and $(O_{yy})_{i, j}^{n-\frac{1}{2}}(h^4)$. To facilitate mathematical analysis below, the following definitions are given.

Definition 2.2. For any mesh function $\{w_{i, j}|0\leqslant i, j\leqslant N\}$, the average operators $\mathcal{A}$ and $\mathcal{B}$ are defined as

Definition 2.3. The maximum norm and $2$-norm errors are defined as

where $w(x_i, y_j, t_M)$ and $w_{{i, j}}^M$ stand for the exact and numerical solutions at the discrete mesh point $(x_i, y_j, t_M)$, respectively.

Definition 2.4. Denote $w_{\max} (t):\triangleq\max\limits_{(x, y)}w(x, y, t)$; the relatively maximum norm and $2$-norm errors in the temporal dimension are defined as

where $w_{\max}^* (t)$ represents the reference solution.

Definition 2.5. The convergence rate is defined as

where $\nu$ stands for $\infty$ or $2$ and $L_\nu (h_1)$ and $L_\nu (h_2)$ are corresponding $||\cdot||_\nu$ norm errors which are closely related to $h_1$ and $h_2$, respectively.

3.   Derivation for the compact difference scheme

In this part, an HOC scheme is derived to border on the solution of the Keller-Segel system given by (1.1)–(1.3). To facilitate the numerical analysis, we rewrite an equivalent form for the system given by (1.1)–(1.3), that is,

where

Equation (3.1) is a nonlinear advection-diffusion-reaction equation. Then, we define

on $\Omega_{h\tau}$, respectively, and employ the FDM to discretize Eqs. (3.1)–(3.4) for the interior points.

3.1.   Fourth-order compact scheme

Now, we focus on Eq. (3.1) at $(x_i, y_j, t_n)$ for $1\leqslant i, j\leqslant N-1$ and $0 < n \leqslant M$; we have

First, the following formulas are applied to discretize the second derivative terms on the right-hand side of Eq. (3.5), i.e.,

where $1\leqslant i, j\leqslant N-1$, and $\delta_x^2$, and $\delta_y^2$ are the central difference operators. Then, for the nonlinear advection terms $\textbf{F}_{x}(x_i, y_j, t_n)$ and $\textbf{G}_{y}(x_i, y_j, t_n)$ on the left-hand side of Eq. (3.5), the following Padé compact schemes ^[39] are applied to compute them, that is,

where $1\leqslant i, j\leqslant N-1$ and $0\leqslant n\leqslant M$. To be consistent with the accuracy of the interior points above, the following fourth-order boundary schemes ^[40] are employed, i.e.,

where

Next, we employ the BDF-4 ^[36] to cope with the temporal derivative on the left-hand side of Eq. (3.5), i.e., for $1\leqslant i, j\leqslant N-1$ and $4\leqslant n\leqslant M,$

Substituting Eqs. (3.6)–(3.12) into Eq. (3.5), using the definitions above, combining Eqs. (3.2)–(3.4) and applying Theorem 2.1 for $4\leqslant n\leqslant M$, we have

where $\textbf{Q}_{{i, j}}^{n} = \mathcal{AB}(O_t)_{i, j}^n(\tau^4) +\{[\mathcal{B}(O_{xx})_{i, j}^n+\mathcal{A}(O_{yy})_{i, j}^n]$+$\mathcal{AB}[(O_x)_{i, j}^n+(O_x)_{0, j}^n+(O_x)_{N, j}^n$+$(O_y)_{i, j}^n+(O_y)_{i, 0}^n+(O_y)_{i, N}^n]\}(h^4)$, $\textbf{Q}_{0, j}^{n} = O_{0, j}^n(h^4)$, $\textbf{Q}_{N, j}^{n} = O_{N, j}^n(h^4)$, $\textbf{Q}_{i, 0}^{n} = O_{i, 0}^n(h^4)$ and $\textbf{Q}_{i, N}^{n} = O_{i, N}^n(h^4)$. Then, there exist positive constants $c_1$, $c_2$, $c_3$, $c_4$ and $c_5$; it holds that

Omitting $\textbf{Q}_{i, j}^n$ in Eqs. (3.13)–(3.16) and replacing $\mathbb{U}$, $\mathbb{F}_x$, $\mathbb{G}_y$ and $\mathbb{R}$ with $\textbf{U}$, $\textbf{F}_x$, $\textbf{G}_y$ and $\textbf{R}$, respectively, the following HOC scheme for solving the chemotaxis system given by (1.1)–(1.3) for $4\leqslant n\leqslant M$ can be obtained:

Remark 3.1. For the function values containing the unknown time steps $\textbf{U}_{i, j}^{n}$ in Eqs. (3.17)–(3.20), the following method is used to compute them, where $k = 1, 2, 3, \cdots$, which stands for the nonlinear cyclic iterative parameter, $\textbf{U}^{n, (*)}$ denotes the approximate convergent solution at the $n$th time step and $\sigma$ represents a small amount. That is, for the $n$th time step, the following is performed:

● Approximate $\textbf{U}^{n, (k)}$ using $\textbf{U}^{n, (k-1)}$ and solve $v_x^{n, (k)}$ and $v_y^{n, (k)}$ using Eqs. (3.3), (3.4) and (3.7);

● Compute the chemotaxis terms $\textbf{F}_x^{n, (k)}$ and $\textbf{G}_y^{n, (k)}$ using Eqs. (3.7)–(3.11);

● Solve the linear system given by (3.17)–(3.20) for $\textbf{U}^{n, (k)}$;

● Set $\textbf{U}^{n, (*)}\leftarrow \textbf{U}^{n, (k)}$; if $\|\textbf{U}^{n, (k)}-\textbf{U}^{n, (k-1)}\| < \sigma$, then proceed to the $(n+1)$th time step.

Remark 3.2. Equations (3.17)–(3.20) demonstrate fully implicit compact scheme including five time steps, and its truncation error is $O(\tau^4+h^4)$, i.e., it has space-time fourth-order accuracy. It is noteworthy that the compactness here only refers to the spatial direction, since no more than nine mesh points are required for the 2D spatial mesh subdomain. However, it is non-compact in the temporal direction.

Remark 3.3. For the HOC scheme given by (3.17)–(3.20), besides the values of $\textbf{U}^{0}$ being known, the values of $\textbf{U}^{1}$, $\textbf{U}^{2}$ and $\textbf{U}^{3}$ are also needed; then, we can use the scheme given by (3.17)–(3.20) to compute the values of $\textbf{U}^{n}(n\geqslant 4)$ in turn. Thus, in the following part, we discuss the computational strategies for solving $\textbf{U}^{1}$, $\textbf{U}^{2}$ and $\textbf{U}^{3}$.

3.2.   Computation for $\textbf{U}^{1}$, $\textbf{U}^{2}$ and $\textbf{U}^{3}$

Next, we consider Eq. (3.1) at $(x_i, y_j, t_{n-\frac{1}{2}})$ for $1\leqslant i, j\leqslant N-1$ and $1\leqslant n \leqslant 3$; we have

First, we employ the Crank-Nicolson (C-N) method to cope with the temporal derivative term in Eq. (3.21), that is,

Then, we employ the following method for the remaining parts, i.e.,

where $\mathbf{\Theta}$ could be $\textbf{U}, \textbf{F}_x, \textbf{G}_y$ and $\textbf{R}$ in Eq. (3.21). Then, Eq. (3.21) becomes

We use Eq. (3.6) to discretize $\textbf{U}_{xx}(x_i, y_j, t_{n})$, $\textbf{U}_{xx}(x_i, y_j, t_{n-1})$, $\textbf{U}_{yy}(x_i, y_j, t_{n})$ and $\textbf{U}_{yy}(x_i, y_j, t_{n-1})$ in the spatial direction of Eq. (3.22), and we use Eqs. (3.7)–(3.11) to compute $\textbf{F}_{x}(x_i, y_j, t_n)$, $\textbf{F}_{x}(x_i, y_j, t_{n-1})$, $\textbf{G}_{y}(x_i, y_j, t_{n})$ and $\textbf{G}_{y}(x_i, y_j, t_{n-1})$ in the spatial direction. According to Eqs. (3.2)–(3.4) and Theorem 2.1, after rearrangement, the following form can be obtained:

where $\textbf{P}_{{i, j}}^{n} = \mathcal{AB}[(O_t)_{i, j}^{n-\frac{1}{2}} +(\tilde{O}_{xx})_{i, j}^{n-\frac{1}{2}}+(\tilde{O}_{yy})_{i, j}^{n-\frac{1}{2}} +(\tilde{O}_{\textbf{F}_x})_{i, j}^{n-\frac{1}{2}} +(\tilde{O}_{\textbf{G}_x})_{i, j}^{n-\frac{1}{2}} +(\tilde{O}_{\textbf{R}})_{i, j}^{n-\frac{1}{2}}](\tau^2) +\frac{1}{2}\{[\mathcal{B}((O_{xx})_{i, j}^n+(O_{xx})_{i, j}^{n-1}) +\mathcal{A}(O_{yy})_{i, j}^n+\mathcal{A}(O_{yy})_{i, j}^{n-1}] +\mathcal{AB}[(O_x)_{i, j}^n+(O_x)_{0, j}^n+(O_x)_{N, j}^n+(O_y)_{i, j}^n+(O_y)_{i, 0}^n+(O_y)_{i, N}^n +(O_x)_{i, j}^{n-1}+(O_x)_{0, j}^{n-1}+(O_x)_{N, j}^{n-1}+(O_y)_{i, j}^{n-1}+(O_y)_{i, 0}^{n-1}+(O_y)_{i, N}^{n-1}]\}(h^4)$. $\textbf{P}_{0, j}^{n} = O_{0, j}^n(h^4)$, $\textbf{P}_{N, j}^{n} = O_{N, j}^n(h^4)$, $\textbf{P}_{i, 0}^{n} = O_{i, 0}^n(h^4)$, and $\textbf{P}_{i, N}^{n} = O_{i, N}^n(h^4)$. Then, there exist positive constants $c_6$, $c_7$, $c_8$, $c_9$ and $c_{10}$ such that

Omitting $\textbf{P}_{i, j}^n$ in Eqs. (3.23)–(3.26) and replacing $\mathbb{U}$, $\mathbb{F}_x$, $\mathbb{G}_y$ and $\mathbb{R}$ with $\textbf{U}$, $\textbf{F}_x$, $\textbf{G}_y$ and $\textbf{R}$, respectively, we can get the following scheme for the initial time steps to solve the chemotaxis system given by (1.1)–(1.3) as follows:

Remark 3.4. The same method in the previous section can be used to compute $(\textbf{F}_{x})_{i, j}^{n}$, $(\textbf{F}_{x})_{i, j}^{n-1}$, $(\textbf{G}_{y})_{i, j}^{n}$, $(\textbf{G}_{y})_{i, j}^{n-1}$, $\textbf{R}_{i, j}^{n}$ and $\textbf{R}_{i, j}^{n-1}$ in Eq. (3.27), as described in Remark 3.1.

Remark 3.5. Equations (3.27)–(3.30) constitute a two-level implicit scheme, and its truncation error is $O(\tau^2+h^4)$.

Remark 3.6. Since the C-N method is used to discretize the time derivative term, the scheme given by (3.27)–(3.30) for the initial time steps is unconditionally stable. On the other hand, the formula (3.12) is a $k$-step method for time integration. In this work, only the fourth-order BDF scheme ($k = 4$) is used, which is widely employed to solve stiff differential equations ^[37]. Therefore, according to ^[36,37,38], the HOC scheme given by (3.17)–(3.20) is $A$-stable.

3.3.   Extrapolation for $\textbf{U}^{1}$, $\textbf{U}^{2}$ and $\textbf{U}^{3}$ in the temporal dimension

As described in Subsections 3.1 and 3.2 above, we use the two-level implicit scheme given by (3.27)–(3.30) to compute $\textbf{U}_{i, j}^{1}$, $\textbf{U}_{i, j}^{2}$ and $\textbf{U}_{i, j}^{3}$ based on a known $\textbf{U}_{i, j}^{0}$. The Richardson extrapolation technique is applied to enhance the accuracy of the scheme given by (3.27)–(3.30) in time and ensure that it is consistent with that of the scheme given by (3.17)–(3.20); the following extrapolation formula is used, that is,

Here, $\widehat{\textbf{U}}_{i, j}^n(h, \tau)$ stands for the extrapolated solution at the $n$th time step. $\textbf{U}_{i, j}^n(h, \tau)$ represents the computed solution obtained by using the scheme given by (3.27)–(3.30) with the temporal step length $\tau$, whereas $\textbf{U}_{i, j}^{2n}(h, \tau/2)$ denotes the computed solution obtained by using the scheme given by (3.27)–(3.30) with $\tau/2$.

Remark 3.7. Equation (3.31) can extrapolate the accuracy in the temporal direction from second to fourth order.

4.   Numerical algorithm

In this part, based on the compact difference schemes proposed above, we will establish the corresponding numerical algorithm. On one hand, because the slow convergence speed of the classical iterative method (for instance, Gauss–Seidel, Jacobi or successive over-relaxation iterations) used to solve the algebraic systems arise out of the fully implicit scheme at each time step, we aim to employ a multigrid method to accelerate the convergence speed. On the other hand, we want to establish a positivity-preserving algorithm to obtain the non-negative solutions of $u$ and $v$ without losing the fourth-order accuracy. Finally, we structure a time advancement algorithm based on these algorithms above to solve the system given by (1.1)–(1.3).

4.1.   Multigrid algorithm

Since the idea of the multigrid method was proposed in the 1930s, until Professor Brandt published his pioneering work ^[41] in 1977, this method was increasingly applied to solve engineering and technical problems. Up to now, the multigrid method has been widely applied in numerous disciplines and engineering fields, such as computational fluid dynamics and computational biology. The multigrid method has been theoretically proved to be a first-rank numerical computational method for linear elliptic problems ^[41,42]. Its convergence rate is independent of the mesh scale, and the computational speed of the existing computational program using the classical iterative method can be increased by employing the multigrid method with $1$ to $2$ orders of magnitude, which is especially suitable for application in large-scale engineering numerical computational problems. Nowadays, the multigrid method is also used more and more to solve non-elliptic problems to speed up the convergence of each time step and satisfy the needs of practical problems ^[43,44].

The multigrid method is implemented by using a cycling algorithm. A multigrid cycle $V$ includes three elements: relaxation, restriction and interpolation operators. We use alternating direction line Gauss–Seidel (ALGS) ^[45] iteration for the relaxation operator, as well as half-weighted and fully weighted restriction operators and a bilinear interpolation operator ^[42]. As the schemes given by (3.27)–(3.30) and (3.17)–(3.20) are both nonlinear, we employ a multigrid full approximation scheme (FAS) ^[42,44] to accommodate nonlinearities. For simplicity, we formally denote the algebraic equations arising from Eqs. (3.27)–(3.30) and (3.17)–(3.20) at each time step as

The multigrid FAS algorithm is implemented recursively. Here, we only give the following two-level FAS $V(\nu_1, \nu_2)$ cycle algorithm based on Eq. (4.1):

Remark 4.1. Here, $\mathcal{L}^{h}$ and $\mathcal{L}^{2h}$ are the difference operators; $\mathcal{I}_h^{2h}$ (half-weighted) and $\overline{\mathcal{I}}_h^{2h}$ (fully weighted) are the restriction operators, which are used to restrict the approximate value and residual from the fine mesh level $h$ to coarse mesh level $2h$; and, $\mathcal{I}_{2h}^{h}$ is the interpolation operator, which is used to transfer the error correction from the coarse to fine mesh levels. $Res^h$ is the residual at the fine mesh level, and $Res^{2h}$ is the residual at the adjacent coarse mesh level; $\nu_1$ and $\nu_2$ represent the pre-smooth and post-smooth numbers, respectively.

Remark 4.2. In this work, we employ fully multi-level algorithms for all calculations. For instance, if the finest grid is $256\times256$, an eight-level algorithm is used. Under such conditions, the coarsest grid level only has $2\times2$ grids.

4.2.   Positivity-preserving algorithm

When the schemes given by (3.27)–(3.30) and (3.17)–(3.20) are employed to approximate the solution of the chemotaxis system, it may generate negative values where blow-up occurs. To obtain the non-negative values of $u$ and $v$ for all time levels, we define the following average solution for each time step:

where $\Delta x_i = x_{i+1}-x_{i-1}$, $\Delta y_j = y_{j+1}-y_{j-1}$, $1\leqslant i, j\leqslant N-1$ and $\Omega_{i, j} = \{[x_{i-1}, x_{i+1}]\times[y_{j-1}, y_{j+1}]|1\leqslant i, j\leqslant N-1\}$. Similar to the numerical integration formula of the one-variable function (i.e., Simpson formula), we structure the numerical integration of the two-variable function $\textbf{U}(x, y)$ on $\Omega_{i, j}$ for each time step, i.e., for $1\leqslant i, j\leqslant N-1$,

Substituting Eq. (4.3) into Eq. (4.2) at each time step, we obtain the following average solutions $\overline{\textbf{U}}$, that is, for $1\leqslant i, j\leqslant N-1$ and $0 < n\leqslant M$, we have

Next, a positivity-preserving limiter (PPL) ^[25,46] is employed to eliminate the negative values of $\textbf{U}_{i, j}^n$, $0\leqslant i, j\leqslant N-1$, $0 < n\leqslant M$, that is,

where $\varpi_{i, j}^n = \min\{\textbf{U}_{i, j}^n, \textbf{U}_{i-1, j}^n, \textbf{U}_{i+1, j}^n, \textbf{U}_{i, j+1}^n, \textbf{U}_{i, j-1}^n\}, \; \; 1\leqslant i, j\leqslant N-1, \; 0 < n\leqslant M.$ For $1\leqslant i, j\leqslant N-1$ and $0 < n\leqslant M$, we can get the positivity-preserving solution $\widetilde{\textbf{U}}_{i, j}^n\geqslant 0$ without losing the fourth-order accuracy obtained by the proposed scheme. Finally, we establish the following algorithm for the $n$th time step:

4.3.   Time advancement algorithm

On the basis of the derivation process described in Section 3 and the proposed algorithms in Subsections 4.1 and 4.2, we employ the difference schemes given by (3.27)–(3.30) and (3.17)–(3.20) to approximate the solution of the original system given by (1.1)–(1.3), which mainly includes the following three steps: first, provide the initial value of $\textbf{U}_{i, j}^{0}$; second, we employ Eqs. (3.7)–(3.11) and (3.27)–(3.31) to compute the values of $\textbf{U}_{i, j}^{1}$, $\textbf{U}_{i, j}^{2}$ and $\textbf{U}_{i, j}^{3}$; finally, we employ Eqs. (3.7)–(3.11) and (3.17)–(3.20) to compute the value of $\textbf{U}_{i, j}^{n}$ $(n\geqslant 4)$. As (3.27)–(3.30) and (3.17)–(3.20) involve nonlinear parts, we establish the following algorithm with a nonlinear iteration tactic to solve them.

5.   Numerical simulation experiments

In this part, we give several numerical experiments to test the various properties of the proposed scheme and algorithm when solving the system given by (1.1)–(1.3). First, in Subsection 5.1 below, the accuracy and stability of the proposed scheme and algorithm in the absence of PPL and in the presence of a PPL when solving the chemotaxis system are tested. Second, in Subsection 5.2 below, we simulate the blow-up phenomena for the system given by (1.1)–(1.3), compare the obtained results from the proposed scheme and algorithm before and after using the positivity-preserving algorithm and verify the mass conservation and energy dissipation of the proposed method.

5.1.   Performance testing of the difference scheme

5.1.1.   Accuracy and stability experiments in the absence of a PPL

We consider the following general type of Eq. (1.1) with source terms to test the accuracy, that is,

We use the periodic boundary conditions, and its analytical solution is $u(x, y, t) = v(x, y, t) = (2+\sin(x+y))e^{-2t}.$ By observing the expression of this analytical solution, we can see that the global solutions of this problem are positive values in the entire physical domain. Therefore, in this computation, we will use the difference scheme proposed in this work to approximate the solution for this problem in the absence of a PPL.

In Table 1, we take $\tau = 0.1h$ and test the convergence of the HOC scheme at the final time $T = 0.1$ for Problem 5.1.1. From this table, the fourth-order accuracy is obtained by using the proposed scheme. The obtained results for the LDG method in ^[3] (the time step size $\tau = 0.001h^2$ is used to solve this problem) represent third-order accuracy, which shows the advantage of our method, that is, we can obtain higher accuracy with a larger time step. Table 2 displays the convergence of the HOC scheme at $T = 1$ and $T = 5$ for Problem 5.1.1, where $\tau = 0.5h$. Table 3 shows the $||\cdot||_2$ errors of $u(x, y, t)$ and $v(x, y, t)$ for different step ratios $\lambda = \tau/h$ for Problem 5.1.1, where $T = 1$ and $\tau = \lambda h$. From these tables, we obtain that the proposed method still converges at a fourth-order rate when computing long times and large step ratios, which indicates that our method has better stability.

5.1.2.   Accuracy and stability experiments in the presence of a PPL

Next, we consider the following general type of Eq. (1.1) with two additional fluxes^[24]:

where $\chi = 1$, and the homogeneous Neumann boundary condition (1.3) is applied.

Case 1: The additional fluxes are $r_1(x, y, t) = -2 \exp(-2t)(\cos^2(x)+\cos(x)\cos(y)+\cos^2(y)-1)$, $r_2(x, y, t) = 0$, with the initial condition $u(x, y, 0) = v(x, y, 0) = \cos(x)+\cos(y), 0\leqslant x, y\leqslant 2\pi$, and the analytical solution^[23,24] ${u(x, y, t) = v(x, y, t) = \exp(-t)(\cos(x)+\cos(y))}.$

Case 2: The additional fluxes are $r_1(x, y, t) = -2e^{-4t}\cos(2x+2y)$, $r_2(x, y, t) = 0$, with the initial condition $u(x, y, 0) = v(x, y, 0) = \cos(x+y), 0\leqslant x, y\leqslant 2\pi$, and the analytical solution $u(x, y, t) = v(x, y, t) = e^{-2t}\cos(x+y)$.

Table 4 shows the convergence of the HOC scheme when we take the time step size $\tau = 0.001$ at $T = 0.1$ when using the proposed method to solve Problem 5.1.2 Case 1: the obtained results are compared with those in ^[24]. According to Table 4, the results obtained via the proposed method in the absence of a PPL can converge at the fourth-order rate. After using the positivity-preserving algorithm, although the errors obtained via the proposed scheme in the presence of a PPL decreases slightly, it still converges at the fourth-order rate and is better than the results in ^[24], as it is only second-order accuracy in ^[24]. In Tables 5 and 6, the results computed by using the proposed method for Problem 5.1.2 are given respectively. And, the convergence of the HOC scheme in the presence and absence of a PPL are compared when $\tau = 0.5h$, $T = 1$ and $T = 5$. By observing these tables, it can be found that the computed results for the proposed method can achieve the fourth-order accuracy in both cases, including with and without a PPL. Finally, we take the spatial meshes $N = 64,128$ and $256$, respectively, at the final time $T = 1$. By using the proposed method to solve Problem 5.1.2 Case 1, we obtain the $||\cdot||_2$ norm errors computed with and without a PPL for different step ratios $\lambda$. The computed results are listed in Table 7. On the basis of this table, with the increase of the step ratio $\lambda$, the proposed method can converge well regardless of PPL existence, which also reflects that the proposed HOC difference scheme has better stability.

Table 8 lists the numerical convergence of HOC scheme Problem 5.1.2 Case 2 in the absence and presence of a PPL, where $T = 0.2$ and $\tau = 0.1h$, and it compares the computed results with those obtained via the method in ^[25]. The results in ^[25] have third-order accuracy before and after using a PPL, while the results computed via the method in this work converge at the fourth-order rate for both cases, including without and with a PPL, and the $||\cdot||_{\infty}$, $||\cdot||_{2}$ errors are better than those in ^[25]. This also reflects the superiority of the HOC difference scheme proposed in this work.

5.2.   Blow up in finite time, mass conservation and non-negativity experiments

In the 2D case of the system given by (1.1)–(1.3), when the initial mass satisfies certain critical values, the solutions will blow up in a finite time. To capture the blow-up time of the solutions, referring to ^[47], the following adaptive technology is used to obtain the optimal time step sizes, i.e.,

The computation is terminated if $||\mathbb{U}||_{\infty}\geqslant 10^{\nu}$, where $\nu$ denotes the maximum order of magnitude of cell density when the actual problem experiences a blow up. And, $t_n = \sum\limits_{m = 0}^{n-1}\tau_m$ is regarded as the approximation of the blow-up time $T$.

5.2.1.   Blow up at the center of the rectangular region

In a rectangular region $[-0.5, 0.5]\times[-0.5, 0.5]$, we first take $\chi = d = 1$ and test the initial-boundary value problem (IBVP) ^[23,24,25] for the chemotaxis system given by (1.1)–(1.3) with the homogeneous Neumann boundary conditions given by (1.3), where its initial conditions are as follow:

The solution of this problem will blow up in finite time at the center of this rectangular region.

(1) Finite-time blow up. Figure 2 plots the numerical solutions of cell density $u(x, y, t)$ as obtained by using the proposed method in the presence of a PPL, to solve the IBVP 5.2.1 at $t = 0$, $t = 5\times10^{-5}$, $t = 10^{-4}$ and $t = 1.2\times10^{-4}$, where the initial time step size $\tau_0 = 5\times10^{-6}$ and the fixed spatial mesh numbers are $100\times 100$. In ^[22], the blow-up time $t$ is approximately equal to $1.21\times10^{-4}$, which is verified by the authors. Based on observation, we can expect to see that Problem 5.2.1 displays blow-up phenomena around $t\thickapprox 1.2\times10^{-4}$ at the center of the rectangular region $[-0.5, 0.5]\times[-0.5, 0.5]$. It shows that the maximum peak values of the cell density $u(x, y, t)$ gradually increase over time and present a very sharp aiguille structure at the central zone. We find that $u(x, y, t)$ is strictly positive during time evolution by using the proposed method in the presence of a PPL. Figure 3 displays the projections of $u(x, y, t)$ on the $xou$ plane for Problem 5.2.1 at two pre-blow-up times $t = 10^{-4}$ and $t = 1.2\times10^{-4}$, with $\tau_0 = 5\times10^{-6}$ and $N = 100$. They intuitively show that the proposed method in this work can achieve positivity-preserving capability.

(2) Non-negativity. To highlight the performance of positivity-preserving capability for our method in the presence of a PPL, we employ the proposed scheme in the absence and presence of a PPL to solve the IBVP 5.2.1 in two different space grids $N = 50$ and $N = 100$ at pre-blow-up times $t = 10^{-4}$ and $t = 1.2\times10^{-4}$, where $\tau_0 = 5\times10^{-6}$. In Figure 4, we can see that $u(x, y, t)$ has negative values in the absence of a PPL, and that all negative values are eliminated after the PPL is used. At the same time, the negative values decrease with the increase of spatial grid numbers, and the solutions become smoother. However, for the same number of grids, the oscillation of the blow-up area increases with time.

(3) Influence of chemotaxis sensitivity coefficient $\chi$. We employ the proposed method to compute the blow-up solution in a finite time for this IBVP 5.2.1 under different chemotaxis sensitivity coefficients $\chi = 2$ and $\chi = 3$ at the same pre-blow-up time $t = 5\times10^{-5}$. The results are shown in Figure 5, where $\tau_0 = 10^{-6}$ and $N = 100$. We find that, for different values of $\chi$, the cell density $u(x, y, t)$ achieves negative values in the absence of a PPL. After using the PPL, the algorithm effectively eliminates the negative values. At the same time, the maximum peak value when $\chi = 3$ is one order of magnitude higher than that when $\chi = 2$. It can be seen that the peak value is gradually increased over $\chi$. Therefore, the sensitivity intensity $\chi$ also affects the value of $u(x, y, t)$ and the blow-up degree of Problem 5.2.1.

(4) Chemoattractant concentration. In Figure 6, we plot the numerical solutions of chemoattractant concentration $v(x, y, t)$ for Problem 5.2.1 when $t = 10^{-4}$, $\tau_0 = 5\times10^{-6}$ and $N = 100$, and we have compared the computed results in the absence and presence of a PPL. Observing Figures 4b and 6, we find that, although we applied the same parameters in the absence of a PPL, $u(x, y, t)$ achieves negative values, but $v(x, y, t)$ does not have negative values, and that the numerical solutions in the presence of a PPL is highly consistent with those in the absence of a PPL. In summary, our method can effectively eliminate these negative values in this computation. It is consistent with the results obtained in the absence of a PPL if the obtained solutions do not have negative values. Meanwhile, Problem 5.2.1 always experiences a finite-time blow-up at the center of the rectangular region under this initial condition.

(5) Accuracy and mass conservation. Since the analytical solution for this IBVP 5.2.1 is very difficult to obtain, to verify the accuracy and mass conservation, we take the fine mesh $800\times800$ as the reference solution. The convergence of the HOC scheme in the presence of a PPL are tested by using the relative errors between the maximum value $\textbf{U}_{\max} (t)$ of $u(x, y, t)$ and the reference solutions $\textbf{U}_{\max}^*(t)$; the results are shown in Table 9. We find that $u(x, y, t)$ can converge at a fourth-order relative error convergence rate. Meanwhile, we plot $u_{max}(t)$, which is the time evolution of the maximum value of $u(x, y, t)$, on the left of Figure 7, and the numerical energy with time on the right of Figure 7, using the proposed scheme in the presence of a PPL to solve this IBVP 5.2.1 with PPL at $\tau_0 = 5\times 10^{-6}$, $N = 50, 80$, $100$. We find that the energy is strictly positive and monotonically decreasing with the evolution of time, which verifies the energy dissipation nature of the original problem. In addition, in Table 10, we give the discrete masses of $u(x, y, t)$ and $v(x, y, t)$ as obtained by using the proposed scheme in the absence and presence of a PPL to solve this IBVP 5.2.1 at different times $T$ with $\tau_0 = 5\times10^{-6}, N = 100$. According to this table, we can find that the proposed scheme in the absence of a PPL well verifies the mass conservation. After using the PPL, it is still conserved before blow up, and the mass conservation is slightly affected near the blow-up time.

5.2.2.   Blow up at the corner of the rectangular region

In the section, we consider the following IBVP ^[24] for the chemotaxis system given by (1.1)–(1.3) with the boundary conditions given by (1.3) for the rectangular region $[-0.5, 0.5]\times[-0.5, 0.5]$. Its initial condition is given as follows:

In the computation, we take $\chi = d = 1$ and employ the proposed scheme in the presence of a PPL to simulate the temporal evolution of $u(x, y, t)$ and $v(x, y, t)$. The solution for this Problem 5.2.2 will blow up in finite time ^[8,9,24] at the corner of the rectangular region $[-0.5, 0.5]\times[-0.5, 0.5]$. Next, we take the initial time step size $\tau_0 = 10^{-3}$ and fix a mesh of $200\times 200$ to compute them.

Figures 8 and 9 respectively display the evolutionary processes of $u(x, y, t)$ and $v(x, y, t)$ over time $t$. It can be seen in these figures that the maximum value of $u(x, y, t)$ gradually moves toward the corner of the rectangular region $[-0.5, 0.5]\times[-0.5, 0.5]$ in the form of a Gaussian profile. When it is close to the corner of the rectangular region, the value of $u(x, y, t)$ increases rapidly, finally blowing up in a finite time. At the same time, $v(x, y, t)$ also increases rapidly at the corner of the rectangular region.

6.   Conclusions

In this work, first, a fourth-order compact difference scheme for solving a 2D Keller-Segel model was derived and the computing strategies for the initial time steps and the nonlinear chemotaxis terms has been presented. Then, a multigrid algorithm was established to improve the computational efficiency and the binary function integration method has been used to establish the positivity-preserving algorithm. Third, the accuracy and reliability of the proposed method in the absence and presence of a PPL have been respectively verified by several numerical examples with flux source terms. And, we simulated the blow-up phenomena for the chemotaxis system in the center and corner of the rectangular region, respectively. Finally, the energy dissipation and mass conservation were also verified by numerical experiments. In summary, compared with the classical high-order numerical methods for solving the chemotaxis system in the literature, this proposed method has the following advantages:

● The proposed HOC scheme is compact in the spatial directions and fully implicit, as no more than nine mesh points are required for the 2D spatial mesh subdomain.

● The truncation error of the proposed HOC scheme is $O(\tau^4+h^4)$, i.e., it has space-time fourth-order accuracy. Thus, we can take large time step sizes to approximate the solution of this chemotaxis system.

● The proposed numerical algorithm can effectively filter the negative values in the solution process to ensure that the values of cell density are not negative at all time steps, while maintaining the original fourth-order accuracy without loss.

● The computed results show that the proposed method is more accurate than most such schemes reported in the literature.

Nevertheless, our method also has some disadvantages. First, the new scheme has five layers in time. When it is used to compute the actual problem, three start-up time steps are required, which is slightly more troublesome for programming. Second, the proposed positivity-preserving algorithm would slightly affect its mass conservation when the actual problem blows up. The above deficiencies are also the driving force for our next work. At the same time, it is also possible to extend our approach to solve more complex chemotaxis and haptotaxis systems, which will be incorporated into our next work.

Acknowledgments

Support of the study was received from the National Natural Science Foundation of China (12161067, 12001015, 12261067), National Natural Science Foundation of Ningxia, China (2022AAC02023) and National Youth Top-Notch Talent Support Program of Ningxia, China.

Conflict of interest

All authors assert no conflict of interest.

Data availability statement

The data that support the findings of this study are available from the authors upon reasonable request.