Convergence of a Lagrangian–Eulerian scheme via weak asymptotic analysis for one-dimensional hyperbolic problems

In this work, we present a fully discrete Lagrangian–Eulerian scheme, and its numerical analysis by a weak asymptotic method, for the computational treatment of first-order hyperbolic partial differential equations given by

where $u=u(x,t):\mathbb{R}\times\mathbb{R}^{+}\to\Omega\subset\mathbb{R}$ and $H:\Omega\to\mathbb{R}$.

To give a brief summary about the previous works that lead us to obtain the new method developed in this paper, we cite first [16]. In that article, an innovative locally conservative scheme was developed for the scalar parabolic convection-dominated convection–diffusion model problem of two-phase, immiscible, incompressible flow in porous media, the Locally Conservative Eulerian–Lagrangian Method (LCELM). This scheme was related to the divergence form of the parabolic flow equation, where the use of a space-time divergence form allowed the localization of the transport and the desired local conservation property could also be localized. Such a scheme was very competitive from a computational point of view for scalar nonlinear transport problems in heterogeneous porous media. To the best of our knowledge, the LCELM procedure was the first work in the literature to introduce a space-time local conservation (by physical and geometric arguments), a region where the mass conservation takes place (locally), which was coined as integral tube along with the so-called integral curves for scalar parabolic convection-diffusion models in porous media transport problems (see Eqs. (5.4a)-(5.4b) in [16]). Moreover, so-called integral curves were associated with points on the boundary (usually vertices) of the finite elements in that framework.

In this work, the fully discrete Lagrangian–Eulerian formulation is based on the new and substantial improvement interpretation of the integral tube, which is now subject to condition $O\left(\frac{H(u)}{u}\right)\propto\left[\frac{\Delta x}{\Delta t}\right]\to 0$ (called no-flow curves [7]), where the quantities $u$ and $H(u)$ come from the scalar initial value problem (1.1) and under a suitable stability estimate that is very effective in computing practice; as a result, we also obtain a weak CFL condition that does not depend on the derivative of the flux function $H(u)$, but only on the introduced no-flow curves. This simple and interesting technique is the key ingredient of the Lagrangian–Eulerian framework that provides information about the local wave propagation speed. This issue is not discussed in [16]. In addition, the no-flow curves reveal to be also a desingularization analysis tool for the construction of computationally stable schemes [6, 5, 7]. In [7], the authors also discussed a new reinterpretation of the Lagrangian–Eulerian no-flow curves as an anti-diffusion term into the viscosity coefficient defined by the quantity $O\left(\frac{H(u)}{u}\right)$, with a distinct identification to the Finite Volume framework. Here we are interested in the design of novel fully-discrete schemes based on the new concept no-flow curves subject to $O\left(\frac{H(u)}{u}\right)\propto\left[\frac{\Delta x}{\Delta t}\right]\to 0$ which is substantially different in theory foundations from the previous and relevant LCELM method.

The idea of how the fully discrete explicit locally conservative Lagrangian–Eulerian scheme works is outlined in the following three basic steps, for each time step, from $t^{n}$ to $t^{n+1}$. First, by using a space-time divergence form of model (1.1) we obtain an exact construction of the set of ODE modeling the so-called no-flow curves (see Figure 1). Second, by integrating the underlying conservation law (in divergence form) over the special region in the space-time domain (see region $D_{j}^{n,n+1}(x,t)$ in Figure 1), where the conservation of the mass flux takes place, we get the Lagrangian–Eulerian conservation relations. Third, by combining steps one and two, we perform an evolution (Lagrangian) step in time, from a primal grid to a nonstaggered grid and turn back by means of the final (Eulerian) projection step with a piecewise linear reconstruction (see Figure 2). A key hallmark of such method is the dynamic tracking forward of no-flow curves subject to $O\left(\frac{H(u)}{u}\right)\propto\left[\frac{\Delta x}{\Delta t}\right]\to 0$, per time step. The scheme is free of local Riemann problem solutions and does not use adaptive space-time discretization. This is a considerable improvement compared to the classical backward tracking over time of the characteristic curves over each time step interval, which is based on the strong form of the problem and that are not reversible for systems in general.

To analyze the properties of the explicit Lagrangian–Eulerian scheme, we use recent improvements on the weak asymptotic analysis (see [2, 3, 10]), which was first defined in a distinct framework by [13]. The weak asymptotic solutions are consistent with the traditional solutions in one-dimensional and multi-D (see also [2, 3, 5]). The weak asymptotic analysis has been used to study the existence of solutions for scalar equations and systems of hyperbolic equations, giving solutions a new meaning (see [2, 3, 11, 14, 12, 13, 15, 21, 22, 23, 25] and the references therein). An interesting aspect of this theory is that it makes it possible to prove the existence (and, for the scalar case, the uniqueness) of a solution by means of numerical methods.

In a previous work [5], we defined a weak asymptotic solution for a scalar equation (1.1) and outlined the proof of stability of the numerical method used. Weak asymptotic methods aim to investigate the nonlinear phenomena that appear in evolutionary equations (see [2, 3, 10]). In addition, they have been used in explicit calculations by several authors to study wave interactions inside the solutions to Riemann or Cauchy problems when these solutions involve nonclassical products of Heaviside functions, $\delta$-Dirac distributions, and even their derivatives. In [2], the authors show how one can construct families of continuous functions which asymptotically satisfy scalar equations with discontinuous nonlinearities and systems having irregular solutions. Through a weak asymptotic analysis, it has been proven, for scalar equations, that the initial value problem is well posed in the $L^{1}$ sense for the approximate solutions constructed. It was possible to demonstrate that the family of solutions obtained from the method satisfies Kruzhkov entropy, which provided a better understanding of the mathematical computations for the construction of accurate numerical schemes satisfying classical and entropy (Kruzhkov) solutions. It turns out that, beside the applications in nontrivial models of hyperbolic conservation laws, the convergent Lagrangian–Eulerian scheme via the weak asymptotic method treat such models with computational efficiency. From this theory, we obtain explicit estimates for the stability condition, and these estimates are numerically implemented and give very good numerical solutions as we see in the Section 4.

The proposed weak asymptotic analysis in this work, together with the Lagrangian–Eulerian approach [4, 5, 6, 7], fits in properly with the classical theory while improving the mathematical computations for the construction of new accurate numerical schemes, given by Proposition 3.1, where convergence, existence, and stability are established. In the classical treatment for proving convergence of a numerical scheme, first it is proved that the approximate solutions, generated by the numerical scheme, are of finite total variation (or bounded variation), actually, for scalar equations are proved that the total variation are not increasing. This condition gives the compact embedding of $BV$ functions in $L^{1}$; one can use the Helly’s selection theorem to show pointwise convergence of the generated sequence of solutions. After, it is proved that the generated sequence is a weak solution for the scalar equations and, finally, it is proved that the solution satisfies an entropy criterion, the most commom used is the Kruzhkov entropy condition, see [18]. The treatment used in the weak asymptotic analysis is very similar to the previous steps. To prove the main properties of the numerical scheme it is proved that the numerical scheme generates a solution with bounded total variation, also it is used a similar argument of Helly’s selection theorem, see Appendix C. Using the asymptotic theory, the solution to the proposed scheme is obtained as a family of functions, $\{u(\cdot,t,\epsilon)\}_{\epsilon}$, bounded in $\mathbb{L}^{1}(\mathbb{S}^{1})$ uniformly in parameter $\epsilon$ of the model (1.1)). The main difference is the kind of limit that is taken in solution in Equation $(\ref{weaksol})$. For instance, our approach allows us to deal with the reconstruction (accomplished by means of robust choices of slope limiters) of variable $u$ of Eq. $(\ref{noLi})$ in the (numerical) flux terms. We also give sufficient conditions for a total variation nonincreasing (Section 3.1) through a suitable semidiscrete formulation of Eq. $(\ref{noLi})$. In addition, we obtain the maximum principle and the entropy (Kruzhkov) solution to model $(\ref{noLi})$ thanks to a suitable interpretation of the approximate solutions provided by the analysis presented in Section 3.2. The weak asymptotic theory handles reconstructions easily and opens new possibilities for the design of new methods for hyperbolic problems under a weak CFL-type stability depending on the no-flow curves associated with the novel Lagrangian–Eulerian approach instead of using estimates on the eigenvalues, the weak asymptotic theory fits very well with the Lagragian-Eulerian scheme and it is a very promissor technique to be used to prove convergence for improved numerical methods of Lagrangian-Eulerian class.

The rest of this paper is structured as follows. Section 2 introduces the Lagrangian–Eulerian scheme. Section 3 presents the convergence of the numerical method via the weak asymptotic analysis. We propose a scheme to find a solution to $(\ref{noLi})$ and prove that the resulting solution converges to the weak solution of $(\ref{noLi})$. We demonstrate that the numerical scheme obtained in Section 2 is compatible with that used in the weak asymptotic analysis. Section 3.1 proves that our scheme has a total variation nonincreasing property for the solution. To complete our analysis, Section 3.2 demonstrates that the obtained solution satisfies the maximum principle and Kruzhkov entropy solution. Appendix A shows evidence that the reconstructions used in this study are Lipschitz continuous. Section 4 shows and discusses a set of representative computational results, including numerical studies with the W1 distance. Section 5 summarizes our concluding remarks.

To construct the fully discrete Lagrangian–Eulerian scheme, we first consider the scalar one-dimensional (1D) Cauchy problem (1.1) in its divergent form

For the sake of clarity and completeness, we sketch in what follows the main three steps of the Lagrangian–Eulerian approach, namely, the construction of the set of ODE modeling the no-flow curves (Section 2.1), the Lagrangian–Eulerian conservation relations (Section 2.2) and the final projection step by piecewise linear reconstruction (Section 2.3), along with an evolution algorithm of the fully discrete nonstaggered Lagrangian-Eulerian method. As in the Lagrangian–Eulerian schemes present in [4, 5, 6, 7], local conservation is obtained by integrating the conservation law over the region in the space-time domain where the conservation of mass flux takes place.

In the classical treatment for proving convergence of a numerical scheme, first it is proved that the approximate solutions, generated by the numerical scheme, are of finite total variation (or bounded variation), actually, for scalar equations are proved that the total variation are not increasing. This condition gives the compact embedding of $BV$ functions in $L^{1}$, one can use the Helly’s selection theorem to show compactness for the approximations and establish pointwise convergence. After, it is proved that the solution is a weak solution for the scalar equations and, finally, it is proved that the solution satisfies an entropy criterion, the most commom used is the Kruzhkov entropy condition, see [18]. In this Section, we show the convergence of our scheme using the weak asymptotic theory, that is very suitable to handle with the Eulerian-Lagrangian approach. From this theory, we obtain explicit estimates for the stability condition, these estimates are implemented and give very good numerical solutions as one can see in the Section 4.

To establish the weak asymptotic solution, we consider the one-dimensional scalar equation in (1.1). In order to avoid boundary conditions in the bounded domain due to numerical purposes, we consider that $x\in\mathbb{S}^{1}=\mathbb{R}/\mathbb{Z}$; variable $t\in\mathbb{R}^{+}$, thus $u=u(x,t):\mathbb{S}^{1}\times\mathbb{R}^{+}\to\Omega\subset\mathbb{R}$, and the flux function $H=H(u):\Omega\to\mathbb{R}$ is assumed to be a locally Lipschitz function in $u$, i.e.,

Assumption 1 - for all $c>0$, $\exists L>0$, such that

The weak asymptotic solution is a sequence of solutions $(u_{\epsilon})_{\epsilon}=(u(x,t,\epsilon))_{\epsilon}$ of class $\mathcal{C}^{1}$ in $t$ and class $L^{\infty}$ and is piecewise continuous in $x$ such that, for all $\psi=\psi(x)\in\mathcal{C}^{\infty}_{c}(\mathbb{R})$ (smooth with compact support), and, for all $t$,

In the weak asymptotic method, a PDE with a special flux (using parameter $\epsilon$) is first proposed. Then, for each fixed $\epsilon$, we obtain an Ordinary Differential Equation (ODE) for variable $t$. Based on the theory of ODEs, we prove the existence and stability of the solution. Finally, we demonstrate that, when taking $\epsilon\to 0$, the limit satisfies (3.2). The idea is for the flux to represent the numerical method, which is why the existence and stability of the PDE with special class may be an extension of the numerical method. In our method, we use an auxiliary function ($f(u)=H(u)/u$) and assume that $u\neq 0$ to avoid technical details (we would like to stress that the convergence of the method can be proven in this case). Moreover, we assume that there exists a $a>0$ such that $u>a>0$. Note that, in this case, $f(u)$ is a locally Lipschitz function in $u$ because for all $c>a$, $\exists\tilde{K}>0$ such that

Notice that we can take $\tilde{K}=L/a+\tilde{M}_{c}/a^{2}$, where $\displaystyle{\tilde{M}_{c}=\max_{a<u\leq c}|H(u)|}$.

For our method, we propose the following ODE:

with initial condition $u_{\epsilon}(x,0)=u_{0}(x,0)$. Here $f^{+}$ and $f^{-}$ are defined in $(\ref{fmm})$ and $u_{{}_{\epsilon-i}}$, $u_{{}_{\epsilon+i}}$, and $u_{{}_{\epsilon}}$ are defined as $\displaystyle{u_{{}_{\epsilon-i}}=u(x-i\epsilon,t,\epsilon)}$, $\displaystyle{u_{{}_{\epsilon+i}}=u(x+i\epsilon,t,\epsilon)}$ and $\displaystyle{u_{{}_{\epsilon}}=u(x,t,\epsilon)}$. Functions ${f^{+}\left(\hat{u}_{{}_{\epsilon-1/2}}\right)}$ and $f^{-}\left(\hat{u}_{{}_{\epsilon-1/2}}\right)$ are evaluated in the middle of each cell, i.e, they are given by following identities ${f^{+}\left(\hat{u}_{{}_{\epsilon-1/2}}\right)=f^{+}\left(\hat{u}_{{}_{\epsilon-1/2}},x-\frac{\epsilon}{2},t\right)}$,
${f^{+}\left(\hat{u}_{{}_{\epsilon+1/2}}\right)=f^{+}\left(\hat{u}_{{}_{\epsilon+1/2}},x+\frac{\epsilon}{2},t\right)}$, ${f^{-}\left(\hat{u}_{{}_{\epsilon-1/2}}\right)=f^{-}\left(\hat{u}_{{}_{\epsilon-1/2}},x-\frac{\epsilon}{2},t\right)},$ and ${f^{-}\left(\hat{u}_{{}_{\epsilon+1/2}}\right)=f^{-}\left(\hat{u}_{{}_{\epsilon+1/2}},x+\frac{\epsilon}{2},t\right)}$.

States $\hat{u}_{{}_{\epsilon-1/2}}$ and $\hat{u}_{{}_{\epsilon+1/2}}$ are obtained from a combination of known states from equations such as $(\ref{umeio})$ and are also determined in the middle of each cell. For the general case, we assume that there is a function $L(x_{-p},\dots,x_{p})$ so that we can define, for instance, $\hat{u}_{{\epsilon+1/2}}$ as:

Function $L$ is generated using the slope limiters (defined in Section 2.4). In addition, we assume the following compatibility condition: for any $x$, $x=L(x,x,\dots,x).$ In Appendix A, we prove that the slope limiters, as well as the reconstructions of type $(\ref{umeio})$, are Lipschitz continuous. In this case, if $u(x+i\epsilon,t,\epsilon)$ are continuous functions

Note that since $f(u)$ and $L$ are Lipschitz continuous, then $f$ applied to $\hat{u}_{{}_{\epsilon-1/2}}$ is also a Lipschitz function.

In the next result, we state the existence and stability result of the solution to $(\ref{ODE})$. As strategy of proof, we apply the Taylor expansion with remainder term to substitute the derivative by a difference defined between $t$ and $t+dt$, where $dt$ represents a time step. We state our result as follows:

Notice that $(\ref{Cflcon})$ is always valid if the CFL condition $(\ref{CFL})$ is satisfied. Indeed, the proof of Proposition 3.1 follows ideas from works [2, 3, 10], however, here is the first time that we prove the result for reconstructions of variable $u(x,t)$. Moreover, with the adaptation of our proof, we are able to obtain new conditions to guarantee the stability of the numerical method. These conditions are not obtained in previous works.

Proof. First, we fix $\epsilon$ and obtain an ODE from (3.3) for variable $t$ as follows:

Notice that $x$ is also a parameter in Eq. $(\ref{xxt})$. We define $F_{\epsilon}:L^{\infty}(\mathbb{S}^{1})\longrightarrow L^{\infty}(\mathbb{S}^{1})$ as

Here $\hat{u}\left(x+{\epsilon/2},t,\epsilon\right)$ is obtained from reconstruction $L(x_{-p},\dots,x_{p})$ as in Eq. $(\ref{ffdse})$:

Since $f(\cdot)$ and reconstruction $L(x_{-p},\dots,x_{p})$ are Lipschitz continuous, so are $f(L(x_{-p},\dots,x_{p}))$. Since flux $F_{\epsilon}$, defined in Eq. $(\ref{fluxeps})$, is a combination of Lipschitz continuous (in a remarkably simple way), we get that $F_{\epsilon}$ is also a Lipschitz function. Thus, based on the classical theory of ODEs in Banach spaces, in the Lipschitz continuous case, there is a local solution to $t\in[0,\delta(\epsilon)]$ for some $\delta(\epsilon)$ that depends on $\epsilon$. For the global solution, since $f$ is bounded (because $H$ is assumed to be Lipschitz continuous), we can extend the solution to $\delta(\epsilon)\to\infty$. From assumption 1, Eq. $(\ref{h1})$, the Lipschitz constants of each $F_{\epsilon}$ can be chosen uniformly on bounded sets $L^{\infty}({\mathbb{S}^{1}})\times[0,\infty)$. To demonstrate the existence of the solution $(\ref{xxt})$ globally (in time), it suffices to prove that, for fixed $\epsilon$, there exists a $c_{\epsilon}(t)<\infty$ such that $||u(\cdot,t,\epsilon)||_{\infty}\leq c_{\epsilon}(t)<\infty$. Here $c_{\epsilon}$ is a continuous function on $[0,\infty)$, which is not uniformly continuous in $\epsilon$. We also have that $H(u)$ is a bounded function; hence, if $u\geq a>0$, then $f$ is also bounded. Let $M$ be defined as

We have that $\displaystyle{|\partial_{t}u(x,t,\epsilon)|\leq\frac{4}{\epsilon}||u(\cdot,t,\epsilon)||_{\infty}M}$ where $\displaystyle{||u(\cdot,t,\epsilon)||_{\infty}=\text{ess}\sup_{x\in{\mathbb{S}^{1}}}|u(x,t,\epsilon)|}$. Solving, we obtain:

From Grönwall formula, we obtain that:

Bound $(\ref{ggr})$ indicates the existence of a global solution to the ODE $(\ref{ODE})$ for each fixed $\epsilon$. However, note that the solutions to the system of ODEs are not uniformly continuous in $\epsilon$. To demonstrate that the solutions to ODEs provide a weak asymptotic solution for $(\ref{noLi})$, we will prove that the solution is $L^{1}$ bounded uniformly with respect to $\epsilon$. To do so, let $T>0$ for $t+dt\leq T$ and $dt>0$. From the Taylor expansion with remainder term (for fixed $\epsilon)$, it follows that $(\ref{ODE})$ can be written as

where $||r(\cdot,t,\epsilon,dt)||_{1}\to 0$ (or $||r(\cdot,t,\epsilon,dt)||_{\infty}\to 0$), uniformly in $t\in[0,T]$ and fixed $\epsilon$ (and not uniformly continuous in $\epsilon$), when $dt\to 0$. This behavior results from the continuous differentiability of map $t\longrightarrow u(\cdot,t,\epsilon)$, $[0,\infty)\longrightarrow L^{\infty}({\mathbb{S}^{1}})$ for fixed $\epsilon$. Since we are interested in obtaining the $L^{1}$ bound, we take the absolute value,

From the definition of $f^{+}$ and $f^{-}$ in Eq. (2.7), we get that if $(\ref{Cflcon})$ is satisfied, then $(\ref{ODE2n})$ is true. Notice that $(\ref{Cflcon})$ is always valid if the CFL condition $(\ref{CFL})$ is satisfied. This proves that the condition $(\ref{CFL})$ provides the method with stability because by integrating $(\ref{ODE2n})$ and the appropriate translations of $\pm\epsilon$, we obtain

Here the remainder value $r_{1}(t,\epsilon,dt)=\int_{\mathbb{S}}|r(\cdot,t,\epsilon,dt)|dx$ is bounded, and $r_{1}(t,\epsilon,dt)\to 0$ when $dt\to 0$, uniformly in $t\in[0,T]$ for each fixed $\epsilon$. Notice that, for each $T>0$ given, we can divide interval $[0,T]$ into $n$ subintervals $[jdt_{n},(j+1)dt_{n}]$, where $dt_{n}=\frac{T}{n}$ and $0\leq j\leq n-1$. Applying this in $(\ref{fdsee})$, we get

Applying recursively for all intervals, we obtain

Notice that

Thus, taking the limit $dt\longrightarrow 0$ and using $r_{1}(t,\epsilon,dt)\longrightarrow 0$, when $dt\longrightarrow 0$ we obtain: