Analysis of a Population Model with Strong Cross-Diffusion in an Unbounded Domain

(1)

Analysis of a Population Model with Strong Cross-Diffusion in an Unbounded Domain

Michael Dreher

Konstanzer Schriften in Mathematik und Informatik Nr. 216, Juni 2006

ISSN 1430–3558

c Fachbereich Mathematik und Statistik

c Fachbereich Informatik und Informationswissenschaft Universit¨at Konstanz

Fach D 188, 78457 Konstanz, Germany Email: preprints@informatik.uni–konstanz.de

WWW: http://www.informatik.uni–konstanz.de/Schriften/

Konstanzer Online-Publikations-System (KOPS) URL: http://www.ub.uni-konstanz.de/kops/volltexte/2007/2236/

URN: http://nbn-resolving.de/urn:nbn:de:bsz:352-opus-22364

(2)

Strong Cross-Diffusion in an Unbounded Domain

Michael Dreher^∗

Abstract

We study a parabolic population model in the full space and prove the global in time existence of a weak solution. This model consists of two strongly coupled diffusion equations describing the population densities of two competing species. The system features intrinsic growth, inter- and intra-specific competition of the species, as well as diffusion, cross-diffusion and self-diffusion, and drift terms related to varying environment quality. The cross-diffusion terms can be large, making the system non-parabolic for large initial data. The method of our proof is a combination of a time semi-discretization, a special entropy symmetrizing the system, and compactness arguments.

2000 Mathematics Subject Classification: 35K55, 35D05, 92D25

Keywords: global existence, weak solution, time semi-discretization, viscous regularization, weighted Sobolev spaces

1 Introduction

Following Shigesada, Kawasaki and Teramoto [17], the time evolution of the population densities of two interacting species can be modeled by the system

∂tuj−div (∇((δj+δj1u1+δj2u2)uj) +τjuj∇U) = (αj−βj1u1−βj2u2)uj, u_j(0, x) =u_j0(x),

)

(1.1) where j = 1,2 and x∈Ω⊂Rⁿ,n= 1,2,3. The function uj =uj(t, x) ≥0 denotes the population density of a species j. The parameters δ_j and δ_ji describe diffusion phenomena: δ_j is the diffusion rate,δ_jj is a self-diffusion rate, andδ_ji forj6=iare the cross-diffusion rates. The parametersτ_j are related to population flows in direction to areas of better environmental quality, which is described by the environment potential U. The coefficient α_j is the intrinsic growth rate, and the parameters β_ji correspond to the inter-specific and intra-specific competition.

In case of Ω 6= Rⁿ, appropriate boundary conditions on u_j have to be added, for instance no-flux boundary conditionsJ_j·ν = 0, whereJ_j =∇((δ_j+δ_j1u₁+δ_j2u₂)u_j) +τ_ju_j∇U, andν is the normal vector on ∂Ω.

∗Fachbereich Mathematik und Statistik, P.O.Box D187, Universit¨at Konstanz, 78457 Konstanz, Germany, michael.dreher@uni-konstanz.de

1

(3)

The system (1.1) can be written in the form

∂_t u₁

u2

−div

A(u₁, u₂) ∇u₁

∇u2

+

τ₁u₁ 0 0 τ2u2

∇U

= f₁

f2

, (1.2)

where Ais the diffusion matrix, A=

δ1+ 2δ11u1+δ12u2 δ12u1

δ₂₁u₂ δ₂+δ₂₁u₁+ 2δ₂₂u₂

. (1.3)

Systems (1.2) with a matrixAas in (1.3) have numerous applications: we mention reaction-diffusion problems [3], where U is an electric potential; or the drift-diffusion equations as in the theory of semi-conductors [14]. And forδ_j = 0,δ₁₂=δ₂₁= 0, we have at hand a degenerate parabolic system as it appears in porous medium problems [10].

The matrixA may not be positive definite for large positive values ofu₁,u₂; and then the standard approach towards a prioriestimates of an energy ofL² type will not work. Additionally, maximum principles are not available because the two equations of (1.2) form a strongly coupled system.

We list some known results:

Steady state solutions in bounded domains Ω were studied, e.g., in [12], [13], [16]. Depending on the ranges of the diffusion, growth and competition parameters, one of the species may be extinct;

or there can be constant steady states; or non-constant steady states are possible and segregation of the species may happen. Numerical simulations of the steady state equations can be found in [8].

If δ12 = 0 or δ21 = 0, then the matrix A has a triangular form, and the general framework of [2]

gives the local well-posedness of initial-boundary value problems for (1.2). For results concerning global existence and global attractors of such weakly coupled systems, we refer to [6] or [15].

In [11], the existence and uniqueness of a local smooth non-negative solution was shown for a one- dimensional domain Ω, and δ_jj = 0, δ₁₂ = δ₂₁ = 1. Under the additional assumption δ₁ = δ₂, this smooth solution turns out to be global in time. For δjj = 0 and small initial data, the global existence and uniqueness of solutions in arbitrary dimensions was proved in [7]. For large self- diffusion coefficients, in the sense of 0 < δ₂₁<8δ₁₁, 0< δ₁₂ < δ₂₂, the existence and uniqueness of a non-negative strict solution for Ω⊂R² was demonstrated in [19]. Under the same assumption on theδ_ji, the global existence of a weak solution in the case of arbitrary dimensions was proved in [8].

The global existence of weak solutions without assumptions on the size of the coefficients, in spatial dimensions up to three, was shown in [5].

In contrast to the above mentioned results, which assumed a bounded domain Ω, the present paper deals with the population model (1.1) in the unbounded domain Ω =Rⁿ,n= 1,2,3.

To be specific, we list our assumptions: the coefficients are supposed to satisfy δ_j, δ_ji, β_ji>0,

αj ≥0, τ_j ∈R.







(1.4) We assume that the initial datau₁₀, u₂₀are positive functions onRⁿand belong to weighted Lebesgue spaces and Orlicz spaces:

u_j0(x)>0, x∈Rⁿ, j= 1,2, (1.5)

u_j0(x)hxi, u_j0(x) lnu_j0(x)∈L¹(Rⁿ), u_j0 ∈L²(Rⁿ), j = 1,2, (1.6)

(4)

where hxi= (1 +|x|²)^1/2.

The environment potential U = (t, x) is a function on (0,∞)×Rⁿ with

∇U ∈C([0,∞), L³(Rⁿ) +L^∞(Rⁿ)), (1.7)

∇U ∈C([0,∞), L²(Rⁿ)), (1.8)

△U ∈C([0,∞), L²(Rⁿ) +L^∞(Rⁿ)). (1.9)

Then a global in time weak solution exists:

Theorem 1.1. Suppose (1.4) and (1.7)–(1.9). Define the entropy density functional e, e(f) =fln(f)−f+ 1, f ≥0.

There is a weight function ̺ =̺(t, x) on (0,∞)×Rⁿ, depending only on the coefficients of (1.1), with the following property:

For each pair of initial functions with (1.5) and (1.6), there is a non-negative solution (u₁, u₂) belonging toL^∞_loc(R₊, L¹(Rⁿ))andL²_loc(R₊, H¹(Rⁿ)), which satisfies (1.1)in the distributional sense.

The entropy of this solution relative to the weight function ̺ exists:

E(t) =

2

X

j=1

Z

Rⁿ

e

uj(t, x)

̺(t, x)

̺(t, x) dx <∞, 0≤t <∞, and we have the a prioriestimate:

kEkL^∞(0,T)+X

j

√uj

2

L²((0,T),L²(Rⁿ,hxidx))+X

j

kujk²L²((0,T),H¹(Rⁿ)) (1.10)

+X

j

√u_j

2

L²((0,T),H¹(Rⁿ))+k√u₁u₂k²_L²_((0,T_),H¹_(Rⁿ₎₎

≤C(T), 0< T <∞.

We give some remarks on the strategy of the proof. Consider first the case of a bounded domain Ω, and assume that sufficiently regular positive solutions u_j exist. Multiplying the equations of (1.1) with lnuj, integrating over Ω, and performing suitable integrations by part, an estimate of the form (1.10) can be derived, where E is defined as in Theorem 1.1, but with ̺ ≡ 1. Of course, this derivation is only formal. To show the existence of non-negative solutions u_j, we introduce wj = lnuj, derive a system of differential equations for the wj, and seek a bound of wj(t,·) in the space H²(Ω). The continuous embedding H²(Ω) ⊂ L^∞(Ω) for n ≤ 3 will then show u_j > 0.

This approach can be made rigorous with a discretization of the time variable as in [9], and a subsequent spatial discretization with finite differences or a Galerkin scheme, and possibly a viscous regularization. This way, an approximate solution can be obtained. The convergence of this sequence of approximate solutions is shown by compactness arguments and the Lions-Aubin Lemma. The limit then is a weak solution of (1.1). However, we have to remark that the uniqueness of such weak solutions and their regularity are delicate questions, see [1].

This strategy will fail in case of Ω = Rⁿ: first, it is natural to assume that u_j(t, x) decays to zero for |x| → ∞, making the standard entropy infinite for all times. Second, since lnuj(t, x) = −∞

at infinity, the partial integrations have to be justified. And ultimately, the above compactness

(5)

arguments no longer hold. We overcome these difficulties by introducing the modified entropy from Theorem 1.1, which compares the functionu_j against an exponentially decaying weight function̺.

The time derivative of this weight function will then reinstate the needed compact embedding.

With minor modifications of the proof, we can also study the case of Ω being the exterior domain of an obstacle, with no-flux boundary conditions on ∂Ω. Finally, we remark that the machinery of our proof works also for the bounded domain case with no-flux boundary conditions. Then we can put̺≡1 andκ0 = 0 and consider the case of vanishing competition rates,βji= 0.

2 Proof of Theorem 1.1

2.1 Construction of an Entropy

For simplicity of notation, we scale the functions u₁ and u₂ by multiplications with appropriate constants in such a way that the constants δ₁₂ and δ₂₁ become both equal to one.

Then we choose a positive constant κ0 with the property that 8κ²₀≤β₁₂+β₂₁, 4δ_jκ²₀≤1,

δ_jj+1

2τ_j²

κ²₀ ≤ 1

16β_jj, j= 1,2. (2.1) Next, we determine a positive number λin such a way that the function κ=κ(x) =hλxi satisfies

k∇xκ(x)kL^∞(Rⁿ)≤κ₀.

Then we set µ(t) = _1+t¹ for 0≤t <∞, and define weight functions

σ =σ(t, x) =−µ(t)κ(x), ̺(t, x) = exp(σ(t, x)), (t, x)∈[0,∞)×Rⁿ. By the choice of the parameters, we have

0< ̺(t, x)<1, |∇^xσ(t, x)| ≤µ(t)κ0 ≤κ0, (t, x)∈[0,∞)×Rⁿ.

For positive real numbers ̺ and v, we put Φ̺(v) =vln

v

̺

−v+̺.

Observe that Φ_̺(v)≥0 for̺, v >0; and Φ_̺(·) has a minimum atv=̺, taking the value zero there.

For two functions u1, u2, taking positive values on [0,∞)×Rⁿ, we write our generalized entropy functional E as

E(t) =

2

X

j=1

Z

Rⁿ

Φ_̺(t,x)(u_j(t, x)) dx, 0≤t <∞. (2.2)

(6)

2.2 The Semi-discretization Scheme

We define a small time step-size h >0 and sett_k =kh, fork= 0,1, . . .. Thinking ofu^k_j =u^k_j(x) as an approximation of u_j(t_k, x), we wish to solve the system

u^k_j −u^k−1_j

h −div

∇(δ_j +δ_j1u^k₁ +δ_j2u^k₂)u^k_j +τ_ju^k_j∇U(t_k,·)

= (α_j−β_j1u^k₁−β_j2u^k₂)u^k_j, forj= 1,2 andk∈N₊. However, it seems hard to prove the existence of a solution to this system.

Instead, we perform an exponential change of the dependent variables and insert a higher order elliptic regularization:

u^k,ε_j −u^k−1,ε_j

h +ε(△⁴+hxi⁸)w^k,ε_j −div

∇(δ_j+δ_j1u^k,ε₁ +δ_j2u^k,ε₂ )u^k,ε_j +τ_ju^k,ε_j ∇U(t_k,·) (2.3)

= (α_j−β_j1u^k,ε₁ −β_j2u^k,ε₂ )u^k,ε_j , j= 1,2, k∈N₊, w^k,ε_j (x) = ln u^k,ε_j (x)

̺(t_k, x)

!

, (2.4)

whereε >0 is a small parameter. Of course, the transformation (2.4) is only valid ifu^k,ε_j (x)>0 for all x∈Rⁿ.

For fixed h >0 and ε >0, we define the discrete entropy E^k,ε =

2

X

j=1

Z

Rⁿ

Φ_̺(t_k_,x)(u^k,ε_j (x)) dx, t_k=kh, k∈N₀.

Lemma 2.1. Suppose ε > 0 and u^k−1,ε_j ∈ L²(Rⁿ) with E^k−1,ε < ∞, and u^k−1,ε_j (x) > 0 almost everywhere in Rⁿ, j= 1,2.

Then the problem (2.3)–(2.4) has a solution w^k,ε₁ , w₂^k,ε ∈H⁴(Rⁿ) with hxi⁴w^k,ε_j (x) ∈L²(Rⁿ). The functions u^k,ε_j belong to H⁴(Rⁿ) and take only positive values on Rⁿ. The functions u^k,ε_j , w^k,ε_j solve (2.3)in the distributional sense:

1 h

Z

Rⁿ

(u^k,ε_j −u^k−1,ε_j )ψdx+ε Z

Rⁿ

w_j^k,ε(△⁴+hxi⁸)ψdx (2.5)

− Z

Rⁿ

((δj+δj1u^k,ε₁ +δj2u^k,ε₂ )u^k,ε_j )△ψdx+τj

Z

Rⁿ

u^k,ε_j (∇U(t_k,·))∇ψdx

= Z

Rⁿ

(αj −βj1u^k,ε₁ −βj2u^k,ε₂ )u^k,ε_j ψdx, j= 1,2, ψ∈C₀^∞(Rⁿ).

Furthermore, if h >0 is small enough, then we have the a prioriestimate E^k,ε−E^k−1,ε+εh

2

X

j=1

△²w^k,ε_j

2

L²(Rⁿ)+

hxi⁴w_j^k,ε

2 L²(Rⁿ)

(2.6)

−hµ^′(kh) 2

2

X

j=1

Z

Rⁿ

κ(x)u^k−1,ε_j (x) dx+1 2h

2

X

j=1

δjj

∇u^k,ε_j

2

L²(Rⁿ)+ h 16

2

X

j=1

βjj

u^k,ε_j

2 L²(Rⁿ)

+h



3

2

X

j=1

δ_j ∇

q u^k,ε_j

2 L²(Rⁿ)

+ 2 ∇

q

u^k,ε₁ u^k,ε₂

2 L²(Rⁿ)

+ 2κ²₀

q

u^k,ε₁ u^k,ε₂

2 L²(Rⁿ)





(7)

≤3hα₀E^k,ε+hC_β ̺^k

2

L²(Rⁿ)+ 3hα₀ ̺^k

L¹(Rⁿ)

+ 2 Z

Rⁿ

̺^k−̺^k−1

dx+hC_τ,δk∇U(t_k,·)k²L²(Rⁿ), where α₀ = max(α₁, α₂) and C_β, C_τ,δ are constants defined below.

Proof. We exploit the Leray-Schauder fixed-point principle. For a parameter ζ ∈ [0,1] and given functionsu^k,ε_j andw^k,ε_j satisfying (2.4) withw_j^k,ε∈H²(Rⁿ), and u^k−1,ε_j ∈L²(Rⁿ) withu^k−1,ε_j (x)>0 on Rⁿ, we look for functionsW_j^k,ε as solutions to

ε(△⁴+hxi⁸)W_j^k,ε =−ζu^k,ε_j −u^k−1,ε_j

h (2.7)

+ζdiv

∇(δ_j+δ_j1u^k,ε₁ +δ_j2u^k,ε₂ )u^k,ε_j +τ_ju^k,ε_j ∇U(t_k,·) +ζ(α_j−β_j1u^k,ε₁ −β_j2u^k,ε₂ )u^k,ε_j , j= 1,2.

Fromw_j^k,ε∈H²(Rⁿ) andn≤3 we getw^k,ε_j ∈L^∞(Rⁿ) and∇w^k,ε_j ∈L⁶(Rⁿ)∩L²(Rⁿ). Then it follows that exp(w^k,ε_j ) ∈L^∞(Rⁿ) and ∇exp(w^k,ε_j )∈L⁶(Rⁿ),∇²exp(w^k,ε_j )∈L²(Rⁿ). As a consequence, the function u^k,ε_j , defined via (2.4), belongs to H²(Rⁿ). Then the right-hand side of (2.7) belongs to L²(Rⁿ), where we have used (1.7) and (1.9).

By the Lax-Milgram theorem, this problem has a unique solution (W₁^k,ε, W₂^k,ε) in the space X⁴ =n

v∈H⁴(Rⁿ) : hxi⁴v(x)∈L²(Rⁿ)o .

For 0 ≤ ζ ≤ 1, we define a mapping, with parameter ζ, S = S(w^k,ε₁ , w₂^k,ε;ζ) = (W₁^k,ε, W₂^k,ε) from H²(Rⁿ) into X⁴, which is compactly embedded intoH²(Rⁿ), as can be seen from a variant of the Arzela-Ascoli theorem, compare also [4].

Clearly, the operatorS(·,·; 0) has a unique fixed point (W₁^k,ε, W₂^k,ε) = (0,0).

Next, we show that a fixed point (w^k,ε₁ , w^k,ε₂ ) to S(·,·;ζ) satisfies an a priori estimate inX⁴, inde- pendent ofζ. The Leray-Schauder fixed-point principle will then guarantee the existence of at least one fixed point of S(·,·; 1).

Let (w₁^k,ε, w^k,ε₂ ) = (W₁^k,ε, W₂^k,ε)∈X⁴ be a solution to (2.7). Multiply (2.7) with w^k,ε_j , integrate over Rⁿ, and sum overj= 1,2:

ζX

j

Z

Rⁿ

(u^k,ε_j −u^k−1,ε_j )w_j^k,εdx+εhX

j

△²w_j^k,ε

2 L²(Rⁿ)+

hxi⁴w^k,ε_j

2 L²(Rⁿ)

(2.8)

−ζhX

j

Z

Rⁿ

w_j^k,εdiv

∇(δj+δj1u^k,ε₁ +δj2u^k,ε₂ )u^k,ε_j +τju^k,ε_j ∇U(tk,·) dx

=ζhX

j

Z

Rⁿ

w_j^k,ε(α_j−β_j1u^k,ε₁ −β_j2u^k,ε₂ )u^k,ε_j dx.

(8)

For simplicity of notation, we set ̺^k(x) = ̺(t_k, x), σ^k(x) = σ(t_k, x), µ^k = µ(t_k), µ^k₁ = µ^′(t_k), U^k(x) =U(t_k, x). Moreover, we fix

α₀= max(α₁, α₂), β₀ = max(β₁₂, β₂₁).

We estimate the first term:

E^k,ε−E^k−1,ε =X

j

Z

Rⁿ

u^k,ε_j w_j^k,ε−u^k,ε_j +̺^k

−

u^k−1,ε_j w_j^k−1,ε−u^k−1,ε_j +̺^k−1 dx

=X

j

Z

Rⁿ

(u^k,ε_j −u^k−1,ε_j )w^k,ε_j dx

+X

j

Z

Rⁿ

u^k−1,ε_j

lnu^k,ε_j −lnu^k−1,ε_j

+ (u^k−1,ε_j −u^k,ε_j ) dx

+X

j

Z

Rⁿ

u^k−1,ε_j (σ^k−1−σ^k) +̺^k−̺^k−1 dx

≤X

j

Z

Rⁿ

(u^k,ε_j −u^k−1,ε_j )w^k,ε_j dx+X

j

Z

Rⁿ

u^k−1,ε_j ln u^k,ε_j u^k−1,ε_j

!

+ 1− u^k,ε_j u^k−1,ε_j

! dx

+X

j

Z

Rⁿ

µ^k₁hκ(x)u^k−1,ε_j +̺^k−̺^k−1 dx.

In the second integral of (2.8), we perform integration by part and exploit Young’s inequality:

−X

j

Z

Rⁿ

w_j^k,εdiv

∇(δ_j+δ_j1u^k,ε₁ +δ_j2u^k,ε₂ )u^k,ε_j +τ_ju^k,ε_j ∇U^k dx

=X

j n

X

l=1

Z

Rⁿ

δ_j(∂_lw^k,ε_j )(∂_lu^k,ε_j ) dx

+X

j n

X

l=1

Z

Rⁿ

(∂_lw_j^k,ε)

(∂_l(δ_j1u^k,ε₁ +δ_j2u^k,ε₂ )u^k,ε_j ) +τ_ju^k,ε_j ∂_lU^k dx

=X

j n

X

l=1

Z

Rⁿ

δ_j ∂_lu^k,ε_j

u^k,ε_j −∂_lσ^k

!

(∂_lu^k,ε_j ) dx+X

j n

X

l=1

Z

Rⁿ

∂_lu^k,ε_j

!

∂_l(u^k,ε₁ u^k,ε₂ ) dx

+X

j n

X

l=1

Z

Rⁿ

∂_lu^k,ε_j

!

δjj∂_l(u^k,ε_j )²+τju^k,ε_j ∂_lU^k dx

= 4X

j

δ_j ∇

q u^k,ε_j

2 L²(Rⁿ)

−2X

j

δ_j Z

Rⁿ

q u^k,ε_j

∇σ^k

∇ q

u^k,ε_j

dx + 4

∇

q

u^k,ε₁ u^k,ε₂

2 L²(Rⁿ)

−4 Z

Rⁿ

q

u^k,ε₁ u^k,ε₂

∇σ^k

∇ q

u^k,ε₁ u^k,ε₂

dx + 2X

j

δjj

∇u^k,ε_j

2

L²(Rⁿ)−2X

j

δjj

Z

Rⁿ

u^k,ε_j

∇σ^k ∇u^k,ε_j dx

+X

j

τj

Z

Rⁿ

(∇u^k,ε_j )

∇U^k

dx−X

j

τj

Z

Rⁿ

u^k,ε_j

∇σ^k ∇U^k dx

(9)

≥4X

j

δ_j ∇

q u^k,ε_j

2 L²(Rⁿ)

+ 4 ∇

q

u^k,ε₁ u^k,ε₂

2 L²(Rⁿ)

+ 2X

j

δ_jj ∇u^k,ε_j

2 L²(Rⁿ)

−X

j

δ_j Z

Rⁿ

u^k,ε_j ∇σ^k

2 dx−X

j

δ_j ∇

q u^k,ε_j

2 L²(Rⁿ)

−2 Z

Rⁿ

u^k,ε₁ u^k,ε₂ ∇σ^k

2 dx

−2 ∇

q

u^k,ε₁ u^k,ε₂

2 L²(Rⁿ)

−X

j

δ_jj Z

Rⁿ

(u^k,ε_j )² ∇σ^k

2 dx−X

j

δ_jj ∇u^k,ε_j

2 L²(Rⁿ)

−1 2

X

j

δjj

∇u^k,ε_j

2

L²(Rⁿ)−



 X

j

τ_j² 2δ_jj + 1



 ∇U^k

2

L²(Rⁿ)−X

j

τ_j² 2

Z

Rⁿ

(u^k,ε_j )² ∇σ^k

2 dx

≥3X

j

δj

∇

q u^k,ε_j

2 L²(Rⁿ)

+ 2 ∇

q

u^k,ε₁ u^k,ε₂

2 L²(Rⁿ)

+1 2

X

j

δjj

∇u^k,ε_j

2 L²(Rⁿ)

−X

j

δ_j(µ^kκ₀)² Z

Rⁿ

u^k,ε_j dx−2(µ^kκ₀)² Z

Rⁿ

u^k,ε₁ u^k,ε₂ dx

−X

j

δjj+1

2τ_j²

(µ^kκ0)² u^k,ε_j

2 L²(Rⁿ)−



 X

j

τ_j² 2δ_jj + 1



 ∇U^k

2 L²(Rⁿ). Finally, we consider the right-hand side of (2.8):

X

j

Z

Rⁿ

w^k,ε_j

α_j −β_j1u^k,ε₁ −β_j2u^k,ε₂ u^k,ε_j dx

=X

j

αj

Z

Rⁿ

Φ_̺k(x)(u^k,ε_j (x)) dx+X

j

αj

Z

Rⁿ

u^k,ε_j −̺^kdx−

2

X

j=1 2

X

i=1

βji

Z

Rⁿ

w_j^k,εu^k,ε_i u^k,ε_j dx

≤α₀E^k,ε+α₀X

j

Z

Rⁿ

u^k,ε_j dx

−β12

Z

Rⁿ

̺^ku^k,ε₂ u^k,ε₁

̺^k ln u^k,ε₁

̺^k

! + 1

!

dx−β21

Z

Rⁿ

̺^ku^k,ε₁ u^k,ε₂

̺^k ln u^k,ε₂

̺^k

! + 1

! dx

−1 2

X

j

βjj

Z

Rⁿ

(̺^k)²



 u^k,ε_j

̺^k

!2

ln



 u^k,ε_j

̺^k

!2

+ 1



 dx +

Z

Rⁿ

(β₁₂u^k,ε₂ +β₂₁u^k,ε₁ )̺^kdx+1 2

X

j

β_jj ̺^k

2

L²(Rⁿ)−(α₁+α₂) ̺^k

L¹(Rⁿ). From the elementary inequality z≤2zlnz−2z+ 4 forz≥0 we then obtain

α₀X

j

Z

Rⁿ

u^k,ε_j dx≤2α₀E^k,ε+ 4α₀ Z

Rⁿ

̺^kdx.

And by Young’s inequality, with a constant C_β = max(2β₀²/β_jj) + (β₁₁+β₂₂)/2, Z

Rⁿ

(β₁₂u^k,ε₂ +β₂₁u^k,ε₁ )̺^kdx+1 2

X

j

β_jj ̺^k

2

L²(Rⁿ) ≤ 1 8

X

j

β_jj u^k,ε_j

2

L²(Rⁿ)+C_β ̺^k

2 L²(Rⁿ).

(10)

We define an auxiliary function L = L(z) = zlnz + 1 for z > 0 and summarize the estimates obtained so far:

ζ

E^k,ε−E^k−1,ε

−ζX

j

Z

Rⁿ

u^k−1,ε_j ln u^k,ε_j u^k−1,ε_j

!

+ 1− u^k,ε_j u^k−1,ε_j

! dx +ζX

j

Z

Rⁿ

−µ^k₁hκ(x)u^k−1,ε_j −̺^k+̺^k−1 dx +εhX

j

△²w_j^k,ε

2

L²(Rⁿ)+

hxi⁴w_j^k,ε

2 L²(Rⁿ)

+ζh



3X

j

δ_j ∇

q u^k,ε_j

2 L²(Rⁿ)

+ 2 ∇

q

u^k,ε₁ u^k,ε₂

2 L²(Rⁿ)

+1 2

X

j

δ_jj ∇u^k,ε_j

2 L²(Rⁿ)





+1 2ζhX

j

β_jj Z

Rⁿ

(̺^k)²L



 u^k,ε_j

̺^k

!2

 dx+ 2ζhκ²₀ Z

Rⁿ

≤ζh



 X

j

δ_j(µ^kκ₀)² Z

Rⁿ

u^k,ε_j dx+ 4κ²₀ Z

Rⁿ



+ 3ζhα₀E^k,ε

−ζh Z

Rⁿ

̺^k β₁₂u^k,ε₂ L u^k,ε₁

̺^k

!

+β₂₁u^k,ε₁ L u^k,ε₂

̺^k

!!

dx +ζhX

j

δ_jj+1

2τ_j²

κ²₀+β_jj 8

u^k,ε_j

2

L²(Rⁿ)+ζhC_β ̺^k

2

L²(Rⁿ)+ 3ζhα₀ ̺^k

L¹(Rⁿ)

+ζhC_τ,δ ∇U^k

2

L²(Rⁿ), C_τ,δ =X

j

τ_j² 2δjj

+ 1.

Noting that L(z)≥ ¹₂z andµ^k≤1, we re-order the terms:

ζ

E^k,ε−E^k−1,ε

+εhX

j

△²w^k,ε_j

2

L²(Rⁿ)+

hxi⁴w^k,ε_j

2 L²(Rⁿ)

+ζX

j

Z

Rⁿ

u^k−1,ε_j −µ^k₁hκ(x)−ln u^k,ε_j u^k−1,ε_j

!

−1 +

1−δ_j(µ^kκ₀)²h u^k,ε_j u^k−1,ε_j

! dx

+ζh



3X

j

δ_j ∇

q u^k,ε_j

2 L²(Rⁿ)

+ 2 ∇

q

u^k,ε₁ u^k,ε₂

2 L²(Rⁿ)

+ 2κ²₀

q

u^k,ε₁ u^k,ε₂

2 L²(Rⁿ)





+ζhX

j

δ_jj 2

∇u^k,ε_j

2

L²(Rⁿ)+ζhX

j

1 8β_jj−

δ_jj+1

2τ_j²

κ²₀

u^k,ε_j

2 L²(Rⁿ)

≤ζh

4κ²₀−1

2(β₁₂+β₂₁) Z

Rⁿ

u^k,ε₁ u^k,ε₂ dx+ 3ζhα₀E^k,ε +ζhC_β

̺^k

2

L²(Rⁿ)+ 2ζ Z

Rⁿ

̺^k−̺^k−1

dx+ 3ζhα₀ ̺^k

L¹(Rⁿ)+ζhC_τ,δ ∇U^k

2 L²(Rⁿ).

(11)

Observe that δjκ²₀ ≤ ¹₄ and (µ^k)² =−µ^k₁, by (2.1). Using 8κ²₀ ≤β12+β21 from (2.1), we can drop the first integral on the right-hand side, and deduce that

ζ

E^k,ε−E^k−1,ε

+εhX

j

△²w^k,ε_j

2

L²(Rⁿ)+

hxi⁴w^k,ε_j

2 L²(Rⁿ)

+ζX

j

Z

Rⁿ

u^k−1,ε_j −µ^k₁hκ(x)−ln u^k,ε_j u^k−1,ε_j

!

−1 +

1 +µ^k₁h 4

u^k,ε_j u^k−1,ε_j

! dx

+ζh



3X

j

δ_j ∇

q u^k,ε_j

2 L²(Rⁿ)

+ 2 ∇

q

u^k,ε₁ u^k,ε₂

2 L²(Rⁿ)

+ 2κ²₀

q

u^k,ε₁ u^k,ε₂

2 L²(Rⁿ)





+1 2ζhX

j

δ_jj ∇u^k,ε_j

2

L²(Rⁿ)+ 1

16ζhX

j

β_jj u^k,ε_j

2 L²(Rⁿ)

≤3ζhα₀E^k,ε+ζhC_β ̺^k

2

L²(Rⁿ)+ 3ζhα₀ ̺^k

L¹(Rⁿ)

+ 2ζ Z

Rⁿ

̺^k−̺^k−1

dx+ζhC_τ,δ ∇U^k

2 L²(Rⁿ).

The a prioriestimate ofw_j^k,ε in the space X⁴ is found if we chooseh so small that 3hα₀ <1 and if we can show that the integral on the left-hand side is positive.

From ln(z)≤z−1 for z >0, we deduce that

−µ^k₁hκ(x)−ln u^k,ε_j u^k−1,ε_j

!

−1 +

1 +µ^k₁h 4

u^k,ε_j

u^k−1,ε_j ≥ −µ^k₁h κ(x)− 1 4

! ,

and this is greater than −¹₂µ^k₁hκ(x) for those xwith u^k,ε_j (x)≤2u^k−1,ε_j (x), sinceκ(x)≥1.

Now we study those xwithu^k,ε_j (x)>2u^k−1,ε_j (x). We have ln(z)≤ηz−1 with η= ¹₂(ln(2) + 1)<1 forz≥2, from which it follows that

−µ^k₁hκ(x)−ln u^k,ε_j u^k−1,ε_j

!

−1 +

1 +µ^k₁h 4

u^k,ε_j

u^k−1,ε_j ≥ −µ^k₁hκ(x) +

1 +µ^k₁h 4 −η

≥ −µ^k₁hκ(x),

under the assumption−µ^k₁h≤4(1−η).

We end up with thea priori estimate ζ

E^k,ε−E^k−1,ε

+εhX

j

△²w^k,ε_j

2

L²(Rⁿ)+

hxi⁴w^k,ε_j

2 L²(Rⁿ)

−ζhµ^k₁ 2

X

j

Z

Rⁿ

κ(x)u^k−1,ε_j (x) dx+ζhX

j

δjj

2 ∇u^k,ε_j

2

L²(Rⁿ)+ζhX

j

βjj

16 u^k,ε_j

2 L²(Rⁿ)

+ζh



3X

j

δ_j ∇

q u^k,ε_j

2 L²(Rⁿ)

+ 2 ∇

q

u^k,ε₁ u^k,ε₂

2 L²(Rⁿ)

+ 2κ²₀

q

u^k,ε₁ u^k,ε₂

2 L²(Rⁿ)



