# Adaptive Numerical Method for Approximation of Traffic Flow Equations

Authors: Neda Najafzadeh; Saeid Abbasbandy; Elyas Shivanian
Journal: Journal of Advanced Transportation
 (2022)

---

## Abstract

For a long time now, traffic equations have been considered, and different modeling has been done for it. In this article, we work on the macroscopic model, especially the most famous light model. Because these models are among the stiff and shocking problems, theoretical methods do not give good answers to these problems. This paper describes a meshless method to solve the traffic flow equation as a stiff equation. In the proposed method, we also use the exponential time differencing (ETD) method and the exponential time differencing fourth-order Runge–Kutta (ETDRK4). The purpose of this new method is to use methods of the moving least squares (MLS) method and a modified exponential time differencing fourth-order Runge–Kutta scheme. To solve these equations, we use the meshless method MLS to approximate the spatial derivatives and then use method ETDRK4 to obtain approximate solutions. In order to improve the possible instabilities of method ETDRK4, approaches have been stated. The MLS method provided good results for these equations due to its high flexibility and high accuracy and has a moving window and obtains the solution at the shock point without any false oscillations. The technique is described in detail, and a number of computational examples are presented.


---

## Body

## 1. Introduction

Car traffic modeling has been around since time immemorial, around the early twentieth century. Many methods have been modeled according to different cases, from microscopic methods that have been modeled individually for each vehicle to dynamic kinetic and macroscopic methods that deal with the relationship between vehicles and even driver behavior, all of which have been modeled in different ways ([1–5]).There are many articles on traffic flow modeling. Existing traffic flow simulations are usually based on mesoscopic modeling; that is, models that are microscopic simulations is mixed with macroscopy, meaning that, in some cases, it is based on aggregation, or macroscopic models that are combined with microscopic modes, and in some cases, it isbased on the individual mode of each vehicle [6–8].There are many traffic models in these categories, such as automated cellular models, car chase models, gas kinetic models, fluid flow models, and hydrodynamics [9–11]. Microscopic models have many applications, such as traffic forecasting, accident detection, and traffic control, which are typically described using a system of nonlinear partial differential equations. Macroscopic traffic models look at traffic flow as a cumulative concept and are commonly used to analyze service level and also supply and demand, during regional planning or surface transportation networks [12]. In 1955, Lighthill and Whitham and, in 1956, Richards showed that equations expressing water flow could explain the flow of car traffic. Their idea was to assume machines to be small particles and their density to be principal quantity. They established one of the most famous models of traffic flow called model LWR. The LWR model is the first model of traffic flow in the form of equations of nonlinear hyperbolic partial derivatives [7, 13].In any case, it makes sense to assume that maintaining the number of vehicles leads to a conservation laws. Because they show severe discontinuities in traffic jams, there is an adaptive relationship between traffic jams and shock waves [6, 14–16]. Fluid dynamic models are often considered the best simulations for traffic flow because they seem to be better able to detect some phenomena, such as the formation and propagation of shocks on roads, and the answers to problems can create strong discontinuities, even with perfectly smooth initial conditions (see [17]). We investigate traffic simulations from the term of macroscopic models that use fluid dynamics and consist of conservation laws. In Section 2, we will explain these equations in more detail.In this article, an adaptive numerical method is proposed based on the method the meshless method and ETD scheme. In recent years, meshless methods including MLS scheme have been used for diverse differential equations.One of the numerical methods examined in this paper and used to discretize time is the “exponential temporal difference” method, which involves the correct integration of dominant equations and then the approximate integral of nonlinear expressions [18, 19]. In principle, this method was in the category of first-order explicit methods, which over time has been extended to explicit and explicit designs. We discuss the sustainability of these methods, which are further explored in Ref. [20]. To improve ETD methods, we consider new and more accurate ETD methods. The method combines ETD scheme with Runge–Kutta [19]. As Cox and Matthews (the originators of this method) recognized, they have a big problem with eigenvalues close to zero, especially when the matrix for the linear part (L) is not diagonal. If these problems are not solved, ETD schemes for PDE equations whose matrix is not diagonal will have many errors and then fail. In this work, we use modified ETD schemes to avoid these problems [21]. Low-order ETD methods are more suitable for application in computational electrodynamics [22]. They are often taken independently [19, 20, 23]. In his article, Iserles stated that in 1928, Filon proposed ideas related to this method in the field of ODE [24, 25]. The best solution to the problems of this method is to use the fourth-order Runge–Kutta formula (ETDRK4) with exponential time difference [19]. Cox and Matthews argued that ETD schemes perform better than implicit-explicit (IMEX) schemes, where linear expressions predominate and perform better than integration-factor (IF) schemes, where nonlinear expressions predominate. We will explain more about this method in Section 3.In recent decades, methods for better solving equations called meshless methods have been proposed [26]. These methods are used to discretize the domain without using a predefined mesh by creating a system of algebraic equations for the entire problem domain. How this method works is that a set of nodes that gets diffused across the domain of the problem is used to indicate its domain, and a set of nodes is used at the domain boundaries to represent its boundaries. The set of these diffuse nodes is called field nodes [27]. Mesh-free methods can be divided into three groups:(1)
Mesh-free methods based on weak forms, such as the element-free Galerkin method (EFG) [28](2)
Mesh-free methods based on strong forms, such as collocation method based on radial basis functions (RBFs) [29, 30](3)
Mesh-free methods based on a combination of weak forms and collocation method [31, 32]In the written works, several methods without weak mesh are stated which are as follows:(1)
Dispersed elements method (DEM) [33](2)
Smooth particle hydrodynamics (SPH) [34](3)
Reproducing kernel particle method (RKPM) [35](4)
Boundary node method (BNM) [36](5)
Partition unity finite element method (PUFEM) [37](6)
Finite sphere method (FSM) [38](7)
Boundary point interpolation method (BPIM) [39](8)
Border radius point interpolation method (BRPIM) [40](9)
Meshless local Petrov-Galerkin(MLPG) [41](10)
Meshless local radius point interpolation (MLRPI) [42–47]One of the meshless methods is moving least squares (MLS) method, in which we use a local interpolation or approximation to express an experimental function with unknown values at some node points [48]. In this paper, the moving least squares (MLS) approximation is used. In Section 4, we discuss this method in detail. In Section 5, we explain about our adaptive method. In Section 6, a number of examples are solved with these methods.


## 2. The Traffic Flow Equation

Traffic modeling and simulation is the subject of research that has attracted much attention today. Travel information technology seeks to alleviate the problems of congestion and travel time that travelers face in most cities around the world. Attempts are being made to answer questions,such as where should the traffic lights or stop signs be installed, what should be the traffic light cycle, and where is the construction site of entrance, exit, and overpass?The aim of this research is to solve the traffic equation in a way that has the least fluctuations and obtain the solution of the shock point with the least error. When solving traffic equations is done by researchers, more general goals such as increasing car traffic and reducing traffic congestion, accidents, and air pollution can be achieved. This is possible through better scheduling of traffic lights and informing drivers in choosing the route based on real and correct information of traffic situations. These issues have led to the issue of modeling and predicting traffic parameters, as one of the main components has attracted considerable scientific interest. For this reason, mathematical modeling and solving traffic flow models as one of the effective methods of traffic management is very important.Traffic models can be found in approximately three models: (1) microscopic models, (2) mesoscopic models, and (3) macroscopic models. There are many models by researchers who study traffic from other perspectives, see Refs. [1–5]. In many cases, attention is paid to a section of road or small sections of the urban network. One of the most widely used models of traffic flow is the macroscopic model that expresses traffic flow as one or a set of nonlinear hyperbolic partial equations. Since there are sudden changes or shocks in the answer to most of these equations, choosing the appropriate numerical method to solve this category of equations is critical.The macro theory looks at traffic flow as a cumulative concept. This method is the most appropriate method to study the steady state of the flow phenomenon based on physical analogies such as heat flow and fluid flow, and therefore, it best explains the device’s potential performance. In these models, traffic flow simulations are performed on a part of the road. The interactions of individual road users are not considered.Macrosimulation models are commonly used to analyze service, supply, and demand levels during regional planning or large-scale transportation networks. Micro models are more detailed than macro models and can therefore be used to assess the effects of the proposed level of improvement on road facilities with a higher degree of accuracy. However, due to the nature and breadth of information that the shredding models simulate, they run slowly compared to the macro models, so they require lengthy calculations.In this paper, we examine models from a macroscopic perspective that are inspired by fluid dynamics and composed of conservation laws. The traffic flow equation is(1)∂ρ∂t+∂f∂x=0.Macroscopic models present the parameters of traffic densityρ, flow rate f, and average speed v, and the basic relationship is f=ρvρx,t, where t is time and x is space. The number of vehicles per unit length is called traffic flow density ρ=n/Δx. Traffic flow rate indicates the number of vehicles passing a certain point per unit of time f=n/Δt.In the model LWR, speed is expressed asv=vρ, s.t. 0≤ρ≤1. Here, it is assumed that the function f has the following conditions:(1)
The functionf is continuous until the second time(2)
The functionf is strictly concave(3)
Letf0=fρmax=0When the number of vehicles gradually increases, the density and therefore the flow also increases. If the number of cars is constantly increased, we will reach a situation where other cars cannot move. In this case, the density is equal to the density of traffic jams and the flow is zero, andρmax is the density of traffic jams.One of the most important discussions of traffic modeling is how to define the velocity function as a function of density. There are many models for defining velocity in terms of density [49].


## 3. The MLS Approximation Scheme

Meshless methods work without using a predetermined meshing to discretize the problem domain. Among the many meshless methods, we chose method 1 because there were good reasons for this choice; that is, the development of a smooth weight function that is a guarantee for the continuity of variable fields, very high flexibility of this method, and high accuracy of this method.MLS combines the moving window with weight functions that have compact support. This method is very accurate even with low-order basis functions and also has very good stability because it uses a local approximation scheme. Calculations are performed in subdomains of the global domain. A very attractive feature of MLS is its flexibility. With the parameters of the weight function, it is easy to provide the appropriate tuning for the MLS and thus control and adjust the data fit.In the traffic equations, shocks inevitably appear and the fitting method is preferred to solve these equations. The least squares (LS) method is the most used one in data fitting.Most LS-based methods are global approximation schemes that are not suitable for large amounts of data, irregular or sparse distribution cases. Therefore, the MLS method, which is a local approximation and has high computational accuracy and stability, was proposed. The moving window introduced in MLS shows a good performance in solving these problems. The compact support weight function indicates that only partial measurement data near the unknown measurement point are involved in the calculation. Then, the moving window in MLS acts as a smooth section. Due to the fact that rigid segmentation causes problems in the way of part selection and fitting discontinuity, segment selection is avoided and thus the continuity and smoothness of the connection is guaranteed.In this section, we explain one of the meshless methods called MLS for the approximation of the functionvx in the domain Λ. Assume that Λ is expressed by computational geometry techniques. In this supposed domain, we assume the nodes xi,i=1,…,N and denote the approximate value of the related with the MLS method at the point of the assumed node i with vi. As explained in section 1, meshless methods use local interpolation or approximation to express an experimental function that at several node points is the value of unknown variable.Here, among the meshless methods, we used the moving least square (MLS) method because it seems to work better for shock problems. We consider a subdomainΛs. In the neighborhood, the point x, which exists within the principal domain of the problem Λ, is represented as the support domain of the MLS approximation for the experimental function at the point x.To approximate the functionv in Λs, on random points of nodes xi,i=1,…,N. We illustrate the approximation of the v function by the method of moving least squares with vhx upon the domain Λs, and we can define ∀x∈Λsas (2)vhx=DTxqx,∀x∈Λs,where DTx is a vector of order m as DTx=d1x,d2x,…,dmx. DTx is a complete monomial basis, and qx is a vector since qx=q1x,q2x,…,qmx. Let qjx,j=1,2,…,m be a function of the coordinates of space x, and djx be monomial in the coordinates of space, and m be the number of basic polynomial functions. To make djx, the Pascal triangle is used and a complete base is certainly preferred. The basic functions are ascertained by(3)DTx=1,m=1,1,x,m=2,1,x,x2,m=3.Here, if we want to shift the origin to a fixed pointxe=X1e,X2e,…,XneT, x in Dx must be replaced by x−xe. Point xe=X1e,X2e,…,XneT is on ℜx, where ℜx represents the influence domain of x. Then, a linear basis is given by(4)Dx=1,x1−x1e,x2−x2e,…,xn−xneT,x∈ℝn,for m=1 and a quadratic basis by(5)D=1,x1−x1e,x1−x1e2T,x∈ℝ,1,x1−x1e,x2−x2e,x1−x1e2,x2−x2e2,x1−x1ex2−x2eT,x∈ℝ2,for m=2. The weight function used in the MLS method is as follows:(6)wix=φx−xih,i=1,…,N,where the function φ is non-negative, γ th times continuously differentiable, and its derivatives are bounded for a higher order than γ and compactly supported. Many weight functions used in meshless methods have these conditions, such as Gaussian, exponential, cubic, and quadratic splines. In this paper, the Gaussian weight function is used, which is defined as follows:(7)wix=exp−di/ci2−exp−rs/ci21−exp−rs/ci2,0≤di≤rs,0,di≥rs,where di=x−xi,ci is a constant that controls the shape of the weight function wi and rs is the size of the support domain. The ci parameter is a constant that controls the shape of the weight function, which can be dependent on the points for more flexibility. This shape parameter, in general, may vary from one subdomain to another. In addition, this shape parameter can be chosen differently in different directions. As an example, a problem with nine nodes that is uniformly distributed is presented in Figure 1(a). The most appropriate weight function for this case is shown in Figure 1(b). In Figure 2(a), the uniform distribution of nodes depends on the coordinate direction, where the grid spacing is different in the x and y directions. To obtain more accurate spatial derivatives, the weight function for this case should be different in each direction in Figure 2(b). In other words, the ci parameter must be larger in the x direction. To read more about weight function and shape parameter, refer to [50].Figure 1
(a) Nodes. (b) Weight function.
(a)(b)Figure 2
(a) Nodes. (b) Weight function.
(a)(b)The support size of the weight functionwi related to node irs is very important. On the one hand, it should be chosen so that it is large enough to cover the required number of nodes in the definition domain of each sample point, and on the other hand, a very small rs It can also cause a large error in calculating the entries of matrix when using Gauss numeric quadrants. Care must also be taken that rs should be small enough, that is, not too small to preserve the local properties of the MLS approximation.As it was saidvhx=DTxqx that we talked about DTx, now we are talking about the coefficients of qx in detail.To define a coefficient vectorqx, qx, we use the L2-norm in such a way that these coefficients are determined by minimizing the weighted discrete L2-norm:(8)Jx=∑i=1nwixDTxiqx−v^i2=D.qx−v^T.W.D.qx−v^.Here,wix is defined as the weight function for node i. Note that for all points in the support domain of wix, wix, the value of the weight function at each point, i.e., wix, wix is positive. xi represents the value of x in node i, and n is the number of nodes in Λs, for each weight function is wix>0. The matrices D and W are defined as(9)D=dTx1dTx2…dTxn,W=diagw1x,…,wnx,v^T=v^1,…,v^n.Note thatv^i,i=1,2,…,n are not node values of the unknown trial function vhx but are artificial node values. There is a linear relationship between qx and v^. This relationship is derived from the stationarity of J in (2) with respect to qx, i.e.,(10)Qxqx=Sxv^.Define the matricesQx and Sx as(11)Qx=DT.W.D,Sx=DT.W.Obtainingqx from the equation (10) and substituting it in (2) yields(12)vhx=μTx.v^=∑i=1nμixvi^,x∈Λs,(13)μTx=DTxQ−1xSx,or(14)μix=∑j=1mdjxQ−1xSxji,where μix is called a shape function of the MLS approximation dependent on the xi node. From (13) and (11), it can sometimes be concluded that μix=0 when wix=0. In practical applications, usually wix is selected such that over the support of xi node is nonzero. In fact, μix=0, for x that is not in the support area of the node xi, is to maintain the local character of the MLS approximation.Note that the sustainability of the linear system is effectively dependent on conditioning of the coefficients matrix. The conditioning of the coefficients matrix can be measured using the condition number.Theorem 1.
If we consider theQx matrix that is produced by the basis Dx (i.e., Qx=DTWxD), then there is a bounded and countable number Cdx,n,m, which is independent of h and such that for the determinate of Qx, we get(15)detQx=Cdx,n,mh2∑i=1mi^,∀x∈Λ.
In addition, there is a constanth0>0, for h≤h0, such that a bounded and computable number Ccx,n,m is independent of h that the norm 2 condition number of Qx is(16)condQx=Ccx,n,mh−2m^,∀x∈Λ.Proof.
A proof is described in [51].
Equation (15) states that(17)detQx⟶0ash⟶0.
Note that from the above discussion we conclude that whenever a sufficiently small nodal distanceh is selected, the instantaneous matrix Qx obtained in MLS approximation is approximately a single matrix. This causes the errors to increase; that is, when the instantaneous matrix Qx has ill conditions or is singular, the answers obtained for the vector of coefficients qx have a large error, and therefore, the function of the shape μix will have a large rounding error.
Referring to (16), we find that the condition number Qx, Qx has an inverse relationship with factor h2m^.That is, as h2m^ decreases, the condition number increases and as a result the rounding error increases because it is very clear that when the matrix conditioning number increases, the degree of ill conditioning of the matrix increases to the point that for singular matrices, the situation worsens and the condition number Infinitely desires. Therefore, we conclude that the instability of the MLS method depends on h and m^ so that by reducing the value of h and increasing the value of m^ too much, the instability of the MLS method increases. So we conclude that the MLS method achieves well and low-error solutions to solve problems when matrix Q in (10) is nonsingular. This happens when the rank D is equal to m.
In order not to have the problem of instability and MLS method being a well-defined method, a necessary condition is that there must be the least at leastm nonzero weight functions (i.e., n>m) for each point x∈Λ and the points determined in Λs do not follow from a regular and special pattern.


## 4. A Modified ETD Scheme

The traffic flow equation combines nonlinear terms with linear terms. To solve these problems and obtain low-error numerical answers to these problems, high-order approximations for place and time must be used. Because we have a nonlinear and stiff combination in these problems, it has not been used much until the second order. The general form of the stiff equations is as follows:(18)vt=ℒv+Nv,t.In this equation,ℒ represents the linear operator and N represents the nonlinear operator. First, the spatial part of the PDE equation is discretized by the method of the previous section, and thus we obtain a system of ODEs as(19)vt=Lv+Nv,t,where L,N are the matrices obtained after using method MLS and approximating spatial derivatives on linear and nonlinear operators, ℒ and N. For this special traffic equation, both matrices are written in Section 5. We begin by multiplying equation (19) through by e−Lt, where the term e−Lt is the integrating factor, then we take the integral from the equation on a single time step of length h from t=tntot=tn+1=tn+h to give(20)vn+1=eLhvn+eLh∫0he−LζNvtn+ζ,tn+ζdζ.In this equation, different orders of ETD schemes can be obtained from how the integrals are approximated. A sequence of recursive formulas is presented that obtains higher order approximations than the multistage type. Andthe generating formula is(21)vn+1=eLhvn+h∑m=0s−1ηm∑k=0m−1kmkNn−k,wherein s represents the order of the method. The coefficients ηm can be obtained from the following recursive relation:(22)Lhη0=eLh−I,Lhηm+1+I=ηm+12ηm−1+13ηm−2⋯+1m+1η0.Whenever we use Runge–Kutta methods in the set of ETD methods and derive its time step with Runge–Kutta, then the method is called ETDRK. In this paper, we use only the fourth-order method and derive it with ETD methods. In this case, the method is called ETDRK4. This is not entirely clear and needs to be explained in the form of a symbolic manipulation system. Here are the formulas for the ETDRK4 method:(23)an=eLh/2vn+eLh/2−ILNvn,tn,bn=eLh/2vn+eLh/2−ILNan,tn+h/2,cn=eLh/2an+eLh/2−IL2Nbn,tn+h/2−Nvn,tn,ϖ=−4−Lh+eLh4−3Lh+Lh2,ϱ=2+Lh+eLh−2+Lh,ς=−4−3Lh−Lh2+eLh4−Lh,vn+1=eLhvn+ϖNvn,tn+2ϱNan,tn+h/2+Nbn,tn+h/2+ςNcn,tn+hh2L3.The coefficients inside the ETDRK4 formula can be rewritten as follows:(24)α=ϖh2L3,β=ϱh2L3,γ=ςh2L3.As mentioned before, we tried to solve the traffic equations with the methods described, and an example of these traffic equations is solved in the next section.In the written codes, for all these formulas, division by matrixL means the inverse of the matrix L is multiplied in the desired expression and to calculate eL, commands exp and expm of MATLAB software are used, if L is diagonal, command exp is used, and if it is a nondiagonal matrix, command expm, expm is used.


## 5. Adaptive Method

In this article, we have tried to solve equation(25)∂ρ∂t+∂fρx,t∂x=0,by a adaptive method. In this section, we explain how to solve this equation according to methods MLS and ETDRK4. To solve equation, we first use the MLS method to approximate spatial derivatives and then use the ETDRK4 method to approximate time derivatives. We divide the main domain of the problem Λ into N parts so that the distance between the nodes is constant, and for each point, we consider a neighborhood Λs that is inside the main domain. To approximate the function ρ, by MLS method, in Λs, on random points of nodes xi,i=1,…,n, from (12) and (14), we get(26)ρhx,t=∑i=1nμixρ^it.By replacing the equation, we have(27)∂ρ∂t+∂f∑i=1nμixρ^it∂x=0.Note thatρ^it=ρ^xi,t in (27) are node values that are fictitious. On the other hand, for the derivative μix, it can be explained as follows.LetCqΛ be the space of functions that are continuously qth differentiable on Λ. If wix∈CqΛ and djx∈CsΛ,i=1,2,…,n,j=1,2,…,m, then μix∈CrΛ with r=minq,s. The k derivatives of μix becomes(28)μi,k=∑j=1mdj,kQ−1Sji+djQ−1S,k+Q,k−1Sji,where we show the derivative inverse of the matrix Q, which depends on xk, with the symbols Q,k−1=Q−1,k, which can be obtained by the formula Q,k−1=−Q−1Q,kQ−1. As mentioned in Section 2, the density function is f=ρVρ. In the first example Vρ=vmax1−ρ, where vmax=1, then we can write (29)∂ρ∂t=−∂ρ∂x+2ρ∂ρ∂x.Now, considering the spatial discretization and using the MLS method, we have(30)∂ρx,t∂t=−∑j=1nμj′xρ^jt+2ρx,t∑j=1nμj′xρ^jt.For xi∈Λ (Λ is the global domain of problem), we have(31)ρi′t=−∑j=1nμj′xiρ^jt+2ρit∑j=1nμj′xiρ^jt.Note that the derivative ofρ is relative to time, and the derivative of ϕ is defined in 20. When we disassemble the PDE space partition, we get a system of ODEs as(32)ρ1′tρ2′t⋮ρn′t=−μ1′x1μ2′x1…μn′x1μ1′x2μ2′x2…μn′x2⋮⋮⋮⋮μ1′xnμ2′xn…μn′xnρ^1tρ^2t⋮ρ^nt+2ρ1tρ2t⋮ρnt°μ1′x1μ2′x1…μn′x1μ1′x2μ2′x2…μn′x2⋮⋮⋮⋮μ1′xnμ2′xn…μn′xnρ^1tρ^2t⋮ρ^nt,where ∘ is the Hadamard product.In this step, to approximate time derivatives by method ETDRK4, we have to divide the flux function into two parts, linear and nonlinear:(33)L=−μ1′x1μ2′x1…μn′x1μ1′x2μ2′x2…μn′x2⋮⋮⋮⋮μ1′xnμ2′xn…μn′xn,and(34)Nρ,t=2ρ1tρ2t⋮ρnt°μ1′x1μ2′x1…μn′x1μ1′x2μ2′x2…μn′x2⋮⋮⋮⋮μ1′xnμ2′xn…μn′xnρ^1tρ^2t⋮ρ^nt.For time discretization, the method ETD, which is described in section, is used. Assuming that we divide the main domain inton parts with equal distances, it is assumed that ρxi,kt,i=1,2,…,n is specified, and our goal is to calculate ρxi,k+1t,i=1,2,…,n. Multiply the equation by e−Lt and then integrate the equation into a single time step of length h, so we get(35)ρik+1=eLhρik+eLh∫0he−Lζ2ρitk+ζ∑j=1nμj′xiρ^jtk+ζdζ,where ρik=ρxi,tk and then we solve it using the formulas of method ETDRK4 that is expressed in (15), but we have problems with the ETDRK4 method and that it becomes unstable for some points. One way to solve this problem is to use a Taylor expansion point with a cut-off point for small eigenvalues so that if the matrix is a linear diagonal operator, the Taylor expansion is used for the diagonal points below the cut, but one of the serious problems with this method is that it possible to extend this method to nondiagonal problems? That is why we went the other way. Another way that can be used is to use the contour integral of the idea of complex analysis. In fact, we have a function to evaluate that has a single point and is analytical elsewhere. Near this singular point, the function has a numerical error. The suggested solution is to use a contour integral on a complex plane that includes z and is separated from z=0 by(36)fz=2πi−1∫Γfss−z−1ds.Now instead ofz-scalar we use the matrix L and apply the same method again, only the expression 1/s−z becomes the matrix sI−L−1 as(37)fL=2πi−1∫ΓfssI−L−1ds,where the contour must be such that it contains the eigenvalues of L. The result of contour integrals of functions in the complex plane is obtained using the trapezoidal rule, which converges exponentially.The contour area can be considered as a circle, and we usually consider 32 points with equal distance in this circle. There are many different choices for the contour region that work well, but the important thing is that the eigenvalues are really restricted byΓ. When L is real, we use only the upper half of the circle for the contour area whose center is the real axis, then we get the real part of the result. The contour integral can be approximated by the mean fs on Γ so that we obtain the mean at points of equal distance along Γ, or again consider only the upper half of Γ, and then consider only that real part.


## 6. Numerical Example

We have described the method implemented in this article in detail in the previous sections. In this section, we will solve two examples of traffic equations and compare their numerical answer form with other methods. As explained, because these equations are in the category of stiff equations, it is very important to display the numerical answer obtained with the least vibrations at the moment of shock and after.


### 6.1. Example 1

Consider the following traffic flow equation:(38)∂ρx,t∂t+∂Fρx,t∂x=0.The density function in this model is expressed as follows:(39)Fρx,t=ρVρ,where V is(40)Vρ=vmax1−ρ,0≤ρ≤1.Here, we assume thatvmax=1, and the initial conditions are as follows:(41)ρx,0=0.25+0.7e−βx2,where β=0.01.The density bulge (or density shock) changes over time. After a while, the shock wave simulation shows that the density of the traffic flow increases and the speed decreases very quickly (drivers who were speeding in light traffic suddenly had to brake). As vehicles exit the shock zone, they gradually accelerate, resulting in a gradual decrease in density. Figure1(a) shows the numerical answer to this problem in the method presented in this paper and rejects the shape of the shock part well. The figure is drawn to t=30.Because these problems do not have an exact answer and they are very difficult to solve, we tried to solve this problem with different methods, and the best answer is obtained from the method mentioned in the article. It is important that the movement of the shock point in the figures after the occurrence of traffic is completely clear, but in some methods, since the occurrence of the shock, the error is very high and does not show the movement of the shock. However, in the method mentioned in the article, as shown in Figure 3(a), it shows the shock point well.Figure 3
A comparison of the solution of equation of Example 1 using the present method: (a) with solution by Mathematica Software and (b) solution by the method of reproduction kernel.
(a)(b)(c)Compared to Figure 3(b), which was obtained by Mathematica Software and did not pass from time 8 due to shock, but in our calculation method, it was obtained up to time 30 and the figure is acceptable. In Figure 3(c), an attempt has been made to solve the equation by the method of reproduction kernel, where unfortunately the error has increased so much after the shock time, which is quite clear in the figure. Also, in Figure4, we compare the moment t=25, which is the time after the shock, with the figure shown in Ref. [52], which shows the routes of individual vehicles and the density observed by each driver. A comparison of the obtained solution using the present method with MLS and without MLS is shown in Figure 5. Table 1 shows a comparison of the obtained solution using the article method with other numerical methods.Figure 4
A comparison of the solution of equation of example 1 using the present method with the solution in Ref. [52].
(a)(b)Figure 5
A comparison of the solution of equation of example 1 using the present method with the solution without MLS.Table 1
A comparison of the solution of equation of example 1 using the present method with other numerical methods.
txWithout MLSMathematicaNew method5−200.26240.26400.2631−100.61360.56000.561300.88720.90890.9184100.46530.46950.4726200.26180.26240.262610−200.2628……0.2633−100.7288……0.653700.8877……0.8535100.4334……0.4476200.2616……0.262515−200.2632……0.2636−100.8444……0.815100.7369……0.7877100.4112……0.4285200.2615……0.262320−200.2599……0.2639−100.9019……0.934700.7410……0.7313100.3976……0.4134200.2621……0.262225−200.2538……0.2646−100.8550……0.938800.5769……0.6844100.3849……0.4016200.2606……0.2620


### 6.2. Example 2

In the second example, we consider another model of the traffic problem in which the density function is the same as in the previous example, but the initial data are piecewise function:(42)ρx,0=1,4≤x≤6,0.5,otherwise.This can be interpreted as cars approaching a congested traffic or encountering a line of cars behind a red light waiting to turn green. This example describes the traffic situation at the traffic light installation site. Here, the traffic light is in positionx=6, and at moment t=0, it has changed from red to green. Figure 6 shows the numerical solutions up to time t=30 computed in the new proposed scheme.Figure 6
Example 2: forN=120 and t=30.


## 6.1. Example 1

Consider the following traffic flow equation:(38)∂ρx,t∂t+∂Fρx,t∂x=0.The density function in this model is expressed as follows:(39)Fρx,t=ρVρ,where V is(40)Vρ=vmax1−ρ,0≤ρ≤1.Here, we assume thatvmax=1, and the initial conditions are as follows:(41)ρx,0=0.25+0.7e−βx2,where β=0.01.The density bulge (or density shock) changes over time. After a while, the shock wave simulation shows that the density of the traffic flow increases and the speed decreases very quickly (drivers who were speeding in light traffic suddenly had to brake). As vehicles exit the shock zone, they gradually accelerate, resulting in a gradual decrease in density. Figure1(a) shows the numerical answer to this problem in the method presented in this paper and rejects the shape of the shock part well. The figure is drawn to t=30.Because these problems do not have an exact answer and they are very difficult to solve, we tried to solve this problem with different methods, and the best answer is obtained from the method mentioned in the article. It is important that the movement of the shock point in the figures after the occurrence of traffic is completely clear, but in some methods, since the occurrence of the shock, the error is very high and does not show the movement of the shock. However, in the method mentioned in the article, as shown in Figure 3(a), it shows the shock point well.Figure 3
A comparison of the solution of equation of Example 1 using the present method: (a) with solution by Mathematica Software and (b) solution by the method of reproduction kernel.
(a)(b)(c)Compared to Figure 3(b), which was obtained by Mathematica Software and did not pass from time 8 due to shock, but in our calculation method, it was obtained up to time 30 and the figure is acceptable. In Figure 3(c), an attempt has been made to solve the equation by the method of reproduction kernel, where unfortunately the error has increased so much after the shock time, which is quite clear in the figure. Also, in Figure4, we compare the moment t=25, which is the time after the shock, with the figure shown in Ref. [52], which shows the routes of individual vehicles and the density observed by each driver. A comparison of the obtained solution using the present method with MLS and without MLS is shown in Figure 5. Table 1 shows a comparison of the obtained solution using the article method with other numerical methods.Figure 4
A comparison of the solution of equation of example 1 using the present method with the solution in Ref. [52].
(a)(b)Figure 5
A comparison of the solution of equation of example 1 using the present method with the solution without MLS.Table 1
A comparison of the solution of equation of example 1 using the present method with other numerical methods.
txWithout MLSMathematicaNew method5−200.26240.26400.2631−100.61360.56000.561300.88720.90890.9184100.46530.46950.4726200.26180.26240.262610−200.2628……0.2633−100.7288……0.653700.8877……0.8535100.4334……0.4476200.2616……0.262515−200.2632……0.2636−100.8444……0.815100.7369……0.7877100.4112……0.4285200.2615……0.262320−200.2599……0.2639−100.9019……0.934700.7410……0.7313100.3976……0.4134200.2621……0.262225−200.2538……0.2646−100.8550……0.938800.5769……0.6844100.3849……0.4016200.2606……0.2620


## 6.2. Example 2

In the second example, we consider another model of the traffic problem in which the density function is the same as in the previous example, but the initial data are piecewise function:(42)ρx,0=1,4≤x≤6,0.5,otherwise.This can be interpreted as cars approaching a congested traffic or encountering a line of cars behind a red light waiting to turn green. This example describes the traffic situation at the traffic light installation site. Here, the traffic light is in positionx=6, and at moment t=0, it has changed from red to green. Figure 6 shows the numerical solutions up to time t=30 computed in the new proposed scheme.Figure 6
Example 2: forN=120 and t=30.


## 7. Conclusion

Traffic issues are very important, and on the other hand, they are very difficult to solve due to problems such as having a shock, and it is very important to choose a good numerical method for traffic flow models.In the present article, a new method has been proposed to identify shock points well and to find a better answer without vibrations in these points. In the method described in the article, first we solved such problems from meshless methods and then with the ETDRK4 method. We chose the MLS method from among the meshless methods. In this method, we used quadratic basic functions. Then, we used the ETDRK4 method to approximate the time derivative. Due to the high flexibility and accuracy of MLS, it is considered to solve these traffic equations, which have shocks. Because they discretize the problem domain without using a predetermined grid, it seems to work well for equations that have shocks at points in the domain. On the other hand, contour integral and other suitable approaches are used for controlling the instability and challenges of ETDRK4. Now, by combining these methods, we get good results for solving these equations. In MLS, parameterc can be changed, and with it, the shape of the weight function can be controlled, and the type of weight function can be different. In ETDRK4, as mentioned, we used contour integral to control instability, but contour integrals are not the only solutions offered for this problem, but the contour integral method attracts us due to its generality for dealing with arbitrary functions. There is considerable flexibility with this procedure. Finally, two examples of traffic equations presented in this way were solved. The results showed that this method can successfully solve traffic problems.



---

*Source: 8208957-2022-11-03.xml*