# A Hybrid Finite Volume Method for Advection Equations and Its Applications in Population Dynamics

Chang Yang, Léon Matar Tine

<figure><img src="image_1.png" /></figure>

## To cite this version:

Chang Yang, Léon Matar Tine. A Hybrid Finite Volume Method for Advection Equations and Its Applications in Population Dynamics. Numerical Methods for Partial Differential Equations, 2017, 33 (4), pp.1114 - 1142. 10.1002/num.22134. hal-01646859

HAL Id: hal-01646859

https://hal.science/hal-01646859

Submitted on 15 Feb 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# A HYBRID FINITE VOLUME METHOD FOR ADVECTION EQUATIONS AND ITS APPLICATIONS IN POPULATION DYNAMICS

CHANG YANG AND LÉON MATAR TINE

**ABSTRACT.** We present in this paper a very adapted finite volume numerical scheme for transport type-equation. The scheme is an hybrid one combining an anti-dissipative method with down-winding approach for the flux [8, 6] and an high accurate method as the WENO5 one [13]. The main goal is to construct a scheme able to capture in exact way the numerical solution of transport type-equation without artifact like numerical diffusion or without "stairs" like oscillations and this for any regular or discontinuous initial distribution. This kind of numerical hybrid scheme is very suitable when properties on the long term asymptotic behavior of the solution are of central importance in the modeling what is often the case in context of population dynamics where the final distribution of the considered population and its mass preservation relation are required for prediction.

**KEYWORDS.** Discontinuity detector; WENO scheme; anti-diffusive method; population dynamics.

**AMS SUBJECT CLASSIFICATIONS.** 65M08, 65M12, 92C40.

## CONTENTS

<table><tr><td>1. Introduction</td><td>1</td></tr><tr><td>2. Biological models</td><td>2</td></tr><tr><td>3. A hybrid finite volume method for advection equations</td><td>4</td></tr><tr><td>4. Numerical results</td><td>13</td></tr><tr><td>5. Conclusion and perspective</td><td>22</td></tr><tr><td>Acknowledgment</td><td>26</td></tr><tr><td>References</td><td>26</td></tr></table>

## 1. INTRODUCTION

In this paper we are interested in finite volume numerical simulations of PDEs of hyperbolic type with transport term. More precisely we are looking to correct the numerical dissipation which appears in the simulations of the asymptotic profile of such PDEs. Indeed this numerical dissipation is an artifact that is inherent to most of existing numerical schemes even for high-order ones. So, as reported first in [3] and confirmed in [8] for the Lifshitz-Slyozov equation which is of transport type, capturing numerically the exact asymptotic profile for transport equation is a real challenge because numerical dissipation smooths out the fronts and leads to an artificial profile. In the context of biology modeling, specially in population dynamics, most of existing models [20, 10, 4] contain at least a transport part which takes into account the growth of the considered population, it is crucial to recover the exact asymptotic profile in order to predict the behavior of the population or to estimate some parameters for instance its growth, division or death rates by using measures that are based on the profile. Some authors address the question consisting to correct this inherent numerical dissipation by establishing adequate schemes such as the WENO (Weighted Essentially Non-Oscillatory) scheme [13], anti-dissipative scheme ADM (Anti Dissipative Method) [8], anti-diffusive flux correction [29], etc. All these schemes define their numerical fluxes in order to minimize at best the artifacts. The WENO schemes in [13] is a fifth-order finite difference one developed
for the resolution of contact discontinuities and is based on the successful ENO (Essentially Non-Oscillatory) schemes [12, 24, 25]. The WENO schemes are very suitable for problems containing both shocks and smooth flow features. Though, they are relatively weak for the resolution of contact discontinuities. It is the difficult challenge in general for the high order methods to maintain at the same time non-oscillatory property and well contact discontinuities [23]. Of course, to overcome this challenge, many efficient strategies are proposed in [11, 30, 23] but all these attempts are unsuccessful for multi-dimensional problems except for few situations [26].

The ADM scheme is based on the approach called limited downwind scheme and is developed by B. Desprès and F. Lagoutière [7]. It is similar to the class of flux limiters by Sweby [27] which allow to keep well the contact discontinuities with non-linear stability. In linear advection, the strategy in ADM scheme is the same with the ultrabee scheme [21] developed by Roe. With the anti-dissipative flux, the scheme is very successful for keeping contact discontinuities even in scalar system cases. Of course the scheme is only first order accurate and is not really appropriate for reconstructing solutions containing both discontinuities and smooth profiles. Indeed for smooth solutions, the ADM scheme shows some “stair steps” looking like oscillations. For more details one can refers to [8, 7, 2].

For the anti-diffusive flux correction [29], the authors use the approach of [7, 2] for the high order finite difference WENO schemes in [1, 13] to maintain a high order accuracy in smooth regions, a non-oscillatory behavior near discontinuities and a sharp contact discontinuity resolution even in larger time. Contrary to the first-order anti-dissipative scheme in [8, 7, 2], Xu and Shu don't use a simple Euler forward time discretization, they perform a high order Runge Kutta time discretization and they introduce some readjustments at each stage of the Runge Kutta methods in order to maintain a sharp contact discontinuity resolution.

We propose here a hybrid finite volume scheme where the construction of the numerical flux appears as a combination of the WENO and ADM fluxes and where the time discretization is a high order Runge-Kutta one with a slight adaptation. This scheme is very suitable for firstly removing the numerical artifacts but also correct the stair treads appearing in ADM scheme. More precisely we use a discontinuity detector at each grid point and use the WENO order 5 scheme when the solution is regular near the point otherwise we apply the ADM. We choose a WENO reconstruction instead of a linear one because in the case where the solution is continuous and change rapidly it is shown that WENO is better. As a validation of this hybrid scheme we use two test cases. The first one is a classical (academic) test case for transport equation where the initial distribution is considered on one hand as a very oscillatory one as given in [13]; in the other hand, we use in 2D the famous Zalesak' disk test that is given in [14]. The second kind of test is based on population dynamics and polymerization process where we consider a population of cells or polymers growing either by nutrients uptake (for cells) or by gain and lost of monomers (for polymers). This application on population dynamics is of great importance because in many cases some predictions on the numerical behavior of the solution allow to investigate inverse problem of estimating relevant parameters of the considered model. So, having a bad numerical reconstruction induces bad parameters estimation.

The paper is organized as follows. In section 2, we recall the biology context on which we focus our study. In section 3, we detail the derivation of our hybrid scheme which is based on a general conservation laws. The section 4 is devoted to the numerical results and comparison of our hybrid method with WENO and ADM schemes.

## 2. BIOLOGICAL MODELS

In cell biology as in physics of particles, the evolution dynamics of a group of cells or macro-particles in cell culture or in a bath of micro-particles plays a central role in the
understanding and the explanation of some physical and biological behaviors. Often the observed quantities evolve by growth process either by nutrients uptake (the example of the micro-organism *Daphnia* which uses the nutrients to grow [20, sect. 4.3.1], [17, 22]) or by earnings micro-particles by polymerization (for polymers modeling [16]).

Lots of conjectures are based on the observation of these quantities, especially on their evolution dynamics in long term. In modeling point of view this long term evolution dynamic is obtained by the analysis of the asymptotic behavior of the considered quantity.

Following the processes taken into account in the model, the asymptotic behavior can either be dependent or be independent of the initial distribution of the considered group of cells or parasites. Indeed, in the case where the considered population evolves only by growth it is proven in Lifshitz-Slyozov equations that the asymptotic behavior depends on the initial distribution [3, 8]. However, when aggregation process or division process is taken into account, the asymptotic behavior is regularized towards a quasi-universal profile as shown in [8] and then it is independent to the initial distribution.

A crucial point linked to the modeling of these phenomena is the numerical simulation which can lead, following the used scheme, to bad conjectures on the behavior of the model. These bad conjectures are resulted from some numerical artifacts caused by numerical dissipation inherent to some standard schemes [8].

In order to apply the hybrid method proposed later in this paper, we consider the following test case where a size-structured cells population model is taken into account and the cells evolve by gain or loss of micro cells. Let denote by $f(t, x, \xi)$ the size density repartition of cells of size $\xi \ge 0$, located at position $x \in \Omega \subset \mathbb{R}^i$ ($i = 1, 2, 3$) at time $t \ge 0$ where $\Omega$ is a smooth bounded domain. Then the model can be written for all $(t, x, \xi) \in \mathbb{R}_+ \times \Omega \times \mathbb{R}_+$ as follows

$$ (2.1) \qquad \begin{cases} \displaystyle \frac{\partial}{\partial t} f(t, x, \xi) + \frac{\partial}{\partial \xi} ((a(\xi)c(t, x) - b)f(t, x, \xi)) = 0, \\ f(t = 0, x, \xi) = f^0(x, \xi). \end{cases} $$

Where $c(t, x)$ is the concentration of micro-organisms (nutrients) and it follows a diffusion equation of this form:

$$ (2.2) \qquad \left\{ \begin{aligned} & \frac{\partial}{\partial t} (c(t, x) + \int_{0}^{\infty} \xi f(t, x, \xi) d\xi) = \Delta_{x} c(t, x), \quad t \ge 0, x \in \Omega, \\ & \frac{\partial}{\partial \nu} c = \nabla c \cdot \nu = 0, \quad \text{on } \partial\Omega. \end{aligned} \right. $$

For this kind of coupling model (2.1)-(2.2), the kinetic coefficients $a(\xi)$, $b$ are interpreted as the rates at which cells gain or loss nutriments (monomers, micro-organisms).

We assume the micro-cells (or monomers) to follow a diffusion equation as depicted in equation (2.2). We endowed this diffusion equation with a homogeneous Neumann boundary condition where $\nu$ is the outward unit vector at point $x \in \partial\Omega$.

The problem (2.1)-(2.2) is a variant of the very known standard Lifshitz-Slyozov system which models the evolution of a population of macro-particles immersed in a bath of monomers [16].

The analytical study of (2.1)-(2.2) concerning the existence, uniqueness and properties of the solution is rigorously done in [9] and the main result is based on the following hypothesis:

**Hypothesis.** The kinetic coefficients $a, b$ are required to satisfy $b = 1$; $a$ is an increasing function with $a(0) = 0$ and $a(+\infty) = +\infty$; $a \in C^0([0, \infty)) \cap C^1((0, \infty))$ and for any $\xi_0 > 0$ there exists $L_{a,0} > 0$ such that $0 \le a'(\xi) \le L_{a,0}$ for $\xi \ge \xi_0 > 0$.
The initial condition satisfy $c(t=0, x) \in L^\infty(\Omega); f(t=0, x, \xi) \in L^\infty(\Omega; L^1((0, \infty), (1+\xi)d\xi))$.

With this previous hypothesis, the authors in [9] prove the following statement on the well-posedness of (2.1)-(2.2):

**Theorem 2.1.** There exists a weak solution $(c, f)$ of (2.1)-(2.2) with, for any $0 < T < \infty$, $c \in L^\infty((0, T) \times \Omega) \cap L^2(0, T; H^1(\Omega))$, $f \in L^\infty((0, T) \times \Omega; L^1((0, \infty), (1+\xi)d\xi))$, $c \in C^0([0, T]; L^2(\Omega)-\text{weak})$, $f \in C^0([0, T]; L^1(\Omega \times (0, \infty)) - \text{weak})$.

In addition they prove thanks to the Neumann boundary condition, the following mass preservation relation:

$$ (2.3) \qquad \frac{d}{dt} \left[ \int_{\Omega} \int_{0}^{\infty} \xi f(t, x, \xi) \, d\xi \, dx + \int_{\Omega} c(t, x) \, dx \right] = 0. $$

The fact that the space variable $x$ acts as a parameter in the size density repartition function $f$ implies that (2.1) is a transport equation. The study of its asymptotic behavior is numerically very challenging.

Indeed, following the chosen model as in (2.1)-(2.2), one needs in the modeling and simulations to recover the evolution dynamics of the considered population such as the time evolution of the total number of individuals (even in asymptotic time), the total mass of the population or the conservation law fulfilled by the model. For the numerical simulations of the evolution dynamics, a very adapted scheme is required in order to capture in exact way the solution of the system without artifact in order to get the right and essential properties. That's the aim to introduce the following hybrid method.

### 3. A HYBRID FINITE VOLUME METHOD FOR ADVECTION EQUATIONS

**3.1. Anti-dissipative method.** In this section, we consider the following advection equation

$$ (3.1) \qquad \frac{\partial f}{\partial t} + \frac{\partial (Vf)}{\partial x} = 0, \quad t \ge 0. $$

where $V(t, x)$ is a given smooth velocity field. Lets consider a regular mesh, with constant spatial mesh size $\Delta x > 0$: the cells are the intervals $[x_{i-1/2}, x_{i+1/2}]$, $i \in \mathbb{N}$ with $x_{-1/2} = 0$, $x_{i+1/2} = (i+1)\Delta x$, and $x_i$ denotes the midpoint of the cell: $x_i = (i+1/2)\Delta x$. We denote by $f_i^n$ the numerical unknown, which is the approximation of $f_i^n = \frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} f(t^{(n)}, x)dx$, where $t^{(0)} = 0 < t^{(1)} < \dots < t^{(n)} < t^{(n+1)}$ are times discretization with a possible variable time step $\Delta t^{(n)} = t^{(n+1)} - t^{(n)}$. We denote by $V_{i-1/2}^n$ the approximations of the velocity at the cell interfaces: namely, we set

$$ V_{i-1/2}^n = V(t^{(n)}, x_{i-1/2}), \quad n \in \mathbb{N}, i \in \mathbb{N}. $$

The finite volume scheme applied to (3.1) gives the following approximation

$$ (3.2) \qquad f_i^{n+1} = f_i^n - \frac{\Delta t^{(n)}}{\Delta x} (V_{i+1/2}^n f_{i+1/2}^n - V_{i-1/2}^n f_{i-1/2}^n). $$

Then the main task is how to define the range of the interface fluxes in order to reconstruct perfect discontinuities if they exist. In other words, the scheme is desired to take the interface fluxes $f_{i+1/2}^n$ so that there is zero numerical dissipation for contact discontinuities. To get this anti-dissipative property, the flux need to be chosen as close as possible to the downwind value of the numerical unknown. The used strategy is the same as the one established by Desprès and Lagoutière [7] and consists to choose the most downwind flux with respect to stability, consistency and positivity constraints of the scheme.

In order to describe the different constraints, let introduce the following useful notations:
• $\nu^n = \frac{\Delta t^{(n)}}{\Delta x}$,

• $m_{i+1/2}^n = \min(f_i^n, f_{i+1}^n)$, and $M_{i+1/2}^n = \max(f_i^n, f_{i+1}^n)$,

• if $V_{i+1/2}^n, V_{i-1/2}^n > 0$:

$$ 
\begin{aligned}
b_{i+1/2}^n &= \frac{1}{\nu^n V_{i+1/2}^n} (f_i^n - \max(f_i^n, f_{i-1}^n)) + \max(f_i^n, f_{i-1}^n) \\
&= \frac{1}{\nu^n V_{i+1/2}^n} (f_i^n - M_{i-1/2}^n) + M_{i-1/2}^n,
\end{aligned}
$$

$$ 
\begin{aligned}
B_{i+1/2}^n &= \frac{1}{\nu^n V_{i+1/2}^n} (f_i^n - \min(f_i^n, f_{i-1}^n)) + \min(f_i^n, f_{i-1}^n) \\
&= \frac{1}{\nu^n V_{i+1/2}^n} (f_i^n - m_{i-1/2}^n) + m_{i-1/2}^n,
\end{aligned}
$$

$$ \mathcal{B}_{i+1/2}^n = \begin{cases} \min\left(B_{i+1/2}^n, m_{i-1/2}^n \frac{V_{i-1/2}^n}{V_{i+1/2}^n} + \frac{f_i^n}{\nu^n V_{i+1/2}^n}\right), & \text{if } m_{i-1/2}^n \ge 0, \\ B_{i+1/2}^n, & \text{otherwise,} \end{cases} $$

• if $V_{i+1/2}^n, V_{i-1/2}^n < 0$:

$$ 
\begin{aligned}
b_{i-1/2}^n &= \frac{1}{\nu^n |V_{i-1/2}^n|} (f_i^n - \max(f_i^n, f_{i+1}^n)) + \max(f_i^n, f_{i+1}^n) \\
&= \frac{1}{\nu^n |V_{i-1/2}^n|} (f_i^n - M_{i+1/2}^n) + M_{i+1/2}^n,
\end{aligned}
$$

$$ 
\begin{aligned}
B_{i-1/2}^n &= \frac{1}{\nu^n |V_{i-1/2}^n|} (f_i^n - \min(f_i^n, f_{i+1}^n)) + \min(f_i^n, f_{i+1}^n) \\
&= \frac{1}{\nu^n |V_{i-1/2}^n|} (f_i^n - m_{i+1/2}^n) + m_{i+1/2}^n,
\end{aligned}
$$

$$ \mathcal{B}_{i-1/2}^n = \begin{cases} \min\left(B_{i-1/2}^n, m_{i+1/2}^n \frac{|V_{i+1/2}^n|}{|V_{i-1/2}^n|} + \frac{f_i^n}{\nu^n |V_{i-1/2}^n|}\right), & \text{if } m_{i+1/2}^n \ge 0, \\ B_{i-1/2}^n, & \text{otherwise,} \end{cases} $$

• if $V_{i+1/2}^n, V_{i-1/2}^n$ do not have the same sign, we set $b_{i+1/2}^n = \mathcal{B}_{i+1/2}^n = f_i^n$ if $V_{i+1/2}^n > 0$ and $b_{i+1/2}^n = \mathcal{B}_{i+1/2}^n = f_{i+1}^n$ if $V_{i+1/2}^n < 0$.

• $\mu_{i+1/2}^n = \max(m_{i+1/2}^n, b_{i+1/2}^n)$, and $\mathcal{M}_{i+1/2}^n = \min(M_{i+1/2}^n, \mathcal{B}_{i+1/2}^n)$.

### 3.1.1. Stability constraints.

Knowing that the flux is nothing else the approximation of the numerical unknown at the interface, it is sought to define in this section a non-empty interval containing the numerical unknown thanks to the following standard Courant-Friedrichs-Levy (CFL) condition:

$$ (3.3) \qquad 0 \le \frac{\Delta t^{(n)}}{\Delta x} \max_i(|V_{i+1/2}^n|) \le 1. $$

From this CFL condition, we consider the case where the velocity is positive in the cell $i$: $V_{i-1/2}^n > 0$ and $V_{i+1/2}^n > 0$, so we have

$$ \frac{1}{\nu^n V_{i+1/2}^n} - 1 \ge 0 $$

and in first hand we multiply this inequality by the non-negative term $f_i^n - \min(f_i^n, f_{i-1}^n)$ and obtain

$$ \left( \frac{1}{\nu^n V_{i+1/2}^n} - 1 \right) \left( f_i^n - \min(f_i^n, f_{i-1}^n) \right) \ge 0. $$
Then we deduce

$$
(3.4) \qquad \frac{1}{\nu^n V_{i+1/2}^n} \left( f_i^n - \min(f_i^n, f_{i-1}^n) \right) + \min(f_i^n, f_{i-1}^n) \ge f_i^n.
$$

In second hand we also multiply by the negative term $f_i^n - \max(f_i^n, f_{i-1}^n)$ and obtain

$$
\left( \frac{1}{\nu^n V_{i+1/2}^n} - 1 \right) \left( f_i^n - \max(f_i^n, f_{i-1}^n) \right) \le 0,
$$

then we deduce

$$
(3.5) \qquad \frac{1}{\nu^n V_{i+1/2}^n} \left( f_i^n - \max(f_i^n, f_{i-1}^n) \right) + \max(f_i^n, f_{i-1}^n) \le f_i^n.
$$

In last hand, assuming that $f_i^n \ge 0 \forall i$ (the positivity constraint is discussed later) we deduce
the obvious relation

$$
(3.6) \qquad \min(f_i^n, f_{i-1}^n) \frac{V_{i-1/2}^n}{V_{i+1/2}^n} + \frac{f_i^n}{\nu^n V_{i+1/2}^n} \ge \frac{f_i^n}{\nu^n V_{i+1/2}^n} \\
\qquad \qquad \qquad \qquad \ge f_i^n.
$$

Owing to (3.4)-(3.6) and the previous notations, one deduce this first non empty interval
where the adequate flux will be chosen

$$
(3.7) \qquad f_i^n \in [b_{i+1/2}^n, \mathcal{B}_{i+1/2}^n] \neq \emptyset.
$$

For the case where the velocity is locally negative mean $V_{i-1/2}^n < 0$ and $V_{i+1/2}^n < 0$, we perform the same reasoning and obtain the non empty interval $f_i^n \in [b_{i-1/2}^n, \mathcal{B}_{i-1/2}^n] \neq \emptyset$.

3.1.2. *Consistency constraints.* For the consistency of the scheme, we write the very definition which is that the numerical flux between cell $i$ and cell $i+1$ belongs necessary to the interval defined by the numerical solutions $f_i^n$ and $f_{i+1}^n$. So the consistency constraint is written as follows

$$
(3.8) \qquad m_{i+1/2}^n \le f_{i+1/2}^n \le M_{i+1/2}^n, \quad \forall i \text{ in the grid.}
$$

3.1.3. *Positivity constraints.* From a positive initial solution, we need to impose a condition on the numerical fluxes in order to ensure the positivity of the numerical approximation of the solution given by (3.2). For this, we shall assume $V_{i-1/2}^n > 0, V_{i+1/2}^n > 0$ and we address the conditions such that

$$
\begin{align*}
f_i^{n+1} \ge 0 &\implies \nu^n (V_{i+1/2}^n f_{i+1/2}^n - V_{i-1/2}^n f_{i-1/2}^n) \le f_i^n \\
&\implies f_{i+1/2}^n \le f_{i-1/2}^n \frac{V_{i-1/2}^n}{V_{i+1/2}^n} + \frac{f_i^n}{\nu^n V_{i+1/2}^n},
\end{align*}
$$

so, obtaining the later inequality is done by imposing the flux to satisfy

$$
(3.9) \qquad f_{i+1/2}^n \le m_{i-1/2}^n \frac{V_{i-1/2}^n}{V_{i+1/2}^n} + \frac{f_i^n}{\nu^n V_{i+1/2}^n},
$$

because by the consistency constraints we know that $m_{i-1/2}^n \le f_{i-1/2}^n$.

**Proposition 3.1.** Assume that the CFL condition (3.3) be satisfied. Then for any $i$ the interval $[\mu_{i+1/2}^n, \mathcal{M}_{i+1/2}^n]$ is non empty. By choosing the fluxes $f_{i+1/2}^n \in [\mu_{i+1/2}^n, \mathcal{M}_{i+1/2}^n]$ for any $i$, then the following assertions hold:

1) The scheme (3.2) is consistent with (3.1).

2) The scheme (3.2) remains positive for positive initial solution: mean if $f_i^n \ge 0$ for any $i$ then $f_i^{n+1} \ge 0$ too.
3) Thanks to stability and consistency constraints at step n and in case $V_{i-1/2}^n > 0$ and $V_{i+1/2}^n > 0$, the discrete solution satisfies the following bounds:

$$ (3.10) \qquad m_{i-1/2}^n - \nu^n (V_{i+1/2}^n - V_{i-1/2}^n) M_{i-1/2}^n \le f_i^{n+1} \le M_{i-1/2}^n - \nu^n (V_{i+1/2}^n - V_{i-1/2}^n) m_{i-1/2}^n. $$

Similar bounds can be obtained if $V_{i-1/2}^n < 0$ and $V_{i+1/2}^n < 0$. More precisely we have

$$ m_{i+1/2}^n - \nu^n (|V_{i-1/2}^n| - |V_{i+1/2}^n|) M_{i+1/2}^n \le f_i^{n+1} \le M_{i+1/2}^n - \nu^n (|V_{i-1/2}^n| - |V_{i+1/2}^n|) m_{i+1/2}^n. $$

**Proof.** The proof of the proposition is essentially based on the gathering of the results (3.4)–(3.6) and (3.8)–(3.9) for the items 1) and 2).

For the last item, we begin by writing the stability constraint:

$$ \frac{1}{\nu^n V_{i+1/2}^n} (f_i^n - M_{i-1/2}^n) + M_{i-1/2}^n \le f_{i+1/2}^n \le \frac{1}{\nu^n V_{i+1/2}^n} (f_i^n - M_{i-1/2}^n) + m_{i-1/2}^n $$

and use the consistency constraint: $m_{i-1/2}^n \le f_{i-1/2}^n \le M_{i-1/2}^n$ to deduce

$$ \frac{1}{\nu^n V_{i+1/2}^n} (f_i^n - M_{i-1/2}^n) + f_{i-1/2}^n \le f_{i+1/2}^n \le \frac{1}{\nu^n V_{i+1/2}^n} (f_i^n - M_{i-1/2}^n) + f_{i-1/2}^n. $$

Multiplying all the terms by $V_{i+1/2}^n$ and adding $-V_{i-1/2}^n f_{i-1/2}^n$ we obtain

$$ \frac{1}{\nu^n} (f_i^n - M_{i-1/2}^n) + (V_{i+1/2}^n - V_{i-1/2}^n) f_{i-1/2}^n \le V_{i+1/2}^n f_{i+1/2}^n - V_{i-1/2}^n f_{i-1/2}^n \le \frac{1}{\nu^n} (f_i^n - M_{i-1/2}^n) + (V_{i+1/2}^n - V_{i-1/2}^n) f_{i-1/2}^n. $$

So, multiplying all the term by $-1$ and rearranging the terms we deduce

$$ m_{i-1/2}^n - \nu^n (V_{i+1/2}^n - V_{i-1/2}^n) f_{i-1/2}^n \le f_i^n - \nu^n (V_{i+1/2}^n f_{i+1/2}^n - V_{i-1/2}^n f_{i-1/2}^n) \le M_{i-1/2}^n - \nu^n (V_{i+1/2}^n - V_{i-1/2}^n) f_{i-1/2}^n, $$

what achieves the proof of the item 3). The proof of the similar bounds in the case where $V_{i-1/2}^n < 0$ and $V_{i+1/2}^n < 0$ is straightforward.

For more details on the proof related to the case of the non conservative transport equation, one can refer to [8].

**Remark 3.2.** In the previous reasoning, we exclude the case where the velocities on the interfaces cell are of different sign. So in the case where $V_{i+1/2}^n < 0$ and $V_{i-1/2}^n > 0$ then it is obvious that the positivity of the numerical solution $f_i^{n+1}$ is obtained without any time step condition. Nevertheless, the case where $V_{i+1/2}^n > 0$ and $V_{i-1/2}^n < 0$ mean the possible empty of the cell from the two sides. So we choose the upwind fluxes and the numerical solution becomes constant in the cell i. Then for this case the non negativity is ensured under a restricted CFL condition $\nu^n \max_i(|V_{i+1/2}^n|) \le \frac{1}{2}$.

**Remark 3.3.** In this above presentation, the numerical flux is designed for general conservation laws, while the one in [8] is more suitable for the transport equations.

For the anti-dissipative strategy, we define the flux $f_{i+1/2}^n$ by solving the following minimization problem:

To minimize $|f_{i+1/2}^n - f_{i+1}^n|$ under the constraint $f_{i+1/2}^n \in [\mu_{i+1/2}^n, \mathcal{M}_{i+1/2}^n]$
which solution is given for instance in the case $V_{i-1/2}^n > 0$ and $V_{i+1/2}^n > 0$ by

$$ (3.11) \qquad \left\{ \begin{array}{ll} f_{i+1/2}^n = \mu_{i+1/2}^n & \text{if } f_{i+1}^n \le \mu_{i+1/2}^n, \\ f_{i+1/2}^n = f_{i+1}^n & \text{if } \mu_{i+1/2}^n \le f_{i+1}^n \le \mathcal{M}_{i+1/2}^n, \\ f_{i+1/2}^n = \mathcal{M}_{i+1/2}^n & \text{if } f_{i+1}^n \ge \mathcal{M}_{i+1/2}^n. \end{array} \right. $$

This kind of anti-dissipative method is very suitable for discontinuous initial solution, which has been shown in [6]. However, it is not suitable for smooth solution. Indeed, it turns the very smooth solutions into a series of step functions with respect to time evolution [6, 8]. The objective of the remaining part is to find an alternative method such that it keeps the shock sharp near discontinuities while having high accuracy in smooth regions. To the end, we propose the following hybrid method.

**3.2. A hybrid method.** We denote $f_{i+1/2}^A$ the flux computed by the limited Downwind scheme and $f_{i+1/2}^W$ the flux computed by a high accurate Upwind scheme, for instance the fifth order WENO reconstruction [13]. So our desirable flux $f_{i+1/2}^H$ by the hybrid method will just be a convex combination of $f_{i+1/2}^A$ and $f_{i+1/2}^W$, i.e.

$$ (3.12) \qquad f_{i+1/2}^H = \omega_{i+1/2}^A f_{i+1/2}^A + \omega_{i+1/2}^W f_{i+1/2}^W $$

where $\omega_{i+1/2}^A + \omega_{i+1/2}^W = 1$, $\omega_{i+1/2}^A, \omega_{i+1/2}^W \ge 0$. Moreover, it is desirable:

* $\omega_{i+1/2}^A = \mathcal{O}(1)$ and $\omega_{i+1/2}^W = o(1)$, near discontinuities,

* $\omega_{i+1/2}^A = \mathcal{O}(\Delta x^M)$ and $\omega_{i+1/2}^W = \mathcal{O}(1)$, in smooth regions,

where $M$ is a large enough number that causes $f_{i+1/2}^W$ to be the dominant term in smooth regions.

According to these considerations, we then need to identify the smoothness of solution. We consider a similar smoothness measurement, which was first proposed in [5]

$$ (3.13) \qquad e_i = \frac{|\hat{f}_i - f_i|}{D}, $$

where $D$ is a scaling value given as

$$ (3.14) \qquad D = \max_i f_i - \min_i f_i. $$

Figure 1 shows the interpolated value, i.e. $\hat{f}_i$ is obtained by the fourth-order interpolation. Note that the point $f_i$ itself, as shown in Figure 1, is not included in the interpolation.

$$ (3.15) \qquad \hat{f}_i = \frac{1}{6}(-f_{i-2} + 4f_{i-1} + 4f_{i+1} - f_{i+2}). $$

Now we need to know whether the corresponding cell $i$ is in the smooth region. For this, let us study the smoothness indicator for the solution of transport equation with initial condition as a step function. We consider the simplest case, with only a transition point

$$ f_j = \begin{cases} 1, & j \le 0, \\ \alpha, & j = 1, \\ 0, & j \ge 2. \end{cases} $$

The smoothness indicator as a function of $\alpha$ is depicted in Figure 2. We see $e_j, j = 0, 1, 2$ are intersected. To keep the contact discontinuity of solution, we need to the three cells $j = 0, 1, 2$ in non-smooth region, or at least two cells in non-smooth region. Thus a necessary condition is $e_j \le \beta$, with the criteria $\beta \le 0.1$ (see the dotted line in Figure 2).
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Interpolation error between the node value $f_i$ and interpolated value $\hat{f}_i$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIGURE 2. The smoothness indicator as a function of $\alpha$.</figcaption></figure>

Then we choose the smoothness measurement at the cell interface by an upwind way

$$ (3.16) \qquad e_{i+1/2} = \begin{cases} e_i, & \text{if } V_{i+1/2} \ge 0, \\ e_{i+1}, & \text{else.} \end{cases} $$

Finally, we consider smooth weights defined in [0, 1], antisymmetric with respect to $\beta$, and vary rapidly around $\beta$ (controlled by the parameter $c$). Such weights can be defined by hyperbolic tangent function as follows

$$ (3.17) \qquad \omega_{i+1/2}^{A} = \frac{1 - \tanh(-c(e_{i+1/2} - \beta))}{2}, \qquad \omega_{i+1/2}^{W} = \frac{1 + \tanh(-c(e_{i+1/2} - \beta))}{2}. $$

It is clear that $\omega_{i+1/2}^A + \omega_{i+1/2}^W = 1$ and $0 \le \omega_{i+1/2}^A, \omega_{i+1/2}^W \le 1$. An example of weight $\omega_{i+1/2}^W$ is illustrated in Figure 3 with parameters $\beta = 0.09, c = 150$.

Let recall that parameters $\beta$ and $c$ have a none negligible influence in the behavior of the numerical solution as described in (3.17).

In smooth regions, using a Taylor expansion analysis, we have:

$$ e_{i+1/2} = \mathcal{O}(\Delta x^4). $$

So the weights are approximately equal to:

$$ \omega_{i+1/2}^{A} \approx 0, \qquad \omega_{i+1/2}^{W} \approx 1. $$

Thus, the Upwind flux with fifth-order WENO reconstruction is activated, i.e.

$$ f_{i+1/2}^{H} \approx f_{i+1/2}^{W}. $$
<figure><img src="image_4.png" /><figcaption>FIGURE 3. The weight $\omega_{i+1/2}^W$ with parameters $\beta = 0.09$, $c = 150$.</figcaption></figure>

For a discontinuity, since the size of the discontinuity does not change as $\Delta x \to 0$, one can conclude

$$e_{i+1/2} = \mathcal{O}(1).$$

So the weights are approximately equal to:

$$\omega_{i+1/2}^A \approx 1, \quad \omega_{i+1/2}^W \approx 0.$$

Thus, the limited Downwind flux is activated, i.e.

$$f_{i+1/2}^H \approx f_{i+1/2}^A.$$

**Remark 3.4.** In the original method [5], the parameter $D$ depended on the index $i$ and it was a convex combination of the global and local scales, i.e.

$$D_i = c_s S_g + (1 - c_s) S_l,$$

where $S_g = \max_i f_i - \min_i f_i$, $i \in \{1, \dots, i_{\max}\}$, $S_l = \max_i f_i - \min_i f_i$, $i \in \{i-2, \dots, i+2\}$. The parameter $c_s$ was chosen as 0.1 or 0.01 in [5]. This choice can highlight the small jumps in the solution. However, in the paper we focus on the major jumps in the solution. Therefore, the global scale $S_g$ seems more appropriate.

3.3. **Time discretization.** The second-order TVD Runge-Kutta time discretization [25] is given by

$$f^{(1)} = f^n + \Delta t L(f^n),$$

$$f^{n+1} = \frac{1}{2} f^n + \frac{1}{2} f^{(1)} + \frac{1}{2} \Delta t L(f^{(1)}),$$

where $L$ is the spatial discretization of $-\frac{\partial(Vf)}{\partial x}$. This is equivalent to

$$f^{(1)} = f^n + \Delta t L(f^n),$$

$$f^{n+1} = f^n + \frac{1}{2} \Delta t L(f^n) + \frac{1}{2} \Delta t L(f^{(1)}).$$

In order to maintain the moving travelling wave solutions for piecewise constant functions containing contact discontinuities as in the Euler forward case, we modify the previous scheme as

$$(3.18) \qquad f^{(1)} = f^n + \Delta t L(f^n),$$

$$(3.19) \qquad f^{n+1} = f^n + \frac{1}{2} \Delta t L(f^n) + \frac{1}{2} \Delta t L'(f^{(1)}),$$
where the operator $L$ is defined for $f^n$ by

$$ (3.20) \qquad L(f^n)_i = -\nu^n (V_{i+1/2} f_{i+1/2}^n - V_{i-1/2} f_{i-1/2}^n) $$

In case where $L'(f^{(1)})$ is chosen equal to $L(f^n)$, it is clear that this time discretization is equivalent to the Euler forward scheme. For example, for the transport equation

$$ \partial_t f + a \partial_x f = 0, $$

when CFL number is equal to 0.4, i.e. $a \frac{\Delta t}{\Delta x} = 0.4$, and the initial condition is a step function

$$ f_j^0 = \begin{cases} 1, & j \le 0, \\ 0, & \text{otherwise.} \end{cases} $$

Moreover, let us consider only the anti-dissipative method described in Section 3.1. The sequence of the numerical solution for the first few time steps are:

after one-time step

$$ f_j^1 = \begin{cases} 1, & j \le 0, \\ 0.4, & j = 1, \\ 0, & \text{otherwise.} \end{cases} $$

after two-time steps

$$ f_j^2 = \begin{cases} 1, & j \le 0, \\ 0.8, & j = 1, \\ 0, & \text{otherwise.} \end{cases} $$

after three-time steps

$$ f_j^3 = \begin{cases} 1, & j \le 1, \\ 0.2, & j = 2, \\ 0, & \text{otherwise.} \end{cases} $$

after four-time steps

$$ f_j^4 = \begin{cases} 1, & j \le 1, \\ 0.6, & j = 2, \\ 0, & \text{otherwise.} \end{cases} $$

after five-time steps

$$ f_j^5 = \begin{cases} 1, & j \le 2, \\ 0, & \text{otherwise.} \end{cases} $$

Clearly $f_j^5 = f_{j-2}^0$, hence the numerical solution repeats itself after four-time steps with a shift. If $a \frac{\Delta t}{\Delta x}$ is not a rational number, the numerical solution will not exactly repeat after finite time steps, but the number of transition points will not be larger than one.

The spatial discretization of the anti-dissipative method is only a first order method [6], thus a second time discretization will not improve accuracy. However, it is necessary to use high order time discretization method when a high order spatial discretization is used. So, in order to implement with our hybrid spatial discretization described in Section 3.2, we use the following discretization for the operator $L$

$$ (3.21) \qquad L(f^n)_i = -\nu^n (V_{i+1/2} f_{i+1/2}^{H,n} - V_{i-1/2} f_{i-1/2}^{H,n}), $$

$$ (3.22) \qquad L'(f^{(1)})_i = -\nu^n (V_{i+1/2} f_{i+1/2}^{H,(1)} - V_{i-1/2} f_{i-1/2}^{H,(1)}), $$
where

$$ (3.23) \qquad f_{i+1/2}^{H,n} = \omega_{i+1/2}^{A,n} f_{i+1/2}^{A,n} + \omega_{i+1/2}^{W,n} f_{i+1/2}^{W,n}, $$

$$ (3.24) \qquad f_{i+1/2}^{H,(1)} = \omega_{i+1/2}^{A,(1)} f_{i+1/2}^{A,n} + \omega_{i+1/2}^{W,(1)} f_{i+1/2}^{W,(1)}. $$

We remark that, in the hybrid $f_{i+1/2}^{H,(1)}$, the anti-dissipative flux is equal to $f_{i+1/2}^{A,n}$ and not to $f_{i+1/2}^{A,(1)}$. So, in the case $\omega_{i+1/2}^{A,n} = \omega_{i+1/2}^{A,(1)} = 1$, the formulae (3.20) for $f^n, f^{(1)}$ and (3.21)-(3.22) are equivalent. Therefore, the second order time discretization is given by (3.18)-(3.19), (3.21)-(3.24).

The third-order TVD Runge-Kutta method in [25] has the following form:

$$ f^{(1)} = f^n + \Delta t L(f^n), $$

$$ f^{(2)} = \frac{3}{4} f^n + \frac{1}{4} f^{(1)} + \frac{1}{4} \Delta t L(f^{(1)}), $$

$$ f^{n+1} = \frac{1}{3} f^n + \frac{2}{3} f^{(2)} + \frac{2}{3} \Delta t L(f^{(2)}), $$

which is equivalent to

$$ f^{(1)} = f^n + \Delta t L(f^n), $$

$$ f^{(2)} = f^n + \frac{1}{4} \Delta t L(f^n) + \frac{1}{4} \Delta t L(f^{(1)}), $$

$$ f^{n+1} = f^n + \frac{1}{6} \Delta t L(f^n) + \frac{1}{6} \Delta t L(f^{(1)}) + \frac{2}{3} \Delta t L(f^{(2)}). $$

In the same strategy, as we did in the second order time discretization, we modify the previous scheme as

$$ (3.25) \qquad f^{(1)} = f^n + \Delta t L(f^n), $$

$$ (3.26) \qquad f^{(2)} = f^n + \frac{1}{4} \Delta t L(f^n) + \frac{1}{4} \Delta t L'(f^{(1)}), $$

$$ (3.27) \qquad f^{n+1} = f^n + \frac{1}{6} \Delta t L(f^n) + \frac{1}{6} \Delta t L'(f^{(1)}) + \frac{2}{3} \Delta t L''(f^{(2)}). $$

The operations $L, L', L''$ are defined as

$$ (3.28) \qquad L(f^n)_i = -\nu^n (V_{i+1/2} f_{i+1/2}^{H,n} - V_{i-1/2} f_{i-1/2}^{H,n}), $$

$$ (3.29) \qquad L'(f^{(1)})_i = -\nu^n (V_{i+1/2} f_{i+1/2}^{H,(1)} - V_{i-1/2} f_{i-1/2}^{H,(1)}), $$

$$ (3.30) \qquad L''(f^{(2)})_i = -\nu^n (V_{i+1/2} f_{i+1/2}^{H,(2)} - V_{i-1/2} f_{i-1/2}^{H,(2)}), $$

where

$$ (3.31) \qquad f_{i+1/2}^{H,n} = \omega_{i+1/2}^{A,n} f_{i+1/2}^{A,n} + \omega_{i+1/2}^{W,n} f_{i+1/2}^{W,n}, $$

$$ (3.32) \qquad f_{i+1/2}^{H,(1)} = \omega_{i+1/2}^{A,(1)} f_{i+1/2}^{A,n} + \omega_{i+1/2}^{W,(1)} f_{i+1/2}^{W,(1)}. $$

$$ (3.33) \qquad f_{i+1/2}^{H,(2)} = \omega_{i+1/2}^{A,(2)} f_{i+1/2}^{A,n} + \omega_{i+1/2}^{W,(2)} f_{i+1/2}^{W,(2)}. $$
Therefore, the third order time discretization is given by (3.25)-(3.33). Here also be carful that in the expressions $f_{i+1/2}^{H,(1)}$ and $f_{i+1/2}^{H,(2)}$ the anti-dissipative fluxes are both equal to $f_{i+1/2}^{A,(n)}$.

## 4. NUMERICAL RESULTS

In this section, we will present several numerical results to depict the behaviors of our hybrid method and its applications in population dynamics.

### 4.1. The classical numerical tests.
Here, we first compare different methods for transport equation with mixed initial condition. Then we perform the classical tests with the 1D free transport equation and the 2D rotation equation to verify the convergence of our hybrid method for the regular and irregular initial data. In the sequel, the modified third order Runge-Kutta method is used for time discretization.

#### 4.1.1. Comparison of different methods.
To compare different methods, we consider the free transport equation

$$ (4.1) \qquad \frac{\partial f}{\partial t} + \frac{\partial f}{\partial x} = 0, \quad x \in [-1, 1], \quad t \ge 0, $$

with a very oscillatory initial condition, given by [13],

$$ (4.2) \qquad f(0,x) = \begin{cases} f_1(x) = \frac{1}{6}[G(x, z - \delta) + G(x, z - \delta) + 4G(x, z)], & \text{if } -0.8 \le x \le -0.6, \\ f_2(x) = 1, & \text{if } -0.4 \le x \le -0.2, \\ f_3(x) = 1 - |10(x - 0.1)|, & \text{if } 0 \le x \le 0.2, \\ f_4(x) = \frac{1}{6}[F(x, z - \delta) + F(x, z - \delta) + 4F(x, z)], & \text{if } 0.4 \le x \le 0.6, \\ 0, & \text{otherwise.} \end{cases} $$

where $G(x, z) = \exp(-\beta(x-z)^2)$, $F(x, a) = \max((1-\alpha^2(x-a)^2)^{1/2}, 0)$ with $a = 0.5$, $z = -0.7$, $\delta = 0.005$, $\alpha = 10$ and $\beta = (\log 2)/36\delta^2$. Moreover, the periodic boundary conditions is considered.

The numerical results is presented in Figure 4. We first see that the WENO method is well adapted for the smooth regions of solution (see the Gaussian function), while it becomes significantly diffuse near contact discontinuity (see the step function). For the refined mesh ($n_x = 400$), we can still observe clearly this diffusion near the step function. The anti-dissipative method adapts perfectly for the step function, however it turns the smooth solution into a stair form, which can not be regularized by refining mesh. The Xu-Shu's method was shown [29] to be a successful method for capturing contact discontinuities. However, by performing the same test with their method, we find it keeps well the contact discontinuity, but it does not adapt very well for the smooth solutions. We see, for the Gaussian function, Xu-Shu's method is better than the anti-dissipative method, but some stair forms can still be observed. Finally, with the hybrid method, both the contact discontinuity and the smooth functions are adapted well on two different meshes.

#### 4.1.2. Convergence for the 1D transport equation.
Let us first consider a smooth solution, where the initial condition is chosen as

$$ f(0, x) = \sin(\pi x), \quad x \in [-1, 1]. $$

The numerical error for different methods is presented in Table 1. We first notice that the anti-dissipative method has only first order of convergence rate. At contrast, the WENO method and the hybrid method have both the fifth order convergence rate. Thus they are much more precise than the anti-dissipative method. Moreover, the hybrid method does not disturb at all the precision of numerical results in smooth solution case.
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Numerical solutions for 1D transport equation (4.1) with initial data (4.2) : The left column for $n_x = 200$ while the right column for $n_x = 400$. $\Delta t$ is chosen such that CFL= 0.2. The final time is $T_{end} = 8$.</figcaption></figure>
<table><caption>TABLE 1. 1D transport equation: Error in L₁-norm and order of convergence r. The final time is $T_{end} = 8$.</caption><thead><tr><th rowspan="2">n<sub>x</sub></th><th colspan="2">200</th><th colspan="2">400</th><th colspan="2">800</th></tr><tr><th>||·||<sub>1</sub></th><th>r</th><th>||·||<sub>1</sub></th><th>r</th><th>||·||<sub>1</sub></th><th>r</th></tr></thead><tbody><tr><td>Anti-dissipative method</td><td>4.31e-2</td><td>1.07</td><td>2.06e-2</td><td>1.07</td><td>1.06e-2</td><td>0.96</td></tr><tr><td>WENO scheme</td><td>1.14e-7</td><td>5.00</td><td>3.57e-9</td><td>5.00</td><td>1.12e-10</td><td>4.99</td></tr><tr><td>Hybrid method</td><td>1.14e-7</td><td>5.00</td><td>3.57e-9</td><td>5.00</td><td>1.13e-10</td><td>4.99</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>FIGURE 5. Numerical solutions for 1D transport equation (4.1) with initial data (4.3): The left column for $T_{end} = 10$ while the right column for $T_{end} = 100$. $n_x = 100$ and $\Delta t$ is chosen such that CFL= 0.2.</figcaption></figure>

We next consider a step function as follows

$$ (4.3) \qquad f(0, x) = \begin{cases} 1, & \text{for } -0.5 \le x \le 0.5, \\ 0, & \text{otherwise.} \end{cases} $$

The results as $T_{end} = 10$ (after 5-time periods) and at $T_{end} = 100$ are shown in Figure 5. We can see the classical fifth order WENO method progressively smears the contact discontinuities, and gets more severely with larger time. While the hybrid method has a sharp resolution for the contact discontinuities, which does not deteriorate noticeably with larger time.

Moreover, the hybrid method does not significantly increase the computational cost. Indeed, in average the WENO scheme and the ADM method take 74 operations and 27 operations respectively. While the hybrid method consists of these two methods (103 operations) and in addition the computation of smoothness indicator (27 operations). So in total the hybrid method takes approximately 130 operations, about 1.76 times of the WENO scheme or 4.81 times of the ADM method. However, considering its important features, we think the additional computational cost of our hybrid method is acceptable.

4.1.3. 2D rotation equation. We can directly extend the hybrid method to two dimensional cases. In fact, the numerical flux can be computed dimension by dimension. For example, in the x direction, the numerical flux can be computed by

$$ (4.4) \qquad f_{i+1/2,j}^{H} = \omega_{i+1/2,j}^{A} f_{i+1/2,j}^{A} + \omega_{i+1/2,j}^{W} f_{i+1/2,j}^{W} $$
where $\omega_{i+1/2,j}^A + \omega_{i+1/2,j}^W = 1$, $\omega_{i+1/2,j}^A, \omega_{i+1/2,j}^W \ge 0$. Then we choose the smoothness measurement at the cell interface by an upwind way

$$ (4.5) \qquad e_{i+1/2,j} = \begin{cases} e_{i,j}^x, & \text{if } V_{i+1/2,j} \ge 0, \\ e_{i+1,j}^x, & \text{else,} \end{cases} $$

where the smoothness measurement $e_{i,j}^x$ is the same as in (3.13)-(3.15) for $j$ fixed. The weights have forms

$$ (4.6) \qquad \omega_{i+1/2,j}^A = \frac{1 - \tanh(-c(e_{i+1/2,j} - \beta))}{2}, \qquad \omega_{i+1/2,j}^W = \frac{1 + \tanh(-c(e_{i+1/2,j} - \beta))}{2}. $$

Similarly, we can compute the numerical flux $f_{i,j+1/2}^H$ in $y$ direction.

To test the efficiency of the hybrid method, we use the famous Zalesak' disk test. The governed equation is the 2D rotation equation

$$ (4.7) \qquad \frac{\partial f}{\partial t} + y \frac{\partial f}{\partial x} - x \frac{\partial f}{\partial y} = 0, \quad x \times y \in [-1, 1] \times [-1, 1], \quad t \ge 0. $$

The initial solution is [14]

$$ (4.8) \qquad f(0, x) = \begin{cases} 1, & \text{if } \sqrt{x^2 + (y - 0.5)^2} \le r_0 \text{ \&} \& (|x| \ge 0.025 \|y\ge 0.75), \\ 1 - \frac{1}{r_0} \sqrt{x^2 + (y + 0.5)^2}, & \text{if } \sqrt{x^2 + (y + 0.5)^2} \le r_0, \\ \frac{1 + \cos \pi \frac{1}{r_0} \sqrt{(x+0.5)^2 + y^2}}{4}, & \text{if } \sqrt{(x+0.5)^2 + y^2} \le r_0, \\ 0, & \text{otherwise,} \end{cases} $$

where the radius $r_0 = 0.3$.

We illustrate the numerical results in Figure 6. We first notice that the initial condition consists of the Zalesak' disk, the conical body and the peak of the hump. Then after one period, the WENO method preserves well the conical body and the peak of the hump, however a clear diffusion appears in the Zalesak' disk. The anti-dissipative method keeps well the shape of the Zalesak' disk, but it destroys completely the others two objects. Finally, we observe that the hybrid method performs well for all these three objects. Therefore, the hybrid method is suitable for both the smooth solution and the irregular solution for the advection equations.

## 4.2. Biological numerical tests.
In this subsection, we will focus on numerical simulations for the polymerization/depolymerization type model. The objective here is to highlight the good performance of our hybrid method on the long term asymptotic behavior of the solution.

### 4.2.1. Space-homogeneous polymerization/depolymerization type model.
Here we are interested on the numerical behavior of the standard Lifshitz-Slyozov equations [16] that we refer as homogeneous in space in comparison to the system (2.1)-(2.2).

This standard Lifshitz-Slyozov equations can be interpreted in the point of view of population dynamics as a model describing a population of cells evolving only by nutrients uptake (without birth, death and division) where the nutrients are characterized by their concentration $c(t)$ which fulfills a mass preservation equation as follows

$$ (4.9) \qquad \begin{cases} \frac{\partial}{\partial t} f(t, \xi) + \frac{\partial}{\partial \xi} ((\xi^{1/3} c(t) - 1) f(t, \xi)) = 0, \\ c(t) + \int_{0}^{\infty} \xi f(t, \xi) d\xi = \rho, \end{cases} $$
<figure><img src="image_7.png" /><figcaption>FIGURE 6. 2D rotation equation: Plot solutions of the linear equation (4.7) with initial data (4.8). Mesh size is $n_x \times n_y = 400 \times 400$, while $\Delta t$ is chosen such that CFL= 0.2. The final time is $T_{end} = 2\pi$.</figcaption></figure>

where $\rho$ is a constant and measures the total initial mass, $\xi$ the cell-size and $f$ the size density repartition.

Despite its simplistic appearance, this model (4.9) is very intriguing when one is interested on the time asymptotic behavior. An interesting discussion is made in [8] and the authors highlight specifically the importance of using an anti-dissipative numerical scheme in order to avoid numerical artifact (diffusion) and then capture the exact asymptotic profile. In order to test our scheme defined in section 3 and show its accuracy, we compare the numerical results with those obtained either by a WENO scheme (see [3]) or by the anti-dissipative scheme (see [8]). Two types of initial distribution functions will be considered: the regular one, which represents the size-density of cells in a normal distribution as in (4.10),

$$ (4.10) \qquad f^0(\xi) = 0.1 \exp(-0.1(\xi - 20)^2). $$
<figure><img src="image_8.png" /><figcaption>FIGURE 7. Polymerization/depolymerization test in homogeneous space: the regular and irregular initial data corresponding to (4.10) and (4.11).</figcaption></figure>

While the irregular one represents that the size-density of cells is concentrated around a spot as in (4.11),

$$ (4.11) \qquad f^0(\xi) = \begin{cases} 0.1, & \text{if } 10 \le \xi \le 30, \\ 0, & \text{else.} \end{cases} $$

These two types of size-density are illustrated in Figure 7.

We first consider the numerical results with the regular data (4.10) (see Figure 8). Our hybrid scheme has the same results as the classical WENO scheme. More precisely, a very smooth profile is formed, which moves towards cells of large size. This observation has a good agreement with the physical behavior of transport equation which does not modify the general shape of the initial distribution function. In biological point of view, the observation show the growth of big size cells at the expense of smaller ones. This competition between big size cells and smaller size ones is due to the fact that at any time there exists a critical size $\xi_{crit} = \frac{1}{c(t)^3}$ such that the velocity vanish; then cells of size $\xi > \xi_{crit}$ grow while cells of size $\xi < \xi_{crit}$ shrink. This competition phenomenon is well known under the name “Ostwald ripening” [18, 19, 15, 28]. At contrast, the anti-dissipative method forms a very sharp front, which is caused by the Downwind flux and it can not be significantly improved by mesh refinement. Even in the zone view for small size particles (see the right column of Figure 8), the anti-dissipative method generates “stairs” looking like oscillations.

Then, with the irregular data (4.11), as pictured at the figure 9, we first notice that the irregular data is regularized by the numerical diffusion by using the WENO scheme. Indeed, I. M. Lifshitz and V. V. Slyozov in their original paper [16] conjecture that the asymptotic solution is in the same form despite the initial data (compare with the left column of Figure 8). However, as pointed out in [8], the irregular data does have an influence for long term solution. Our scheme as the anti-dissipative one preserves well the front propagation at long time (no numerical dissipation) what is physically very important when dealing with transport equation and can be interpreted in biology as the fact that very stiff localized population of cells must remain stiff localized if they are subject to only transport process. Moreover, our scheme corrects well the “stairs” looking like oscillations in the numerical solution of the anti-dissipative method (see the right column of Figure 9).

### 4.2.2. Non-space-homogeneous polymerization/depolymerization type model.
Here we are interested on the numerical behavior of the non-homogeneous system (2.1)-(2.2) with the choice
<figure><img src="image_9.png" /><figcaption>FIGURE 8. Polymerization/depolymerization test in homogeneous space: *Plot solutions of the equation (4.9) with the regular initial data (4.10).* Mesh size is $n_x = 800$, $\Delta t = 1/8$, $CFL \approx 0.12$, $\rho = 41$.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIGURE 9. Polymerization/depolymerization test in homogeneous space: *Plot solutions of the equation (4.9) with the irregular initial data (4.11). Mesh size is $n_x = 800$, $\Delta t = 1/8$, CFL $\approx 0.12$, $\rho = 41$.*</figcaption></figure>
$a(\xi) = \xi^{1/3}$, $b = 1$, that we rewrite as follows

$$ (4.12) \qquad \left\{ \begin{array}{l} \displaystyle \frac{\partial}{\partial t} f(t, x, \xi) + \frac{\partial}{\partial \xi} ((\xi^{1/3} c(t, x) - 1) f(t, x, \xi)) = 0, \quad t \ge 0, x \in \Omega, \xi \ge 0, \\[1.5ex] \displaystyle \partial_t \left( c(t, x) + \int_0^\infty \xi f(t, x, \xi) d\xi \right) = \Delta_x c(t, x), \quad t \ge 0, x \in \Omega. \end{array} \right. $$

The diffusion equation is endowed with homogeneous Neumann boundary condition

$$ \partial_{\nu} c = \nabla c \cdot \nu = 0, \text{ on } \partial \Omega. $$

Here we choose $\Omega \subset \mathbb{R}$ means $1D$ in space-variable. The system (4.12) describes the immersion of a population of cells in a culture of micro-organisms (nutrients) that are subjected to diffusion equation while de size density repartition of cells is parametrized by the space position $x$.

The numerical simulations are performed in the slab $x \in [0, 60]$. The size variable is truncated to $\xi \in [0, 100]$. The initial concentration is defined by

$$ (4.13) \qquad c(t=0, x) = 0.5 \mathbb{I}_{x \in [20, 40]} $$

As we did in the previous test, two types of size-density are considered: a regular one

$$ (4.14) \qquad f(t=0, x, \xi) = 0.01 \exp(-0.2(\xi-30)^2) \mathbb{I}_{x \in [20, 40]} $$

and an irregular one

$$ (4.15) \qquad f(t=0, x, \xi) = 0.01 \mathbb{I}_{x \in [20, 40]} \times \mathbb{I}_{\xi \in [30, 35]} $$

The initial data is presented in Figure 10.

Let us briefly describe the algorithm for the system (4.12). The third order Runge-Kutta method is used as the time discretization for the advection equation of the size-density of macro-particles. In the reconstruction of numerical flux, both the WENO scheme and our hybrid scheme are applied for the purpose of comparison. Then for the diffusion equation of monomers, the crank-Nicolson method is used for time discretization. A classical second-order finite difference discretization method is applied for the lapacian operator, and the integral is approximated by the classical Simpson's rule.

Firstly, in the case with the regular data (4.14), there are almost no difference from the size-density distribution between the two numerical methods (see Figure 11). The initial distribution is in a plaque form at $t = 0$, then it generates a moon shape caused by the diffusion of the concentration of monomers. In evolution in time, we see that the solution is very smooth in size direction.

Secondly, we consider the irregular data (4.15). The micro-organisms (nutrients) concentration $c(t, x)$ and the mass of the cells $\int_0^\infty \xi f(t, x, \xi) d\xi$ are presented in Figure 12. From these two quantities, we see there are no difference at different times. Moreover, the total mass for both methods preserves well in whole time evolution (see Figure 13). However, we observe clearly that numerical dissipation appears with the WENO scheme, while our hybrid method preserves a sharp profile (see Figure 14).

In general, having a scheme able to reconstruct accurately the profile of the solution is of paramount importance in population dynamics; indeed this profile is often used for parameter estimation and comparison with experimental measurements for prediction purposes such as prediction of the evolution of bacterial populations followed by genetic algorithms or the assessment of the average time for the balance of a product (micro-organisms or monomers in our case) in the considered mixture. That's the case in this paper because following the initial distribution, the expected solution does not have the same profile and thus does not lead to same predictions.
<figure><img src="image_11.png" /><figcaption>(a) Regular initial size distribution function</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 10. Polymerization/depolymerization test in non-homogeneous space: the initial data corresponding to (4.13)-(4.15).</figcaption></figure>

## 5. CONCLUSION AND PERSPECTIVE

In this paper we have proposed a hybrid finite volumes scheme based on the flux convex combination between the Anti-Dissipative Method (ADM) [6, 8] and the fifth order WENO method [13]. The obtained numerical results show a good accuracy in reconstructing the solution of transport type equation even in case of discontinuous initial data. Indeed the simulations in figure 8 show results that are as good as the ones in WENO scheme and better than ADM scheme which fail for regular initial distribution. In reverse, for irregular initial distribution, the hybrid scheme show better numerical results than de ADM scheme in the sense that it advects very well the solution without “stairs” and it shows also better results than WENO scheme which develop numerical diffusion for irregular distribution as depicted in figure 9. So, when the WENO order 5 scheme fails because of numerical diffusion artifact, the hybrid method remains anti-dissipative and when ADM scheme develops “stairs” like oscillations, the hybrid scheme corrects them. This property is very suitable for long term asymptotic behavior of the solution of population dynamics, as presented in the numerical simulations for the polymerization/depolymerization type models.
<figure><img src="image_13.png" /><figcaption>(a) WENO Scheme</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b) Hybrid Method</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(c) $x = 30, t = 20$</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(d) WENO Scheme</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(e) Hybrid Method</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(f) $x = 30, t = 80$</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(g) WENO Scheme</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(h) Hybrid Method</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(i) $x = 30, t = 200$</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(j) WENO Scheme</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(k) Hybrid Method</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(l) $x = 30, t = 350$</figcaption></figure>

FIGURE 11. Polymerization/depolymerization test in non-homogeneous space: Plot the size distribution function of the equation (4.12) with the regular initial data (4.13), (4.14) at different time. Mesh size is $n_x \times n_\xi = 100 \times 800$, $\Delta t = 0.1$, CFL $\approx 0.13$.
<figure><img src="image_25.png" /><figcaption>FIGURE 12. Polymerization/depolymerization test in non-homogeneous space: Evolution of the monomers concentration $c(t,x)$ (first row) and the mass of the marco-particles $\int_{0}^{\infty} \xi f(t,x,\xi) d\xi$ (second row) corresponding to the equation (4.12) with the irregular initial data (4.13), (4.15). Mesh size is $n_x \times n_\xi = 100 \times 800$, $\Delta t = 0.1$, CFL $\approx 0.13$.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>FIGURE 13. Polymerization/depolymerization test in homogeneous space: Mass conservation property corresponding to (2.3). Mesh size is $n_x \times n_\xi = 100 \times 800$, $\Delta t = 0.1$, CFL $\approx 0.13$.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>(a) WENO Scheme</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>(b) Hybrid Method</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>(c) $x = 30, t = 20$</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>(d) WENO Scheme</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>(e) Hybrid Method</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>(f) $x = 30, t = 80$</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>(g) WENO Scheme</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>(h) Hybrid Method</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>(i) $x = 30, t = 200$</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>(j) WENO Scheme</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>(k) Hybrid Method</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>(l) $x = 30, t = 350$</figcaption></figure>

<caption>FIGURE 14. Polymerization/depolymerization test in non-homogeneous space: *Plot the size distribution function of the equation (4.12) with the irregular initial data (4.13), (4.15) at different times. Mesh size is $n_x \times n_\xi = 100 \times 800$, $\Delta t = 0.1$, CFL $\approx 0.13$.*</caption>
ACKNOWLEDGMENT

Chang Yang is supported by National Natural Science Foundation of China (Grant No. 11401138) and Heilongjiang Postdoctoral Scientific Research Development Fund (No. LBH-Q14080).

Léon M. Tine is supported by the “BQR” funding set up by Lyon 1 University for the new
associate professors (DR\CG\CD\GG/2014-096).

REFERENCES

[1] D. BALSARA, C.-W. SHU, Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy, *Journal of Computational Physics*, 160 (2000), pp. 405–452.

[2] F. BOUCHUT, An anti-diffusive entropy scheme for monotone scalar conservation law, *Journal of Scientific Computing*, 21 (2004), pp. 1–30.

[3] J. A. CARRILLO AND T. GOUDON, A Numerical Study on Large-Time Asymptotics of the Lifshitz-Slyozov System, *Journal of Scientific Computing*, 20(1) (2004), pp. 69–113.

[4] J. M. CUSHING, An Introduction to Structured Population Dynamics. *Society for Industrial & Applied Mathematics, U.S.*, (1987).

[5] H. M DARIAN, V. ESFAHANIAN, K. HEJRANFAR, A shock-detecting sensor for filtering of high-order compact finite difference schemes, *Journal of computational physics*, 230 (2011), pp. 494–514.

[6] B. DESPRÉS, F. LAGOUTIÈRE, Un schéma non linéaire anti-dissipatif pour l'équation d'advection linéaire. (French) [A nonlinear anti-diffusive scheme for the linear advection equation], *C. R. Acad. Sci. Paris Sér. I Math.*, 328 (1999), no. 10, pp. 939–944.

[7] B. DESPRÉS, F. LAGOUTIÈRE, Contact discontinuity capturing schemes for linear advection, compressible gas dynamics, *Journal of Scientific Computing*, 16 (2001), pp.479–524.

[8] T. GOUDON, F. LAGOUTIÈRE AND L. TINE, Simulations of the Lifshitz-Slyozov equations: the role of coagulation terms in the asymptotic behavior, *Mathematical Methods in the Applied Sciences*, 23(7) (2013), pp. 1177–1215.

[9] T. GOUDON, F. LAGOUTIÈRE AND L. TINE, The Lifshitz-Slyozov equations with space-diffusion of monomers, *Kinetic and Related Models*, 5 (2012), no.2, pp. 325–355.

[10] M. L. GREER, L. PUJO-MENJOUET, GL. WEBB, A mathematical analysis of the dynamics of prion proliferation, *J. Theor. Biol.*, 242 (2006), pp. 598–606.

[11] A. HARTEN, ENO schemes with subcell resolution, *Journal of Computational Physics*, 83 (1989) pp. 148-184.

[12] A. HARTEN, B. ENGQUIST, S. OSHER AND S. CHAKRAVATHY, Uniformely high order accurate essentially non-oscillatory schemes: III, *Journal of Computational Physics*, 71 (1987), pp. 231–303.

[13] G.-S. JIANG AND C.-W. SHU, Efficient implementation of weighted ENO schemes, *Journal of computational physics*, 126 (1996), pp. 202–228.

[14] R.J. LEVEQUE, Finite volume methods for hyperbolic problems, Cambridge university press, 2002.

[15] J. LI, C. GUO, Y. MA, Z. WANG AND J. WANG, Effect of initial particle size distribution on the dynamics of transient Ostwald ripening: A phase field study. *Acta Materialia*, 90 (2015), pp. 10–26.

[16] I. M. LIFSHITZ, V. V. SLYOZOV, The kinetics of precipitation from supersaturated solid solutions. *J. Phys. Chem. Solids*, 19 (1961), pp. 35–50.

[17] J. A. J. METZ AND O. DIEKMANN, The dynamics of physiologically structured populations. *L.N. in biomathematics 68*, Springer., (1986).

[18] B. NIETHAMMER, R. L. PEGO, On the initial-value problem in Lifshitz-Slyozov-Wagner theory of Ostwald ripening. *SIAM J. Math. Anal.*, 31(3) (2000), pp. 467–485.

[19] B. NIETHAMMER, The LSW model for Ostwald ripening with kinetic undercooling. *Proceedings of the Royal Society of Edinburgh, Section: A Mathematics*, 130 (2000), pp. 1337–1361.

[20] B. PERTHAME, Transport equations in biology. Frontiers in mathematics. *Birkhauser Verlag, Basel*, (2007).

[21] P. L. ROE, SOME CONTRIBUTION TO THE MODELLING OF DISCONTINUOUS FLOWS, *Lectures in Applied Mathematics*, 22 (1985), pp. 163–193.

[22] A. M. DE ROOS, A gentle introduction to physiologically structured population models. In: Structured-population models in marine terrestrial, and freshwater Systems. *Chapman & Hall, Population and Community Biology Series*, Vol. 18, (1997).

[23] C.-W. SHU, High order finite difference and finite volume WENO schemes and discontinuous Galerkin methods for CFD, *International Journal of Computational Fluid Dynamics*, 17 (2003), pp. 107–118.

[24] C.-W. SHU AND S. OSHER, Efficient implementation of essentially non-oscillatory shock-capturing schemes, *Journal of Computational Physics*, 77 (1988), pp. 439–471.
[25] C.-W. SHU AND S. OSHER, Efficient implementation of essentially non-oscillatory shock capturing schemes II, *Journal of Computational Physics*, 83(1989), pp. 32–78.

[26] K. SIDDIQI, B. KIMIA, C.-W. SHU, Geometric shock-capturing ENO schemes for subpixel interpolation computation and curve evolution, *Graphical Models and Image Processing* (CVGIP:GMIP), 59 (1997) pp. 278–301.

[27] P. K. SWEBY, HIGH RESOLUTION SCHEMES USING FLUX LIMITERS FOR HYPERBOLIC CONSERVATION LAWS, *SIAM Journal on Numerical Analysis*, 21 (1984), pp. 995–1011.

[28] J. SVOBODA AND F.D. FISCHER, Generalization of the Lifshitz-Slyozov-Wagner coarsening theory to non-dilute multi-component systems. *Acta Materialia*, 79 (2014), pp. 304–314.

[29] Z. XU AND C.-W. SHU, Anti-diffusive flux corrections for high order finite difference WENO schemes, *Journal of Computational Physics*, 205(2005), pp. 458–485.

[30] H. YANG, An artificial compression method for ENO schemes: the slope modification method, *Journal of Computational Physics*, 89 (1990), pp. 125–160.

DEPARTMENT OF MATHEMATICS, HARBIN INSTITUTE OF TECHNOLOGY, 92 WEST DAZHI STREET, NAN GANG DISTRICT, HARBIN 150001, CHINA

*E-mail address*, C. Yang: *yangchang@hit.edu.cn*

UNIVERSITÉ DE LYON, UNIVERSITÉ LYON 1, CNRS UMR 5208, INSTITUT CAMILLE JORDAN 43 BLVD DU 11 NOVEMBRE 1918, F-69622 VILLEURBANNE-CEDEX, FRANCE. & INRIA TEAM DRACULA, INRIA CENTER GRENOBLE RHONE-ALPES, FRANCE.

*E-mail address*, L.M. Tine: *leon-matar.tine@univ-lyon1.fr*