A mixed system modeling two-directional pedestrian flows

. In this article, we present a simpliﬁed model to describe the dy-namics of two groups of pedestrians moving in opposite directions in a corridor. The model consists of a 2 × 2 system of conservation laws of mixed hyperbolic-elliptic type. We study the basic properties of the system to understand why and how bounded oscillations in numerical simulations arise. We show that Lax-Friedrichs scheme ensures the invariance of the domain and we investigate the existence of measure-valued solutions as limit of a subsequence of approximate solutions.


1.
Introduction. Systems of conservation laws that are not everywhere hyperbolic in the phase space arise naturally in the modeling of physical phenomena. Two wellknown examples are the two-fluid single-pressure model for two-phase flow [22], and a model for three-phase porous medium flow that has been widely used in petroleum reservoir simulation [1]. Another model arises in modeling two-directional traffic flows [3,4]. These models display an elliptic region in the phase space. An elliptic state in the solution space is a state where the Jacobian matrix of the vectorvalued flux function has complex eigenvalues. The set of all elliptic points forms the elliptic region. This type of systems has been addressed since several decades now, see [10,17,18] for a general overview. Nevertheless, their solutions have not been completely understood yet. Rapidly oscillating but bounded numerical approximations suggest that solutions could be defined in the framework of Young measures [13].
In this article, we consider a mixed type system of conservation laws describing two populations of pedestrians moving in opposite directions, adapted from [2,4].
be the governing equations, together with the initial conditions where t ∈ R + , x ∈ R. The flux function is therefore given by where f (u, v) = u(1 − u − v), and U = (u, v) are densities of the two groups of pedestrians that take values in As announced, the main feature of system (1) lies in the hyperbolicity loss for certain density values. Indeed, the Jacobian of the flux exhibits complex eigenvalues in an elliptic region of the phase space Ω. It was suggested in [3] that oscillations arising in the elliptic region could be related to the lane formation phenomenon observed in groups of pedestrians moving in opposite directions [14,21]. We aim to investigate solutions properties in relation with the modeled pedestrian dynamics, also relying on numerical simulations. In particular, we will study the solutions of (1), (2) corresponding Riemann-like initial data U 0 = (u 0 , v 0 ) of the form A similar problem was addressed by Vinod [24], who considered a slightly different version of model (1) including a parameter β > 0, which sensibly changes the solutions behavior. The paper is organized as follows. Section 2 contains the basic analytical study of the models properties. In Section 3, we introduce a Lax-Friedrichs finite volume scheme and we prove an L ∞ bound on the corresponding approximate solutions, ensuring the convergence towards Young measures. In Section 4 we give examples of weak solutions in distributional sense. The study two examples of initial data generating persistent oscillations is presented in Section 5 and we give conclusions in Section 6.
2. Basic analytical study. This section is devoted to the study of the basic properties of system (1), and to the identification of the wave types appearing in the solutions of the corresponding Riemann problem (1)-(4). This study does not pretend to be exhaustive, the problem being non classical and still not completely understood. In particular, we cannot give any global existence result for weak solutions, and their uniqueness is not expected. Some examples of solutions displaying the described features will be showed through numerical computations in Section 4.
First of all, we compute the Jacobian of the flux (3): and its characteristic polynomial The discriminant of (5) is MIXED SYSTEM FOR TWO-DIRECTIONAL PEDESTRIAN FLOWS 3 Figure 1. The domain Ω and the elliptic region E delimited by the curve 4 + 14uv − 12u − 12v + 9u 2 + 9v 2 = 0.
which is negative in the region Fig. 1. Therefore the eigenvalues can take complex values : and corresponding eigenvectors are :

PAOLA GOATIN AND MATTHIAS MIMAULT
The gradient of the eigenvalues are: For seek of clarity, we recall here the notion of distributional solution.
Weak solutions in distributional sense, when they exist, consist of a combination of rarefactions and shock waves, which are described below.
where V (ξ) satisfieṡ Let R i (U L ) denote the solution of (8). The integral curves of r i (U ) are illustrated on Fig. 2. The arrows indicate directions of increase of the corresponding eigenvalue. The direction is reversed across the fognals F i , see Sec. 2.3.
Vector fields r i (U ), i = 1, 2, oriented so that ∇λ i · r i (U ) ≥ 0, are given in Fig. 3, which give an idea of the orientation of rarefaction curves.

Shocks. The discontinuous function
is a weak solution of (1)-(4) called shock wave if and only if it satisfies the Rankine-Hugoniot relation for some speed s = s(U L , U R ) ∈ R. Given any point U L ∈ Ω, the Hugoniot locus [5] is defined by A Hugoniot locus is composed of three disjoint branches. When U L lies in the hyperbolic region, one of them has a loop closing at U L , which crosses the elliptic region. Another branch can also cross a different hyperbolic region but the three branches can never be in Ω at the same time. Those informations are visible on Fig. 4. Since the system is not genuinely non-linear (and not even hyperbolic), the  . Eigenvectors fields of r 1 (left) and r 2 (right) outside the elliptic region E, oriented so that ∇λ i · r i ≥ 0. entropy admissible branch can be selected using the Liu-Oleinik condition [20]: for each U ∈ H(U L ) between U L and U R Anyway, each branch section of H(U L ) belonging to Ω lies in a region were the corresponding field is genuinely non-linear, as shown in Fig. 4. Therefore, if U R belongs to the same branch than U L , condition (11) coincides with the usual Lax 2.3. Fognals and umbilic points. Due to the mixed character of system (1), weak solutions can display other type of discontinuities. Combinations of contact discontinuities moving with zero speed appear along fognals : We define (13) F coincides with the line u + v = 1 and crosses ∂E and ∂Ω at points called umbilic points. One can see orientation's change on the eigenvector field in Fig. 3 and observe that the wave type on the u and the v-axis changes at points 2 3 , 0 , 0, 2 3 and 1 2 , 1 2 . In particular, we will call crossing shocks the discontinuities satisfying see [16].
MIXED SYSTEM FOR TWO-DIRECTIONAL PEDESTRIAN FLOWS 7 3. Numerical scheme. We take a space step ∆x and a time step ∆t subject to a CFL condition which will be specified later in (17). For j ∈ Z and n ∈ N, let x j+1/2 = j∆x be the cells interfaces, x j = (j−1/2)∆x the cells centers and t n = n∆t the time mesh. We want to construct a finite volume approximate solution of (1)- [. We use the following Lax-Friedrichs scheme: with the numerical flux F = (F 1 , F 2 ) defined by for α ≥ 1. We prove by induction that the domain Ω is invariant for (15)-(16) for any initial data U 0 ∈ Ω the approximate solutions computed by scheme (15)-(16) satisfy the following uniform bounds: Proof. We proceed by induction: assuming that u n j ≥ 0, v n j ≥ 0 and u n j + v n j ≤ 1 for all j ∈ Z, we show that the same holds for u n+1 j and v n+1 j .
To prove positiveness, we focus on the u component, the procedure being similar for v. Dropping the index n, we compute By

PAOLA GOATIN AND MATTHIAS MIMAULT
and we compute (dropping again the index n)

Using the hypothesis that
therefore concluding the proof.
The uniform L ∞ -bound provided by Lemma 3.1 ensures the convergence towards Young measures, which are weak- * measurable maps ν : R + × R → P(R 2 ), where P(R 2 ) denote the space of probability measures on R 2 .
Frid and Liu [13] provide an explicit formula for computing the probability measure ν satisfying (18) in the case of Riemann-type initial data (4). For any space step ∆x fixed, let U k = U ∆x/k , k ∈ N, be the sequence of approximate solutions obtained dividing ∆x by k. Then, for any h ∈ C(R 2 ; R) it holds for almost every x/t ∈ R. (For the proof, see [13,Appendix A.2].) For numerical purposes, we will use the following discretized version of formula (19) to compute the moments of interest: where [·] denotes the integer part. Relying on Young measures, DiPerna [8] introduced the concept of measurevalued solutions. Definition 3.3. Let P(R 2 ) denote the space of probability measures on R 2 . A measure-valued solution of (1) is a measurable map ν : In order to prove that the Young measure obtained as limit of finite volume approximations is indeed a measure-valued solution, one needs some weak BV estimates that allow to pass to the limit in the discrete conservation equation, replacing the strong convergence argument used in the proof of Lax-Wendroff theorem [19]. Weak BV estimates have been derived in the case of monotone schemes for scalar conservation laws in several space-dimensions [7,9], but also for entropy stable methods applied to hyperbolic systems [11,12]. Unfortunately, in the present case the system's entropy is not convex, and these techniques don't apply.
4. Distributional solutions. In this section, we present some numerical computations illustrating the principal features of weak (distributional) solutions of problem (1)-(4). From [15] we know that, if U L , U R ∈ Ω \Ē, then corresponding weak solutions must satisfy U (x, t) ∈ Ω \ E a.e. In particular, U consists of a combination of rarefactions and shock waves, as illustrated by the following numerical tests, where we have taken ∆x = 0.001, α = 1, CF L = 0.9, and we display the approximate solution U ∆x at time t = 1, both in the x-u, v plane (Figures 5-9 and 11-15, left) and in the phase-plane u-v ( Figures 5-9 and 11-15, right). Test 1. We consider initial data U L = (0.2, 0.1) and U R = (0.1, 0.2). The solution showed in Fig. 5 consists of a rarefaction of the first family joining U L with an intermediate state on ∂E, followed by a rarefaction of the second family to U R . This is a limit situation, since if the Lax curves do not intersect in the same connected region, the structure of the solution becomes more complex, as illustrated by the following examples.
Test 2. We consider initial data U L = (0.2, 0.1) and U R = (0.1, 0.3). The solution showed in Fig. 6 consists of a shock of the first family joining U L with a state U 1 = (u 1 , 0) on the u-axis, followed by a crossing shock between U 1 and the state (1, 0), a contact discontinuity form (1, 0) to (0, 1) with zero speed, another crossing shock from (0, 1) to a point U 2 = (0, v 2 ) and a 2-shock from U 2 to U R .
Test 3. We consider initial data U L = (0.2, 0.1) and U R = (0.1, 0.8). The solution showed in Fig. 7 consists of a shock of the first family joining U L with a state U 1 = (u 1 , 0) on the u-axis, followed by a crossing shock between U 1 and the state (1, 0), a standing contact discontinuity form (1, 0) to a point U 2 = (u 2 , 1 − u 2 ) ∈ ∂Ω and a 2-shock from U 2 to U R .
Test 4. We consider initial data U L = (0.2, 0.1) and U R = (0.85, 0.1). The solution showed in Fig. 8 consists of a shock of the first family joining U L with a   We observe a 1-shock, a crossing shock, a contact discontinuity and a 2-shock. state U 1 = (u 1 , 0) on the u-axis, followed by a crossing shock between U 1 and a state U 2 = (u 2 , 0) and a 2-shock from U 2 to U R . We remark that the crossing shock is sharply captured. Test 5. We consider initial data U L = (0.2, 0.1) and U R = (0.75, 0.1). The solution showed in Fig. 9 consists of a shock of the first family joining U L with a state U 1 in the interior of the domain, followed by a crossing shock between U 1 and a state U 2 superposed to a 2-shock from U 2 to U R . Note that this composite wave could be replaced by a 2-shock joining directly U 1 to U R , see [16]. Indeed, we have U 2 ∈ H(U 1 ), U R ∈ H(U 2 ), U R ∈ H(U 1 ), and s(U 1 , U 2 ) = s(U 2 , U R ) = s(U 1 , U R ). Anyway, the two solutions are identical as L 1 functions. Figure 9. Solution to (1)-(4) with U L = (0.2, 0.1) and U R = (0.75, 0.1). We observe a 1-shock, a crossing shock and a 2-shock. Remark 1. The configurations displayed in Tests 2-3 are unrealistic from the modeling point of view, because they result in a complete blocking of one or both groups of pedestrians, represented by the vacuum regions delimited by the standing contact discontinuities. In reality, such stuck situations never occur in normal conditions, and the flows always organize so that few people manage to pass, even if the resulting capacity can be very reduced [6]. System (1) must be seen as a toy model, whose understanding can give some insight for more realistic approaches.
If one or both values of the Riemann initial data U 0 = (U L , U R ) belong to E, we can still observe distributional solutions in some cases. In accordance to [15], since E is convex, if U L ∈ E or U R ∈ E and U is a weak solutions, then if U (t, x) ∈ E for some t, x, than U (t, x) ∈ {U L , U R }. Indeed, the initial state belonging to E will be connected through a shock to some U M ∈ Ω \ E. In particular, given any point Fig. 10. To get further information on the Young measures associated to the above initial data, we compute their means (or expected values) and variance relying on formula (20). In particular, we get  and U R = (0.1, 0.6) ∈ E. We observe a 1-shock, a crossing shock and a contact discontinuity in the hyperbolic region, followed by a shock joining directly U R (the last state in the hyperbolic region belongs to H((1, 0)) ∩ H(U R )).
Aiming at giving an estimation of the accuracy of the results, we also compute the conservation errors and the convergence rates, which are given in Tables 1, 2. The conservation error is given by the formula The L 1 -convergence rate is defined by where the L 1 -error is computed at final time t = 1 as .
(In the above expressions, the integrals are intended component by component.)    Table 2. Convergence rates (26) at time t = 1.
We observe that conservation errors remain quite stable, while the L 1 -convergence rate is lower than one. Indeed, the very long computing times prevented us to reach sharper approximations. 6. Conclusions. The 2 × 2 system (1) is a simplified model for the motion of two groups of people walking in opposite directions along a corridor. This situation is known for displaying characteristic patterns as lane formation [14,21]. The system consists of two conservation laws of mixed hyperbolic-elliptic type. Though the solution configurations are not always realistic, the instabilities observed for densities in the elliptic region could be related to these auto-organization phenomena that result in a transition from a mixture to separate phases.
Following [13], we suggest that the corresponding Cauchy problem should be recast in the framework of Young measures. This allows to compute the expected values and the variance for densities and fluxes, and to estimate the flow characteristics.