A weighted essentially non-oscillatory numerical scheme for a multi-class LWR model Mengping Zhanga, Chi-Wang Shub, George C.K. Wongc and S.C. Wongc a Department of Mathematics, University of Science and Technology of China, Hefei, Anhui 230026, P.R. China b Division of Applied Mathematics, Brown University, Providence, RI 02912, USA c Department of Civil Engineering, The University of Hong Kong, Pokfulam Road, Hong Kong, P.R. China ABSTRACT In this paper we present a high order weighted essentially non-oscillatory (WENO) scheme for solving a multi-class extension of the Lighthill-Whitham-Richards (LWR) model. We first review the multi-class LWR model and present some of its analytical properties. We then present the WENO schemes, which were originally designed for computational fluid dynamics problems and for solving hyperbolic conservation laws in general, and demonstrate how to apply these to the present model. We found through numerical experiments that the WENO method is vastly more efficient than the low order Lax-Friedrichs scheme, yet both methods converge to the same solution of the physical model. It is especially interesting to observe the small staircases in the solution which are completely missed out, because of the numerical viscosity, if a lower order method is used without a sufficiently refined mesh. To demonstrate the applicability of this new, efficient numerical tool, we study the multi-class model under different parameter regimes and traffic stream models. We consider also the convergence of the multi-class LWR model when the number of classes goes to infinity. We 1 Report Documentation Page Form Approved OMB No. 0704-0188 Public reporting burden for the collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to a penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number. 1. REPORT DATE 3. DATES COVERED 2002 2. REPORT TYPE 00-00-2002 to 00-00-2002 4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER A weighted essentially non-oscillatory numerical scheme for a multi-class 5b. GRANT NUMBER LWR model 5c. PROGRAM ELEMENT NUMBER 6. AUTHOR(S) 5d. PROJECT NUMBER 5e. TASK NUMBER 5f. WORK UNIT NUMBER 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION Brown University,Division of Applied Mathematics,182 George REPORT NUMBER Street,Providence ,RI,02192 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR’S ACRONYM(S) 11. SPONSOR/MONITOR’S REPORT NUMBER(S) 12. DISTRIBUTION/AVAILABILITY STATEMENT Approved for public release; distribution unlimited 13. SUPPLEMENTARY NOTES 14. ABSTRACT 15. SUBJECT TERMS 16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF 18. NUMBER 19a. NAME OF ABSTRACT OF PAGES RESPONSIBLE PERSON a. REPORT b. ABSTRACT c. THIS PAGE 44 unclassified unclassified unclassified Standard Form 298 (Rev. 8-98) Prescribed by ANSI Std Z39-18 show that the solution converges to a smooth profile without staircases when the number of classes increases. Key Words: multi-class LWR model, traffic flow, weighted essentially non-oscillatory scheme, Lax-Friedrichs scheme, Godunov scheme 1. INTRODUCTION Lighthill and Whitham (1955) and Richards (1956) independently proposed a simple continuum model, now known as the LWR model, to describe the characteristics of traffic flow. In this model, a traffic stream model (a relationship between the traffic state variables of flow, speed and density, e.g. Greenshields (1934)) is supplemented by the continuity equation of vehicles, and the resultant partial differential equation presumably could be solved to obtain the density as a function of space and time. For a specific form of Greenshields’ traffic stream model, the solution can be obtained analytically (Wong and Wong, 2002). Although aiming to provide a coarse representation of traffic behavior, the LWR model is capable of reproducing qualitatively a remarkable amount of real traffic phenomena such as shock formation. Nevertheless, there are still some puzzling traffic phenomena observed on the highway, such as the two-capacity or reverse-lambda state in the fundamental diagram, hysteresis of traffic flow and platoon dispersion, that this simple LWR model cannot address or explain. Recently, a multi-class model, extended from the LWR model, with heterogeneous drivers was formulated (MCLWR model) (Wong and Wong, 2001). Let there be M classes of road users with different speed choice behaviors in response to the same traffic density when 2 traveling on a highway section. It means that for a given total density, there exists a distribution of equilibrium speeds by different user classes. It is expected that the variation around the mean speed (averaged over all user classes) decreases when traffic density increases, due to the tighter interactions between road users. Let q (x,t), k (x,t) and m m u (x,t) be, respectively, the flow, density and speed of user class m in the space-time m domain. The total density on a highway section can then be obtained as M (cid:1) k(x,t) (cid:1) k (x,t). (1) m m(cid:1)1 The flow, density and speed variables of a particular class are subject to the following definitional relationship, q (x,t)(cid:1)u (x,t)(cid:3)k (x,t), (cid:2)m(cid:1)1,2,(cid:1),M . (2) m m m From the law of conservation of vehicles, each user class should satisfy the following continuity equation, (cid:3)k (x,t) (cid:3)q (x,t) m (cid:4) m (cid:1)0, (cid:2)m(cid:1)1,2,(cid:1),M , (3) (cid:3)t (cid:3)x which describes the conservation of vehicles at any location at any time along a topographically homogeneous highway section without intermediate entrances or exits. 3 The core of the present extension is to assume that the choice of speed of a particular user class is not only affected by the presence of this user class, but also by all other user classes on the highway. A general form of speed-density relationship can be written as u (x,t) (cid:1)U (k ,k ,(cid:1),k ), (cid:2)m (cid:1)1,2,(cid:1),M . (4) m m 1 2 M For the isotropic case, the above relationship would take a simpler functional form as u (x,t) (cid:1)U (k), (cid:2)m (cid:1)1,2,(cid:1),M , (5) m m where k is the total density determined by equation (1). Combining the above equations, the problem can be formulated into a set of partial differential equations, (cid:3)k (x,t) (cid:1)M (cid:3)k (x,t) m (cid:4) c (x,t) n (cid:1)0, (cid:2)m(cid:1)1,2,(cid:1),M , (6) (cid:3)t mn (cid:3)x n(cid:1)1 where (cid:3)U c (cid:1)U (cid:5) (cid:4)k m , (cid:2)m,n(cid:1)1,2,(cid:1),M , (7) mn m mn m (cid:3)k n is the kinematic wave speed of user class m in response to the presence of class n users, and (cid:4) (cid:1)1 if m = n; and (cid:2) (cid:1)0 if m (cid:1) n. Note that the problem stipulated in equation (6) mn mn reduces to the original LWR model when M = 1 (i.e. homogeneous users). The problem becomes one of solving the set of differential equations (6), or better still the conservation 4 form (3) with q defined by equation (2) and u defined by equation (5), subject to certain m m initial spatial and time boundary conditions. Although the problem formulation is straight forward, it was found that the model is capable of producing the desired properties of a macroscopic traffic flow model and it explains many puzzling phenomena, such as the two- capacity or reverse-lambda state, hysteresis, and platoon dispersion, but it would not be subject to other deficiencies such as wrong-way travel (Daganzo, 1995). In Wong and Wong (2001), the MCLWR model was solved by a first order Lax-Friedrichs finite difference scheme. Although this finite difference scheme is commonly used to solve the original LWR model (Lebacque, 1984; Michalopoulos et al., 1984), it is argued that this first order Lax-Friedrichs scheme may produce smeared solutions near discontinuities due to excessive numerical viscosity. The effect of numerical viscosity will diminish with mesh refinement, but it will be very costly to solve a very refined mesh. More recently, Lebacque (1996) successfully applied the Godunov scheme, introduced by Godunov (1959), to solve the LWR model. The Godunov scheme is subject to smaller numerical viscosity, but it requires a Riemann solver as its building block, which is very difficult, if not impossible, to develop for the MCLWR model. This is because the multi-class model does not seem to be either genuinely nonlinear or linearly degenerate (LeVeque, 1992). Nevertheless, it is important to note that, even though for first order methods the Godunov scheme is more accurate than the Lax-Friedrichs scheme, this difference diminishes dramatically when higher order schemes are considered (Shu, 1998). Both Godunov and Lax-Friedrichs schemes converge to the same physical solution of the model with a sufficiently refined mesh. This can be proved for the scalar and some system cases, and can be observed for more complex systems (LeVeque, 1992; Shu, 1998). 5 This paper presents the solution of the MCLWR model by a weighted essentially non- oscillatory (WENO) scheme (Jiang and Shu, 1996). The WENO scheme is a very robust numerical scheme and is found to be very useful in computational fluid dynamics as well as in other applications. The numerical results from WENO are compared with those obtained from the first order Lax-Friedrichs method. In the special case when all the eigenvalues of the kinematic wave matrix of the system (6) are positive, the Godunov solver becomes the simple upwind solver. In this special case we have verified that the first order Godunov solver converges faster than the first order Lax-Friedrichs solver, but slower than the fifth order WENO solver, while all three converge to the same physical solution. In Section 2, some analytical properties of the MCLWR model are presented for a 2-class model. The WENO scheme for the MCLWR is given in Section 3. Section 4 compares the convergence characteristics of the numerical schemes, shows the numerical solutions for different congestion regimes, and studies the asymptotic case of an infinite number of classes. 2. SOME ANALYTICAL PROPERTIES OF THE MCLWR MODEL 2.1 Hyperbolicity of the system Let k (cid:1)Col(k ,k ,(cid:1),k ) be a column vector containing all M density variables in the 1 2 M MCLWR model. For a smooth solution k (meaning that k has at least first order continuous derivatives in x and t), the set of partial differential equations (PDEs) (3) can be rewritten as (cid:2)k (cid:2)k (cid:3)A(k) (cid:1)0, (8) (cid:2)t (cid:2)x 6 where A(k)(cid:2)(cid:1) q(k) is a kinematic wave matrix containing all the elements of c in k mn equation (7), and q(cid:2)Col(q ,q ,(cid:1),q ) is a column vector of the flow fluxes of all classes. 1 2 M However, when the solution k becomes discontinuous (containing shocks or other discontinuities), the two systems (3) and (8) are not equivalent. Thus to be on the safe side, one should always use conservative schemes to solve (3) directly. The system (3) is called hyperbolic if the eigenvalues of the kinematic wave matrix A(k) are all real and there is a complete, linearly independent set of eigenvectors. Hyperbolic systems are mathematically well-posed, meaning that their solutions depend continuously on the initial conditions. This can be proved for the linear case and also for some nonlinear cases. An important issue in verifying the reasonableness of a model is to check if it is hyperbolic. We confirm that the system (3) is hyperbolic for the practical choices of traffic stream models and their parameters. This verification is performed analytically for the 2-class case and numerically for the general M-class case to be discussed in Section 4. For the 2-class case (M=2), the kinematic wave matrix is a 2(cid:1)2 matrix given by (cid:6)U (k)(cid:8)k U(cid:7)(k) k U(cid:7)(k) (cid:3) A (cid:9) (cid:4)(cid:5) 1 k U(cid:7)1(k)1 U (k)1(cid:8)1k U(cid:7)(k)(cid:1)(cid:2), (9) 2 2 2 2 2 where U (k) and U (k) are defined in equation (5). The two eigenvalues of the kinematic 1 2 wave matrix are thus given by (cid:1) (cid:2) (cid:8) (cid:1) U (k)(cid:7)k U(cid:6)(k)(cid:7)U (k)(cid:7)k U(cid:6)(k)(cid:5) D 2, (10) 1,2 1 1 1 2 2 2 where 7 (cid:1)(cid:1) (cid:2) (cid:1) (cid:2)(cid:2) D(cid:4) U (k)(cid:1)kU(cid:2)(k) (cid:3) U (k)(cid:1)k U(cid:2)(k) 2 (cid:1)4k k U (k)U (k). (11) 1 1 1 2 2 2 1 2 1 2 Clearly D > 0 which implies that the two eigenvalues are both real and distinct. Thus the system is hyperbolic. We could obtain analytical formulas for the eigenvalues of the kinematic wave matrix for the 3-class or even the 4-class case (M = 3 or 4), however the formulas are quite complex and it is not easy to see whether the eigenvalues are always real. At any rate, this approach would not work for the multi-class case with M > 4, as no analytical formulas for the eigenvalues would be available. We thus resort to implementing a numerical eigenvalue solver to verify a posteriori that the eigenvalues are always real for all the test cases in Wong and Wong (2001) and in this paper. Indeed through all the numerical tests, non-real eigenvalues never appear for the kinematic wave matrix. Although this is not a rigorous proof that the MCLWR model is always hyperbolic, it at least gives validity to the numerical experiments in Wong and Wong (2001) and in this paper since the models under all these cases are hyperbolic. 2.2 First order traveling waves In this section, we apply a linearization approach to demonstrate the traveling wave properties of the MCLWR model (Whitham, 1974). To simplify the analysis, we also consider the simple 2-class system and assume a modified Greenshields’ form of traffic stream model, 8 u (cid:9)u (cid:4)(cid:6)1(cid:8) k1(cid:7)k2 (cid:1)(cid:3) and u (cid:9)u (cid:4)(cid:6)1(cid:8) k1(cid:7)k2 (cid:1)(cid:3), (12) 1 f1(cid:4) (cid:1) 2 f2(cid:4) (cid:1) (cid:5) kjam (cid:2) (cid:5) kjam (cid:2) where, for Class 1 and Class 2 traffic, u and u are the traffic speeds, k and k are the 1 2 1 2 densities, and u and u are the free-flowing speeds, respectively, while k is the jam f1 f2 jam density of the highway. The system of differential equations can be written as (cid:1) (cid:2) (cid:1) (cid:2) (cid:2)k (cid:2) u k (cid:2)k (cid:2) u k 1 (cid:3) 1 1 (cid:1) 0 and 2 (cid:3) 2 2 (cid:1) 0. (13) (cid:2)t (cid:2)x (cid:2)t (cid:2)x For small perturbations of densities, r and w, around the steady state densities k and k 1 2 respectively, we can write k (cid:2)k (cid:1)r and k (cid:2)k (cid:1)w. (14) 1 1 2 2 Substituting equation (14) into equation (13) and neglecting higher order terms, we can show that (cid:1)r u (cid:1)r u k (cid:1)w (cid:4) f1 (k (cid:3)2k (cid:3)k ) (cid:2) f1 1 (15) (cid:1)t k jam 1 2 (cid:1)x k (cid:1)x jam jam and (cid:1)w u (cid:1)w u (cid:1)r (cid:4) f2 (k (cid:3)2k (cid:3)k ) (cid:2) f2 k . (16) (cid:1)t k jam 2 1 (cid:1)x k 2 (cid:1)x jam jam Eliminating w from equations (15, 16), we have 9