1 Introduction

Following the seminal works of ^[1] and ^[2],the higher-order model has become the focus of research on traffic flow problems (see also ^{[3, 4, 5, 6, 7, 8]}). The model has been used to describe the behaviors of averaged quantities,such as the density ρ(x,t) and the average velocity v(x,t) of vehicles at location x and time t. Viewed as a continuum,traffic flow can be described by the equations of continuity and acceleration. The continuity equation reads

and the acceleration of higher-order models takes the general form as follows: In (2),the left-hand side denotes the traffic acceleration with the first term expressing the relaxation to the equilibrium state (ρ,v_e(ρ)),the v_x and ρx terms being related to the physical forces of pressure with p′(ρ) > 0 and c₀(ρ) > 0,and the final term representing the “viscosity” with η = η(ρ,v) > 0.

There are extensive literatures on traveling waves within the higher-order model. For the inviscid model (with η(ρ,v) = 0),the weak solution theory was used to analytically determine a shock and a transition layer,which constituted the upstream and downstream fronts, respectively,of a cluster or stop-and-go wave solutions^{[7, 9, 10]}. For the viscid model with a small viscosity η ＜ 1,the boundary layer method was used to demonstrate a similar solution profile^{[11, 12]},which became to be identical to that in the corresponding inviscid model when the coefficient η vanished^[10]. In another development,the long-wavelength analysis^{[13, 14, 15, 16, 17, 18, 19]} and the Taylor expansion^{[20, 21, 22, 23]} were used to derive the Korteweg-de Vries equation (KdV) equation by assuming/adding a sufficiently large viscosity in/to the car-following or the higher-order model.

The present paper uses the Taylor expansion to seek the traveling wave in a higher-order traffic flow model with large diffusion. Although this issue has been intensively investigated, we significantly improve the KdV solution by assuming a non-constant viscosity coefficient η(ρ,v),which effectively balances the other terms in the expansion. The resultant wave speed of the KdV solution is indicated to be either negative or positive,which corresponds to the first or second characteristic field,respectively. In the previous studies,the KdV wave speed can only be negative and correspond to the first characteristic field as the viscosity coefficient is assumed to be constant^[20]. Moreover,the present paper is the first to discuss the validity of approximation. To compare the analytical and numerical solutions,which is not considered in the previous studies,a finite element scheme is designed to test the accuracy and credibility of the approximate KdV solution.

The rest of this paper is organized as follows. In Section 2,the Taylor expansion is used to derive the approximate KdV solution,and the validity of the approximation is analyzed by the physical boundary and the stability of the solution. In Section 3,a finite element scheme is designed to provide general solutions to the system and,in particular,to test the stability and accuracy of the approximate KdV solution. Section 4 concludes the paper. 2 KdV solution

Let c₀(ρ) be constant,as assumed by ^{[1, 2, 3]}. Then,the conservation forms of (1) and (2) can be uniformly rewritten as

which is applied for discussion in this section and for numerical simulation in Section 3. Here, we denote the flow by q = ρv and the flow-density relationship (the fundamental diagram) by q_e(ρ) = ρv_e(ρ). We note that (4) takes a standard conservation form with viscosity (See ^[10] and the references therein). 2.1 Approximate KdV solution

We consider the traveling wave solutions for the system of (3) and (4): ρ = ρ(z),v = v(z), and q = q(z),where z = x − at,and the wave speed a is constant. With η = η(ρ,v) = μρ(v − a)⁻¹,the substitution gives

where μ is constant. Replacing qz by aρz (see (5)) in (6),we have where λ₁(ρ,q) and λ₂(ρ,q) are the first and second characteristic speeds of the system,respectively, given by We have the following assumptions for the approximation of (7).

Assumption (A1) The flow q is in the vicinity of the equilibrium flow q_e(ρ),namely,there is a △z such that q(ρ) ≈ q_e(ρ(z +△z)),and thus the flow is approximated by

Here,a₁ and a₂ are the parameters depending on △z and ρz,respectively.

Assumption (A2) There exists a constant density ρ such that

and thus for the small perturbation ρ(z),q_e(ρ) is approximated by the first three terms in the Taylor expansion with Apply (9) and (10) of Assumptions (A1) and (A2) to (7). By omitting the higher-order terms and balancing the resultant equation,we have which together with (5),(10),and (11) suggests For the derivation of a KdV solution,we assume a₁(ρ) = 0,which means that a = λ₁(ρ,q_e(ρ)) or a = λ₂(ρ,q_e(ρ)). In this case,(13) can be transformed into still assuming a traveling wave solution U = U(X − αT ). The standard solitary wave solution of (14) reads which,in return,gives an approximate traveling wave solution for the system of (3) and (4) as follows: 2.2 Validity of approximate KdV solution

The relations between (x,t) and (X,T ) are indicated in Table 1 as the main properties of the resultant KdV solution. In Case I,we assume that a₂ < 0,and μ > 0 (see (12)),by which the transformation implies that a = −α < 0. Thus,we can only choose a = λ₁(ρ,q_e(ρ)). In this case,the requirement η = μρ(v − a)⁻¹ > 0 is satisfied,and we can roughly assume that q′ _e(ρ^*) ≈ a < 0. In Case II,we assume that a₂ > 0,and μ < 0,by which the transformation implies that a = α > 0. Moreover,the requirement of η = μρ(v − a)⁻¹ > 0 implies that a > v. Therefore,we can only choose a = λ₂(ρ,q_e(ρ)) and roughly assume q′ _e(ρ^*) ≈ a > 0,which implies q′′ e (ρ^*) < 0 for an applicable fundamental diagram. In either case,we can determine the minimal and maximal densities ρm and ρM of the solution (16). The indicated sign of ρ^* in the table ensures that ρm and ρM are within the physical solution region [0,ρ_jam],and ρ^* is reasonably within (ρ_m,ρ_M). We note that the wave profile is open to the bottom if the amplitude 3a/q′′ e (ρ^*) > 0; otherwise,it is open to the top.

The above arguments give rise to considerable restrictions on ρ^*,which appears to suggest a reasonable approximation. However,as an approximate solution of the present system,the profile of (16) might be unstable. For the space interval [−L,L] with periodic boundary conditions, the traveling wave of (16) that is initially symmetric to x = 0 is supposed to propagate in the interval. In this case,the integral average of (16) is independent of t,which can be easily indicated as

where we have Χ = m (if the profile is open to the bottom) or Χ = M (otherwise). We argue that the stability of the profile is related to the average ρ (see (17)). Consider an equilibrium solution (ρ,q_e(ρ)),which exactly solves the system of (3) and (4) in the same [−L,L] with the same periodic boundary condition and the same unchanged average ρ. It is easy to see that,if the equilibrium solution is stable,then the corresponding profile of (16) is expected to evolve into the former solution,and thus is not stable. In other words,to derive a stable KdV solution profile of (16),the average ρ would essentially ensure an unstable equilibrium solution (ρ,q_e(ρ)). Moreover,the profile of (16) is occupied by a sufficiently long constant region in the ring,which is very close to an equilibrium state (ρ_Χ,q_e(ρ_Χ)). For the stable evolution of the constant part of (16),ρ_Χ would essentially ensure the stability of the equilibrium state (ρ_Χ,q_e(ρ_Χ)). For Case I,assume that (ρ_c1,ρ_c2) is the only unstable region of an equilibrium solution,and the requirements are also listed in Table 1. For Case II,the requirements cannot be satisfied unless we relieve the assumption such that the fundamental diagram q_e(ρ) is (very close to) a linear function in the involved density region. In this case,the approximate KdV solution is very close to a linear wave,and we require that ρ_M < ρ_c1 .

Finally,we note that the two requirements in Case I contradict for L = +∞,which leads to ρ_Χ = ρ. This suggests that,to ensure the two requirements,L would not be too large. 3 Finite element scheme and numerical simulation

For a system like (3) and (4),discontinuous numerical schemes are essential for the resolution of shock profiles when the diffusion is small or vanishes^[10]. These discontinuous schemes are still valid for a large diffusion coefficient η(ρ,v),but they would be inefficient due to a very rigorous Courant-Friedrichs-Lewy condition for numerical stability. In this case,the finite element method is recommended because it is free from the requirement for the resolution of shock and a rigorous stability condition. 3.1 Finite element scheme

The system of (3) and (4) or the like can be rewritten in the following vector form:

For simplicity,we only consider x ∈ ^{[0, 1]}. The interval [−L,L] can be properly transformed by scaling and translation. The discretization of (18) by the element method can be then proceeded by the following steps.

Step 1 Multiply (18) with an arbitrary function w(x) ∈ H,and integrate the resultant equation over ^{[0, 1]}. Then,applying the Green formula to the second and third terms leads to the following variational form of (18):

where q(x) ≡ wε(u)ux −wf(u),H = H ×H,H = {v(x,·)|v(x,·) ∈ H¹(^{[0, 1]}),v(0,·) = v(1,·), v_x(0,·) = v_x(1,·)},and H¹(^{[0, 1]}) is the Sobolev space. It is also required that u ∈ H,and we note that the variational form of (19) will surely lead to a conservative scheme.

Step 2 Divide ^{[0, 1]} into a set of elements Ii = [x_i−1,x_i]|^K _i=1 with △x = x_i − x_i−1. Then, approximate H by the space L,where L is the set of piece-wise linear functions with the following bases:

where i=1,2,· · ·,K−1.

Step 3 Express each component of u as the linear combination of {'i(x)}^K _i=1,which has been granted the periodic boundary conditions u(0,t) = u(1,t) and u_x(0,t) = u_x(1,t) in the sense of distribution. Replace each component of w with {'i(x)} in (19) for i = 1,2,· · ·,K. Then,an ODE system of u_i(t) is derived. The Euler time discretization of the system gives

where c₁ ≡ 1/△x,c₂ ≡ 1/(△x)²,c₃ ≡ 1/(△x)³,and the periodic boundary conditions u₀ = u_K and u_K+1 = u₁ are implied. Moreover, where i=1,2,· · ·,K. 3.2 Numerical simulation

Although the extended finite element scheme can be used to solve the higher-order model with large diffusion,it is only used to examine the accuracy of the approximate KdV solution of (15),which is also analytically described in Table 1. For certainty,we assume

in (4),by which the present system of (3) and (4) is essentially the same as the Kerner- Konh¨auser model^[3]. Here,we assume that the free flow velocity v_f =30 m/s,and the maximal density ρ_jam = 0.2 vehicle/m. In accordance with the analysis in subsection 2.2,the interval [−L,L] (L=2 500 m) and the periodic boundary conditions are also assumed for computation. Moreover,we assume μ=0.1v_f²L in (4).

An equilibrium solution (ρ₀,v_e(ρ₀)) is linearly stable if and only if the kinetic speed q′ _e(ρ₀) is between the two characteristic speeds,λ₁(ρ₀,v_e(ρ₀)) and λ₂(ρ₀,v_e(ρ₀)) (i.e.,(8)),which turns out to be

Accordingly,we can solve the two critical densities ρ_c1 and ρ_c2 in Table 1^[2]. Under the above assumptions,the critical densities ρ_c1 and ρ_c2 and the valid range of ρ (as indicated in Table 1) only depend on c₀.

The present study indicates that the requirements in Table 1 allow for very narrow ranges of c₀ and ρ. Choosing c₀ = 0.71vf and c₀ = 0.60,by which ρ_c1 and ρ_c2 are determined,we derive two “good” solutions which satisfy the requirements for Cases I(a) and I(b) in Table 1. These parameter values are given in the first and the third lines of Table 2. To test the stability and accuracy of the resultant KdV solution of (16) through numerical simulation,a solution profile symmetric to x= 0 is taken as the initial state. Although it is considerably distorted,the profile of the numerical solution at t=3 600 s becomes steady and is thus very close to a traveling wave. Moreover,as predicted in subsection 2.2,the bottom part falls within the stable region with ρ<ρ_c0 . Figure 1 shows the profiles of the KdV and numerical solutions,in which the density is scaled by ρjam,and the interval [−L,L] is transformed into ^{[0, 1]}.

We repeat the above process,only slightly changing the value of ρ such that all the requirements in Table 1 are satisfied except that ρ falls within a stable region with ρ < ρ_c1 or ρ > ρ_c2 . These two situations correspond to the second and fourth lines of Table 2,and the KdV solution profile and the developed solution profile are shown in Fig. 2. As predicted in subsection 2.2, the KdV solution profile appears to converge with an equilibrium constant solution.

4 Conclusions

large and variable viscosity coefficient in the higher-order traffic flow model is assumed to derive an approximate KdV solution. By the validity and stability of the approximate solution, considerable restrictions on the related solution parameters can be prescribed,as shown in Table 1. The present study indicates that the restrictions give rise to a very narrow range of solution parameters (e.g.,ρ) or even model parameters (e.g.,c₀). Besides other restrictions, the restriction on the average density ρ of the approximate KdV solution is also proven to be essential in the analysis in subsection 2.2 and the numerical simulation in subsection 3.2,in which the KdV solution profile becomes unstable when ρ falls within an unstable region of the corresponding equilibrium solution. Finally,the analysis and the extended finite element scheme can be used to study similar approximations in higher-order traffic flow models.