Casimir Emako
email: emako@ann.jussieu.fr

Min Tang
email: tangmin@sjtu.edu.cn

WELL-BALANCED AND ASYMPTOTIC PRESERVING SCHEMES FOR KINETIC MODELS

Keywords: Mathematics Subject Classification. 65M08, 35Q20, 35Q92 Kinetic models, well-balanced, asymptotic preserving

In this paper, we propose a general framework for designing numerical schemes that have both well-balanced (WB) and asymptotic preserving (AP) properties, for various kinds of kinetic models. We are interested in two different parameter regimes, 1) When the ratio between the mean free path and the characteristic macroscopic length ε tends to zero, the density can be described by (advection) diffusion type (linear or nonlinear) macroscopic models; 2) When ε = O(1), the models behave like hyperbolic equations with source terms and we are interested in their steady states. We apply the framework to three different kinetic models: neutron transport equation and its diffusion limit, the transport equation for chemotaxis and its Keller-Segel limit, and grey radiative transfer equation and its nonlinear diffusion limit. Numerical examples are given to demonstrate the properties of the schemes.

Introduction

Transport equations are important since they arise in many important applications, ranging from neutron transport, radiative transfer, semiconductor device simulation to E.coli chemotaxis. These models describe particles that travel freely for a certain distance and then change their directions by scattering, interacting with background media or tumbling. The mean free path (the average distance a particle travels between two successive velocity changes) is an important parameter, when it is small compared to the typical length scales, various macroscopic model can be derived asymptotically or analytically [START_REF] Chalub | Kinetic models for chemotaxis and their drift-diffusion limits[END_REF][START_REF] Edward W Larsen | Asymptotic solutions of numerical transport problems in optically thick, diffusive regimes[END_REF][START_REF] Othmer | The diffusion limit of transport equations. II. Chemotaxis equations[END_REF].

Let the dimensionless parameter ε denote the ratio of the mean free path and the typical length scale. Numerical solutions to the transport equation are challenging when ε is small, since it requires the numerical resolution of the small scale. To develop a multi-scale scheme whose stability and convergence are independent of ε refers to the asymptotic preserving (AP) property. When ε = O(1), the models behave like hyperbolic equations with source terms. When the source terms in the system become stiff, the usual numerical methods may give poor approximations to the steady state solutions [START_REF] Botchorishvili | Equilibrium schemes for scalar conservation laws with stiff sources[END_REF][START_REF] Leveque | Balancing source terms and flux gradients in high-resolution godunov methods: The quasi-steady wave-propagation algorithm[END_REF]. To maintain the steady states or to achieve them in the long time limit with an acceptable level of accuracy refers to the well balanced (WB) property.

In this work, we present a general framework to build schemes that have both properties. Two different kinetic models are considered to illustrate the idea of our proposed framework Transport equation for Chemotaxis and its Keller-Segel limit Bacteria undergo run and tumble process as mentioned in [START_REF] Berg | coli in Motion[END_REF][START_REF] Saragosti | Mathematical description of bacterial traveling pulses[END_REF][START_REF] Saragosti | Directional persistence of chemotactic bacteria in a traveling concentration wave[END_REF]. During the run phase, bacteria move along a straight line and change their directions during the tumble phase. This individual motion is called the velocity jump process and modeled by the so-called Othmer-Dunbar-Alt model [START_REF] Alt | Biased random walk models for chemotaxis and related diffusion approximations[END_REF][START_REF] Othmer | Models of dispersal in biological systems[END_REF] which reads (1.1)

       ∂ t f ε + v ε • ∇ x f ε = 1 ε 2 1 |V | V 1 + εφ(v ′ • ∇ x S) f ε (v ′ )dv ′ -1 + εφ(v • ∇ x S) f ε (v) , ∂ t S ε -D∆S ε + αS ε = βρ ε (x, t) := 1 |V | V f ε dv.
Here f ε (x, v, t) is the amount of cells with velocity v ∈ V (V is a sphere) at position x and time t > 0, φ(u) is a decreasing function in u, S ε (x, t) is the concentration of the chemical substance. The parameters D, α, β are positive constants and ε is the ratio of the average run distance between two successive tumbles and the characteristic macroscopic length. When ǫ → 0, the solution of (1.1) f ε (x, v, t) tends to ρ 0 (x, t), which is independent of v and satisfies the following Keller-Segel type equation [START_REF] Chalub | Kinetic models for chemotaxis and their drift-diffusion limits[END_REF][START_REF] Emako | Existence and diffusive limit of a two-species kinetic model of chemotaxis[END_REF][START_REF] Hyung | Drift-diffusion limits of kinetic models for chemotaxis: a generalization[END_REF][START_REF] Keller | Model for chemotaxis[END_REF][START_REF] Othmer | The diffusion limit of transport equations. II. Chemotaxis equations[END_REF] :

(1.2)    ∂ t ρ = 1 3 ∆ρ + ∇ 1 |V | V vφ(v∂ x S)dv ρ , ∂ t S -D∆S + αS = βρ.
Many numerical schemes have been proposed to study (1.1) [START_REF] José | An asymptotic preserving scheme for the diffusive limit of kinetic systems for chemotaxis[END_REF][START_REF] Filbet | Numerical simulations of kinetic models for chemotaxis[END_REF][START_REF] Gosse | Asymptotic-preserving and well-balanced schemes for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF].

Gray radiative transfer equation and its nonlinear diffusion limit: The gray radiative transfer equation concerns photon transport and its interaction with the background material. It has wide application in astrophysics and inertial confinement fusion. The system for the radiative intensity I and the material temperature is

(1.3)      1 c ∂ t I + 1 ε v • ∇ x I = σ ε 2 1 |V | acT 4 -I + q(v), C v ∂ t T = σ ε 2 (ρ -acT 4 ),
where ρ := V Idv, σ(x, T ) is the opacity and a, c, C v are positive constants represent the radiation constant, light speed and heat capacity respectively. Similar to the chemotaxis case, when ε goes to zero, the radiative intensity I approaches a Planckian at the local temperature 1 |V | acT 4 and the photon temperature satisfies the following nonlinear diffusion equation:

(1.4) a∂ t T 4 + C v ∂ t T = ∇ x ( ac 3σ ∇ x T 4 ) + V q dv.
It is important to preserve the steady state temperature distribution when t → ∞.

A scheme for such problems is AP if it preserves the limiting equation (1.2) or (1.4) when ε → 0 at the discrete level. AP schemes were first designed for the neutron transport equation to use unresolved cells to capture the macroscopic diffusion limit model [START_REF] Larsen | Asymptotic solutions of numerical transport problems in optically thick,diffusive regimes ii[END_REF][START_REF] Edward W Larsen | Asymptotic solutions of numerical transport problems in optically thick, diffusive regimes[END_REF]. It has been successfully extended to a lot of applications, we refer to the review paper [START_REF] Shi | Asymptotic preserving (ap) schemes for multiscale kinetic and hyperbolic equations: a review[END_REF] for more discussions. In the literature, many papers have been devoted to build AP schemes for diffusion limit of the transport equations, for example [START_REF] Adams | Discontinuous finite element transport solutions in thick diffusive problems[END_REF][START_REF] José | An asymptotic preserving scheme for the diffusive limit of kinetic systems for chemotaxis[END_REF][START_REF] Di Russo | Global existence of smooth solutions to a two-dimensional hyperbolic model of chemotaxis[END_REF][START_REF] Mieussens | On the asymptotic preserving property of the unified gas kinetic scheme for the diffusion limit of linear kinetic models[END_REF].

WB schemes are developed for hyperbolic equations with source terms. They have been proposed for various applications, including the Saint-Venant system with a source term, the non-isothermal nozzle flow equations, etc. By balancing the numerical flux with the source term, WB schemes can capture the steady state solutions exactly or with at least a second order accuracy [START_REF] Shi | An efficient method for computing hyperbolic systems with geometrical source terms having concentrations[END_REF]. A very limited list of such references includes but not limits to quasi-steady scheme [START_REF] Leveque | Balancing source terms and flux gradients in high-resolution godunov methods: The quasi-steady wave-propagation algorithm[END_REF], kinetic schemes [START_REF] Perthame | A kinetic scheme for the saint-venant system with a source term[END_REF][START_REF] Xu | A well-balanced gas-kinetic scheme for the shallow-water equations with source terms[END_REF], central schemes [START_REF] Kurganov | Central-upwind schemes for the saint-venant system[END_REF], interface scheme [START_REF] Shi | An efficient method for computing hyperbolic systems with geometrical source terms having concentrations[END_REF], and other references there in [START_REF] Georgij Bispen | Imex large time step finite volume methods for low froude number shallow water flows[END_REF].

Recently designing schemes that are either AP or WB, or have both properties attracts a lot of interests, both for kinetic chemotaxi model [START_REF] Gosse | Asymptotic-preserving and well-balanced schemes for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF][START_REF] Gosse | A well-balanced scheme for kinetic models of chemotaxis derived from onedimensional local forward-backward problems[END_REF][START_REF] Natalini | A well-balanced numerical scheme for a one dimensional quasilinear hyperbolic model of chemotaxis[END_REF] and grey radiative transfer equation [START_REF] Gosse | Transient radiative transfer in the grey case: Well-balanced and asymptoticpreserving schemes built on case's elementary solutions[END_REF][START_REF] Sun | Asymptotic preserving unified gas kinetic scheme for grey radiative transfer equations[END_REF]. Generally, it is very hard to have a scheme that have both properties and this is the goal of our present paper.

The paper is organized as follows: In section 2, the scheme framework composed of two steps: prediction step and steady problem step, is described and we show that the AP and WB properties can be achieved. Section 3 and 4 are respectively devoted to the construction of AP and WB schemes for the transport equation for chemotaxis and the gray radiative transfer equation. In the prediction step, we first extend the unified gas approach for neutron transport equation [START_REF] Mieussens | On the asymptotic preserving property of the unified gas kinetic scheme for the diffusion limit of linear kinetic models[END_REF] to the chemotaxis model (1.1) and construct a new AP scheme for E.coli chemotaxis, while, for the gray radiative transfer equation, the scheme in [START_REF] Sun | Asymptotic preserving unified gas kinetic scheme for grey radiative transfer equations[END_REF] is employed. Then, in the steady problem step, we use the numerical results obtained by an AP scheme for the steady state equation of (1.1) or (1.3) to modify the numerical flux. The performances of the proposed schemes are presented in section 5 and we conclude with some discussion in section 6.

The scheme framework and its WB and AP properties

In this part, we introduce a general framework of designing WB and AP schemes, while the details of the discretization are given in section 3 and section 4. In the subsequent part of the paper, we consider the one dimensional problem and use a uniform grid with

x i = i∆x, i ∈ Z, t n = n∆t, n ∈ N.
Extensions to the two dimensional case are straightforward.

To illustrate the idea, we consider the following simplest one dimensional neutron transport equation (2.1)

∂ t f + 1 ε v∂ x f = σ T ε 2 (ρ -f ) -σ a ρ + q, with ρ := 1 2 1
-1 f dv. When ε → 0, the solution of the above equation tends to ρ 0 which satisfies the following diffusion equation

∂ t ρ 0 -∂ x 1 3σ T ∂ x ρ 0 + σ a ρ 0 = q.
The formal derivation is standard by substituting the Chapman-Enskog expansion

f = f (0) + εf (1) + ε 2 f (2) + • • • into (2.1
). When ε is small, the solution to (2.1) f can be approximated by

(2.2) f = ρ (0) - ε σ T v∂ x ρ (0) + O(ε 2 ).
We are interested in ε ranging from O(1) to very small, the numerical scheme writes:

(2.3)

f n+1 i -f n i ∆t + v ε∆x (1 -α)f n i + αf n+1 i -f n i-1/2 = 0, v > 0, f n+1 i -f n i ∆t + v ε∆x f n i+1/2 -(1 -α)f n i + αf n+1 i = 0, v < 0,
where f n i are the approximation of the average f (x, t n ) in the interval

(x i-1/2 , x i+1/2 ) = (x i -∆x/2, x i + ∆x/2) and (2.4) α = min 1, ∆t ε .
The determination of f n i+1/2 consists of two steps: the prediction step and the steady problem step.

Prediction step: The prediction step is to evolve the equation (2.1) for one time step by an AP scheme. Starting from f n i obtained from the nth iteration, the predictions f n+1 i and ρn+1

i = 1 2 1 -1 f n+1 i
dv can be found by any scheme that is AP. The requirement that, when ε → 0, the scheme preserves the diffusion limit at the discrete level can only be achieved when

(2.5) f n+1 i ≈ ρn+1 i + O(ε)
and ρn i is a discretization of the limiting diffusion equation. Here in this paper we use the unified gas approach developed in [START_REF] Mieussens | On the asymptotic preserving property of the unified gas kinetic scheme for the diffusion limit of linear kinetic models[END_REF][START_REF] Xu | A gas-kinetic {BGK} scheme for the Navier-Stokes equations and its connection with artificial dissipation and godunov method[END_REF][START_REF] Sun | Asymptotic preserving unified gas kinetic scheme for grey radiative transfer equations[END_REF] for various kinetic models.

Steady problem step:

The second step is devoted to the computation of the steady state problem. On each cell [x i , x i+1 ], we solve the following stationary problem:

(2.6) v • ∂ x f = σ T ε ρ -f -εσ a ρ + εq together with inflow boundary conditions (2.7) f (x i , v) = (1 -α)f n i + α f n+1 i , v > 0, f (x i+1 , v) = (1 -α)f n i+1 + α f n+1 i+1 , v < 0.
Here f n+1 i is determined by the prediction step. Then f n i+1/2 in (2.3) is given by the outflow of the above steady state problem in each cell such that

(2.8) f n i+1/2 (v) = f (x i+1 , v), v > 0 f n i+1/2 (v) = f (x i , v), v < 0.
First of all, when ǫ is at O(1) or ∆x, ∆t ≪ ε, i.e. α = ∆t/ε, we show that (2.3) is a consist discretization for (2.1). We only consider the case when v > 0. Since f n i-1/2 is obtained from the steady problem step, it can be approximated by (1

-∆t/ε)f n i-1 + ∆t/ε f n+1 i-1 + ∆x∂ x f n (x i-1 , v) =(1 -∆t/ε)f n i-1 + ∆t/ε f n+1 i-1 + ∆x v σ T ε ρn i-1 -f n i-1 -εσ a ρn i-1 + εq =f n i-1 + ∆x v σ T ε ρ n i-1 -f n i-1 -εσ a ρ n i-1 + εq + O(∆t 2 ) + ∆x v O(∆x).
From the CFL condition, ∆t < ǫ max |v| ∆x, (2.3) can then be approximated by

f n+1 i -f n i ∆t + v ε∆x f n i -f n i-1 = σ T ε 2 ρ n i-1 -f n i-1 -σ a ρ n i-1 + q + O(∆x).
The derivation for v < 0 is the same. Moreover, we have the following proposition: Proposition 2.1. When AP schemes are used in the prediction step and steady problem, the numerical scheme given in (2.3) has both WB and AP properties.

Proof. AP property: When ε tends to zero, since AP schemes are used for the prediction step, we have (2.5). From (2.4), when ε is small, α = 1, then the boundary conditions for the steady state problem (2.6)

on each cell becomes f (x i , v) = ρn+1 i + O(ε), v > 0, f (x i+1 , v) = ρn+1 i+1 + O(ε), v < 0.
When ε is small, at the leading order, the boundary conditions for the steady state problem (2.6) are isotropic. Therefore we have f = ρ + O(ε)

and f n i+1/2 = f (x i+1 , v) = ρn+1 i+1 + O(ε), v > 0, f n i+1/2 = f (x i , v) = ρn+1 i + O(ε), v < 0.
Then from (2.3),

f n+1 i = fi-1 2 + ε∆x v∆t f n i 1 + ε∆x v∆t = f n i-1 2 + O(ε) = ρn+1 i + O(ε), v > v 0 , f n+1 i = fi+ 1 2 -ε∆x v∆t f n i 1 -ε∆x v∆t = f n i+ 1 2 + O(ε) = ρn+1 i + O(ε), v < -v 0 .
with v 0 small positive value away from 0. Therefore, we get

ρ n+1 i = ρn+1 i + O(ε).

Since ρn+1

i satisfy the macroscopic equation at the discrete level, i.e. the leading order of f n+1 i evolves according to a discrete diffusion equation, which gives the AP property of the proposed scheme (2.3).

WB property:

To prove the WB property, we assume that at time t n , f n i = f (x i , v) where f (x, v) is the solution to the steady state equation (2.6) with given inflow boundary conditions on the whole computational domain. We need to show that when ∆x, ∆t are small enough,

f n+1 i = f n i + O(∆x∆t).
For the resolved case under consideration, α = ∆t ε .

Since the scheme at the prediction step is consistent, consider the simplest first order method, we have

f n+1 i = f n i + O(∆x, ∆t). The inflow boundary conditions in (2.7) becomes f (x i , v) = (1 -α)f n i + α f n+1 i = f n i + O(∆x∆t), v > 0, f (x i+1 , v) = (1 -α)f n i+1 + α f n+1 i+1 = f n i+1 + O(∆x∆t), v < 0.
Since the steady state problem is solved in a cell of size ∆x, by the simplest Taylor expansion, the out flow can be approximated by at least second order accuracy, so that

f n i+1/2 = f (x i+1 , v) + O(∆x∆t, ∆x 2 ), v > 0; f n i+1/2 = f (x i , v) + O(∆x∆t, ∆x 2 ), v < 0.
Thus, (2.3) yields (2.9)

f n+1 i = fi-1 2 + ε∆x v∆t f n i 1 + ε∆x v∆t = f (x i , v) + O(∆x∆t, ∆x 2 ), v > 0, f n+1 i = fi+ 1 2 -ε∆x v∆t f n i 1 -ε∆x v∆t = f (x i , v) + O(∆x∆t, ∆x 2 ), v < 0.
Due to the CFL condition, O(∆t) = O(∆x), the solution approximate the steady state f with a formally second order accuracy O(∆x 2 ).

Remark 2.2. The framework in this section does not depend on specific AP methods for the time evolutionary problem and the steady state problem, we can choose any scheme in the literature that has been proved to be AP for various kinetic models.

Remark 2.3. The accuracy of preserving the steady state problem can be improved by using higher order method in space and time, or we can repeat once the steady problem step by replacing f n+1 i from the prediction step by f n+1 i obtained in (2.9).

The chemotaxis kinetic model

In this section, we apply the framework in section 2 to the chemotaxis kinetic model (1.1). Two specific AP schemes are respectively chosen for the time evolutionary problem and steady state problem. We first extend the unified gas kinetic scheme (UGKS) in [START_REF] Mieussens | On the asymptotic preserving property of the unified gas kinetic scheme for the diffusion limit of linear kinetic models[END_REF][START_REF] Xu | A unified gas-kinetic scheme for continuum and rarefied flows[END_REF][START_REF] Xu | A gas-kinetic {BGK} scheme for the Navier-Stokes equations and its connection with artificial dissipation and godunov method[END_REF] to get an AP scheme for the time evolutionary problem, then solve the steady state equation by extending the scheme in [START_REF] Shi | A uniformly second order numerical method for the one-dimensional discrete-ordinate transport equation and its diffusion limit with interface[END_REF] which was originally designed for the isotropic neutron transport equation. In [START_REF] José | An asymptotic preserving scheme for the diffusive limit of kinetic systems for chemotaxis[END_REF], the authors have proposed AP schemes for the chemotaxis kinetic model, our approach is different. The details are described below.

As in [START_REF] Mieussens | On the asymptotic preserving property of the unified gas kinetic scheme for the diffusion limit of linear kinetic models[END_REF], the UGKS is a finite volume approach of integrating (1.1) over

[x i-1/2 , x i+1/2 ] × [t n , t n+1 ] × V . Let f n i = 1 ∆x x i+1/2 x i-1/2 f (x, v, t n ) dx, ρ n i = 1 |V | V f n i dv.
ρ n+1 i -ρ n i ∆t + F n i+1/2 -F n i-1/2 ∆x = 0, (3.2) f n+1 i -f n i ∆t + Φ n i+1/2 -Φ n i-1/2 ∆x = 1 ε 2 ρ n+1 i -f n+1 i + 1 ε 1 |V | V φ(vσ i+1/2 )f n i -φ(vσ i+1/2 )f n i .
Here the numerical flux

(3.3) Φ n i+1/2 ≈ 1 ε∆t t n+1 t n vf (x i+1/2 , v, t) dt, F n i+1/2 ≈ 1 |V | V 1 ε∆t t n+1 t n vf (x i+1/2 , v, t) dt dv.
Since one dimensional problem is considered here, we use discrete ordinate method for the velocity discretization. The most crucial step for UGKS is to determine Φ n i+1/2 and F n i+1/2 . The details are as below:

3.1. Determine Φ n i+1/2 , F n i+1/2 .
In one dimensional space, the chemotaxis model (1.1) can be rewritten as follows:

(3.4) ∂ t f ε + 1 + εφ(v∂ x S ε ) ε 2 f ε + v ε ∂ x f ε = 1 ε 2 T 1 f, where (T 1 f )(x, t) := 1 |V | V 1 + εφ(v∂ x S) f (x, v, t)dv.
To approximate the flux F n i+1/2 in (3.3), we first approximate ∂ x S by a piecewise constant function such that

∂ x S ≈ ∂ x S(x i+1/2 ) := σ i+1/2 , x ∈ [x i , x i+1 ).
It is important to note that σ i+1/2 approximate ∂ x S in the interval [x i , x i+1 ) while f n i is the density average over the cell [x i-1/2 , x i+1/2 ). This choice is important to preserve the advection term in the limiting Keller-Segel model when ε becomes small. On the interval [x i , x i+1 ), multiplying both sides of (3.4) by exp

(1+εφ(vσ i+1/2 ) ε 2 t yields d dt f (x + v ε t, v, t) exp (1 + εφ(vσ i+1/2 ) ε 2 t = T 1 f (x, t) ε 2 exp (1 + εφ(vσ i+1/2 ) ε 2 t .
Integrating the above equation over (t n , t) yields (3.5)

f (x i+1/2 , v, t) = f (x i+1/2 - v ε (t -t n ), v, t n ) exp - (1 + εφ(vσ i+1/2 ) ε 2 (t -t n ) + 1 ε 2 t t n T 1 f (x i+1/2 - v ε (t -s), s) exp - (1 + εφ(vσ i+1/2 ) ε 2 (t -s) ds.
This is an exact expression for f (x i+1/2 , v, t) that will be used to determine Φ n i+1/2 , F n i+1/2 in (3.3). At this stage, we need to approximate f (x, v, t n ) and T 1 f (x, t) on the right hand side of (3.5). f is approximated by a piecewise constant function and T 1 f by a piecewise linear function as follows:

f (x, v, t n ) = f n i , x < x i+1/2 , f n i+1 , x > x i+1/2 , T 1 f (x, t) = T 1 f n i+1/2 + δ L T 1 f n i+1/2 (x -x i+1/2 ), x < x i+1/2 , T 1 f n i+1/2 + δ R T 1 f n i+1/2 (x -x i+1/2 ), x > x i+1/2 .
Here

T 1 f n i+1/2 , δ L T 1 f n i+1/2 , δ R T 1 f n i+1/2
are defined by:

                   T 1 f n i+1/2 := 1 |V | V - (1 + εφ(vσ i+1/2 ))f n i+1 + 1 |V | V + (1 + εφ(vσ i+1/2 ))f n i , δ L T 1 f n i+1/2 := T 1 f n i+1/2 -T 1 f n i ∆x/2 , δ R T 1 f n i+1/2 := T 1 f n i+1 -T 1 f n i+1/2 ∆x/2 , with V + = V ∩ R + and V -= V ∩ R -.
Substituting the above approximations into (3.5) yields an expression for f (x i+1/2 , v, t) such that: For v > 0, (3.6)

f (x i+1/2 , v, t) = f n i exp - (1 + εφ(vσ i+1/2 ) ε 2 (t -t n ) + T 1 f n i+1/2 1 + εφ(vσ i+1/2 ) 1 -exp - (1 + εφ(vσ i+1/2 ) ε 2 (t -t n ) + vε δ L T 1 f n i+1/2 (1 + εφ(vσ i+1/2 )) 2 • 1 + 1 + εφ(vσ i+1/2 ) ε 2 (t -t n ) exp - (1 + εφ(vσ i+1/2 ) ε 2 (t -t n ) -1 ,
and for v < 0, (3.7)

f (x i+1/2 , v, t) = f n i+1 exp - (1 + εφ(vσ i+1/2 ) ε 2 (t -t n ) + T 1 f n i+1/2 1 + εφ(vσ i+1/2 ) 1 -exp - (1 + εφ(vσ i+1/2 ) ε 2 (t -t n ) + vε δ R T 1 f n i+1/2 (1 + εφ(vσ i+1/2 )) 2 • 1 + 1 + εφ(vσ i+1/2 ) ε 2 (t -t n ) exp - (1 + εφ(vσ i+1/2 ) ε 2 (t -t n ) -1 .
Then the flux Φ n i+1/2 (v) in (3.3) can be approximated by

(3.8) Φ i+1/2 (v) = Avf n i+1 + BvT 1 f n i+1/2 + Cv 2 δ R T 1 f n i+1/2 , for v < 0, Φ i+1/2 (v) = Avf n i + BvT 1 f n i+1/2 + Cv 2 δ L T 1 f n i+1/2 , for v > 0,
where the coefficients A(v, ε, ∆t), B(v, ε, ∆t), C(v, ε, ∆t) can be determined explicitly such that (3.9)

A(v, ε, ∆t) : = ε ∆t 1 + εφ(vσ i+1/2 ) 1 -exp - 1 + εφ(vσ i+1/2 ) ε 2 ∆t , B(v, ε, ∆t) : = 1 ε(1 + εφ(vσ i+1/2 )) - ε ∆t(1 + εφ(vσ i+1/2 )) 2 1 -exp - 1 + εφ(vσ i+1/2 ) ε 2 ∆t , C(v, ε, ∆t) : = 2ε 2 ∆t(1 + εφ(vσ i+1/2 )) 3 1 -exp - 1 + εφ(vσ i+1/2 ) ε 2 ∆t - 1 (1 + εφ(vσ i+1/2 )) 2 1 + exp - 1 + εφ(vσ i+1/2 ) ε 2 ∆t . Furthermore, F n i+1/2 in (3.3) is given by (3.10) F n i+1/2 = 1 |V | V - Avf n i+1 dv + 1 |V | V + Avf n i dv + 1 |V | T 1 f n i+1/2 V vBdv + 1 |V | δ R T 1 f n i+1/2 V - Cv 2 dv + 1 |V | δ L T 1 f n i+1/2 V + Cv 2 dv.
and we complete the construction of the scheme. We can prove that the scheme is AP by employing similar approach as in [START_REF] Mieussens | On the asymptotic preserving property of the unified gas kinetic scheme for the diffusion limit of linear kinetic models[END_REF].

AP property.

In this part, we give a formal derivation of the AP property for the UGKS proposed in (3.1)-(3.2). When ε goes to zero, asymptotic expansions of A, B, C given in (3.9) read

A = O(ε), B = 1 ε -φ(vσ i+1/2 ) + O(ε), C = -1 + O(ε).
The leading order term of (3.2) yields f n+1 i = ρ n+1 i + O(ε) and we only need to show that (3.1) satisfies the equation for ρ in (1.2), at the discrete level. Suppose that

f n i = ρ n i + O(ε), then              T 1 f n i+1/2 = 1 2 ρ n i + ρ n i+1 + O(ε), δ L T 1 f n i+1/2 = ρ n i+1 -ρ n i ∆x + O(ε), δ R T 1 f n i+1/2 = ρ n i+1 -ρ n i ∆x + O(ε).
We deduce that the expansion of F n i+1/2 reads:

F n i+1/2 = - ρ n i + ρ n i+1 2|V | V vφ(vσ i+1/2 )dv - ρ n i+1 -ρ n i 3∆x + O(ε).
Therefore,

F n i+1/2 -F n i-1/2 ∆x = - ρ n i+1 -2ρ n i + ρ n i-1 3(∆x) 2 + - 1 |V | V vφ(vσ i+1/2 )dv ρ n i + ρ n i+1 2 + 1 |V | V vφ(vσ i-1/2 )dv ρ n i + ρ n i-1 2 + O(ε).
In the limit of ε → 0, the discretization (3.1) becomes

ρ n+1 i -ρ n i ∆t = ρ n i+1 -2ρ n i + ρ n i-1 3(∆x) 2 + 1 |V | V vφ(vσ i+1/2 )dv ρ n i + ρ n i+1 2 - 1 |V | V vφ(vσ i-1/2 )dv ρ n i + ρ n i-1 2 .
which is a consistent discretization of the equation for ρ in (1.2). The proposed scheme is AP.

3.3.

Steady state problem for the chemotaxis kinetic model. To solve the steady state problem, we start from the most used discrete ordinate method [START_REF] Chandrasekhar | Radiative transfer[END_REF][START_REF] Kubelka | New contributions to the optics of intensely light-scattering materials[END_REF].

The discrete ordinate method is to discretize the velocity space

V = [-1, 1] by a quadrature set {µ m , ω m } (m ∈ V m = {-N, • • • , -1, 1, • • • , N }).
In one dimensional case, the most used and well accepted quadrature is the Gaussian quadrature. In order to preserve the diffusion limit equation, the points µ m and weights ω m satisfy

µ -m = -µ m , ω -m = ω m , m ∈ 1 • • • N, m ω m = 2, m ω m µ 2 m = 2 3 .
The discrete ordinate method for the steady state chemotaxis kinetic model on each cell [x i , x i+1 ) writes (3.11)

µ m ∂ x f m = 1 2ε n∈V ω n 1 + εφ(µ n σ i+1/2 ) f n - 1 ε 1 + εφ(µ m σ i+1/2 ) f m , m ∈ V m .
This is a linear ODE system with constant matrix coefficient. Together with the inflow boundary conditions, the exact solution can be obtained analytically. We use the following procedure to construct the general solution on each interval [x i , x i+1 ). We seek eigenfunctions of the form

(3.12) l m exp(- ζ ε x).
By substituting this form into (3.11), one obtains

(3.13) 1 + εφ(µ m σ i+1/2 ) -µ m ζ l m = 1 2 n∈V ω n 1 + εφ(µ m σ i+1/2 ) l m .
The above equation holds for all m ∈ V m which indicates that ζ is the eigenvalue of the matrix

(3.14) E = U -1 I - 1 2 W I + εφ(U σ i+1/2 ) .
Here U = Diag(µ m ) m∈Vm and φ(U σ i+1/2 ) = Diag φ(µ m σ i+1/2 ) m∈Vm are diagonal matrixes, W is a rank 1 matrix whose rows are (ω m ) m∈Vm . Due to (3.13), ζ is the root of

(3.15) 1 2 m∈Vm 1 + εφ(µ m σ i+1/2 ) ω m 1 + εφ(µ m σ i+1/2 ) -µ m ζ = 1
and the eigenvector associated to ζ can be given by (3.16)

l m = 1 1 + εφ(µ m σ i+1/2 ) -µ m ζ , m ∈ V m .
It is easy to check that when

ε = 0, ζ = 0 is a simple root if φ(u) = -φ(-u). Besides if ν m = 1+εφ(µmσ i+1/2 ) µm (m ∈ V m ) are different from each other, since
the sign of ν m are determined by µ m the functional 1 2 m∈Vm 1+εφ(µmσ i+1/2 ) ωm 1+εφ(µmσ i+1/2 )-µmζ tends to plus or minus infinity when ζ approaches ν m from left or right. Therefore, (3.15) has 2N simple roots when ε is away from zero. However when ε → 0, E becomes U -1 (I -W/2) whose eigenvalues admit a zero double root [START_REF] Chandrasekhar | Radiative transfer[END_REF][START_REF] Shi | A uniformly second order numerical method for the one-dimensional discrete-ordinate transport equation and its diffusion limit with interface[END_REF]. The general solutions to (3.11) are different for these two different cases, we discuss them separately in the subsequent part.

f m (x) = ζn<0 B n l n m exp - ζ n ε (x -x i+1 ) + ζn>0 A n l n m exp - ζ n ε (x -x i ) ,
The coefficients B n , A n can be determined by the inflow boundary conditions

f m (x i+1 ) = f R m , m < 0, f m (x i ) = f L m , m > 0.
Then the outflows are given by

f m (x i ) = ζn<0 B n l n m exp ζ n ε ∆x + ζn>0 A n l n m , m < 0, f m (x i+1 ) = ζn<0 B n l n m + ζn>0 A n l n m exp - ζ n ε ∆x , m > 0.
Case II: We consider the limiting case when E = U -1 (I -W/2), the eigenvalues for this matrix have been proved in [START_REF] Shi | A uniformly second order numerical method for the one-dimensional discrete-ordinate transport equation and its diffusion limit with interface[END_REF] to have the following propperty Theorem 3.1. The equation

1 2 m∈Vm ω m 1 -µ m ζ = 1
has 2(N -1) simple roots appear in positive/negative paris while 0 is a double root.

We set the double zero eigenvalue ζ N and ζ -N and rearrange eigenvalues ζ n from the lowest to the highest.

ζ 1-N < • • • < ζ -1 < 0 < ζ 1 < • • • < ζ N -1 .
For eigenvalues different of zero, eigenfunctions are exponential functions as in (3.12). Let us find eigenfunctions corresponding to the zero eigenvalue. By using the same token as in [START_REF] Shi | A uniformly second order numerical method for the one-dimensional discrete-ordinate transport equation and its diffusion limit with interface[END_REF], we are looking for solution of the form α m x + β m .

Injecting α m x + β m into (3.11) gives two solutions

1 1 + εφ(µ m σ i+1/2 ) , x 1 + εφ(µ m σ i+1/2 ) - εµ m 1 + εφ(µ m σ i+1/2 ) 2 .
Finally the general solution of (3.11) writes

f m (x) = c N 1 + εφ(µ m σ i+1/2 ) + c -N 1 + εφ(µ m σ i+1/2 ) x - εµ m 1 + εφ(µ m σ i+1/2 ) + 1≤|n|≤N -1 c n l n m exp - ζ n ε x ,
where c n (1 ≤ |n| ≤ N ) are constants to be determined by the boundary conditions. Since ε can become very small, overflow may occur when evaluating exp -ξ n x/ε . In order to have c n at the same scale as f m (x), we rewrite the general solution in the interval (x i , x i+1 ) into (3.17)

f m (x) = c 1 1 + εφ(µ m σ i+1/2 ) + c 2 1 + εφ(µ m σ i+1/2 ) (x -x i ) - µ m ε 1 + εφ(µ m σ i+1/2 ) + n<0 B n l n m exp - ζ n ε (x -x i+1 ) + n>0 A n l n m exp - ζ n ε (x -x i ) ,
with c 1 , c 2 , B n , A n constants. Then we can determine c 1 , c 2 , B n , A n by the inflow boundary condition, so that the outflow.

The radiative transport equation

For the radiative transport equation (1.3), the same as in the previous section, we can choose two specific AP schemes respectively for the time evolutionary problem and steady state problem. For the time evolutionary problem, the UGKS scheme proposed in [START_REF] Sun | Asymptotic preserving unified gas kinetic scheme for grey radiative transfer equations[END_REF] is employed, while similar discretization as in subsection 3.3 can be used to solve the steady state problem. We sketch the discretization used in [START_REF] Sun | Asymptotic preserving unified gas kinetic scheme for grey radiative transfer equations[END_REF] in the subsequent part and for more details and the proof of AP, one can refer to [START_REF] Sun | Asymptotic preserving unified gas kinetic scheme for grey radiative transfer equations[END_REF].

4.1.

The AP UGKS for the grey radiative transport equation. The idea is to first introduce ψ = acT 4 and rewrite the system (1.3) into (4.1)

     1 c ∂ t I + 1 ε v∂ x I = σ ε 2 1 |V | ψ -I + q(x, v), ε 2 ∂ t ψ = β(ψ)σ(ρ -ψ),
where ρ := V Idv is the total radiation intensity and β is a function of ψ given by

β(ψ) = 4ac C v ψ ac 3/4
.

Then we take the integral with respect to v in (4.1) and get a system for ρ, ψ, with an advection term depending on V vI dv. Similar as the UGKS for the chemotaxis model, the first step is to get a prediction for ρ n+1 by (4.2)

       ρ n+1 i = ρ n i + ∆t ∆x F n i-1/2 -F n i+1/2 + σc∆t ε 2 (ψ n+1 i -ρ n+1 i ) + c∆t V q i dv, ψ n+1 i = ψ n i + σβ n+1 i ∆t ε 2 ρ n+1 i -ψ n+1 i ,
where

F n i+1/2 := c |V | V 1 ε∆t t n+1 t n vI(x i+1/2 , v, t) dt dv, q i (v) := 1 ∆x x i+1/2
x i-1/2 q(x, v)dx.

After discretizing the velocity space by discrete-ordinate method, the numerical flux F n i-1/2 and source term are approximated by

F n i-1/2 ≈ A N m=-N ω m µ m (I n- i-1/2,m 1 µm>0 + I n+ i-1/2,m 1 µm<0 ) + 2D 3∆x (ψ n+1 i -ψ n+1 i-1 ) +B N m=-N ω m µ 2 m (δ x I n i-1,m 1 µm>0 + δ x I n i,m 1 µm<0 ) + |V |ε 2 C σ m∈Vm ω m µ m q i,m , V q i dv ≈ m∈Vm ω m q i,m = m∈Vm ω m q i (µ m ).
Here 1 µm<0 , 1 µm>0 are the characteristic functions, δ x I n i,m are the approximations to ∂ x I(x i , µ m , t n ) by the slope limiters and I n± i-1/2,m are given by

I n- i-1/2,m = I n i-1,m + ∆x 2 δ x I n i-1,m , I n+ i-1/2,m = I n i,m - ∆x 2 δ x I n i,m .
The coefficients A, B, C, D depend on ε and ∆t, which will be specified later.

After obtaining the macroscopic variables ρ n+1 i , ψ n+1 i from (4.2), I n+1 i,m are computed by (4.3)

I n+1 i,m = I n i,m + ∆t ∆x Φ n i-1/2,m -Φ n i+1/2,m + σc∆t ε 2 ψ n+1 i |V | -I n+1 i,m + c∆tq i,m ,
where Φ n i-1/2,m are given by

Φ n i-1/2,m =Aµ m I n- i-1/2,m 1 µm>0 + I n+ i-1/2,m 1 µm<0 + Cµ m ψ n+1 i-1/2 +D µ 2 m δ x ψ n+1,L i-1/2 1 µm>0 + µ 2 m δ x ψ n+1,R i-1/2 1 µm<0 +B µ 2 m δ x I n i-1,m 1 µm>0 + µ 2 m δ x I n i,m 1 µm<0 +Eδ t ψ n+1 i-1/2 + |V |ε 2 C σ µ m q i,m , with ψ n+1 i+1/2 , δ x ψ n+1,R i-1/2 , δ x ψ n+1,L i-1/2 , δ t ψ n+1 i-1/2 being defined by ψ n+1 i+1/2 = 1 2 (ψ n+1 i + ψ n+1 i+1 ). δ x ψ n+1,L i-1/2 = ψ n+1 i-1/2 -ψ n+1 i-1 ∆x/2 , δ t ψ n+1 i-1/2 = ψ n+1 i-1/2 -ψ n i-1/2 ∆t , δ x ψ n+1,R i-1/2 = ψ n+1 i -ψ n+1 i-1/2 ∆x/2 . δ x ψ n+1,L i-1/2 and δ x ψ n+1,R i-1/2
are the same and we use different notations in order to be consistent with notations in [START_REF] Sun | Asymptotic preserving unified gas kinetic scheme for grey radiative transfer equations[END_REF]. In [START_REF] Sun | Asymptotic preserving unified gas kinetic scheme for grey radiative transfer equations[END_REF], different notations are used for δ x ψ n+1,R i-1/2 and δ x ψ n+1,L i-1/2 in order to consider different approximations for ψ n+1 i-1/2 . Some of them can lead to the case when δ x ψ n+1,R i-1/2 and δ ψ n+1,L i-1/2 are different. Let ν = cσ ε 2 , the coefficients A, B, C, D, E appeared in all above formula are given by

A = c ε∆tν (1 -exp -(ν∆t)) , C = c 2 σ |V |ε 3 ∆tν ∆t - 1 ν (1 -exp (-ν∆t)) , D = - c 3 σ |V |ε 4 ∆tν 2 ∆t(1 + exp (-ν∆t)) - 2 ν (1 -exp (-ν∆t)) , B = - c 2 ε 2 ν 2 ∆t (1 -exp (-ν∆t) -ν∆t exp (-ν∆t)) , E = c 2 σ |V |ε 3 ν 3 ∆t 1 -exp (-ν∆t) -ν∆t exp (-ν∆t) - 1 2 (ν∆t) 2 .

4.2.

Steady states for the radiative transport equation. We use the discrete ordinate method similar as in subsection 3.3 and refer to [START_REF] Gosse | Transient radiative transfer in the grey case: Well-balanced and asymptoticpreserving schemes built on case's elementary solutions[END_REF] for more details. The one dimensional discrete ordinate steady state problem for (4.1) reads (4.4)

µ m ∂ x I m = σ ε ( m ω m I m -I m ) + εq m .
We use the same AP scheme as in subsection 3.3 and the above equation falls into the category of Case II. The only difference is that there exists a source term q m , which can be easily built into the scheme by constructing an approximated particular solution. The general solution with source term q m can be given by the summation of the general solution of the homogeneous problem and the approximated particular solution. Then the coefficients in front of the general solutions of the homogeneous equation can be determined by the continuity of solution at the cell edge. For more details, we refer to [START_REF] Shi | A uniformly second order numerical method for the one-dimensional discrete-ordinate transport equation and its diffusion limit with interface[END_REF].

Numerical simulations

In this section, we apply the well-balanced and AP scheme designed in Section 2 to two test cases: kinetic chemotaxis and radiative transport models. We impose specular boundary conditions for f ε and Dirichlet conditions for S ε . The initial condition is composed of two bumps in x located at ±0.65 such that f 0 (x, v) = 5 exp -10(x -0.65) 2 -10(x + 0.65) 2 exp -20(v -0.5) 2 -20(v + 0.5) 2 .

To ensure the stability of the numerical scheme, the time step ∆t is chosen as follows ∆t = ∆x 2 , for ε < ∆x 2 , ǫ∆x, else.

The requirement for the time step can be improved if the limiting discretization is implicit, the strategy of constructing an implicit limiting discretization has been discussed in [START_REF] Sun | Asymptotic preserving unified gas kinetic scheme for grey radiative transfer equations[END_REF]. The goal of our present paper is to illustrate the idea of the AP-WB framework and thus we keep using the simple scheme.

5.1.1. AP property. In order to verify the AP property of our scheme, we show that the convergence order of the scheme is independent of ε. In figure 1, we plot e ∆x := f ∆x (t) -f 2∆x (t) 1 with respect to ∆x for different values of ε and uniform convergence with respect to ǫ can be observed numerically. This guarantees the AP property.

5.1.2. Model convergence in ε. As predicted by the theoretical analysis, when ε → 0,

f ε → ρ 0 1 v∈V ,
where ρ 0 is the solution to the limiting Keller-Segel equation. To verify the above convergence, the total densities ρ ε are displayed in Figure 2 for different values of ε ranging from 10 -1 to 10 -6 . As a comparison, we plot the solution of the limiting model in Figure 2 and can observe that ρ ε get closer to ρ 0 as ε goes to zero. The convergence order in ε can be seen in Figure 3, where we have plotted f ε (x, v) -ρ ε (x)1 v∈V 2 for different ε in the logarithm scale. The numerical results indicate that the order of convergence is close to 1. 5.1.3. WB property. In order to check the WB property, we calculate the long time behavior of the kinetic scheme for ε = 1. The steady state is characterised by the vanishing of the macroscopic flux J. A way to measure whether the steady state is reached is to look at the evolution of the residues r n at a given time t n defined by

r n := m∈V ω m |f n+1 (µ m ) -f n (µ m )| 2 .
Figure 4 shows the decrease of the residues and its stabilization when the steady state is reached. In addition, the flux J is of order ∆x 2 at this steady state as displayed in Figure 5. When ε = 1, we can have analytical nonzero steady state for this choice of q(x, v). Though q(x, v) in (5.1) can take negative values, which is non-physical, it can be

AP property.

In this part, we show that the convergence order of the scheme is independent of ε. In Figure 7, we plot e ∆x := max x∈(0,1) I ∆x (x, t) -I 2∆x (x, t) 1 with respect to ∆x for different values of ε. Uniform convergence can be observed which confirms the AP property of the scheme.

5.2.2.

Model convergence with respect to ε. The theoretical analysis indicates that

I ε → ψ 0 2 1 v∈V , ε → 0,
where ψ 0 = ac(T 0 ) 4 and T 0 is the solution to the limiting equation (1.4). To verify the above convergence, the total densities ρ ε are displayed in Figure 7 for different values of ε ranging from 10 -1 to 10 -6 . As a comparison, we plot the solution of the limiting model in Figure 2 and can observe that I ε get closer to ψ 0 /2 as ε goes to zero.

The convergence order in ε can be seen in Figure 8, where we have plotted I ε (x, v) -ρ ε 1 v∈V 2 , for different ε in the logarithm scale. The numerical results indicate that the order of convergence is close to 0.6.

WB property.

For the source term q as in (5.1), when ε = 1, the steady state is given by I(x, v) = |v|x, ρ = x. The time evolution of the residues r n is given in Figure 9, it decrease and stabilize at T = 1. In Figure 10, the captured steady state is plotted.

Conclusion

A general framework of developing AP and WB schemes is proposed in this paper. Two parameter regimes are considered: the (advection) diffusion type (linear or nonlinear) limit when ε << 1 and the long time behavior of the hyperbolic equations with source terms when ε = O(1). The framework is composed of two steps. We first calculate a prediction with an AP scheme for the time evolutionary problem and then use the prediction to set up the inflow boundary conditions for the steady state problem in each cell. The steady state problem in each cell is solved by an AP scheme for the steady state equation, then the obtained outflows in each cell are used for modifying the numerical flux to achieve WB. The details of the discretizations for transport equations for chemotaxis and gray radiative transfer are given. The numerical results in section 5 verify the AP and WB properties of our proposed scheme. In particular, the AP schemes for the time evolutionary transport equation for chemotaxis and the corresponding steady state equation are new.

One may argue about the computational cost of our proposed framework, the advantage is that, any AP schemes for the time evolutionary problem and steady state problem are applicable, thus we can choose those simple ones. First of all, the key principle underlying most WB schemes consists of using values constructed from a local discrete equilibrium, therefore if we are interested in unknown steady states, it is hard to avoid the steady problem step. On the other hand, in the prediction step, we can use those AP scheme that can implement the implicit stiff terms explicitly as in [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF], or at least more efficient than using the Newton type solvers for nonlinear algebraic systems. (2.3) can be updated explicitly afterwards, the total computational cost is composed of some explicit calculations plus solving the steady state problems.

The proposed framework can be extended to the hyperbolic system with relaxations, we will pursue more applications in the future.

  ρ n+1 i and f n+1 i are updated as follows (3.1)

  Case I: If (3.15) has 2N different simple roots, we can find 2N linearly independent eigenfunctions of the form l m exp(-ζ ε x). Then the general solution to (3.11) can be written as

5. 1 .

 1 Chemotaxis model. The first simulation concerns the chemotaxis model(1.1). Simulations are set on x ∈ [-1, 1] and v ∈ [-1, 1]. We take N v = 32 points for the velocity and N x = 500 points for the space discretization. Parameters in (1.1) are chosen as in[START_REF] Gosse | A well-balanced scheme for kinetic models of chemotaxis derived from onedimensional local forward-backward problems[END_REF] such that χ S = 1, D = 15, β = 60, α = 3.

Figure 1 .

 1 Figure1. Convergence order with respect to ∆x. Here the log-log plot of the mesh sizes with respect to the L 1 norm of the probability density error in all directions is displayed. The mesh sizes are 1/20, 1/40, 1/80, 1/160 and e ∆x := max x∈(-1,1) f ∆x (x, t)-f 2∆x (x, t) 1 are plotted for different ε.

6 Figure 2 . 2 0Figure 3 .

 6223 Figure 2. Plot of the macroscopic density for different values of ε. ε = 0 is the solution to the macroscopic Keller-Segel equation.

Figure 4 .Figure 5 . 5 . 2 .

 4552 Figure 4. The evolution of residues with respect to time.

Figure 6 .

 6 Figure 6. Convergence order with respect to ∆x. Here the log-log plot of the mesh sizes with respect to the L 1 norm of the probability density error in all directions is displayed. The mesh sizes are 1/10, 1/20, 1/40, 1/80 and e ∆x := max x∈(0,1) I ∆x (x, t) -I 2∆x (x, t) 1 are plotted for different ε.

ρ 1 ε = 0 ε = 10 - 1 ε = 10 - 2 ε = 10 - 5 ε

 10101102105 Plot of macroscopic density rho for radiative equation at T=0.= 10 -6

Figure 7 .

 7 Figure 7. Total density for different values of ε. ε = 0 corresponds to the limiting nonlinear diffusion model (1.4).

1 0Figure 8 .

 18 Figure 8. Convergence order in ε of the grey radiative transport model (1.3).

Plot of the residues r in the logarithm scaleFigure 9 .Figure 10 .

 910 Figure 9. Evolution of residues with respect to time.