¹

Keywords. Conservation laws on network, Chebyshev spectral method, Filters, Super Spectral Viscosity

1 Introduction↩︎

Modeling transport on networks is an important issue arising several contexts, such as for managing crowd evacuations in case of natural disasters (see [1]), to study the evolution of cracks in peridynamics and in fluid-dynamic models, (see for instance [2], [3]), to investigate infiltration in soils (see [4]–[6]) and for traffic management [7], [8], etc.

Several mathematical models have been proposed to deepen such issue, and, in particular, models consisting of systems of conservation laws are often used to describe such situations. They are characterized by the fact that their solutions may develop discontinuities in finite time even when the initial conditions are analytical (see [9]). Moreover, the problem is further complicated by the fact that, when conservation laws are applied on network, there is the need to establish a suitable transmission condition at the junction, (see [10]).

From a numerical point of view, systems of conservation laws are treated by using finite volume schemes, Glimm schemes (see [11]–[14]) or also mimetic finite difference methods (see [15]). They are first order monotone, consistent and stable schemes. A different approach consists in considering high-order spectral schemes. Indeed they are efficient method with exponential accuracy when applied to problems whose solutions have a certain amount of regularity. When solutions develop discontinuities in the domain, the exponential accuracy is lost and spectral methods show the Gibbs phenomenon, that is the formation of spurious oscillations in the profile of the solution. Such phenomenon can appear both at the boundaries of the domain or close to discontinuities. In the first case, one can overcome the phenomenon by choosing non-uniform mesh grid, with collocation points denser at the boundaries or by imposing periodic boundary conditions (see for instance [16]–[18]).

Other strategies have to be adopted in order to recover global high-order accuracy when damping appears in proximity of shocks. This represents an interesting challenge and a good framework for a machine learning approach to the problem. Indeed, in [19], the authors propose a projection-based reduced order model based on the Fourier collocation method for compressible flows and use a neural network to limit oscillations.

A different approach to recover spectral accuracy and to remove Gibbs phenomenon is based on the introduction of dissipative or filtered mechanism. Following this idea, in this paper we proposed a filtered Chebyshev spectral method to study numerically a system of conservation laws on network. The choice of using Chebyshev polynomials, instead of the Fourier ones, is related to the fact that Chebyshev interpolation does not require any periodicity assumption at the boundaries, and is also able to avoid Gibbs phenomenon at the boundaries thanks to the use of Chebyshev Gauss-Lobatto points. We filter the Chebyshev modes in the frequency space to reduce oscillations near discontinuities. We couple the method with a monotone transmission condition at the junction and we reformulate the methods by adding a super spectral viscosity in order to prove the convergence of the semi-discrete method, and, thanks to several simulations we show the gain of one order of convergence with respect to the implementation of a finite volume scheme.

The paper is organized as follows. In Section 2 we state the problem and provide a briefly review on conservation laws on network. In Section 3, we describe the filtered Chebyshev collocation method and prove its convergence. In Section 4 we construct the fully discrete method . Section 5 is devoted to the numerical simulations. Finally, Section [sec:concl] concludes the paper.

2 Statement of the problem↩︎

We consider a model of conservation laws describing the flow of a mass density on a star shaped network consisting of two incoming edges, one outgoing edge and one junction. The generalization to network consisting of more edges is trivial. Moreover, due to the property of finite speed of propagation in hyperbolic equations, it is enough to consider networks having a unique junction.

Each edge is parameterized by \(x\in\Omega_h\subset\mathbb{R}\), for \(h\in\{1,2,3\}\), so that the network is given by \(\Gamma=\Pi_{h=1}^3 \Omega_h\subset\mathbb{R}^3\) and the junction belongs to \(\bigcap_{h=1}^3\partial \Omega_h\), (see Figure 1 for a basic representation of the network). Let denote by \(\bar{x}\) the location of the junction.

On each edge, we consider an hyperbolic conservation law of the form \[\label{eq:cl} \frac{\partial u_h}{\partial t}(x,t) + \frac{\partial }{\partial x} f_h(u_h(x,t)) =0,\qquad t>0,\quad x\in\Omega_h,\quad h=1,2,3,\tag{1}\] where \(u_h\) is the density and \(f_h\) is the flux on each edge. We assume that the edges have a common maximal density \(u_{\max}>0\) and the fluxes are bell-shaped, Lipschitz and nonlinear non-degenerate, namely we require that \(f_h\) satisfies the following properties

• \(f_h\in Lip([0,u_{\max}];\mathbb{R}_+)\), with \(\left\| f_h \right\|_\infty\le L_h\), for certain \(L_h>0\),

• \(f_h(0)=f_h(u_{\max})=0\),

• there exists \(u_{h,c}\in(0,u_{\max})\) such that \(f_h'(u)\left(u_{h,c}-u\right)>0\) for a.e. \(u\in[0,u_{\max}]\),

• \(f_h'\) is not constant on any non-trivial subinterval of \([0,u_{\max}]\),

for every \(h=1,2,3\).

Under such assumptions on the flux, we have the existence of the strong traces for the solution (see [20], [21]).

To complete the model, we add the initial conditions \[\label{eq:initcond} u_h(x,0)=u_{h,0}(x),\quad x\in\Omega_h,\qquad h=1,2,3,\tag{2}\] and we assume that \(u_{h,0}\in L^2(\Omega_h)\cap L^\infty(\Omega_h)\), for \(h=1,2,3\).

We also require the conservation of the total density at the junction, namely for a.e. t>0, it holds \[\label{eq:consjunc} f_1\left(u_1\left(\bar{x}^-,t\right)\right)\;+ f_2\left(u_2\left(\bar{x}^-,t\right)\right) = f_3\left(u_3(\bar{x}^+,t)\right).\tag{3}\]

The well-posedness of the Cauchy problem is proved in [13] in the framework of admissible solution at the junction in terms of vanishing viscosity germ, (see [22]–[24]). In their work, the authors propose a finite volume scheme to approximate numerically the solution on the network and prove its convergence. While a validation of such a scheme can be found in [12].

The construction of an admissible weak solution is related with the definition of a Riemann solver at the junction, that is an operator which associates to any constant initial condition on the whole network a self-similar weak solution of 1 -2 .

Following [12], [13], we find that the selection of an admissible Riemann solver at the junction can be understood as a collection of initial boundary value problems of the form \[\label{eq:IBVP} \begin{cases} \frac{\partial }{\partial t} u + \frac{\partial}{\partial x}f_h\left(u_h\right)=0&\qquad (x,t)\in\Omega_h\times [0,T],\\ u_h(x,0)=u_{h,0}(x)&\qquad x\in\Omega_h\\ u_h(\bar{x},t)=u_{h,b}(t),&\qquad t\in[0,T], \end{cases}\qquad h=1,2,3,\tag{4}\] where the boundary functions have to be fixed in order to ensure the conservation at the junction 3 .

In [25], Bardos-LeRoux-Nédélec reformulate the conservation condition 3 in terms of traces of the solution and by using the Godunov flux.

The Godunov flux is the function which associates to any couple \((a,b)\in[0,u_{\max}]^2\) the value \(f_h(u_h(\bar{x}^-,t))=f_h(u_h(\bar{x}^+,t))\), where \(u_h\) is the Kruzhkov entropy solution to the Riemann problem (see [25]–[27]) \[\label{eq:Riemannpb} \begin{cases} \frac{\partial}{\partial t} u_h + \frac{\partial}{\partial x} f_h(u_h) = 0,&\quad x\in\mathbb{R},\, t>0,\\ u_h(x,0)=\begin{cases} a,\quad x<\bar{x},\\ b,\quad x>\bar{x}, \end{cases} &\quad x\in\mathbb{R}. \end{cases}\tag{6}\]

As proved in [25], the conservation condition 3 is equivalent to ask \[\label{eq:BCleft} f_h\left(u_h\left(\bar{x}^-,t\right)\right)=G_h\left(u_h\left(\bar{x}^-,t\right),\;u_{h,b}(t)\right),\quad h=1,2,\tag{7}\] and \[\label{eq:BCright} f_h\left(u_h\left(\bar{x}^+,t\right)\right)=G_h\left(u_{h,b}(t),\;u_h\left(\bar{x}^+,t\right)\right),\quad h=3.\tag{8}\]

If the system 1 admits a solution, then there exists \(u_b:\mathbb{R}_+\to[0,u_{\max}]\) such that (see [22], [28]) \[\label{eq:existub} u_{h,b}(t)=u_b(t),\qquad \text{for a.e. }t>0,\text{and for all }h=1,2,3.\tag{9}\]

Using 7 -8 -9 , the conservation condition 3 becomes \[\label{eq:BCcriterion} G_1\left(u_1\left(\bar{x}^-, t\right),\;u_b(t)\right) + G_2\left(u_2\left(\bar{x}^-, t\right),\;u_b(t)\right) = G_3\left(u_b(t),\;u_3\left(\bar{x}^+, t\right)\right).\tag{10}\]

As a consequence, studying the model 1 -2 -3 is equivalent to study \[\label{eq:model} \begin{cases} \frac{\partial }{\partial t} u + \frac{\partial}{\partial x}f_h\left(u_h\right)=0&\qquad (x,t)\in\Omega_h\times [0,T],\\ u_h(x,0)=u_{h,0}(x)&\qquad x\in\Omega_h\\ u_h(\bar{x},t)=u_{b}(t),&\qquad t\in[0,T], \end{cases}\qquad h=1,2,3,\tag{11}\] where for every \(t>0\), we have to compute the boundary value \(u_b(t)\) by solving the non-linear equation 10 .

3 Chebyshev collocation method↩︎

In what follows, we construct a spectral approximation of the solution to 11 -10 by using the discrete Chebyshev expansion.

We fix the total number of collocation points for each edge of the network \(N>0\) and we assume \(\Omega=[-1,1]^3\). So, the junction is located at \(\bar{x}_{in}=1\) for the incoming edges and at \(\bar{x}_{out}=-1\) for the outgoing one. By an affine linear transformation we can always consider the product of more general intervals.

We discretize each interval by the Chebyshev-Gauss-Lobatto points (CGL) defined as \(x_j=\cos\left(\frac{\pi j}{N}\right)\), for \(j=0,\dots,N\). They are non uniform mesh points, simmetrically distibuted with respect the average value of the reference domain. Due to their geometry, CGL points allows us to avoid Gibb’s phenomenon at the boundaries. However, as we will see in the next section, in order to remove such phenomenon near discontinuities, we need to introduce a filtered procedure.

On each edge, we seek an approximation of the solution \(u_h\) given by \[\label{eq:Chebyapprox} u_h^N(x,t)=\sum_{k=0}^N \tilde{u}_{h,k}(t) \;T_k(x),\qquad h=1,2,3,\tag{12}\] where for every \(k=0,\dots,N\), \(T_k(x)=\cos\left(k\arccos(x)\right)\) are the Chebyshev polynomials of the first order, and \(\tilde{u}_{h,k}(t)\) are the Chebyshev coefficients given by \[\label{eq:Chebycoef} \tilde{u}_{h,k} (t) = \frac{1}{\gamma_k}\sum_{n=0}^N u_h(x_n,t) T_k(x_n) w_n,\qquad h=1,2,3,\tag{13}\] where \[\label{eq:gamma} \gamma_k=\begin{cases} \pi&\qquad k=0,\,\text{or}\,k=N\\ \frac{\pi}{2}&\qquad k=1,\dots,N-1 \end{cases}\tag{14}\] and \[\label{eq:w} w_n=\begin{cases} \frac{\pi}{2N}&\qquad k=0,\,\text{or}\,k=N\\ \frac{\pi}{N}&\qquad k=1,\dots,N-1 \end{cases}\tag{15}\] Chebyshev polynomials are orthogonal with respect to the weight function \(w(x)=1/sqrt{1-x^2}\), moreover, thanks to their definition, Chebyshev polynomials can be seen, under a suitable change of variables, as trigonometric cosine functions. Therefore, we can compute them efficiently by using the Fast Fourier Transform algorithm (FFT).

In order to deduce a semidiscrete formulation of the model, we have to replace \(u_h(x,t)\) in 11 by its Chebyshev approximation \(u_h^N(x,t)\), for every \(h=1,2,3\). \[\label{eq:50}\frac{\partial}{\partial t} u_h^N(x,t) + \frac{\partial}{\partial x} I_N\left(f_h\left(u_h^N(x,t)\right)\right)=0,\tag{16}\] where \(I_N\) is the unique interpolation operator such that \(I_N u_h^N(x_j,t)=u_h(x_j,t)\), for \(j=0,\dots,N\) and \(h=1,2,3\).

We want to study 16 in the frequency space. To do this, we need to recall for the relationship between Chebyshev coefficients of a function and Chebyshev coefficients of its partial derivative of any order. It is easy to verify that the Chebyshev coefficients of \(\frac{\partial\;u_h^N}{\partial x}(x,t)\) are given by \[\tilde{u}^{\,'}_{h,k}=\sum_{n=0}^N D_{kn}\tilde{u}_{h,n},\] where \(D\) is the derivative matrix, whose entries are given by \[D_{k,n}=\begin{cases} n-1&\qquad k=1\,\text{and}\,n\,\text{even}\\ 2(n-1)&\qquad k>1\,\text{and}\,n>k\,\text{and} n+k\,\text{odd}\\ 0&\qquad \text{otherwise} \end{cases}\] Moreover, we can relate the coefficients of a product of two functions to those of these two functions.

We replace \(u_h(x,t)\) by \(u_h^N(x,t)\) in the model 11 , for every \(h=1,2,3\). Thus, if we make use of the previous results, at each collocation point, we can reformulate the model 16 in the frequency space as follows \[\label{eq:discmethodinterior} \frac{d}{dt}\tilde{u}_{h,k}(t) + \sum_{n=0}^N\sum_{m=0}^N\sum_{i=0}^N P_{knm}D_{mi}\tilde{f'}_{h,n}(t)\;\tilde{u}_{h,i}(t)=0, \qquad k=0,\dots,N,\quad h=1,2,3.\tag{17}\] where \(\tilde{f'}_{h,n}\) is the \(n\)-th Chebyshev coefficients of \(f'_h (u_h^N(x,t))\).

Using the same argument as before, we can approximate the initial conditions in the frequency space as follows \[\label{eq:discmethodinitcond} \tilde{u}_{h,k}(0) = \tilde{u_0}_{h,k},\qquad k=0,\dots N,\quad h=1,2,3.\tag{18}\]

To close the discrete model, we need to compute and approximate the boundary conditions 10 in the frequency space.

For every \(t>0\) we have to find a zero of a scalar nonlinear function. It is easy to prove that the nonlinear function in 10 is monotone and continuous but not everywhere differentiable since so is the Godunov flux. As a consequence, \(u_b(t)\) can be computed by using Regula Falsi method.

Once the boundary value has been found, we impose \[\label{eq:discmethodboundary} \sum_{k=0}^N \tilde{u}_{1,k}(t)=\sum_{k=0}^N \tilde{u}_{2,k}(t)=\sum_{k=0}^N (-1)^k\; \tilde{u}_{3,k}(t)=u_b(t),\qquad t>0,\tag{19}\] so that, in the physical space we have \[\label{eq:bcphysical} u_1\left(x_N,t\right)=u_2\left(x_N,t\right)=u_3\left(x_0,t\right)=u_b(t),\qquad t>0.\tag{20}\]

3.1 Filtering for Chebyshev method↩︎

Due to their high-order accuracy, spectral methods, when applied to nonlinear problem, may cause artificial numerical oscillations, known as Gibbs phenomenon. Such phenomenon can reduce drastically the accuracy of the method.

The formulation and the analysis of high-order accurate methods able to do not generate Gibbs phenomenon is an interesting challenge. Several results has been obtained in order to better understand the problem.

Filtering is a strategy to limit oscillations in the solution [29]. It has been introduced in the framework of signal process with the aim to reduce the noises in a process. In our context, we implement a filter to make a specific modification of the Chebyshev coefficients in order to have a control of the oscillations in the profile of the solution.

In particular, a low-order filter dampens high-frequency components, while a high-order filters are able to localize oscillations close to discontinuities.

In what follows we consider an exponential filter function defined as \[\label{eq:sigma} \sigma(\eta,p)=\begin{cases} \exp\left(-\beta\;\eta^p\right),&\quad 0\le\eta\le1\\ 0,&\quad \text{otherwise} \end{cases}\tag{21}\] where \(p\) refers to the order of the filter, and \(\beta\) is the damping rate. To ensure that the exponential filter 21 satisfies Definition 2, we choose \(\beta=-\log(\varepsilon_M)\), where \(\varepsilon_M\) is the machine accuracy. In Figure 2 we show the behavior of exponential filters of different orders.

To apply the filtering strategy to our model, we modify the \(k\)-th Chebyshev coefficient defined in 13 by pre-multiplication with \(\sigma\left(\frac{k}{N}\right)\), namely the \(k\)-th filtered coefficients becomes \[\label{eq:FChebycoef} \tilde{u}_{h,k}^{\,\,\sigma}(t)=\sigma\left(\frac{k}{N}\right)\;\tilde{u}_{h,k}(t).\tag{22}\] Thus, the semi-discrete method 17 becomes \[\label{eq:filterdiscmethodinterior} \frac{d}{dt}\tilde{u}_{h,k}^{\,\,\sigma} (t) + \sum_{n=0}^N\sum_{m=0}^N\sum_{i=0}^N P_{knm}D_{mi}\tilde{f'}_{h,n}^{\,\,\sigma} (t)\;\tilde{u}_{h,i}^{\,\,\sigma}(t)=0, \qquad k=0,\dots,N,\qquad h=1,2,3.\tag{23}\]

3.2 Super Spectral Viscosity↩︎

When spectral methods are applied to nonlinear hyperbolic equations in conservation form, the problem of obtaining an entropy satisfying solution arises. Unmodified spectral methods do not converge to the correct entropy solution if the solution contains shocks [30].

Spectral calculations can be stabilized by using exponential filters on the conserved variable at each time step. An equivalent approach consists in applying a Super Spectral Viscosity (SSV) in the discretized model. Indeed, a SSV method converges to the correct entropy solution and maintains the spectral accuracy, when the solution is uniformly bounded.

In this section we introduce the Super Spectral Viscosity in the semi-discrete scheme and we show the equivalence of such approach with the implementation of an exponential filter. As we will see in the next section, the introduction of the SSV method allows us to prove the convergence of the semi-discrete scheme.

We defined the Super Spectral Viscosity for the Chebyshev collocation method as follows \[\label{eq:SSV} SV(u_h^N) = \frac{\varepsilon(-1)^{s+1}}{N^{2s-1}}\left(\sqrt{1-x^2}\frac{\partial}{\partial x}\right)^{2s} \left(u_h^N\right) = \frac{\varepsilon(-1)^{s+1}}{N^{2s-1}}Q^{2s} \left(u_h^N\right), \quad h=1,2,3,\tag{24}\] where \(\varepsilon>0\) is the super viscosity coefficient, and \(s\) is an integer growing with \(N\).

Thus, the Chebyshev collocation method can be written in the following way \[\label{eq:ssvmethod} \frac{\partial}{\partial t} u_h^N(x,t) + \frac{\partial}{\partial x} I_N\left(f_h\left(u_h^N\right)\right) = \frac{\varepsilon(-1)^{s+1}}{N^{2s-1}}\left(\sqrt{1-x^2}\frac{\partial}{\partial x}\right)^{2s} \left(u_h^N\right),\quad h=1,2,3.\tag{25}\]

The reformulation of the filtered Chebyshev method by using Super Spectral Viscosity, allows us to prove the convergence of the semi-discrete scheme.

3.3 Convergence of the semi-discrete scheme↩︎

We consider the spectral approximation 25 of 11 , with initial condition \(u_{0,h}^N\), \(h=1,2,3\). For simplicity, we focus on the case \(s=1\), but the results of this section can be easily generalized to \(s>1\).

Due to the Chebyshev approximation properties, we can assume \[\label{eq:assump1} \begin{align} &u_{0,h}^N \to u_{0,h}\quad\text{in}\quad L^p_{loc}(\Omega_h),\quad 1\le p <\infty,\quad\text{as}\quad N\to\infty,\\ &\left\| u_{0,h}^N \right\|_{L^\infty\left(\Omega_h\right)}\le \left\| u_{0,h} \right\|_{L^\infty\left(\Omega_h\right)},\quad \left\| u_{0,h}^N \right\|_{L^2\left(\Omega_h\right)}\le \left\| u_{0,h} \right\|_{L^2\left(\Omega_h\right)},\quad N>0. \end{align}\tag{28}\]

Additionally, the maximum principle for parabolic equations applied to 25 ensures \[\label{eq:assump2} \left\| u_{h}^N \right\|_{L^\infty\left(\Omega_h\times (0,\infty)\right)} \le \left\| u_{0,h} \right\|_{L^\infty\left(\Omega_h\right)},\quad N>0,\quad h=1,2,3.\tag{29}\]

We prove the following preliminary Lemmas.

We can now prove the convergence of the semi-discrete scheme.

4 Fully-discrete scheme↩︎

Once the Cauchy problem on network 11 has been discretized in space, we need to integrate in time the obtained system of ODEs 17 . In this paper, we construct the fully-discrete scheme by implementing the Midpoint method.

We fix the time step \(\Delta t>0\) satisfying the CFL condition \[\label{eq:CFL} \Delta t \max_h \left(L_h\right)\le \frac{\min_{j=0,\dots,N-1} \left|x_{j+1}-x_j\right|}{2},\tag{39}\] and for \(s\in\mathbb{N}\) we define the time mesh \(t^s=s\Delta t\). At each time step \(t^s\), \(\tilde{u}_{h,k}^{\, s}\) represents an approximation of the solution \(\tilde{u}_{h,k}(t^s)\), for every \(k=0,\dots,N\) and \(h=1,2,3\).

Let us introduce the following notation for the second term of the left hand side of 23 \[\label{eq:discflux} F_{h,k}^\sigma\left(t^s,\tilde{u}_h^\sigma(t^s)\right) = \sum_{n=0}^N\sum_{m=0}^N\sum_{i=0}^N P_{knm}D_{mi}\tilde{f'}_{h,n}^{\,\sigma}(t^s)\;\tilde{u}_{h,i}^{\,\sigma}(t^s)=0, \qquad k=0,\dots,N,\qquad h=1,2,3.\tag{40}\]

The fully-filtered discrete scheme is given by \[\label{eq:fullymethod} \tilde{u}_{h,k}^{\,s+1} = \tilde{u}_{h,k}^{\,s} - \Delta t F_{h,k} \left(t^s+\frac{\Delta t}{2}, \tilde{u}_{h} \left(t^s+\frac{\Delta t}{2}\right)\right),\qquad k=0,\dots, N,\qquad h=1,2,3.\tag{41}\]

4.1 A 2D Chebyshev collocation method for the fully discretized model↩︎

In order to obtain a second order scheme in time, one can use a spectral Chebyshev method both for space and time discretization, as shown in [32].

The method consists in looking for an approximation of the solution on the whole network in the form of finite combination of the product of Chebyshev polynomials in space and time variables.

More precisely, for simplify the notation, we assume that the time interval under consideration is \([-1,1]\), so that the initial condition corresponds to \(t=-1\). We consider as full domain of computation the product \(\tilde{\Omega}=\Omega\times[-1,1]\) and we look for \(u^N_h(x,t)\approx u_h(x,t)\), \(h=1,2,3\) given by \[\label{eq:chebser2d} u^N_h(x,t)=\sum_{j=0}^{N_x}\sum_{k=0}^{N_t} \tilde{u}_{h,jk}\;T_j(x)\;T_k(t)\tag{42}\] where \(N_x\) and \(N_t\) represent the total number of collocation points in space and time, respectively, and \(N=(N_x,N_t)\).

The coefficients \(\tilde{u}_{h,jk}\), \(j=0,\dots,N_x\) , \(k=0,\dots,N_t\), and \(h=1,2,3\), are the Chebyshev coefficients of the discrete Chebyshev expansion and when the grid points are the Chebyshev Gauss-Lobatto points \((x_n,t_m)=(\cos(n\pi/N_x),\cos(m\pi/N_t))\), then their expression is given by \[\label{eq:fcoef} \tilde{u}_{h,jk}=\frac{1}{\gamma_j\gamma_k}\sum_{n=0}^{N_x}\sum_{m=0}^{N_t} u_h(x_n,t_m)\;T_j(x_n)T_k(t_m)w_n w_m,\tag{43}\] where \(\gamma\) and \(w\) are defined in 14 and 15 , respectively.

For the purpose of our work, we show here the expression of the expansion of the first order derivative in time and in space: \[\label{eq:ptchebser} \begin{align} \frac{\partial }{\partial t}u^N_h(x,t)&=\sum_{j=0}^{N_x}\sum_{k=0}^{N_t}\sum_{\ell=0}^{N_t} D_{k\ell}^{(t)}\;\tilde{u}_{h,j\ell}\;T_j(x)\;T_k(t),\\ \frac{\partial }{\partial x}u^N_h(x,t)&=\sum_{j=0}^{N_x}\sum_{k=0}^{N_t}\sum_{i=0}^{N_x} D_{ji}^{(x)}\;\tilde{u}_{h,ik}\;T_j(x)\;T_k(t) \end{align}\tag{44}\] where \(D\) is the derivative matrix, and the superscripts \((t)\) and \((x)\) in such a matrix denote the differentiation with respect to the temporal or spatial coordinates, respectively.

Moreover, it is possible to find a compact expression for the product of two functions as linear combination of Chebyshev polynomials in the same way as in Lemma 1. We have that the \(ij\)-th Chebyshev coefficient of the product between two functions \(f\) and \(g\) is given by \[\label{eq:cheb2dproduct} \tilde{\left(fg\right)}_{ij}=trace\left(P^{(x)}_i\tilde{f}\;\tilde{g}^{\, T}\;P^{(t)}_j\right),\qquad i=0,\dots, N_x,\quad j=0,\dots,N_t+1,\tag{45}\] where \(\tilde{f}\) and \(\tilde{g}\) are the \((N_x+1)\times (N_t+1)\) matrix of Chebyshev coefficients of \(f\) and \(g\), while \(P^{(x)}\) and \(P^{(t)}\) are the \((N_x+1)\times(N_x+1)\times(N_x+1)\) and \((N_t+1)\times(N_t+1)\times(N_t+1)\) tensors of Lemma 1, respectively.

The advantage of implementing such scheme is related to the fact that even in this case we can benefit of the implementation of the 2D-FFT to compute the Chebyshev coefficients \(\tilde{u}_{h,jk}\) in 42 .

If we replace \(u_{h}(x,t)\) with its Chebyshev approximation \(u_h^N(x,t)\), we find \[\begin{align} \label{eq:2Dapprox} 0=&\frac{\partial}{\partial t} u_h^N(x,t) + \frac{\partial }{\partial x} I_N\left(f_h(u_h^N(x,t))\right)\\ =&\sum_{j=0}^{N_x}\sum_{k=0}^{N_t}\left(\sum_{\ell=0}^{N_t} D_{k\ell}^{(t)}\;\tilde{u}_{h,j\ell} + a_{h,jk}\right)\;T_j(x)\;T_k(t), \end{align}\tag{46}\] where for every \(h=1,2,3\), \(j=0,\dots,N_x\) and \(k=0,\dots,N_t\), \(a_{h,jk}\) is given by \[\label{eq:a} a_{h,jk}=trace\left(P_j^{(x)}\tilde{(f')}_h\left(\sum_{\ell =0}^{N_x} D^{(x)}_{j\ell}\tilde{u}_{h,\ell k}\right)^{\,T}\;P_k^{(t)}\right).\tag{47}\]

The unknowns \(\tilde{u}_{h,jk}\) have to satisfy the conservation law in the interior of the integration domain with special care for imposing the initial and boundary conditions. We have \[\begin{align} \label{eq:initcond2D} u^N_{h}(x,-1)&=\sum_{j=0}^{N_x}\sum_{k=0}^{N_t} (-1)^k \tilde{u}_{h,jk} T_j(x)\\ &=\sum_{j=0}^{N_x} \tilde{u_0}_{h,j} T_j(x)\qquad\qquad\qquad \qquad\qquad \qquad \qquad\qquad h=1,2,3.\\ &=u_{h,0}(x) \end{align}\tag{48}\] The boundary conditions for the incoming edges can be obtained by solving the system \[\label{eq:BCincom2D} \sum_{j=0}^{N_x} \tilde{u}_{h,jk} = \tilde{(u_b)}_k\qquad h=1,2,\quad k=0,\dots,N_t.\tag{49}\] While the boundary condition for the outgoing edge writes \[\label{eq:BCoutg2D} \sum_{j=0}^{N_x} (-1)^k \;\tilde{u}_{h,jk} = \tilde{(u_b)}_k\qquad h=3,\quad k=0,\dots,N_t.\tag{50}\]

We can notice that, in order to find the solution on the whole network, we need to impose the boundary condition for every time step \(t_k\), \(k=0,\dots,N_t\). However, the proposed Riemann solver used to select a unique weak solution to 1 , requires to compute the boundary values implicitly at each iteration. So, in our analytical framework, a 2D Chebyshev spectral method is not suitable to get the admissible weak solution.

Such scheme can be applied successfully to each edge, but it needs a further analytical investigation to be couple with a different transmission condition. This represents an ongoing study, which will be presented in a future work.

5 Numerical Simulations↩︎

In this section we present some simulations in order to show the properties of the proposed filtered Chebyshev spectral method applied to the network.

More in details, we make a validation of the scheme by comparing the approximated solution provided by the method with an exact solution on network computed explicitly in [12]. Then, we compare the filtered Chebyshev scheme with the finite volume scheme proposed in [12] in terms of CPU cost.

Finally, in order to verify that, at least far from the junction, even 2D-Chebyshev method constructed in 4 performs well, in the framework of hyperbolic conservation laws.

5.1 Validation of the scheme on network↩︎

We fix \(N>0\) and discretize \(\Omega=[-1,1]^3\) by using the CGL points. The validation of the filtered Chebyshev method on network is made by comparing the obtained approximated solution with the explicit solution computed in [12], with flux \(f_h(u) = u\left(1-u\right)\), for every \(h=1,2,3\), and with initial condition \[u_{1,0}(x) = \chi_{[-1/2,0]}(x),\quad u_{2,0}(x)=\frac{3}{4} \chi_{[-1/4,0]}(x),\quad u_{3,0}(x)=0.\]

Figure 3 compare the profile of the numerical solution on each edge with the exact one, showing a good agreement and the capability of the spectral method to capture well the shock without showing the appearance of spurious oscillations near discontinuities. The result is also confirmed in terms of convergence rate.

We introduce the relative \(L^1\)-error on network as follows \[E^s_{L^1}=\frac{\sum_{h=1}^3 \sum_{j=0}^N \left|u^N_h(x_j,t^s)-u^\ast_h(x_j,t^s)\right|}{\sum_{h=1}^3\sum_{j=0}^N \left|u^\ast(x_j,t^s)\right|},\] where \(u^\ast\) is the reference solution.

Table 1 shows the relative error \(E^{s}_{L^1}\) corresponding to the incoming edges of the network for different values of the total number of mesh points \(N\) at time \(t^s=1/6\) and the associated convergence rate. The table shows also a comparison between the error and the convergence rate obtained with the finite volume scheme (FVS) proposed in [12].

We can conclude that thanks to spectral method we can gain one order of accuracy, even in presence of shocks in the profile of the solution.

5.2 Computation of the execution time for a Riemann problem at the junction↩︎

We analyze here the performance of the method applied to a Riemann problem on network in terms of the computational cost required to complete the simulation and compare the results with the ones obtained by using finite volume schemes.

On each edge we consider the same flux \(f_h(u)=u\left(1-u\right)\) and fix constant initial conditions: \[u_{1,0}(x) = \frac{1}{4},\quad u_{2,0}(x) =\frac{2}{3},\quad u_{3,0}(x) = \frac{4}{5}.\]

In Table 2 we find that, even if the method allows us to have a better accuracy with respect to the implementation of (FVS), it takes much more time to complete the simulation. The reason of such result could be related to the fact that, although FFT allows us to reduce the computation cost to obtain Chebyshev coefficients, the method has to take into account the cost of multiplications in the right hand side of 23 .

5.3 Implementation of the 2D-Chebyshev spectral method↩︎

In this section, we perform a simulation on a degenerate network consisting in one single edge without any junction. This particular situation does not require the computation of the boundary data to establish the transmission condition, thus represents a good framework in order to validate the 2D-Chebyshev method described in 4.1.

We consider the same flux as in the previous simulations and we fix \([-1,1]\) as time domain, so that the initial condition has to be considered at time \(t=-1\). We take \[u_0(x)=\chi_{[-3/4,-1/4]}(x) -\frac{2}{3}\left(x-\frac{5}{2}\right)\chi_{[-1/4,1/2]}(x)\] as initial condition.

In Figure 4, we make a comparison between the profile of the solution at time \(t=1\) computed by means of the filtered Chebyshev method and the 2D-Chebyshev method. As expected, We can observe a perfect agreement in the results.

A difference in the performance can be found in terms of execution time. Indeed, midpoint scheme seems more competitive, as shown in Table 3.

6 Conclusions↩︎

The paper focuses on the implementation of a filtered Chebyshev spectral scheme to solve numerically a system of conservation laws on network. The modification of the Chebyshev coefficients by multiplication with a filter function allows us to avoid the Gibbs phenomenon near shock discontinuity. As a consequence, the method can benefit of the high-accuracy properties of spectral methods. We prove the convergence of the semi-discretized method and we perform some simulations in order to make a comparison between the proposed scheme and the implementation of a finite volume scheme. Our tests show the gain of one order in the convergence rate, even if the filtered Chebyshev method requires an higher computational cost. We also investigate the performance of a Chebyshev method to discretize the problem both in the space and in the time variable. This approach is applied only to a single edge is able to approximate well the solution of the problem, but it seems expensive from a computational point of view. Moreover, it need to be coupled with a different Riemann solver, because it requires to know in advance the values of the boundary data for every time. This issue represents the starting point for a further investigation.

Acknowledgments↩︎

The author is member of the INdAM Research group GNCS. She has been supported by REFIN Project, grant number D1AB726C funded by Regione Puglia, by INdAM - GNCS 2023 Project, grant number CUP\(\_\)E53C22001930001 and by PNRR MUR - M4C2, grant number N00000013 - CUP D93C22000430001. She also acknowledges the partial support of “Finanziamento giovani ricercatori 2022” funded by GNCS-INdAM.

References↩︎