Two algorithms for a fully coupled and consistently macroscopic PDE-ODE system modeling a moving bottleneck on a road

In this paper we propose two numerical algorithms to solve a coupled PDE-ODE system which models a slow vehicle (bottleneck) moving on a road together with other cars. The resulting system is fully coupled because the dynamics of the slow vehicle depends on the density of cars and, at the same time, it causes a capacity drop in the road, thus limiting the car flux. The first algorithm, based on the Wave Front Tracking method, is suitable for theoretical investigations and convergence results. The second one, based on the Godunov scheme, is used for numerical simulations. The case of multiple bottlenecks is also investigated.


Introduction
In this paper we focus on a system modeling a slow vehicle (e.g. large truck or bus) which moves in a car traffic flow. We assume that the slow vehicle's velocity is influenced by the car density and, at the same time, it causes a capacity drop in the road at the position where it is located. Therefore we use the terminology moving bottleneck.
Mathematically, the problem gives rise to a coupled PDE-ODE system. Car density ρ solves the TWO ALGORITHMS FOR A PDE-ODE SYSTEM MODELING A MOVING BOTTLENECK 2 following conservation law with t ∈ [0, T ], ρ ∈ [0, ρ max ], (t, x) ∈ R + × R, f : [0, ρ max ] × R → R the flux function,ρ the initial datum and y(t) the position of the moving bottleneck, which is assumed to solve the following Cauchy problem ẏ(t) = w(ρ(t, y)), with y 0 the position at initial time t = 0 and w ≥ 0 the velocity of the moving bottleneck. An example of f and w will be given in the next Section, see equations (2.1)-(2.2). It is well known that the function ρ(t, ·) may be discontinuous, then solutions to both (1.1) and (1.2) are to be intended in generalized sense. The case of multiple moving bottlenecks will be also considered.
Relevant literature. Let us mention a few results related to our work. Modelling of vehicular traffic with fluid dynamic approach started with the seminal works of Lighthill and Whitham [38] and Richards [41]. The attention on this subject continued in the last years, both for the proposal of more accurate models, see, e.g., [3,16,18], and for the treatment of complex networks, see, e.g., [12,13,29,30] and the book [24]. To deal with the case of different vehicles, such as motorcycles, cars, buses, multi-population models were also considered, see [4,43]. The problem of determining the trajectory of a single car given the density of cars was addressed in [19], and the relative numerics in [11], taking advantage of vector techniques. In this case the observed car is assumed not to influence significantly the traffic flow, thus the system is not fully coupled. Here we consider the more complicated situation in which the position of the single vehicle influences the whole traffic flow.
The coupled system (1.1)-(1.2) was proposed in [32]. More precisely, that paper introduced a notion of solution based on weak solutions for the PDE and Filippov solutions for the ODE, and proved existence of such a solution under suitable assumptions.
Recently, another fully coupled system for a moving bottleneck problem was studied in [20,21,22,37] from both the theoretical and numerical point of view (see also [42] for an extension to the second order ARZ model) . Furthermore, in [40] an optimal control problem using the maximal speed of the coordinated vehicle as control variable is stated. It is worth to point out the difference with the model we consider here: the model in [20,21], following the lines of [39], assumes that the moving bottleneck reduces the maximal density attainable in the road, together with its maximal flux (i.e. the capacity of the road). Moreover, it assumes that the moving bottleneck moves at constant speed for low car density. Instead, in the model we consider here the moving bottleneck does not reduce the maximal density attainable in the road, while it only reduces the maximal flux. Moreover, its velocity is always dependent on the car density. These features make the model fully consistent with a pure macroscopic point of view.
Let us also recall that the problem of modeling bottlenecks, moving or not, was addressed in [1,17,23] and the case of multiple moving bottlenecks was already considered in [33].
Finally, the terminology "moving bottleneck" is known also in transport engineering literature: see, among others, [26,34,36,39]. In these papers the problem of moving bottlenecks, both with a given 3 path [26,36,39] and with a mutual interaction with the surrounding flow [34], is studied, but without a detailed description of the resulting schemes and of their convergence.
Goal. We propose two different numerical algorithms, the first allowing theoretical investigations, while the second allowing an actual implementation on a computer.
The first algorithm, called WFT-ODE, combines the Wave Front Tracking (WFT) method [6,31] to solve (1.1) with an exact method to solve (1.2). The moving bottleneck position is traced taking into consideration interactions with (shock or rarefaction) waves. Focusing on bounded variation data, we establish a convergence result for the approximate solution towards a solution to (1.1)-(1.2). We are also able to provide a convergence rate for y in terms of total variation of initial datum of density, namely T V(ρ). To achieve the estimates, we measure the distance among successive approximations by means of tangent vectors. The latter represents distances among discontinuities (shocks or rarefaction shocks) for car density, and the difference of positions for moving bottleneck.
The second algorithm, called GOF, is a fractional step method combining the classical Godunov method [28] for (1.1) with an exact method for (1.2). There are two reasons for using the Godunov method: the easy implementation and the importance of this method in the transportation literature [35,36].
The GOF method is then extended to the case of multiple bottlenecks. For that, we consider two models. The first allows slow vehicles to pass each other, while the second inhibits such a possibility. The latter is more challenging from theoretical and numerical point of view, thus we focus on this case. Both methods are useful to deal, for instance, with the case of many buses on the same urban route [25].
Paper organization. In Section 2 we provide basic definitions and preliminary results. The theoretical scheme WFT-ODE is introduced in Section 3, while convergence results are given in Section 4. The second scheme GOF is described in Section 5. Then we deal with multiple bottlenecks in Section 6. Finally, numerical simulations are reported in Section 7.

Basic definitions and preliminary results
In this section we provide the definition of solution to the system (1.1)-(1.2) and we detail the choice of the flux f and the velocity w. We also give a description of Wave Front Tracking (WFT) algorithm for the case of a flux function depending on the space variable, as it occurs for (1.1). Finally, a result on possible interactions of the moving bottleneck with fronts in the WFT solution is proved. The latter will be extensively used in the rest of the paper.
We denote by BV(R) the space of functions with bounded variation. For every T > 0, the definition of solution to the system (1.1)-(1.2) reads as follows.

] is defined as the smallest interval containing a and b.
We assume hereafter that the maximal density of cars is ρ max =1. We denote by w(ρ) the velocity of the moving bottleneck and we assume that for some constant w max > 0. We denote by v(ρ, x − y) the velocity of cars and we assume that where ϕ ∈ C 1 (R; (0, +∞)) and it is strictly decreasing in (−β, 0] and strictly increasing in [0, β) for some β > 0. We also assume that there exist two positive constants v,v such that and that v is reached for ζ = 0 andv is reached for any ζ ∈ R\[−β, β], see Fig. 1(a). Note that, in particular ϕ ′ (ζ) is compactly supported and, for any fixed density ρ, the minimal velocity of cars v is taken when they overtake the moving bottleneck, i.e. for x = y. We finally assume v > w max .
The last assumption is crucial for the model: it means that cars can overtake the bottleneck, see Fig.  1

Preliminaries for Wave-Front Tracking algorithms
As it is well known, a solution to the Cauchy problem of a conservation law can be constructed using the WFT algorithm. Roughly speaking, first a BV initial datum is approximated by a piecewise constant function via sampling (thus with a smaller BV norm). For small times, a piecewise constant weak solution is obtained by piecing together solutions to Riemann problems, where rarefactions are replaced by a fan of rarefaction shocks. Then, when waves interact, a new Riemann problem is generated and solved again approximating rarefactions by rarefaction shocks' fans and so on. For details we refer the reader to [6] for the general theory and to [24] for the case of networks.
In our case the flux does not depend only on ρ, but also on x and t (via y(t)), so the WFT algorithm needs to be modified. Moreover, we will evolve, at the same time, the position of the bottleneck y(t) according to equation (1.2). In order to perform the construction we first need a preliminary result comparing the speed of the bottleneck with that of neighboring waves of the WFT approximation of ρ.
Proof. Recall that the discontinuity travels with speed: We divide the proof in two cases: 1. The wave is a shock wave. 2. The wave is a rarefaction-shock wave.
2. In the case of rarefaction-shock we have ρ r = ρ ℓ − δ ν , where 0 < δ ν < 1 ν is a parameter that will be defined by the WFT algorithm. We only prove the inequality in (2.9), then (2.10) follows because ρ ℓ > ρ r . The inequality in (2.9) can be rewritten as Reasoning as above, we have Now the worst case occurs when ρ r = ρ min , ρ ℓ = ρ min + δ ν : Then we conclude as in step 1.
The WFT algorithm will need to be further split because of source effect due to the x dependence of the flux. Indeed, consider a general equation of the type where we assume (H1) g is Lipschitz continuous in both variables; for every ρ, g(ρ, ·) ∈ C 2 (R); max ρ∈[0,1] R |∂ xx g(ρ, x)|dx < +∞; for every x, g(·, x) is concave.
Therefore, formally deriving in space the flux function g depending both on the density of cars ρ(t, x) and on the position of cars x, one gets: so that equation (2.13) can be written as a conservation law with source term In our case we have and the source term reads (see Fig. 2). Note that the assumptions (H1) are satisfied and that the x dependence in the flux is indeed effective only for ζ ∈ (−β, β), that is when cars are close to the bottleneck's position. Equation (2.13) can be solved applying an operator splitting method in which one alternates between solving the homogeneous conservation law and the ordinary differential equation

The semi-discrete WFT-ODE scheme and its properties
We are now ready to define the WFT-ODE scheme. In this section we describe a semi-discrete scheme for the coupled dynamics (1.1)-(1.2). The idea is to combine the WFT method for the conservation laws with exact solution to the ODE of the moving bottleneck. Indeed, since the approximate solution given by WFT is piecewise constant (in time and space), we can get a simple ODE with piecewise constant (in time) right-hand side for the moving bottleneck.
Let us now describe in detail the WFT-ODE scheme. Fix T > 0, then for every ν ∈ N and ∆t > 0, we define an approximate solution (ρ ν , y ν ) on the time interval [0, T ] via the following steps.
Step 2. At time t = 0, we solve all Riemann problems at each discontinuity point ofρ ν for the homogeneous equation (2.16) and approximate every rarefaction wave with a rarefaction fan, formed by TWO ALGORITHMS FOR A PDE-ODE SYSTEM MODELING A MOVING BOTTLENECK 9 rarefaction shocks (non-entropic shocks) of strength less than δ ν travelling with speed given by the Rankine-Hugoniot condition. More precisely: Assume, for simplicity, that y 0 ∈]s j , s j+1 [ for some j (i.e. y 0 s i for every i = 1, . . . , N.) For small times we can get an approximate solution by solving the equations: for the discontinuities s i of ρ ν and the equation: for the bottleneck, which gives y(t) = y 0 + tw(ρ(t, s j +)).
If two wave fronts s i and s i+1 interacts, then we solve the new Riemann problem. Because of and so on for all the possible interactions with ρ ν waves up to time ∆t. Finally, possibly slightly modifying the wave speeds, we can assume that y ν (∆t) s i (∆t) for any j.
Step 3. To take into account the source term ρ ν (∆t, x) is updated in order to include the effect of the source term. This is done by means of one step of the explicit Euler method approximating the ODE (2.17) by: After this modification the function ρ source ν is no longer piecewise constant and we need to apply a specific sampling to ensure accurate estimates of continuous dependence.
Step 4. Letȳ = y(∆t). We update ρ source ν by sampling only on the interval First we insert the points y j =ȳ − β + j 2β 2 ν , j = 1, . . . , 2 ν . Moreover, if there exists i such that no y j is in the interval ]x i , x i+1 [, then we insert the additional point y i = x i +x i+1 2 . Again possibly slightly modifying the wave speeds we can assume that x i y j for every i and j. We define ρ sampled by sampling ρ source on subintervals generated by the partition: The rationale behind this definition is the following. To obtain neat estimates of increase in jumps size of shocks and rarefaction shocks we need always to sample to the left and right of the discontinuities which were generated by the WFT step. In the sequel, we shall refer to these new generated (at positive time) waves as s b waves, being generated by the source term due to the bottleneck.
We can restart the procedure at step 2 in the interval [∆t, 2∆t] and so on.

TWO ALGORITHMS FOR A PDE-ODE SYSTEM MODELING A MOVING BOTTLENECK 10
It is important to note that, while the sampling procedure does not increase the total variation, the updating step (3.1) can do it. For the first time step, the increment can be estimated easily as follows. Let us denote by ω(x; ρ ν , ∆t) the contribution given by the source term at time ∆t, i.e.
Recalling the assumption (H1), we have Then, at every time step the additional total variation is bounded by max ρ ω ′ (·; ρ, ∆t) L 1 , and at final time T , namely after T ∆t time steps, the total increment is bounded by T max ρ ∂ xx g(ρ, ·) L 1 . More precisely we have the following: Bounds on the BV norm and on number of waves and interactions are directly obtained as in the scalar case, when the flux is not depending on x, provided this dependence does not create resonance effects among waves. Indeed, we can rewrite a scalar conservation law with a x-dependent flux g(ρ, x) as a 2 × 2 system as follows: The above system is strictly hyperbolic, and therefore standard WFT analysis may be carried out, provided ∂ ρ g(ρ, z) 0. In our case, g(ρ, x; y(t)) = ρ(1 − ρ)ϕ(x − y(t)) and the above condition corresponds to the requirement that the bottleneck's position y(t), where the x-dependence is effective as noticed before, travels with a speed uniformly different with the one of the closest ρ waves. This is precisely the statement of Lemma 1. The construction thus relies on Lemma 1, which in turn is applicable if ρ ν >¯v −w max 2v−w max , see (2.5). Due to the source term −∂ x g this is not guaranteed for all times even if the initial data satisfies (2.5). However we have the following: Then for all positive times T such that: the condition (2.5) of Lemma 1 is verified for ∆t sufficiently small.
Notice that forv → w max we have η → 0, thus T → ∞. In other words, if the maximal speed of cars is close the that of the moving bottleneck then the WFT-ODE scheme is defined for large times.
Proof. Notice that solutions to the ODE linked to the source term (2.17) . This gives the estimate (3.3). Since Step 3 is an approximation of the source term, the estimate remains valid for ∆t sufficiently small.

Convergence of the WFT-ODE scheme
In this section we prove the convergence of the WFT-ODE scheme. The BV estimates on the WFT scheme allows to pass to the limit by standard arguments. On the other side, the convergence of the the bottleneck position requires additional work, in the same spirit of [11]. However, our case is more complex due to the complete coupling of the PDE for car density and the ODE for the bottleneck position. As in [11], to achieve this goal we use the technique of generalized tangent vectors [7,8], to the WFT approximate solutions ρ ν and to the bottleneck position y ν .
More precisely, after introducing the tangent vector technique, we consider ρ ν+1 as obtained from ρ ν by shifts of discontinuities, which generate generalized tangent vectors. Then we need to estimate: 1) The increase in the norm of tangent vector to ρ ν due to the WFT algorithm.
2) The increase in the norm of tangent vector to ρ ν due to the source term.
3) The increase in the norm of tangent vector to y ν due to interactions with waves of ρ ν .

Generalized tangent vectors
The technique is based on the idea of considering L 1 as a Finsler manifold with generalized tangent vectors with appropriate norm. We first consider the subspace of piecewise constant functions and "generalized tangent vectors" consisting of two components ((v θ , ξ θ ), η θ ), where ξ θ describes the infinitesimal displacement of discontinuities and the scalar η θ is the infinitesimal shift of the car trajectory. Let us take a family of piecewise constant functions {ρ θ } θ∈[0,1] , each of which has the same number of jumps, say at the points s θ and also the quantities Note that v θ is not defined if x = s θ k , k = 1, . . . , N. Indeed, the contribution of the jumps to the tangent vector is given by {ξ θ k } k , which take into account the presence of the discontinuities. Moreover, v θ solve the linearized equation along the trajectory, see [6] while the shifts change only at interactions times or, as detailed in next remark, due to Step 3 in the WFT-ODE scheme.

Remark 1.
It is worth observing that, as a consequence of the shift η θ of the moving bottleneck's position, the waves s b generated in Step 4 of the scheme will be all shifted by the same quantity.
Then we say that the path γ : θ → (ρ θ , y θ ) admits tangent vectors ((v θ , ξ θ ), η θ ) ∈ T ρ θ=L 1 (R, R) × R N × R. Note that in general such path is not differentiable w.r.t. the usual differential structure of L 1 . One can compute the L 1 -length of the path γ in the following way:

TWO ALGORITHMS FOR A PDE-ODE SYSTEM MODELING A MOVING BOTTLENECK 12
According to (4.1), in order to compute the L 1 -length of a path γ, we integrate the norm of its tangent vector which is defined as follows: where ∆ρ θ k is the jump across the discontinuity s θ k . Notice that this is not the usual L 1 length of a path since there is the additional component η.
Let us introduce the following definition.
Definition 2. We say that a continuous map γ :
(possibly) except for θ = 0, 1, where the number of jumps changes. Finally we have: Lemma 4. The initial dataρ ν andρ ν+1 can be connected by a regular path with shifts whose norm is bounded by 2 −(ν+1) .
Since ρ θ is piecewise constant, v θ (x) is equal to 0 (where it is defined), then v θ L 1 = 0. At time t = 0, we also have η θ = 0 since y ν (0) = y ν+1 (0). By (4.1), we have Now we consider the evolution in time of γ, denoted by γ t , in the sense that we perform the wave front tracking forρ ν andρ ν+1 . It is easy to prove that γ t is still regular (some more point of not differentiability will arise because of interactions order). We aim at proving that for an appropriate constant K > 0. Then uniqueness and Lipschitz continuous dependence of solutions to Cauchy problems is straightforwardly achieved passing to the limit on the WFT-ODE approximate solutions.

Estimates on tangent vectors to ρ ν due to WFT
To estimate the shifts of waves at interactions times between waves, we can use the estimate [24, Lemma 2.7.2] (originally proved in [5]).
Lemma 5. Consider two waves with speeds λ 1 and λ 2 respectively, that interact together producing a wave with speed λ 3 . If the first wave is shifted by (the tangent vector of the shift) ξ 1 and the second wave by (the tangent vector of the shift) ξ 2 , then the (tangent vector of the) shift of the resulting wave at time of interaction is given by Moreover we have

4)
where ∆ρ i are the signed strengths of the corresponding waves.
Proof. The interacting waves satisfy the ODE: From the regularity of ϕ we have that x i are smooth functions of time, thus we can linearize them around the interaction time and apply Lemma 2.7.2 of [24].
The wave shifts change also when the waves are located within the influence zone of the moving bottleneck, more precisely we have:

TWO ALGORITHMS FOR A PDE-ODE SYSTEM MODELING A MOVING BOTTLENECK 14
Proof. Let x(t; z), t ∈ [t 1 , t 2 ] be the unique solution to the ODE: Letx be the position of the wave at time t 1 and consider the solution to shifted initial position x(t;x+ξ), t ∈ [t 1 , t 2 ]. Note that, since λ = ϕ(ζ)∆(ρ(1−ρ)) ∆ρ , the solution x(t; ·) satisfies uniformly in t ∈ [t 1 , t 2 ], with C = ϕ ′ ∞ = max R |ϕ ′ | the uniform Lipschitz constant for λ, and then we also obtain the estimate

Estimates on tangent vectors to ρ ν due to the source
As specified above, our WFT algorithm includes a sampling procedure on the interval [ȳ − β,ȳ + β] where ρ source is not piecewise constant due to the action of the source term; see (3.1). The so-defined piecewise constant approximation ρ updated has new tangent vectors which are estimated in the following lemma.
Proof. Due to Step 4 in our WFT algorithm (cfr. Section 3), a wave (ρ l , ρ r ) of ρ source located at a point in the region of influence of the bottleneck will be located in ρ updated at the same point, that is, ξ − = ξ + , indicating by ± the values before and after the source action and the sampling. However, the jump will be changed by this procedure, namely (ρ l − ∆t ∂ x ϕ ρ l (1 − ρ l ), ρ r − ∆t ∂ x ϕ ρ r (1 − ρ r )). Therefore and in addition

Estimates on tangent vectors to y ν
The tangent vector to y ν , i.e. η, changes only at interaction times of the moving bottleneck with waves. Both the moving bottleneck and the waves satisfy ODEs with smooth right-hand-side, thus by linearization we find the same interaction estimate as for the case treated in [11]. Proposition 1. Let t * be an interaction time between the moving bottleneck and a wave (ρ ℓ , ρ r ) and indicate by η ± the moving bottleneck shift after and before the interaction respectively. Then we have: if the wave is a shock; (4.12) where µ is given by Lemma 1.

Convergence of WFT-ODE
To prove convergence of the WFT-ODE scheme, we estimate the increase of the tangent vectors η to the moving bottleneck position y ν . We fix ∆t and ν and start introducing some notation.

Definition 3. We denote by (T V) p j the total variation of waves of ρ ν ( j∆t) which are to the right of the moving bottleneck (thus can potentially interact with it), and by (T V) i j the total variation of waves of ρ ν which interact with the moving bottleneck in the time interval [ j∆t, ( j + 1)∆t]. Similarly we denote by V p j the sum of the norms of tangent vectors to waves of ρ ν ( j∆t) which are to the right of the moving bottleneck (thus can potentially interact with it), and by V i j the norms of tangent vectors to waves of ρ ν which interact with the moving bottleneck in the time interval [ j∆t, ( j + 1)∆t].
We denote by η j the tangent vector to y ν at time j∆t.
With these notations we have the following estimates. From Lemma 2, we get: (4.14) From Lemmas 2, 5, 6 and 7, and Remark 1, we get: More precisely, Lemma 2 guarantees that the total variation of new waves due to Step 3 is bounded by K ϕ ∆t, thus the new contribution to the tangent vector norm is bounded by K ϕ ∆t η j+1 by Remark 1; Lemma 5 guarantees that the norm of tangent vector does not increase for waves interactions; Lemma 6 shows that the magnification due to permanence in the moving bottleneck influence zone is bounded by e ϕ ′ ∞ ∆t ; Lemma 7 provides the bound (1 + ϕ ′ ∞ ∆t) on magnification due to the source term. Finally, from Proposition 1, η may increase by interacting with a wave (ρ ℓ , ρ r ) having shift ξ with multiplicative factor at most (1+ w max µ |ρ ℓ −ρ r |) for any rarefaction wave and additive factor w max µ |ξ| |ρ ℓ −ρ r | for any wave. The worst case scenario happens when first all shocks interact and then all rarefactions, this is bounded by the estimate:

TWO ALGORITHMS FOR A PDE-ODE SYSTEM MODELING A MOVING BOTTLENECK 16
We can rewrite the estimates (4.14)-(4.16) as: (4.17) where K = max{K ϕ , 2 ϕ ′ ∞ , w max µ }. We now have the following: (4.17). Let µ j be the solution to the system: On the other side this may not achieve the maximal increase in V p j+1 because of the multiplicative term (1 + K∆t). However, such maximal increase of V p j+1 on the interval [0, T ] (thus after all time steps) is bounded by e KT , which is the term appearing on the right-hand side of (4.18) (second equation). Therefore, with this corrected multiplicative term, the increase is maximized. Finally, η 0 = 0 thus the initial data for (4.18) are given by µ 1 = Ke KT V 0 and V 1 = K∆t µ 1 .
We are now ready to estimate η j . First notice that we have: and finally: Thus we obtain the following: Proposition 2. Consider initial conditionsρ ∈ BV(R), y(0) ∈ R and time horizon T > 0. Using the notation of Lemma 3, assume that η < 1 and T satisfy (3.3). Let (ρ ν , y ν ) be the approximate solutions computed by the WFT-ODE scheme, then: where K 1 depends only on the total variation ofρ, T , w max , ϕ ∞ , ϕ ′ ∞ , µ defined in (2.6) and K ϕ defined in (3.2).
Proof. From (4.19), we have that η ν can be estimated in terms of the total variation ofρ, T , the constant K and V 0 . But V 0 is estimated by the total variation ofρ and the initial shifts, thus we conclude by Lemma 4.

Godunov-ODE-FS scheme
Here we introduce a numerical scheme called Godunov-ODE-FS (GOF), which is based on fractional step method combining Godunov scheme for the PDE and exact solution for the ODE. The dynamics of (1.1) and (1.2) are thus solved separately at each iteration. The use of Godunov scheme is motivated by its easy implementation and its connection with the modelling of vehicular traffic problems, see [24,35].

Godunov scheme for the PDE
Following the ideas described in Section 2.1 for the WFT algorithm, we report in the following the modified Godunov scheme for a conservation law with source term (see (2.15)). We first introduce a numerical grid, denoting by ∆x the space mesh size and by ∆t the time mesh size. Moreover ; ρ ℓ , ρ r the self-similar entropy solution of the unique Riemann problem defined on C l m (with discontinuity point x m− 1 2 ) and let us define the numerical flux F as F(ρ ℓ , ρ r , t, x) = g(W(0; ρ ℓ , ρ r ), t, x). First, we replace the initial datumρ(x) by a piecewise constant approximation, Then, we alternate a single step of the classical Godunov scheme with a single step of the Euler scheme ρ l+1 m = ρ * m + ∆t(−∂ x g(ρ * m , t l , x m )).

GOF scheme
Before describing the GOF scheme, we point out that shock waves solutions to (1.1) have velocities depending on the moving bottleneck position y. The modified Godunov scheme can be thus used, provided that, in each discretization cell, the shock speed λ(ρ r , ρ ℓ , ζ) defined in (2.8) does not change sign as a function of ζ. This is established by next Lemma.
The situation described by Lemma 9 is depicted in Fig. 5.
The discretization grid is defined as for the Godunov scheme and the value of density on grid points is called ρ l m . The moving bottleneck position will be defined for times t l and denoted by y l .
Step 0. We approximate the initial datumρ as for the Godunov scheme and set y 0 = y(0).
Step 1. We assume that y l and ρ l m were defined by previous step. We keep fixed the position of the moving bottleneck, thus set y(s) = y l for s ∈ [t l , t l+1 ]. Define W m = W m (t, x) to be the self similar solution to the Riemann problem defined in x m− 1 2 (for the bottleneck position fixed). This solution is formed by a wave, whose velocity varies in time but mantains its sign constant as proved by Lemma 9. We can use Gauss-Green formula as for Godunov scheme and define Step 2. There exists m such that y l ∈ [x m−1− 1 2 , x m− 1 2 ), see Fig. 6. Using the computed density ρ l+1 , the initial velocity of the moving bottleneck is w(ρ l+1 m−1 ) up to the interaction time with the line x = x m− 1 2 then it travels with velocity w(ρ l+1 m ). More precisely, we compute the interaction time as .

Multiple bottlenecks
In this section we introduce a model describing the evolution of car traffic flow along a road in presence of P > 1 moving bottlenecks. In this case, we face the event that two slower vehicles can occupy the same position, which could lead to a non realistic dropping of flux capacity. To solve this problem, we introduce two modelling choices: (A) Two moving bottlenecks can overtake each other, thus we re-define the flux function. (B) Two moving bottlenecks cannot overtake each other. To achieve that, we re-define the velocity function of moving bottlenecks.

Model A
Let us consider the following PDE-ODE coupled system

TWO ALGORITHMS FOR A PDE-ODE SYSTEM MODELING A MOVING BOTTLENECK 20
where y i = y i (t) is the position of the i-th moving bottleneck, and ρ = ρ(t, x) is the density on the road.

Speeds of moving bottlenecks are defined as
Now we admit the presence of many slow-moving vehicles in the car traffic flow, with possibly different characteristics among each other. Therefore, in general, functions w i (·) and ϕ i (·) are different and P can be arbitrarily large. In view of our definition of the function Φ in the flux, when two slow vehicles are far enough, we go back to the case of a single moving bottleneck. On the other side, if two slow vehicles occupy the same position, the flux capacity drop is the largest among the two.

Model B
We now introduce the alternative modelling choice in (B), which avoid overtaking between slow vehicles, again described by a coupled PDE-ODE system with y 1,0 < . . . < y P,0 . The flux function is defined as in (6.1)-(6.2), but this time the function Φ is given by We define differently also the velocities of moving bottlenecks. More precisely, to avoid overtaking among them, and indeed to avoid the superposition of the zones where functions ϕ i take values less than

TWO ALGORITHMS FOR A PDE-ODE SYSTEM MODELING A MOVING BOTTLENECK 21
v, we suitable define the velocity functions w i (ρ, y i , y i+1 ). First let ω i = ω i max (1 − ρ), with 0 < ω i max < v i , and then back recursively on i set (6.5) and complete the definition for y i+1 −y i ∈ [β i+1 +β i , 2(β i+1 +β i )] by smooth monotone interpolation. This choice introduces a follow-the-leader flavour in the microscopic dynamic and avoids the modelling problems described above. A bottleneck at position y i moves with velocity ω i (ρ), depending only on the density of cars along the road, provided the distance from the vehicle ahead in position y i+1 is sufficiently large, namely, larger than 2(β i+1 + β i ). If such distance decreases, then y i starts to decelerate and eventually move with the same velocity as the vehicle ahead when y i+1 − y i = β i+1 + β i .

GOF algorithm for multiple bottlenecks
We introduce a Godunov-ODE-FS algorithm for model (B), while model (A) can be treated in a manner completely similar to the case of a single bottleneck. Notice that for model (B) slow moving vehicles can influence each other, thus the numerics is more involved. The discretization grid is defined as for the Godunov scheme and the value of density on grid points is called ρ l m . The moving bottlenecks' positions will be defined for times t l and denoted by y l i , i = 1, . . . , P. For consistency we suppose Notice that the P-th bottleneck plays the role of the leader and it is not influenced by positions of other slow vehicles. Thus its trajectory can be traced as in Section 5.2. Furthermore, because of definition of w i , the i-th bottleneck's trajectory is only influenced by the (i + 1)-th bottleneck trajectory. More precisely, in the scheme we let the i-th bottleneck proceed with velocity ω i as long as the distance with the (i + 1)-th bottleneck is larger than β i+1 + β i , otherwise we let it proceed with the velocity of the (i + 1)-th bottleneck. Let us first note that, as for the case of a single bottleneck, the velocity of any wave inside each cell change in time, but its sign does not. Indeed, let λ(ρ r , ρ ℓ , x, y) be the speed of the wave, then, using notation of Lemma 9, we have and, since Φ(x, y) > 0, it follows that sgn(λ(ρ r , ρ ℓ , x, y)) = sgn(λ(ρ r , ρ ℓ )).
We are now ready to introduce the GOF algorithm.
Step 0. We approximate the initial datumρ as for the Godunov scheme and set y 0 i = y i (0).
Step 1. We assume that y l i , for i = 1, . . . , P, and ρ l m were defined by previous step. We keep fixed the positions of moving bottlenecks, thus set y i (s) = y l i for s ∈ [t l , t l+1 ]. Let W m = W m (t, x) be the self similar solution to the Riemann problem defined in x m− 1 2 (for bottlenecks' positions fixed.) This solution is formed by waves, whose velocity sign is constant. Setting y l = y(t l ) and using the Gauss-Green formula, the scheme is expressed in the integral formulation as Step 2. We compute y l+1 P using the density ρ l+1 obtained at Step 1, as for the case of a single bottleneck. More precisely, there exists m such that y l P ∈ [x m−1− 1 2 , x m− 1 2 ); y P moves with velocity ω P (ρ l+1 m−1 ) until the interaction time ∆t in with the line x = x m− 1 2 , provided that ∆t in ≤ ∆t. After that time, y P moves with speed ω P (ρ l+1 m ). Finally, if ∆t ≤ ∆t in then we set y l+1 P = y l P + ∆t ω P (ρ l+1 m−1 ), Step 3. We compute y l+1 i by backward recursion on i = P−1, . . . , 1. More precisely, we define a trajectory y i (t) for every t ∈ [t l , t l+1 ] and set y l+1 . Now, fix i and assume to have defined y j (t) for j ≥ i + 1 and t ∈ [t l , t l+1 ]. Let m be such that y l i ∈ [x m−1− 1 2 , x m− 1 2 ) and define ∆t in the time at which the line y l i + (t − t l )ω i (ρ l+1 m−1 ) intersects the vertical line x = x m− 1 2 . If ∆t ≤ ∆t in we set otherwise we simply set y i (t) =ỹ i (t) for every t ∈ [t l , t l+1 ].

Numerical tests
In this section we show some simulations obtained by models (1.1)-(1.2) and (6.3). In particular, we focus on interactions of bottlenecks with different type of waves, namely shocks and rarefactions. For both models we implement the corresponding GOF scheme.

Model (1.1)-(1.2)
To compute the approximate solution to (1.1)-(1.2) we use the numerical algorithm presented in Section 5.2. We set the space and the time mesh size of the grid to be ∆x = 0.02 and ∆t = 0.01. Setting parameters, we assume w(ρ) = w max (1 − ρ) with w max = 0.4 and wherev = 1, v = 0.6 and β = 0.1. In Fig. 8 a moving bottleneck starting from y 0 = 0.5 interacts with a rarefaction wave centred in x = 1.4, generated by the Riemann datum ρ ℓ = 0.9 for x ≤ 1.4 and ρ r = 0.45 for x > 1.4. The vehicle is initially moving with low speed (w = 0.04) due to the high density until it interacts with the first characteristic of the rarefaction; after that it accelerates since the density is decreasing.
In Fig. 9 a bottleneck starting from y 0 = 0.5 is interacting with a backward moving shock starting from x = 1.4. The latter is generated by the Riemann datum ρ ℓ = 0. x > 1.4, so that the speed of the shock, before it interacts with the bottleneck, is λ = −0.2. On a following time interval, during the interaction with the bottleneck (from around t = 1 until around t = 2.5 in Fig. 9), the speed of the shock almost vanishes, because of the capacity drop caused by the bottleneck.   10 shows a shock and a rarefaction which interact with each other and with the bottleneck. The rarefaction wave is centred in x = 0.6 while the shock is starting from x = 1.6.

Model (6.3)
Simulations for the system in (6.3) are provided by means of the numerical algorithm defined in Section 6.3. Interactions of three bottlenecks and different type of waves are considered. We set the space and time mesh size to be ∆x = 0.02 and ∆t = 0.01. Velocity functions for slow moving vehicles are ω i (ρ) = ω i max (1 − ρ), with ω 1 max = 0.49 and ω i max = 0.4 i = 2, 3. We set again where, as before,v = 1. In Fig. 11, we set v i = 0.5 and β i = 0.25 for i = 1, 2, 3. The initial datum for the density is piecewise constant with single discontinuity at x = 2.5 and values ρ ℓ = 0.9 and ρ r = 0.1. This gives rise to a rarefaction wave. Each slow vehicle in the road is interacting with the rarefaction wave, but in different ways. The second vehicle, starting from y 2 (0) = 1.5, is accelerating, but not as much as the first one, starting from y 3 (0) = 2. This is due to differences in the rate of density decrease. The last one, starting from y 1 (0) = 1, has maximal velocity w 1 max higher than the others. However, it is not capable of strong acceleration, since its initial distance w.r.t. the second is equal to β 1 + β 2 = 0.5, thus its speed is bounded by that of the second vehicle. In Fig. 12, we simulate an homogeneous case, i.e. all bottlenecks have the same characteristics. We set v i = 0.5 and β i = 0.25 for i = 1, 2, 3. The initial datum is piecewise constant with ρ ℓ = 0.85 for any x ≤ 3.5 and ρ r = 0.95 for any x > 3.5, thus generating a shock with negative speed λ = −0.8. Initial  data for bottlenecks are y 1,0 = 1, y 2,0 = 2 and y 3,0 = 3. Starting from t = 0, due to capacity dropping, the density behind each slow vehicle increases giving rise to a queue of cars. This phenomenon generates further waves travelling and interacting with bottlenecks as well as the main shock. When a bottleneck cross the shock, it enters a region with higher density, so that its velocity decreases. On the other hand, during the interaction with each slow vehicle, the speed of the shock decreases substantially, as