Unformatted text preview: Notes for Mathematical Modeling, Spring 2011
Esteban G. Tabak 1 Traffic Flow In our first modeling exercise, we will build and study a mathematical model for traffic flow. There are various possible approaches to modeling traffic. First, we can distinguish between discrete and continuous models. In a discrete model, one follows individual cars, studying their position and speed as functions of time, and modeling their evolution in terms of interactions with the other cars, traffic lights, etc. In a continuous model, one looks at the road from far enough that the behavior of individual cars becomes unimportant, and deals instead with average concepts, such as the number of cars per mile and the cars mean speed. This is analogous to the distinction, in the study of fluid motion, between the molecular and the continuous viewpoints: the former is more accurate for the description of behavior in the nano-scales, while the latter is far more practical for a macroscopic description. A second distinction is between stochastic and deterministic models. In a stochastic model, the evolution of the traffic is not fully specified by the model, but depends on random variables with specified probability distributions, such as the probability of a car exceeding a certain speed, of a number of cars to cluster in a section of the road, etc. In a deterministic model, the evolution is dictated exclusively by the model's laws. In our first few classes, we will describe traffic flow through continuous, deterministic models. In the process, we will learn something about the modeling process and on a class of partial differential equations (hyperbolic conservation laws), and unveil some simple reasons for a number of mysteries that we observe daily on the road, such as the sudden formation of seemingly causeless traffic congestions, and their eventual dissipation into thin air. 1.1 Variables In describing the traffic through a road from a macroscopic perspective, a sensible choice for independent variables is that of the distance along the road, that we will denote x, and the time t. The main macroscopic dependent variables, on the other hand, are and Q, the car density and flux: number of cars per unit length of the road, and number of cars crossing a section per unit time. A related variable is U , the mean velocity of cars. Thus we seek three functions 1 of two variables: (x, t), Q(x, t) and U (x, t). The three dependent variables are related by the identity Q = U , (1) which follows from the following considerations: The number of cars crossing a section at speed u in the time interval dt equals the number of cars with that speed in the spatial interval dx = u dt. Hence, for a uniform car speed, (1) follows. To take the average over the various car speeds in the general case, we need to introduce the density in phase space r(x, t, u) , so that u2 r(x, t, u) du
u1 represents the amount of cars with speeds between u1 and u2 per unit length of the road at time t and in position x. Clearly (x, t) = r(x, t, u) du
0 and Similarly, we may introduce a flux q(x, t, u) in phase space. Then Q= q(x, t, u) du = r(x, t, u) u du = U ,
0 0 r(x, t, u) u du 0 U = . r(x, t, u) du 0 justifying (1). 1.2 Equations Our "physics" will be limited to a simple, rather intuitive fact: in a road segment with no entrances or exits, the number of cars does not change. The differential form of the equation for car conservation is t + Q x = 0 . It follows from the integral identity x2 ((x, t2 ) - (x, t1 )) dx =
x1 (2) t2 t1 (Q(x1 , t) - Q(x2 , t)) dt , (3) which equates the change in the number of cars in a given segment of the road to the difference between the cars entering and those leaving the segment through 2 its two ends. The way the derivation goes is to rewrite (3), by invoking the fundamental theorem of calculus, in the form x 2 t2 (t + Qx ) dt dx = 0 ,
x1 t1 which must hold for any values of (x1,2 , t1,2 ). Hence, if the integrand is continuous, (2) follows. Equation (2) is a single constraint for two unknowns; it alone cannot determine the evolution of the traffic. In order to close the system, we need another constraint. The simplest one follows from the realization that the typical car speed depends on the level of traffic congestion, i.e. the car density. Turning this qualitative statement into a sharp functional relation, we may propose that Q(x, t) = Q((x, t)) , (4) where Q() is a prescribed function, that depends on the characteristics of the road. Under this assumption, equation (2) becomes t + Q () x = 0 , (5) a closed equation for the single variable . Since the mean speed U should clearly be a decreasing function of , reaching zero at the fully congested level max , then the flux Q(), computed from (1), must be concave, with zeros at both = 0 and = max , and achieving its maximum value in between, say at = . Figure 1 plots one such function, corresponding to the simplest choice Q = (1 - ) . (6) 1.3 Characteristics We will solve the partial differential equation (5) for (x, t) using the fact that there are certain lines in the (x, t) plane along which it becomes an ordinary differential equation. To see this, consider the time derivative of along a prescribed curve x = X(t): d dX (X(t), t) = t + x . dt dt If we pick the curve X(t) in such a way that dX = Q ((X(t), t) , dt then (5) yields d (X(t), t) = 0 . dt 3 (8) (7) Constituitive equation Q() 0.25 0.2 Q = (1 - ) 0.15 0.1 0.05 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 1: Simple constitutive relation for traffic flow Equations (7, 8) are the characteristic equations for (5), and the family of curves X(t, ) satisfying (7) , its characteristics (Here is any parameter, such as the initial value X(t = 0), distinguishing the family members.) It follows that the density is constant along characteristics, and the characteristic curves X(t) are straight lines. The characteristics are often described as directions along which information propagates. In our case, this information is just the value of the solution (x, t). It is worth noting that, for traffic flow, the characteristics are always slower than the cars: dX dQ d( U ) dU = = =U + U, dt d d d since U is a decreasing function of . This fact makes driving doable, since information about the state of the traffic always comes from the cars ahead, not from those behind! A consequence of the characteristics being straight is that the solution at any point (x, t) equals that at (x - Q ((x, t) t, 0). This allows us to write the following implicit solution to (5) subject to initial data (x, 0) = 0 (x) : (x, t) = 0 (x - Q ((x, t)) t) (9) 1.4 Shocks From the construction of characteristics in the section above, summarized in the characteristic equations (7, 8) or alternatively in the implicit solution (9), we see that the information yielding a solution (x, t) to (5) "travels" along 4 straight lines in the (x, t) plane. These straight lines, however, are not generally parallel to each other, since Q () adopts different values for different values of . Hence characteristics will eventually intersect, yielding contradictory information about the solution (x, t). What becomes of the solution after the first such intersection occurs? The concavity of Q() for traffic flow implies that Q is in fact a decreasing function of . Hence lower traffic densities will travel faster than higher ones. Consequently, any nonconstant initial profile will deform as time progresses, with the areas of low density catching up with those of higher density ahead of them, and getting farther away from those behind them. This phenomenon is illustrated in figure 2, which displays various snapshots of a numerical simulation of (5), with the constitutive relation given by (6) (The algorithm used will be described a few subsections below), and initial data 0 (x) = 1 (1.5 + sin(x - )). 4 Notice that at approximately t = 2, the solution seems to deform to the point of developing a vertical slope at x 3.6. This corresponds to the first crossing of characteristics. This breaking time and position could have been predicted from the implicit solution (9): differentiating this, one obtains x (x, t) = 1+ (x0 ) 0 (x ) Q ( (x )) t 0 0 0 0 , with x0 = x - Q ((x, t)) t. Hence the slope x (x, t) blows up at the minimum positive value of -1 t= . 0 (x0 ) Q (0 (x0 )) Since Q is negative, blow up takes place at a location of positive slope, agreeing with the picture of a region of low density catching up with a more congested one. In our case, Q = -2 for all values of , so blow up occurs at the point of maximal slope = 1/4, achieved at x0 = . The corresponding time is t = 2, 0 at position x = +Q ((, 0)2 = 3.64, in agreement with the numerical results of figure (2). After the blowup time, the implicit solution given by (9) is multivalued; typically three characteristics converge at each point of the overlap region, each carrying a different value for . The profile corresponding to the initial data above, evaluated at t = 4, is presented in figure (3). Even though we shall see examples later where multivalued solutions of this kind make perfect physical sense, this is certainly not the case of traffic flow, where only one car density can exist at each point of the road. A possible solution is to discard two of the three solutions in the region of characteristic overlap, keeping only one. At the point where the solutions switches from a characteristic branch to another, the solution will be discontinuous. An example of this construction is shown in figure (3) with a dotted line; the locus of discontinuity is called a shock wave, a name originating in gas dynamics. In traffic flow, it corresponds to a 5 0.8 0.7 0.6 0.5 0.4 0.3 t=0 t=1 0.2 t=2 0.1 0 1 2 3 x 4 5 6 7 Figure 2: Snap shots of a numerical solution phenomenon that we frequently observe: a car driving relatively fast in light traffic is suddenly forced to slow down and enter a more congested segment of the road, with no apparent reason explaining the sharp transition between the two. Yet it remains to determine where to position the shock, and more generally which meaning to give to a solution to (2) once it becomes discontinuous. The answer lies in the integral formulation (3). Notice that this simpleminded statement of car conservation does not depend on either or Q being continuous. Hence we may think of (3) as more fundamental than (2), and call weak solutions of (2) those functions, not necessarily continuous, that satisfy (3) for all choices of x1,2 , t1,2 . Notice first that (3) can be generalized, switching from rectangles to arbitrary regions in the (x, t) plane: ( dx - Q dt) = 0 , (10) where is the closed boundary of the domain . Equation (10) can be derived in many ways: physically, by noticing that it represents the generalized flux of cars across , and so just states that cars are neither created nor destroyed; geometrically, by approximating the surface by a union of rectangles and applying (3); and analytically, by invoking Green's theorem in conjunction with (2). To see that (10) really stands for car flux in (x, t) space, think of an observer moving at speed c. The number of cars that the moving observer would encounter during a small time interval t equals (c - U ) t, i.e.: the car density times the distance traveled with respect to the moving traffic. But this equals c t - U t = x - Q t, justifying (10). 6 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1 2 3 4 x 5 6 7 8 Figure 3: A multivalued solution We shall see next that the integral formulation (10) allows us to compute the speed at which a shock wave propagates. To this end, consider the location of a shock in the (x, t) plane, and draw the boundary of a thin domain , including a small enough segment of the shock that its speed c = dx/dt, as well as the values - and + to the left and right of it, can be approximated by constants (see figure(4)). When applied to this domain, (10) yields c= dx Q(+ ) - Q(- ) = . dt + - - (11) This relation between the speed of the shock and the jumps of and Q across it. is called the RankineHugoniot condition, or simply the jump condition across shocks. Is the shock speed c all we need to compute weak solutions to (2) near shocks? The answer is positive if the characteristics on both sides of the shock are converging toward the shock, or at least not emanating from it. For then the information about - and + is updated by the characteristic flow on each smooth side of the shock, while the speed of the shock itself is given by (11). Notice that this convergence of characteristics is what gave rise to the shock to start with, with areas of lower car density catching up with more congested ones. Hence we should expect the condition Q (- ) Q(+ ) - Q(- ) Q (+ ) + - - (12) to apply to realizable shocks. This is the Lax condition for shocks, that will reappear below relating to issues of uniqueness and stability. Notice that the 7 3 2.5 2 shock 1.5 x t t 1 0.5 0 0 0.2 0.4 0.6 0.8 x 1 1.2 1.4 1.6 1.8 Figure 4: Control volume around a shock leading to the RH conditions downward concavity of Q() for traffic flow makes this condition equivalent to - + , i.e., realizable shocks always have the more congested area ahead. Notice that the jump condition (11) does not depend just on the differential equation (5), but on the integral principle (10) from which it derives. To see that this is not determined uniquely by the differential equation alone, take any smooth function F (), and rewrite (5) in the form Ft + G x = 0 , where G = F () Q () d. The shock speed for (13) is given by c= G+ - G- Q+ - Q- = + . + - F- F - - (13) Hence different forms of writing the same differential equation give rise to different weak solutions. It is crucial, therefore, to know which is the conservation principle from which the equation has been derived. 1.5 The Riemann Problem Equation (2) is an example of a conservation law, referring specifically in this case to the conservation of cars. Other phenomena modeled through systems of conservation laws include gas dynamics and river flow. For such systems, one of the basic building blocks for both theory and numerics is the Riemann problem: a set of initial conditions consisting of two constant states separated by a discontinuity. For traffic flow, this problem is given by equations (2) and 8 (4), subject to the initial condition - (x, 0) = + for x < 0 for x > 0 . (14) Notice that both the PDE (2) and the initial condition (14) are invariant under the stretching x x , y y . This implies that if (x, t) is a solution to the Riemann problem, so is ( x, t), for any value of . Hence for any two points (x1 , t1 ), (x2 , t2 ), such that x2 /x1 = t2 /t1 , the solution will adopt the same value. In other words, (x, t) = R(x/t) . If R( = x/t) is a smooth function, plugging it into (5) yields (Q (R) - ) R () = 0 , with alternative solutions R = constant and Q (R) = . (15) (16) Since Q is a decreasing function of R, the solution given by (16) is a decreasing function of = x/t. Hence, if - > + , we can build a solution to the Riemann problem (14) by combining (15) and (16) in the following way: - for x/t Q (- ) -1 (x, t) = (17) Q (x/t) for Q (- ) < x/t < Q (+ ) + for x/t Q (+ ) . An example of such a solution, called a rarefaction, is plotted in figure 5a. Since in our example Q (R) = 1 - 2R is a linear function, the resulting rarefaction is linear in = x/t too. If - < + , on the other hand, no such matching of solutions is possible. However, in this case, we may fit a shock wave between - and + , with speed given by (11): + )-Q(- - for x/t Q(+ -- ) (x, t) = (18) + )-Q(- + for x/t > Q(+ -- ) . This discontinuous solution, plotted in figure 5b, completes the solution of the Riemann problem for traffic flow. 9 a: A rarefaction 0.8 0.8 b: A shock 0.7 0.7 0.6 0.6 0.5 0.5 0.4 -1 0 =x/t 1 2 0.4 0.3 0.3 0.2 -2 0.2 -2 -1 0 =x/t 1 2 Figure 5: Solution to the Riemann problem for (2), with constitutive relation Q = (1 - ). a: A rarefaction. b: A shock wave. 1.6 Stability of shocks and uniqueness of weak solutions The attentive reader may be slightly troubled by the sudden end of the previous section. There, we found two qualitatively different solutions to the Riemann problem, depending on the relative sizes of - and + : a rarefaction and a shock wave. Clearly, the rarefaction only fits the case when - > + ; there is no way to match it to the outer, constant solution otherwise. But couldn't the shock solution (18) apply to this case as well? Indeed, there is nothing in the discontinuous solution (18) that requires that - < + . Hence we have a problem of nonuniqueness: both (17) and (18) are solutions to the Riemann problem (14), when - > + . In a way, by extending our definition of a solution to (2) to include weak solutions, i.e. piece-wise smooth functions satisfying the integral constraint (10), we gave ourselves too much freedom: now we have more solutions than we need; certainly at least two for the Riemann problem with - > + . Is there a way (and a need) to eliminate all but one of them? We shall argue below that indeed there is a need and a way, and that the only valid solution to (14) with - > + is in fact the rarefaction (17). Let us start with some plausibility arguments. We must recall the reason why we were forced to generalize the notion of solution to (2): because different characteristics may catchup with each other and cross, thus yielding multi valued solutions. But, for this to happen, the characteristic starting behind needs to be faster than the one ahead; for traffic flow, this means that it needs to have a lower density . Hence we should not expect a shock to form with 10 - > + : it would not arise naturally from smooth initial data. This is a first argument in favor of rejecting shocks with the "wrong" ordering of - and + . Yet, in posing the Riemann problem (14), we just started with a discontinuity in the initial data, which could well have had - > + . Is there any reason why such discontinuity, once there, could not survive as a shock? To see that there is one reason, recall our logic behind the introduction of the Lax condition (12): that the information required to evolve the smooth solutions on each side of the shock should be available; i.e., that the characteristics carrying this information should be converging toward the shock, not diverging from it. When this condition is violated, as is the case of (18) with - > + , it is not clear where the information comes from, that keeps the solution at the two sides of the shock at their constant values. This is a second argument in favor of rejecting the "rarefaction shock". Here comes a third, related argument, based on stability. When one writes a differential equation as a model for a physical process, it is often a desired property that it should lead to well posed problems: small changes in the data, or small thermal perturbations along the way, should not give rise to dramatically different solutions. Yet this is exactly the case of the shock solution (18) when - < + . For consider, instead of (14), the slightly perturbed initial condition - for x < - x for - < x < + , (x, 0) = (19) + for x > + . where is an arbitrarily small number. This smooth, monotonically increasing initial condition, will evolve along diverging characteristics, hence never breaking, and approaching a profile arbitrarily close to (17) as 0. This is, therefore, the solution of choice, if one is to require wellposedness of the initial value problem. An alternative wording of the argument above will state that a shock wave not satisfying the Lax condition (12) will be unstable to small perturbations: if the shock should be slightly smoothed at any time, it will immediately expand into a fully blown rarefaction. Shocks satisfying (12), on the other hand, will recompose shortly after being pushed apart, due to the "compressive" effect of the converging characteristics. A condition such as (12), that restores unicity to the weak formulation of a system of conservation laws, is often denoted an entropy condition. The origin of this name is to be found in gas dynamics, where physical shocks create entropy, while the unstable, unphysical ones manage to dissipate it, contrary to the second principle of thermodynamics. 1.7 Numerics: Godunov's method We now switch from theoretical considerations to the problem of simulating the traffic flow equation (2) numerically. Three main challenges associates with solving systems of conservation laws numerically are the following: 11 The solutions one is after are typically not smooth, due to the presence of shock waves. Hence methods that approximate derivatives by assuming some local smoothness, such as finite differences or spectral collocation, are not wellsuited to the task. The speed of propagation of shocks depends not just on the differential equations being solved, but on the integral form from which they derive (and these, as shown above, are not uniquely determined by the differential equations.) Hence the numerics needs to know about these integral principles. Information propagates along particular directions, the system's characteristics. If the numerics does not mimic this information flow correctly, instabilities ensue. One numerical method that addresses successfully all the issues above, is due to Godunov, who devised it originally in the context of gas dynamics. This is the method that has been used to produce the numerical results in figure (2), and is the one that we shall choose to describe and use in these notes. A first ingredient of the method is that it is in conservation form, hence enforcing the integral formulation of the problem: the number of cars is preserved exactly by the numerics. The way this is achieved is through the introduction of a discrete spatial grid with mesh size x (that we shall for simplicity assume uniform), and a discretized dependent variable j approximating not quite (xj ), as in most numerical methods, but rather its average over a grid cell: xj +x/2 1 j (x) dx . x xj -x/2 Every time step t, j is updated by an equation of the form n+1 = n + j j t Qj-1/2 - Qj+1/2 , x (20) where Qj+1/2 is an approximation to the flux Q(x + x/2) averaged over the time interval (t, t + t). Notice that, independently of the accuracy of this approximation to Q, the form of (20), leading to telescopic sums, makes the total mass x j
j be exactly conserved for all times (all changes are due to fluxes at the boundaries of the domain.) Hence the only task left to complete the algorithm is to decide how to compute the approximate fluxes Qj+1/2 . What Godunov proposed, at this stage, is to compute essentially the exact fluxes corresponding to an approximate profile (x, tn ) consistent with the local numerical averages n . He chose the simplest j such profile: a piecewise constant function, defined by (x, tn ) = n j for xj - x/2 < x < xj + x/2 12 0.6 0.5 , j 0.4 0.3 0.2 0.1 0 1 2 3 x 4 5 6 Figure 6: Approximation of by a piecewise constant function (see figure 6.) In order to compute the fluxes Q at the cell interfaces, it is enough to notice that, locally, these are the solutions to Riemann problems: two initially constant states separated by a discontinuity. Since information travels along characteristics at a finite speed, nonneighboring cells will not interact with each other, provided that we pick t small enough that characteristics coming from two adjacent interfaces do not intersect. As we saw in the prior subsection, the solution to the Riemann problem will be constant along lines of constant x/t; for our procedure, we are interested in Q((x/t = 0)), the flux at the interface between cells. This completes the description of the (first order) Godunov method, applied to single conservation laws. Its Matlab implementation traffic.m can be found in the class' website. 1.8 An empirical validation of the theory The simple kinematic model for traffic flow described above illustrates the power of the theory of partial differential equations at its best: using just a handful of intuitive facts about traffic cars don't suddenly appear or vanish; the lighter the traffic, the faster the cars one develops a model with rich mathematical structure characteristics, weak solutions, shock waves, entropy conditions capable of explaining a mysterious fact of everyday driving: the sudden appearance of traffic congestions, and their smooth dissolution into nothingness. Here's another empirical fact that I have noticed while driving in the streets of New York and Buenos Aires, and that the theory explains rather well: You are in a long traffic queue in a two-way street say Park Avenue, 13 waiting for the green light. When the light finally turns green, the first car in the queue moves, then the second, and so on: a wave propagates backward through the queue, telling the cars when to start. On the other hand, you see on your left the first cars on the other lane, going the opposite way, as they start, accelerate, and soon reach your point in the queue. My observation is please confirm this with your own experience! that the time at which the first car of the opposite lane reaches your position in the queue agrees almost exactly with your starting time. Is there an explanation for this peculiar coincidence? In the model that we have developed above, the backward starting wave in the queue is the leading edge of a rarefaction wave, with speed given by c = Q (max ), where max is the car density at maximal congestion. The first cars on the opposite lane, on the other hand, move at speed -U (0) the minus sign because of our other way perspective, the zero density because they move into an area with no cars. Then our observation is equivalent to the statement Q (max ) = -U (0), linking the flow at the two extreme opposite conditions: maximal congestion and no traffic at all. Is there a reason why these two should agree? Consider the simplest model for the speed of cars: maximal speed when the road is empty, decreasing linearly to zero at maximal congestion: U () = U (0) max - . max The corresponding flux function is given by Q() = U () = U (0) with max - , max Q (max ) = -U (0) as the observation suggests. Then the observation is consistent with the simplest rule for driving speed as a function of traffic conditions. That it should apply to cities like New York and Buenos Aires, where people drive following complex, rather idiosyncratic rules, remains nonetheless a mystery! 14 ...
View Full Document
This note was uploaded on 03/25/2011 for the course MATH 250 taught by Professor Tabak during the Spring '11 term at NYU.
- Spring '11