Having been studied carefully the space and time mesh sizes and patterns or schematic diagrams of all these schemes, another but a new scheme has been developed and named as one-dimensional explicit Tolesa numerical scheme. This new scheme has been visualized with the unique schematic diagram and is illustrated in Figure 5. The implication of the advection Equation (1) together with initial value problem in the context of Lighthill Whitham and Richards (LWR) traffic flow model is implemented. The LWR model is the well-known model and describes the traffic flow using a partial differential equation based on the conservation law of the vehicles in traffic [7] [8] [9] [10] . In this study, the single lane traffic flow with constant speed and linear density-speed relationship are considered as a test case. The simulations of traffic flow model with constant speed will also be done by using the mentioned one-dimensional explicit numerical schemes including Tolesa scheme, and linear density-speed relationship will be simulated by using Tolesa scheme.

In this research work we present and discuss some of the one-dimensional explicit numerical schemes available in the literature and the newly propose scheme. Specifically, the paper is organized in three more sections that follow this Introduction. Section 2 is devoted to carry out finite difference method. The finite difference approximations of the first order hyperbolic partial differential equation using one-dimensional explicit numerical schemes are presented. Section 3 reports about macroscopic continuum traffic flow depend mainly on three quantities flux, speed and density, and present some cases for speed-density relationship. Finally section 4 deals with the numerical simulation. In this section, we present the discretization of the macroscopic continuum traffic flow model, specifically Lighthill-Whitham and Richards Traffic flow model presented in section 3, using finite difference schemes presented in Section 2.

2. Overviews of Finite Difference Method

In this study, it is mainly focused to carry out finite difference method (FDM). The FDM was introduced by Euler in 18^{th} century and has been greatly regarded as the easiest method and widely used to solve simple geometrical problems [11] . The FDM is classically obtained by approximating the derivatives appearing in the partial differential equation by a Taylor expansion up to some given order which will give the order of the scheme. The application of FDM in solving a PDE is to transform a calculus problem into an algebraic problem by discretizing the continuous physical domain into a number of cells or intervals, and approximating the individual exact partial derivatives in the PDE by algebraic finite difference approximations.

Let the spatial and temporal domains be divided into N and M cells respectively. The index $\left({x}_{j},{t}^{n}\right)$ represents a grid point where ${x}_{j}$ spatial and ${t}^{n}$ temporal grids intersect. Here $j=0,1,2,\cdots ,N$ and $n=0,1,2,\cdots ,M$. Thus, the discretization of space-time domain is assumed as ${x}_{j}=j\Delta x$ and ${t}^{n}=n\Delta t$. The quantities $\Delta x>0$ and $\Delta t>0$ respectively denote the spatial and temporal step sizes. These are also known as the increments between two consecutive spatial and temporal nodes respectively.

Let the quantity ${u}_{j}^{n}$ denote an approximate value of the function $u\left(x,t\right)$ at the grid point or space-time location $\left({x}_{j},{t}^{n}\right)$. That is ${u}_{j}^{n}\approx u\left({x}_{j},{t}^{n}\right)$. Thus, the function $u\left(x,t\right)$ can now be replaced with a discrete set of point-wise approximate values $\left\{{u}_{j}^{n}\right\}$. These approximations are called the finite difference approximations. The finite difference approximations of the first order hyperbolic partial differential Equation (1) using one-dimensional explicit numerical schemes are presented in the following subsection.

One Dimensional Explicit Numerical Schemes

In this subsection, the one-dimensional explicit numerical schemes with Upwind, Downwind, FTCS and Lax-Friedrichs schemes including Tolesa scheme are presented to approximate the spatial and temporal partial derivatives of the advection Equation (1) in different views as shown in Figures 1-5. Also the orders of accuracy and numerical stability have been determined in which the space and time mesh sizes are restricted. The approximations and the schematic diagrams of the mentioned schemes are explained sequentially.

In one-dimensional explicit upwind scheme, the first order spatial and temporal partial derivatives are approximated respectively as $\left(\partial u/\partial x\right)\approx \left[\left({u}_{j}^{n}-{u}_{j-1}^{n}\right)/\Delta x\right]+O\left(\Delta x\right)$ and $\left(\partial u/\partial t\right)\approx \left[\left({u}_{j}^{n+1}-{u}_{j}^{n}\right)/\Delta t\right]+O\left(\Delta t\right)$. In this view, the stencil is given as shown in Figure 1 and the approximations of the finite difference form of Equation (1) can be expressed as

${u}_{j}^{n+1}={u}_{j}^{n}-v\left(\Delta t/\Delta x\right)\left({u}_{j}^{n}-{u}_{j-1}^{n}\right)+O\left(\Delta t,\Delta x\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}v>0$ (2)

The scheme (2) is a first order accurate in both space and time, $O\left(\Delta t,\Delta x\right)$. This scheme is stable if the condition $v\left(\Delta t/\Delta x\right)\le 1$ is satisfied.

In one-dimensional explicit downwind scheme, the first order spatial and temporal partial derivative are approximated respectively as

$\left(\partial u/\partial x\right)\approx \left[\left({u}_{j}^{n}-{u}_{j+1}^{n}\right)/\Delta x\right]+O\left(\Delta x\right)$ and $\left(\partial u/\partial t\right)\approx \left[\left({u}_{j}^{n+1}-{u}_{j}^{n}\right)/\Delta t\right]+O\left(\Delta t\right)$.

In view of this the finite difference approximations form of Equation (1) can be expressed as

${u}_{j}^{n+1}={u}_{j}^{n}-v\left(\Delta t/\Delta x\right)\left({u}_{j+1}^{n}-{u}_{j}^{n}\right)+O\left(\Delta t,\Delta x\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}v<0$ (3)

Equation (3) is a first order accurate in both space and time $O\left(\Delta t,\Delta x\right)$. This scheme is unconditionally unstable.

Figure 1. Stencil of One-dimensional explicit upwind numerical scheme.

Figure 2.Stencil of one-dimensional explicit downwind scheme.

In one-dimensional explicit FTCS scheme, as its name implies, the first order temporal and spatial partial derivatives are obtained respectively by taking a first order forward finite differencing in time and a central differencing in space. These are approximated respectively as $\left(\partial u/\partial t\right)\approx \left[\left({u}_{j}^{n+1}-{u}_{j}^{n}\right)/\Delta t\right]+O\left(\Delta t\right)$ and $\left(\partial u/\partial x\right)\approx \left[\left({u}_{j+1}^{n}-{u}_{j-1}^{n}\right)/2\Delta x\right]+O\left(\Delta {x}^{2}\right)$. Then, the approximations of the finite difference form of Equation (1) can be expressed as in Equation (4) and the schematic diagram is given in Figure 3.

${u}_{j}^{n+1}={u}_{j}^{n}-\frac{v\Delta t}{2\Delta x}\left({u}_{j+1}^{n}-{u}_{j-1}^{n}\right)+O\left(\Delta t,\Delta {x}^{2}\right)$ (4)

Equation (4) is a first order accurate in time and second order accurate in space, i.e., $O\left(\Delta t,\Delta {x}^{2}\right)$. This scheme is unconditionally unstable.

In one-dimensional explicit Lax-Friedrichs scheme, the term ${u}_{j}^{n}$ in FTCS scheme is replaced by its average value $\frac{1}{2}\left({u}_{j+1}^{n}+{u}_{j-1}^{n}\right)$. In Lax-Friedrich scheme the finite difference form of Equation (1) can be expressed as (5) and its stencil is viewed as in Figure 4.

${u}_{j}^{n+1}=\frac{1}{2}\left({u}_{j+1}^{n}+{u}_{j-1}^{n}\right)-v\left(\Delta t/2\Delta x\right)\left({u}_{j+1}^{n}-{u}_{j-1}^{n}\right)+O\left(\Delta t,\Delta {x}^{2}\right)$ (5)

Equation (5) is a first order accurate in time and second order accurate in space, i.e., $O\left(\Delta t,\Delta {x}^{2}\right)$. This scheme is conditionally stable, if the condition $v\left(\Delta t/\Delta x\right)\le 1$ is satisfied.

Having been studied carefully the space and time mesh sizes, and the patterns of all these schemes a new scheme has been developed and named as one-dimensional explicit Tolesa numerical scheme. The one-dimensional explicit Tolesa numerical scheme is still another alternative numerical scheme to solve advection equation and can be applied to traffic flows model. The schematic diagram of Tolesa evolution scheme is shown in Figure 5 and the application of this scheme to the advection Equation (1) is straightforward. In view of this, the first order temporal and spatial partial derivatives are approximated respectively as Equations (6) and (7).

Figure 3. Stencil of one-dimensional explicit FTCS scheme.

Figure 4. Stencil of one-dimensional explicit Lax-Friedrichs scheme.

Figure 5. Stencil of one-dimensional explicit Tolesa scheme.

$\left(\partial u/\partial t\right)=\frac{{u}_{j}^{n+1}-{u}_{j}^{n+\left(1/2\right)}}{\left(1/2\right)\Delta t}$ (6)

$\left(\partial u/\partial x\right)=\left(\frac{1}{\Delta x}\right)\left[{u}_{j+\left(1/2\right)}^{n+\left(1/2\right)}-{u}_{j-\left(1/2\right)}^{n+\left(1/2\right)}\right]$ (7)

Let the term ${u}_{j}^{n+\left(1/2\right)}$ on the left hand side of (6) be expressed as the average value as

${u}_{j}^{n+\left(1/2\right)}=\frac{1}{2}\left[{u}_{j+\left(1/2\right)}^{n+\left(1/2\right)}+{u}_{j-\left(1/2\right)}^{n+\left(1/2\right)}\right]$ (8)

In (8), the terms on the left hand side can be expanded, as in view of (5), as

${u}_{j+\left(1/2\right)}^{n+\left(1/2\right)}=\frac{1}{2}\left({u}_{j+1}^{n}+{u}_{j}^{n}\right)-v\frac{\Delta t}{2\Delta x}\left({u}_{j+1}^{n}-{u}_{j}^{n}\right)$ (9)

${u}_{j-\left(1/2\right)}^{n+\left(1/2\right)}=\left(1/2\right)\left({u}_{j}^{n}+{u}_{j-1}^{n}\right)-v\left(\Delta t/2\Delta x\right)\left({u}_{j}^{n}-{u}_{j-1}^{n}\right)$ (10)

Using the three expansions (8)-(10) in the two Equations (6)-(7) and on substituting them in (1), the advection equation reduces to the form as

$\begin{array}{c}{u}_{j}^{n+1}=\left(1/4\right)\left({u}_{j+1}^{n}+2{u}_{j}^{n}+{u}_{j-1}^{n}\right)-\alpha \left({u}_{j+1}^{n}-{u}_{j-1}^{n}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+{\alpha}^{2}\left({u}_{j+1}^{n}-2{u}_{j}^{n}+{u}_{j-1}^{n}\right)+O\left(\Delta t,\Delta {x}^{2}\right)\end{array}$ (11)

Here in (11), the notation $\alpha =v\left(\Delta t/2\Delta x\right)$ is used. Thus, (11) is the newly proposed one-dimensional explicit Tolesa scheme. Equation (11) is a first order accurate in time and second order accurate in space, $O\left(\Delta t,\Delta {x}^{2}\right)$.

Local Truncation Error Local truncation error represents the difference between an exact differential equation and its finite difference representation at a point in space and time. Local truncation error provides a basis for comparing local accuracies of various difference schemes. Accordingly the local truncation error for Tolesa scheme (11) will be

${u}_{t}+c{u}_{x}+\frac{\Delta t}{4}{u}_{tt}-\frac{\Delta {x}^{2}}{4\Delta t}{u}_{xx}+\frac{\Delta {t}^{2}}{6}{u}_{ttt}+\frac{c\Delta {x}^{2}}{6}{u}_{xxx}=0$

$\Rightarrow {u}_{t}+c{u}_{x}+O\left(\Delta t,\Delta {x}^{2}\right)=0$

by using Taylor’s expansion. Therefore, Tolesa scheme is first-order in time and second-order in space.

The fundamental properties that every finite difference approximation of a partial differential equation should possess are consistency, convergence and stability.

Consistence: The notion of consistency addresses the problem of whether the finite difference approximation is really representing the partial differential equation. We say that a finite difference approximation is consistent with a differential equation if the finite difference equations converge to the original equation as the time and space grids are refined. Hence, if the truncation error goes to zero as time and space grids are refined we conclude that the scheme is consistent. For the explicit solution to the advection equation, the truncation error is,

${u}_{t}+v{u}_{x}+\frac{\Delta t}{4}{u}_{tt}-\frac{\Delta {x}^{2}}{4\Delta t}{u}_{xx}+\frac{\Delta {t}^{2}}{6}{u}_{ttt}+\frac{v\Delta {x}^{2}}{6}{u}_{xxx}+\text{Higherorder}=0$

${u}_{t}+v{u}_{x}+{e}_{l}+\text{Higherorder}=0$ ,

where

${e}_{l}=\frac{\Delta t}{4}{u}_{tt}-\frac{\Delta {t}^{-1}\Delta {x}^{2}}{4}{u}_{xx}+\frac{\Delta {t}^{2}}{6}{u}_{ttt}+\frac{v\Delta {x}^{2}}{6}{u}_{xxx}+\text{Higherorder}$.

Thus as $\Delta x\to 0$ and $\Delta t\to 0$ , then ${e}_{l}=0$ , hence the Tolesa scheme is consistent with partial differential Equation (1) as long as $\Delta {t}^{-1}\Delta {x}^{2}\to 0$.

Stability Analysis: A finite difference scheme is stable if the scheme do not allows the growth of error in the solution with different time level. Stability analysis is a useful tool for checking validity of a given numerical scheme [12] . There are many approaches to analyze whether a finite difference scheme is stable or unstable. In this paper, we will consider the Von Neumann stability analysis for presented finite difference schemes. The basic idea of this analysis is given by defining the discrete Fourier transform of u as (12). Let it be assumed that the solution can be seen as eigenmodes [13] which at each grid point have the form

${u}_{j}^{n}={\xi}^{n}{\text{e}}^{ipj\Delta x}$ (12)

Here in (12), $\xi =\xi \left(p\right)$ is a complex number dependent on p and it works as an amplification factor; p is a real spatial wave number; $i=\sqrt{-1}$ is an imaginary number. Equation (12) shows the time dependence of a single eigenmode. The differential equations are said to be stable if $\left|\xi \left(p\right)\right|\le 1$.

Also, the Courant-Friedrichs-Lewy CFL criteria for stability say that $\left|\xi \right|\le 1$ if and only if $\left|\frac{v\Delta t}{\Delta x}\right|\le 1$. CFL is necessary condition for stability [14] .

Substituting (12) in (11) and solving the expression for $\xi \left(p\right)$ is obtained:

$\begin{array}{c}{\left|\xi \right|}^{2}=|(\frac{1}{2}\mathrm{cos}\left(p\Delta x\right)-8{\alpha}^{4}\mathrm{cos}\left(p\Delta x\right)+\left(\frac{1}{4}+2{\alpha}^{2}+4{\alpha}^{4}\right){\mathrm{cos}}^{2}\left(p\Delta x\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left(\frac{1}{4}-2{\alpha}^{2}+4{\alpha}^{4}\right))+4{\alpha}^{2}{\mathrm{sin}}^{2}\left(p\Delta x\right)|\end{array}$

From Von-Neumann stability analysis, Tolesa method is stable when ${\left|\xi \right|}^{2}\le 1$.

${\left|\xi \right|}^{2}=\left|\left(\frac{1}{2}-8{\alpha}^{4}\right)\mathrm{cos}\left(p\Delta x\right)+\left(\frac{1}{4}-2{\alpha}^{2}+4{\alpha}^{4}\right){\mathrm{cos}}^{2}\left(p\Delta x\right)+\left(\frac{1}{4}+2{\alpha}^{2}+4{\alpha}^{4}\right)\right|\le 1$

After some computations, and considering different cases which satisfy equality and the CFL condition, then we decided that the $\alpha \le 1/2$.

$\alpha =\left|\frac{v\Delta t}{2\Delta x}\right|\le 1/2\Rightarrow \left|v\right|\frac{\Delta t}{\Delta x}\le 1$ (13)

Thus, the explicit Tolesa scheme is conditionally stable (13).

Convergence: A numerical scheme is convergent if the computed solution of the discretized equation leads to the exact solution of the differential equation as the time and grid spacing lead to zero. The computed solution ${u}_{j}^{n}$ must approach the exact solution ${u}_{exact}$ of the differential equation at any point ${x}_{j}=j\Delta x$ and ${t}^{n}=n\Delta t$ when $\Delta x$ and $\Delta t$ lead to zero while keeping ${x}_{j}$ and ${t}^{n}$ constant. In other hand, the error

${e}_{j}^{n}={u}_{j}^{n}-{\left({u}_{exact}\right)}_{j}^{n}$

satisfying the following convergence condition.

${\mathrm{lim}}_{\Delta x,\Delta t\to 0}\left|{e}_{j}^{n}\right|\to 0$ , at fixed ${x}_{j}=j\Delta x$ and ${t}^{n}=n\Delta t$. Hence, the explicit Tolesa scheme is convergent; we can see this from the simulations result.

3. Lighthill-Whitham and Richards Traffic Flow Model (LWR Model)

In this section, macroscopic continuum traffic flow model is introduced and analyzed. Macroscopic continuum traffic flow depends mainly on three quantities: traffic density, traffic flow or flux and traffic velocity [15] . The number of vehicles on a highway per unit length is defined as traffic density and is denoted by $\rho \left(x,t\right)$. The traffic flow rate or flux is defined as the number of vehicles passing through a given point x at time t and is denoted by $q\left(x,t\right)$. Here $x\in R$ and $t\in \left[0,\infty \right)$. In this study, highway is considered as a unidirectional roadway of finite length with no entrances and exits.

The well-known Lighthill-Whitham and Richards (LWR) model describes the traffic flow using a partial differential equation constructed based on the conservation law of the vehicles in traffic. In this model the traffic flow is represented using a first order hyperbolic partial differential equation and is put as

$\frac{\partial \rho \left(x,t\right)}{\partial t}+\frac{\partial q\left(x,t\right)}{\partial x}=0$ (14)

The flux can also be expressed in terms of the traffic density and the traffic speed as

$q\left(x,t\right)=v\left(x,t\right)\rho \left(x,t\right)$ (15)

In view of (15), the traffic flow model (14) with initial condition takes the form as

$\{\begin{array}{l}\frac{\partial \rho}{\partial t}+\frac{\partial \left(v\rho \right)}{\partial x}=0\\ \rho \left(x,{t}_{0}\right)={\rho}_{0}\left(x\right)\end{array}$ (16)

Equation (16) is called an initial value problem IVP of the macroscopic traffic flow model.

In this study, Equation (16) has been considered in two different cases depending on the speed-density relationship as constant speed, and linear speed-density relationship. Also assuming that traffic flux and speed are expressed as a function of density $q=q\left(\rho \right)$ and $v=v\left(\rho \right)$ respectively.

Case 1: In this case Equation (16) can be expressed as (17) by considering constant speed v and the flux as a function of density $q\left(\rho \right)$.

$\{\begin{array}{l}\frac{\partial \rho}{\partial t}+v\frac{\partial \rho}{\partial x}=0\\ \rho \left(x,{t}_{0}\right)={\rho}_{0}\left(x\right)\end{array}$ (17)

The analytical solution of the form (17) has been calculated using the method of characteristics in implicit form [16] as follows.

$\rho \left(x,t\right)={\rho}_{0}\left(x-vt\right)$ (18)

Case 2: Linear speed-density function $v\left(\rho \right)={v}_{\mathrm{max}}\left(1-\frac{\rho}{{\rho}_{\mathrm{max}}}\right)$ in 1935, Greenshields [17] proposed what was perhaps the first traffic flow model. According to his observations made using photographic methods, Greenshields postulated that there existed a linear relationship between speed and density. Then traffic flux $q\left(\rho \right)=\left({v}_{\mathrm{max}}\left(\rho -\frac{{\rho}^{2}}{{\rho}_{\mathrm{max}}}\right)\right)$. In this case Equation (16) can be expressed as