Web page maintained by Evans M. Harrell, II, harrell@math.gatech.edu.
This section indicates a common method to model the flow of traffic. The importance of the problem for us is that it is a nonlinear system for which we have some intuition and, more important, it introduces the notion of a "shock" curve.
We consider the flow of cars on a highway under the assumptions that cars do not enter or leave the highway. We take the x-axis along the highway and assume that the traffic flows in the positive direction. Suppose that u(t,x) is the density representing the number of cars per unit length at the point x on the highway at time t. Thus,
I(a,b, )u(t,x) dx
is the number of cars between a and b at time t. Let f(t,x) be the flow of cars per unit time. Then
I(r,s, )f(t,x) dt
is the number of cars going past x from time r to time s.
We assume a conservation law which states that the change in the total amount of a physical quantity contained in any region of space must be equal to the flux of that quantity across the boundary of that region. The total change is what goes in minus what goes out. (For this example, the implication is no wrecks! That is, nothing enters a region and stays there.) In this case and recalling that u is the "density" of cars, the time rate of change of the total number of cars in any segment a <= x <= b of the highway is given by
Total change = F(d ,dt)I(a,b, )u(t,x) dx = I(a,b, )F([[partialdiff]]u,dt) dx.
This rate of change must be equal to the net flux across a and b given by
f(t,a) - f(t,b)
which measures the flow of cars entering the segment at "a" minus the flow of cars leaving the segment at "b." Thus we have the conservation equation
I(a,b, )F([[partialdiff]]u,dt) dx = f(t,a) - f(t,b).
But, by the fundamental theorem of relating integral and differential calculus, the right side can be re-written as
f(t,a) - f(t,b) = - I(a,b, )F([[partialdiff]] ,[[partialdiff]]x)f(t,x) dx
so that
I(a,b, ) F([[partialdiff]] ,[[partialdiff]]t)u(t,x)dx = - I(a,b, )F([[partialdiff]] ,[[partialdiff]]x) f(t,x) dx
Since the integrands are continuous and this equation holds for all [a,b], it follows that the integrands must be equal so that we have the partial differential equation
F([[partialdiff]]u,[[partialdiff]]t) + F([[partialdiff]]f,[[partialdiff]]x) = 0. (6.1)
If we add the hypothesis that the flow depends on the density alone: f(t,x) = f(u(t,x)) and that g = f[[minute]] then
F([[partialdiff]]u,[[partialdiff]]t) + g(u)F([[partialdiff]]u,[[partialdiff]]x) = 0. (6.2)
This is a non-linear equation such as we considered in the last section. We add an initial value which should be a model of an initial distribution of traffic. Call this initial distribution j(x), so that
u(0,x) = j(x) .
In each of the following examples, remember that the characteristics for equation (6.2) are given by x(t) = g(j([[xi]])) t + [[xi]]. On these curves, u is given through the implicit relation
x = g(u(t,x)) t + u(0,x). (6.3)
In the discussions below, we take j to be given by
j(x) = BLC{(A(1 if x < 0, 1-x if 0 < x < 1, 0 if 1 < x )) (6.4)
for each of the examples. This function could be considered as a model for traffic stopped at t <= 0 for a traffic light which turns green at t = 0. For negative x, the density is one. For x > 1 and t <= 0, the density is zero. The transition region provides a smooth connection. The function j has graph as suggested here:
Figure 6.1
Example 6.1
For this first example, we take f(u) = u/2. Because f[[minute]](u) = 1/2, the partial differential equation (6.2) becomes
F([[partialdiff]]u,[[partialdiff]]t) + 1/2 F([[partialdiff]]u,[[partialdiff]]x) = 0.
Characteristics are x = t/2 + [[eta]] for all real numbers [[eta]] and the characteristics fill a half plane with straight, parallel lines. In this case, the solution for the equation is u(t,x) = j(x - t/2) whose graph follows. It will be of value to give this picture a physical -- a traffic -- intrepretation. A way to visualize this traffic interpretation is to choose x, visualize standing at x and let time move forward. Or, What does u(.8,x) represent and what is its graph?
Figure 6.2
Example 6.2
In this second example, we take f(u) = u (1-u). Then f[[minute]](u) = g(u) = 1-2 u. With [[phi]] as specified above, characteristics are given by
If 1 < [[xi]], then j([[xi]]) = 0, g(j([[xi]]))= g(0) = 1, and x = t + [[xi]].
If [[xi]] < 0, then j([[xi]]) = 1, g(j([[xi]]))= g(1) = -1, and x = -t + [[xi]].
If 0 < [[xi]] < 1, then j([[xi]]) = 1-[[xi]], g(j([[xi]]))= g(1-[[xi]]) = -1 + 2[[xi]], and x = (-1+2[[xi]])t + [[xi]].
Here are graphs of the characteristics. Note they do not cross!
Figure 6.3
We should be able to give an explicit formula for u. If t and x lie below the line t + x = 0, then j([[eta]]) = 1. If t and x lie above the line t + x = 1, then j([[eta]]) = 0. Between these two x = [[eta]] (2t + 1) + t. Or,
[[eta]] = F(x+t,2t + 1),
and, j([[eta]]) = 1 - [[eta]] = 1 - F(x+t,2t + 1) = F(t + 1 - x,1+2t).
U(t, x)=BLC{(A(1 if x < -t, ,(t+1-x)/(1+2t) if -t < x < t +1, ,0 if t +1 < x))
We should be able to draw the graph.
A "traffic" interpretation of the graph is suggested in the next example, which pulls from this example what is important.
Figure 6.4
Example 6.3
We make the assumption that the flow rate f depends on t and x only through the density. That is, f = f(u) for some function f. This assumption seems to be reasonable in the sense that the density of cars surrounding a given car effects the speed of that car. The functional relation between f and u depends on many factors, including speed limits, road characteristics, and driving experience of the drivers. (In this regard, a knowledge by the driver of very elementary physics can be important.)
We consider here the particular relation f = u. v(u) where v is the average local velocity of cars. In view of this relation,
F([[partialdiff]]f,[[partialdiff]]x) = v(u) F([[partialdiff]]u,[[partialdiff]]x) + u v[[minute]](u) F([[partialdiff]]u,[[partialdiff]]x)
the differential equation (6.1) reduces to the nonlinear equation
F([[partialdiff]]u,[[partialdiff]]t) + g(u)F([[partialdiff]]u,[[partialdiff]]x) = 0.
where g(u) = f[[minute]](u) = v(u) + u v[[minute]](u).
This graph in Figure 6.5 shows how v might depend on u. For this picture, v(u) = 1-u2.
Figure 6.5 Graph of v(u)
The local velocity v(u) is a decreasing function of u: the more cars around, the slower the traffic ... except in Atlanta. Also, v(u) has a finite maximum value vmax at u = 0 and decreases to zero at u = umax = um. For the value u = um, the cars are "bumper to bumper." Since f = u v(u), then f(u) = 0 when u = 0 and when u = um. A graph of f might be graph in Figure 6.6. For this illustration, f(u) = u.v(u) = u.(1-u2).
Figure 6.6 Graph of f(u)
For our example, f is an increasing function of u until it reaches a maximum value fmax for uM = .5 and then decreases to zero at um = 1.
From (6.3) we have a wave propagation velocity g(u) = v(u) + u v[[minute]](u). Since v[[minute]](u) < 0, g(u) < v(u) and the wave velocity is less than the car velocity. That is, waves propagate backwards through the stream of cars and drivers start to stop before a disruption ahead. The function f in an increasing function in [0,uM], and f is a decreasing function in [uM,um] and attains a maximum at uM. Hence, g(u) = f[[minute]](u) is positive in [0,uM], zero at uM, and negative in [uM,um]. The result is that in [0,uM] the waves propagate forward relative to the highway, are stationary at uM, and travel backwards in [uM,um].
The solution with initial value [[phi]] is
u(t,x) = [[phi]]([[xi]]), x = [[xi]] + t f([[xi]])
where
f([[xi]]) = g([[phi]]([[xi]])).
Since g[[minute]] = f[[minute]][[minute]](u) < 0, the g(u) is a decreasing function of u. This means that breaking occurs at the left due to a formation of shock at the back. Waves propagate slower than the cars, so drivers enter such a local density increase from behind, and they must decelerate rapidly through the shock. They speed up slowly as they get out from the crowded area.
Please be aware that when cars are loaded onto the train for shipping under the English Channel, there is no worry that a "shock" will ripple through the tunnel due to the lack of a smooth solution for a pde! In case the cars are loaded on a train, while f(t,x) may be u(t,x).v(u(t,x)), v is a simple constant function. Alternatively, on I-85, why not have everyone put cruise control on 68? Then any distribution can move down the highway. The trouble comes when one of the cruise controls is not accurate. Then the density changes; and the following example becomes appropriate.
Example 6.4
Let u(t,x) be the density of cars on the I-85 and d be the of a density of a bumper-to-bumper parking lot. Suppose that
v(t,x) = 68 B(1 - F(u(t,x),d)).
Think of the straight line graph of v as a function of u. Note the values of v when u = 0 or d. This example seems to fit Atlanta for when u = d/2, so that there is one car length between vehicles, the velocity is 34 miles per hour. Too fast!!
As in Example 6.3, we take the flow rate to be u(t,x).v(t,x), so that
ut + (u v)x = 0.
Now,
(u v)x = [u 68 B(1 - F(u(t,x),d)) ]x = 68 B(1 - F(2 u,d)) ux.
Thus, we have the equation
ut + 68 B(1 - F(2 u,d)) ux = 0.
The associated 3 characteristic equations have solution
T(t) = t, X(t) = 68 B(1 - F(2 f(x),d)) t + [[eta]], z(t) = f([[eta]])
where the initial density is u(0,x) = f(x).
With some choices of f(x), the characteristic lines are parallel: for example if f(x) is constant. With other choices for f(x), the lines do not cross for t > 0; in particular, take f(x) to be decreasing. Finally, here is a distribution where they do cross. This latter is the case if f is increasing on some interval. At the point of crossing, this method of finding solutions needs further investigating.
We illustrate these three choices with different f's.
Illustration 1: Take d = 20 (cars/100 yards), f1(x) = 2. Then
X(t) = 68.4.t /5 + x.
Characteristic lines are parallel, and cars move down the road uniformly.
Illustration 2: With d = 20, take
f2(x) = F(3,200) F(550 - 3 x, x + 50) .
This function predicts 5 cars in the first 50 yards and about 2 in the second 50 yards. The following syntax verifies these calculations.
* f2:=x->3*(550-3*x)/(200*(x+50)):
eq1:=evalf(int(f2(x),x=0..50));
eq2:=evalf(int(f2(x),x=50..100));
The function X would be given by
X(t) = 68 (1 - F(f2(x),10) ) t + x.
To simplify the coefficient of t,
* mx:=simplify(68*(1-f2(x)/10));
and get
X(t) = F(119,500) F( 287 x + 14050 , x+50 ) t + x.
To see that the slopes of these lines are increasing as x increases, note the derivative of the slope as a function of x:
* diff((287*x+14050)/(x+50),x);
* simplify(");
F(d F( 287 x + 14050 , x+50 ) ,dx) = F(300,(x+50)2).
This is positive for all x; so that the slope must be increasing. Hence the characteristic lines do not cross.
To think of how the solution surface looks, we could solve this equation for [[eta]]:
X(t,x) = [[eta]].
This would give x as a function of t, and [[eta]]. We would plot
u(t,[[eta]]) = f(x(t,[[eta]])).
Solving the equation X(t,x) = [[eta]] for [[eta]] is really a mess.
Here's an alternate idea. Since over each characteristic line, the function is constant, we plot the surface parametrically:
u = { {t, X(t,x), f(x)}: 0 <= t, 0 <= x} }.
Here is the MAPLE realization of that idea.
* sol:=solve(X(t,eta)=xx,eta);
* simplify(sol[1]);
* simplify(sol[2]);
* plot3d([t,X(t,x),f2(x)],t=0..0.1,x=0..2,
orientation=[-145,45],axes=NORMAL);
A portion of the solution surface for Illustration 2.
Illustration 3: With d = 20, take
f3(x) = - F(47,100) + F(54,5) F(1,1+(x-2)2)
This initial distribution is designed so that there are about 3 cars in the first mile and 8 in the second.
* a:=-47/100: b:=108/10:
* f3:=x-> a + b/(1+(2-x)^2);
* eq1:=int(f3(x),x=0..1);
eq2:=int(f3(x),x=1..2);
* evalf(eq1); evalf(eq2);
One has only to look at the a graph of the initial distribution f on the interval [0,2] to see that it is increasing.
* plot(f3(x),x=0..2);
Graph of u(0,x) for 0 <= x <= 2.
The characteristic lines should cross.
* X:=(t,x)->68*(1 - f3(x)/10)*t + x;
* plot({[t,X(t,0),t=0..0.1],[t,X(t,.5),t=0..0.1],[t,X(t,1),t=0..0.1]});
Graphs of Characteristic lines
There are now two interesting plots to consider. This first is a plot of the initial distribution at different times. We see how the initial distribution moves forward in part and holds back, in the process destroying the property of being a function.
* plot({[X(0,x),f3(x),x=0..2],[X(.01,x),f3(x),x=0..2],
[X(.02,x),f3(x),x=0..2],[X(.03,x),f3(x),x=0..2],
[X(.04,x),f3(x),x=0..2]});
Parametric graphs of {X(t,x),f3(x)} for selected t's
Finally, here is the graph of u(t,x) drawn parametrically. One could have predicted how the graph would look from the previous.
* plot3d([t,X(t,x),f3(x)],t=0..0.03,x=0..2,
orientation=[-145,45],axes=NORMAL);
3D plot of the multivalued function u(t,x)
* a:='a': b:='b':
* f3:=x-> a + b/(1+(2-x)^2);
* eq1:=int(f3(x),x=0..1);
* eq2:=int(f3(x),x=1..2);
* solve({eq1=3,eq2=8},{a,b});
* assign(");
* plot(f3(x),x=0..2);
* evalf(a);evalf(b);
* a:=-47/100; b:=108/10;
* evalf(eq1); evalf(eq2);
* X:=(t,x)->68*(1 - f3(x)/10)*t + x;
* diff(1-f3(x)/10,x);
* solve(X(0,x)=X(2,x),x);
* plot({[t,X(t,0),t=0..0.1],[t,X(t,.5),t=0..0.1],[t,X(t,1),t=0..0.1]});
* plot({[X(0,x),f3(x),x=0..2],[X(.01,x),f3(x),x=0..2],
[X(.02,x),f3(x),x=0..2],[X(.03,x),f3(x),x=0..2],
[X(.04,x),f3(x),x=0..2]});
* plot3d([t,X(t,x),f3(x)], t=0..0.03, x=0..2,
orientation=[-145,45],axes=NORMAL);
Example 6.5
In this example, we assume that f(u) = u2. If [[xi]] < 0 the characteristics are given by x = [[xi]] + 2t. If 0 < [[xi]] < 1, the characteristics are x = [[xi]] + 2(1-[[xi]])t. And, if 1 <= [[xi]] then the characteristics are constant [[xi]].
Characteristic curves can be sketched. Note that characteristics begin to cross if t > 1/2. One can provide a solution u for t < 1/2.
u(t,x) = BLC{(A(1 if x < t, (1-x)/(1-2t) if t < x < 1, 0 if 1 < x)).
What happens if t > 1/2 when the characteristics begin to cross? We will see in the next section that the following procedure works for getting the equation for the "shock" in the t,x-plane: re-write the equation in the form
F([[partialdiff]]u,[[partialdiff]]t) + F([[partialdiff]]f,[[partialdiff]]x) = 0.
Form the quotient by taking the numerator to be
f(the limit of the solution from the right) -
f(the limit of the solution from the left) .
Take the denominator to be
the limit of the solution from the right -
the limit of the solution from the left.
In this case, taking u to be 0 on the right of the shock and 1 on the left, we have the slope of the shock = (0 - 1)/(0 - 1) = 1. Further, the shock starts at the overlap of the t = 1/2, x = 1. Hence the shock follows the equation
t = (x - 1)/2.
These ideas will be worked out in the next section.
Exercises
(6.1) Solve these partial differential equations:
a. ut + 3 ux = -u/2, with u(0,x) = [[pi]]/2 - arctan(x).
b. ut + 3 ux = -1/(1+x2), with u(0,x) = 1.
(6.2) Verify that if uy + g(u) ux = 0 with u(0,x) = f(x) then u(t,x) is given implicitly by u(t,x) = f(x - g(u(t,x)). t ).
(6.3) In example 6.3, take v(u) = 1 - u2 and [[phi]](x) as given by equation (6.4). Solve the resulting PDE.