Convergence of a finite-volume scheme

# Convergence of a finite-volume scheme for a degenerate cross-diffusion model for ion transport

Clément Cancès Inria, Univ. Lille, CNRS, UMR 8524 - Laboratoire Paul Painlevé, F-59000 Lille Claire Chainais-Hillairet Univ. Lille, CNRS, UMR 8524, Inria - Laboratoire Paul Painlevé, F-59000 Lille Anita Gerstenmayer Institute for Analysis and Scientific Computing, Vienna University of Technology, Wiedner Hauptstraße 8–10, 1040 Wien, Austria  and  Ansgar Jüngel Institute for Analysis and Scientific Computing, Vienna University of Technology, Wiedner Hauptstraße 8–10, 1040 Wien, Austria
July 7, 2019
###### Abstract.

An implicit Euler finite-volume scheme for a degenerate cross-diffusion system describing the ion transport through biological membranes is analyzed. The strongly coupled equations for the ion concentrations include drift terms involving the electric potential, which is coupled to the concentrations through the Poisson equation. The cross-diffusion system possesses a formal gradient-flow structure revealing nonstandard degeneracies, which lead to considerable mathematical difficulties. The finite-volume scheme is based on two-point flux approximations with “double” upwind mobilities. It preserves the structure of the continuous model like nonnegativity, upper bounds, and entropy dissipation. The degeneracy is overcome by proving a new discrete Aubin-Lions lemma of “degenerate” type. Under suitable assumptions, the existence and uniqueness of bounded discrete solutions, a discrete entropy inequality, and the convergence of the scheme is proved. Numerical simulations of a calcium-selective ion channel in two space dimensions indicate that the numerical scheme is of first order.

###### Key words and phrases:
Ion transport, finite-volume method, gradient flow, entropy method, existence of discrete solutions, convergence of the scheme, calcium-selective ion channel.
###### 2000 Mathematics Subject Classification:
65M08, 65M12, 35K51, 35K65, 35Q92.
The authors have been supported by the Austrian-French Program Amadée of the Austrian Exchange Service (ÖAD). The work of the first and second authors is supported by the LABEX CEMPI (ANR-11-LABX-0007-01). The third and fourth authors acknowledge partial support from the Austrian Science Fund (FWF), grants P27352, P30000, F65, and W1245

## 1. Introduction

The ion transport through biological channels plays an important role in all living organisms. On a macroscopic level, the transport can be described by nonlinear partial differential equations for the ion concentrations (or, more precisely, volume fractions) and the surrounding electric potential. A classical model for ion transport are the Poisson-Nernst-Planck equations [25], which satisfy Fick’s law for the fluxes. However, this approach does not include size exclusion effects in narrow ion channels. Taking into account the finite size of the ions, one can derive from an on-lattice model in the diffusion limit another set of differential equations with fluxes depending on the gradients of all species [10, 27]. These nonlinear cross-diffusion terms are common in multicomponent systems [23, Chapter 4]. In this paper, we propose an implicit Euler finite-volume discretization of the resulting cross-diffusion system. The scheme is designed in such a way that the nonnegativity and upper bound of the concentrations as well as the entropy dissipation is preserved on the discrete level.

More specifically, the evolution of the concentrations and fluxes of the th ion species is governed by the equations

 (1) ∂tui+divFi=0,Fi=−Di(u0∇ui−ui∇u0+u0uiβzi∇Φ)in Ω, t>0,

for , where is the concentration (volume fraction) of the electro-neutral solvent, is a diffusion coefficient, is the (scaled) inverse thermal voltage, and the charge of the th species. Observe that we assumed Einstein’s relation which says that the quotient of the diffusion and mobility coefficients is constant, and we call this constant . The electric potential is determined by the Poisson equation

 (2) −λ2ΔΦ=n∑i=1ziui+fin Ω,

where is the (scaled) permittivity constant and is a permanent background charge density. Equations (1) and (2) are solved in a bounded domain .

In order to match experimental conditions, the boundary is supposed to consist of two parts, the insulating part , on which no-flux boundary conditions are prescribed, and the union of boundary contacts with external reservoirs, on which the concentrations are fixed. The electric potential is prescribed at the electrodes on . This leads to the mixed Dirichlet-Neumann boundary conditions

 (3) Fi⋅ν=0on ΓN, ui=¯¯¯uion ΓD,i=1,…,n, (4) ∇Φ⋅ν=0on ΓN, Φ=¯¯¯¯Φon ΓD,

where the boundary data and can be defined on the whole domain . Finally, we prescribe the initial conditions

 (5) ui(⋅,0)=uIiin Ω, i=1,…,n.

The main mathematical difficulties of equations (1) are the strong coupling and the fact that the diffusion matrix , defined by for and is not symmetric and not positive definite. It was shown in [10, 22] that system (1) possesses a formal gradient-flow structure. This means that there exists a (relative) entropy functional with the entropy density

 h(u)=n∑i=0∫ui¯¯¯uilogs¯¯¯uids+βλ22|∇(Φ−¯¯¯¯Φ)|2,

where and , such that (1) can be formally written as

 ∂tui=div(n∑j=1Bij∇wj),

where , for provide a diagonal positive definite matrix, and are the entropy variables, defined by

 ∂h∂ui=wi−¯¯¯¯wi,where wi=loguiu0+βziΦ,¯¯¯¯wi=log¯¯¯ui¯¯¯u0+βzi¯¯¯¯Φ,i=1,…,n.

We refer to [20, Lemma 7] for the computation of .

The entropy structure of (1) is useful for two reasons. First, it leads to bounds for the concentrations. Indeed, the transformation to entropy variables can be inverted, giving with

 ui(w,Φ)=exp(wi−βziΦ)1+∑nj=1exp(wj−βzjΦ),i=1,…,n.

Then is positive and bounded from above, i.e.

 (6) u∈D={u∈(0,1)n:n∑i=1ui<1}.

This yields bounds without the use of a maximum principle. Second, the entropy structure leads to gradient estimates via the entropy inequality

 dHdt+12∫Ωn∑i=1Diu0ui|∇wi|2dx≤C,

where the constant depends on the Dirichlet boundary data. Because of

 (7) n∑i=1u0ui∣∣∣∇loguiu0∣∣∣2=4n∑i=1u0|∇u1/2i|2+4|∇u1/20|2+|∇u0|2,

we achieve gradient estimates for and . Since may vanish locally, this does not give gradient bounds for , which expresses the degenerate nature of the cross-diffusion system. As a consequence, the flux has to be formulated in the terms of gradients of and only, namely

 (8) Fi=−Di(u1/20∇(u1/20ui)−3u1/20ui∇u1/20+u0uiβzi∇Φ).

The challenge is to derive a discrete version of this formulation. It turns out that (24) below is the right formulation in our context (assuming vanishing drift parts).

Our aim is to design a numerical approximation of (1) which preserves the structural properties of the continuous equations. This suggests to use the entropy variables as the unknowns, as it was done in our previous work [20] with simulations in one space dimension. Unfortunately, we have not been able to perform a numerical convergence analysis with these variables. The reason is that we need discrete chain rules in order to formulate (7) on the discrete level and these discrete chain rules seem to be not easily available. Therefore, we use the original variables for the numerical discretization. Interestingly, we are still able to prove that the scheme preserves the nonnegativity, upper bound, and entropy inequality. However, the upper bound comes at a price: We need to assume that all diffusion coefficients are the same. Under this assumption, solves a drift-diffusion equation for which the (discrete) maximum principle can be applied. It is not surprising that the bound can be shown only under an additional condition, since cross-diffusion systems usually do not allow for a maximum principle.

The key observation for the numerical discretization is that the fluxes can be written on each cell in a “double” drift-diffusion form, i.e., both and have the structure , where is the diffusion term and is the drift term. We discretize and by using a two-point flux aproximation with “double” upwind mobilities.

Our analytical results are stated and proved for no-flux boundary conditions on . Mixed Dirichlet-Neumann boundary conditions could be prescribed as well, but the proofs would become even more technical. The main results are as follows.

• We prove the existence of solutions to the fully discrete numerical scheme (Theorem 1). If the drift part vanishes, the solution is unique. The existence proof uses a topological degree argument in finite space dimensions, while the uniqueness proof is based on the entropy method of Gajewski [18], recently extended to cross-diffusion systems [20, 28].

• Thanks to the “double” upwind structure, the scheme preserves the nonnegativity and upper bound for the concentrations (at least if for all ). Moreover, convexity arguments show that the discrete entropy is dissipated with a discrete entropy production analogous to (7) (Theorem 2). The proof of the discrete entropy inequality only works if the drift term vanishes, since we need to control a discrete version of the sum from below; see the discussion after Theorem 2.

• The discrete solutions converge to the continuous solutions to (1) as the mesh size tends to zero (Theorem 3). The proof is based on a priori estimates obtained from the discrete entropy inequality. The compactness is derived from a new discrete Aubin-Lions lemma, which takes into account the nonstandard degeneracy of the equations; see Lemma 10 in the appendix.

• Numerical experiments for a calcium-selective ion channel in two space dimensions show the dynamical behavior of the solutions and their large-time asymptotics to the equilibrium. The tests indicate that the order of convergence in the norm is one.

In the literature, there exist some results on finite-volume schemes for cross-diffusion systems. An upwind two-point flux approximation similar to our discretization was recently used in [1] for a seawater intrusion cross-diffusion model. A two-point flux approximation with a nonlinear positivity-preserving approximation of the cross-diffusion coefficients, modeling the segregation of a two-species population, was suggested in [4], assuming positive definiteness of the diffusion matrix. The Laplacian structure of the population model (still for positive definite matrices) was exploited in [24] to design a convergent linear finite-volume scheme, which avoids fully implicit approximations. A semi-implicit finite-volume discretization for a biofilm model with a nonlocal time integrator was proposed in [26]. Finite-volume schemes for cross-diffusion systems with nonlocal (in space) terms were also analyzed; see, for instance, [3] for a food chain model and [2] for an epidemic model. Moreover, a finite-volume scheme for a Keller-Segel system with additional cross diffusion and discrete entropy dissipation property was investigated in [7]. All these models, however, do not include volume filling and do not possess the degenerate structure explained before.

The paper is organized as follows. The numerical scheme and the main results are presented in Section 2. In Section 3, the existence and uniqueness of bounded discrete solutions are shown. We prove the discrete entropy inequality and further a priori estimates in Section 4, while Section 5 is concerned with the convergence of the numerical scheme. Numerical experiments are given in Section 6 in order to illustrate the order of convergence and the long time behavior of the scheme. For the compactness arguments, we need two discrete Aubin-Lions lemmas which are proved in the appendix.

## 2. Numerical scheme and main results

### 2.1. Notations and definitions

We summarize our general hypotheses on the data:

{labeling}

(A44)

Domain: ( or ) is an open, bounded, polygonal domain with , .

Parameters: , , , and , .

Background charge: .

Initial and boundary data: , satisfy , and , in for , and .

For our main results, we need additional technical assumptions:

{labeling}

(A44)

, i.e., we impose no-flux boundary conditions on the whole boundary.

The diffusion constants are equal, for .

The drift terms are set to zero, .

###### Remark 1 (Discussion of the assumptions).

Assumption (A1) is supposed for simplicity only. Mixed Dirichlet-Neumann boundary conditions can be included in the analysis (see, e.g., [20]), but the proofs become even more technical. Mixed boundary conditions are chosen in the numerical experiments; therefore, the numerical scheme is defined for that case. Assumption (A2) is needed for the derivation of an upper bound for the solvent concentration. Indeed, when for all , summing (1) over gives

 ∂tu0=Ddiv(∇u0−u0w∇Φ),where w=βn∑i=1ziui.

On the discrete level, we replace by an upwind approximation. This allows us to apply the discrete maximum principle showing that and hence with defined in (6). Finally, Assumption (A3) is needed to derive a discrete version of the entropy inequality. Without the drift terms, the upwinding value does not depend on the index of the species, which simplifies some expressions; see Remark 2. ∎

For the definition of the numerical scheme for (1)-(2), we need to introduce a suitable discretization of the domain and the interval . For simplicity, we consider a uniform time discretization with time step , and we set for , where , are given and . The domain is discretized by a regular and admissible triangulation in the sense of [16, Definition 9.1]. The triangulation consists of a family of open polygonal convex subsets of (so-called cells), a family of edges (or faces in three dimensions), and a family of points associated to the cells. The admissibility assumption implies that the straight line between two centers of neighboring cells is orthogonal to the edge between two cells and . The condition is satisfied by, for instance, triangular meshes whose triangles have angles smaller than [16, Examples 9.1] or Voronoi meshes [16, Example 9.2].

We assume that the family of edges can be split into internal and external edges with and . Each exterior edge is assumed to be an element of either the Dirichlet or Neumann boundary, i.e. . For given , we define the set of the edges of , which is the union of internal edges and edges on the Dirichlet or Neumann boundary, and we set .

The size of the mesh is defined by . For with , we denote by the Euclidean distance between and , while for , we set . For a given edge , the transmissibility coefficient is defined by

 (9) τσ=m(σ)dσ,

where denotes the Lebesgue measure of .

We impose a regularity assumption on the mesh: There exists such that for all and , it holds that

 (10) d(xK,σ)≥ζdσ.

This hypothesis is needed to apply discrete functional inequalities (see [6, 16]) and a discrete compactness theorem (see [19]).

It remains to introduce suitable function spaces for the numerical discretization. The space of piecewise constant functions is defined by

 HT={v:¯¯¯¯Ω→R:∃(vK)K∈T⊂R, v(x)=∑K∈TvK1K(x)}.

The (squared) discrete norm on this space is given by

 (11) ∥v∥21,T=∑σ=K|L∈Eintτσ(vK−vL)2+∑K∈Tm(K)v2K.

The discrete norm is the dual norm with respect to the scalar product,

 (12) ∥v∥−1,T=sup{∫Ωvwdx:w∈HT, ∥w∥1,T=1}.

Then

 ∣∣∣∫Ωvwdx∣∣∣≤∥v∥−1,T∥w∥1,Tfor v,w∈HT.

Finally, we introduce the space of piecewise constant in time functions with values in ,

 HT,△t={v:¯¯¯¯Ω×[0,T]→R:∃(vk)k=1,…,N⊂HT, v(x,t)=N∑k=1vk(x)1(tk−1,tk)(t)},

equipped with the discrete norm

 ∥v∥1,T,△t=(N∑k=1△t∥vk∥21,T)1/2.

For the numerical scheme, we introduce some further definitions. Let with values on the Dirichlet boundary (). Then we introduce

 (13) DK,σ(ui)=ui,K,σ−ui,K, whereui,K,σ=⎧⎪ ⎪⎨⎪ ⎪⎩ui,L for % σ∈Eint, σ=K|L,¯¯¯ui,σ for σ∈EDext,K,ui,K for σ∈ENext,K,¯¯¯ui,σ=1m(σ)∫σ¯¯¯uids.

The numerical fluxes should be consistent approximations to the exact fluxes through the edges . We impose the conservation of the numerical fluxes for edges , requiring that they vanish on the Neumann boundary edges, for . Then the discrete integration-by-parts formula becomes for

 ∑K∈T∑σ∈EKFK,σuK=−∑σ∈EFK,σDK,σ(u)+∑σ∈EDextFK,σuK,σ.

When , this formula simplifies to

 (14) ∑K∈T∑σ∈EKFK,σuK=∑σ=K|L∈EintFK,σ(uK−uL).

### 2.2. Numerical scheme

We need to approximate the initial, boundary, and given functions on the elements and edges :

 uIi,K =1m(K)∫KuIi(x)dx, fK =1m(K)∫Kf(x)dx, ¯¯¯ui,σ =1m(σ)∫σ¯¯¯uids, ¯¯¯¯Φσ =1m(σ)∫σ¯¯¯¯Φds,

and we set and .

The numerical scheme is as follows. Let , , , and be given. Then the values are determined by the implicit Euler scheme

 (15) m(K)uki,K−uk−1i,K△t+∑σ∈EKFki,K,σ=0,

where the fluxes are given by the upwind scheme

 (16) Fki,K,σ=−τσDi(uk0,σ%DK,σ(uki)−uki,σ(DK,σ(uk0)−ˆuk0,σ,iβziDK,σ(Φk))),

where is defined in (9),

 (17) uk0,K=1−n∑i=1uki,K,uk0,σ=max{uk0,K,uk0,L}, (18) uki,σ={uki,Kif Vki,K,σ≥0,uki,K,σif Vki,K,σ<0,,ˆuk0,σ,i={uk0,Kif ziDK,σ(Φk)≥0,uk0,K,σif ziDK,σ(Φk)<0,,

and is the “drift part” of the flux,

 (19) Vki,K,σ=DK,σ(uk0)−ˆuk0,σ,iβziDK,σ(Φk)

for . Observe that we employed a double upwinding: one related to the electric potential, defining , and another one related to the drift part of the flux, . The potential is computed via

 (20) −λ2∑σ∈EKτσDK,σ(Φk)=m(K)(n∑i=1ziuki,K+fK).

We recall that the numerical boundary conditions are given by and for .

We denote by , the functions in associated to the values and , respectively. Moreover, when dealing with a sequence of meshes and a sequence of time steps , we set , .

###### Remark 2 (Simplified numerical scheme).

When Assumptions (A1)-(A3) hold, the numerical scheme simplifies to

 (21) m(K)uki,K−uk−1i,K△t+∑σ∈EK,intFki,K,σ=0, (22) Fki,K,σ=−τσD(uk0,σ(uki,L−uki,K)−uki,σ(uk0,L−uk0,K)),

where and are defined in (17), and the definition of simplifies to

 uki,σ={uki,Kif uk0,K−uk0,L≤0,uki,Lif uk0,K−uk0,L>0.

In the definition of , the upwinding value does not depend on anymore such that

 (23) n∑i=0uki,σ=1+max{uk0,K,uk0,L}−min{uk0,K,uk0,L}=1+|uk0,K−uk0,L|.

This property is needed to control the sum from below in the proof of the discrete entropy inequality; see (35). Finally, we are able to reformulate the discrete fluxes such that we obtain a discrete version of (8) (without the drift part):

 (24) Fi,K,σ=τσD{u1/20,σ(u1/20,Kui,K−u1/20,Lui,L)−ui,σ(u1/20,K−u1/20,L)(u1/20,σ+2u1/20,K+u1/20,L2)}.

This formulation is needed in the convergence analysis. ∎

### 2.3. Main results

Since our scheme is implicit and nonlinear, the existence of an approximate solution is nontrivial. Therefore, our first result concerns the well-posedness of the numerical scheme.

###### Theorem 1 (Existence and uniqueness of solutions).

Let (H1)-(H4) and (A2) hold. Then there exists a solution to scheme (15)-(20) satisfying and, if the initial data lie in , . If additionally Assumptions (A1) and (A3) hold, the solution is unique.

Assumption (A2) is needed to show that is nonnegative. Indeed, summing (15) and (16) over , we obtain

 m(K)uk0,K−uk−10,K△t=−∑σ∈EKτσ(uk0,σDK,σ(n∑i=1Diuki)−n∑i=1Diuki,σVki,K,σ).

Under Assumption (A2), it follows that , and we can apply the discrete minimum principle, which then implies an bound for . This bound allows us to apply a topological degree argument; see [13, 14]. For the uniqueness proof, we additionally need Assumption (A3), since we use the entropy method of Gajewski [18], and it seems that this method cannot be applied to cross-diffusion systems including drift terms [28]. The idea is to prove first the uniqueness of , which solves a discrete nonlinear equation, and then to show the uniqueness of for by introducing a semimetric for two solutions and and showing that it is monotone in , such that a discrete Gronwall argument implies that .

The second result shows that the scheme preserves a discrete version of the entropy inequality.

###### Theorem 2 (Discrete entropy inequality).

Let Assumptions (H1)-(H4) and (A1)-(A3) hold. Then the solution to scheme (21)-(22) constructed in Theorem 1 satisfies the discrete entropy inequality

 (25) Hk−Hk−1△t+Ik≤0,

with the discrete entropy

 (26) Hk=∑K∈Tm(K)n∑i=0(uki,K(loguki,K−1)+1)

and the discrete entropy production

 Ik xx+4((uk0,K)1/2−(uk0,L)1/2)2+(uk0,K−uk0,L)2).

Assumption (A3) is required to estimate the expression . In the continuous case, this sum equals . On the discrete level, this identity cannot be expected since the value of depends on the upwinding value; see (18). If the drift part vanishes, the upwinding value does not depend on , as mentioned in Remark 2, and we can derive the estimate ; see Section 4.1. Note that the entropy production is the discrete counterpart of (7).

The main result of this paper is the convergence of the approximate solutions to a solution to the continuous cross-diffusion system.

###### Theorem 3 (Convergence of the approximate solution).

Let (H1)-(H4) and (A1)-(A3) hold and let and be sequences of admissible meshes and time steps, respectively, such that and as . Let be the solution to (21)-(22) constructed in Theorem 1. Then there exist functions , satisfying ,

 u1/20, u1/20ui∈L2(0,T;H1(Ω)),i=1,…,n, u1/20,m→u1/20, u1/20,mui,m→u1/20uistrongly in L2(Ω×(0,T)),

where is a weak solution to (1), (3)-(5) (with ), i.e., for all and ,

 (27) ∫T0∫Ωui∂tϕdxdt+∫ΩuIiϕ(⋅,0)dx=D∫T0∫Ωu1/20(∇(u1/20ui)−3ui∇u1/20)⋅∇ϕdxdt.

The compactness of the concentrations follows from the discrete gradient estimates derived from the entropy inequality (25), for which we need Assumption (A3). By the discrete Aubin-Lions lemma [17], we conclude the strong convergence of the sequence . The difficult part is to show the strong convergence of , since there is no control on the discrete gradient of . The idea is to apply a discrete Aubin-Lions lemma of “degenerate” type, proved in Lemma 10 in the appendix.

## 3. Existence and uniqueness of approximate solutions

### 3.1. L∞ bounds and existence of solutions

In order to prove the existence of solutions to (15)-(20), we first consider a truncated problem. This means that we truncate the expressions in (18); more precisely, we consider scheme (15), (16), and (20) with

 uk0,K=1−n∑i=1(uki,K)+,uk0,σ=max{0,uk0,K,uk0,K,σ}, (28) ˆuk0,σ,i={(uk0,K)+%ifziDK,σ(Φk)≥0,(uk0,K,σ)+if ziDK,σ(Φk)<0, uki,σ={(uki,K)+if Vki,K,σ≥0,(uki,K,σ)+if Vki,K,σ<0,

where for and . We show that this truncation is, in fact, not needed if the initial data are nonnegative. In the following let (H1)-(H4) hold.

###### Lemma 4 (Nonnegativity of uki).

Let be a solution to (15), (16), (20), and (28). Then for all , , and . If and then also for all , .

###### Proof.

We proceed by induction. For , the nonnegativity holds because of our assumptions on the initial data. Assume that for all . Then let for some and assume that . The scheme writes as

 (29) m(K)uki,K−uk−1i,K△t=∑σ∈EKτσDi(uk0,σDK,σ(uki)−uki,σVki,K,σ).

By assumption, . If , we have and if , it follows that . Hence, the right-hand side of (29) nonnegative. However, the left-hand side is negative, which is a contradiction. We infer that and consequently, for all . When the initial data are positive, similar arguments show the positivity of for . ∎

We are able to show the nonnegativity of only if the diffusion coefficients are the same. The reason is that we derive an equation for by summing (15) for , and this gives an equation for only if for all .

###### Lemma 5 (Nonnegativity of uk0).

Let Assumption (A2) hold and let be a solution to (15), (16), (20), and (28). Then for all , . If and then also for all , .

###### Proof.

Again, we proceed by induction. The case follows from the assumptions. Assume that for all . Then let for some and assume that . Summing equations (15) from , we obtain

 m(K)uk0,K−uk−10,K△t (30) ≥−D∑σ∈EKτσn∑i=1βziˆuk0,σ,iDK,σ(Φk),

since and by construction and because of the minimality property of . The remaining expression is nonnegative:

However, the left-hand side of (30) is negative, by induction hypothesis, which gives a contradiction. ∎

Lemmas 4 and 5 imply that we may remove the truncation in (28). Moreover, by definition, we have such that or, if the initial and boundary data are positive,