Numerical study of vanishing and spreading dynamics of chemotaxis systems with logistic source and a free boundary

TL;DR

This paper develops a numerical method to simulate the spreading and vanishing behavior of chemotaxis models with free boundaries, validating theoretical results and exploring factors influencing invasion dynamics.

Contribution

It introduces a front fixing finite difference scheme for chemotaxis models with free boundaries, addressing challenges in tracking moving fronts and nonlinear chemical interactions.

Findings

Numerical results confirm the vanishing-spreading dichotomy.

Simulations demonstrate the influence of initial conditions on spreading.

The method shows stability and positivity in solutions.

Abstract

The current paper is to investigate the numerical approximation of logistic type chemotaxis models in one space dimension with a free boundary. Such a model with a free boundary describes the spreading of a new or invasive species subject to the influence of some chemical substances in an environment with a free boundary representing the spreading front (see Bao and Shen [1], [2]). The main challenges in the numerical studies lie in tracking the moving free boundary and the nonlinear terms from chemical. To overcome them, a front fixing framework coupled with finite difference method is introduced. The accuracy of the proposed method, the positivity of the solution, and the stability of the scheme are discussed.The numerical simulations agree well with theoretical results such as the vanishing spreading dichotomy, local persistence, and stability. These simulations also validate some conjectures in our future theoretical studies such as the dependence of the vanishing-spreading dichotomy on the initial solution u0, initial habitat h0, the moving speed {\nu} and the chemotactic sensitivity coefficients \c{hi}1,\c{hi}2.

Equations

$$\{\begin{array}{l}{u}_t=u_{xx}-{\chi }_1(u{v}_{1,x}{)}_x+{\chi }_2(u{v}_{2,x}{)}_x+u(a(t,x)-b(t,x)u),0<x<h(t)\\ 0={\partial }_{xx}{v}_1-{\lambda }_1{v}_1+{\mu }_{1}u,0<x<h(t)\\ 0={\partial }_{xx}{v}_2-{\lambda }_2{v}_2+{\mu }_{2}u,0<x<h(t)\\ h^{\prime }(t)=-\nu u_x(t,h(t))\\ u_x(t,0)=v_{1,x}(t,0)=v_{2,x}(t,0)=0\\ u(t,h(t))=v_{1,x}(t,h(t))=v_{2,x}(t,h(t))=0\\ h(0)=h_0,u(x,0)=u_0(x),0\le x\le h_0,\end{array}$$

$$a_{\inf }:=\underset{t\in \mathbb{R},x\in [0,\infty )}{\inf }a(t,x)>0,b_{\inf }:=\underset{t\in \mathbb{R},x\in [0,\infty )}{\inf }b(t,x)>0.$$

$$C_{\mathrm{unif}}^b({\mathbb{R}}^{+})=\{u\in C({\mathbb{R}}^{+})|u(x)\text{is uniformly continuous and bounded on}{\mathbb{R}}^{+}\}$$

$$C_{\mathrm{unif}}^b(\mathbb{R})=\{u\in C(\mathbb{R})|u(x)\text{is uniformly continuous and bounded on}\mathbb{R}\}$$

$$\displaystyle M=\min\Big{\{}$$

$$\displaystyle\qquad\frac{1}{\lambda_{1}}\big{(}(\chi_{2}\mu_{2}\lambda_{2}-\chi_{1}\mu_{1}\lambda_{1})_{+}+\chi_{2}\mu_{2}(\lambda_{1}-\lambda_{2})_{+}\big{)}\Big{\}}$$

$$\displaystyle K=\min\Big{\{}$$

$$\displaystyle\quad\frac{1}{\lambda_{1}}\Big{(}|\chi_{1}\mu_{1}\lambda_{1}-\chi_{2}\mu_{2}\lambda_{2}|+\chi_{2}\mu_{2}|\lambda_{1}-\lambda_{2}|\Big{)}\Big{\}}.$$

$$M\le {\chi }_2{\mu }_2.$$

$$M_0=\frac{a_{\sup }}{b_{\inf }+{\chi }_2{\mu }_2-{\chi }_1{\mu }_1-M}$$

$$m_{0}=\frac{a_{\inf}\big{(}b_{\inf}-(1+\frac{a_{\sup}}{a_{\inf}})\chi_{1}\mu_{1}+\chi_{2}\mu_{2}-M\big{)}}{(b_{\inf}-\chi_{1}\mu_{1}+\chi_{2}\mu_{2}-M)(b_{\sup}-\chi_{1}\mu_{1}+\chi_{2}\mu_{2})}.$$

$$H(a,b)=\mathrm{cl}\{(a(t+\cdot ,\cdot ),b(t+\cdot ,\cdot ))|t\in \mathbb{R}\}$$

$$u_0\in C^2[0,h_0],u_0(x)\ge 0\mathrm{for}x\in [0,h_0],\mathrm{and}{u}_0^{\prime }(0)=0,u_0(h_0)=0,$$

$$0\le h^{\prime }(t)\le 2\nu M_1{C}_0,$$

$$0\le u(t,x;t_0,u_0,h_0)\le \max \{\Vert u_0{\Vert }_{\infty },M_0\}\forall t\in [t_0,\infty ),x\in [0,h(t;t_0,u_0,h_0))$$

$$\underset{t\to \infty }{\mathrm{limsup}}\underset{x\in [0,h(t;t_0,u_0,h_0))}{\sup }u(t_0+t,x;t_0,u_0,h_0)\le M_0,$$

$$h_{\infty }(t_0,u_0,h_0)=\underset{t\to \infty }{\lim }h(t;t_0,u_0,h_0).$$

$$\underset{t\to \infty }{\lim }\Vert u(t,\cdot ;t_0,u_0,h_0){\Vert }_{C([0,h(t;t_0,u_0,h_0)])}=0.$$

$$\underset{t\to \infty }{\mathrm{liminf}}\underset{0\le x\le L}{\inf }u(t,x;u_0,h_0)>0.$$

$$\{\begin{array}{l}{v}_t=v_{xx}+a(t,x)v,0<x<l\\ v_x(t,0)=v(t,l)=0.\end{array}$$

$$\underset{t\to \infty }{\mathrm{liminf}}\underset{0\le x\le L}{\inf }u(t,x;t_0,u_0,h_0)\ge {\tilde{m}}_0,$$

$$\underset{t\to \infty }{\mathrm{liminf}}\underset{0\le x\le L}{\inf }u(t,x;t_0,u_0,h_0)>m_0,$$

$$\underset{t\to \infty }{\lim }\underset{0\le x\le L}{\sup }|u(t,x;t_0,u_0,h_0)-u^{\ast }(t,x;a,b)|=0.$$

$$\underset{t\to \infty }{\lim }\underset{0\le x\le L}{\sup }|u(t,x;u_0,h_0)-u^{\ast }(t)|=0,$$

$$u^{{}^{\prime }}=u(a(t)-b(t)u)$$

$${\partial }_{xx}{v}_1$$

$${\partial }_{xx}{v}_2$$

$$u_t=u_{xx}+(-{\chi }_1{v}_{1x}+{\chi }_2{v}_{2x})u_x+(-{\chi }_1{\lambda }_1{v}_1+{\chi }_2{\lambda }_2{v}_2+a)u+({\chi }_1{\mu }_1-{\chi }_2{\mu }_2-b)u^2.$$

$$z(t,x)=\frac{x}{h(t)},w(t,z)=u(t,x),V_1(t,z)=v_1(t,x),V_2(t,z)=v_2(t,x).$$

$$\frac{{\partial }^2{V}_1}{\partial z^2}\cdot \frac{1}{h^2(t)}-{\lambda }_1{V}_1+{\mu }_{1}w=0,0<z<1,$$

$$\frac{{\partial }^2{V}_2}{\partial z^2}\cdot \frac{1}{h^2(t)}-{\lambda }_2{V}_2+{\mu }_{2}w=0,0<z<1.$$

$$V_{1,z}(t,0)=V_{2,z}(t,0)=0,$$

Full text

Key words. Chemoattraction-repulsion system, nonlinear parabolic equations, free boundary problem, spreading-vanishing dichotomy, front fixing, finite difference, invasive population.

AMS subject classifications. 35R35, 35J65, 35K20, 78M20, 92B05.

1 Introduction

The current paper is to study, in particular, numerically, the spreading and vanishing dynamiccs of the following attraction-repulsion chemotaxis system with a free boundary and logistic source,

[TABLE]

where $\nu>0$ in (1.1) is a positive constant, $a,b,\chi_{i}$, $\lambda_{i}$, and $\mu_{i}$ ($i=1,2$) are nonnegative constants, and $a(t,x)$ and $b(t,x)$ satisfy the following assumption,

(H0) $a(t,x)$* and $b(t,x)$ are bounded $C^{1}$ functions on $\mathbb{R}\times[0,\infty)$, and*

[TABLE]

Biological backgrounds of (1.1) are discussed in the paper ([1]). The free boundary condition in (1.1) is also derived in [1] based on the consideration of “population loss” at the front which assumes that the expansion of the spreading front is evolved in a way that the average population density loss near the front is kept at a certain preferred level of the species, and for each given species in a given homogeneous environment, this preferred density level is a positive constant determined by their specific social and biological needs, and the environment.

One of the first mathematical models of chemotaxis were introduced by Keller and Segel ([13], [14]) to describe the aggregation of certain type of bacteria in 1970. Since their publications, considerable progress has been made in the analysis of various particular case of chemotaxis (Keller-Segel) model on both bounded and unbounded fixed domain (see [3], [5], [6], [9], [10], [12], [25], [30], [31], [33], [34], [35], [36], [37], [38], [39], [40], and the references therein). Among the fundamental problems in studying chemotaxis model are the existence of nonnegative solutions which are globally defined in time or blow up at a finite time and the asymptotic behavior of time global solutions.

Du and Lin studied the population invasion represented by Fisher-KPP free boundary problem in 2010 [7]. The breaking difference between the asymptotic behaviours of Fisher-KPP with a free boundary and on the fixed or fixed unbounded domain is the vanishing-spreading dichotomy, which is well supported by some empirical evidences, for example, the introduction of several bird species from Europe to North America in the 1900s was successful only after many initial attempts (see [23],[29]).

Compared to the studying chemotaxis model on fixed bounded or fixed unbounded domain and the asymptotic behaviour of Fisher-KPP equation with a free boundary, the central problems in studying system (1.1) are the existence of nonnegative solutions which are globally defined in time, the vanishing-spreading dichotomy, local persistence, local stability, and so on.

To state the main results of the current paper, we first recall some theoretical results proved in [2] which will all be validated in our numerical simulations.

Let

[TABLE]

with norm $\|u\|_{\infty}=\sup_{x\in\mathbb{R}^{+}}|u(x)|$, and

[TABLE]

with norm $\|u\|_{\infty}=\sup_{x\in\mathbb{R}}|u(x)|$. Define

[TABLE]

and

[TABLE]

Let (H1)- (H3) be the following standing assumptions.

(H1) $b_{\inf}>\chi_{1}\mu_{1}-\chi_{2}\mu_{2}+M$.

(H2) b_{\inf}>\Big{(}1+\frac{a_{\sup}}{a_{\inf}}\Big{)}\chi_{1}\mu_{1}-\chi_{2}\mu_{2}+M.

(H3) $b_{\inf}>\chi_{1}\mu_{1}-\chi_{2}\mu_{2}+K$.

Note that

[TABLE]

Hence $b_{\inf}\geq\chi_{1}\mu_{1}$ implies (H1). In the case $\chi_{2}=0$, we can choose $\lambda_{2}=\lambda_{1}$, and then $M=0$ and $K=\chi_{1}\mu_{1}$. Hence (H1) becomes $b_{\inf}>\chi_{1}\mu_{1}$, (H2) becomes $b_{\inf}>(1+\frac{a_{\sup}}{a_{\inf}})\chi_{1}\mu_{1}$, and (H3) becomes $b_{\inf}>2\chi_{1}\mu_{1}$. In the case $\chi_{1}=0$, we can also choose $\lambda_{1}=\lambda_{2}$, and then $M=\chi_{2}\mu_{2}$ and $K=\chi_{2}\mu_{2}$. Hence (H1) (resp.(H2), (H3)) becomes $b_{\inf}>0$. Biologically, (H1), (H2), and (H3) indicate that the chemo-attraction sensitivity is relatively small with respect to logistic damping.

When (H1) holds, we put

[TABLE]

and

[TABLE]

Note that if (H2) holds, then $m_{0}>0$.

Let

[TABLE]

with open compact topology, where the closure is taken under the open compact topology.

The main results of the paper [2] are stated in the following. The first result is on the global existence of nonnegative solutions of (1.1).

*Global existence [2, Theorem 1.2]:

If (H1) holds, then for any $t_{0}\in\mathbb{R}$, and any $h_{0}>0$ and any function $u_{0}(x)$ on $[0,h_{0}]$ satisfying*

[TABLE]

(1.1) has a unique globally defined solution $(u(t,x;t_{0},u_{0},h_{0})$, $v_{1}(t,x;t_{0},u_{0},h_{0})$, $v_{2}(t,x;t_{0},u_{0},h_{0})$, $h(t;t_{0},u_{0},h_{0}))$ with $u(t_{0},x;t_{0},u_{0},h_{0})=u_{0}(x)$ and $h(t_{0};t_{0},u_{0},h_{0})=h_{0}$. Moreover,

[TABLE]

[TABLE]

and

[TABLE]

*where $M_{1}$ is a big enough constant and $C_{0}=\max\{\|u_{0}\|_{\infty},M_{0}\}$. *

Assume (H1). For any given $t_{0}\in\mathbb{R}$, and any given $h_{0}>0$ and $u_{0}(\cdot)$ satisfying (1.6), by the nonnegativity of $u(t,x;t_{0},u_{0},h_{0})$, $h^{{}^{\prime}}(t;t_{0},u_{0},h_{0})\geq 0$ for all $t>t_{0}$. Hence $\lim_{t\to\infty}h(t;t_{0},u_{0},h_{0})$ exists. Put

[TABLE]

We say vanishing occurs if $h_{\infty}(t_{0},u_{0},h_{0})<\infty$ and

[TABLE]

We say spreading occurs if $h_{\infty}(t_{0},u_{0},h_{0})=\infty$ and for any $L>0$,

[TABLE]

For given $l>0$, consider the following linear equation,

[TABLE]

Let $[\lambda_{\min}(a,l),\lambda_{\max}(a,l)]$ be the principal spectrum interval of (1.10) (see Definition 2.1 [2]). Let $l^{*}>0$ be such that $\lambda_{\min}(a,l)>0$ for $l>l^{*}$ and $\lambda_{\min}(a,l^{*})=0$ (see Lemma 2.2, 2.3 [2] for the existence and uniqueness of $l^{*}$).

The second result is about the spreading and vanishing dichotomy scenario in (1.1).

Spreading-vanishing dichotomy [2, Theorem 1.3]: Assume that (H1) holds. For any given $t_{0}\in\mathbb{R}$, and $h_{0}>0$ and $u_{0}(\cdot)$ satisfying (1.6), we have that either (i) vanishing occurs and $h_{\infty}(t_{0},u_{0},h_{0})\leq l^{*}$; or (ii) spreading occurs.

For given $t_{0}\in\mathbb{R}$, and $h_{0}>0$ and $u_{0}(\cdot)$ satisfying (1.6), if spreading occurs, it is interesting to know whether local uniform persistence occurs in the sense that there is a positive constant $\tilde{m}_{0}$ independent of the initial data such that for any $L>0$,

[TABLE]

and whether local uniform convergence occurs in the sense that $\lim_{t\to\infty}u(t,x;t_{0},u_{0},h_{0})$ exists locally uniformly. We have the following result along this direction.

Persistence and convergence [2, Theorem 1.4]:

Assume that (H1) holds and that $h_{0}>0$ and $u_{0}(\cdot)$ satisfy (1.6).

(i) (Local uniform persistence) For any given $t_{0}\in\mathbb{R}$, if $h_{\infty}(t_{0},u_{0},h_{0})=\infty$ and (H2) holds, then for any $L>0$,

[TABLE]

where $m_{0}$ is as in (1.5).

(ii) (Local uniform convergence) Assume that (H3) holds, and that for any $(\tilde{a},\tilde{b})\in H(a,b)$, there has a unique strictly positive entire solution $(u^{*}(t,x;\tilde{a},\tilde{b}),v_{1}^{*}(t,x;\tilde{a},\tilde{b})$, $v_{2}^{*}(t,x;\tilde{a},\tilde{b}))$. Then for any given $t_{0}\in\mathbb{R}$, if $h_{\infty}(t_{0},u_{0},h_{0})=\infty$, there are $\chi_{1}^{*}>0,\chi_{2}^{*}>0$ such that to any $0\leq\chi_{1}\leq\chi_{1}^{*}$, $0\leq\chi_{2}\leq\chi_{2}^{*}$, for any $L>0$,

[TABLE]

(iii) (Local uniform convergence) Assume that (H3) holds, and that $a(t,x)\equiv a(t)$ and $b(t,x)\equiv b(t)$. Then for any given $t_{0}\in\mathbb{R}$, if $h_{\infty}(t_{0},u_{0},h_{0})=\infty$, then for any $L>0$,

[TABLE]

where $u^{*}(t)$ is the unique strictly positive entire solution of the ODE

[TABLE]

(see [11, Lemma 2.5] for the existence and uniqueness of strictly positive entire solutions of (1.12)).

Remark 1.1**.**

Biologically, the invasion or spreading of the population is depending on the initial solution, initial habitat, the moving speed $\nu$ ([7]). When the spreading happens, the local persistence and convergence can be guaranteed in Fisher-KPP equation with a free boundary, furthermore, there is a asymptotic spreading speed such that $\lim_{t\to\infty}\frac{h(t)}{t}=c^{*}>0$ (see [7], [16]).

Compared to the vanishing-spreading dichotomy in Fisher-KPP equation with a free boundary, the chemotaxis system (1.1) do not have comparison principle which leads the following interesting open problems but has positive answers in Fisher-KPP case.

For given $u_{0}(\cdot)$ and $0<h_{0}<l^{*}$, whether there is $\nu^{*}>0$ such that for $0<\nu\leq\nu^{*}$, vanishing occurs, and for $\nu>\nu^{*}$, spreading occurs. 2. 2.

For given $\nu>0$, $\phi(\cdot)$, and $0<h_{0}<l^{*}$, whether there is $\sigma^{*}>0$ such that for $u_{0}=\sigma\phi$ with $\sigma\leq\sigma^{*}$, vanishing occurs, and for $u_{0}=\sigma\phi$ with $\sigma>\sigma^{*}$, spreading occurs. 3. 3.

Whether there is a spreading speed $c^{*}>0$ such that $\lim_{t\to\infty}\frac{h(t)}{t}=c^{*}$ as long as the spreading occurs.

Remark 1.2**.**

Vanishing-spreading result [2, Theorem 1.2] indicate there is a separating value $l^{*}$, which is independent of the chemotactic sensitivity coefficients $\chi_{1},\chi_{2}$, such that in the vanishing scenario the limiting moving boundary $h_{\infty}<l^{*}$ and when the initial habitat $h_{0}>l^{*}$ the spreading guaranteed. The dependence of the dynamics of the system on the chemotactic sensitivity coefficients is another important and interesting questions [28], [38]. We also have the following question in this direction.

4.

If the asymptotic spreading speed exist, whether the limit $\lim_{t\to\infty}\frac{h(t)}{t}=c^{*}$ depends on the chemotactic sensitivity coefficient $\chi_{1}$ and $\chi_{2}$.

The objective of the current paper is to study the numerical effect of the parameters $\sigma,\nu,\chi_{1},\chi_{2}$ on the vanishing and spreading dynamics in the system (1.1) which will give us directions in the theoretical studies. For simplicity, we only consider the constant logistic coefficients in the system, where $a(t,x)=a,b(t,x)=b$. In general, it is always difficult to handle the attraction term in chemotaxis system which may lead to convection dominant in the system [4], [19], [20], [26]. However, in the system (1.1), we have extra numerical challenges in efficiently and accurately handling the moving boundaries [24]. These two challenges require us to construct a new numerical algorithm in the numerical study.

Thanks to the maximum principle in the elliptic equation and death damping coefficient $b$ in the parabolic equation, the chemoattraction term can be controlled by the magnitude of the population density which has a global bound by the suppressing of the death rate. The front-fixing method has been successfully applied to solve one dimensional free boundary problem [15], [21], [22], [24] which changes the moving boundary to fixed domain, and is the main concern in our numerical studies of (1.1). Combined with the front-fixing method, finite difference in parabolic and finite volume method in the elliptic equations in the system (1.1), we construct a new algorithm in the numerical study, as a by-product we also obtain the consistency, monotonicity of the moving boundary, positivity of the solution and stability results.

Our numerical experiments validate the vanishing and spreading dichotomy in the numerical scheme of system (1.1) which is similar to Fisher-KPP equation with a free boundary and give evidences to our conjectures that:

For given $u_{0}(\cdot)$ and $0<h_{0}<l^{*}$, there is $\nu^{*}>0$ such that for $0<\nu\leq\nu^{*}$, vanishing occurs, and for $\nu>\nu^{*}$, spreading occurs. Which means in order to spread to the half space $\mathbb{R}^{+}$, the moving speed $\nu$ should be large enough and otherwise the population will be extinct. 2. 2.

For given $\nu>0$, $\phi(\cdot)$, and $0<h_{0}<l^{*}$, there is $\sigma^{*}>0$ such that for $u_{0}=\sigma\phi$ with $\sigma\leq\sigma^{*}$, vanishing occurs, and for $u_{0}=\sigma\phi$ with $\sigma>\sigma^{*}$, spreading occurs. Biologically, large initial population density helps the establishment and spreading of invasion which is an indirect evidence to the early birds introduction problem in 1900s [23], [29]. 3. 3.

There is a spreading speed $c^{*}>0$ such that $\lim_{t\to\infty}\frac{h(t)}{t}=c^{*}$ as long as the spreading occurs, which is independent of chemotactic sensitivity coefficient $\chi_{1}$ and $\chi_{2}$. Chemical $v_{1}$ and $v_{2}$ are produced by the species and the density is close to zero near the spreading front. In such case the decisive effect of the spreading speed should not depend on the chemotactic sensitivity coefficients $\chi_{1}$ and $\chi_{2}$.

The rest of this paper is organized in the following way. In section 2, we first use Landau transformation to transfer the moving boundary to a fixed domain, then we use the finite difference, finite volume, and iteration method to approximate the continuous chemotaxis system. We also prove the monotonicity of the moving boundary, the positivity and stability of the discrete solutions. In section 3, we study the numerical spreading-vanishing dichotomy in (1.1) which validates our theoretical results (Vanishing-spreading dichotomy, local persistence and convergence). Our simulations also indicate the dependence or independence of the vanishing-spreading dichotomy on parameters $\nu,u_{0},h_{0},\chi_{1},\chi_{2}$ and so on. In section 4, some future works are briefly discussed.

2 Numerical approximation of the free boundary problem

In this section, we study the numerical approximation of system (1.1) with constant logistic coefficients $a(t,x)=a,b(t,x)=b$. First through the well-known Landau transformation (see [15]), we convert (1.1) into a fixed spatial domain problem. In such a way the length of the moving boundary is included as another variable to be solved apart from the population density. Then we solve the converted new problem on the basis of finite difference and finite volume method. There is a circulation that each time when time variable increase we solve the elliptic equations first and using the forward differential method to find the solution of the parabolic equation for the new time.

From the elliptic equations in (1.1), we know that

[TABLE]

Combining (2.1), (2.2) and the first equation in (1.1), we have

[TABLE]

Then we introduce the Landau transformation,

[TABLE]

Under this substitution the elliptic equations (2.1), (2.2) take the form:

[TABLE]

The elliptic boundary conditions are

[TABLE]

Equation (2.3) takes the form:

[TABLE]

Let $G(t)$ denote $h^{2}(t)$ and multiply it on both sides of the above equation.

[TABLE]

Boundary conditions and Stefan condition take the form

[TABLE]

and

[TABLE]

The initial conditions in (1.1) become:

[TABLE]

and the initial function $u_{0}(x)$ is changed into $w_{0}(z)$ which maintains:

[TABLE]

Under the transformation, our aim is to solve the nonlinear parabolic partial differential system (2.5), (2.6), (2.10) in the fixed domain $(0,\infty)\times(0,1)$ for the variables $(t,z)$.

The following is the process according to the theory of the finite difference method. First we consider the time and space discretization $\tau=\triangle t$, $h=\triangle z=1/M$, which means the interval (0,1) is divided into $M$ equal cells

[TABLE]

and the mesh points $(t^{n},z_{j})$, with $t^{n}=n\tau,n\geq 0,z_{j}=jh,0\leq j\leq M$. For abbreviation, the approximate value of $w(t^{n},z_{j})$ can be denoted by $w^{n}_{j}$, the approximate values of $V_{1}(t^{n},z_{j})$ and $V_{2}(t^{n},z_{j})$ can be denoted by $V_{1,j}^{n}$ and $V_{2,j}^{n}$. Besides we write $g^{n}$ for the value of $G(t^{n})$. Let us consider the central approximation of the spatial derivatives,

[TABLE]

[TABLE]

By using the forward approximation of the time derivative, we get

[TABLE]

Let us apply the approximation (2.15) on the elliptic equations (2.5), (2.6), then we get

[TABLE]

We mainly focus on solving the values of $V_{1}$, the relevant results about $V_{2}$ can be obtained in a similar way. For (2.15), $1\leq j\leq M-1$, we can get $M-1$ equations for $V_{1}$ and there are two other equations we need for the boundary. In order to achieve a higher $O(h^{2})$ accuracy (see [18]), we use the idea of finite volume method to handle the boundary conditions.

We get the equation for the left boundary:

[TABLE]

where

[TABLE]

and $z_{1/2}=\frac{z_{0}+z_{1}}{2}$.

The equation for the right boundary is similar:

[TABLE]

where

[TABLE]

and $z_{M-1/2}=\frac{z_{M-1}+z_{M}}{2}$.

By now we have $M+1$ equations which is enough to form a system of linear algebraic equations for $V_{1}$. It is tridiagonal and there exists the unique solutions $V_{1,0}^{n},V_{1,1}^{n},\cdots,V_{1,M}^{n}$. Similarly, we can acquire the system of linear algebraic equations about $V_{2}$ and its corresponding solutions $V_{2,0}^{n},V_{2,1}^{n},\cdots,V_{2,M}^{n}$.

With the information of $V_{1}$ and $V_{2}$, we can concentrate on solving the parabolic equation (2.10). From (2.15), (2.16), and (2.17), (2.10) is approximated by

[TABLE]

for $n\geq 0,\;0\leq j\leq M-1$. Because of the initial conditions (2.14), we can assume a fictitious value $w_{-1}^{n}$ at the point $(t^{n},-h)$, and then

[TABLE]

Considering the Stefan condition (2.12), according to (2.17) and three points backward spatial approximation of $\frac{\partial w}{\partial z}(t,1)$, we obtain

[TABLE]

Because of (2.13), it can also be written as:

[TABLE]

Let us replace $g^{n+1}$ with (2.25) in (2.22), we get the explicit scheme:

[TABLE]

for $n\geq 0,\;0\leq j\leq M-1,$ where

[TABLE]

The solution of (2.10) is classical (see [2, Lemma 3.1]), we can obtain a optimal error estimates for the numerical approximation. Consider (2.10), (2.11) and (2.12), we denote that

[TABLE]

From (2.22), we let

[TABLE]

for $n\geq 0,\;0\leq j\leq M-1$. Also from (2.23) and (2.24), let

[TABLE]

Then we have the following error estimates for the parabolic equation.

Proposition 2.1**.**

(Error estimates) Under conditions of global existence ([2, Theorem 1.2]), let $w$ be the solution of Equation (2.30), (2.31), (2.32) and $w_{h}$ be the numerical solution of (2.33), (2.34), (2.35). Then we have the following error estimates:

[TABLE]

The proof of Proposition 2.1.

At the point $(t^{n},z_{j})$, we have the error estimate

[TABLE]

By using Taylor’s expansion at the point $(t^{n},z_{j})$, we have the following result

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Similarly, we have

[TABLE]

Using estimates (2.43) and (2.44), and the linear combinations we obtain

[TABLE]

Then from (2.40),(2.41),(2.42), and (2.45), we have

[TABLE]

By the central difference approximation, we have the error estimate:

[TABLE]

Using the three points derivative formula (2.24), we obtain

[TABLE]

∎

From the above results, we have the following error estimate for the transferred system (1.1).

Theorem 2.1** (Error estimate).**

Let $L_{h}(w,V_{1},V_{2},g)$ denote all the finite difference scheme of the transferred system (1.1) including the boundary approximations, and the corresponding continuous scheme as $L(w,V_{1},V_{2},g)$, then $L_{h}(w,V_{1},V_{2},g)$ is consistent with $L(w,V_{1},V_{2},g)$ and the local truncation error is

[TABLE]

Proof.

Combining Proposition 2.1 for the parabolic equation estimate and the elliptic equation estimate above (central finite difference for the elliptic equation and finite volume scheme for the boundary), we have the local truncation error estimate (2.49). ∎

Before proving the positivity of the density function $w$ and the monotonicity property of the spreading front $g^{n}$, we need the following maximum principle of the discrete elliptic equations.

Proposition 2.2**.**

Consider equations (2.18) and (2.19), we have the following discrete maximum principle

[TABLE]

where $w^{n}_{\min},w^{n}_{\max}$ are the minimum and maximum solution of Equation (2.22).

Proof.

We only study the discrete maximum principle of Equation (2.18), the proof for Equation (2.19) is the same. Rewrite Equation (2.18) as the following

[TABLE]

Let $j=j^{*}(j^{*}\neq 0,M)$ be the point where $V_{1}$ takes the maximum at $t^{n}$ and denote it as $V_{1,max}^{n}$, then

[TABLE]

Because $V_{1,max}^{n}-V_{1,j^{*}-1}^{n}\geq 0$, $V_{1,max}^{n}-V_{1,j^{*}+1}\geq 0$, and $2h\mu_{1}G(t)>0$, we have

[TABLE]

where $w_{max}^{n}$ is the maximum of $w$ at $t^{n}$. We then conclude that

[TABLE]

Similarly we can obtain

[TABLE]

∎

The moving boundary $h(t)$ is monotonicity in the system (1.1) which is also preserved in the discretized one.

Theorem 2.2** (Monotonicity of spreading front).**

Let $\chi_{1},\chi_{2}$ be small enough and $\tau$ satisfy

[TABLE]

then $w^{n}_{M-1}>0,$ and $g^{n}$ is monotonicity in $n$ with $n\geq 0,C=e^{aT}w_{M-1}^{0}$.

Proof.

We prove the positivity and monotonicity of the free boundary $g^{n}$ by using the induction principle on the index $n$. For $n=0$, from the initial condition of $w_{0}(z)$, we have $w_{j}^{0}>0,\;0\leq j\leq M-1,w_{M}^{0}=0$. Additionally, by Hopf lemma the left derivative of $w^{0}$ at $z_{M}$ is negative and hence the corresponding difference approximation of Equation (2.24) with small enough $h$ has

[TABLE]

which is equivalent to $(w_{M-2}^{0}-4w_{M-1}^{0})<0$. So combined with $g^{0}>0$ and (2.24) with small enough $h$, we have

[TABLE]

By using Taylor’s expansion on the left of $z_{M}=1$ at $t^{n}$, we have

[TABLE]

Let $w=w_{M-2}^{n}$, then we have

[TABLE]

combined with $w_{M}^{n}=0$, we obtain

[TABLE]

Plugging (2.57) into (2.24), we have

[TABLE]

If $w_{M-1}^{n}$ is positive, then from (2.58) $g^{n}$ is positive and increasing with $n$. By using (2.57) and let $j=M-1$, from (2.26) we have

[TABLE]

In order to preserve the stability in the forward approximation of parabolic equation, one needs the requirement such that $\frac{\tau}{h^{2}}<\frac{1}{2}$, and from the discrete maximum principle the bound of $V_{1}^{n}$ and $V_{2}^{n}$ are controlled by $w^{n}$ which is bounded by iteration. Assume $\chi_{1}$ and $\chi_{2}$ are small enough, we have the approximation

[TABLE]

In the following, we prove that $w_{M-1}^{n}$ has a up bound independent of $n$ (n depends on T). By induction, we assume that $g^{n}>g^{n-1}>\cdots>g^{1}>g^{0}$. With small enough parameters $\tau,h,\chi_{1},\chi_{2}$, and $g^{n}>g^{0},z_{M-1}<1$ we have

[TABLE]

where $\phi_{M-1}^{n}=1+\tau(a-bw_{M-1}^{n})-\frac{\tau}{h^{2}}\frac{1}{g^{0}}\nu w_{M-1}^{n}$. From (2.60) we obtain $w_{M-1}^{n+1}<(1+\tau a)w_{M-1}^{n}$ and

[TABLE]

Denote the total time from $t^{0}$ to $t^{n+1}$ as $T$ and we have $T=(n+1)\tau$, then

[TABLE]

which leads to $w_{M-1}^{n}<e^{aT}w_{M-1}^{0}$. Furthermore, if the time step $\tau$ satisfies

[TABLE]

we obtain $\phi^{n}_{M-1}>0$ and $w^{n+1}_{M-1}>0$ where $C=e^{aT}w_{M-1}^{0}$. Then combined Equation (2.58) with the positivity of $w^{n}_{M-1}$, we conclude the monotonicity of the spreading front $g^{n}$. ∎

Theorem 2.3** (Positivity and boundedness of the discrete solution).**

In Equation (2.33), let $\tau$ satisfy

[TABLE]

with small $\chi_{1},\chi_{2}$ and let $M$ large enough which is equivalent to small $h$, then the solution of (2.33) satisfies $w_{j}^{n}\geq 0$ with uniform up-bound, for $0\leq j\leq M,n\geq 0$.

Proof.

From Equations (2.26), (2.27), (2.28), (2.29), we can rewrite $w_{j}^{n+1}$ as

[TABLE]

where

[TABLE]

If $a_{j}^{n},b_{j}^{n},c_{j}^{n}$ are positive, the positivity of $w^{n}_{j},n\geq 0,0\leq j\leq M$ can be proved by induction. First, we consider

[TABLE]

From the discrete maximum principle, regularity of $w$ at $z_{M}=1(w^{n}_{M}=0),$ by taking large enough $M$ or small enough $h$ and small enough chemotactic sensitivity $\chi_{1},\chi_{2}$, then plugging (2.57) into (2.67) and combined with the positivity of $w^{n}_{M-1}$ and $0\leq z_{j}<1$, we have

[TABLE]

with negative $S_{1}$, or stronger condition

[TABLE]

with positive $S_{1}$. Both conditions lead to $a_{j}^{n}>0$.

Next, we study the positivity of $c_{j}^{n}$.

[TABLE]

Plugging (2.57) into (2.72) and combine with the positivity of $w^{n}_{M-1}$ and small enough $\chi_{1},\chi_{2}$, and $h$, the positivity of $c_{j}^{n}$ can be guaranteed.

Thirdly, we investigate the positivity of

[TABLE]

In order to obtain $b_{j}^{n}>0$, we need

[TABLE]

(While $\chi_{1},\chi_{2}$ are small enough, $\tau(\frac{2}{g^{n}h^{2}}-S_{2}-S_{3}w_{j}^{n})>0$.) By the discrete maximum principle of the elliptic equation $V_{1}$ and $V_{2}$, the monotonicity of $g^{n}$ $(g^{0}<g^{n}),$ we can improve the requirement (2.74) to

[TABLE]

where $|w^{n}_{max}|=\max_{0\leq j\leq M}|w^{n}_{j}|$, which leads to the positivity of $b_{j}^{n}$.

The positivity of $w_{j}^{n},0\leq j\leq M,n\geq 0$ is followed from the positivity of $a_{j}^{n},b_{j}^{n},c_{j}^{n}$.

For the uniform up-bound, we use the induction method again. By using discrete maximum principle of $V_{1}^{n}$ and $V_{2}^{n}$, whose values are controlled by maximum of $w^{n}$, $S_{3}<0$ and positivity of $a_{j}^{n},b_{j}^{n},c_{j}^{n}$, we have estimate in (2.26) such that

[TABLE]

Because of the existence of up bound for $w^{n}$, there exist a constant $U$ such that $|w^{n}_{max}|\leq U$. Furthermore, we have the following uniform up bound estimate

[TABLE]

Denote the total time $T=n\tau$, we have

[TABLE]

and the requirement of (2.75) can be sharped to

[TABLE]

which guarantee the positivity of the discrete solution $w^{n}$ for $0\leq j\leq M,n\geq 0$. ∎

Before discussing about the stability of our algorithm and for the sake of clarity in the presentation we specify the concept of stability we use below. We recall the definition of the supremum norm of a vector $x=(x_{1},x_{2},\cdots,x_{n})^{T}$ in $\mathbb{R}^{n}$ as $\|x\|_{\infty}=\max(|x_{1}|,|x_{2}|,\cdots,|x_{n}|).$

Definition 2.1**.**

In Equation (2.33), the numerical scheme is said to be $\|\cdot\|_{\infty}$ stable in the domain $[0,T]\times[0,1]$, if for every partition with $T=N\tau,Mh=1$ it hold true that:

[TABLE]

where $w^{n}=[w_{0}^{n},w_{1}^{n},\cdots,w_{M}^{n}]^{T}$is the vector solution of the scheme at $t^{n}$, $K$ is a constant independent of $h,\tau,n$.

From the results in Theorem 2.3 and Estimate (2.78), we have the uniform stability result.

Theorem 2.4** (Stability of discrete solution).**

In Equations (2.18), (2.19) and (2.33), let $\tau$ satisfy

[TABLE]

with small $\chi_{1},\chi_{2}$ and let $M$ large enough which is equivalent to small $h$, then the discrete solution of (1.1) is $\|\cdot\|_{\infty}$ stable.

Remark 2.1**.**

In the investigation, (2.80) is a strong requirement and it is a sufficient condition for the positivity and stability. However, our simulations indicate larger $h,\tau$ can also guarantee the positivity and stability of the scheme.

3 Numerical experiments

In this section, we study the numerical simulations of the free boundary problem of (1.1) in the case that $a(t,x)=2,b(t,x)=1$. Our theoretical results indicate there exists a critical value $l^{*}=\frac{\pi}{2}\sqrt{\frac{1}{a}}$, which is independent of the chemotactic sensibility coefficients $\chi_{1},\chi_{2}$, such that spreading of the species is guaranteed for $h_{0}\geq l^{*}$ and vanishing happens for $h_{0}<l^{*}$. Furthermore, in order to obtain local convergence and persistency, chemotactic sensitivity must be small enough.

Compared to Fisher-KPP free boundary problems, $\chi_{1}=\chi_{2}=0$ in the system (1.1), even if $h_{0}<l^{*}$, the spreading is guaranteed under condition $\nu>\nu^{*}>0$ or $\sigma>\sigma^{*}>0$, where $u_{0}(x)=\sigma\phi(x)$ and $\nu^{*}$ is an unknown threshold depending on $u_{0}$ (see [7, Theorem 3.9]). Because of the lack of comparison principle in Chemotaxis system (1.1), the existence of $\nu^{*}$ and $\sigma^{*}$ is still an open problem. However, our numerical simulation do validate such existence. Furthermore, when the spreading happens in Fisher-KPP free boundary problem, we have the following asymptotic spreading speed result:

[TABLE]

where $c^{*}$ is depending on the logistic coefficient $a(t,x)$ (see [7], [8], and [17]). Such result is also confirmed in our numerical simulations and the theoretical analysis will be under our investigation in the future.

In the following, we show different simulation results depending on different parameter selections. All parameters satisfy conditions (H1)-(H3).

3.1 Numerical vanishing-spreading dichotomy

Compared to the vanishing-spreading dichotomy in Fisher-KPP free boundary problems, vanishing happens when the initial habitat $h_{0}$, initial solution $u_{0}(x)=\sigma\phi(x)$, and moving speed $\nu$ are small enough; spreading happens when either $h_{0},u_{0}(x),$ or $\mu$ is big enough. However, because of the lack of comparison principle, similar results are still open in the chemotaxis system (1.1). The following numerical simulations validate these similar results in chemotaxis system.

Exmple 3.1**.**

In the logistic chemotaxis model (1.1), let $h_{0}=2.5>l^{*}=1.11,u_{0}=\cos(\pi x/2h_{0})$ and $(\chi_{1},\chi_{2},\nu,\lambda_{1},\lambda_{2},\mu_{1},\mu_{2})=(0.02,0.01,0.01,2,1,2,1)$. Figure 2 and 2 show the system has spreading tendency and the asymptotic speed $\frac{h(t)}{t}$ converges to a constant.

Exmple 3.2**.**

In the logistic chemotaxis model (1.1), let $h_{0}=0.5<l^{*}=1.11,u_{0}=\cos(\pi x/2h_{0})$, and $(\chi_{1},\chi_{2},\nu,\lambda_{1},\lambda_{2},\mu_{1},\mu_{2})=(2,1,0.8,1,2,1,2)$. Figure 4 and 4 show the system has vanishing tendency and the asymptotic speed $\frac{h(t)}{t}$ converges to zero.

Exmple 3.3**.**

Increase the initial habitat $h_{0}$ in Example 3.2 such that $h_{0}=2.5>l^{*}=1.11$, and let $(\chi_{1},\chi_{2},\nu,\lambda_{1},\lambda_{2},\mu_{1},\mu_{2})=(2,1,0.8,1,2,1,2)$. Figure 6 and 6 show the system has spreading tendency and the asymptotic speed $\frac{h(t)}{t}$ converges to a positive constant.

3.2 Critical value of the moving speed $\nu^{*}$

In the following, we study the influence of the moving speed $\nu$ on the vanishing-spreading dichotomy. Compared to Fisher-KPP free boundary problem, there exists a critical value $\nu^{*}$ such that vanishing happens when $\nu<\nu^{*}$ and spreading happens in other direction. However, because of the lack of comparison principle, whether there exists a critical value $\nu^{*}$ in chemotaxis free boundary problem is still an open question, but numerically we find the existence of such $\nu^{*}$ by dichotomy method.

Exmple 3.4**.**

In the logistic chemotaxis model (1.1), let $h_{0}=1.0<l^{*}=1.11,u_{0}=\cos(\pi x/2h_{0})$ and $(\chi_{1},\chi_{2},\nu,\lambda_{1},\lambda_{2},\mu_{1},\mu_{2})=(0.2,0.1,2,1,2,1,2)$. Figure 8 and 8 show the system has a spreading tendency and the asymptotic speed $\frac{h(t)}{t}$ converges to a positive constant.

Exmple 3.5**.**

We only change the moving speed to a smaller positive number $\nu=0.01$ in Exmple 3.4. Figure 10 and 10 show the system (1.1) has a tendency of vanishing.

Exmple 3.6**.**

In the logistic chemotaxis model (1.1), let $h_{0}=1.0<l^{*}=1.11,u_{0}=\cos(\pi x/2h_{0})$ and $(\chi_{1},\chi_{2},\lambda_{1},\lambda_{2},\mu_{1},\mu_{2})=(0.2,0.1,1,2,1,2)$. By the dichotomy method, the simulations indicate the critical $\nu^{*}$ is between $0.05$ and $0.025$ (see Figure 12 and 12).

3.3 Vanishing-spreading dichotomy dependence on $u_{0}$

In the following, we study the vanishing-spreading dependence on initial solution $u_{0}(x)$.

Exmple 3.7**.**

We first investigate the system (1.1) with large initial solution $u_{0}(x)=4\cos(\pi x/2h_{0})$ with the following parameters $(h_{0},\chi_{1},\chi_{2},\nu,\lambda_{1},\lambda_{2},\mu_{1},\mu_{2})=(1,0.2,0.1,0.8,1,2,1,2)$. The system has a tendency of spreading and converges to the constant $\frac{a}{b}=2$ (see Figure 14,14).

Exmple 3.8**.**

In the case of small initial solution $u_{0}=0.01\cos(\pi x/2h_{0})$ and with fixed other parameters as in Exmple 3.7, the system has a tendency of vanishing (see Figure 16, 16).

In summary, the system (1.1) has tendency of vanishing when the initial solution $u_{0}(x)$ is small and the moving boundary converges to a constant less than $l^{*}$. In the case of large initial solution, the system has tendency of spreading whose spreading speed also converges to a constant which is similar to the Fisher-KPP free boundary problems.

3.4 The dependence of asymptotic spreading speed $\frac{h(t)}{t}$ on parameters

In the following, we mainly focus on the spreading speed’s dependence on parameters $u_{0},\nu$, and the chemotactic sensitivity coefficients $\chi_{1},\chi_{2}$ when spreading happens.

We first consider the dependence of the spreading speed on the moving speed $\nu$ with small $u_{0}$ and $h_{0}>l^{*}$.

Exmple 3.9**.**

*With a large $\nu$, let the parameters in the system (1.1) be the following

$(h_{0},u_{0},\chi_{1},\chi_{2},\nu,\lambda_{1},\lambda_{2},\mu_{1}.\mu_{2})=(2,\cos(\pi x/2h_{0}),0.2,0.1,2,1,2,1,2)$, the system has tendency of spreading to the half line $\mathbb{R}^{+}$ (see Figure 18).*

Exmple 3.10**.**

With a even smaller moving speed $\nu=0.01$, let other parameters are the same as in Exmple 3.9. The system (1.1) has tendency of spreading and the spreading speed is smaller compared to the system with larger $\nu$ (see Figure 20, 20).

The above simulations indicate that the spreading happens when the moving speed $\nu$ is small and we also can conclude that the asymptotic spreading speed depends on the moving speed $\nu$. Figure 12 and 12 also indicate that when the spreading happens, the spreading speed is a increasing function of moving speed $\nu$.

In the following, we use a table to compare the spreading speeds with different initial solution $u_{0}$, initial habitat $h_{0}$, and the effects of the chemotactic sensitivity $\chi_{1},\chi_{2}$. These simulations indicate the asymptotic spreading speed is independent of the parameters $u_{0},h_{0},\chi_{1},\chi_{2}$.

Exmple 3.11**.**

Let $(h_{0},\chi_{1},\chi_{2},\nu,\lambda_{1},\lambda_{2},\mu_{1}.\mu_{2})=(2,0.2,0.1,0.8,1,2,1,2)$ and $u_{0}=\sigma\cos(\pi x/2h_{0}).$ With different $\sigma$ choices, we have the following spreading speed data.

[TABLE]

The table indicate the spreading speed converges to a constant near $0.7$ with different choices of $\sigma$. The smaller $\sigma$ has a slower spreading speed at the beginning, but converges to a similar constant to other larger choice of $\sigma$ as time increase. For the larger $\sigma$, the initial spreading speed decrease first and then converges to a constant which is independent of the initial $u_{0}$.

In Fisher-KPP free boundary problem, the spreading speed is also independent of the initial habitat when spreading happens. Because of the lack of comparison principle in chemotaxis system, this result is still open. Our simulation indicate this results should be also exist in chemotaxis system.

Exmple 3.12**.**

Let the parameters in the system (1.1) as $(h_{0},\chi_{1},\chi_{2},\nu,\lambda_{1},\lambda_{2},\mu_{1}.\mu_{2})=(2,0,0,\\ 0.8,1,2,1,2)$ and $u_{0}=\sigma\cos(\pi x/2h_{0})$. we have the following spreading speed data.

[TABLE]

Exmple 3.13**.**

*Let the parameters in the system (1.1) as $u_{0}=\cos(\pi x/2h_{0})$ and

$(\chi_{1},\chi_{2},\nu,\lambda_{1},\lambda_{2},\mu_{1}.\mu_{2})$ = $(0.2,0.1,0.8,1,2,1,2)$. With different $h_{0}$ choices, we have the following spreading speeds table:*

[TABLE]

Fisher-KPP free boundary problem [7] is a special case of the chemotaxis free boundary problem with $\chi_{1}=0,\chi_{2}=0$, which are fully investigated and its asymptotic spreading speed $\frac{h(t)}{t}$ is only depending on $a(t,x)$. Our numerical simulations indicate the speed may be independent of the chemotactic sensitivity $\chi_{1},\chi_{2}$.

Exmple 3.14**.**

*Let the parameters in the system (1.1) as $u_{0}=\cos(\pi x/2h_{0})$ and

$(\chi_{1},\chi_{2},\nu,\lambda_{1},\lambda_{2},\mu_{1}.\mu_{2})=(0,0,0.8,1,2,1,2)$. With different $h_{0}$ choices, we have the following spreading speeds table:*

[TABLE]

3.5 Asymptotic behaviors with large chemotactic sensitivities

The long time behaviors of the system is depend on the choices of parameters [32]. In order to guarantee the convergence of the system to the constant $a/b=2$, we need chemotactic sensitivity coefficient to be small enough, and if not, the system may converges to some other constants. The following simulations indicate such result.

Exmple 3.15**.**

Let $h_{0}=2.5>l^{*}=1.11$, $u_{0}=\cos(\pi x/2h_{0})$, and $(\chi_{1},\chi_{2},\nu,\lambda_{1},\lambda_{2},\mu_{1}.\mu_{2})=(0.2,0.1,0.8,2,1,2,1)$. Simulation indicates the system has spreading tendency but does not converge to any constant (see Figure 21).

Exmple 3.16**.**

Compared to Exmple 3.15, fix other parameters and let $(\lambda_{1},\lambda_{2},\mu_{1},\mu_{2})=(2,1,1,2)$, the simulation indicate the system converge to the constant $a/b=2$ (see Figure 23).

The corresponding spreading speed $\frac{h(t)}{t}$ tends to a positive constant 0.69 (see Figure 23),which is similar to Fisher-KPP free boundary problems. It is an evidence that the spreading speed is independent of the ratio of $a/\lambda_{1}$ [28].

Exmple 3.17**.**

Compared to Exmple 3.15, fix other parameters and let $\lambda_{1}<\lambda_{2}$ such that $(\lambda_{1},\lambda_{2},\mu_{1},\mu_{2})=(1,2,2,1)$, we have the following spreading result which does not converge to $a/b=2$.

4 Conclusion

In this paper we construct a numerical scheme to approximate the continuous logistic type chemotaxis system with a free boundary. The scheme’s error estimate, positivity preserving, the monotonicity of the free boundary and stability are investigated. Numerical simulations validate some proved theoretical results such as vanishing-spreading dichotomy, persistency and stability. Compared to Fisher-KPP free boundary problem, the dependence of the vanishing-spreading dichotomy on initial solution $u_{0}$ and initial habitat $h_{0}$ are still open problems both theoretically and numerically. Furthermore, the existence of the asymptotic spreading speed $\frac{h(t)}{t}$ and its dependence on the parameters are also open problems. All these continuous and discrete dynamical questions should be investigated in the future.

Acknowledgement

The authors appreciate Prof. Wenxian Shen at Auburn University for introducing this interesting topic and many valuable and insightful discussions and help.