Simulaciones numéricas de un modelo macroscópico de tráfico vehicular con relajación

Numerical simulations of a macroscopic traffic Flow model with relaxation

DOI: https://dx.doi.org/10.17981/ingecuc.20.2.2024.18

Artículo de Investigación Científica.

Fecha de Recepción: 02/10/2023, Fecha de Aceptación: 15/12/2023.

Carlos Arturo Vega-Fuentes

Universidad del Norte. Barranquilla, (Colombia)

cvega@uninorte.edu.co

To cite this paper

C. Vega-Fuentes, “Numerical simulations of a macroscopic traffic Flow model with relaxation,” INGE CUC, vol. 20, no. 2, 2024. DOI: https://dx.doi.org/10.17981/ingecuc.20.2.2024.18

INTRODUCCIÓN

Currently, the study of traffic flow models constitutes a very active line of research due to the growing number of vehicles in urban areas and the imperative need to improve traffic conditions. From a mathematical point of view, models to describe and simulate vehicular traffic can be classified into two large groups [1]: microscopic models and macroscopic models. The former describe the behavior (movement) of each vehicle based on the behavior of the vehicle or vehicles in front of it. Although microscopic models provide a very realistic description of traffic flow, they usually require considerable computational effort to obtain numerical results when there is a high volume of traffic because such models involve the individual behavior of each vehicle simultaneously. On the other hand, macroscopic models do not study the individual behavior of each vehicle but rather describe traffic as continuous flow. Macroscopic models were first introduced by Lighthill and Whitham [2] and by Richards [3]. This first model, called the LWR model, is described by the following first order differential equation obtained from the continuity equation for density of vehicles:

${\partial }_t\rho +{\partial }_x\left(\rho v\right)=0$ (1)

where ρ(x,t) and v(x,t) denote, respectively, the local density of vehicles (veh/m) and is the average speed (m/s) at position and at time . The LWR model is useful for characterizing equilibrium traffic flows in which speed is a function of traffic density, i.e., all drivers adjust their speed depending on the spacing between vehicles. However, the LWR model fails to describe other traffic phenomena such as accelerations and decelerations or non-equilibrium conditions, for example, instantaneous speed changes. To overcome these difficulties, the so-called second-order models have been proposed, which are described by a system of two partial differential equations, one of which includes a relaxation term describing drivers’ speed adjustments due to traffic conditions ahead of them. One of the most famous second-order macroscopic models is the Aw-Rascle model [4]-[5]. This model is described by the scalar equation (1) together with an equation describing the mass flow variations due to the road conditions in front of the driver. More precisely, the model in conservative form is described by the following system of partial differential equations

$\left\{\begin{array}{c}{\partial }_t\rho +{\partial }_x\left(\rho v\right)=0\\ {\partial }_t\left(\rho w\right)+{\partial }_x\left(v\rho w\right)=\frac{1}{T}g(\rho ,w)\end{array}\right.$ (2)

where w=v+P(ρ) with P(ρ) a given increasing function describing the anticipation of road conditions in front of the drivers, g(ρ,w) is a relaxation which considers possible entrances or exits on the road and the parameter T is a relaxation time.

By introducing the variable z=ρw (called generalized momentum), the system (2) can be written as the one-dimensional balance law

${\partial }_{t}U+{\partial }_{x}F\left(U\right)=S\left(U\right)$ (3)

where U=(ρ,z) is the vector of conservative variables, F(U)=(ρv,vρw)=(z-ρP(ρ),z2⁄ρ-zP(ρ)) is the flux function and y S(U)=(0,g(ρ,w)⁄T) is the source term. Characteristic speeds are the eigenvalues of the Jacobian matrix of the flux function and are given by

$\begin{array}{c}{\lambda }_1=\frac{z}{\rho }-P\left(\rho \right)-\rho \mathrm{P\text{'}}\left(\rho \right)\\ {\lambda }_2=\frac{z}{\rho }-P\left(\rho \right)\end{array}$ (4)

Since P(ρ) is an increasing function, then the characteristic speeds are distinct and, consequently, the system (3) is hyperbolic.

Since, in general, it is not possible to explicitly find an analytical solution of the system (2), the alternative is to find numerical approximations. However, the development of efficient numerical methods is a complex task due to the presence of the relaxation parameter in the source term given by (2) imposes some degree of stiffness on the system. This implies that when performing the time discretization using standard explicit methods, the time step is very restrictive which affects the efficiency of the method. A very popular alternative that avoids this restriction is the implicit-explicit Runge-Kutta (IMEX-RK) schemes, in which the hyperbolic part is treated explicitly and the source term implicitly [6]. Due to the nonlinearity of the source term in (2), if such a term is treated implicitly, then it is necessary to solve nonlinear algebraic equations at each time step, thus increasing the computational effort. On the other hand, systems of type (3) usually admit steady state solutions in which the balance between the gradient of the flux function and the source term must be maintained. However, most IMEX schemes do not satisfy this requirement [7]. A very attractive technique for temporal discretization to overcome the disadvantages of the aforementioned schemes are the high-order Runge-Kutta exponential methods introduced in [8] and recently improved in [9]. A detailed description of these methods and the schemes used for spatial discretization will be given in the next section.

METHODOLOGY

Since in general it is not possible to find an exact or analytical solution of the system (3), it is imperative to approximate the solution numerically. For this purpose, a discretization of the spatial derivative of the flux function will first be performed using the robust fifth-order weighted essentially non-oscillatory (WENO) method [10]-[11]. This method is especially useful for approximating hyperbolic partial differential equations whose solutions may be discontinuous as is the case of the model considered in this work.

We start by discretizing the computational domain [a,b] where a and b are the ends of the portion of the road considered for the study. For simplicity, we take a uniform mesh, i.e., we perform a regular partition of [a,b] into M cells $[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}]$ with mesh width h = (b - a)/M and centers at the points xj = a + jh, j = 1,2, ..., M. Thus, a semi-discrete scheme for approximating the solution of the system (3) has the form

$\frac{d{U}_j\left(t\right)}{dt}=-\frac{1}{h}\left({\hat{f}}_{j+\frac{1}{2}}-{\hat{f}}_{j-\frac{1}{2}}\right)+S_j$ (5)

Here, Uj(t) denotes the numerical approximation of U(xj, t), ${\hat{f}}_{j+\frac{1}{2}}$ is a numerical flux at the interface $x_{j+\frac{1}{2}}$ and $S_j=\left(0,s_j\right)$ is the discretization of the source term, where

sj=g(ρj,zj )/T.

Note that the first term on the right-hand side in (5) corresponds to the discretization of the spatial derivative ${\partial }_{x}f\left(U\right)$ and such a vector will be denoted by (Lj(U), Μj(U)). To compute ${\hat{f}}_{j+\frac{1}{2}}$ , the widely used fifth order WENO method is used. For the sake of clarity in presentation we next describe in detail the construction of ${\hat{f}}_{j+\frac{1}{2}}$ . Due to nonlinearity of the model, it is necessary to split the flux into positive and negative parts, that is, F(U) = F+ (U)+ F- (U). In this paper we use the Lax-Friedrichs flux splitting $F^{\pm }\left(U\right):=\frac{1}{2}\left(F\left(U\right)\pm \gamma U\right)$ . Here, γ is the so-called viscosity coefficient or maximal characteristic velocity which is calculated from the eigenvalues of the Jacobian matrix given in (4) by using the expression

$\Gamma =\max {}_j\left\{\left|{\lambda }_1\left(U_{j+\frac{1}{2}}\right)\right|,\left|{\lambda }_2\left(U_{j+\frac{1}{2}}\right)\right|\right\}$ (6)

where $U_{j+\frac{1}{2}}=\frac{1}{2}\left(U_j+U_{j+1}\right)$ .

Next, the numerical fluxes ${\hat{f}}_{j+\frac{1}{2}}^{+}$ and ${\hat{f}}_{j+\frac{1}{2}}^{-}$ are constructed by using the smoothness indicators ${\mathrm{IS}}_p^{+}$ y ${\mathrm{IS}}_p^{-}$ for p = 0,1,2, which are given by the following equations

${\mathrm{IS}}_0^{+}=\frac{13}{12}{\left(f_{j-2}^{+}-2f_{j-1}^{+}+f_j^{+}\right)}^2+\frac{1}{4}{\left(f_{j-2}^{+}-4f_{j-1}^{+}+3f_j^{+}\right)}^2$

${\mathrm{IS}}_1^{+}=\frac{13}{12}{\left(f_{j-1}^{+}-2f_j^{+}+f_{j+1}^{+}\right)}^2+\frac{1}{4}{\left(f_{j-1}^{+}-f_{j+1}^{+}\right)}^2$ (7)

${\mathrm{IS}}_2^{+}=\frac{13}{12}{\left(f_j^{+}-2f_{j+1}^{+}+f_{j+2}^{+}\right)}^2+\frac{1}{4}{\left(3f_j^{+}-4f_{j+1}^{+}+f_{j+2}^{+}\right)}^2$

and

${\mathrm{IS}}_0^{-}=\frac{13}{12}{\left(f_{j+1}^{-}-2f_{j+2}^{-}+f_{j+3}^{-}\right)}^2+\frac{1}{4}{\left(3f_{j+1}^{-}-4f_{j+2}^{-}+f_{j+3}^{-}\right)}^2$

${\mathrm{IS}}_1^{-}=\frac{13}{12}{\left(f_j^{-}-2f_{j+1}^{-}+f_{j+2}^{-}\right)}^2+\frac{1}{4}{\left(f_j^{-}-f_{j+2}^{-}\right)}^2$ (8)

${\mathrm{IS}}_2^{-}=\frac{13}{12}{\left(f_{j-1}^{-}-2f_j^{-}+f_{j+1}^{-}\right)}^2+\frac{1}{4}{\left(f_{j-1}^{-}-4f_j^{-}+3f_{j+1}^{-}\right)}^2$

Next, the smoothness indicators are used to compute the weights ${\sigma }_p^{+}$ y ${\sigma }_p^{-}$ for p = 0,1,2:

${\sigma }_p^{+}=\frac{{\theta }_p^{+}}{{\theta }_0^{+}+{\theta }_1^{+}+{\theta }_2^{+}};{\sigma }_p^{-}=\frac{{\theta }_p^{-}}{{\theta }_0^{-}+{\theta }_1^{-}+{\theta }_2^{-}}$ (9)

where

${\theta }_p^{+}=\frac{C_p}{{\left(\delta +{\mathrm{IS}}_p^{+}\right)}^2};{\theta }_p^{-}=\frac{C_p}{{\left(\delta +{\mathrm{IS}}_p^{-}\right)}^2}$ (10)

The coefficients Cp for p = 0,1,2, are the ideal weights [12] which are given by $C_0=\frac{1}{10},C_1=\frac{6}{10},\; and{C}_2=\frac{3}{10}$ . The positive parameter δ was originally introduced to avoid null denominator. However, as it is pointed out in [12], large values of may cause oscillations, while small values may cause an order drop. In this paper, this parameter is taken as δ=10-40 which corresponds to the recommended value in [12].

Finally, the numerical flow is formed by adding the positive and negative parts, that is,

${\hat{f}}_{j+\frac{1}{2}}={\hat{f}}_{j+\frac{1}{2}}^{+}+{\hat{f}}_{j+\frac{1}{2}}^{-}$ (11)

So far only the spatial discretization of (3) has been considered, leaving the time variable continuous. For the discretization of the time derivative in scheme (5) a type of exponential Runge-Kutta method (abbreviated RK-ERK resulting from a combination of a standard third-order Runge-Kutta method and a modified exponential Runge-Kutta method) of third-order accuracy recently proposed in [9] will be used. The description of the method is summarized as follows. The computations of an RK-ERK scheme necessary to advance from Un=(ρn, zn ) which corresponds to the approximate value of the solution at time $t^n=n\Delta t$ to $U^{n+1}=\left({\rho }^{n+1},z^{n+1}\right)$ which corresponds to the approximate value of the solution at time are given by

${\rho }^1={\alpha }_{10}{\rho }^n+{\beta }_{10}\Delta tL({\rho }^n,z^n)$

$z^1=z^n+\frac{{\beta }_{10}\Delta t\left[M({\rho }^n,z^n)+s({\rho }^n,z^n)\right]}{A_1}$

${\rho }^2={\alpha }_{20}{\rho }^n+{\alpha }_{21}{\rho }^1+{\beta }_{21}\Delta tL({\rho }^1,z^1)$

$z^2=\frac{{\alpha }_{20}{z}^n}{A_2}+\frac{e{}^{{\beta }_{10}\mu \Delta t}\left\{{\alpha }_{21}{z}^1+{\beta }_{21}\Delta t\left[M({\rho }^1,z^1)+s({\rho }^1,z^1)+\mu z^1\right]\right\}}{A_2}$

${\rho }^{n+1}={\alpha }_{30}{\rho }^n+{\alpha }_{31}{\rho }^1+{\alpha }_{32}{\rho }^2+{\beta }_{32}\Delta tL({\rho }^2,z^2)$

$z^{n+1}=\frac{{\alpha }_{30}{z}^n}{A_3}+\frac{e^{{\beta }_{10}\mu \Delta t}{\alpha }_{31}{z}^1}{A_3}+\frac{e^{A\mu \Delta t}\left({\alpha }_{32}{z}^2+{\beta }_{32}\Delta t\left[M({\rho }^2,z^2)+s({\rho }^2,z^2)+\mu z^2\right]\right)}{A_3}$

The optimal values of the coefficients are given in Table 1 (more details can be found in [9]).

α10=0	α21=0.3313107066925596	β10=0.7071933376925014
α20=0.6686892933074404	α31=0.2039576138780898	β21=0.4178047564915065
α30=0.3487419430256090	α32=0.4473004430963011	β32=0.5640754637100439

Table 1. Optimal coefficients for the RK-ERK method.

In addition,

$A={\alpha }_{21}{\beta }_{10}+{\beta }_{21},A_1={\alpha }_{10}+{\beta }_{10}\mu \Delta t,A_2={\alpha }_{20}+e^{{\beta }_{10}\mu \Delta t}\left({\alpha }_{20}+{\beta }_{10}\mu \Delta t\right)$

$A_3={\alpha }_{30}+{\alpha }_{31}{e}^{{\beta }_{10}\mu \Delta t}+e^{A\mu \Delta t}\left({\alpha }_{32}+{\beta }_{32}\mu \Delta t\right)$

and the parameter μ is calculated dynamically based on Theorem 2.3 proved in [8], namely

$\mu =\max \left\{\frac{-s({\rho }^n,z^n)}{z^n},0\right\}$ (12)

This choice of guarantees that the scheme satisfies the sign property, i.e., if $U^n\ge 0then{U}^{n+1}\ge 0$ . Such a property is especially important because it avoids unphysical negative densities. in the approximate solution. Another remarkable advantage of the RK-ERK method is that it is not necessary to take a restrictive time step [9], unlike other explicit methods for which the presence of the stiffness source term forces to take a very small time step to guarantee stability. Furthermore, the RK-ERK method preserves the steady state solutions, i.e., if $\left(L_j\left(U^n\right),M_j\left(U^n\right)\right)+\left(0,s_j\left(U^n\right)\right)=0then{U}^{n+1}=U^n$ .

NUMERICAL RESULTS

In this section, three examples are presented to obtain numerical solutions of the model (2) with Riemann-type initial conditions, i.e., piecewise constant initial data. In all examples a spatial mesh with M = 200 is taken. The units of length are given in meters (m) and time in seconds (s).

Example 1.

In this example, taken from [5] the relaxation function is g(ρ,z) = ρ(V(ρ) - v) = ρ(V(ρ) + P(ρ)) - z, where V(p) is a function describing the dependence of the velocity with respect to the density for an equilibrium situation. This equilibrium velocity is decreasing with respect to density and is taken as the following fitting function for experimental data

$V\left(\rho \right)=v_m\cdot \frac{\pi }{2+\arctan \left(\frac{11(\rho -0.22)}{\rho -1}\right)}⁄\frac{\pi }{2+\arctan (11\cdot 0.22)}$ (13)

where vm is the normalized maximum speed. The initial conditions are

$U(x,0)=\left\{\begin{array}{cc}{U}_l=\left({\rho }_l,z_l\right)& six\le x_d\\ U_r=\left({\rho }_r,z_r\right)& six>x_d\end{array}\right.$ (14)

where x(d ) = 0 is the point where the discontinuity is located with left and right states given by Ul = (0.05,-0.2971) and Ur = (0.05,-0.2746) respectively. As in [5], we also take P(ρ) = 2 ln(ρ) and the computational domain [-250,150]. In Fig 1 the evolution of the density and momentum of the traffic flow ρv is shown (first column) along with the profiles of the solutions at t = 1s (second column) with relaxation time T = 0.2s. These results are contrasted with the numerical solutions for the same problem but without relaxation term (T = ∞) presented in Fig 2. In the latter case, the solutions at t = 1s for both density and momentum of the traffic flow consist of a rarefaction wave followed by a contact discontinuity.

Fig. 1. Solutions evolution for density and momentum of the traffic Flow (left) with the corresponding solutions at t = 1 s (right) for example 1 with relaxation.

Source: author.

Example 2.

In this second scenario, the same relaxation function of example 1 is considered with P(ρ)=ρ2 , computational domain [0,20], relaxation time T = 0.3s and initial conditions

$U(x,0)=\left\{\begin{array}{cc}\left(0.5,0.15\right)& \mathrm{si}x\le 10\\ \left(0.15,0.95\right)& \mathrm{si}x>10\end{array}\right.$ (15)

describing the situation of faster vehicles downstream followed by slower ones from behind. The numerical solutions corresponding to the time evolution are shown in Fig 3 along with the density and velocity profiles at t = 8s.

Fig. 2. Solutions evolution for density and momentum of the traffic Flow (left) with the corresponding solutions at t = 1 s (right) for example 1 without relaxation.

Source: author.

Example 3.

For this test case, we take P(ρ)=ρ, the relaxation time is T = 1, the computational domain is the interval and the relaxation function is taken as $g(\rho ,z)=\left(0.1\left(z-\rho P\left(\rho \right)\right)+2\rho \right)\left(\frac{z}{\rho }-9\right)$ The aim of this numerical experiment, whose parameters are taken from [13], is to use the numerical scheme described in this paper to simulate the double Riemann problem with initial data

$U(x,0)=\left\{\begin{array}{cc}{U}_1=\left({\rho }_1,z_1\right)& \mathrm{si}x<0\\ U_2=\left({\rho }_2,z_2\right)& \mathrm{si}0<x<L\\ U_3=\left({\rho }_3,z_3\right)& \mathrm{si}x>L\end{array}\right.$ (16)

where ${\rho }_1=\frac{13}{4},{\rho }_2=7,{\rho }_3=4,z_1=\frac{299}{8},z_2=63and{z}_3=36$ . This situation corresponds to the traffic light problem where in there is a red traffic light which a switches to green. Double Riemann problems are demanding tests for testing numerical methods because their solutions may involve shock waves, contact discontinuities and generalized rarefaction waves. In [13] the exact solution of this problem expressed in terms of the density and velocity was obtained, namely

$\left(\rho (x,t),v(x,t)\right)=\{\begin{array}{cc}\left({\rho }_1,v_1\right)& x<x_s\left(t\right)\\ \left(\hat{\rho }(x,t),\hat{v}(x,t)\right)& x_s\left(t\right)<x\le x_r\left(t\right)\\ \left({\rho }_1,v_1\right)& x\ge x_r\left(t\right)\end{array}$ (17)

Here,${\hat{x}}_s\left(t\right)=\left(v_1-{\rho }_1\right)t-2L\frac{\sqrt{t}}{\sqrt{t_c}}+L,{\hat{x}}_r\left(t\right)=\left(v_3-{\rho }_3\right)t+L,\hat{\rho }(x,t)=0.5\left(v_1+{\rho }_1-\frac{x-L}{t}\right),\hat{v}(x,t)=v_2+{\rho }_2-\hat{\rho }(x,t)$ . The solution given by (17) is valid for t > tc where $t_c=\frac{L}{v_1-v_2}$ is the critical time and compared with the approximate solution in Fig 4

Fig. 3. Solutions evolution for density and velocity (left) with the corresponding solutions at t = 8 s for example Ej. 2.

Source: autor.

In Fig 4 it can be observed that the numerical solution agrees very well with the analytical solution which corresponds to an initial shock and a generalized rarefaction wave propagating backward until they meet.

Fig. 4. Comparison between exact and approximate solutions for double Riemann problem at t = 5 s.

Source: author.

Next, Table 2 shows the relative errors E1 (M) and E2 (M) in L1-norm for different values of M (number of points of spatial mesh). These errors are calculated according to the expressions given in (18):

$E_1\left(M\right)=\frac{\sum j_1^M\left|{\rho }_j\right(t)-\rho (x_j,t\left)\right|}{\sum j_1^M\left|\rho \right(x_j,t\left)\right|},E_2\left(M\right)=\frac{\sum j_1^M\left|v_j\right(t)-v(x_j,t\left)\right|}{\sum j_1^M\left|v\right(x_j,t\left)\right|}$ (18)

where ρj (t) is the numerical approximation of the exact solution ρ(xj,t) at the discrete point (xj,t) given by (17). Note that although the solution has discontinuities, the error decreases consistently as the mesh is refined.

Table 2. Relative errors for density and velocity at final time t = 5 s.

M	E1 (M)	E2 (M)
100	6.128×10-3	7.340×10-3
200	4.978×10-3	5.119×10-3
400	1.859×10-3	1.866×10-3
800	8.751×10-4	8.628×10-4
1600	4.669×10-4	4.287×10-4

Source: author.

Fig 5 displays the errors versus mesh width are shown in logarithmic scale in order to observe the behavior of the convergence rate (estimation of the slope of the straight lines) which is approximately equal to 1.

Fig. 5. Approximate -errors for density and velocity in logarithmic scale.

Source: author.

CONCLUSIONS

In this paper numerical approximations for the solutions of a second-order traffic flow model were computed using the well-known fifth-order WENO scheme for spatial discretization and a sign-preserving Runge-Kutta exponential method for the temporal discretization obtaining a explicit scheme that does not require taking restrictive time steps for stability. In numerical experiments it was observed that different kinds of waves are correctly captured by the method and the errors exhibit a consistent behavior.

REFERENCES

[1] F. Kessels. (2018, Aug. 21). Traffic Flow Modelling: Introduction to Traffic Flow Theory Through a Genealogy of Models. (First ed.) [Online]. Available: https://doi.org/10.1007/978-3-319-78695-7_1

[2] M.J. Lighthill, G.B. Whitham, “On kinematic waves II: a theory of traffic flow on long crowded roads”, Proc. R Soc Lond A Math Phys Sci., vol. 229, no. 1178, pp. 317–345, May. 1955. https://doi.org/10.1098/rspa.1955.0089

[3] P.I. Richards, (1956) “Shock waves on the highway”, Oper. Res., vol. 4, no. 1, pp. 42–51, Feb. 1956. https://www.jstor.org/stable/167515

[4] A. Aw and M. Rascle, “Resurrection of second order models of traffic flow”, SIAM J. Appl. Math., vol. 60, no. 3, pp. 916–938, Feb. 2000. https://doi.org/10.1137/S0036139997332099

[5] A. Aw, A. Klar, T. Materne and M. Rascle, “Derivation of continuum traffic flow models from microscopic follow-the-leader models”, SIAM J. Appl. Math., vol. 63, no. 1, pp. 259–278, Aug. 2002. https://doi.org/10.1137/S0036139900380955

[6] L. Pareschi, G. Russo, “Implicit-Explicit Runge-Kutta Schemes and Applications to Hyperbolic Systems with Relaxation”, J. Sci. Comput., vol. 25, Nov. 2005. https://doi.org/10.1007/s10915-004-4636-4

[7] A. Chertock, S. Cui, A. Kurganov and T. Wu, “Steady state and sign preserving semi-implicit Runge-Kutta methods for ODES with stiff damping term”, SIAM J. Numer. Anal., vol. 53, no. 4, pp. 2008-2029, Aug. 2015. https://doi.org/10.1137/151005798

[8] J. Du, Y. Yang, “Third- order conservative sign-preserving and steady-state preserving time integrations and applications in stiff multispecies and multireaction detonations”, J. Comput. Phys.. vol. 395, pp. 489-510, Oct. 2019. https://doi.org/10.1016/j.jcp.2019.06.040

[9] R. Yang, Y. Yang and Y. Xing, “High order sign-preserving and well-balanced exponential Runge-Kutta discontinuous Galerkin methods for the shallow water equations with friction”, J. Comput. Phys., vol. 444, pp. 1-23, Nov. 2021. https://doi.org/10.1016/j.jcp.2021.110543

[10] C. W. Shu, “Essentially non-oscillatory and weighted essentially non-oscillatory schemes”, Acta Numerica, vol. 29, pp 701-762, May. 2020. https://doi.org/10.1017/S0962492920000057

[11] G. Jiang, C.W. Shu. “Efficient implementation of weighted ENO schemes”. J. Comput. Phys.. vol. 126, pp. 202-228, 1996. https://doi.org/10.1006/jcph.1996.0130

[12] A. Henrick, T. Aslam, J. Powers. “Mapped weighted essentially non-oscillatory schemes: Achieving optimal order near critical points”. J. Comput. Phys. Vol. 207, pp. 542-567, 2005. https://doi.org/10.1016/j.jcp.2005.01.023

[13] A. Jannelli, N. Manganaro, A. Rizzo. “Riemann problems for nonhomogeneous Aw-Rascle model ” Commun. Nonlinear Sci. Numer. Simul., vol. 118, pp. 1-13, 2023. https://doi.org/10.1016/j.cnsns.2022.107010

The Author received the M.Sc. degree in Mathematics from Universidad Nacional de Colombia (Sede Medellín) and Universidad del Norte and the Ph.D. degree in Mathematical Engineering from Universidad de Concepción, Chile, in 2004 and 2010, respectively. Since 2011, I have been a full profesor at the Universidad del Norte in Barranquilla – Colombia.