Vortex formation for a non-local interaction model with Newtonian repulsion and superlinear mobility

J.A. Carrillo; D. Gómez-Castro; J.L. Vázquez

doi:10.1515/anona-2021-0231

Open Access Published by De Gruyter February 25, 2022

Vortex formation for a non-local interaction model with Newtonian repulsion and superlinear mobility

J.A. Carrillo , D. Gómez-Castro and J.L. Vázquez

From the journal Advances in Nonlinear Analysis

https://doi.org/10.1515/anona-2021-0231

Abstract

We consider density solutions for gradient flow equations of the form u_t = ∇ · (γ(u)∇ N(u)), where N is the Newtonian repulsive potential in the whole space ℝ^d with the nonlinear convex mobility γ(u) = u^α, and α > 1. We show that solutions corresponding to compactly supported initial data remain compactly supported for all times leading to moving free boundaries as in the linear mobility case γ(u) = u. For linear mobility it was shown that there is a special solution in the form of a disk vortex of constant intensity in space u = c₁t⁻¹ supported in a ball that spreads in time like c₂t^1/d, thus showing a discontinuous leading front or shock. Our present results are in sharp contrast with the case of concave mobilities of the form γ(u) = u^α, with 0 < α < 1 studied in [10]. There, we developed a well-posedness theory of viscosity solutions that are positive everywhere and moreover display a fat tail at infinity. Here, we also develop a well-posedness theory of viscosity solutions that in the radial case leads to a very detailed analysis allowing us to show a waiting time phenomena. This is a typical behaviour for nonlinear degenerate diffusion equations such as the porous medium equation. We will also construct explicit self-similar solutions exhibiting similar vortex-like behaviour characterizing the long time asymptotics of general radial solutions under certain assumptions. Convergent numerical schemes based on the viscosity solution theory are proposed analysing their rate of convergence. We complement our analytical results with numerical simulations illustrating the proven results and showcasing some open problems.

Keywords: nonlinear mobility equations; conservation laws; viscosity solutions; shock conditions

MSC 2010: 35L65; 35D40; 65M25

1 Introduction

We are interested in the family of equations of the form

{ut=∇⋅(γ(u)∇v)(0,+∞)×ℝd,−Δv=u(0,+∞)×ℝd,u=u0t=0,

where the function γ(u) is called the mobility. They all correspond to gradient flows with nonlinear mobility of the Newtonian repulsive interaction potential in dimension d ≥ 1

ℱ(u)=12∫ℝdN(u)u dx ,

with N(u) the Newtonian repulsive potential [8], as they can be written in the form

(1.1) {ut+∇⋅(γ(u)w)=0(0,+∞)×ℝd,w=−∇δℱδu(0,+∞)×ℝd.

We will consider nonnegative data and solutions. The linear case γ(u) = u is well-known in the literature as a model for wave propagation in superconductivity or superfluidity, cf. Lin and Zhang [16], Ambrosio, Mainini, and Serfaty [4, 5], Bertozzi, Laurent, and Léger [7], Serfaty and Vazquez [18]. In that case the theory leads to uniqueness of bounded weak solutions having the property of compact space support, and in particular there is a special solution in the form of a disk vortex of constant intensity in space u = c₁t⁻¹ supported in a ball that spreads in time like c₂t^1/d, thus showing a discontinuous leading front or shock. This vortex is the generic attractor for a wide class of solutions.

We want to concentrate on models with nonlinear mobility of power-like type γ(u) = u^α, α > 0. The sublinear concave 0 < α < 1 range was studied in our previous paper [10]. For nonnegative data the study provides a theory of viscosity solutions for radially symmetric initial data that are positive everywhere, and moreover display a fat tail at infinity. In particular the standard vortex of the linear mobility transforms into an explicit selfsimilar solution that reminds of the Barenblatt solution for the fast diffusion equation. A very detailed analysis is done for radially symmetric data and solutions via the corresponding mass function that satisfies a first-order Hamilton-Jacobi equation.

The present paper contains the rest of the analysis of power-like mobility for convex superlinear cases when γ(u) = u^α, and α > 1. Again, we perform a fine analysis in the case of radially symmetric solutions by means of the study of the corresponding mass function. The theory of viscosity solutions for the mass function still applies. As for qualitative properties, let us stress that in this superlinear parameter range α > 1 solutions recover the finite propagation property and the existence of discontinuity fronts (shocks). We analyse in detail how the stable asymptotic solution goes from the fat tail profile of the sublinear case α < 1 to the shock profile of the range α > 1 when passing through the critical value α = 1.

Another important aspect of the well-posedness theory that we develop for viscosity solutions with radially symmetric initial data, is that the classical approach based on optimal transport theory for equations of the form (1.1) developed in [3, 8, 13] fail for convex superlinear mobilities as described in [8] since the natural associated distance is not well-defined [13]. Therefore, our present results are the first well-posedness results for gradient flows with convex superlinear power-law mobilities, even if only for radially symmetric initial data. Finally, let us mention that we still lack of a well-posedness theory for gradient flow equations of the form (1.1) with convex superlinear mobilities for general initial data, possibly showing that the vortex-like solutions are generic attractors of the flow.

We also highlight how different convex superlinear mobilities are with respect to the linear mobility case by showing the property of an initial waiting time for the spread of the support for radial solutions that we are able to characterize. Indeed, let u₀ be radial and supported in a ball: supp u0=B¯R . We prove that there is finite waiting time at r = R if and only if

(1.2) lim supr→R− (R−r)1−αα∫r<|x|<Ru0(x) dx=C<+∞.

The waiting time phenomenon is typical of slow diffusion equations like the Porous Medium Equation [20] or the p-Laplacian equation. In our class of equations it does not occur for the whole range 0 < α ≤ 1. We are able to estimate the waiting time in terms of the limit constant C in (1.2).

We combine the theory with the computational aspect: we identify suitable numerical methods and perform a detailed numerical analysis. Indeed, we construct numerical finite-difference convergent schemes and prove convergence to the actual viscosity solution for radially symmetric solutions based on the mass equation. By taking advantage of the connection to nonlinear Hamilton-Jacobi equations, we obtain monotone numerical schemes showing their convergence to the viscosity solutions of the problem with a uniform rate of convergence, see Theorem 5.6 in constrast with the case of concave sublinear mobilities in [10].

The paper is structured as follows. We start by constructing some explicit solutions and developing the theory for radially symmetric initial data by using the mass equation in Sections 2 and 3 respectively. We construct the general viscosity solution theory for radially symmetric initial data in Section 4 showing the most striking phenomena for convex superlinear mobilities: compactly supported free boundaries determined by sharp fronts and the waiting time phenomena. Section 5 is devoted to show the convergence of monotone schemes for the developed viscosity solution theory with an explicit convergence rate. The numerical schemes constructed illustrate the sharpness of the waiting time phenomena result, and allow us to showcase interesting open problems in Section 6. A selection of figures illustrates a number of salient phenomena. We provide videos for some interesting situations as supplementary material in [1].

2 Explicit solutions

The aim of this section is to find some important explicit solutions of

(P) {ut=∇⋅(uα∇v)(0,+∞)×ℝd−Δv=u(0,+∞)×ℝdu=u0t=0

Notice that, as in [10], we still have that, for C ≥ 0

(2.1) u¯(t)=(C+αt)−1α

is a solution of the PDE. The repulsion potential v diverges quadratically at infinity. For C = 0 we recover the Friendly Giant solution, introduced in [10].

2.1 Self-similar solution

The algebraic calculations developed in [10] still work, we get self-similar solutions of the form

U(t,x)=t−1αF(|x|t−1αd),

with the self-similar profile

(2.2) F(|y|)=(α+(ωd|y|dα)−αα−1)−1α.

Let us remark that, for α > 1, we have F(0) = 0 and F(+∞)=α−1α (see Figure 1 for a sketch of the self-similar profiles depending on α). This is different to the case 0 < α < 1, where F(0) is a positive constant and F decays at infinity. In our present case α > 1, the self-similar solutions have infinite mass, whereas for 0 < α < 1 the self-similar solutions have finite mass.

Fig. 1

Self-similar profiles when d = 1

For 0 < α < 1 these self-similar give the typical asymptotic behaviour as t → +∞. For α > 1 we will show this is no longer the case, for finite mass solutions.

2.2 Vortices

The vortex solutions defined as

(2.3) u(t,x)={(u0−α+αt)−1αωd|x|d<S(t)=M/(u0−α+αt)−1α0otherwise

are local weak solutions of (P).

Remark 2.1

Notice that, for t→−u0−α/α , the vortex collapses to the Dirac delta of mass M.

This solution was also constructed by characteristics and the Rankine-Hugoniot condition in [10, Section 5.2]. However, in that case the Lax-Oleinik condition of incoming characteristics failed. In our present setting for α > 1, this shock-type solutions are entropic, and we will prove that they are indeed viscosity solutions of the mass equation (3.1). We will prove that, for α > 1, they now have significant relevance. In particular, they describe the asymptotic behaviour as t → +∞. Notice that the the radius of the support S(t) of this kind of solutions solves an equation of type

(2.4) dSdt=M(u0−α+αt)−1+1α.

There are also complementary vortex solutions:

(2.5) u(t,x)={0|x|d<a(c0−α+αt)−1α|x|d>a

which are stationary (and solve the mass problem (3.1) by characteristics). This type of solution belongs to a theory of solutions in L^∞, but not L¹.

3 Mass of radial solutions

In [10] we proved that the mass of a radial solution

m(t,r)=∫Bru(t,x) dx

written in volume coordinates ρ = dω_dr^d is a solution of the equation

(3.1) mt+m(mρ)α=0.

3.1 Characteristics for the mass equation

In [10] we computed the generalised characteristics for the mass equation in the case of sublinear mobility, α < 1. The algebra for characteristics still works

(3.2) ρ=ρ0+αm(ρ0)u0(ρ0)α−1t ,

and the solutions behave like

(3.3) u(t,ρ)=(u0(ρ0)−α+αt)−1α ,

and

m(t,ρ)=m0(ρ0)(1+αu0(ρ0)αt)1−1α .

Remark 3.1

Notice that the the generalised characteristics are not the level sets of m.

These solutions are well defined, for a given u₀, as long as the characteristics do not cross.

Proposition 3.2

Let u₀ be non-decreasing and C¹. Then, there is a classical global solution of the mass equation, given by characteristics.

Proof

Let P_t(ρ₀) = ρ₀ + αm(ρ₀)u₀(ρ₀)^α−1t. Clearly P_t(0) = 0. If u₀ is non-decreasing, dPtdρ0≥1 , and hence it is invertible. We construct

(3.4) u(t,ρ)={((u0(Pt−1(ρ)))−α+αt)−1αif u0(Pt−1(ρ))≠0,0if u0(Pt−1(ρ))=0.

It is immediate to see that u is continuous and C¹.

For 0 < α < 1, in [10] we developed a theory of classical solutions for non-increasing initial data. In that case, rarefaction fan tails appeared, which gave rise to classical solutions of the mass equation. In our present case α > 1, it seems that the good data is radially non-decreasing, but this is not possible in an L¹ ∩ L^∞ theory, unless a jump is introduced.

3.2 The Rankine-Hugoniot condition

We will prove in Section 4.3 that solutions with a jump, given by a Rankine-Hugoniot condition, will be the correct “stable” solutions. As in [10], shocks (i.e. discontinuities) propagate following a Rankine-Hugoniot condition. If S(t) is the position of the shock, we write the continuity of mass condition

m(t,S(t)−)=m(t,S(t)+).

Taking a derivative and applying the equation (3.1) we have that

(3.5) dSdt(t)=m(t,S(t))u(t,S(t)+)α−u(t,S(t)−)αu(t,S(t)+)−u(t,S(t)−).

Remark 3.3

Notice that, in the case of the vortex the Rankine-Hugoniot condition determines precisely the support. In particular, we have u(t, S(t)⁺) = 0 and u(t,S(t)−)=(u0−α+αt)−1α so (3.5) is precisely (2.4). In fact, the vortex is simply a cut-off of the Friendly Giant with a free-boundary determined by the Rankine-Hugoniot condition.

3.3 Local existence of solutions by characteristics

Theorem 3.4

Let 0 ≤ u₀ ∈ L¹(ℝⁿ) ∩ L^∞(ℝⁿ) be radial and such that u0α−1 is Lipschitz. Then, there exists a small time T > 0 and a classical solution of the mass equation given by characteristics defined for t ∈ [0, T].

Proof

The solution given by (3.2) and (3.3) is well defined as long as the characteristics cover the whole space, and do not cross. This is equivalent to P_t(ρ₀) = ρ₀ + αm(ρ₀)u₀(ρ₀)^α−1t being a bijection [0, +∞) → [0, +∞). Again, we construct (3.4). Since P_t(0) = 0, it suffices to prove that dPtdρ0≥c0>0 . We take the derivative explicitly and find that

dPtdρ0=1+αt(dm0dρ0u0α−1+m0ddρ0(u0α−1))=1+αt(u0α+m0ddρ0(u0α−1)).

Due to the hypothesis

L:=supρ0≥0|u0α+m0ddρ0(u0α−1)|<∞,

and we have that dPtdρ0≥1−α tL which is strictly positive if t ≤ T < 1/(αL). Since u0α−1 is Lipschitz and bounded, then 1−α TL≤dPtdρ0≤C is Lipschitz. Hence, Pt−1 is Lipschitz in ρ₀. Also, it is easy to see that Pt−1 is continuous in t. Since it is immediate to check that u is continuous by composition, we have that m is of class C¹, and the proof is complete.

Corollary 3.5

(Waiting time). Let 0 ≤ u₀ ∈ L¹(ℝⁿ) ∩ L^∞(ℝⁿ) be radial and such that u0α−1 is Lipschitz. Then, there is a short time T > 0 such that, if supp u0⊂BR¯ , then any classical solution of the mass equation satisfies supp u(t,⋅)⊂BR¯ .

Proof

Notice that, if u₀ is compactly supported, then outside the support the characteristics are given by P_t(ρ₀) = ρ₀. As long as the solution is given by characteristics, if supp u0⊂BR¯ for ρ > ω_dR^d, then u(t, ρ) = 0.

Remark 3.6

This effect of preservation of the support for a finite time is known as waiting time. In Section 3.4 we will show this holds true as long as the solution is smooth. In Section 4.5 we show that this waiting time effect must be finite. This will lead us to show that solutions must lose regularity.

Remark 3.7

Notice that higher regularity of u₀ is preserved by characteristics. Taking a derivative

dudρ=((u0(Pt−1(ρ)))−α+αt)−1−1αu0−1−αdu0dρ0dPt−1dρ.=(1+αtu0(Pt−1(ρ))α)−1−1αdu0dρ01dPtdρ0=(1+αtu0(Pt−1(ρ)))α)−1−1α1+αt(u0α+m0ddρ0(u0α−1))du0dρ0.

It is easy to see that, for small time, if u₀ is smooth enough, then u is of class 𝒞¹.

Remark 3.8

The condition u0α−1 Lipschitz is sharp. Let us take, for ε > 0

(3.6) u0(ρ)=(c0−ρ)+1−εα−1

and let us show that characteristics cross for all t > 0. Looking at the characteristics for ρ₀ = c₀ − δ with δ positive but small (so that m₀(ρ) > M/2) we have that

ρ=ρ0+αm0(ρ)(c0−ρ0)+1−εt≥c0−δ+αM2δ1−εt.

But for δ<(α Mt2)1ε we have ρ > c₀. But this is not possible, since it must have crossed the characteristic ρ = c₀ coming from ρ₀ = c₀. No solutions by characteristics can exist.

The crossing of characteristics will immediately lead to the formation of shock waves. The shock waves will be led by a Rankine-Hugoniot condition as above.

Remark 3.9

As shown in Remark 3.8, with initial datum (3.6) we cannot expect solutions by characteristics. We could potentially paste solutions by characteristics on either side of a shock. We will show that this is the case, and we will show that solutions with bounded and compactly supported initial data will indeed produce a propagating shock at the end of their support, possibly with a waiting time (see the main results in Section 4).

3.4 Explicit Ansatz with waiting time

For fixed mass M and prescribed support of u = m_ρ we can construct local classical solutions with waiting time T. We will prove that

(3.7) m(t,ρ)=(Mαα−1−α1α−1(c0−ρ)+αα−1(T−t)1α−1)α−1α, if t<T and ρ>(c0−α1αM(T−t)1α)+

is a classical solution of mt+mραm=0 . We represent this function in Figure 2. We will extend this function by zero for ρ≤(c0−α1αM(T−t)1α)+ to construct a viscosity subsolution of the mass equation.

Fig. 2

Ansatz solution

The intuition to construct this kind of explicit Ansatz in “separated variables” is well known in the context of nonlinear PDEs of power type (see, e.g., [20]). One starts with a general formula of the type

m(t,ρ)=(M1β−c(c0−ρ)γ(T−t)δ)β

and we match the exponents through the scaling properties of the equation. By taking β=α−1α , γ=αα−1>1 and δ=1α−1 this gives

mt+mmρα=(M1β−c(c0−ρ)γ(T−t)δ)β−1(c0−ρ)γ(T−t)δ+1[−cδβ+(cγβ)α].

Hence, the sign is that of −cδβ + (cγβ)^α, which in our case is precisely −cα+cα. .

Remark 3.10

We have also checked the following properties

m ∈ 𝒞¹ in its domain of definition, since γ ≥ 1, the matching at ρ = c₀ is guaranteed by the explicit formula for m_ρ.
The domain of definition of m depends on the value of T, and ρ = 0 may not be contained in the domain.
It is easy to check that the function u = m_ρ associate to the Ansatz satisfies
(3.8) u(t,ρ)=cm(t,ρ)1α−1(T−t)1α−1(c0−ρ)+1α−1.
Notice that u^α−1 is a Lipschitz function of ρ. The easiest way to check this expression is by using
mαα−1=Mαα−1−c(c0−ρ)+αα−1(T−t)1α−1,
taking a derivative in ρ, and solving for m_ρ.
Conversely, notice that the condition u0α−1 in Theorem 3.4 is sharply satisfied by
u0(ρ)∼(c0−ρ)+1α−1.
If this is the behaviour of u₀ then close to ρ = c₀ we have that
(m0αα−1)ρ=αα−1m01α−1(m0)ρ=αα−1m1α−1u0∼αα−1M1α−1(c0−ρ)+1α−1
and, integrating in [ρ, c₀]
Mαα−1−m0αα−1∼M1α−1(c0−ρ)+αα−1.
Solving for m, we precisely recover the Ansatz
m0∼(Mαα−1−M1α−1(c0−ρ)+αα−1)α−1α.

3.5 A change of variable to a Hamilton-Jacobi equation

The change in variable m=θα−1α leads to the equation

(3.9) θt+(α−1α)α−1θρα=0.

This equation is of classical Hamilton-Jacobi form, and falls in the class studied by Crandall and Lions in the famous series of papers (see, e.g. [12]). Letting w = θ_ρ we recover a Burguer's conservation equation

(3.10) wt+(α−1α)α−1(wα)ρ=0.

The theory of existence and uniqueness of entropy solutions for this problem is well know (see, e.g. [6, 9, 15] and related results in [14, 17]). Notice that the relation between u and w is somewhat difficult

w(t,ρ)=dρd[(∫0ρu(t,σ) dσ)α−1α].

Remark 3.11

Notice that, for α < 1 this change of variable does not make any sense since as ρ → 0 we have that m → 0 so θ → +∞. We would be outside the L¹ ∩ L^∞ framework.

4 Viscosity solutions of the mass equation

4.1 Existence, uniqueness and comparison principle

The Crandall-Lions theory of viscosity solutions developed in [10] for α < 1, works also for the case α ≥ 1 without any modifications. Since m_ρ = u ≥ 0 we can write the equation as

(4.1) mt+(mρ)+αm=0.

Then, the Hamiltonian H(z, p₁, p₂) = (p₂)^αz is defined and non-decreasing everywhere. We write the initial and boundary conditions

(4.2) {mt+(mρ)+αm=0t,ρ>0m(t,0)=0t>0m(0,ρ)=m0(ρ)ρ>0.

The natural setting is with m₀ Lipschitz (i.e. m_ρ = u ∈ L^∞) and bounded (i.e. u ∈ L¹). We introduce the definition of viscosity solution for our problem and some notation.

Definition 4.1

Let f : Ω ⊂ ℝⁿ → ℝ. We define the Fréchet subdifferential and superdifferential

D−u(x)={p∈ℝn:lim infy→xu(y)−u(x)−p(y−x)|y−x|≥0}D+u(x)={p∈ℝn:lim supy→xu(y)−u(x)−p(y−x)|y−x|≤0}.

Definition 4.2

We say that a continuous function m ∈ 𝒞([0, +∞)²) is a:

viscosity subsolution of (4.2) if
(4.3) p1+(p2)+αm(t,ρ)≤0, ∀(t,ρ)∈ℝ+2and(p1,p2)∈D+m(t,ρ).
and m(0, ρ) ≤ m₀(ρ) and m(t, 0) ≤ 0.
a viscosity supersolution of (4.2) if
p1+(p2)+αm(t,ρ)≥0, ∀(t,ρ)∈ℝ+2 and (p1,p2)∈D−m(t,ρ).
and m(0, ρ) ≥ m₀(ρ) and m(t, 0) ≥ 0.
a viscosity solution of (4.2) if it is both a sub and supersolution.

The main results in [10] show the comparison principle and the well-posedness result for viscosity solutions of (4.2).

Theorem 4.3

Let m and M be uniformly continuous sub and supersolution of (4.2) in the sense of Definition 4.2. Then m ≤ M.

We will denote by BUC the space of bounded uniformly continuous functions.

Theorem 4.4

If m₀ ∈ BUC([0, +∞)) be non-decreasing such that m₀(0) = 0. Then, there exists a unique bounded and uniformly continuous viscosity solution. Furthermore, we have that

(4.4) ∥m(t,⋅)∥∞=limρ→+∞m(t,ρ)=∥m0∥∞, ∥mρ(t,⋅)∥∞≤∥(m0)ρ∥∞.

If m₀ is Lipschitz, then so is m.

4.2 The vortex

We next show that the mass associated to vortices (2.3) (see Figure 3) are viscosity solutions if and only if α ≥ 1. In [10] we showed that the mass associated to vortex solutions are not of viscosity type for α < 1. Intuitively, this is another instance of degenerate diffusion versus fast diffusion-like behaviour, i.e. compactly supported versus fat tails.

Fig. 3

Representation of (4.5) where c₀ = L = 1

Theorem 4.5

The mass associated to the vortex solution (2.3)

(4.5) m(t,ρ)=min{(c0−α+αt)−1αρ,c0L}.

is a viscosity solution for α ≥ 1.

Proof

Let us fix a t₀, ρ₀ > 0. If we are not on the edge, i.e. (c0−α+αt0)−1αρ0≠c0L , then m is a classical solution and it satisfies the viscosity formulation.

Let us, therefore, look at a point such that

(4.6) (c0−α+αt0)−1αρ0=c0L.

It is clear that no smooth function ϕ ≤ m can be tangent to m at (t₀, ρ₀), since the derivative of mρ(t0,ρ0−)>mρ(t0,ρ0+) . Therefore, the definition of viscosity supersolution is immediately satisfied, since D⁻m(t₀, ρ₀) = ∅.

To check that m is a viscosity subsolution, let us take (p₁, p₂) ∈ D⁺m(t₀, ρ₀). Due to the explicit formula p2∈[0,(c0−α+αt)−1α] , and since it is non-increasing in t, p₁ ≤ 0.

If p₂ = 0 then, since p₁ ≤ 0, (4.3) is satisfied. Assume p₂ ≠ 0. There exists ϕ ∈ 𝒞¹ such that ϕ ≥ m, ϕ(t₀, ρ₀) = m(t₀, ρ₀) and ϕ_t(t₀, ρ₀) = p₁ and

(4.7) φρ(t0,ρ0)=p2≤(c0−α+αt)−1α.

Let us take a look at the level sets. Since ϕ ≥ m, we have that {ϕ < c₀L} ⊂ {m < c₀L}. Therefore ∂{ϕ < c₀L}, is contained in the region {m<c0L}¯ . Since ϕ and m coincide at (t₀, ρ₀), then ∂{ϕ < c₀L} and ∂{m < c₀L} are tangent at (t₀, ρ₀) (see Figure 4).

Fig. 4

Relation between level sets. In the figure, M = c₀L.

Since p₂ > 0, ∂{ϕ < c₀L} can be parametrised by a curve (t, ρ^*(t)) defined for t ∈ (t₀ − ε, t₀ + ε) such that ϕ(t, ρ^*(t)) = c₀L and ρ^*(t₀) = ρ₀. Hence, taking a derivative and evaluating at (t₀, ρ₀)

φt(t0,ρ0)+dρ*dt(t0)φρ(t0,ρ0)=0.

On the other hand, ∂{m < c₀L} can be parametrised as (t, S(t)) where

S(t)=c0L(c0−α+αt)1α

and hence

(4.8) dSdt=c0L(c0−α+αt)−1+1α.

Notice that the sign of α − 1 decides the convexity or concavity of the matching curve. Since the level sets are tangent at (t₀, ρ₀), the derivatives coincide and we have

dρ*dt(t0)=dSdt(t0)=c0L(c0−α+αt0)−1+1α.

Hence, we have that

p1=φt(t0,ρ0)=−dρ*dt(t0)φρ(t0,ρ0)=−c0L(c0−α+αt)−1αp2.

Applying (4.7) and that α ≥ 1 we have that

p1+m(t0,ρ0)p2α=−c0L(c0−α+αt)−1αp2+c0Lp2α=c0Lp2(p2α−1−(c0−α+αt)−1α)≤0,

which is precisely (4.3). This completes the proof.

Remark 4.6

The notion of viscosity solution can be extended, by approximation, to cover non-negative finite measures as initial data u₀. Notice that these vortex solutions concentrate as t↘−c0−α/α to the Heaviside function m(t, ρ) → c₀LH₀(ρ). As we point out above, this proves that if

(4.9) u0=Mδ0 (i.e. m0=MH0),

then the solution is a cut-off of the Friendly Giant (2.1)

u(t,ρ)={(αt)−1αρ<M(αt)1α0ρ>M(αt)1α,

and hence

(4.10) m(t,ρ)={(αt)−1αρρ≤M(αt)1αMρ>M(αt)1α.

For every ε > 0 this is a viscosity solution of (4.2) in 𝒞((ε, +∞) × ℝ₊). Notice that m is of self-similar form

(4.11) m(t,ρ)=MG(ρM(αt)1α), where G(y)={yy≤11y>1.

In fact, by translation invariance, it is possible to show that if u₀ = Mδ_c₀, i.e.

m0(ρ)=MHc0(ρ),

then

(4.12) m(t,ρ)={0ρ<c0(αt)−1α(ρ−c0)ρ∈[c0,c0+M(αt)1α)Mρ>c0+M(αt)1α.

is a viscosity solution of the mass equation defined for t > 0. With self-similar form

(4.13) m(t,ρ)=MG(ρ−c0M(αt)1α), where G(y)={0y≤0yy∈(0,1)1y>1.

4.3 Monotone non-decreasing data with final cut-off

As in [10], the theory of existence and uniqueness is written in terms of m, but we take advantage of the intuition from the conservation law (P) for u, to construct explicit solutions through characteristics. Notice that taking a derivative in (4.2) we can write

ut+(uαm)ρ=0.

Afterwards, we check that the constructed solution fall into our viscosity theory for m.

Since it is suggested by (3.2) that characteristics do not cross if u₀ is non-decreasing, let us look for solutions with initial data

(4.14) u0(ρ)={u˜0(ρ)ρ<S00ρ>S0., where u˜0 in continuous and non-decreasing in [0,S0].

When m₀ and u₀ are non-decreasing, it clear that characteristics with foot on ρ₀ ∈ [0, S₀] do not cross. If ũ₀ ≢ 0, then u0(S0−)=u˜0(S0)>0 and there is a shock starting at S₀ which will propagate following the Rankine-Hugoniot condition (3.5).

We construct the viscosity solution as follows. The characteristic of foot ρ₀ = S₀ is precisely

S¯(t)=S0+αMu˜0(S0)α−1t.

For ρ<S¯(t) we can go back through the characteristics with P_t defined above. Let us define

(4.15) u˜(t,ρ)={0if u˜0(Pt−1(ρ))=0(u˜0(Pt−1(ρ))−α+αt)−1αif u˜0(Pt−1(ρ))>0 and ρ<S¯(t)0if ρ>S¯(t).

The shock is given by

(4.16) {dSdt=Mu˜(t,S(t))α−1S(0)=S0.

Finally we define

(4.17) u(t,ρ)={u˜(t,ρ)ρ<S(t)0ρ>S(t).

Solving explicitly is not possible. However, we can prove that

Proposition 4.7

Let α ≥ 1, (4.14). Then, the mass of (4.17) is a viscosity solution of (4.2) and S(t)≤S¯(t) .

Proof

Since α ≥ 1, we have that

dSdt=m(t,S(t))u(t,S(t)−)α−1=m(t,S(t))u˜(t,S(t))α−1≤αMu˜0(S0)α−1=dS¯dt.

where M = sup m₀. Also S(0)=S¯(0) . Hence the shock is slower than the last characteristic (i.e. S(t)≤S¯(t) ) and this implies that there are no outgoing characteristics so the Lax-Oleinik condition is satisfied.

Now, in order to check that it is a viscosity solution, we can repeat the proof of Theorem 4.5, replacing (4.8) by (4.16).

4.4 Two Dirac deltas

Let us now consider that

(4.18) u0=m1δρ1+m2δρ2.

where 0 ≤ ρ₁ < ρ₂. Then the initial mass m₀ is discontinuous, and this creates some technical difficulties. We will show the viscosity solution for (4.2) is the primitive of

(4.19) u(t,ρ)={0ρ<ρ1(αt)−1αρ∈[ρ1,S1(t)]0ρ∈[S1(t),ρ2]((ρ−ρ2αm1t)−αα−1+αt)−1αρ∈[ρ2,S2(t)]0ρ>S2(t)

where

(4.20) S1(t)=ρ1+m1(αt)α−1α,

and

(4.21) S2(t)=ρ2+αm1K−1(m2αm1)t1α,

with

(4.22) K(τ)=∫0τ(s−αα−1+α)−1αds.

We have the following estimates for the function 𝒦⁻¹: for τ ≤ s₀ there exists c(s₀), C(s₀) > 0 such that

(4.23) c(s0)ταα−1≤K−1(τ)≤C(s0)ταα−1.

This solution is defined for all t such that S₁(t) ≤ ρ₂, i.e. t≤((ρ2−ρ1)/m1)αα−1/α. . For large t, S₁ would need to be computed from another further Rankine-Hugoniot condition. We will only use the value for t small, so this computation is enough for our purposes.

4.4.1 Approximation by viscosity solutions

We will prove there is an explicit solution defined for some T > 0 which is of viscosity type for t > 0. We will approximate the initial data by

(4.24) u0ε,δ={m1ερ∈[ρ1,ρ1+ε]m2ε(ρ−(ρ2−δ))δρ∈[ρ2−δ,ρ2]m2ερ∈[ρ2,ρ2+ε]0otherwise,

for ε and δ small enough. The ε-regularisation is used to approximate the Dirac deltas at the level of u. The δ-regularisation is used to resolve the appearance of a rarefaction wave at ρ₂ due to a gap in the characteristics. Since viscosity solutions are stable by passage to the limit, we only need to show that our approximating solution are viscosity solutions.

The first part of the solutions does not notice the δ-regularisation. We take ε small enough so that ρ₁ + ε < ρ₂ − 2ε. For ρ < ρ₂ we reconstruct a vortex type solution following Section 4.2 with an initial gap

(4.25) uε(t,ρ)={0ρ<ρ1((m1ε)−α+αt)−1αρ∈[ρ1,S1ε(t)]

where the first shock is given by

(4.26) S1ε(t)=ρ1+ε1+m1(((m1ε)−α+αt)α−1α−(m1ε)1−α).

Solutions in this form are defined for t ∈ [0, T_ε) such that S1ε(Tε)=ρ2−δ . We leave to the reader to check that T_ε does not tend to zero with ε → 0.

For the second part, the characteristics with foot ρ₀ ∈ [ρ₂ − δ, ρ₂] are given by

(4.27) ρ=ρ0+α(m1+m22(ρ0+δ−ρ2)2εδ)(m2εδ(ρ2+δ−ρ))α−1t.

On the other hand, if ρ₀ ∈ [ρ₂, ρ₂ + ε] we have

(4.28) ρ=ρ0+α(m1+m22εδ+m2ε(1−(ρ2+ε−ρ0)))(m2ε)α−1t.

Notice that in both cases u is given by

u(t,ρ)=(u0(ρ0)−α+αt)−1α.

By mass conservation we have a further shock starting from ρ₂ + ε given by a Rankine-Hugoniot condition

dS2ε,δdt(t)=(m1+m2)uε,δ(t,S2ε(t)−)α−1.

The first part of solution is of viscosity type, by an argument analogous to Section 4.2 and the second part have a monotone non-decreasing datum with final cut-off as in Section 4.3. We are reduced now to pass to the limit as ε and δ tend to 0.

4.4.2 Passage to the limit as δ → 0

For [0, ρ₂ − δ] the solution did not depend on δ, so there is no work needed. Applying a similar argument as in [10] we can pass to the limit in (4.27). The characteristics with foot in [ρ₂, ρ₂ + δ] collapse to a rarefaction fan at ρ₂ of the form

(4.29) ρ=ρ2+m1η0α−1t, η0∈[0,m22ε].

By inverting η₀ in (4.29) we recover the solution

uε(t,ρ)=(η0−α+αt)−1α=((ρ−ρ2αm1t)−αα−1+αt)−1α.

The other characteristics are for foot ρ₀ ∈ [ρ₂, ρ₂ + ε] and, by passing analogously to the limit in (4.28), we have

ρ=ρ0+α(m1+m2ε(1−(ρ2+ε−ρ0)))(m2ε)α−1t.

Since u0ε is non-decreasing the characteristics do not cross. The Rankine-Hugoniot condition is now

dS2εdt(t)=(m1+m2)uε(t,S2ε(t)−)α−1.

4.4.3 Passage to the limit as ε → 0

Passing to the limit we end up only with the rarefaction fan characteristics and recover (4.19) where the first shock is given by (4.20) and the second shock, S₂ which defines the support, is a solution of the ODE

(4.30) {dS2dt(t)=(m1+m2)((S2(t)−ρ2αm1t)−αα−1+αt)−α−1α,S2(0)=ρ2.

Notice that this equation is singular at t = 0 but it can be rewritten as

dS2dt(t)=(m1+m2)t−α−1α(t1α−1(S2(t)−ρ2αm1)−αα−1+α)−α−1α.

Since α−1α∈(0,1) the Cauchy problem is well-posed. Alternatively, one can write S₂ implicitly as the only value such that

∫ρ2S2(t)u(t,ρ) dρ=m2.

In other words,

(4.31) ∫ρ2S2(t)((ρ−ρ2αm1t)−αα−1+αt)−1αdρ=m2.

This solution is defined for 0 ≤ t < T where

T=1α(ρ2−ρ1m1)αα−1.

By scaling analysis on the integral, we can give an algebraic expression of S₂(t). We apply the change of variables ρ=ρ2+αm1st1α to deduce

m2=∫ρ2S2(t)((ρ−ρ2αm1t)−αα−1+αt)−1αdρ=∫0S2(t)−ρ2αm1t1α((st−α−1α)−αα−1+αt)−1αt1ααm1ds=αm1∫0S2(t)−ρ2αm1t1α(s−αα−1+α)−1αds=αm1K(S2(t)−ρ2αm1t1α).

Hence, we recover (4.21). To show (4.23) we simply indicate that, for s ≤ s₀

s−αα−1≤s−αα−1+α≤C(s0)s−αα−1

4.5 Waiting time

4.5.1 Existence of waiting time

We turn the explicit solution in (3.7) into a viscosity subsolution by extending it by zero, that is we define m(t, ρ) as

(4.32) m_(t,ρ):={0if ρ<c0−α1αM(T−t)1α(Mαα−1−α1α−1(c0−ρ)αα−1(T−t)1α−1)α−1αif ρ∈(c0−α1αM(T−t)1α,c0)Mif ρ≥c0.

Proposition 4.8

m(t, ρ) is a viscosity subsolution of mt+mmρα=0 .

Proof

It is clear that 0 is a solution of mt+mmρα=0 . So is the second part for ρ>c0−c1αM(T−t)1α , as we have checked in Section 3.4. At the matching point ρ=c0−c1αM(T−t)1α , we have that

mρ(t,ρ−)=0, mρ(t,ρ+)=c(0+)−1α(c0−ρ)1α(T−t)1α=+∞

This corner does not allow any smooth ϕ to be tangent from above at this point, and hence the condition of viscosity subsolution is trivially satisfied.

We will denote by c₀ = max supp u₀, where u₀ = (m₀)_ρ, that coincides with the boundary of m₀ = M in the sense that

(4.33) m0(ρ)<M for ρ<c0 and m0(ρ)=M for ρ≥c0.

Corollary 4.9

Let m₀ ∈ BUC([0, +∞)) and let c₀ = max supp u₀. If

(4.34) lim supρ→c0−M−m0(ρ)(c0−ρ)αα−1<+∞,

then there is waiting time as in Corollary 3.5.

Proof

First, we prove that

supρ∈[0,c0]M−m0(ρ)(c0−ρ)αα−1<+∞.

Let ρ_k be such that

M−m0(ρk)(c0−ρk)αα−1→supρ∈[0,c0]M−m0(ρ)(c0−ρ)αα−1.

If the supremum were infinite, since M − m₀(ρ_k) is bounded, then we have that ρ_k → c₀. This results in

limkM−m0(ρk)(c0−ρk)αα−1≤lim supρ→c0−M−m0(ρ)(c0−ρ)αα−1<+∞

leading to a contradiction.

Therefore, there exists C > 0 such that for all ρ ∈ [0, c₀]

M−m0(ρ)(c0−ρ)αα−1≤C.

In particular, we have that

m0(ρ)≥M−C(c0−ρ)αα−1.

We can apply the convexity of the function f(x)=xαα−1 to show that

m0(ρ)α−1α≥Mαα−1−αα−1M1α−1C(c0−ρ)αα−1=Mαα−1−α1α−1T1α−1(c0−ρ)+αα−1=m_(0,ρ)α−1α,

for a well chosen T (see, e.g. Figure 6). Therefore, applying the comparison principle Theorem 4.3 then m ≥ m, and thus m has waiting time.

Fig. 5

Solutions with u₀ given by two characteristics. Computed with the numerical scheme for m in Section 5 reproducing the exact solution up to approximation error. The function u = m_ρ is recovered by numerical differentiation.

Fig. 6

The explicit Ansatz viscosity subsolution (4.32) guarantees existence of waiting time. The subsolution is represented from the explicit solution, whereas u is computed through the numerical scheme in Section 5. See a movie simulation in the supplementary material [1, Video 1].

4.5.2 Non-existence of waiting time

Theorem 4.10

Let m₀ ∈ BUC([0, +∞)) and let c₀ = max supp u₀. If

(4.35) lim supρ→c0−M−m0(ρ)(c0−ρ)αα−1=+∞,

then there is no waiting time.

Proof

To prove there is no waiting time, we want to show that S(t) > c₀ for any t > 0. In order to do this, we will construct a sequence of supersolutions m¯k with Sk(t)=maxsupp(m¯k)ρ(t,⋅) and times t_k ↘ 0 such that S_k(t_k) > c₀. This ensures that, for some k we have 0 < t_k < t and hence S(t) ≥ S_k(t) ≥ S_k(t_k) > c₀.

Let us consider a sequence d_k such that

dk→lim supρ→c0−M−m0(ρ)(c0−ρ)αα−1.

There exists ρ_k ↗ c₀ such that

(4.36) M−m0(ρk)≥dk(c0−ρk)αα−1.

We construct the viscosity supersolutions m¯k with initial derivative

uk,0=m0(ρk)δ0+(M−m0(ρk))δρk.

It is clear that m¯k(0,ρ)≥m(0,ρ) . By using the comparison principle Theorem 4.3, m¯k≥m .

Now we apply the theory of Section 4.4. We will select t_k > 0 such that S_k(t) ≥ c₀ + ε_k for t ≥ t_k with εk=c0−ρk2>0 . Using (4.21), (4.22) and (4.36) we have that

Sk(t)=ρk+αm0(ρk)K−1(M−m0(ρk)αm0(ρk))t1α≥ρk+αm0(ρk)K−1(dk(c0−ρk)αα−1αm0(ρk))t1α.

Due to our choice of ρ_k, it is sufficient that

ρk+αm0(ρk)K−1(dk(c0−ρk)αα−1αm0(ρk))t1α≥c0+εk.

Solving for t, we have that

(4.37) t≥((c0−ρk)+εkαm0(ρk)K−1(dk(c0−ρk)αα−1αm0(ρk)))α.

We know that dk(c0−ρk)αα−1≤M−m0(ρk)→0 , therefore we need to study 𝒦⁻¹ close to 0. Going back to (4.23) there exists a constant C > 0 such that

((c0−ρk)+εkαm0(ρk)K−1(dk(c0−ρk)αα−1αm0(ρk)))α≥C((c0−ρk)+εkαm0(ρk)(dk(c0−ρk)αα−1αm0(ρk))α−1α)α=3α2αCdkα−1=:tk.

Due to the hypothesis of the theorem d_k → +∞ and hence t_k → 0.

See a movie simulation of one of the mass supersolutions interacting with a solution without waiting time in the supplementary material [1, Video 2].

Remark 4.11

Notice that if the lim sup is finite, the previous proof can be adapted to show that the supersolutions m¯k give an upper bound of the waiting time.

Remark 4.12

As pointed out in Corollary 3.5, the spatial support of classical solutions does not change in time. Taking c₀ = max supp u₀, we construct the supersolution m¯ with initial derivative

u¯0=m0(c02)δ0+(M−m0(c02))δc02.

This supersolution shows that the support of u must move after a finite time and therefore that the solution is no longer classical.

4.6 Asymptotic behaviour

We give first a general result of asymptotic behaviour in mass variable.

Theorem 4.13

Assume that u₀ ∈ L^∞(0, ∞) has compact support, M = ‖u₀‖_L¹, m be the viscosity solution with initial data m₀ and S(t) = inf{ρ : m(t, ρ) = M}. Then S(t)∼M(αt)1α with estimate

(4.38) 0≤S(t)M(αt)1α−1≤S(0)M(αt)−1α.

Furthermore, m has the asymptotic profile in rescaled variable y=ρM(αt)1α with an asymptotic estimate

(4.39) supy≥ε|m(t,M(αt)1αy)MG(y)−1|→0, as t→+∞ where G(y)={yy≤11y>1

for any ε > 0.

Proof

By Remark 4.6 we take as super and subsolution m¯ and m with initial data

m¯0(ρ)=MH0(ρ), m_0(ρ)=MHS(0)(ρ).

Hence m¯≥m≥m_ . Due to the explicit form of m¯ and m, we have that

M(αt)1α≤S(t)≤S(0)+M(αt)1α.

Due to the self-similar form of m¯ and m given in Remark 4.6, the result is proven.

Remark 4.14

Notice that for m₀ = 0 in [0, a], we have m(t, ρ) = 0 in [0, a] so the supremum of y ≥ 0 is always 1. If u₀ is continuous and u₀(0) > 0, then the supremum can be taken for y ≥ 0.

Let us discuss the asymptotic behaviour when the datum is monotone non-decreasing with final cut-off. We recall (4.14)–(4.17) and define

U(t,ξ)=u(t,(αt)1αξ)(∥u0∥L∞−α+αt)−1α.

Since the solution is constructed by characteristics we have that

U(t,ξ)={0if u0((αt)−1αξ)=0(η0(t,ξ)−α+αt∥u0∥L∞−α+αt)−1αif u0((αt)−1αξ)>0 and (αt)1αξ≤S(t)0if (αt)1αξ>S(t)

where η₀(t, ξ) ∈ (0, ‖u₀‖_L^∞). Due to (4.38), as t → +∞ we have that

U(t,ξ)→{1if ξ∈(0,M)0if ξ∈(M,+∞).

The value at 0 depends on whether u₀(0) = 0.

5 A numerical scheme

In the pioneering paper by Crandall and Lions [11], the authors developed a theory of monotone schemes for finite differences of Hamilton-Jacobi equations, where solutions are shown to converge to the viscosity solution. They study equations of the form

(5.1) mt+H(mρ)=0.

For these equations it is natural to develop only explicit methods. However, for our case m_t + H(m_ρ)m = 0, we will see that it more natural, and probably more stable, to do an explicit-implicit approximation of the non-linear term H(m_ρ)m. In fact, since the nonlinear term is linear in m, we can solve for the implicit step in an explicit manner. More precisely, considering an equispaced discretization

(5.2) tn=htn ρj=hρj.

We select the following finite-difference schemes

Mjn+1−Mjnht+(Mjn−Mj−1nhρ)αMjn+1=0

which can be written as

(5.3) Mjn+1=Mjn1+ht(Mjn−Mj−1nhρ)+α=G(Mjn,Mj−1n).

Here, G is given by

G(p,q)=p1+htH(p−qhρ), where H(s)=s+α.

Notice that the method depends only on the parameter ht/hρα . Taking derivatives we have that

∂G∂p=1+htH(p−qhρ)−hthρH′(p−qhρ)p(1+htH(p−qhρ))2, ∂G∂q=phtH′(p−qhρ)hρ(1+htH(p−qhρ))2≥0.

Then G is non-decreasing in p under the simple CFL condition:

(5.4) hthρH′(p−qhρ)p≤12.

Since the denominator in G is larger than 1, we have that G(p, q) ≤ p. This is immediately translated to a maximum principle for Mjn

(5.5) Mjn+1≤Mjn≤∥m0∥∞.

For α ≥ 1 we have two options to obtain a CFL condition. We can check whether the numerical derivative is bounded (this can be done for some methods, see Section 5.2) or cut-off the equation by a nice value. For m₀ fixed, since m_ρ ≤ ‖(m₀)_ρ‖_L^∞ due to (4.4), the equation (5.1) where H(s)=s+α is equivalently to itself with

(5.6) H(s)=(max{s,∥(m0)ρ∥∞})+α.

We write this cut-off to ensure monotonicity. Nevertheless, once the method is monotone, Lemma 5.3 ensures that the cut-off part is not reached. Hence, this cut-off is purely technical.

For α ≥ 1 this new H given by (5.6) satisfies

0≤H′(s)≤α∥(m0)ρ∥∞α−1.

Therefore, (5.4) can be taken as

(CFL) hthρ≤12α∥(m0)ρ∥∞α−1∥m0∥∞.

We propose the scheme

(M) {Mjn+1=Mjn1+htH(Mjn−Mj−1nhρ)if j>0,n≥0Mj0=m0(hρj)if j≥0M0n=0if n>0.

Remark 5.1

As we pointed out in [10], for 0 < α < 1 this method is not monotone. This was fixed by regularising H. For δ > 0 we take

(5.7) Hδ(s)=(s++δ)α−δα.

By including the boundary and initial condition, we constructed the method

(M_δ) {Mjn+1=Mjn1+htHδ(Mjn−Mj−1nhρ)if j>0,n≥0Mj0=m0(hρj)if j≥0M0n=0if n>0.

with this regularisation we know that Hδ′(s)≤αδα−1 so we have a CFL condition

(CFL_δ) hthρ≤δ1−α2α∥m0∥∞.

In [10] we made δ to converge to 0 with h_t and h_ρ, showing the convergence of the numerical solutions.

5.1 Properties of monotone methods

The following properties of (M) when G is monotone in each variable are a classical matter (see the original result in [11] and the presentation and references in [2]). We just briefly sketch them for completeness.

Lemma 5.2

Let α ≥ 1, m₀ ≥ 0 and bounded and consider the sequence Mjn constructed by (M) and assume (CFL). We have the following properties:

1. Mjn+1=G(Mjn,Mj−1n) where G is non-decreasing.
2. Mjn≤∥m0∥∞
3. If m₀ ≥ 0 is non-decreasing then:
- (a) 0≤Mjn≤Mj+1n for all n, j
- (b) There is mass conservation in the numerical scheme
  M∞n+1=limj→+∞Mjn+1=limj→+∞Mjn=M∞n.

Proof

We have shown this above.
This is true for Mj0 by construction, and hence for all n, due to the previous item.
1. We proceed by induction in n. For time n = 0 this is true due to the monotonicity of m₀. Assume Mjn≤Mj+1n for all j. Since G is monotone in each coordinate
  Mj+1n+1=G(Mj+1n,Mjn)≥G(Mjn,Mjn)≥G(Mjn,Mj−1n)=Mjn+1.
2. Since the sequence is non-decreasing and bounded, it has a limit. Furthermore limj(Mjn−Mj−1n)=0 . Hence, since H_δ(0) = 0 we have that
  M∞n+1=limj→+∞Mjn+1=limj→+∞Mjn1+htH(Mjn−Mj−1nhρ)=M∞n.

Notice the biggest advantage of the method (M) is that it preserves the space monotonicity of m and the total mass, as it should be for a mass equation.

5.2 Convergence of the numerical scheme (M) to the viscosity solution

In order to construct a convergent scheme, it is better to work with a single parameter. For h > 0 we define

hρ=h, ht=h2α∥(m0)ρ∥∞α−1∥m0∥∞.

so that (CFL) is satisfied. We now allow M jn to be constructed from (M). For t_n ≤ t < t_n+1 and ρ_j ≤ ρ < ρ_j+1 we write the piecewise linear interpolation of the discrete values as

(5.8) mh(t,ρ)={M jn+(ρ−ρj)M j+1n−M jnhρ+(t−tn)M jn+1−M jnhtif ρ≤tM j+1n+1−(ρj+1−ρ)M j+1n+1−M jn+1hρ−(tn+1−t)M j+1n+1−M j+1nhtif ρ>t

This construction ensures that

∂mh∂ρ={Uj+1nif ρ<tUj+1n+1if ρ>t,∂mh∂t={−H(Ujn)Mjn+1if ρ<t−H(Uj+1n)Mj+1n+1if ρ>t

where Ujn is the numerical space derivative

(5.9) Ujn=Mjn−Mj−1nhρ≥0,

and the numerical time derivative is given by the relation

Mjn+1−Mjnht=−H(Ujn)Mjn+1≤0.

The strategy of the proof is the following. We will show that these space and time numerical derivatives are uniformly bounded, and hence m^h is uniformly continuous, non-decreasing in ρ and non-increasing in t. We can then apply the Ascoli-Arzelá precompactness theorem and show there exists a convergent subsequence. We will prove the limit is the viscosity solution.

If we subtract (M) for j and j − 1 we recover an equation for the numerical derivative Ujn

(5.10) Ujn+1−Ujnht+H(Ujn)Mjn+1−H(Uj−1n)Mj−1n+1hρ=0.

Notice that the natural scaling for this equation is h_t/h_ρ.

5.2.1 Boundedness of the numerical derivative

Since (5.10) is a numerical approximation by a monotone method of the nonlinear conservation law (P), we can expect a maximum principle.

Lemma 5.3

Let 0 ≤ m₀ be uniformly Lipschitz continuous, bounded and non-decreasing, Mjn be given by (M), that (CFL) holds and let Ujn given by (5.9). Then, Ujn≥0 and

(5.11) supjUjn+1≤supjUjn ∀n≥0.

Remark 5.4

Once this is proven, the cut-off (5.6) is not needed.

Proof

That Ujn≥0 follows form Lemma 5.2. We write

0=Ujn+1−Ujnht+Mj−1n+1H(Ujn)−H(Uj−1n)hρ+H(Ujn)Mjn−Mj−1nhρ=Ujn+1−Ujnht+Mj−1n+1H(Ujn)−H(Uj−1n)hρ+H(Ujn)Ujn.

Solving for Ujn+1 , using the fact that H is non-decreasing and Ujn≥0 , we have that

Ujn+1≤Ujn−hthρMj−1n+1(H(Ujn)−H(Uj−1n))=Ujn−hthρMjn+1H′(ξjn)Ujn+hthρMjn+1H′(ξjn)Uj−1n=(1−hthρMjn+1H′(ξjn))Ujn+hthρMjn+1H′(ξjn)Uj−1n.

Due to (CFL) we have that the coefficients in front of Ujn and Uj−1n are non-negative. Hence

Ujn+1≤(1−hthρMjn+1H′(ξjn))supj Ujn+hthρMjn+1H′(ξjn)sup jUjn=supj Ujn.

And this holds for every j so the result is proved.

5.2.2 Convergence via Ascoli-Arzelá. Existence of a viscosity solution

Theorem 5.5

Let us m₀ ∈ W^1,∞(0, +∞) and non-decreasing, (CFL), Mjn constructed by (M) and m^h be given by (5.8). Then, m^h is a family of uniformly continuous functions. Then, for every P > 0

(5.12) mh→m in C([0,P]×[0,T]) as hρ→0.

where m is a viscosity solution of (4.2). Furthermore, (4.4) holds.

Proof

First, we notice that m^h satisfies (4.4). Due to (5.8), we have that

By the Ascoli-Arzelá theorem there is a subsequence that converges uniformly in [0, P] × [0, T].

We will show every convergent subsequence converges to the same limit m, and hence the whole sequence converges. We still denote by h the indices of the convergent subsequences.

Let m^h be a subsequence converging in [0, T] × [0, P]. We check that it is a viscosity subsolution, and the proof of viscosity supersolution is analogous. Let (t₀, ρ₀) ∈ (0, T) × (0, P) and ϕ ∈ C² be such that m − ϕ has a strict local maximum at (t₀, ρ₀) and m(t₀, ρ₀) = ϕ(t₀, ρ₀). We can modify ϕ outside a bounded neighbourhood of (t₀, ρ₀), so that m − ϕ attains a unique global maximum at (t₀, ρ₀), for h large enough m^h − ϕ attains a global maximums in [0, T] × [0, P] at an interior points (t_h, ρ_h), and (t_h, ρ_h) → (t₀, ρ₀) as h → 0. Our argument is a variation of [19, Lemma 1.8].

Let B ⊂ [0, T]×[0, P] be a small open ball around (t₀, ρ₀) where the maximum is global. Let ε = − inf_B(m− ϕ)/2. Define U = {m − ϕ > −ε} ∩ B which is a open and bounded neighbourhood of (t₀, ρ₀). We modify ϕ so that is greater than m + ε outside U. With the modification, m − ϕ attains a unique global maximum at (t₀, ρ₀).

Let h be small enough so that |mh−m|<ε2 in [0, T] × [0, P]. We have that

max[0,t0+1]×[0,ρ0+1]\U(mh−φ)<max[0,t0+1]×[0,ρ0+1]\U(m−φ)+ε2≤−ε2.

On the other hand

mh(t0,ρ0)−φ(t0,ρ0)>m(t0,ρ0)−φ(t0,ρ0)−ε2=−ε2.

Therefore the maximum over [0, T]×[0, P] is attained at some (t_h, ρ_h) ∈ U. The sequence (t_h, ρ_h) is bounded, and therefore as a convergent subsequence. Let (t₁, h₁) be its limit. We have that

mh(th,ρh)−φ(th,ρh)≥mh(t,ρ)−φ(t,ρ) ∀(t,ρ)∈[0,T]×[0,P].

Passing to the limit, since the maximum is unique, we have that (t₁, ρ₁) = (t₀, ρ₀). Since all convergent subsequences share a limit, the whole sequence converges.

For such small values of h, let us define

nh=⌊thht⌋−1, jh=⌈ρhhρ⌉.

Since m^h − ϕ has a global maximum in [0, T] × [0, P], we have that

mh(th,ρh)−φ(th,ρh)≥mh(t,ρ)−φ(t,ρ).

Evaluating on the points of the mesh, we get that

Mjn≤φ(tn,ρj)−φ(th,ρh)+mh(th,ρh).

Since m^h is increasing in ρ and decreasing in t and the fact that G is non-decreasing, we recover

mh(th,ρh)≤mh((nh+1)ht,jhhρ)=Mjhnh+1=G(Mjhnh,Mjh−1nh)≤G(φ(tnh,ρjh)−φ(th,ρh)+mh(th,ρh),φ(tnh,ρjh−1)−φ(th,ρh)+mh(th,ρh))=φ(tnh,ρjh)−φ(th,ρh)+mh(th,ρh)1+htH(φ(tnh,ρjh)−φ(tnh,ρjh−1)hρ)

due to the definition of G. Since ϕ is smooth, for h small enough the denominator is positive and hence

φ(th,ρh)−φ(tnh,ρjh)ht+H(φ(tnh,ρjh)−φ(tnh,ρjh−1)hρ)mh(th,ρh)≤0.

Adding and subtracting ϕ(t_{n_h+1}, ρ_j)/h_ρ on both sides

φ(tnh+1,ρjh)−φ(tnh,ρjh)ht+H(φ(tnh,ρjh)−φ(tnh,ρjh−1)hρ)mh(th,ρh)≤φ(tnh+1,ρjh)−φ(th,ρh)th−tnhth−tnhht.

Clearly t_h − t_{n_h} ≥ 0 and, since ϕ is of class C¹ and m is non-increasing in t, we have that

limh→0φ(tnh+1,ρjh)−φ(th,ρh)th−tnh=φt(t0,ρ0)≤0,

Therefore, as h → 0, we conclude that

φt(t0,ρ0)+H(φρ(t0,ρ0))m(t0,ρ0)≤0,

for any ϕ such that m − ϕ has a global maximum at (t₀, ρ₀). Therefore, m is a viscosity subsolution.

5.3 Rate of convergence

Theorem 5.6

Let α ≥ 1 and let h_t and h_ρ satisfy (CFL). Let m₀ be Lipschitz continuous and bounded and let m be the viscosity solution of (4.2) and Mjn be constructed by (M). Then, for any T > 0

supj≥00≤n≤T/ht|m(tn,ρj)−Mjn|≤Chρ13.

where C does not depend on h_ρ.

Remark 5.7

The original paper by Crandall and Lions allows, by a longer and more involved argument, to prove estimates of the form O(ht) with H continuous, but requiring that the function defining the method is Lipschitz continuous.

Proof

For convenience, in the proof we denote N = ⌊T/h_t⌋. Our aim is to prove that

σ=supj≥00≤n≤T/ht(m(tn,ρj)−Mjn)≤Chρ13.

The argument can be analogously repeated for the infimum. If σ ≤ 0 there is nothing to prove. Assume that σ > 0. Let L be the Lipschitz constant of m₀. Due to (4.4), it is also the Lipschitz constant of m.

We begin by indicating there exist n₁, j₁ such that

m(t1,ρ1)−Mj1n1≥3σ4, where t1=htn1 and ρ=hρj1.

We define

Φ(t,htn,ρ,hρj)=m(t,ρ)−Mjn−ϕ(t,htn,ρ,hρj)

where, for ε, λ > 0 we define

ϕ(t,s,ρ,ξ)=(|ρ−ξ|2+|t−s|2ε2+λ(t+s))

Then the maximum is at t_ε ∈ [0, T], ρ_ε ∈ [0, +∞) and t_ε = h_tn_ε with n_ε ∈ {0, · · · , N}, ξ_ε = h_ρj_ε with j_ε ∈ ℕ ∪ {0}. Again this function is continuous and

Defined over a bounded set in t_ε and s_ε.
If ρ → +∞ and j remains bounded then Φ → −∞ (analogously in ρ bounded and j → +∞).
If ρ → +∞ and j → +∞ then
lim supρ,j→+∞Φ≤m∞−m∞≤0,

so there exists a point of maximum (t_ε, h_tn_ε, ρ_ε, h_ρj_ε) such that

Φ(tε,sε,ρε,ξε)≥Φ(t,s,ρ,ξ) ∀(t,s,ρ,ξ).

In particular

Φ(tε,sε,ρε,ξε)≥Φ(t1,t1,ρ1,ρ1)=m(t1,ρ1)−Mj1n1−2λt1.

Taking

(5.13) λ=σ8(1+T)

we have

Φ(tε,ρε,nε,jε)≥σ2.

In particular,

(5.14) m(tε,ρε)−Mjεnε≥σ2+ϕ(tε,htnε,ρε,hρjε)>0

Step 1. Variables collapse. As Φ(t_ε, s_ε, ρ_ε, ξ_ε) ≥ Φ(0, 0, 0, 0) = 0, we have

|ρε−ξε|2+|tε−sε|2ε2+λ(tε+sε)≤m(tε,ρε)−Mjεnε≤2∥m0∥∞.

Therefore, we obtain

|ρε−ξε|2+|tε−sε|2≤2∥m0∥∞ε2, and ρε2+ξε2≤2∥m0∥∞ε.

This implies that, as ε → 0, the variable doubling collapses to a single point.

Step 2. For ε small enough, t_ε, ρ_ε, n_ε, j_ε > 0. Since m is Lipschitz continuous

σ2<m(tε,ρε)−Mjεnε=m(tε,ρε)−m(0,ρε)+m(0,ρε)−m(0,ξε) +m(0,ξε)−Mjε0+Mjε0−Mjεnε≤Ltε+L|ρε−ξε|,

using the fact that m(0,ξε)=Mjε0 and Mjn is decreasing in n. If ε is small enough

(5.15) ε<σ4L2∥m0∥∞,

we have Lt_ε > σ/4 and hence t_ε > 0. Analogously for ρ_ε > 0.

If n_ε = 0 then

σ2<m(tε,ρε)−Mjε0=m(tε,ρε)−m(0,ρε)+m(0,ρε)−m(0,ξε)+m(0,ξε)−Mjε0≤Ltε+L|ρε−ξε|=L|tε−nε|+L|ρε−ξε|≤L2∥m0∥∞ε.

This is a contradiction if (5.15) holds. An analogous contradiction holds if j_ε = 0.

Step 3. An inequality for m via viscosity. We check that

(t,ρ)↦m(t,ρ)−ϕ(t,sε,ρ,ξε)=m(t,ρ)−ψ(t,ρ)

has a maximum at (t_ε, ρ_ε). Hence, since m is a viscosity subsolution, we have

∂ϕ∂t(tε,sε,ρε,ξε)+H(∂ϕ∂ρ(tε,sε,ρε,ξε))m(tε,ρε)≤0.

Computing the derivatives

(5.16) 2(tε−sε)ε2+λ+H(2(ρε−ξε)ε2)m(tε,ρε)≤0.

Step 4. An inequality for M applying that G is monotone. As before, the function

(n,j)↦Mjn−(−ϕ(tε,htn,ρε,hρj))=Mjn−ψ(j,n)

has a minimum at (n_ε, j_ε). In particular, we obtain

Mjn≥Mjεnε−ψ(jε,nε)+ψ(j,n).

Since G is monotone, it yields

Mjεnε=G(Mjεnε−1,Mjε−1nε−1)≥G(Mjεnε−ψ(jε,nε)+ψ(jε,nε−1)︸S1,Mjεnε−ψ(jε,nε)+ψ(jε−1,nε−1)︸S2).

Similarly to the proof of Theorem 5.5, for h small, one can rewrite the previous inequality as

Mjεnε−S1ht+H(S1−S2hx)Mjεnε≥0.

Hence, we recover

(5.17) ψ(jε,nε)−ψ(jε,nε−1)ht+H(ψ(jε,nε−1)−ψ(jε−1,nε−1)hx)Mjεnε≥0.

We could compute this explicitly, but it is sufficient and clearer to apply the intermediate value theorem to deduce

−∂ϕ∂s(tε,s¯ε,ρε,ξε)+H(−∂ϕ∂ξ(tε,sε−ht,ρε,ξ¯ε))Mjεnε≥0

where s¯ε∈(sε−ht,sε) and ξ¯ε∈(ξε−hρ,ξε) . Hence, we conclude that

(5.18) 2(tε−s¯ε)ε2−λ+H(2(ρε−ξ¯ε)ε2)Mjεnε≥0.

Step 5. An estimate for σ. Substracting (5.16) from (5.18) we have that

σ4(1+T)≤sε−s¯εε2+H(2(ρε−ξ¯ε)ε2)Mjεnε−H(2(ρε−ξε)ε2)m(tε,ρε)≤(H(2(ρε−ξ¯ε)ε2)−H(2(ρε−ξε)ε2))Mjεnε+H(2(ρε−ξ¯ε)ε2)(Mjεnε−m(tε,ρε)).

Notice that the second term is non-positive due to (5.14). We now use the Lipschitz continuity of H, which holds for the cut-off given by (5.6), and we obtain that

σ8≤C|2(ξ¯ε−ξε)ε2|≤Chρε2∥m0∥L∞.

Step 6. A first choice of ε. We take ε = Cσ where C is chosen so that (5.15) hold. Then, we have that

σ3≤Chρ.

where C(T, α) is independent of h_t, h_x or δ. This completes the proof.

Remark 5.8

Notice that we do not use the equation until step 4 and that the Lipschitz continuity of m₀ plays a key role. However, the homogeneous boundary conditions do not.

Remark 5.9

Notice that we recover the exponent h13 from the Lipschitz continuity of H. If H is only α-Hölder continuous as in [10], then the rate of convergence is given by hα1+2α .

6 Numerical results

6.1 Asymptotics as t → +∞

Through numerical experiments, we see that the vortex seems to be the asymptotic solution also in u variable. In Figure 7 we represent the asymptotic state of the two-bump initial data constructed explicitly for small times in Section 4.4. We recall that why the computations in Section 4.4 are only valid for small time is that the first bump reaches the second bump, and we did not compute the first shock after this happens. However, as we see in Figure 7, the first bump "eats" the second bump (possibly in finite time), and we recover the vortex. Since u₀(0) = 0, we have that u(t, 0) = 0 so the vortex cannot be reached in the supremum norm. Notice also that if u₀(0) ≠ 0, then u(0,t)=(u0(0)−α+αt)−1α . Nevertheless, the simulation suggest convergence in all L^p norms for 1 ≤ p < ∞.

Fig. 7

Asymptotic behaviour of u in rescaled variables. See a movie simulation in the supplementary material [1, Video 3].

It is an open problem to determine if the first singularity catches the boundary front in finite time for these particular solutions.

6.2 Comparison of the waiting time

It is interesting to compare the behaviour of different powers u0(ρ)=(β+1)(1−ρ)+β which have total mass M = 1. There is waiting time if β≥1α−1 (see Corollary 4.9 and Theorem 4.10). We will work with α = 2. Since the masses are ordered, the waiting time for u₀ = 2(1 − ρ)₊ is shorter than that of u0=3(1−ρ)+2. . It is interesting to notice that the solution for u0=3(1−ρ)+2 develops a singularity at the interior of the support, before the support starts moving.

6.3 Level sets of a solution with and without waiting time

In Section 4.6 we showed that t1α is the asymptotic behaviour of the support of u for compactly supported u₀. For instance, if u₀ is a Dirac δ function at S(0) of mass M, we have shown that the support is [S(0),S(0)+M(αt)1α] . However, for solutions with waiting time, we do not know what is the behaviour of the support for t small. We illustrate an example when u₀ = (1 − ρ)₊ for α = 2 in Figure 9 (cf. Figure 6). This initial datum produces a solution with waiting time due to Corollary 4.9, which by Theorem 3.4 is initially given by the generalised characteristics. However, as pointed out in Remark 3.1 the characteristics are not the level sets of m. Notice that the level sets of m are not straight even for t small. For comparison, we show a solution not given by characteristics (therefore not a classical solution) and without waiting time (by Theorem 4.10) which we represent in Figure 10.

Fig. 8

Behaviour of two different powers with waiting time. See a movie simulation in the supplementary material [1, Video 4].

Fig. 9

Level sets of the numerical solution with u₀ = (1 − ρ)₊ for α = 2, and a uniform mesh in space of equispaced grid h_ρ = 1e − 3. In Figure 6 the reader may find a comparison with the mass subsolution with explicit Ansatz.

Remarks and open problems

We have constructed a theory of radial solutions and proved well-posedness of the mass formulation. Uniqueness in terms of the u variable is an open problem.
Is there a non-radial theory? This seems to be a very difficult problem.
Is there asymptotic convergence to the vortex solution in the u variable in general?
An interesting problem is to construct a theory for infinite mass solutions.
In the two bump solution, is there actually convergence to the vortex in finite time? The numerical experiments suggest so. The ODE for S₁ can be written explicitly from the Rankine-Hugoniot condition, and the question is whether S₁(t) = S₂(t) for some t > 0.

Fig. 10

Level sets of the numerical solution not given by characteristics presented in Remark 3.8, which does not have waiting time due to Theorem 4.10.

Acknowledgments

JAC was partially supported by EPSRC grant number EP/P031587/1. The research of JAC and DGC was supported by the Advanced Grant Nonlocal-CPD (Nonlocal PDEs for Complex Particle Dynamics: Phase Transitions, Patterns and Synchronization) of the European Research Council Executive Agency (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 883363). The research of DGC and JLV was partially supported by grant PGC2018-098440-B-I00 from the the Ministerio de Ciencia, Innovación y Universidades of the Spanish Government. JAC, DGC, and JLV were partially supported by Grant RED2018-102650-T funded by MCIN/AEI/10.13039/501100011033. JLV was an Honorary Professor at Univ. Complutense.

Conflict of interest statement.
Authors state no conflict of interest.

References

[1] Videos supplementary material. https://figshare.com/projects/Vortex_formation_for_a_non-local_interaction_model_with_Newtonian_repulsion_and_superlinear_mobility/83687.Search in Google Scholar

[2] Yves Achdou, Guy Barles, Hitoshi Ishii, and Grigory L. Litvinov. Hamilton-Jacobi Equations: Approximations, Numerical Analysis and Applications, volume 2074 of Lecture Notes in Mathematics. Springer Berlin Heidelberg, Berlin, Heidelberg, 2013.10.1007/978-3-642-36433-4Search in Google Scholar

[3] Luigi Ambrosio, Nicola Gigli, and Giuseppe Savaré. Gradient flows in metric spaces and in the space of probability measures. Lectures in Mathematics ETH Zürich. Birkhäuser Verlag, Basel, second edition, 2008.Search in Google Scholar

[4] Luigi Ambrosio, Edoardo Mainini, and Sylvia Serfaty. Gradient flow of the Chapman-Rubinstein-Schatzman model for signed vortices. Annales de l’Institut Henri Poincare (C) Non Linear Analysis, 28(2):217–246, 2011.10.1016/j.anihpc.2010.11.006Search in Google Scholar

[5] Luigi Ambrosio and Sylvia Serfaty. A gradient flow approach to an evolution problem arising in superconductivity. Communications on Pure and Applied Mathematics, 61(11):1495–1539, 2008.10.1002/cpa.20223Search in Google Scholar

[6] Philippe Bénilan and Stanislzav Kružkov. Conservation laws with continuous flux functions. Nonlinear Differential Equations and Applications, 3(4):395–419, 1996.10.1007/BF01193828Search in Google Scholar

[7] Andrea L. Bertozzi, Thomas Laurent, and Flavien Léger. Aggregation and spreading via the Newtonian potential: the dynamics of patch solutions. Math. Models Methods Appl. Sci., 22(suppl. 1):1140005, 39, 2012.10.1142/S0218202511400057Search in Google Scholar

[8] J. A. Carrillo, S. Lisini, G. Savaré, and D. Slepčev. Nonlinear mobility continuity equations and generalized displacement convexity. J. Funct. Anal., 258(4):1273–1309, 2010.10.1016/j.jfa.2009.10.016Search in Google Scholar

[9] José Carrillo. Entropy solutions for nonlinear degenerate problems. Archive for Rational Mechanics and Analysis, 147(4):269–361, 1999.10.1007/s002050050152Search in Google Scholar

[10] José A. Carrillo, David Gómez-Castro, and Juan Luis Vázquez. A fast regularisation of a Newtonian vortex equation. 2022. To appear in Annales de l’Institut Henri Poincaré C, Analyse non linéaire.10.4171/aihpc/17Search in Google Scholar

[11] M. G. Crandall and P. L. Lions. Two Approximations of Solutions of Hamilton-Jacobi Equations. Mathematics of Computation, 43(167), 1984.10.1090/S0025-5718-1984-0744921-8Search in Google Scholar

[12] Michael G Crandall and Pierre Louis Lions. Viscosity solutions of Hamilton-Jacobi equations. Transactions of the American Mathematical Society, 277(1):1–1, 1983.10.1016/S0304-0208(08)71044-2Search in Google Scholar

[13] Jean Dolbeault, Bruno Nazaret, and Giuseppe Savaré. A new class of transport distances between measures. Calc. Var. Partial Differential Equations, 34(2):193–231, 2009.10.1007/s00526-008-0182-5Search in Google Scholar

[14] Miguel Escobedo, Juan Luis Vazquez, and Enrike Zuazua. Asymptotic behaviour and source-type solutions for a diffusion-convection equation. Archive for Rational Mechanics and Analysis, 124(1):43–65, 1993.10.1007/BF00392203Search in Google Scholar

[15] S N Kružkov. First Order Quasilinear Equations in Several Independent Variables. Mathematics of the USSR-Sbornik, 10(2):217–243, 1970.10.1070/SM1970v010n02ABEH002156Search in Google Scholar

[16] Fang-Hua Lin and Ping Zhang. On the hydrodynamic limit of Ginzburg-Landau vortices. Discrete and Continuous Dynamical Systems, 6(1):121–142, 2000.10.3934/dcds.2000.6.121Search in Google Scholar

[17] Pierre-Louis Lions, Panagiotis Souganidis, and Juan Vázquez. The Relation Between the Porous Medium and the Eikonal Equations in Several Space Dimensions. Revista Matemática Iberoamericana, pages 275–310, 1987.10.4171/RMI/51Search in Google Scholar

[18] Sylvia Serfaty and Juan Luis Vázquez. A mean field equation as limit of nonlinear diffusions with fractional Laplacian operators. Calculus of Variations and Partial Differential Equations, 49(3–4):1091–1120, 2014.10.1007/s00526-013-0613-9Search in Google Scholar

[19] Hung Vinh Tran. Hamilton-Jacobi equations: viscosity solutions and applications. Lecture notes available from author at http://www.math.wisc.edu/~hung/lectures.html. Accessed: 2019-10-21.Search in Google Scholar

[20] Juan Luis Vázquez. The interfaces of one-dimensional flows in porous media. Transactions of the American Mathematical Society, 285(2):717–737, feb 1984.10.21236/ADA132862Search in Google Scholar

Received: 2021-04-28

Accepted: 2022-01-04

Published Online: 2022-02-25

This work is licensed under the Creative Commons Attribution 4.0 International License.