A novel numerical approach for solving fractional order differential equations using hybrid functions

1.   Introduction

The Keller-Segel equations have been widely used to model chemotaxis. These equations and in particular, the properties of their solutions have been intensively investigated, see for example [9,26,32,36]. For adapted and expanded versions of the Keller-Segel model we refer to [7,13,14,29,32,33,12]. Moreover, in recent literature improved flux-limited Keller-Segel models, taking into account the finite speed of propagation, have been developed in [1]. For surveys and extended reference lists, see, for example [2,3].

Our starting point is the classical kinetic chemotaxis equation [11]. Scaling it with the so called diffusive scaling leads to the Keller-Segel equation, see [11]. In general, the derivation of Keller-Segel type models, including flux-limited diffusion models and Fokker-Planck type models, from underlying kinetic or microscopic models is discussed for example in [2,12,11]. In particular, using moment closure approaches one may obtain macroscopic equations intermediate between kinetic and Keller-Segel equations, see the above mentioned references or [19,28,6]. Linear moment closure models called the P1-model (in the radiative transfer literature) or the Cattaneo equations are in many situations resonable approximations of the kinetic equation [6]. However, a major drawback of these approximations is that they do not guarantee the positivity of the cell density. Thus, in this paper we consider a non-linear closure method to deal with this problem. Additionally, we also introduce half-moment non-linear closures for the kinetic equations and obtain associated hydrodynamic macroscopic equations similar to [17].

To model chemotaxis on a network, the crucial point is to define suitable coupling conditions. In previous work [5,6] coupling conditions for the scalar Keller-Segel equation and kinetic equation have been discussed. An analytical investigation can be found in [10]. These investigations deal with the case of a network problem for a parabolic system of equations. Conservation of mass and continuity of density and chemo-attractant at the nodes are used as coupling conditions. Moreover, in [8,25] the hyperbolic-parabolic Cattaneo equations on a graph are studied. For general work on kinetic equations on networks, for example for traffic flow equations, we refer to [18] or [27].

In the present work we will derive the coupling conditions for nonlinear half-moment models from the coupling condition for the kinetic model. These coupling conditions guarantee on the one hand the conservation of mass through nodes, and on the other hand, they satisfy a positivity condition for the density. In the present case, as in [8,25], the hyperbolic equations require another approach to the coupling conditions as in the Keller-Segel context. However, in the diffusive limit, when the scaling parameter goes to $0$, all coupling conditions converge to those of the Keller-Segel model, i.e. the conservation of mass through nodes and a continuity condition, compare for example [5].

Numerical methods for these nonlinear models have to deal with issues like the numerical moment realizability and (as in the linear case) with a large range of the scaling parameter describing different regimes. Thus, it is desirable to use an asymptotic preserving well-balanced (APWB) scheme. These numerical methods work uniformly with respect to the scaling parameter and one obtains in particular a scheme suitable for computations near the diffusive limit. We refer to [22] for the basic idea of such schemes. Combining this with a numerical scheme to solve the coupling problem we obtain results for the full network problem.

The paper is organized as follows. In section 2 we discuss the kinetic model and its nonlinear hydrodynamic full-and half-moment approximations which are derived using non-linear closure function. Section 3 considers the coupling conditions for the kinetic equations and derived conditions for the macroscopic models. Section 4 and 5 contain details about the numerical schemes on an interval and the numerical treatment of the coupling conditions. Finally, the numerical results for tripod and more general networks are shown in Section 6.

2.   Kinetic and macroscopic models for chemotaxis

2.1.   Kinetic equation

Denote with $f=f(x,t,v)$ the cell distribution at position $x\in\mathbb{R}$, moving with velocity $v\in V\equiv [-1,\space 1]$ at the time $t\in [0,\space T]$, $\rho(x,t)=\int_Vf(x,t,v)dv$ is the macroscopic density of cells, $m=m(x,t)$ is the density of the chemoattractant. Let $\alpha \in \mathbb{R}^+$ denote the rate of turning into the favourable direction given by the chemoattractant, let $\lambda \in \mathbb{R}^+$ be the turning rate of the cells and $D, \gamma_\rho, \gamma_m \in$ $\mathbb{R}^+$ be the diffusivity, production and decay rate of the chemoattractant, respectively. As in [6], we consider the kinetic equation

Here, we have defined

We assume

in order to guarantee positivity of the turning kernel

for all $v$, compare [11]. We introduce a parameter $\epsilon\in [0, \infty)$ and scale the kinetic equation (1) in the following way

Condition (2) turns into $\lambda\geq \epsilon \alpha$ or

This equation will be used as the starting point for the macroscopic models considered in the next section.

The kinetic equation (3) can be transformed using the even-and odd-parities [30]

for positive velocities $(v\geq 0)$. This leads to the system

When $\epsilon\to 0$, the limit of the kinetic equation (3) is the Keller-Segel equation with flux-limited chemoattractant.

The procedure is shortly revisited. Start with the scaled version of the kinetic equation (6). Integrating the equation for $r$ over $v \ge 0$ we obtain

Moreover we observe $r= \frac{\rho}{2} + \mathcal{O}(\epsilon^2)$ and $j = \frac{\alpha \overline{\partial_x m} }{2 \lambda} v \rho - \frac{v}{\lambda}\partial_x r +\mathcal{O}(\epsilon^2)$. Plugging the expansions into (8) gives the diffusion limit of the kinetic equation, the modified or flux-limited Keller-Segel equation

In addition, from the scaled version for kinetic chemotaxis equation (3), using moment closure approaches one may obtain macroscopic, hyperbolic equations intermediate between kinetic and Keller-Segel equations, see [19,29]. There are two kinds of moment models for chemotaxis, full-moment models and half-moment models, that are corresponding to full-and half-moment closure methods respectively.

2.2.   Full-moment non-linear models

Consider the following averaged quantities

where $\rho$ represents the density and $q$ the flow of the cells. The fact that $f(x,t,v)\geq 0$, it leads to the realizability conditions for full-moment models, see [16,37]:

This leads to the equations

In order to close the model (11), a moment closure is specified by using a suitable approximation of $f$ in terms of a function of $\rho$ and $q$ to obtain expressions for $P$. A linear closure with $f(x,t,v) = \frac{1}{2}\rho(x,t) + \frac{3}{2} vq(x,t)$ leads to the P1 model for chemotaxis as in [8,25,6]. However this approach does not guarantee the realizability condition (10) and the cell density $f$ can be negative.

Therefore, one may assume that the cell density $f(x,t,v)$ is approximated by an exponential function

with $a \in \mathbb{R}^+$ and $b \in \mathbb{R}$. This closure function is chosen using an entropy minimization principle, see for example [23]. This leads to

where

guaranteeing the realizability condition. The function $h(u)$ is called the Eddington factor. Without refering to an underlying closure function the function $h(u)$ can also be chosen directly such that the realizability condition is fulfilled. In [35,23] models with different forms for this function have been proposed, for example, the Kershaw closure $h(u)=\frac{1}{3}+\frac{2}{3}u^2$ or the Levermore-Lorentz Eddington factor $h(u)=\frac{1}{3}+\frac{2u^2}{2+\sqrt{4-3u^2}}$. The resulting full-moment macroscopic model for chemotaxis derived by the nonlinear closure (12) reads

Scaling and denoting $q_{\epsilon}:=\frac{q}{\epsilon}$ it can be rewritten as

Under the realizability condition, the model (14) is hyperbolic and converges to the Keller-Segel equations when $\epsilon\to 0$.

2.3.   Half-moment non-linear models

We will construct a macroscopic, half-moment model for (3) by applying a half-moment nonlinear closure, which was introduced in [17], see [6] for a linear version. Consider the following quantities

Using $f(x,t,v)\geq 0$, one obtains the conditions $\rho^{\pm}\geq 0$, $\pm q^{\pm}\geq 0$, $P^{\pm}\geq 0$ and

These conditions again guarantee the existence of a non-negative density for these equations [16,37].

Integrating the kinetic equation (3) leads to the half-moment model

The nonlinear closure for the half-moment models is derived similarly as for the full-moment case. One uses the ansatz

Plugging (18) in (15) leads to

such that

Condition (19) guarantees the realizability condition (16). We note that an explicit maximum-entropy Eddington factor approximating the above closure is obtained by using the Kershaw closure [34]

Finally, the half-moment model for chemotaxis reads

Moreover, by introducing new quantities

the half-moment model can be rewritten as

with

The half-moment model (23) is hyperbolic and has again the Keller-Segel equations as macroscopic diffusive limit when $\epsilon$ goes to $0$.

3.   Asymptotic preserving well-balanced (APWB) schemes for the nonlinear models on an interval

We concentrate in the following on the discretization of the equations for the movements of the cell density. The equation for the chemoattractant $m$ is always solved with a forward central difference scheme as in [5]. We construct an APWB numerical scheme for the nonlinear chemotaxis model following closely the strategies developed in [21] for the linear Cattaneo model and [4,15] for radiation hydrodynamics. Moreover, we refer to [24,20,31] for further details on AP and well-balanced schemes.

3.1.   Framework for APWB scheme

First, let us recall the framework of the APWB scheme in [24,20,21] for the Cattaneo model

with $c(x)=\frac{\alpha}{3}\frac{x}{\sqrt{1+x^2}}$, $a,\lambda,\alpha>0$, $0<\epsilon\leq\frac{\lambda}{\alpha}$. The limit problem as $\epsilon \rightarrow 0$ is

By introducing the diagonal variables $U = \frac12\rho + \frac{\epsilon}{2\sqrt{a}} q_{\epsilon}, V = \frac12\rho -\frac{\epsilon}{2\sqrt{a}} q_{\epsilon}$, the system (25) is rewritten as

Using

the well-balanced Godunov scheme for (26) is given as

with $l: = \frac{\sqrt{a}}{\epsilon}\frac{\Delta t}{\Delta x}$ and $\partial_xm^n_{i+\frac12}: = \frac{m^n_{i+1}-m^n_i}{\Delta x}$. The APWB numerical solution (28) can be rewritten explicitly as

It is easy to check that for small enough values of $\epsilon\leq\sqrt{a}\frac{\Delta t}{\Delta x}\leq 1$, the APWB scheme (29) is positivity preserving for $U$ and $V$, if the following parabolic CFL restriction holds

Furthermore, as mentioned in [21], (28) converges to a centered discretization of the Keller-Segel equation, when $\epsilon\to 0$.

3.2.   APWB scheme for the nonlinear full-moment model

Based on the above considerations and the schemes developed in the above cited references we consider the following relaxation system for (14)

For $\beta\to 0$ the system converges to (14). Note that for any value of $\beta$ the system (30) converges to the Keller-Segel equation with a diffusion coefficient equal to $a$.

System (30) is split into two into separate sub systems. The right hand side treated with a simple implicit numerical scheme and we obtain for $\beta=0$

which is used as the initial condition for the system given by the left hand side. By introducing the new variables

the left hand side of system (30) can be rewritten as

Now, we can apply the APWB schemes presented in 3.1 for each subsystem in (33). Remark that for $a:=\max\Big\vert h\left(\epsilon\frac{q_{\epsilon}}{\rho}\right)\Big\vert$ the quantities in (32) are nonnegative and the scheme converges formally to the Keller-Segel equations as $\epsilon \rightarrow 0$.

Remark 1. We note that under similar conditions as in the linear case, it can be shown, that the scheme preserves positivity in $U$ and $V$. However, this does not mean that we can prove that the density $\rho$ remains positive. In the linear case, this is not an issue to be investigated, since the density $\rho$ in the Cattaneo model is not necessarily positive. In contrast, in the nonlinear case, it is an important open issue to be investigated.

3.3.   APWB scheme for half-moment M1 model

Following the same procedure as for the full-moment model in the previous section, we consider the following relaxation system for the half-moment model (23)

with $h^{\pm}:= h^{\pm}\left(\frac{\epsilon q_{\epsilon}\pm\hat{q}}{\rho\pm\epsilon\hat{\rho}}\right)$ and $P=P\left(\rho,q_{\epsilon},\hat{\rho},\hat{q}\right),\hat{P} = {P}\left(\rho,q_{\epsilon},\hat{\rho},\hat{q}\right)$ defined in section 2.3. To obtain the correct asymptotic behaviour using a splitting scheme, we introduce a parameter $\sigma < \lambda$, and rewrite the system (34) as

Using again an implicit scheme for the right hand side of (35) and setting $\beta=0$, one obtains for the first step of the splitting scheme

This relaxation scheme ensures that system (35) has the correct asymptotic limit as $\epsilon\to 0$. Thus (35) converges to the Keller-Segel equation (9). Moreover, the solution in the relaxation part (36) preserves the realizability condition i.e., $\vert\epsilon q_{\epsilon}^{(1)}\vert\leq \rho^{(1)}$.

The transport part in (35) includes four independent subsystems with (36) as the initial data. Introducing the characteristic variables $(U_1,V_1)$, $(\bar{U}_1,\bar{V}_1)$, $(U_2,V_2)$, $(\bar{U}_2,\bar{V}_2)$,

we use the APWB scheme for the left hand side of (35). With $a:=\max\left(\Big\vert\frac{w}{\rho}\Big\vert\right)$ we have again the positivity for the characteristic variables (37). Moreover, the half-moment model (34) converges formally to the Keller-Segel equation (9) when $\beta\to 0$ and $\epsilon\to 0$.

4.   Coupling conditions for chemotaxis models on networks

In this section we state a hierarchy of coupling conditions for the kinetic and the nonlinear half-moment model proposed in the previous section. In [6] coupling conditions have been derived for the kinetic model and linear macroscopic approximations. These conditions can be used as well for the nonlinear half moment model discussed in the previous section. Finding conditions for the nonlinear full moment model is much harder since the signs of the eigen values can change with the states of the system. First the coupling conditions for the kinetic model are recalled to derive thereof the nonlinear half moment coupling conditions.

4.1.   Coupling condition for kinetic equations [6]

Consider the kinetic equation

As stated in [6], the coupling conditions for (38) should assign a value to all $f(v)$ for $v\in[0,1]$. A possible choice is

where $f_i^+=f_i(v)$ and $f_i^-=f_i(-v)$ for $v\in[0,1]$ and $i = 1,\dots,N$. In order to conserve the total mass in the system the matrix $A\in \mathbb{R}^{N\times N}$ has to fulfill

Moreover, positivity requires $a_{i,j}\geq 0$, $i,j=1,\dots,N$. If all edges are treated equally, the coupling conditions for three edges are

This can be generalized to the cases $N>3$ by choosing $a_{i,j}=\frac{1}{N-1}$ $i\neq j$ and $a_{i,i}=0$.

Remark 2. Proving existence, uniqueness or stability of the kinetic equation with the above coupling conditions on the graph are interesting further issues, which would deserve a deeper analysis.

4.2.   Coupling condition for the nonlinear half-moment model

Since the kinetic coupling conditions (39) are linear and independent of $v$ we obtain directly by integrating the coupling conditions for the half moment model (21)

Recall that all properties of (39) are inherited by the above coupling conditions. For example the total mass in the system is conserved since

In particular, the coupling conditions for the half-moment model yield the positivity of the densities, i.e. $\rho^+_i\geq 0\ \forall i=1,\dots,N$ if $\rho^-_j\geq 0\ \forall j=1,\dots,N$.

4.3.   Numerical coupling of the nonlinear half-moment model

Here only the case of a junction coupling three edges will be treated, but all procedures easily extend to arbitrary junctions. Consider the left hand side of (35) which includes four independent subsystems of eight variables

For the additional numerical quantities we need further coupling conditions. Note that for $i = 1,...,N$ from the relaxation part of (35) we have: $z_i = q_i,\hat{z}_i = \hat{q}_i$ and $w_i = P_i,\hat{w}_i = \hat{P}_i$: Thus, conditions for the states $z$ and $\hat z$ are found from those prescribed for $q^\mp$. To obtain conditions for $w$ and $\hat w$ we have to find conditions for $P_i^\pm$. These are easily derived from the kinetic coupling conditions (39). Multiply (39) by $v^2$ and integrate over positive velocities to obtain

Inserting (22) into (40), we express the coupling conditions in the variables of (41)

and with (42) we obtain

This leads to

Since on each edge there are four waves with positive speed, we need exactly twelve equations (43), (45) for three edges, as provided.

Remark 3. We have not investigated the positivity preserving property of the scheme for the density $\rho$, see Remark 1 for the case of a single line. It would be even more interesting to investigate this property on the graph taking into account the discretization of the coupling conditions.

5.   Numerical tests

In this section we investigate the proposed models in several numerical test cases for different values of $\epsilon$. For simplicity we use for the nonlinear half-moment model the explicit Kershaw closure, i.e. $h^{\pm}$ given by (20). In all the considered examples we will use the following values for the parameters $\lambda =\alpha =1$, $D = 1$, $\gamma_{\rho}=1$ and $\gamma_m=0.1$. The spatial resolution is $\Delta x = 0.02$ and the time step is chosen according to the CFL condition. At end points where no coupling conditions are imposed zero Neumann boundary conditions are applied. For the kinetic model we discretize the velocity space $V=[-1,1]$ with $N_v=50$ cells.

5.1.   Numerical solutions on an interval

Test 1. We consider the interval $[0,2]$ and $\Delta x=0.005$. As initial conditions for the kinetic equation we use

with

The remaining initial values for the hydrodynamic equations can be derived from the kinetic initial condition. At $t=0$ no chemoattractant $m(x,0)=0$ is present.

In the Figures 2 and 3 the densities at $t=0.2$ for the values $\epsilon=1,\ 0.5,\ 0.1$ and $\epsilon=\ 10^{-6}$ are shown. For the kinetic model we observe that diffusion depends strongly on the value of $\epsilon$ as expected. The closest approximation to the kinetic equations are the half moment models, in particular, the nonlinear half-moment model. For the half moment $P_1$ model for $\epsilon=1$ and $\epsilon=0.5$ one can clearly observe the four waves generated by the advective part. In the full moment P1 model only two waves are used to approximate the kinetic solution. The flux-limited Keller-Segel equation is evolving too fast for large $\epsilon$. In general, the half-moment models are clearly superior to the full moment models and the nonlinear models provide better approximations than the linear models. For small $\epsilon$ all models converge to the solution of the Keller-Segel equation.

Test 2. In the second test we investigate the positivity preserving of the nonlinear moment models which is not guaranteed in the case of the linear models. We use non-equilibrium initial conditions

for the kinetic equation (3). $a$ and $b$ are chosen such that the initial data of the half-moment models are given by

with

The initial conditions for the full-moment models are

The solutions are plotted in Figure 4, for $\Delta x = 0.02$ at $t=0.5$. Even though the initial conditions satisfy the realizability conditions $\rho^{\pm}\geq\vert q^{\pm}\vert$ (for half-moment models) and $\rho\geq\vert q\vert$ (for full-moment models), only the nonlinear full-and half-moment models preserve the positivity of the cell density.

5.2.   Numerical solutions on a tripod network

This test case focuses on the coupling conditions for the non-linear half moment model proposed in section 4. We study a junction connecting three outgoing edges, as depicted in Figure 1. The results of the nonlinear half-moment model with the coupling conditions derived in section 4 are compared with the results of the kinetic, linear half and full moment and Keller-Segel model with coupling conditions described in [6]. On each of the edges we consider the interval $[0,1]$. As initial conditions we choose

which is consistent with the following values for the kinetic equation

All other quantities are initially zero.

In Figures 5-8 the densities at $t=0.2$ for different values of $\epsilon$ are shown. The results are similar to those on a single interval. The nonlinear half moment model is again the best approximation to the kinetic result. As the value of $\epsilon$ decreases all solutions approach the solution of the Keller-Segel equation.

5.3.   Numerical solutions on a larger network

In this last test case we consider a larger network of $31$ edges and $23$ nodes as shown in Figure 9. The length of the short edges is $0.5$, for the longer ones we have $1$ and $\sqrt{2}$ respectively. Note that this network does not only contain three way junctions, but also nodes connecting up to five edges. At the open ends of the network Dirichlet boundary conditions are imposed. For the kinetic model we prescribe $f(0,v,t)=\frac12$ for $v>0$. The boundary values for the other models can be derived thereof. As initial conditions all values are set to zero except the density in the edges at the outer boundaries, which is set to $1$. For the numerical computation we used $15$, $30$ and $42$ cells for the spatial discretization of the edges, respectively. The velocity space in the kinetic model is resolved with $30$ points. The scaling parameter is chosen as $\epsilon=1$.

In Figure 9 the density at time $t=5$ for all four models is shown. There is an inflow from both sides of the network. As observed in the previous tests, the states of the Keller-Segel model propagate faster, such that the network is filled at an earlier time. The solution of the kinetic model and the one of the half moment models almost coincide.

In Figure 10 the evolution of the total mass in the network up to $T=30$ is shown. As before, we observe that the Keller-Segel model fills the network much faster, than the other three. The values of the half moment and of the kinetic model almost coincide. Concerning the computation times, the least expensive model is the full moment $P_1$ model with approximately $15 \%$ of the computation time of the kinetic model. The Keller-Segel model needs $18\%$, the half moment $P_1$ model needs $22\%$ and the half moment $M_1$ model needs $23 \%$ of the kinetic computation time. The overall CPU time for the kinetic computation with $30$ spatial and $30$ velocity cells per edge has been approximately $210$ min on an Intel Core i5-3230M with 2.6 GHz.

6.   Conclusions

In this paper we developed nonlinear half-moment models for chemotaxis models which approximate very efficiently the solution of a network problem governed by kinetic equations. An APWB scheme is investigated for these nonlinear moment models. We derive from coupling conditions for the kinetic models corresponding conditions for the macroscopic models. Note that this procedure also might be applied to other kinetic models. In the numerical tests we investigated the dynamics for different coupling conditions and for varying values of $\epsilon$. A simulation on a larger network showed the applicability to more complicated structures and differences in behavior of macroscopic and kinetic models. Finally, we remark that there are several open issues. Most important would be analytical investigations of the kinetic network problem and a thorough investigation of the positivity preserving property of the APWB scheme for the density $\rho$ in the nonlinear half moment model on the graph.