High-order explicit numerical methods for the one- and two-dimensional haptotaxis problems of tumor invasion

Abstract

High-precision numerical methods are developed for biomathematical models that describe various biotaxis problems of tumor invasion, facilitating an in-depth exploration of the underlying evolutionary mechanism of tumor invasion. In this paper, we construct a high-order explicit numerical method for solving haptotaxis models of tumor invasion. To incorporate the specific characteristics of the initial and no-flux boundary conditions for the haptotaxis models, we employ a high-order compact finite difference method to discretize the spatial derivatives, thereby obtaining a series of semi-discrete ordinary differential systems to approximate the solutions of these models. The strong stability-preserving (SSP) Runge-Kutta method is utilized in the semi-discrete ordinary differential systems to obtain the third-order accuracy in time. A high-order numerical integration method is used to design a positivity-preserving algorithm with spatially fourth-order accuracy guaranteed. The designed explicit method not only has high accuracy and good stability, but also has obvious positivity-preserving effect, which is more accurate and reliable than the numerical methods for this problem in the existing literature. To validate the performances of the proposed methods, several numerical examples are provided, and all sorts of complicated biotaxis dynamics for the haptotaxis models are simulated.

Introduction

Malignant tumor, commonly known as Cancer, pose a serious threat to human life and health worldwide. Especially in recent years, the global number of new cancer cases has been increasing annually [1]. Therefore, it is of far-reaching significance for the treatment and prevention of the tumor diseases to gain a deeper understanding of tumor growth and their underlying evolutionary mechanisms. The main biological characteristic of malignant tumor is invasion and metastasis. Tumor invasion is a dynamic process involving protein synthesis and decomposition, whereby tumor cells (TCs) secrete matrix degrading enzyme (MDE) to degrade the extracellular matrix (ECM) around the tumor and provide space for TCs to metastasize to distant tissues [2]. Tumor invasion almost exists whole procedure of tumor growth and is also the key to cancer metastasis.

Over the past thirty years, numerous biomathematical models [2–7] have been established to describe the complex process of tumor invasion. In this paper, we focus on the classical haptotaxis model for tumor invasion proposed by Anderson et al. [6] in 2000, which considers the mechanism of random motion and haptotaxis directional motion of tumor invasion. Haptotaxis refers to the directed movement of tumor cells along the gradient of non-diffusible ECM density. The haptotaxis model describes the interaction for the density $u(\mathbf{x},t)$ of TCs, the density $v(\mathbf{x},t)$ of ECM and the concentration $w(\mathbf{x},t)$ of MDE in the tumor microenvironment. Its general form is

1

$\{\begin{array}{ll}{u}_t=D_1\mathrm{△}u-\mathrm{\nabla }\cdot (\chi u\mathrm{\nabla }v)+F_1(u,v),& (\mathbf{x},t)\in \mathrm{\Omega }\times (0,T],\\ v_t=F_2(u,v,w),& (\mathbf{x},t)\in \mathrm{\Omega }\times (0,T],\\ w_t=D_2\mathrm{△}w+F_3(u,v,w),& (\mathbf{x},t)\in \mathrm{\Omega }\times (0,T],\end{array}$ with the initial data

2

$[u(\mathbf{x},0),v(\mathbf{x},0),w(\mathbf{x},0)]=[u_0(\mathbf{x}),v_0(\mathbf{x}),w_0(\mathbf{x})],\mathbf{x}\in \mathrm{\Omega },$ and the no-flux boundary

3

$D_1\mathrm{\nabla }u\cdot \mathbf{n}-\chi u\mathrm{\nabla }v\cdot \mathbf{n}=0,D_2\mathrm{\nabla }w\cdot \mathbf{n}=0,(\mathbf{x},t)\in \partial \mathrm{\Omega }\times (0,T],$ in which, $\mathrm{\Omega }\subset {\mathbf{R}}^d(d⩾1)$ is a smoothly convex bounded region, and ∂Ω is its boundary. $D_1$ and $D_2$ stand for the diffusion coefficients of TC and MDE, respectively. $-\mathrm{\nabla }\cdot (\chi u\mathrm{\nabla }v)$ denotes diffusion by haptotaxis. $\chi \in {\mathbf{R}}^{+}$ is the haptotaxis sensitivity function. n is the outward normal vector defined on ∂Ω, $T\in {\mathbf{R}}^{+}$ is the physical terminal time. In addition, $F_1$, $F_2$ and $F_3$ have the forms as in Table 1.

Table 1. Various expressions of $F_i(i=1,2,3)$

Situations	$F_1$	$F_2$	$F_3$
Case 1 [6, 8–12]	0	−ηwv	αu − βw
Case 2 [6, 13, 14]	${\mu }_{1}u(1-u-v)$	−ηwv	αuv − βw
Case 3 [6, 8, 15]	${\mu }_{1}u(1-u-v)$	−ηwv	αu(1 − w)−βw
Case 4 [6, 12, 16]	${\mu }_{1}u(1-u-v)$	$-\eta wv+{\mu }_{2}v(1-u-v)$	αu − βw
Case 5 [6, 12, 16–20]	${\mu }_{1}u(1-u-v)$	$-\eta wv+{\mu }_{2}v(1-u-v)$	αu(1 − w)−βw

In Table 1, α, β, η, ${\mu }_1$, and ${\mu }_2$ are positive constants. $F_1(u,v)$ stands for the proliferation by the interaction between TC and ECM which abides by the logistic growth law of the competition for space. $F_2(u,v,w)$ indicates that the ECM is degraded by MDE, and it satisfies logistic remodeling growth. $F_3(u,v,w)$ represents the natural degradation of MDE and the proliferation by the interaction between TC, ECM and MDE. Moreover, $u_0(\mathbf{x})$, $v_0(\mathbf{x})$ and $w_0(\mathbf{x})$ are known smooth functions.

Since many mathematical models for tumor invasion were proposed, more and more researchers were attracted to study them in depth to reveal various biotaxis phenomena in which TCs invade body tissues. These models are a class of cross-coupled nonlinear differential dynamical systems. They have complex structures and strong nonlinearity such that it is difficult to obtain their exact solutions. Therefore, in recent years, many researchers tend to develop various numerical methods for solving them, such as finite volume methods (FVMs) [5, 21], finite element methods (FEMs) [8, 13, 17, 19, 20, 22], meshless methods [9], finite difference methods (FDMs) [11, 12, 23, 24] and spectral methods [25], etc. For instance, Kolbe et al. [21] used the FVM to establish a numerical method with second-order accuracy in space and conditional stability for solving one-dimensional (1D) and two-dimensional (2D) chemotaxis-haptotaxis models for tumor invasion, as well as to design a space-time mesh refinement adaptive algorithm with a positivity-preserving property. Niño-Celis et al. [13] constructed two numerical methods to solve the 2D haptotaxis model for tumor invasion by using the FEM, their truncation error terms are $O(\tau +h^2)$. Dehghan and Narimani [9] used Galerkin meshless and moving least squares methods to derive two numerical schemes to solve the 2D haptotaxis model for tumor invasion, and their truncation error terms are also $O(\tau +h^2)$. Kolev et al. [11] used the FDM to derive a numerical method for solving the 1D haptotaxis model, but its truncation error term is $O(\tau +h^2)$, and is conditionally stable. Recently, Kolev et al. [12] improved the scheme in [11] by using the non-standard FDM, and obtained an unconditionally positivity-preserving difference scheme for solving the 1D haptotaxis model.

With an in-depth understanding of the tumor invasion mechanisms, high-precision and high-resolution have become a decisive factor in the development of technology for improving the effectiveness and reliability of numerical simulation. However, most of the above numerical methods are of low-precision [8, 9, 11–13, 17, 19, 20, 23]. For the haptotaxis model (1)-(3), there are almost no high-precision numerical methods, only first-order accuracy in time [8, 9, 11–13, 20, 23], as well as first-order [8] or second-order accuracy [9, 11–13, 17, 20, 21, 23, 24] in space. It is well known that low-precision methods on coarse grids often yield inaccurate results. To achieve the accuracy required for practical applications, the computational grid must be sufficiently refined, which inevitably increases computational cost. While high-precision methods provide obvious advantages like smaller errors and higher resolution. On the other hand, the unknown functions in model (1)-(3), such as density and concentration, must remain non-negative due to their biomedical significance.

Researchers [12, 13] have tried to construct positivity-preserving numerical schemes, but it is quite challenging. For instance, the method proposed in [12] was proven to preserve non-negativity for ECM density and MDE concentration, but not for TC density. Meanwhile, the authors in [21, 24] use the monotone central limiter minmod function to construct a kind of positivity-preserving algorithm, but the resolution of the method is low and only has second-order accuracy. Usually, to match the proposed methods, it is necessary to maintain the original high-precision without loss while ensuring the non-negative numerical solution. In addition, most of the existing difference schemes [11, 12, 23, 24] are non-compact, meaning they need to increase the number of mesh points in their template to improve the computational accuracy. This leads to the increase of the bandwidth in the algebraic equations generated by the schemes, and the use of points outside the boundary, posing significant difficulties for computing the mesh points in the vicinity of the boundary. Conversely, high-order compact FDM [26, 27] has the advantages of small discrete subdomain, high precision and easy boundary processing. It can achieve high accuracy with relatively fewer grid points, and has been widely used in the numerical solutions of various partial differential equations [27–30].

For this purpose, we intend to construct a kind of high-precision, high-resolution and positivity-preserving numerical method for the 1D and 2D haptotaxis models of tumor invasion in this paper. The design strategy leverages the fourth-order compact difference method, the advantages of the boundary conditions of the original problem (1)-(3), the strong stability-preserving (SSP) Runge-Kutta method [31, 32], and the high-order numerical integration idea. The highlights of this work are as follows:

1. The proposed method for solving the haptotaxis problems (1)-(3) of tumor invasion achieves fourth-order accuracy in space and third-order accuracy in time.

2. A class of positivity-preserving algorithms with spatially fourth-order accuracy guaranteed for 1D and 2D haptotaxis problems (1)-(3) are established, respectively.

3. The dynamic behaviors of the solutions to haptotaxis models, including various forms of TCs proliferation terms, ECM degradation and remodeling, MDE attenuation and proliferation are numerically simulated.

The layout of the rest of this paper is as follows: In Sec. 2, we will provide a fourth-order compact difference computational strategy to approximate the derivative in space for the model (1)-(3). In Sec. 3, we deduce a semi-discrete approximation system in space, and a high-order explicit strategy for this system in time. And we devote to the high-order positivity-preserving strategies in Sec. 4. In Sec. 5, we will carry out a number of numerical experiments to verify the accuracy, reliability and non-negativity of the proposed methods. The conclusion is provided in Sec. 6.

High-order compact difference discretization in the spatial direction

Discretization of the spatial derivative terms for 1D problem

Firstly, we consider the 1D haptotaxis model for tumor invasion (1)-(3) as follows, where ${\mathrm{\Omega }}_1=\{x|a<x<b;a,b\in \mathbf{R}\}$, $\partial {\mathrm{\Omega }}_1$ is its boundary.

In order to construct high-order numerical method for solving these spatial derivative terms at the right end of the sign of equality in system (4a)-(4e), such as $u_{xx}$, $w_{xx}$ and $\chi {(u{v}_x)}_x$, we denote $p_1=D_1{u}_x-\chi u{v}_x$, $p_2=D_2{w}_x$, and rewrite the Eqs. (4a) and (4c) as follows, $\begin{array}{rlr}& u_t={(p_1)}_x+F_1(u,v),p_1=D_1{u}_x-\chi u{v}_x,& (x,t)\in {\mathrm{\Omega }}_1\times (0,T],\\ & w_t={(p_2)}_x+F_3(u,v,w),p_2=D_2{w}_x,& (x,t)\in {\mathrm{\Omega }}_1\times (0,T].\end{array}$

Next, we divide the physical domain ${\mathrm{\Omega }}_1$ by uniform mesh $N\in {\mathbf{Z}}^{+}$ in 1D space, as shown in Fig. 1. The discrete mesh points $x_i=a+ih$, $0⩽i⩽N$, where $h=(b-a)/N$ means that the spatial mesh step size. Then, we define ${\mathrm{\Gamma }}_h=\{x_i|0⩽i⩽N\}$, and some mesh functions u, v and w on ${\mathrm{\Gamma }}_h$. To save space, we only use u as an example, i.e., for $\mathrm{\forall }u\in {\mathcal{U}}_h=\{u(x_i)|0⩽i⩽N\}$ defined on ${\mathrm{\Gamma }}_h$, denote $u_{i\pm 1}:=u(x_{i\pm 1})$, and $u_{x^{\varsigma }}:=\frac{{\partial }^{\varsigma }u}{\partial x^{\varsigma }}$, $\varsigma \in {\mathbf{Z}}^{+}$, then define ${\delta }_x{u}_i=\frac{1}{2h}(u_{i+1}-u_{i-1})$, and present the following fourth-order symmetrical compact scheme to compute internal points by referring to [26, 33] and references therein. That is, we assume $u(x)\in C^5([x_{i-1},x_{i+1}])$, and denote $\vartheta (\xi )={(1-\xi )}^3(\frac{1}{3}+\xi )$, then have

5

$\begin{array}{rl}& \frac{1}{6}{u}_x(x_{i-1})+\frac{2}{3}{u}_x(x_i)+\frac{1}{6}{u}_x(x_{i+1})={\delta }_x{u}_i+\frac{h^4}{48}\sum _{k=1}^2{\int }_0^1\vartheta (\xi )u_{x^5}(x_i+{(-1)}^{k-1}\xi h)d\xi ,\\ & 1⩽i⩽N-1.\end{array}$

[See PDF for image]

Figure 1

1D space mesh discrete stencil

In order to cooperate with the use of Eq. (5), we give the following fourth-order boundary interpolation formula matching Eq. (5).

Theorem 1

Denote vector$\mathbf{I}{(\xi )}_{1\times 4}=[\frac{1}{18}(19-43\xi ){(1-\xi )}^3,-\frac{3}{2}{(1-\xi )}^4,\frac{13}{18}{(1-\xi )}^4,-\frac{5}{36}{(1-\xi )}^4]$, ${\mathbf{J}}_{1\times 5}^l=[-\frac{53}{36},\frac{43}{18},-\frac{3}{2},\frac{13}{18},-\frac{5}{36}]$, ${\mathbf{J}}^r=-{\mathbf{J}}^l$, then it holds,

6

$\begin{array}{rl}& \frac{2}{3}{u}_x(x_0)+\frac{1}{3}{u}_x(x_1)\\ & =\frac{1}{h}\sum _{i=0}^4{J}_{i+1}^l{u}_i+\frac{h^4}{24}\sum _{k=1}^4{\int }_0^1{k}^5{I}_k(\xi )u_{x^5}(x_0+k\xi h)d\xi ,u(x)\in C^5([x_0,x_4]),\end{array}$

7

$\begin{array}{rl}& \frac{2}{3}{u}_x(x_N)+\frac{1}{3}{u}_x(x_{N-1})\\ & =\frac{1}{h}\sum _{i=0}^4{J}_{i+1}^r{u}_{N-i}+\frac{h^4}{24}\sum _{k=1}^4{\int }_0^1{k}^5{I}_k(\xi )u_{x^5}(x_N-k\xi h)d\xi ,u(x)\in C^5([x_{N-4},x_N]).\end{array}$

Proof

Applying directly the Taylor expansion with Peano remainder, we can expand $u_x(x_1)$, $u_1$, $u_2$, $u_3$ and $u_4$ at $x_0$, $u_x(x_{N-1})$, $u_{N-1}$, $u_{N-2}$, $u_{N-3}$, $u_{N-4}$ at $x_N$, respectively, such as $\begin{array}{rl}& u_x(x_1)={(u_x)}_0+h{(u_{x^2})}_0+\frac{h^2}{2}{(u_{x^3})}_0+\frac{h^3}{6}{(u_{x^4})}_0+\frac{h^4}{6}{\int }_0^1{(1-\xi )}^3{u}_{x^5}(x_0+\xi h)d\xi ,\\ & u(x_k)=\sum _{l=0}^4\frac{{(kh)}^l}{l!}{(u_{x^l})}_0+\frac{{(kh)}^5}{24}{\int }_0^1{(1-\xi )}^4{u}_{x^5}(x_0+k\xi h)d\xi ,k=1,2,3,4.\end{array}$ As well as Hirt enlightenment method, we can refer to the proof of Theorem 2.1 in [34], and easily deduce Eqs. (6) and (7). Therefore, it is omitted here. □

Then, according to the Eq. (5) and Theorem 1, the linear equations composed of (5), (6) and (7) are used to solve $u_x$, $v_x$ and $w_x$, respectively, i.e.,

8

$A_1{U}_x=\frac{1}{h}{A}_{2}U+{\mathcal{T}}_1(h^4),A_1{V}_x=\frac{1}{h}{A}_{2}V+{\mathcal{T}}_2(h^4),A_1{W}_x=\frac{1}{h}{A}_{2}W+{\mathcal{T}}_3(h^4),$ where $\begin{array}{rl}& A_1={\left[\begin{array}{ccccccc}\frac{2}{3}& \frac{1}{3}& 0& \cdots & 0& 0& 0\\ \frac{1}{6}& \frac{2}{3}& \frac{1}{6}& \cdots & 0& 0& 0\\ 0& \frac{1}{6}& \frac{2}{3}& \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& \cdots & \frac{2}{3}& \frac{1}{6}& 0\\ 0& 0& 0& \cdots & \frac{1}{6}& \frac{2}{3}& \frac{1}{6}\\ 0& 0& 0& \cdots & 0& \frac{1}{3}& \frac{2}{3}\end{array}\right]}_{(N+1)\times (N+1)},\\ & A_2={\left[\begin{array}{ccccccc}-\frac{53}{36}& \frac{43}{18}& -\frac{3}{2}& \cdots & 0& 0& 0\\ -\frac{1}{2}& 0& \frac{1}{2}& \cdots & 0& 0& 0\\ 0& -\frac{1}{2}& 0& \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& \cdots & 0& \frac{1}{2}& 0\\ 0& 0& 0& \cdots & -\frac{1}{2}& 0& \frac{1}{2}\\ 0& 0& 0& \cdots & \frac{3}{2}& -\frac{43}{18}& \frac{53}{36}\end{array}\right]}_{(N+1)\times (N+1)},\\ & U_x={[u_x(x_0),u_x(x_1),u_x(x_2),\dots ,u_x(x_{N-2}),u_x(x_{N-1}),u_x(x_N)]}^{⊺},\\ & U={[u_0,u_1,u_2,\dots ,u_{N-2},u_{N-1},u_N]}^{⊺},\\ & \begin{array}{rl}{\mathcal{T}}_1(h^4)& =\frac{h^4}{48}{\int }_0^1\{\sum _{i=1}^{N-1}\sum _{k=1}^2\vartheta (\xi )u_{x^5}(x_i+{(-1)}^{k-1}\xi h)\\ & +2\sum _{k=1}^4{k}^5{I}_k(\xi )[u_{x^5}(x_0+k\xi h)+u_{x^5}(x_N-k\xi h)]\}d\xi ,\end{array}\end{array}$ the forms of $V_x$, V, $W_x$, W, ${\mathcal{T}}_2(h^4)$ and ${\mathcal{T}}_3(h^4)$ are similar to $U_x$, U, and ${\mathcal{T}}_1(h^4)$, respectively. Therefore, by multiplying the inverse matrix $A_1^{-1}$ of $A_1$ on both sides of Eq. (8), respectively, and omitting the truncation error terms $A_1^{-1}{\mathcal{T}}_1(h^4)$, $A_1^{-1}{\mathcal{T}}_2(h^4)$, and $A_1^{-1}{\mathcal{T}}_3(h^4)$, we can obtain the following linear equations with the fourth-order accuracy to approximate the $u_x$, $v_x$ and $w_x$, i.e.,

9

$U_x=\frac{1}{h}{A}_1^{-1}{A}_{2}U,V_x=\frac{1}{h}{A}_1^{-1}{A}_{2}V,W_x=\frac{1}{h}{A}_1^{-1}{A}_{2}W.$

On the basis of the values $u_x$, $v_x$ and $w_x$ computed by Eq. (9), we can present the $p_1(x_i)$ and $p_2(x_i)$ in the interior points $x_i$, $1⩽i⩽N-1$, i.e.,

10

$p_1(x_i)=D_1{u}_x(x_i)-\chi u(x_i)v_x(x_i),p_2(x_i)=D_2{w}_x(x_i),1⩽i⩽N-1.$

Noting that the natural no-flux boundary condition (4e), we obtain

11

$p_1(x_0)=D_1{u}_x(x_0)-\chi u(x_0)v_x(x_0)=0,p_2(x_0)=D_2{w}_x(x_0)=0,$

12

$p_1(x_N)=D_1{u}_x(x_N)-\chi u(x_N)v_x(x_N)=0,p_2(x_N)=D_2{w}_x(x_N)=0.$ Hence, $p_1(x_i)$ and $p_2(x)$, $i=0,1,2,\dots ,N$, have been obtained.

Finally, it is logical to come to the stage of solving ${(p_1)}_x(x)$ and ${(p_2)}_x(x)$. According to Eq. (5) and Theorem 1, combining with Eqs. (10)-(12), we have $A_1{(P_1)}_x=\frac{1}{h}{A}_2{P}_1+{\mathcal{T}}_4(h^4),A_1{(P_2)}_x=\frac{1}{h}{A}_2{P}_2+{\mathcal{T}}_5(h^4),$ where $\begin{array}{rl}& {(P_1)}_x={[{(p_1)}_x(x_0),{(p_1)}_x(x_1),{(p_1)}_x(x_2),\dots ,{(p_1)}_x(x_{N-2}),{(p_1)}_x(x_{N-1}),{(p_1)}_x(x_N)]}^{⊺},\\ & P_1={[{(p_1)}_0,{(p_1)}_1,{(p_1)}_2,\dots ,{(p_1)}_{N-2},{(p_1)}_{N-1},{(p_1)}_N]}^{⊺},\\ & {\mathcal{T}}_4(h^4)=\frac{h^4}{48}{\int }_0^1\{\sum _{i=1}^{N-1}\sum _{k=1}^2\vartheta (\xi ){(p_1)}_{x^5}(x_i+{(-1)}^{k-1}\xi h)\\ & +2\sum _{k=1}^4{k}^5{I}_k(\xi )[{(p_1)}_{x^5}(x_0+k\xi h)+{(p_1)}_{x^5}(x_N-k\xi h)]\}d\xi ,\end{array}$ the forms of ${(P_2)}_x$, $P_2$ and ${\mathcal{T}}_5(h^4)$ are similar to ${(P_1)}_x$, $P_1$ and ${\mathcal{T}}_4(h^4)$, respectively. Similarly, omitting the truncation error terms $A_1^{-1}{\mathcal{T}}_4(h^4)$ and $A_1^{-1}{\mathcal{T}}_5(h^4)$, we can obtain the following fourth-order accuracy strategies to approximate the ${(P_1)}_x$ and ${(P_2)}_x$, i.e.,

13

${(P_1)}_x=\frac{1}{h}{A}_1^{-1}{A}_2{P}_1,{(P_2)}_x=\frac{1}{h}{A}_1^{-1}{A}_2{P}_2.$

Discretization of the spatial derivative terms for 2D problem

Consider the following 2D haptotaxis equations for tumor invasion model (1)-(3), where ${\mathrm{\Omega }}_2=\{(x,y)|a<x<b,c<y<d;a,b,c,d\in \mathbf{R}\}$, $\partial {\mathrm{\Omega }}_2$ is its boundary.

Combined with the structural characteristics of the haptotaxis system (14a)-(14e), to approximate the spatial derivative terms in this system, such as $u_{xx}$, $u_{yy}$, $\chi {(u{v}_x)}_x$, $\chi {(u{v}_y)}_y$, $w_{xx}$ and $w_{yy}$, we denote $p_1=D_1{u}_x-\chi u{v}_x$, $q_1=D_1{u}_y-\chi u{v}_y$, $p_2=D_2{w}_x$, $q_2=D_2{w}_y$, and rewrite the Eqs. (14a) and (14c) as follows,

15

$\begin{array}{rl}& u_t={(p_1)}_x+{(q_1)}_y+F_1(u,v),p_1=D_1{u}_x-\chi u{v}_x,q_1=D_1{u}_y-\chi u{v}_y,\\ & (x,y,t)\in {\mathrm{\Omega }}_2\times (0,T],\end{array}$

16

$\begin{array}{rl}& w_t={(p_2)}_x+{(q_2)}_y+F_3(u,v,w),p_2=D_2{w}_x,q_2=D_2{w}_y,\\ & (x,y,t)\in {\mathrm{\Omega }}_2\times (0,T],\end{array}$ In what follows, we will construct the fourth-order accuracy strategies to approximate $u_x$, $u_y$, $v_x$, $v_y$, $w_x$, $w_y$, $p_1$, $p_2$, ${(p_1)}_x$, ${(p_2)}_x$, $q_1$, $q_2$, ${(q_1)}_y$ and ${(q_2)}_y$. First, we divide the computational domain ${\mathrm{\Omega }}_2$ by uniform mesh $N_x\times N_y$ in 2D space, where $N_x,N_y\in {\mathbf{Z}}^{+}$, as shown in Fig. 2. Analogously, denote $h_x=(b-a)/N_x$, $h_y=(d-c)/N_y$, and $(x_i,y_j)$ stands for the discrete mesh points, where $x_i=a+i{h}_x$, $0⩽i⩽N_x$, $y_j=c+j{h}_y$, $0⩽j⩽N_y$. Define ${\mathrm{\Gamma }}_{h_x\times h_y}=\{(x_i,y_j)|0⩽i⩽N_x,0⩽j⩽N_y\}$, its discrete boundary $\partial {\mathrm{\Gamma }}_{h_x\times h_y}=\{(0,j),(N_x,j)|0⩽j⩽N_y\}\cup \{(i,0),(i,N_y)|0⩽i⩽N_x\}$. Next, we define mesh functions u, v and w on ${\mathrm{\Gamma }}_{h_x\times h_y}$, use u as an example, i.e., for $\mathrm{\forall }u\in {\mathcal{U}}_{h_x\times h_y}=\{u(x_i,y_j)|0⩽i⩽N_x,0⩽j⩽N_y\}$ is defined on ${\mathrm{\Gamma }}_{h_x\times h_y}$, denote $u_{i\pm 1,j}:=u(x_{i\pm 1,j})$, and $u_{i,j\pm 1}:=u(x_{i,j\pm 1})$, we have ${\delta }_x{u}_{i,j}=\frac{1}{2h_x}(u_{i+1,j}-u_{i-1,j})$, and ${\delta }_y{u}_{i,j}=\frac{1}{2h_y}(u_{i,j+1}-u_{i,j-1})$.

[See PDF for image]

Figure 2

2D space mesh discrete stencil

Based on the methods of Refs. [26, 33, 34] and Theorem 1, we can obtain the following fourth-order symmetric compact interior points schemes and the corresponding boundary computational strategies for solving the first-order spatial derivative terms in 2D space. That is, for $1⩽i⩽N_x$, $1⩽j⩽N_y$, denoting $\vartheta (\xi )={(1-\xi )}^3(\frac{1}{3}+\xi )$, and assuming $u(x,y)\in C^{5,0}([x_{i-1},x_{i+1}]\times [y_0,y_{N_y}])$, then

17

$\frac{1}{6}{u}_x(x_{i-1,j})+\frac{2}{3}{u}_x(x_{i,j})+\frac{1}{6}{u}_x(x_{i+1,j})={\delta }_x{u}_{i,j}+\frac{h_x^4}{48}\sum _{k=1}^2{\int }_0^1\vartheta (\xi )u_{x^5}(x_i+{(-1)}^{k-1}\xi h_x,y_j)d\xi .$ Assuming $u(x,y)\in C^{0,5}([x_0,x_{N_x}]\times [y_{i-1},y_{i+1}])$, then

18

$\frac{1}{6}{u}_y(x_{i,j-1})+\frac{2}{3}{u}_y(x_{i,j})+\frac{1}{6}{u}_y(x_{i,j+1})={\delta }_y{u}_{i,j}+\frac{h_y^4}{48}\sum _{k=1}^2{\int }_0^1\vartheta (\xi )u_{y^5}(x_i,y_j+{(-1)}^{k-1}\xi h_y)d\xi .$

Theorem 2

Denote vector$\mathbf{I}{(\xi )}_{1\times 4}=[\frac{1}{18}(19-43\xi ){(1-\xi )}^3,-\frac{3}{2}{(1-\xi )}^4,\frac{13}{18}{(1-\xi )}^4,-\frac{5}{36}{(1-\xi )}^4]$, ${\mathbf{J}}_{1\times 5}^l=[-\frac{53}{36},\frac{43}{18},-\frac{3}{2},\frac{13}{18},-\frac{5}{36}]$, ${\mathbf{J}}^r=-{\mathbf{J}}^l$, and we have

I. If$u(x,y)\in C^{5,0}([x_0,x_4]\times [y_0,y_{N_y}])$, then

19

$\frac{2}{3}{u}_x(x_0,y_j)+\frac{1}{3}{u}_x(x_1,y_j)=\frac{1}{h_x}\sum _{i=0}^4{J}_{i+1}^l{u}_{i,j}+\frac{h_x^4}{24}\sum _{k=1}^4{\int }_0^1{k}^5{I}_k(\xi )u_{x^5}(x_0+k\xi h_x,y_j)d\xi ,$

II. If$u(x,y)\in C^{5,0}([x_{N_x-4},x_{N_x}]\times [y_0,y_{N_y})$, then

20

$\begin{array}{rl}\frac{2}{3}{u}_x(x_{N_x},y_j)+\frac{1}{3}{u}_x(x_{N_x-1},y_j)& =\frac{1}{h_x}\sum _{i=0}^4{J}_{i+1}^r{u}_{N_x-i,j}\\ & +\frac{h_x^4}{24}\sum _{k=1}^4{\int }_0^1{k}^5{I}_k(\xi )u_{x^5}(x_{N_x}-k\xi h_x,y_j)d\xi ,\end{array}$

III. If$u(x,y)\in C^{0,5}([x_0,x_{N_x}]\times [y_0,y_4])$, then

21

$\frac{2}{3}{u}_y(x_i,y_0)+\frac{1}{3}{u}_y(x_i,y_1)=\frac{1}{h_y}\sum _{j=0}^4{J}_{j+1}^l{u}_{i,j}+\frac{h_x^4}{24}\sum _{k=1}^4{\int }_0^1{k}^5{I}_k(\xi )u_{y^5}(x_i,y_0+k\xi h_y)d\xi ,$

IV. If$u(x,y)\in C^{0,5}([x_0,x_{N_x}]\times [y_{N_y-4},y_{N_y}])$, then

22

$\begin{array}{rl}\frac{2}{3}{u}_y(x_i,y_{N_y})+\frac{1}{3}{u}_y(x_i,y_{N_y-1})& =\frac{1}{h_y}\sum _{j=0}^4{J}_{j+1}^r{u}_{i,N_y-j}\\ & +\frac{h_y^4}{24}\sum _{k=1}^4{\int }_0^1{k}^5{I}_k(\xi )u_{y^5}(x_i,y_{N_y}-k\xi h_y)d\xi .\end{array}$

Proof

The proof of this theorem is similar to that of Theorem 2.1 in [34], which is omitted here. □

Now, let us solve $u_x$, $u_y$, $v_x$, $v_y$, $w_x$, $w_y$, $p_1$, ${(p_1)}_x$, $p_2$, ${(p_2)}_x$, $q_1$, ${(q_1)}_y$, $q_2$ and ${(q_2)}_y$, by the above definitions, Eqs. (17) and (18), as well as Theorem 2, for $0⩽i⩽N_x$, $0⩽j⩽N_y$, we have

23

$\begin{array}{rl}& A_1{U}_{xj}=\frac{1}{h_x}{A}_2{U}_j+{({\mathcal{T}}_1)}_{xj}(h_x^4),A_1{V}_{xj}=\frac{1}{h_x}{A}_2{V}_j+{({\mathcal{T}}_2)}_{xj}(h_x^4),\\ & A_1{W}_{xj}=\frac{1}{h_x}{A}_2{W}_j+{({\mathcal{T}}_3)}_{xj}(h_x^4),\end{array}$

24

$\begin{array}{rl}& U_{yi}{B}_1=\frac{1}{h_y}{U}_i{B}_2+{({\mathcal{T}}_1)}_{yi}(h_y^4),V_{yi}{B}_1=\frac{1}{h_y}{V}_i{B}_2+{({\mathcal{T}}_2)}_{yi}(h_y^4),\\ & W_{yi}{B}_1=\frac{1}{h_y}{W}_i{B}_2+{({\mathcal{T}}_3)}_{yi}(h_y^4),\end{array}$ where $\begin{array}{rl}& B_1=A_1^{⊺},B_2=A_2^{⊺},\\ & U_{xj}={[u_x(x_0,y_j),u_x(x_1,y_j),u_x(x_2,y_j),\dots ,u_x(x_{N_x-2},y_j),u_x(x_{N_x-1},y_j),u_x(x_{N_x},y_j)]}^{⊺},\\ & U_j={[u_{0,j},u_{1,j},u_{2,j},\dots ,u_{N_x-2,j},u_{N_x-1,j},u_{N_x,j}]}^{⊺},\\ & U_{yi}=[u_y(x_i,y_0),u_y(x_i,y_1),u_y(x_i,y_2),\dots ,u_y(x_i,y_{N_y-2}),u_y(x_i,y_{N_y-1}),u_y(x_i,y_{N_y})],\\ & U_i=[u_{i,0},u_{i,1},u_{i,2},\dots ,u_{i,N_y-2},u_{i,N_y-1},u_{i,N_y}],\\ & {({\mathcal{T}}_1)}_{xj}(h_x^4)=\frac{h_x^4}{48}{\int }_0^1\{\sum _{i=1}^{N_x-1}\sum _{k=1}^2\vartheta (\xi )u_{x^5}(x_i+{(-1)}^{k-1}\xi h_x,y_j)\\ & +2\sum _{k=1}^4{k}^5{I}_k(\xi )[u_{x^5}(x_0+k\xi h_x,y_j)+u_{x^5}(x_{N_x}-k\xi h_x,y_j)]\}d\xi ,\\ & {({\mathcal{T}}_1)}_{yi}(h_y^4)=\frac{h_y^4}{48}{\int }_0^1\{\sum _{j=1}^{N_y-1}\sum _{k=1}^2\vartheta (\xi )u_{y^5}(x_i,y_i+{(-1)}^{k-1}\xi h_y)\\ & +2\sum _{k=1}^4{k}^5{I}_k(\xi )[u_{y^5}(x_i,y_0+k\xi h_y)+u_{y^5}(x_i,y_{N_y}-k\xi h_y)]\}d\xi ,\end{array}$ the forms of $V_{xj}$, $V_j$, $V_{yi}$, $V_i$, $W_{xj}$, $W_j$, $W_{yi}$, $W_i$, ${({\mathcal{T}}_2)}_{xj}(h_x^4)$, ${({\mathcal{T}}_2)}_{yi}(h_y^4)$, ${({\mathcal{T}}_3)}_{xj}(h_x^4)$ and ${({\mathcal{T}}_3)}_{yi}(h_y^4)$ are similar to $U_{xj}$, $U_j$, $U_{yi}$, $U_i$, ${({\mathcal{T}}_1)}_{xj}(h_x^4)$ and ${({\mathcal{T}}_1)}_{yi}(h_y^4)$, respectively. Consequently, multiplying the $A_1^{-1}$ on both sides of Eq. (23), the $B_1^{-1}$ on both sides of Eq. (24), respectively, and omitting the truncation error terms $A_1^{-1}{({\mathcal{T}}_1)}_{xj}(h_x^4)$, $A_1^{-1}{({\mathcal{T}}_2)}_{xj}(h_x^4)$, $A_1^{-1}{({\mathcal{T}}_3)}_{xj}(h_x^4)$, $B_1^{-1}{({\mathcal{T}}_1)}_{yi}(h_y^4)$, $B_1^{-1}{({\mathcal{T}}_2)}_{yi}(h_y^4)$, and $B_1^{-1}{({\mathcal{T}}_3)}_{yi}(h_y^4)$, we have $\begin{array}{rl}& U_{xj}=\frac{1}{h_x}{A}_1^{-1}{A}_2{U}_j,V_{xj}=\frac{1}{h_x}{A}_1^{-1}{A}_2{V}_j,W_{xj}=\frac{1}{h_x}{A}_1^{-1}{A}_2{W}_j,0⩽j⩽N_y,\\ & U_{yi}=\frac{1}{h_y}{U}_i{B}_2{B}_1^{-1},V_{yi}=\frac{1}{h_y}{V}_i{B}_2{B}_1^{-1},W_{yi}=\frac{1}{h_y}{W}_i{B}_2{B}_1^{-1},0⩽i⩽N_x.\end{array}$ Thus, we can also obtain the linear equations with the fourth-order accuracy to approximate the $u_x$, $u_y$, $v_x$, $v_y$, $w_x$ and $w_y$ as follows,

25

$U_x=\frac{1}{h_x}{A}_1^{-1}{A}_{2}U,V_x=\frac{1}{h_x}{A}_1^{-1}{A}_{2}V,W_x=\frac{1}{h_x}{A}_1^{-1}{A}_{2}W,$

26

$U_y=\frac{1}{h_y}U{B}_2{B}_1^{-1},V_y=\frac{1}{h_y}V{B}_2{B}_1^{-1},W_y=\frac{1}{h_y}W{B}_2{B}_1^{-1},$ where $U_x={\left({(u_x)}_{i,j}\right)}_{(N_x+1)\times (N_y+1)}$, $U_y={\left({(u_y)}_{i,j}\right)}_{(N_x+1)\times (N_y+1)}$, $V_x={\left({(v_x)}_{i,j}\right)}_{(N_x+1)\times (N_y+1)}$, $V_y={\left({(v_y)}_{i,j}\right)}_{(N_x+1)\times (N_y+1)}$, $W_x={\left({(w_x)}_{i,j}\right)}_{(N_x+1)\times (N_y+1)}$, $W_y={\left({(w_y)}_{i,j}\right)}_{(N_x+1)\times (N_y+1)}$, $U={(u_{i,j})}_{(N_x+1)\times (N_y+1)}$, $V={(v_{i,j})}_{(N_x+1)\times (N_y+1)}$, and $W={(w_{i,j})}_{(N_x+1)\times (N_y+1)}$. Furthermore, we can obtain the following values of $p_1$, $p_2$, $q_1$ and $q_2$ at the interior point $(x_i,y_j)$, where $1⩽i⩽N_x-1$, $1⩽j⩽N_y-1$, i.e.,

27

$p_1(x_i,y_j)=D_1{u}_x(x_i,y_j)-\chi u(x_i,y_j)v_x(x_i,y_j),p_2(x_i,y_j)=D_2{w}_x(x_i,y_j),$

28

$q_1(x_i,y_j)=D_1{u}_y(x_i,y_j)-\chi u(x_i,y_j)v_y(x_i,y_j),q_2(x_i,y_j)=D_2{w}_y(x_i,y_j).$ Combining with the boundary condition (14e), we can compute the values of $p_1$, $p_2$, $q_1$ and $q_2$ at the boundary points $\partial {\mathrm{\Gamma }}_{h_x\times h_y}$ as follows,

29

$p_1(x_0,y_j)=D_1{u}_x(x_0,y_j)-\chi u(x_0,y_j)v_x(x_0,y_j)=0,p_2(x_0,y_j)=D_2{w}_x(x_0,y_j)=0,$

30

$p_1(x_{N_x},y_j)=D_1{u}_x(x_{N_x},y_j)-\chi u(x_{N_x},y_j)v_x(x_{N_x},y_j)=0,p_2(x_{N_x},y_j)=D_2{w}_x(x_{N_x},y_j)=0,$

31

$q_1(x_i,y_0)=D_1{u}_y(x_i,y_0)-\chi u(x_i,y_0)v_y(x_i,y_0)=0,q_2(x_i,y_0)=D_2{w}_y(x_i,y_0)=0,$

32

$q_1(x_i,y_{N_y})=D_1{u}_y(x_i,y_{N_y})-\chi u(x_i,y_{N_y})v_y(x_i,y_{N_y})=0,q_2(x_i,y_{N_y})=D_2{w}_y(x_i,y_{N_y})=0.$

Finally, on account of Eqs. (17)-(18), as well as Theorem 2, combining with Eqs. (27)-(32), we have $\begin{array}{rl}& A_1{(P_1)}_{xj}=\frac{1}{h_x}{A}_2{(P_1)}_j+{({\mathcal{T}}_4)}_{xj}(h_x^4),A_1{(P_2)}_{xj}=\frac{1}{h_x}{A}_2{(P_2)}_j+{({\mathcal{T}}_5)}_{xj}(h_x^4),0⩽j⩽N_y,\\ & {(Q_1)}_{yi}{B}_1=\frac{1}{h_y}{(Q_1)}_i{B}_2+{({\mathcal{T}}_4)}_{yi}(h_y^4),{(Q_2)}_{yi}{B}_1=\frac{1}{h_y}{(Q_2)}_i{B}_2+{({\mathcal{T}}_5)}_{yi}(h_y^4),0⩽i⩽N_x,\end{array}$ where $\begin{array}{rl}& {(P_1)}_{xj}={[{(p_1)}_x(x_0,y_j),{(p_1)}_x(x_1,y_j),\dots ,{(p_1)}_x(x_{N_x-1},y_j),{(p_1)}_x(x_{N_x},y_j)]}^{⊺},\\ & {(P_1)}_j={[{(p_1)}_{0,j},{(p_1)}_{1,j},{(p_1)}_{2,j},\dots ,{(p_1)}_{N_x-2,j},{(p_1)}_{N_x-1,j},{(p_1)}_{N_x,j}]}^{⊺},\\ & {(Q_1)}_{yi}=[{(q_1)}_y(x_i,y_0),{(q_1)}_y(x_i,y_1),\dots ,{(q_1)}_y(x_i,y_{N_y-1}),{(q_1)}_y(x_i,y_{N_y})],\\ & {(Q_1)}_i=[{(q_1)}_{i,0},{(q_1)}_{i,1},{(q_1)}_{i,2},\dots ,{(q_1)}_{i,N_y-2},{(q_1)}_{i,N_y-1},{(q_1)}_{i,N_y}],\\ & {({\mathcal{T}}_4)}_{xj}(h_x^4)=\frac{h_x^4}{48}{\int }_0^1\{\sum _{i=1}^{N_x-1}\sum _{k=1}^2\vartheta (\xi ){(p_1)}_{x^5}(x_i+{(-1)}^{k-1}\xi h_x,y_j)\\ & +2\sum _{k=1}^4{k}^5{I}_k(\xi )[{(p_1)}_{x^5}(x_0+k\xi h_x,y_j)+{(p_1)}_{x^5}(x_{N_x}-k\xi h_x,y_j)]\}d\xi ,\\ & {({\mathcal{T}}_4)}_{yi}(h_y^4)=\frac{h_y^4}{48}{\int }_0^1\{\sum _{j=1}^{N_y-1}\sum _{k=1}^2\vartheta (\xi ){(q_1)}_{y^5}(x_i,y_i+{(-1)}^{k-1}\xi h_y)\\ & +2\sum _{k=1}^4{k}^5{I}_k(\xi )[{(q_1)}_{y^5}(x_i,y_0+k\xi h_y)+{(q_1)}_{y^5}(x_i,y_{N_y}-k\xi h_y)]\}d\xi ,\end{array}$ the forms of ${(P_2)}_{xj}$, ${(P_2)}_j$, ${(Q_2)}_{yi}$, ${(Q_2)}_i$, ${({\mathcal{T}}_5)}_{xj}(h_x^4)$ and ${({\mathcal{T}}_5)}_{yi}(h_y^4)$ are similar to ${(P_1)}_{xj}$, ${(P_1)}_j$, ${(Q_1)}_{yi}$, ${(Q_1)}_i$, ${({\mathcal{T}}_4)}_{xj}(h_x^4)$ and ${({\mathcal{T}}_4)}_{yi}(h_y^4)$, respectively. Analogously, omitting the truncation error terms $A_1^{-1}{({\mathcal{T}}_4)}_{xj}(h_x^4)$, $A_1^{-1}{({\mathcal{T}}_5)}_{xj}(h_x^4)$, $B_1^{-1}{({\mathcal{T}}_4)}_{yi}(h_y^4)$, and $B_1^{-1}{({\mathcal{T}}_5)}_{yi}(h_y^4)$, we can derive the following fourth-order accuracy schemes to approximate ${(p_1)}_x$, ${(p_2)}_x$, ${(q_1)}_y$ and ${(q_2)}_y$, i.e.,

33

${(P_1)}_x=\frac{1}{h_x}{A}_1^{-1}{A}_2{P}_1,{(P_2)}_x=\frac{1}{h_x}{A}_1^{-1}{A}_2{P}_2,$

34

${(Q_1)}_y=\frac{1}{h_y}{Q}_1{B}_2{B}_1^{-1},{(Q_2)}_y=\frac{1}{h_y}{Q}_2{B}_2{B}_1^{-1},$ in which, $\begin{array}{rl}& {(P_1)}_x={\left({\left({(p_1)}_x\right)}_{i,j}\right)}_{(N_x+1)\times (N_y+1)},{(Q_1)}_y={\left({\left({(q_1)}_y\right)}_{i,j}\right)}_{(N_x+1)\times (N_y+1)},\\ & {(P_2)}_x={\left({\left({(p_2)}_x\right)}_{i,j}\right)}_{(N_x+1)\times (N_y+1)},{(Q_2)}_y={\left({\left({(q_2)}_y\right)}_{i,j}\right)}_{(N_x+1)\times (N_y+1)},\\ & P_1={\left({(p_1)}_{i,j}\right)}_{(N_x+1)\times (N_y+1)},P_2={\left({(p_2)}_{i,j}\right)}_{(N_x+1)\times (N_y+1)},\\ & Q_1={\left({(q_1)}_{i,j}\right)}_{(N_x+1)\times (N_y+1)},Q_2={\left({(q_2)}_{i,j}\right)}_{(N_x+1)\times (N_y+1)}.\end{array}$

Computational approach in the temporal direction

Base on the dissociation for the spatial derivative terms above, the semi-discrete schemes of the systems (4a)-(4e) and (14a)-(14e) are equivalent to the following ODE systems, i.e.,

35

$\{\begin{array}{ll}\frac{dU}{dt}=L(U),& t\in (0,T],\\ \frac{dV}{dt}=L(V),& t\in (0,T],\\ \frac{dW}{dt}=L(W),& t\in (0,T].\end{array}$ According to Eqs. (9)-(13), we can obtain the forms of $L(U)$, $L(V)$, and $L(W)$ in 1D space as follows, $\{\begin{array}{l}L(U)=\frac{1}{h^2}{A}_1^{-1}{A}_2(D_1{A}_1^{-1}{A}_{2}U-\chi U{A}_1^{-1}{A}_{2}V)+F_1(U,V),\\ L(V)=F_2(U,V,W),\\ L(W)=\frac{1}{h^2}{A}_1^{-1}{A}_2(D_2{A}_1^{-1}{A}_{2}W)+F_3(U,V,W),\end{array}$ where $U={[u_0,u_1,\dots ,u_N]}^{⊺}$, $V={[v_0,v_1,\dots ,v_N]}^{⊺}$, and $W={[w_0,w_1,\dots ,w_N]}^{⊺}$. In the same way, base on the Eqs. (25)-(34), we can also obtain the following forms of $L(U)$, $L(V)$, and $L(W)$ in 2D space, $\{\begin{array}{l}L(U)=\frac{1}{h_x^2}{A}_1^{-1}{A}_2(D_1{A}_1^{-1}{A}_{2}U-\chi U{A}_1^{-1}{A}_{2}V)\\ +\frac{1}{h_y^2}(D_{1}U{B}_2{B}_1^{-1}-\chi UV{B}_2{B}_1^{-1})B_2{B}_1^{-1}+F_1(U,V),\\ L(V)=F_2(U,V,W),\\ L(W)=\frac{1}{h_x^2}{A}_1^{-1}{A}_2(D_2{A}_1^{-1}{A}_{2}W)+\frac{1}{h_y^2}(D_{2}W{B}_2{B}_1^{-1})B_2{B}_1^{-1}+F_3(U,V,W),\end{array}$ where $U={(u_{i,j})}_{(N_x+1)\times (N_y+1)}$, $V={(v_{i,j})}_{(N_x+1)\times (N_y+1)}$, and $W={(w_{i,j})}_{(N_x+1)\times (N_y+1)}$.

In what follows, we discuss the high-accuracy discrete strategies to solve the system (35). Firstly, divide the temporal interval $[0,T]$ by uniform mesh $M\in {\mathbf{Z}}^{+}$, the temporal step-length $\tau =T/M$, and $t_n=n\tau$, $0⩽n⩽M$. In order to obtain an accurate and stable approximation method in the temporal direction, we adopt the optimal nonlinear third-order SSP Runge-Kutta method [31, 32], namely

36

${\mathrm{\Theta }}^{(1)}={\mathrm{\Theta }}^n+\tau L({\mathrm{\Theta }}^n),$

37

${\mathrm{\Theta }}^{(2)}=\frac{3}{4}{\mathrm{\Theta }}^n+\frac{1}{4}{\mathrm{\Theta }}^{(1)}+\frac{1}{4}\tau L({\mathrm{\Theta }}^{(1)}),$

38

${\mathrm{\Theta }}^{n+1}=\frac{1}{3}{\mathrm{\Theta }}^n+\frac{2}{3}{\mathrm{\Theta }}^{(2)}+\frac{2}{3}\tau L({\mathrm{\Theta }}^{(2)}),$ where $\mathrm{\Theta }\in \{U,V,W\}$, and $L({\mathrm{\Theta }}^n)$ is the discrete result of all terms in Eq. (35) at time $t_n$.

Fourth-order positivity-preserving approaches

Since $u(\mathbf{x},t)$, $v(\mathbf{x},t)$ and $w(\mathbf{x},t)$ in the tumor invasion haptotaxis system (1)-(3) represent the densities of TC and ECM, the concentration of MDE, respectively, their values cannot be negative in the biological sense. When the newly constructed computational strategies are applied to the numerical approximation systems (1)-(3), these local regions may produce small negative values due to the mathematical properties of the constructed schemes, such as stability and reliability, as well as the inherent properties of the problem itself, such as large gradient solutions or sharp structural solution regions. Therefore, in this section, to keep consistency with the fourth-order accuracy of the spatial direction, we will use the idea of high-order numerical integration to establish positivity-preserving algorithms with fourth-order accuracy without loss at all time-steps.

On the one hand, for the 1D problem, we define the average solution for the local computational domain ${({\mathrm{\Upsilon }}_h)}_i=[x_{i-1},x_{i+1}],1⩽i⩽N-1$ on each time step as follows,

39

${\bar{\mathrm{\Theta }}}_i^n=\frac{1}{{\mathrm{\Delta }}_h{x}_i}{\int }_{{({\mathrm{\Upsilon }}_h)}_i}\mathrm{\Theta }(x,t_n)dx,$ where $\mathrm{\Theta }\in \{U,V,W\}$, $1⩽n⩽M$, ${\mathrm{\Delta }}_h{x}_i=x_{i+1}-x_{i-1}$, $1⩽i⩽N-1$. Based on the idea of numerical integration and Simpson formula, we have

40

$\begin{array}{rl}{\int }_{{({\mathrm{\Upsilon }}_h)}_i}\mathrm{\Theta }(x,t_n)dx& =\frac{1}{6}{\mathrm{\Delta }}_h{x}_i[\mathrm{\Theta }(x_{i-1},t_n)+4\mathrm{\Theta }(x_i,t_n)+\mathrm{\Theta }(x_{i+1},t_n)]\\ & -\frac{{\mathrm{\Delta }}_h{x}_i}{180}{\left(\frac{{\mathrm{\Delta }}_h{x}_i}{2}\right)}^4{\mathrm{\Theta }}^{(4)}(\eta ,t_n),\end{array}$ where $\eta \in (x_{i-1},x_{i+1})$, $1⩽i⩽N-1$, $1⩽n⩽M$. Substituting Eq. (40) into Eq. (39), we get

41

${\bar{\mathrm{\Theta }}}_i^n=\frac{1}{6}[\mathrm{\Theta }(x_{i-1},t_n)+4\mathrm{\Theta }(x_i,t_n)+\mathrm{\Theta }(x_{i+1},t_n)]-\frac{{({\mathrm{\Delta }}_h{x}_i)}^4}{2880}{\mathrm{\Theta }}^{(4)}(\eta ,t_n),\eta \in (x_{i-1},x_{i+1}).$ Eliminating the error remainder $-\frac{{({\mathrm{\Delta }}_h{x}_i)}^4}{2880}{\mathrm{\Theta }}^{(4)}(\eta ,t_n)$ in Eq. (41), and substituting the values of $\mathrm{\Theta }(x_{i-1},t_n)$, $\mathrm{\Theta }(x_i,t_n)$ and $\mathrm{\Theta }(x_{i+1},t_n)$ into Eq. (41) to compute via the proposed method (9)-(13) and (36)-(38), we can derive the following discrete type linear weighting compact formula with fourth-order accuracy to approximate the average solution ${\mathrm{\Theta }}_i^n$, i.e.,

42

${\bar{\mathrm{\Theta }}}_i^n\approx \frac{1}{6}({\mathrm{\Theta }}_{i-1}^n+4{\mathrm{\Theta }}_i^n+{\mathrm{\Theta }}_{i+1}^n),1⩽i⩽N-1,1⩽n⩽M.$ And then, referring to the linear scaling limiter methods in [35, 36], the following positivity-preserving limiter (PPL) is structured to effectively eliminate the negative values generated during the computational process, namely,

43

${\tilde{\mathrm{\Theta }}}_i^n=\lambda {\mathrm{\Theta }}_i^n+(1-\lambda ){\bar{\mathrm{\Theta }}}_i^n,1⩽i⩽N-1,1⩽n⩽M,$ where $\lambda =min\left\{\right|\frac{{\bar{\mathrm{\Theta }}}_i^n}{{\bar{\mathrm{\Theta }}}_i^n-{({\mathrm{\Theta }}_{min})}_i^n}|,1\},{({\mathrm{\Theta }}_{min})}_i^n=min\{{\mathrm{\Theta }}_{i-1}^n,{\mathrm{\Theta }}_i^n,{\mathrm{\Theta }}_{i+1}^n\}.$ By using the fourth-order accuracy linear weighting compact formula and PPL, we can get the positivity-preserving solution $\tilde{\mathrm{\Theta }}⩾0$ without losing the fourth-order accuracy obtained by the proposed numerical method (9)-(13).

On the other hand, we consider the 2D problem for this model (1)-(3). Analogously, we define the following average solutions on each time step for local areas ${({\mathrm{\Upsilon }}_{h_x\times h_y})}_{i,j}=[x_{i-1},x_{i+1}]$$\times [y_{j-1},y_{j+1}]$, $1⩽i⩽N_x-1$, $1⩽j⩽N_y-1$, i.e., ${\bar{\mathrm{\Theta }}}_{i,j}^n=\frac{1}{{\mathrm{\Delta }}_{h_x}{x}_i{\mathrm{\Delta }}_{h_y}{y}_j}{\iint }_{{({\mathrm{\Upsilon }}_{h_x\times h_y})}_{i,j}}\mathrm{\Theta }(x,y,t_n)dxdy,$ where $1⩽n⩽M$, ${\mathrm{\Delta }}_{h_x}{x}_i=x_{i+1}-x_{i-1}$, ${\mathrm{\Delta }}_{h_y}{y}_j=y_{j+1}-y_{j-1}$, $1⩽i⩽N_x-1$, $1⩽j⩽N_y-1$. According to the derivation method of the above 1D problem and the idea of numerical integration of multivariate functions, for the 2D problem, we can also obtain the following fourth-order accuracy linear weighted compact approximation formula, namely,

44

$\begin{array}{rl}& {\bar{\mathrm{\Theta }}}_{i,j}^n\approx \frac{1}{6}({\mathrm{\Theta }}_{i+1,j}^n+{\mathrm{\Theta }}_{i-1,j}^n+2{\mathrm{\Theta }}_{i,j}^n+{\mathrm{\Theta }}_{i,j+1}^n+{\mathrm{\Theta }}_{i,j-1}^n),\\ & 1⩽i⩽N_x-1,1⩽j⩽N_y-1,1⩽n⩽M.\end{array}$ Thus, we can establish the following PPL to filter the local negative values of ${\mathrm{\Theta }}_{i,j}^n$ for 2D space, i.e.,

45

${\tilde{\mathrm{\Theta }}}_{i,j}^n=\lambda {\mathrm{\Theta }}_{i,j}^n+(1-\lambda ){\bar{\mathrm{\Theta }}}_{i,j}^n,$ where $\lambda =min\left\{\right|\frac{{\bar{\mathrm{\Theta }}}_{i,j}^n}{{\bar{\mathrm{\Theta }}}_{i,j}^n-{({\mathrm{\Theta }}_{min})}_{i,j}^n}|,1\},{({\mathrm{\Theta }}_{min})}_{i,j}^n=min\{{\mathrm{\Theta }}_{i,j}^n,{\mathrm{\Theta }}_{i-1,j}^n,{\mathrm{\Theta }}_{i+1,j}^n,{\mathrm{\Theta }}_{i,j+1}^n,{\mathrm{\Theta }}_{i,j-1}^n\}.$

Numerical experiments

In this part, we carry out several numerical experiments: on the one hand, to verify the accuracy and reliability of the proposed methods, on the other hand, to simulate and analyze the complex biological evolution dynamics of tumor invasion. Noting that the computation and simulation of all the following examples are programmed in Fortran language and run on the desktop computer with Intel(R) Core(TM) i5-3470S CPU @ 2.90 GHz 2.90 GHz, RAM 8.00 GB. For the convenience of testing, we define the $L_{\mathrm{\infty }}$ and $L_2$ norm errors in 1D and 2D spaces, as well as convergence rates, as follows $\begin{array}{rl}1Dspace:& L_{\mathrm{\infty }}(\mathrm{\Theta })=\underset{\begin{array}{c}0⩽i⩽N\end{array}}{max}|{\mathrm{\Theta }}_i^M-\mathrm{\Theta }(x_i,t_M)|,\\ & L_2(\mathrm{\Theta })=\sqrt{h\sum _{i=0}^N{[{\mathrm{\Theta }}_i^M-\mathrm{\Theta }(x_i,t_M)]}^2},\\ 2Dspace:& L_{\mathrm{\infty }}(\mathrm{\Theta })=\underset{\begin{array}{c}0⩽i⩽N_x\\ 0⩽j⩽N_y\end{array}}{max}|{\mathrm{\Theta }}_{i,j}^M-\mathrm{\Theta }(x_i,y_j,t_M)|,\\ & L_2(\mathrm{\Theta })=\sqrt{h_x{h}_y\sum _{i=0}^{N_x}\sum _{j=0}^{N_y}{[{\mathrm{\Theta }}_{i,j}^M-\mathrm{\Theta }(x_i,y_j,t_M)]}^2},\\ & Rat{e}_{\mathrm{\infty }}=\frac{log[L_{\mathrm{\infty }}({\mathrm{\Theta }}_{m_1})/L_{\mathrm{\infty }}({\mathrm{\Theta }}_{m_2})]}{log(m_1/m_2)},or\\ & Rat{e}_2=\frac{log[L_2({\mathrm{\Theta }}_{m_1})/L_2({\mathrm{\Theta }}_{m_2})]}{log(m_1/m_2)},\end{array}$ in which, $\mathrm{\Theta }\in \{u,v,w\}$, $\mathrm{\Theta }(x_i,t_M)$, $\mathrm{\Theta }(x_i,y_j,t_M)$ and ${\mathrm{\Theta }}_i^M$, ${\mathrm{\Theta }}_{i,j}^M$ stand for the exact and numerical solutions at $(x_i,t_M)$ and $(x_i,y_j,t_M)$, respectively. $m_1$ and $m_2$ are the mesh step sizes.

Accuracy test without PPL for 1D problem

Since the original system (1)-(3) is difficult to obtain its analytical solution, to verify the accuracy of the proposed method without PPL, we add the additional fluxes (AFs) $g_k$$(k=1,2,3)$ to the right end of each equation of the original system (1)-(3), such that it satisfies the following preset analytical solutions, $\begin{array}{rl}& u(x,t)=exp(-t-\frac{x^2}{\epsilon }),v(x,t)=1-\frac{1}{2}exp(-t-\frac{x^2}{\epsilon }),\\ & w(x,t)=\frac{1}{2}exp(-t-\frac{x^2}{\epsilon }),0⩽x⩽1,t⩾0.\end{array}$$u(x,0)$, $v(x,0)$ and $w(x,0)$ are derived from the analytical solutions above, and the boundary condition (4e) is used. The AFs $g_k$$(k=1,2,3)$ for the original system (1)-(3) with the initial-boundary condition under Cases 1, 4, and 5 are given out in Table 2.

Table 2. Expressions of $g_k(k=1,2,3)$ for Example 5.1 Cases 1, 4, and 5

	AFs	Expressions	AFs
Case 1	$g_1$	$[-1-D_1(-\frac{2}{\epsilon }+\frac{4x^2}{{\epsilon }^2}\left)\right]exp(-t-\frac{x^2}{\epsilon })+\chi (\frac{1}{\epsilon }-\frac{4x^2}{{\epsilon }^2})exp(-2t-\frac{2x^2}{\epsilon })$
	$g_2$	$\frac{1}{2}(1+\eta )exp(-t-\frac{x^2}{\epsilon })-\frac{1}{4}\eta exp(-2t-\frac{2x^2}{\epsilon })$
	$g_3$	$[-\frac{1}{2}-D_2(-\frac{1}{\epsilon }+\frac{2x^2}{{\epsilon }^2})-\alpha +\frac{\beta }{2}]exp(-t-\frac{x^2}{\epsilon })$	$g_3$	Case 4
Case 5	$g_1$	$[-1-D_1(-\frac{2}{\epsilon }+\frac{4x^2}{{\epsilon }^2}\left)\right]exp(-t-\frac{x^2}{\epsilon })+\left[\chi \right(\frac{1}{\epsilon }-\frac{4x^2}{{\epsilon }^2})+\frac{{\mu }_1}{2}]exp(-2t-\frac{2x^2}{\epsilon })$	$g_1$
	$g_2$	$\frac{1}{2}(1+\eta +{\mu }_2)exp(-t-\frac{x^2}{\epsilon })-\frac{1}{4}(\eta +{\mu }_2)exp(-2t-\frac{2x^2}{\epsilon })$	$g_2$
	$g_3$	$[-\frac{1}{2}-D_2(-\frac{1}{\epsilon }+\frac{2x^2}{{\epsilon }^2})-\alpha +\frac{\beta }{2}]exp(-t-\frac{x^2}{\epsilon })+\frac{\alpha }{2}exp(-2t-\frac{2x^2}{\epsilon })$

In this computation, we take parameters $\alpha =0.1$, $\beta =0.5$, $\eta =10$, $\epsilon =0.01$, $D_1=D_2=0.001$, $\chi =0.01$, ${\mu }_1=0.1$, and ${\mu }_2=0.5$, and $\tau =h^2$, which was used to compute this problem in [11, 12], and $\tau =h^{4/3}$ is used in this paper. The $L_{\mathrm{\infty }}$ and $L_2$ norm errors and convergence rates, which are obtained separately to solve Example 5.1 Cases 1, 4 and 5 with $T=0.5$ by using the proposed method, are listed in Tables 3, 4, 5 and 6, and the results are compared with those obtained in [11, 12]. It can be seen that under the same computational parameters, the errors obtained by our method is 1-3 orders of magnitude smaller than that obtained by the methods in [11, 12]. At the same time, the optimal convergence rates are 2 in [11, 12], while our method can achieve the fourth-order convergence rate. Therefore, our method is more accurate than the existing literature to solve the haptotaxis system (1)-(3).

Table 3. $L_{\mathrm{\infty }}$ error and convergence rate of the proposed method for Example 5.1 Case 1 with $T=0.5$

	N	Ref. [12]		Method I [11]		Method II [11]		Present
		$L_{\mathrm{\infty }}$	$Rat{e}_{\mathrm{\infty }}$	$L_{\mathrm{\infty }}$	$Rat{e}_{\mathrm{\infty }}$	$L_{\mathrm{\infty }}$	$Rat{e}_{\mathrm{\infty }}$	$L_{\mathrm{\infty }}$	$Rat{e}_{\mathrm{\infty }}$
u	40	1.0428E-02		1.9061E-01		1.1751E-02		4.8382E-03
	80	4.3296E-03	1.2692	5.6042E-02	1.7660	2.6471E-03	2.1504	9.3388E-05	5.6951
	160	2.0103E-03	1.1057	1.4692E-02	1.9315	7.9880E-04	1.7285	8.2895E-06	3.4939
	320	1.0099E-03	0.9932	3.7840E-03	1.9571	2.2980E-04	1.7974	4.9402E-07	4.0687
v	40	1.1052E-03		4.6583E-03		4.1509E-03		4.8239E-04
	80	2.3443E-04	2.2371	1.3207E-03	1.8185	1.0695E-03	1.9565	1.3189E-05	5.1928
	160	4.5897E-05	2.3527	3.2872E-04	2.0064	2.7186E-04	1.9760	7.0235E-07	4.2311
	320	2.5551E-05	0.8450	7.4628E-05	2.1391	6.8697E-05	1.9846	3.9200E-08	4.1633
w	40	9.0058E-04		3.3970E-03		3.1735E-03		3.6448E-04
	80	1.7304E-04	2.3798	9.9194E-04	1.7759	7.8768E-04	2.0104	9.3151E-06	5.2901
	160	3.4840E-05	2.3123	2.6873E-04	1.8841	2.0158E-04	1.9497	6.1887E-07	3.9119
	320	1.9679E-05	0.8241	6.1197E-05	2.1346	5.1159E-05	1.9793	3.5131E-08	4.1388

Table 4. $L_2$ error and convergence rate of the proposed method for Example 5.1 Case 1 with $T=0.5$

	N	Ref. [12]		Method I [11]		Method II [11]		Present
		$L_2$	$Rat{e}_2$	$L_2$	$Rat{e}_2$	$L_2$	$Rat{e}_2$	$L_2$	$Rat{e}_2$
u	40	2.5253E-03		3.5289E-02		2.1007E-03		7.8299E-04
	80	1.2091E-03	1.0625	8.2914E-03	2.0895	5.2685E-04	1.9945	1.8822E-05	5.3785
	160	6.0491E-04	0.9992	1.5066E-03	2.4603	1.6415E-04	1.6823	1.2010E-06	3.9701
	320	3.0864E-04	0.9708	2.6537E-04	2.5052	4.5271E-05	1.8584	6.3541E-08	4.2404
v	40	3.3042E-04		1.1465E-03		1.0840E-03		8.6792E-05
	80	6.6619E-05	2.3103	2.8717E-04	1.9972	2.6672E-04	2.0230	2.4396E-06	5.1529
	160	1.5641E-05	2.0906	7.1723E-05	2.0014	6.6591E-05	2.0019	1.3130E-07	4.2157
	320	8.7707E-06	0.8346	1.6988E-05	2.0780	1.6623E-05	2.0021	7.7055E-09	4.0908
w	40	2.1003E-04		6.5958E-04		6.6917E-04		6.2106E-05
	80	4.0051E-05	2.3906	1.7702E-04	1.8976	1.5760E-04	2.0861	1.5154E-06	5.3570
	160	1.0450E-05	1.9383	4.5882E-05	1.9479	3.8783E-05	2.0227	8.9930E-08	4.0747
	320	5.9918E-06	0.8241	1.0342E-05	2.1519	9.6126E-06	2.0124	5.1941E-09	4.1138

Table 5. Error and convergence rate of the proposed method for Example 5.1 Case 4 with $T=0.5$

	N	Ref. [12]				Present
		$L_{\mathrm{\infty }}$	$Rat{e}_{\mathrm{\infty }}$	$L_2$	$Rat{e}_2$	$L_{\mathrm{\infty }}$	$Rat{e}_{\mathrm{\infty }}$	$L_2$	$Rat{e}_2$
u	40	8.746E-03		2.316E-03		4.1016E-03		6.7519E-04
	80	3.234E-03	1.4354	9.954E-04	1.2183	7.5162E-05	5.7700	1.7238E-05	5.2917
	160	1.793E-03	0.8510	5.089E-04	0.9679	7.7012E-06	3.2868	1.1261E-06	3.9361
	320	9.878E-04	0.8618	2.594E-04	0.9720	4.3944E-07	4.1313	5.7826E-08	4.2835
v	40	1.515E-03		4.158E-04		7.5836E-04		1.2686E-04
	80	4.899E-04	1.6290	1.322E-04	1.6530	1.6521E-05	5.5205	3.2876E-06	5.2701
	160	2.582E-04	0.9239	6.806E-05	0.9580	1.1134E-06	3.8913	1.9463E-07	4.0783
	320	1.373E-04	0.9111	3.562E-05	0.9342	6.4205E-08	4.1161	1.1096E-08	4.1326
w	40	1.000E-03		2.236E-04		3.5416E-04		6.0477E-05
	80	2.557E-04	1.9674	4.939E-05	2.1787	9.2230E-06	5.2630	1.5039E-06	5.3296
	160	7.122E-05	1.8441	1.337E-05	1.8850	6.0909E-07	3.9205	8.8829E-08	4.0815
	320	2.740E-05	1.3779	6.757E-06	0.9847	3.4308E-08	4.1501	5.1074E-09	4.1204

Table 6. Error and convergence rate of the proposed method for Example 5.1 Case 5 with $T=0.5$

	N	Ref. [12]				Present
		$L_{\mathrm{\infty }}$	$Rat{e}_{\mathrm{\infty }}$	$L_2$	$Rat{e}_2$	$L_{\mathrm{\infty }}$	$Rat{e}_{\mathrm{\infty }}$	$L_2$	$Rat{e}_2$
u	40	9.056E-03		2.357E-03		4.1718E-03		6.8510E-04
	80	3.389E-03	1.4179	1.016E-03	1.2147	7.6201E-05	5.7747	1.7322E-05	5.3056
	160	1.800E-03	0.9128	5.176E-04	0.9724	7.7494E-06	3.2977	1.1308E-06	3.9372
	320	1.153E-03	0.9105	2.642E-04	0.9710	4.4430E-07	4.1245	5.8266E-08	4.2785
v	40	1.339E-03		3.178E-04		7.1698E-04		1.2055E-04
	80	4.281E-04	1.6448	1.154E-04	1.6881	1.5900E-05	5.4949	3.1603E-06	5.2535
	160	2.347E-04	0.8672	6.034E-05	0.9353	1.0409E-06	3.9331	1.8527E-07	4.0924
	320	1.293E-04	0.8605	1.897E-05	0.9298	6.0021E-08	4.1162	1.0606E-08	4.1266
w	40	8.604E-04		2.078E-04		3.1558E-04		5.4683E-05
	80	2.144E-04	2.0061	4.726E-05	2.1362	8.5858E-06	5.1999	1.3951E-06	5.2927
	160	5.045E-05	2.0875	1.140E-05	2.0516	5.2775E-07	4.0240	7.9326E-08	4.1364
	320	1.466E-05	1.7829	4.089E-06	1.4792	2.9785E-08	4.1472	4.6058E-09	4.1063

On the other hand, we use the proposed numerical methods to solve the original system (1)-(3) directly. We compute the difference $E({\mathrm{\Theta }}_N)(=|{\mathrm{\Theta }}_{i,2N}-{\mathrm{\Theta }}_{i,N}|)$ between the numerical solutions obtained with the spatial mesh numbers N and 2N, as well as compute the relative error $E_{\nu }^r$ and the convergence $Rat{e}_{\nu }^r(={log}_2[E_{\nu }^r({\mathrm{\Theta }}_{N/2})/E_{\nu }^r({\mathrm{\Theta }}_N)])$ of the consecutive meshes. In which, $\mathrm{\Theta }\in \{u,v,w\}$, $\nu \in \{\mathrm{\infty },2\}$, ${\mathrm{\Theta }}_{i,N}$ stands for the numerical solution at $x_i$ with space mesh N at nth time step. With the same calculation parameters as above, the AFs $g_k=0$$(k=1,2,3)$ in the original system (1)-(3) with Cases 1, 4, and 5. Table 7 shows the convergence rates of the relative errors for solving the original problem (1)-(3) of Example 5.1 Cases 1, 4, and 5 with $T=0.5$, $\tau =h^2$, and compares these results with ones in [12]. It is evident that the convergence rates obtained by the method proposed in this paper is higher than those obtained in [12], which shows that our method is more advantageous than the method in [12].

Table 7. Convergence rates of the relative errors of the proposed method to solve the original problem for Example 5.1 Cases 1, 4, and 5

	N	Case 1				Case 4				Case 5
		Ref. [12]		Present		Ref. [12]		Present		Ref. [12]		Present
		$Rat{e}_{\mathrm{\infty }}^r$	$Rat{e}_2^r$	$Rat{e}_{\mathrm{\infty }}^r$	$Rat{e}_2^r$	$Rat{e}_{\mathrm{\infty }}^r$	$Rat{e}_2^r$	$Rat{e}_{\mathrm{\infty }}^r$	$Rat{e}_2^r$	$Rat{e}_{\mathrm{\infty }}^r$	$Rat{e}_2^r$	$Rat{e}_{\mathrm{\infty }}^r$	$Rat{e}_2^r$
u	40	1.2572	1.3348	2.7573	3.3310	1.2543	1.2854	2.6552	3.2388	0.8920	0.9621	2.6942	3.2715
	80	1.0495	1.1364	6.3984	6.4474	1.1381	1.0644	6.5801	6.6042	0.9432	1.0631	6.4924	6.5639
	160	1.0634	1.0054	4.8572	4.5994	0.9332	0.9456	4.8290	4.5723	0.9648	0.9823	4.9239	4.6229
v	40	1.2543	1.2932	4.1804	4.5701	1.2142	1.3256	5.5699	5.4263	1.1544	1.2832	5.5415	5.3842
	80	1.6836	1.4838	5.4257	5.0755	1.0073	1.0172	4.0486	4.1818	0.8321	1.1844	4.0225	4.1871
	160	0.8753	0.8645	3.9823	4.0416	0.8244	0.8401	3.9450	4.0910	0.8005	0.9187	3.9529	4.0885
w	40	1.7131	1.9730	4.9556	5.4212	1.8562	1.8783	4.9388	5.4042	1.5329	1.6621	4.6835	5.1581
	80	0.9567	1.1561	5.3866	5.6782	1.4374	1.3562	5.4346	5.7108	1.4233	1.6832	5.5952	5.8643
	160	0.8452	0.8534	3.2874	3.9148	1.1253	0.9734	3.2587	3.8995	1.3852	1.3562	3.2866	3.9168

Accuracy test with PPL for 1D problem

Next, to test our methods still have high accuracy after adding the designed PPL or bound-preserving limiter, similar to Example 5.1, we give the analytical solution in the following form, $u(x,t)=v(x,t)=w(x,t)=exp(-t)cosx,0⩽x⩽2\pi ,t⩾0.$

Now we consider two cases of the haptotaxis system (1)-(3), i.e., Cases 2 and 3. The AFs $g_k$$(k=1,2,3)$ are listed respectively in the following Table 8 for Cases 2 and 3 on account of the analytical solution above.

Table 8. Expressions of $g_k(k=1,2,3)$ for Example 5.2 Cases 2 and 3

	AFs	Expressions	AFs
Case 2	$g_3$	$(-1+D_2+\beta )exp(-t)cosx-\alpha exp(-2t){cos}^{2}x$
	$g_1$	$(-1+D_1-{\mu }_1)exp(-t)cosx+(2{\mu }_1{cos}^{2}x-\chi cos2x)exp(-2t)$	$g_1$	Case 3
	$g_2$	−exp(−t)cosx + ηexp(−2t)cos2x	$g_2$
		$(-1+D_2-\alpha +\beta )exp(-t)cosx+\alpha exp(-2t){cos}^{2}x$	$g_3$

In this computation, we take $\alpha =\beta =\eta =1$, $D_1=0.05$, $D_2=0.1$, $\chi =0.001$, ${\mu }_1=0.5$. The $L_{\mathrm{\infty }}$, $L_2$ norm errors and convergence rates, which obtained separately to solve Example 5.2 Cases 2 and 3 with $\tau =0.5h^{4/3}$, $T=1$, are listed in Tables 9 and 10. Observing the data in these tables, we find that although the error is slightly increased after adding the PPL, it still has fourth-order convergence rate. It shows that the positivity-preserving algorithm designed in this paper maintains the accuracy of the proposed difference scheme, and is consistent with the previous analysis.

Table 9. Error and convergence rate of the proposed method for Example 5.2 Case 2 with $T=1$

	N	Without PPL				With PPL
		$L_{\mathrm{\infty }}$	$Rat{e}_{\mathrm{\infty }}$	$L_2$	$Rat{e}_2$	$L_{\mathrm{\infty }}$	$Rat{e}_{\mathrm{\infty }}$	$L_2$	$Rat{e}_2$
u	20	4.2675E-04		4.4860E-04		4.2547E-04		4.4952E-04
	40	3.0462E-05	3.8083	3.2092E-05	3.8052	3.0481E-05	3.8031	3.2104E-05	3.8076
	80	1.9125E-06	3.9935	2.0138E-06	3.9942	1.9129E-06	3.9941	2.0140E-06	3.9946
	160	1.2197E-07	3.9709	1.2843E-07	3.9709	1.2198E-07	3.9710	1.2843E-07	3.9710
	320	7.6579E-09	3.9935	8.0627E-09	3.9935	7.6580E-09	3.9935	8.0627E-09	3.9936
v	20	8.9431E-05		1.0813E-04		8.9350E-05		1.0815E-04
	40	6.0367E-06	3.8889	7.1729E-06	3.9140	6.0393E-06	3.8870	7.1737E-06	3.9142
	80	3.7538E-07	4.0074	4.4936E-07	3.9966	3.7542E-07	4.0078	4.4937E-07	3.9967
	160	2.3745E-08	3.9827	2.8454E-08	3.9811	2.3746E-08	3.9828	2.8454E-08	3.9812
	320	1.4872E-09	3.9969	1.7853E-09	3.9944	1.4872E-09	3.9970	1.7853E-09	3.9944
w	20	5.0742E-05		6.8433E-05		5.1492E-05		6.8442E-05
	40	3.5762E-06	3.8267	4.4358E-06	3.9474	3.5843E-06	3.8446	4.4425E-06	3.9454
	80	2.3382E-07	3.9350	2.9205E-07	3.9249	2.3399E-07	3.9372	2.9220E-07	3.9263
	160	1.5031E-08	3.9594	1.8605E-08	3.9725	1.5035E-08	3.9600	1.8609E-08	3.9729
	320	9.4789E-10	3.9871	1.1724E-09	3.9881	9.4800E-10	3.9873	1.1725E-09	3.9883

Table 10. Error and convergence rate of the proposed method for Example 5.2 Case 3 with $T=1$

	N	Without PPL				With PPL
		$L_{\mathrm{\infty }}$	$Rat{e}_{\mathrm{\infty }}$	$L_2$	$Rat{e}_2$	$L_{\mathrm{\infty }}$	$Rat{e}_{\mathrm{\infty }}$	$L_2$	$Rat{e}_2$
u	16	1.2096E-03		1.2483E-03		1.2131E-03		1.2510E-03
	32	8.7114E-05	3.7955	9.0076E-05	3.7927	8.7188E-05	3.7984	9.0119E-05	3.7952
	64	5.5078E-06	3.9833	5.6900E-06	3.9846	5.5091E-06	3.9842	5.6907E-06	3.9852
	128	3.5075E-07	3.9730	3.6220E-07	3.9735	3.5077E-07	3.9732	3.6222E-07	3.9737
	256	2.2099E-08	3.9884	2.2817E-08	3.9886	2.2099E-08	3.9885	2.2817E-08	3.9887
v	16	6.4822E-04		6.3553E-04		6.4903E-04		6.3559E-04
	32	4.7349E-05	3.7751	4.6119E-05	3.7845	4.7367E-05	3.7763	4.6119E-05	3.7847
	64	3.0273E-06	3.9672	2.9490E-06	3.9670	3.0276E-06	3.9676	2.9491E-06	3.9670
	128	1.9312E-07	3.9705	1.8803E-07	3.9712	1.9313E-07	3.9706	1.8803E-07	3.9712
	256	1.2174E-08	3.9876	1.1851E-08	3.9878	1.2174E-08	3.9877	1.1851E-08	3.9878
w	16	8.8553E-04		9.3469E-04		8.8830E-04		9.3668E-04
	32	6.3672E-05	3.7978	6.7583E-05	3.7898	6.3741E-05	3.8008	6.7623E-05	3.7920
	64	4.0337E-06	3.9805	4.2800E-06	3.9810	4.0349E-06	3.9816	4.2806E-06	3.9816
	128	2.5693E-07	3.9727	2.7250E-07	3.9733	2.5695E-07	3.9730	2.7251E-07	3.9735
	256	1.6188E-08	3.9883	1.7167E-08	3.9885	1.6188E-08	3.9884	1.7167E-08	3.9886

Non-negative test and evolutionary solution over time for 1D problem

Next, we aim to use the established compact difference scheme and algorithm to verify the non-negativity, and to simulate and analyze the temporal evolutionarily phenomena of the solution in the complex process of TC invades into the surrounding tissues.

In this computation, we take the same parameter values with [12], i.e., $\alpha =0.1$, $\beta =0.1$, $\eta =1$, $\epsilon =0.01$, $D_1=D_2=0.001$, $\chi =0.01$, and $\tau =h^2$, a fixed space mesh number $N=160$. Figure 3 displays respectively that the time-varying curves of TC, ECM and MDE obtained by the method proposed in this paper without and with using the PPL to simulate system (1)-(3), where $T=1$. From Fig. 3, noticing that TC and MDE appeared small negative values in some local areas without using PPL, and the density of ECM does not generate negative values. After using the PPL, the local negative values are effectively kicked out, which is also the original intention of designing the positivity-preserving algorithm. The experiment shows that the positivity-preserving algorithm designed in this paper cannot only ensure high resolution, but also have a good non-negative guarantee property.

[See PDF for image]

Figure 3

Comparisons of TC, ECM and MDE obtained by the proposed method with/without PPL

On the other hand, we take $\tau =h^{4/3}$, and use the constructed scheme with PPL to solve the tumor invasion haptotaxis model (1)-(3), as well as to simulate the change of TC, ECM and MDE over time. Fig. 4 displays respectively that the space-time evolution for TC, ECM and MDE at $T=5$. It can be seen that TC, ECM and MDE interact in the tumor microenvironment and undergo subtle changes. TC adhere to the surrounding ECM, produce and secrete MDE, and gradually destroy the surrounding healthy tissues, moving in the gradient direction of the ECM. Exactly as Fig. 4, at the $x=0$ boundary areas, the initial value of TC continues to spread with the spatial region, and gradually migrates to the degraded ECM region over time. This is consistent with the conclusion obtained by simulation in [12]. In addition, the solution of the haptotaxis system (1)-(3) with time tends to a constant steady-state solution $(0,1,0)$, which also verifies the conclusion obtained in theoretical research, which exists the unique globally uniformly bounded solution and is no finite-time blow-up phenomenon in the 1D space of the haptotaxis system (1)-(3) for tumor invasion [4, 6, 14, 18].

[See PDF for image]

Figure 4

Space-time evolution for TC, ECM and MDE

Accuracy test without PPL for 2D problem

Next, we carry out the accuracy test for the 2D problems of the haptotaxis system (1)-(3). Under in the absence of PPL, we add the AFs $g_k$$(k=1,2,3)$ with the following analytical solution [13] to the right end of each equation for the original system (1)-(3) in 2D space, i.e., $\begin{array}{rlr}& u(x,y,t)=exp(-t)(\frac{1}{6}cos(2\pi x)+\frac{1}{6}cos(2\pi y)+\frac{1}{2}),& 0⩽x,y⩽1,t⩾0,\\ & v(x,y,t)=1-u(x,y,t),& 0⩽x,y⩽1,t⩾0,\\ & w(x,y,t)=exp(-t)(\frac{1}{2}cos(2\pi x)-\frac{1}{2}sin(2\pi y)+\pi y+1),& 0⩽x,y⩽1,t⩾0.\end{array}$$u(x,y,0)$, $v(x,y,0)$ and $w(x,y,0)$ are derived from the analytical solution above, and the boundary condition (14e) is used, but note that $\frac{\partial u}{\partial \mathbf{n}}=\frac{\partial v}{\partial \mathbf{n}}=\frac{\partial w}{\partial \mathbf{n}}=0$. Table 11 listed respectively that the AFs $g_k$$(k=1,2,3)$ for the system (1)-(3) in 2D space under Cases 2 and 4.

Table 11. Expressions of $g_k(k=1,2,3)$ for Example 5.4 Cases 2 and 4

	AFs	Expressions	AFs
Case 2	$g_3$	$\frac{1}{36}exp(-t)\{18(\beta -1+4D_2{\pi }^2)(cos(2\pi x)-sin(2\pi y))+\alpha (cos(2\pi x)+cos(2\pi y)+3\left)\right[exp(-t)(cos(2\pi x)+cos(2\pi y)+3)-6]+36(\beta -1)(\pi y+1)\}$
	$g_2$	$\frac{1}{12}exp(-t)\left\{\right[2-\eta exp(-t)(cos(2\pi x)-sin(2\pi y)+2\pi y+2)\left]\right(cos(2\pi x)+cos(2\pi y))+3\eta (2-exp(-t))(cos(2\pi x)-sin(2\pi y)+2\pi y+2)+6\}$	$g_2$	Case 4
	$g_1$	$\frac{1}{18}exp(-t)\{(12D_1{\pi }^2+6\chi {\pi }^{2}exp(-t)-3)[cos(2\pi x)+cos(2\pi y)]+2\chi {\pi }^{2}exp(-t)[cos(4\pi x)+cos(4\pi y)+2cos(2\pi x)cos(2\pi y)]-9\}$	$g_1$
		$\frac{1}{6}exp(-t)\{3(\beta -1+4{\pi }^2{D}_2)[cos(2\pi x)-sin(2\pi y)]-\alpha [cos(2\pi x)+cos(2\pi y)+3]+6(\beta -1)(\pi y+1)\}$	$g_3$

In this computation, we preset $\alpha =\beta =\eta =1$, $\chi =0.001$, ${\mu }_1=0.5$, ${\mu }_2=0.1$, $a=c=0$, and $b=d=1$. Table 12 displays that the errors, convergence rate and CPU time of the proposed method in time for Example 5.4 Cases 2 and 4, where $D_1=0.0002$, $D_2=0.0001$, $T=1$, $h_x=h_y=200$. Observe the Table 12, we find that our methods proposed can obtain third-order accuracy in the temporal direction for both cases. Table 13 shows that the errors, convergence rates and CPU time of the proposed method in space for Example 5.4 Cases 2 and 4, in which, $D_1=0.05$, $D_2=0.01$, $T=1$, $\tau =2\times {10}^{-4}$, which shows that our method has fourth-order accuracy in space. To test the spatial-temporal convergence rate of our methods for solving Example 5.4, we take the spatial-temporal mesh ratio $\tau /h_x^{4/3}=\tau /h_y^{4/3}=0.1$ for this computation on the equidistant grid in the x− and y− directions, and list the results obtained for Cases 2 and 4 in Table 14, where $T=1$. Meanwhile, we take $\tau =min\{0.1h_x^{4/3},0.2h_y^{4/3}\}$ for this computation on the non-equidistant grid in the x− and y− directions, and exhibit the results obtained for Cases 2 and 4 in Table 15, where $T=1$. The above numerical experiments show that our numerical methods can obtain the fourth-order convergence rates regardless of the equidistant or non-equidistant grids in space.

Table 12. Errors, convergence rates and CPU time of the proposed method in time for Example 5.4 Cases 2 and 4 with $T=1$, $h_x=h_y=200$

	M	Cases 2					Cases 4
		$L_{\mathrm{\infty }}$	$Rat{e}_{\mathrm{\infty }}$	$L_2$	$Rat{e}_2$	CPU time	$L_{\mathrm{\infty }}$	$Rat{e}_{\mathrm{\infty }}$	$L_2$	$Rat{e}_2$	CPU time
u	40	7.7958E-08		5.0289E-08		8.915	8.4017E-08		5.5302E-08		8.690
	50	3.9461E-08	3.0513	2.5470E-08	3.0486	10.792	4.2579E-08	3.0458	2.8046E-08	3.0427	11.065
	60	2.2640E-08	3.0474	1.4633E-08	3.0399	13.215	2.4451E-08	3.0425	1.6127E-08	3.0351	13.034
	70	1.4156E-08	3.0463	9.1670E-09	3.0338	14.471	1.5298E-08	3.0420	1.0109E-08	3.0299	15.241
	80	9.4205E-09	3.0496	6.1171E-09	3.0295	17.512	1.0186E-08	3.0455	6.7486E-09	3.0262	16.937
v	40	1.4885E-06		7.3242E-07		8.915	1.4616E-06		7.1469E-07		8.690
	50	7.5642E-07	3.0336	3.7155E-07	3.0414	10.792	7.4266E-07	3.0341	3.6255E-07	3.0416	11.065
	60	4.3556E-07	3.0274	2.1370E-07	3.0337	13.215	4.2760E-07	3.0279	2.0851E-07	3.0339	13.034
	70	2.7331E-07	3.0232	1.3399E-07	3.0285	14.471	2.6830E-07	3.0235	1.3073E-07	3.0286	15.241
	80	1.8260E-07	3.0201	8.9464E-08	3.0247	17.512	1.7925E-07	3.0204	8.7290E-08	3.0248	16.937
w	40	1.1984E-06		6.8781E-07		8.915	9.3458E-07		6.0403E-07		8.690
	50	6.1172E-07	3.0137	3.5116E-07	3.0127	10.792	4.7645E-07	3.0193	3.0799E-07	3.0185	11.065
	60	3.5328E-07	3.0112	2.0283E-07	3.0104	13.215	2.7492E-07	3.0159	1.7775E-07	3.0151	13.034
	70	2.2215E-07	3.0094	1.2756E-07	3.0088	14.471	1.7276E-07	3.0137	1.1171E-07	3.0128	15.241
	80	1.4866E-07	3.0081	8.5365E-08	3.0077	17.512	1.1555E-07	3.0123	7.4728E-08	3.0111	16.937

Table 13. Errors, convergence rates and CPU time of the proposed method in space for Example 5.4 Cases 2 and 4 with $T=1$, $\tau =2\times {10}^{-4}$

	$N_x\times N_y$	Cases 2					Cases 4
		$L_{\mathrm{\infty }}$	$Rat{e}_{\mathrm{\infty }}$	$L_2$	$Rat{e}_2$	CPU time	$L_{\mathrm{\infty }}$	$Rat{e}_{\mathrm{\infty }}$	$L_2$	$Rat{e}_2$	CPU time
u	8 × 8	9.4069E-04		3.8795E-04		1.885	9.4176E-04		3.8863E-04		1.986
	16 × 16	1.1295E-04	3.0580	4.2786E-05	3.1807	6.850	1.1321E-04	3.0564	4.2882E-05	3.1800	7.001
	32 × 32	5.5611E-06	4.3442	1.9614E-06	4.4472	25.409	5.5742E-06	4.3441	1.9659E-06	4.4471	26.954
	64 × 64	2.5715E-07	4.4347	9.3630E-08	4.3888	102.592	2.5777E-07	4.4346	9.3835E-08	4.3890	109.772
	128 × 128	1.2792E-08	4.3293	5.1566E-09	4.1825	419.716	1.2823E-08	4.3293	5.1671E-09	4.1827	442.0740
v	8 × 8	4.7109E-04		1.6134E-04		1.885	4.8053E-04		1.7379E-04		1.986
	16 × 16	2.5113E-05	4.2295	9.6246E-06	4.0672	6.850	3.0271E-05	3.9886	1.1488E-05	3.9191	7.001
	32 × 32	1.4939E-06	4.0713	5.2310E-07	4.2016	25.409	1.7466E-06	4.1153	6.1679E-07	4.2192	26.954
	64 × 64	8.7728E-08	4.0899	3.0648E-08	4.0932	102.592	9.9808E-08	4.1292	3.5380E-08	4.1238	109.772
	128 × 128	5.2770E-09	4.0552	1.8751E-09	4.0308	419.716	5.9088E-09	4.0782	2.1437E-09	4.0447	442.074
w	8 × 8	9.3848E-04		4.0295E-04		1.885	1.0274E-03		4.3550E-04		1.986
	16 × 16	8.5876E-05	3.4500	2.9191E-05	3.7870	6.850	1.0580E-04	3.2797	3.4064E-05	3.6763	7.001
	32 × 32	4.1842E-06	4.3592	1.5113E-06	4.2716	25.409	5.0308E-06	4.3943	1.7566E-06	4.2773	26.954
	64 × 64	2.0859E-07	4.3262	8.3785E-08	4.1730	102.592	2.4650E-07	4.3511	9.5805E-08	4.1966	109.772
	128 × 128	1.1107E-08	4.2311	4.9905E-09	4.0694	419.716	1.2980E-08	4.2473	5.6469E-09	4.0846	442.074

Table 14. Errors and convergence rates of the proposed method for Example 5.4 Cases 2 and 4 with $T=1$, $\tau =0.1h_x^{4/3}=0.1h_y^{4/3}$

	$N_x\times N_y$	Cases 2				Cases 4
		$L_{\mathrm{\infty }}$	$Rat{e}_{\mathrm{\infty }}$	$L_2$	$Rat{e}_2$	$L_{\mathrm{\infty }}$	$Rat{e}_{\mathrm{\infty }}$	$L_2$	$Rat{e}_2$
u	10 × 10	5.9323E-04		2.5227E-04		5.9435E-04		2.5279E-04
	20 × 20	4.4052E-05	3.7513	1.6153E-05	3.9651	4.4155E-05	3.7507	1.6190E-05	3.9648
	40 × 40	2.0588E-06	4.4193	7.2504E-07	4.4776	2.0637E-06	4.4192	7.2669E-07	4.4776
	80 × 80	9.6922E-08	4.4088	3.6332E-08	4.3188	9.7155E-08	4.4088	3.6409E-08	4.3190
	160 × 160	4.9600E-09	4.2884	2.0695E-09	4.1339	4.9722E-09	4.2883	2.0737E-09	4.1341
v	10 × 10	2.0816E-04		6.6956E-05		2.0829E-04		7.5996E-05
	20 × 20	1.0299E-05	4.3372	3.7482E-06	4.1589	1.2319E-05	4.0797	4.4807E-06	4.0841
	40 × 40	6.0011E-07	4.1011	2.0827E-07	4.1697	6.9249E-07	4.1529	2.4361E-07	4.2010
	80 × 80	3.5450E-08	4.0814	1.2432E-08	4.0663	3.9985E-08	4.1143	1.4292E-08	4.0913
	160 × 160	2.1438E-09	4.0475	7.6564E-10	4.0213	2.3925E-09	4.0628	8.7395E-10	4.0315
w	10 × 10	5.3931E-04		1.8939E-04		6.2987E-04		2.1380E-04
	20 × 20	3.2315E-05	4.0608	1.1284E-05	4.0690	4.0175E-05	3.9707	1.3200E-05	4.0176
	40 × 40	1.5851E-06	4.3495	5.8927E-07	4.2592	1.8961E-06	4.4052	6.8123E-07	4.2763
	80 × 80	8.0462E-08	4.3002	3.3608E-08	4.1320	9.4770E-08	4.3225	3.8270E-08	4.1539
	160 × 160	4.4599E-09	4.1732	2.0283E-09	4.0505	5.1561E-09	4.2001	2.2904E-09	4.0626

Table 15. Errors and convergence rates of the proposed method for Example 5.4 Cases 2 and 4 with $T=1$, $\tau =min\{0.1h_x^{4/3},0.2h_y^{4/3}\}$

	$N_x\times N_y$	Cases 2				Cases 4
		$L_{\mathrm{\infty }}$	$Rat{e}_{\mathrm{\infty }}$	$L_2$	$Rat{e}_2$	$L_{\mathrm{\infty }}$	$Rat{e}_{\mathrm{\infty }}$	$L_2$	$Rat{e}_2$
u	5 × 8	3.7679E-03		1.9717E-03		3.7723E-03		1.9742E-03
	10 × 16	3.5202E-04	3.4200	1.8489E-04	3.4147	3.5308E-04	3.4174	1.8526E-04	3.4136
	20 × 32	2.4859E-05	3.8239	1.1541E-05	4.0017	2.4936E-05	3.8237	1.1567E-05	4.0015
	40 × 64	1.1618E-06	4.4192	5.1721E-07	4.4799	1.1655E-06	4.4192	5.1833E-07	4.4799
	80 × 128	5.5050E-08	4.3995	2.6034E-08	4.3123	5.5229E-08	4.3994	2.6087E-08	4.3124
v	5 × 8	1.7675E-03		9.6389E-04		1.9886E-03		1.0833E-03
	10 × 16	1.2603E-04	3.8099	4.8452E-05	4.3142	1.4956E-04	3.7329	5.9071E-05	4.1969
	20 × 32	6.8896E-06	4.1932	3.1205E-06	3.9567	8.1467E-06	4.1984	3.8213E-06	3.9503
	40 × 64	3.8498E-07	4.1616	1.6988E-07	4.1992	4.4779E-07	4.1853	2.0391E-07	4.2281
	80 × 128	2.2422E-08	4.1018	9.8570E-09	4.1072	2.5734E-08	4.1211	1.1689E-08	4.1247
w	5 × 8	3.7602E-03		2.3562E-03		4.0946E-03		2.5435E-03
	10 × 16	3.0104E-04	3.6428	1.4604E-04	4.0121	3.2575E-04	3.6519	1.6839E-04	3.9169
	20 × 32	2.1488E-05	3.8084	9.6729E-06	3.9162	2.5352E-05	3.6836	1.1250E-05	3.9038
	40 × 64	1.1179E-06	4.2646	5.0316E-07	4.2649	1.3088E-06	4.2758	5.8017E-07	4.2773
	80 × 128	5.6197E-08	4.3142	2.8186E-08	4.1579	6.5691E-08	4.3164	3.2231E-08	4.1699

Accuracy test with PPL for 2D problem

Now, we consider three cases of the haptotaxis system (1)-(3) for tumor invasion, i.e., Cases 1, 3 and 5 in the following form, under the designed PPL or bound-preserving limiter, $u(x,y,t)=v(x,y,t)=w(x,y,t)=exp(-2t)cos(2x)cos(2y),a⩽x⩽b,c⩽y⩽d,t⩾0.$$u(x,y,0)$, $v(x,y,0)$ and $w(x,y,0)$ are derived from the analytical solution above, and the boundary condition (14e) is used, and $\frac{\partial u}{\partial \mathbf{n}}=\frac{\partial v}{\partial \mathbf{n}}=\frac{\partial w}{\partial \mathbf{n}}=0$. In which, $\alpha =0.25$, $\beta =0.5$, ${\mu }_1=0.75$ and ${\mu }_2=1$. $g_k(k=1,2,3)$ see Table 16 for details.

Table 16. Expressions of $g_k(k=1,2,3)$ for Example 5.5 Cases 1, 3 and 5

	AFs	Expressions	AFs
Case 5	$g_2$	$(-2-{\mu }_2)exp(-2t)cos(2x)cos(2y)+(\eta +2{\mu }_2)exp(-4t){cos}^2(2x){cos}^2(2y)$
	$g_3$	$(8D_2-2-\alpha +\beta )exp(-2t)cos(2x)cos(2y)+\alpha exp(-4t){cos}^2(2x){cos}^2(2y)$	$g_3$	Case 3
	$g_1$	$(8D_1-2-{\mu }_1)exp(-2t)cos(2x)cos(2y)+2{\mu }_{1}exp(-4t){cos}^2(2x){cos}^2(2y)-4\chi exp(-4t)(cos(4x){cos}^2(2y)+{cos}^2(2x)cos(4y))$	$g_1$
Case 1	$g_2$	(ηexp(−2t)cos(2x)cos(2y)−2)exp(−2t)cos(2x)cos(2y)	$g_2$
	$g_1$	$(8D_1-2)exp(-2t)cos(2x)cos(2y)-4\chi exp(-4t)(cos(4x){cos}^2(2y)+{cos}^2(2x)cos(4y))$
	$g_3$	$(8D_2-2-\alpha +\beta )exp(-2t)cos(2x)cos(2y)$

First, we take $\eta =1$, $\chi =0.001$, $D_1=0.0001$ and $D_2=0.0002$, and utilize the proposed scheme with PPL to solve the Cases 1 and 5 for Example 5.5 in $[a,b]\times [c,d]=[0,\pi ]\times [0,\pi ]$. Table 17 lists that the errors, convergence rates and CPU time obtained by using our methods to solve the Case 5 for Example 5.5, where $T=1$, $\tau =0.5h_x^{4/3}=0.5h_y^{4/3}$, and the equidistant meshes in x− and y− directions are used. The results obtained by using the constructed scheme with PPL and without PPL for solving the Cases 1 of Example 5.5 are listed Table 18, in which, $T=1$, $\tau =min\{0.5h_x^{4/3},0.3h_y^{4/3}\}$, and the non-equidistant meshes in x− and y-directions are used. Second, to test the performance of the proposed method more abundantly, we take $\eta =5$, $\chi =0.005$, $D_1=0.0005$ and $D_2=0.0002$, and use the proposed method without and with PPL on the non-equidistant meshes in x− and y-directions for solving Example 5.5 Cases 3 for in $[a,b]\times [c,d]=[-\frac{\pi }{2},\frac{\pi }{2}]\times [0,\pi ]$. The computed results are shown in Table 19, where $T=1$, $\tau =min\{0.5h_x^{4/3},0.3h_y^{4/3}\}$. By observing the data in these tables, it is not difficult to find that although the computational errors obtained by using PPL is slightly increased, the numerical method proposed in this paper can still converge at the fourth-order rate after using PPL.

Table 17. Errors, convergence rates and CPU time of the proposed method for Example 5.5 Cases 5 with $T=1$, $\tau =0.5h_x^{4/3}=0.5h_y^{4/3}$

	$N_x\times N_y$	Without PPL			With PPL
		$L_2$	$Rat{e}_2$	CPU time	$L_2$	$Rat{e}_2$	CPU time
u	32 × 32	5.0426169E-05		0.291	5.0427794E-05		0.300
	40 × 40	2.1373612E-05	3.8466	0.625	2.1373917E-05	3.8467	0.631
	64 × 64	3.4638031E-06	3.8719	3.041	3.4638129E-06	3.8719	3.040
	80 × 80	1.4486893E-06	3.9065	6.722	1.4486914E-06	3.9065	6.388
	128 × 128	2.2541664E-07	3.9584	31.398	2.2541672E-07	3.9584	32.747
v	32 × 32	1.5061632E-04		0.291	1.5061797E-04		0.300
	40 × 40	6.4024154E-05	3.8337	0.625	6.4024376E-05	3.8338	0.631
	64 × 64	1.0453656E-05	3.8559	3.041	1.0453659E-05	3.8560	3.040
	80 × 80	4.3849862E-06	3.8933	6.722	4.3849866E-06	3.8933	6.388
	128 × 128	6.8403557E-07	3.9530	31.398	6.8403558E-07	3.9530	32.747
w	32 × 32	6.7276900E-06		0.291	6.7279747E-06		0.300
	40 × 40	2.8475382E-06	3.8530	0.625	2.8475909E-06	3.8531	0.631
	64 × 64	4.6228157E-07	3.8681	3.041	4.6228326E-07	3.8682	3.040
	80 × 80	1.9357159E-07	3.9012	6.722	1.9357194E-07	3.9012	6.388
	128 × 128	3.0093663E-08	3.9603	31.398	3.0093678E-08	3.9603	32.747

Table 18. Errors, convergence rates and CPU time of the proposed method for Example 5.5 Cases 1 with $T=1$, $\tau =min\{0.5h_x^{4/3},0.3h_y^{4/3}\}$

	$N_x\times N_y$	Without PPL					With PPL
		$L_{\mathrm{\infty }}$	$Rat{e}_{\mathrm{\infty }}$	$L_2$	$Rat{e}_2$	CPU time	$L_{\mathrm{\infty }}$	$Rat{e}_{\mathrm{\infty }}$	$L_2$	$Rat{e}_2$	CPU time
u	20 × 32	2.5462E-06		1.8833E-06		0.310	3.0012E-03		1.6197E-03		0.325
	40 × 64	3.4886E-07	2.8676	1.6143E-07	3.5443	3.004	3.4853E-04	3.1062	1.4289E-04	3.5028	3.270
	50 × 80	1.5262E-07	3.7049	6.4078E-08	4.1407	6.411	1.5631E-04	3.5937	7.0323E-05	3.1771	7.136
	80 × 128	2.4256E-08	3.9134	8.5719E-09	4.2800	31.455	2.4200E-05	3.9691	7.4167E-06	4.7859	34.565
	100 × 160	9.7733E-09	4.0737	3.2466E-09	4.3509	67.214	9.2496E-06	4.3101	2.5877E-06	4.7187	73.327
v	20 × 32	2.9863E-07		2.2964E-07		0.310	3.0459E-03		1.6396E-03		0.325
	40 × 64	2.3769E-08	3.6512	1.5694E-08	3.8711	3.004	3.6227E-04	3.0718	1.4738E-04	3.4757	3.270
	50 × 80	9.9983E-09	3.8807	6.4736E-09	3.9684	6.411	1.6571E-04	3.5050	7.2586E-05	3.1741	7.136
	80 × 128	1.5475E-09	3.9697	9.9209E-10	3.9908	31.455	2.8077E-05	3.7772	8.4152E-06	4.5845	34.565
	100 × 160	6.2895E-10	4.0349	4.0638E-10	3.9998	67.214	1.1618E-05	3.9544	3.1406E-06	4.4170	73.327
w	20 × 32	6.3652E-07		3.7663E-07		0.310	2.8648E-03		1.5446E-03		0.325
	40 × 64	6.8654E-08	3.2128	2.7545E-08	3.7733	3.004	3.1714E-04	3.1752	1.3016E-04	3.5689	3.270
	50 × 80	2.9430E-08	3.7960	1.0854E-08	4.1734	6.411	1.3848E-04	3.7133	6.4874E-05	3.1204	7.136
	80 × 128	4.4478E-09	4.0204	1.4286E-09	4.3145	31.455	1.9786E-05	4.1399	6.1271E-06	5.0207	34.565
	100 × 160	1.7338E-09	4.2220	5.3780E-10	4.3782	67.214	7.1862E-06	4.5388	2.0517E-06	4.9029	73.327

Table 19. Errors, convergence rates and CPU time of the proposed method for Example 5.5 Cases 3 with $T=1$, $\tau =min\{0.5h_x^{4/3},0.3h_y^{4/3}\}$

	$N_x\times N_y$	Without PPL					With PPL
		$L_{\mathrm{\infty }}$	$Rat{e}_{\mathrm{\infty }}$	$L_2$	$Rat{e}_2$	CPU time	$L_{\mathrm{\infty }}$	$Rat{e}_{\mathrm{\infty }}$	$L_2$	$Rat{e}_2$	CPU time
u	20 × 32	5.3059E-05		2.2165E-05		0.311	3.5610E-03		1.9527E-03		0.334
	40 × 64	6.5666E-06	3.0144	1.8955E-06	3.5476	3.117	4.4692E-04	2.9942	1.8854E-04	3.3725	3.101
	50 × 80	2.7850E-06	3.8439	7.3088E-07	4.2707	6.582	1.9886E-04	3.6289	9.7470E-05	2.9566	6.654
	80 × 128	3.7878E-07	4.2447	8.5802E-08	4.5579	31.832	2.5151E-05	4.3994	8.6731E-06	5.1475	33.673
	100 × 160	1.3177E-07	4.7317	2.9555E-08	4.7763	66.690	8.0507E-06	5.1049	2.7162E-06	5.2029	69.677
v	20 × 32	3.5291E-05		1.8104E-05		0.311	3.2121E-03		1.7227E-03		0.334
	40 × 64	2.8379E-06	3.6364	1.2656E-06	3.8385	3.117	3.8993E-04	3.0422	1.5706E-04	3.4552	3.101
	50 × 80	1.2185E-06	3.7890	5.2672E-07	3.9284	6.582	1.7872E-04	3.4962	7.6094E-05	3.2476	6.654
	80 × 128	1.9224E-07	3.9289	8.1441E-08	3.9719	31.832	2.9980E-05	3.7984	8.8943E-06	4.5671	33.673
	100 × 160	7.8497E-08	4.0139	3.3471E-08	3.9849	66.690	1.2286E-05	3.9977	3.2922E-06	4.4539	69.677
w	20 × 32	7.1914E-06		2.9322E-06		0.311	2.9140E-03		1.5724E-03		0.334
	40 × 64	9.1536E-07	2.9739	2.5654E-07	3.5147	3.117	3.2703E-04	3.1555	1.3438E-04	3.5486	3.101
	50 × 80	3.9377E-07	3.7803	9.9922E-08	4.2256	6.582	1.4267E-04	3.7174	6.6644E-05	3.1429	6.654
	80 × 128	5.6392E-08	4.1350	1.2007E-08	4.5083	31.832	1.9866E-05	4.1947	6.2049E-06	5.0511	33.673
	100 × 160	2.0626E-08	4.5073	4.1919E-09	4.7159	66.690	7.0718E-06	4.6287	2.0462E-06	4.9714	69.677

Non-negative test and evolution dynamics simulation of solution for 2D problem

Case 1 [9]: $\begin{array}{rl}& u(x,y,0)=exp\left(\frac{{(x-0.5)}^2+{(y-0.5)}^2}{ϵ}\right),0⩽x,y⩽1,\\ & v(x,y,0)=\{\begin{array}{ll}1-0.5u(x,y,0),& x^2+y^2⩽0.25;\\ 1-0.5|sin(\pi (x-0.5))|-0.5exp(-\frac{0.25}{ϵ}),& otherwise\end{array}\\ & w(x,y,0)=0.5u(x,y,0),0⩽x,y⩽1.\end{array}$

Case 2 [13]: $\begin{array}{rlr}& u(x,y,0)=exp(-400{(x-0.5)}^2-400{(y-0.5)}^2),& 0⩽x,y⩽1,\\ & v(x,y,0)=1-u(x,y,0),& 0⩽x,y⩽1,\\ & w(x,y,0)=0.5u(x,y,0),& 0⩽x,y⩽1.\end{array}$

In this computation, firstly, we give parameters $D_1=D_2=0.001$, $\alpha =0.1$, $\beta =0$, $\eta =10$, $\chi =0.005$ [9, 13], where $ϵ=0.0025$, $\tau =h_x^{4/3}=h_y^{4/3}$ for Case 1, and $\tau =min\{0.5h_x^{4/3},0.3h_x^{4/3}\}$, ${\mu }_1=2$ for Case 2. Then, the proposed method at the absence of PPL and at the present of PPL is used to solve and simulate the haptotaxis system (1)-(3) Cases 1 and 2. Figures 5 and 6 show the comparison of the numerical solutions for TC, ECM and MDE obtained by the proposed method with and without PPL to solve the Cases 1 and 2 for Example 5.6, where $h_x=h_y=1/20$, $T=1$. Figures 7 and 8 show the 1D profiles of the numerical solutions of TC and MDE along $x=0.5$, respectively. By observe, we find that the numerical solutions of TC and MDE have small negative values in the local region without using PPL, and these negative values are effectively kicked out after using PPL, which verifies the non-negative properties of the proposed method. Figures 9 and 10 show the time evolution of the solutions for simulating the tumor invasion haptotaxis system (1)-(3) Cases 1 and 2, where we use PPL in constructed scheme, respectively, where the snapshots of the numerical solutions for TC, ECM and MDE at $T=0,1,2,5,10,15,20$ are shown, respectively. It can be seen that the density of TC is attracted by MDE and gradually diffuses to its surrounding healthy tissues, while ECM is degraded by MDE and gradually destroyed by TC. This also fully reflects the rich dynamic behavior of the interaction between TC, ECM and MDE in the tumor microenvironment over time.

[See PDF for image]

Figure 5

Numerical solutions of TC, ECM, and MDE by the proposed method with/without PPL for Example 5.6 Cases 1

[See PDF for image]

Figure 6

Numerical solutions of TC, ECM, and MDE by the proposed method with/without PPL for Example 5.6 Cases 2

[See PDF for image]

Figure 7

1D profile of TC, and MDE along $x=0.5$ by using the proposed method with/without PPL for Example 5.6 Cases 1

[See PDF for image]

Figure 8

1D profile of TC, and MDE along $x=0.5$ by using the proposed method with/without PPL for Example 5.6 Cases 2

[See PDF for image]

Figure 9

Contour plots of the evolution for TC, ECM and MDE over time for Example 5.6 Cases 1, with $Nx=Ny=80$

[See PDF for image]

Figure 10

Contour plots of the evolution for TC, ECM and MDE over time for Example 5.6 Cases 2, with $Nx=Ny=50$

Conclusions

This paper focuses on the exploration of high-precision numerical methods for 1D and 2D haptotaxis models describing tumor invasion biotaxis phenomena in the biomedical field. The fourth-order compact finite difference method and the third-order SSP Runge-Kutta method are used to discretize the primitive differential system (1)-(3). To obtain their numerical solutions that maintain non-negative properties of the unknown function variables in such systems, a positivity-preserving algorithm with fourth-order accuracy is designed for the discrete function variables at each time step by using the idea of high-order numerical integration. Finally, numerical experiments are conducted to verify the accuracy and non-negativity of the proposed methods. And the dynamic behaviors of the solution to the haptotaxis models, which incorporates various forms of TC proliferation terms, ECM degradation and remodeling, MDE attenuation and proliferation are numerically simulated and analyzed. The methods proposed in this paper are efficient and easy to be extended to the higher dimensional problems, as well as the chemotaxis-haptotaxis models for tumor invasion and other complex mathematical models [37, 38].

Acknowledgements

Thanks for the valuable and constructive suggestions put forward by anonymous reviewers, which has greatly improved the quality of this paper. Thanks to the second and third authors of this paper for providing project funds support for the research on this topic.

Author contributions

LZ: Writing-original draft, investigation, formal analysis, programming. WG: Investigation, methodology, writing, review, revision and editing, visualization, validation, project administration, funding acquisition. YG: Conceptualization, methodology, investigation, data curation, writing draft, review and editing, visualization, project administration, funding acquisition, supervision. All authors read and approved the final manuscript.

Funding information

This work is supported by the Guangdong Basic and Applied Basic Research Foundation, China (2022A1515110436), the National Natural Science Foundation of China (12301622, 12161067), and National Science Foundation of Ningxia, China (2022AAC02023, 2023AAC03002).

Data availability

All the experimental data are true and reliable, which are obtained from the calculation of the method proposed in this paper.

Declarations

Competing interests

The authors declare that they have no competing of interests regarding the publication of this paper.

Publisher’s note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Qiu, H.; Cao, S.; Xu, R. Cancer incidence, mortality, and burden in China: a time-trend analysis and comparison with the United States and United Kingdom based on the global epidemiological data released in 2020. Cancer Commun.; 2021; 41, 10 pp. 1037-1048. [DOI: https://dx.doi.org/10.1002/cac2.12197]

2. Liotta, L.A.; Clair, T. Checkpoint for invasion. Nature; 2000; 405, 6784 pp. 287-288. [DOI: https://dx.doi.org/10.1038/35012728]

3. Gatenby, R.A.; Gawlinski, E.T. A reaction-diffusion model of cancer invasion. Cancer Res.; 1996; 56, 24 pp. 5745-5753.

4. Chaplain, M.A.J.; Lolas, G. Mathematical modelling of cancer cell invasion of tissue: the role of the urokinase plasminogen activation system. Math. Models Methods Appl. Sci.; 2005; 15, 11 pp. 1685-1734.2180714 [DOI: https://dx.doi.org/10.1142/S0218202505000947]

5. Giesselmann, J.; Kolbe, N.; Lukáčová-Medvid’ová, M.; Sfakianakis, N. Existence and uniqueness of global classical solutions to a two dimensional two species cancer invasion haptotaxis model. Discrete Contin. Dyn. Syst., Ser. B; 2018; 23, 10 pp. 4397-4431.3927826

6. Anderson, A.R.A.; Chaplain, M.A.J.; Newman, E.L.; Steele, R.J.C.; Thompson, A.M. Mathematical modelling of tumour invasion and metastasis. J. Theor. Med.; 2002; 2, 2 pp. 129-254.

7. Dzyubak, L.; Dzyubak, O.; Awrejcewicz, J. Nonlinear multiscale diffusion cancer invasion model with memory of states. Chaos Solitons Fractals; 2023; 168, 4534662 [DOI: https://dx.doi.org/10.1016/j.chaos.2022.113091] 113091.

8. Khaled-Abad, L.J.; Salehi, R. Application of weak Galerkin finite element method for nonlinear chemotaxis and haptotaxis models. Appl. Math. Comput.; 2021; 409, 4274891 126436.

9. Dehghan, M.; Narimani, N. An element-free Galerkin meshless method for simulating the behavior of cancer cell invasion of surrounding tissue. Appl. Math. Model.; 2018; 59, pp. 500-513.3788594 [DOI: https://dx.doi.org/10.1016/j.apm.2018.01.034]

10. Preziosi, L. Cancer Modelling and Simulation; 2003; 1 USA, Chapman Hall/CRC: [DOI: https://dx.doi.org/10.1201/9780203494899]

11. Kolev, M.K.; Koleva, M.N.; Vulkov, L.G. Two positivity preserving flux limited, second-order numerical methods for a haptotaxis model. Numer. Methods Partial Differ. Equ.; 2013; 29, 4 pp. 1121-1148.3053860 [DOI: https://dx.doi.org/10.1002/num.21748]

12. Kolev, M.K.; Koleva, M.N.; Vulkov, L.G. An unconditional positivity-preserving difference scheme for models of cancer migration and invasion. Mathematics; 2022; 10, 1 131. [DOI: https://dx.doi.org/10.3390/math10010131]

13. Niño-Celis, V.; Rueda-Gómez, D.A.; Villamizar-Roa, É.J. Convergence and positivity of finite element methods for a haptotaxis model of tumoral invasion. Comput. Math. Appl.; 2021; 89, pp. 20-33.4230499 [DOI: https://dx.doi.org/10.1016/j.camwa.2021.02.007]

14. Marciniak-Czochra, A.; Ptashnyk, M. Boundedness of solutions of a haptotaxis model. Math. Models Methods Appl. Sci.; 2010; 20, 3 pp. 449-476.2647028 [DOI: https://dx.doi.org/10.1142/S0218202510004301]

15. Enderling, H.; Anderson, A.R.A.; Chaplain, M.A.J.; Munro, A.J.; Vaidya, J.S. Mathematical modelling of radiotherapy strategies for early breast cancer. J. Theor. Biol.; 2006; 241, 1 pp. 158-171.2250157 [DOI: https://dx.doi.org/10.1016/j.jtbi.2005.11.015]

16. Gerisch, A.; Chaplain, M.A.J. Mathematical modelling of cancer cell invasion of tissue: local and non-local models and the effect of adhesion. J. Theor. Biol.; 2008; 250, 4 pp. 684-704.2945039 [DOI: https://dx.doi.org/10.1016/j.jtbi.2007.10.026]

17. Ganesan, S.; Lingeshwaran, S. Galerkin finite element method for cancer invasion mathematical model. Comput. Math. Appl.; 2017; 73, 12 pp. 2603-2617.3650000 [DOI: https://dx.doi.org/10.1016/j.camwa.2017.04.006]

18. Meral, G.; Surulesu, C. Mathematical modelling, analysis and numerical simulations for the influence of heat shock proteins on tumour invasion. J. Math. Anal. Appl.; 2013; 408, 2 pp. 597-614.3085056 [DOI: https://dx.doi.org/10.1016/j.jmaa.2013.06.017]

19. Aswin, V.S.; Manimaran, J.; Chamakuri, N. Space-time adaptivity for a multi-scale cancer invasion model. Comput. Math. Appl.; 2023; 146, pp. 309-322.4620077 [DOI: https://dx.doi.org/10.1016/j.camwa.2023.07.005]

20. Ganesan, S.; Lingeshwaran, S. A biophysical model of tumor invasion. Commun. Nonlinear Sci. Numer. Simul.; 2017; 46, pp. 135-152.3573862 [DOI: https://dx.doi.org/10.1016/j.cnsns.2016.10.013]

21. Kolbe, N.; Kat’uchová, J.; Sfakianakis, N. et al. A study on time discretization and adaptive mesh refinement methods for the simulation of cancer invasion: the urokinase model. Appl. Math. Comput.; 2016; 273, pp. 353-376.3427758

22. Irum, S.; Almakayeel, N.; Deebani, W. Computational insights into tumor invasion dynamics: a finite element approach. Math. Comput. Simul.; 2025; 229, pp. 814-829.4819021 [DOI: https://dx.doi.org/10.1016/j.matcom.2024.10.026]

23. Zhou, Y.; Chen, J.-W.; Dai, X.-N. et al. Numerical simulation of avascular tumor growth based on p27 gene regulation. Appl. Math. Mech.; 2013; 34, 3 pp. 327-338.3077556 [DOI: https://dx.doi.org/10.1007/s10483-013-1673-6]

24. Sulman, M.; Nguyen, T. A positivity preserving adaptive moving mesh method for cancer cell invasion models. J. Comput. Appl. Math.; 2019; 361, pp. 487-501.3955141 [DOI: https://dx.doi.org/10.1016/j.cam.2019.05.011]

25. Dai, F.; Liu, B. Optimal control and pattern formation for a haptotaxis model of solid tumor invasion. J. Franklin Inst.; 2019; 356, 16 pp. 9364-9406.4025628 [DOI: https://dx.doi.org/10.1016/j.jfranklin.2019.08.039]

26. Lele, S.K. Compact finite difference schemes with spectral-like resolution. J. Comput. Phys.; 1992; 103, 1 pp. 16-42.1188088 [DOI: https://dx.doi.org/10.1016/0021-9991(92)90324-R]

27. Orszag, S.A.; Israeli, M. Numerical simulation of viscous incompressible flows. Annu. Rev. Fluid Mech.; 1974; 6, 1 pp. 281-318. [DOI: https://dx.doi.org/10.1146/annurev.fl.06.010174.001433]

28. Peng, B.; Guo, X.; Zu, Y.; Tian, Z. A physics-preserving pure streamfunction formulation and high-order compact solver with high-resolution for three-dimensional steady incompressible flows. Phys. Fluids; 2023; 35, 4 [DOI: https://dx.doi.org/10.1063/5.0140054] 043104.

29. Kumar, N.; Dubey, R.K. A new development of sixth order accurate compact scheme for the Helmholtz equation. J. Appl. Math. Comput.; 2020; 62, 1 pp. 637-662.4056918 [DOI: https://dx.doi.org/10.1007/s12190-019-01301-x]

30. Wang, S.; Ge, Y.; Liu, S. Numerical solutions of the nonlinear wave equations with energy-preserving sixth-order finite difference schemes. Comput. Math. Appl.; 2024; 168, pp. 100-119.4756439 [DOI: https://dx.doi.org/10.1016/j.camwa.2024.05.028]

31. Gottlieb, S.; Shu, C.-W.; Tadmor, E. Strong stability-preserving high-order time discretization methods. SIAM Rev.; 2001; 43, 1 pp. 89-112.1854647 [DOI: https://dx.doi.org/10.1137/S003614450036757X]

32. Gottlieb, S.; Shu, C.-W. Total variation diminishing Runge-Kutta schemes. Math. Comput.; 1998; 67, 221 pp. 73-85.1443118 [DOI: https://dx.doi.org/10.1090/S0025-5718-98-00913-2]

33. Gao, G.-H.; Sun, Z.-Z. Compact difference schemes for heat equation with Neumann boundary conditions (II). Numer. Methods Partial Differ. Equ.; 2013; 29, 5 pp. 1459-1486.3092316 [DOI: https://dx.doi.org/10.1002/num.21760]

34. Zhang, L.; Peng, J.; Ge, Y. et al. High-accuracy positivity-preserving finite difference approximations of the chemotaxis model for tumor invasion. J. Comput. Biol.; 2024; 31, 12 pp. 1224-1258.4840153 [DOI: https://dx.doi.org/10.1089/cmb.2023.0316]

35. Liu, X.-D.; Osher, S. Nonoscillatory high order accurate self-similar maximum principle satisfying shock capturing schemes I. SIAM J. Numer. Anal.; 1996; 33, pp. 760-779.1388497 [DOI: https://dx.doi.org/10.1137/0733038]

36. Qiu, C.; Liu, Q.; Yan, J. Third order positivity-preserving direct discontinuous Galerkin method with interface correction for chemotaxis Keller-Segel equations. J. Comput. Phys.; 2021; 433, 4218537 [DOI: https://dx.doi.org/10.1016/j.jcp.2021.110191] 110191.

37. Yan, C.; Xu, J.; Wang, S. et al. Numerical study of convective heat transfer to supercritical CO2 in vertical heated tubes. Int. Commun. Heat Mass Transf.; 2022; 137, [DOI: https://dx.doi.org/10.1016/j.icheatmasstransfer.2022.106242] 106242.

38. Song, C.; Fang, R.-R.; Zhang, H.-L. et al. The exact solutions to a new type space reverse nonlocal Lakshmanan-Porserzian-Daniel equation. Nonlinear Dyn.; 2024; 112, 1 pp. 591-599. [DOI: https://dx.doi.org/10.1007/s11071-023-09057-7]