Numerical Schemes For The Traffic Model

Chapter 3 Numerical schemes for the traffic model In this chapter, we present a finite difference scheme of the first order partial differential equation (PDE) based on [10], [12]. The finite difference scheme was devised and first tested by L. F. Pichardson and it was later improved by H. Liebmann. In this scheme the partial differential coefficient is replaced by their finite difference approximations. On the basis of this study, we present the discretization of the macroscopic traffic model by finite difference formula which leads to formulate an explicit upwind difference scheme (E.U.D.S) for the numerical solution of the traffic model as a first order quasilinear PDE. We also establish the stability condition and well-posed-ness of the explicit scheme for traffic flow model.

3.1 Linear Finite Difference schemes 3.1.1 Basics of discretization We present finite difference approximations to obtain solutions of general nonlinear scalar conservation laws

∂ρ ∂f ( ρ ) + =0 ∂t ∂x

for x ∈ℜ, t > 0

− − − −(3.1)

In the beginning, we suppose that the characteristics speed satisfy λ = f ′( ρ ) > 0 for all states that occur in the solution of (3.1). Thus we also provide initial data

ρ ( x,0 ) = ρ 0 ( x )

− − − −(3.2)

and data at the left-hand boundary

ρ ( a, t ) = ρ a ( t )

− − − −(3.3)

The equations (3.1), (3.2) and (3.3) lead to an initial boundary value problem (IBVP) and the IBVP is well-posed if λ = f ′( ρ ) > 0 . In order to approximate solution to (3.1), we will discretize spatial domain ( a, b )

a=x

−

1 2

< x 1 < ... < x 2

l−

3 2

<x

l−

1 2

=b

and the temporal domain ( 0, T ) 0 = t 0 < t 1 < . . . < t N −1 < t N = T

Numerical schemes for the traffic flow model 23

We will define the computational grid cells to be the intervals

spatial grid size,

∆x i ≡ x

i+

1 2

−x

i−

 x 1 , x 1  i− i+  2 2

 with  

the

1 2

1

We also define the temporal grid size, ∆t n + 2 ≡ t n +1 − t n We will integrate the differential equation (3.1) over the space-time rectangle  x 1 , x 1  i− i+ 2 2 

 × t n , t n +1  

(

)

for an illustration of spatial and temporal discretization figure-

3.1,

x

x

i+

1 2

we

x

i+

1 2

t       f  ρ  x 1 , t  dt − ∫ f  ρ  x 1 , t  dt     i− i+ tn   2    2 

t n +1

n +1

∫ ρ ( x, t )dx = ∫ ρ ( x, t )dx + ∫

i−

n +1

n

x

1 2

i−

tn

1

get,

2

No approximations involve with this equation.

x

1 − 2

x

x

1 i− 2

x

1 i+ 2

a

1 l− 2

b

T

t

N

t n +1 tn

0

Fig 3.1: Spatial and temporal discretization

Numerical schemes for the traffic flow model 24

We will discretize the conservation law by approximating the conserved densities x

1 ρ ≈ ∆x i n i

f

and the flux functions

n+ i+

1 2

1 2

1

i+

∫ ρ ( x, t )dx 2

n

x

i−

1 2

t n +1

1

≈

n+

∆t

1 2

∫

tn

   f  ρ  x 1 , t  dt   i+   2 

We will require our discretization to satisfy the conservative difference n+

1 n+   n + 12 ∆t 2 f − f  i+ 1 1  i − ∆x i  2 2 

ρ in +1 = ρ in − x

Initially, we define

1 ρ = ∆x i 0 i

i+

f

−

1 2

1 2

=

1 2

∫ ρ ( x ) dx

i−

1 2

x

−

1 2

=a

we define

t n +1

1 ∆t

− − − − (3.1)

− − − −(3.2)

0

x

Further, at the left-hand boundary n+

1 2

n+

∫ f ( ρ ( t ) )dt a

1 2

tn

3.1.2 Explicit schemes for linear advection equation To demonstrate the scheme, we consider the linear advection equation as follows: ∂ρ ∂ρ +c =0 ∂t ∂x

− − − −(3.3)

where f ( ρ ) = cρ and c > 0 . Then equation (3.3) gives the explicit upwind difference scheme for linear advection equation (3.3) as n+

ρ in +1

c∆t = ρ in − ∆x i

1 2

[ρ

n i

− ρ in−1

]

− − − −(3.4)

Numerical schemes for the traffic flow model 25

3.2 Explicit Schemes for nonlinear equation We consider a nonlinear equation

∂ρ ∂f ( ρ ) + =0 ∂t ∂x

− − − −(3.5)

where is f ( ρ) nonlinear. Then the equation (3.5) provides the explicit upwind difference scheme is of the form n+

ρ in +1

1

[

∆t 2 n = ρ in − f i − f i −n1 ∆x i

]

3.3 Stability of the scheme Courant, Friedrichs and Lewy, in their fundamental paper (1928) on difference methods for partial differential equations, formulated a necessary condition now known as the CFL condition for the convergence of a difference approximation in terms of the concept of a domain of dependence. For the model problem (3.3), we have seen the solution is ρ( x, t ) = ρ 0 ( x − ct ) , where the function ρ0 is determined by the initial

(

n conditions. The solution at the point x i , t

)

is obtained by drawing the characteristic

(

)

n through this point back to where it meets the initial line at Q , x i − ct ,0 - Fig 3.2.

Now the scheme (3.4) can be written as the averaging process

ρ in +1 = (1 − γ ) ρ in + γ ρin−1 n+

where γ i

1 2

n+

c∆t := ∆x i

1 2

− − − −(3.6)

is the dimensionless courant number.

t P

Q

Fig 3.2: Typical domain of dependence Numerical schemes for the traffic flow model 26

x

n +1 The value of ρi depends on the values of ρ at two points on the previous time

level: each of these depends on two points on the time level t n −1 and so on. As

ρin +1 depends on data given in a triangle with

illustrated in Fig. 3.2, the value of

(

)

n +1 vertex x i , t , and ultimately on data at the points on the initial line

x i −n −1 , x i −n , . . . , x i −1 , x i . n +1 So, ρi depends on data given at all points of the triangle. This triangle is called the

(

)

n +1 n +1 domain of dependence of ρi , or of the point x i , t , for this particular numerical

scheme. The corresponding domain of dependence of the partial differential equation is the

(

)

n +1 characteristic path drawn back from x i , t to the initial line.

The CFL condition then states that for a convergent scheme the domain of dependence of the partial differential equation must lie within the domain of dependence of the numerical scheme. It also gives a restriction on the size of the time step, for the condition that the characteristic must lie within the triangle of dependence requires that n+

c∆t ∆x i

1 2

1

n+ ≤ 1 , that is γ i 2 ≤ 1 . The expression on the R.H.S of equation (3.6) would

become a convex combination.

3.4 Error analysis of the scheme We recall the scheme (3.6)

ρ in +1 = (1 − γ ) ρ in + γ ρin−1 ∆t where γ := V max

n+

− − − −(3.7)

1 2

∆x i

This can be interpreted as follows. In fig-3.2, the characteristic through the point

P = ( x i , t n +1 ) meets the previous line t = t n at the point condition

(

)

(

n n must lie between the points A = x i , t and B = x i −1 , t

Numerical schemes for the traffic flow model 27

)

Q , which by the CFL

P

t B

Q

A

C

n +1

n

x

j

Fig 3.3: Error Analysis of the scheme Moreover the exact solution ρ( x, t ) is constant along the characteristic, so that

ρ( P ) = ρ( Q ) . Knowing an approximate numerical solution at all the points on the line

t n , we can therefore interpolate the value of ρ( Q ) and use this to give the required value ρi

n +1

of

x

(

. If we use linear interpolation, approximating ρ x, t n

)

by a linear function

determined by the approximations at the two points A and B , we obtain (3.6)

exactly when V max

is constant because AQ = γ∆x and QB = (1 − γ ) ∆x ; when varies

smoothly this still gives a good approximation. Notice also that all the coefficients in (3.7) are non-negative so that a maximum principle applies, provided that γ

≤1

at all mesh points. We can therefore obtain an

error bound for the linear, variable coefficient problem. We must first consider more carefully what domain is given, and what boundary conditions; although the physical problem may be given on the whole line, for all values of

x , a numerical solution must

be confined to a finite region. Suppose, for example, that the region of interest is

0 ≤ x ≤ X , so that we have boundaries at x = 0 and x = X . Since the partial differential equation is hyperbolic and first order, we will usually have only one boundary condition. The direction of the characteristics shows that we need a boundary condition at x = 0 , we therefore have just one boundary condition. The exact solution of the partial differential equation would when be determined by drawing the

Numerical schemes for the traffic flow model 28

characteristic backward from the point P , until it reaches either the initial line t = 0 , or a boundary on which a boundary condition is given.

(

n The truncation error of the scheme is defined as usual and expansion about x i , t

)

gives, if ρ is sufficiently smooth,

ρ in +1 − ρ in ρ in − ρ in−1 Ti := + Vmax ∆t ∆x n

n

n

 1 1     ~  ρ t + ∆tρ tt + . . . + Vmax  ρ x − ∆xρ xx + . . . 2 2  i    i =

1 ( ∆tρ tt − Vmax ∆xρ xx ) + . . . 2

− − − −(3.8)

Even if V max is constant, we still find T In = −

1 (1 − γ )Vmax ∆xρ xx + . . .; 2

hence generally the method is first order accurate. Suppose the difference scheme is applied for i =1, 2, . . . , I , at the points x i = i∆x with I∆x = X , and the boundary

(

)

n n n n n value Ρ0 = ρ 0, t is given. Then for the error ei = Ρi − ρ i we have as usual

ein +1 = (1 − γ ) ein + γein−1 − ∆tTi n

− − − −(3.9)

n and e 0 = 0 , from which we deduce that if 0 ≤ γ ≤1 at all points

E n +1 := max ein +1 ≤ E n + ∆t max Ti n i

i

If we suppose that the truncation error is bounded, so that Ti n ≤T

for all i and

n

− − − −(3.10 )

in the domain, the usual induction argument shows that

E n ≤ n∆tT ≤ t F T

− − − −(3.11)

0 0 if Ρi = ρ ( x i ) .

This result is sufficient to prove first order convergence of the upwind scheme along a refinement path which satisfies the CFL condition everywhere, provided that the solution has bounded second derivatives.

Numerical schemes for the traffic flow model 29

3.5 Numerical scheme for the traffic flow model Based on the discussion in the previous sections we investigate a finite difference scheme for our considered traffic flow model. As mentioned at section 2.1, traffic model yields

  ρ ∂ρ ∂  + ρ.V max 1 −    ρ max ∂t ∂x   

  

2

  = 0  

− − − −(3.12 )

where x ∈( a, b ) ; t ∈( 0, T ) The initial and boundary conditions are given by

ρ ( x,0 ) = ρ 0 ( x )

and

ρ ( a, t ) = ρ a ( t )

− − − −(3.13)

The equations (3.12) and (3.13) produce an initial boundary value problem (IBVP) which is well-posed for left hand boundary if velocity, V max > 0 and for the right hand boundary if velocity, V max < 0. Equation (3.12) may be written as

∂ρ ∂q( ρ ) + =0 ∂t ∂x

− − − −(3.14)

   ρ  where q( ρ ) =  ρ .V max 1 −    ρ max  

  

2

   

Firstly, we would like to derive an appropriate numerical scheme for the equation (3.14). In the case of linear part of the equation (3.12), f

n+ i+

1 2

1 2

= V max ρ .

In order to develop the scheme, we discretize the space and time. The discretization of

∂ρ( x, t ) is obtained by first order forward difference in time and the discretization of ∂t ∂ρ( x, t ) is obtained by first order backward difference in space. ∂x The possible finite difference approximations for

∂ρ ∂ρ and : ∂t ∂x

Forward difference in time: From Taylor’s series we write

ρ ( x , t + h ) = ρ ( x, t ) + h

∂ρ ( x, t ) h 2 ∂ 2 ρ ( x, t ) + + ... ∂t 2! ∂t 2

Numerical schemes for the traffic flow model 30

⇒ ∴

∂ρ ( x, t ) ρ ( x, t + h ) − ρ ( x, t ) = − o( h ) ∂t h

∂ρ ( x, t ) ρ ( x, t + h ) − ρ ( x, t ) ≈ ∂t h

− − − −(3.15)

Similarly, backward difference in space

∂ρ ( x, t ) ρ ( x, t ) − ρ ( x − k , t ) ≈ ∂x k

− − − −(3.16 )

We assume the uniform grid spacing with step size h and k for time and space respectively t n +1 = t n + h and xi +1 = xi + k . n We also write ρi for ρ( x, t ) in equation (3.15) and (3.16).

Now equation (3.14) takes the form

ρ in +1 − ρ in ∆t

n+

1 2

q in − q in−1 + =0 ∆ xi n+

⇒ ρ in +1

1

[

∆t 2 n = ρ in − q i − q in−1 ∆x i

  ρn   i q = ρ V where max 1 −   ρ   max n i

n i

   

2

]

− − − −(3.17)

   

This is the explicit upwind difference scheme for the equation (3.12). Therefore, equation (3.17) leads the desired numerical scheme for the traffic model.

3.6 Stencil of the scheme Stencil visualizes the flow of the scheme. The stencil for the E.U.D.S (3.4) is presented below:

ρin +1

ρin−1

ρin

Fig 3.4: Stencil of the scheme

Numerical schemes for the traffic flow model 31

3.7 Stability and well-posed-ness of the numerical scheme We recall the traffic flow model as

∂ρ ∂f ( ρ ) + =0 ∂t ∂x



where f ( ρ ) = ρV ( ρ ) = ρVmax 1 −



ρ2 2 ρ max

   

As mentioned at 3.3, we conclude that the stability condition of the scheme for traffic

flow model is γ = V max ∆t

∆x i

n+

1 2

≤1.

 ρ2  Now we have, f ( ρ ) = ρV ( ρ ) = ρVmax 1 − 2 ρ max 

   

 ρ3  = V max  ρ − 2 ρ max 

   

 3ρ 2  ∴ f ′( ρ ) = V max 1 − 2 ρ max 

   

We know that the explicit upwind difference scheme is well-posed iff f ′( ρ ) > 0.

 3ρ 2  V 1 − It is evident that max  2 ρ max 

 >0  

Here, V max is essentially positive.



So, 1 −



⇒

3ρ 2 2 ρ max

 >0  

3ρ 2 <1 2 ρ max

2 ⇒ 3ρ 2 < ρ max

⇒ρ <

1 3

ρ max

⇒ ρmax > 3 ρ

This is the condition of well-posed ness.

Numerical schemes for the traffic flow model 32

With the condition of well-posed ness, we see that

 3ρ 2 f ′( ρ ) = Vmax 1 − 2 ρ max 

  ≤ V max  

The stability condition of the EUDS for the nonlinear problem is given by n+

λ ∆t γ= ∆x i

1 2

f ′( ρ ) ≤ λ

≤ 1,

− − − −(3.18)

Therefore, the stability condition (3.17) becomes

γ=

V max ∆t ∆x i

n+

1 2

≤1

1

i.e. ∆t n + 2 ≤ V ∆x max i with this condition of stability, the condition of well-posed ness can be incorporated by

ρ max = c * max( ρ 0k ) ; k

c≥ 3

Numerical schemes for the traffic flow model 33