Two algorithms for a fully coupled and consistently macroscopic PDE-ODE system modeling a moving bottleneck on a road

1.   Introduction

In this paper we focus on a system modeling a slow vehicle (e.g. large truck or bus) which moves in a car traffic flow. We assume that the slow vehicle's velocity is influenced by the car density and, at the same time, it causes a capacity drop in the road at the position where it is located. Therefore we use the terminology moving bottleneck.

Mathematically, the problem gives rise to a coupled PDE-ODE system. Car density $\rho$ solves the following conservation law

with $t \in [0, T]$, $\rho \in [0, \rho_{max}]$, $(t, x) \in \mathbb{R}^+\times\mathbb{R}$, $f:[0, \rho_{max}]\times\mathbb{R}\to\mathbb{R}$ the flux function, $\bar\rho$ the initial datum and $y(t)$ the position of the moving bottleneck, which is assumed to solve the following Cauchy problem

with $y_0$ the position at initial time $t = 0$ and $w\geq 0$ the velocity of the moving bottleneck. An example of $f$ and $w$ will be given in the next Section, see equations (2.1)–(2.2).

It is well known that the function $\rho(t, \cdot)$ may be discontinuous, then solutions to both (1.1) and (1.2) are to be intended in generalized sense. The case of multiple moving bottlenecks will be also considered.

Relevant literature. Let us mention a few results related to our work. Modelling of vehicular traffic with fluid dynamic approach started with the seminal works of Lighthill and Whitham ^[38] and Richards ^[41]. The attention on this subject continued in the last years, both for the proposal of more accurate models, see, e.g., ^[3,16,18], and for the treatment of complex networks, see, e.g., ^{[12,13,29,30]} and the book ^[24]. To deal with the case of different vehicles, such as motorcycles, cars, buses, multi-population models were also considered, see ^[4,43].

The problem of determining the trajectory of a single car given the density of cars was addressed in ^[19], and the relative numerics in ^[11], taking advantage of vector techniques. In this case the observed car is assumed not to influence significantly the traffic flow, thus the system is not fully coupled. Here we consider the more complicated situation in which the position of the single vehicle influences the whole traffic flow.

The coupled system (1.1)–(1.2) was proposed in ^[32]. More precisely, that paper introduced a notion of solution based on weak solutions for the PDE and Filippov solutions for the ODE, and proved existence of such a solution under suitable assumptions.

Recently, another fully coupled system for a moving bottleneck problem was studied in ^{[20,21,22,37]} from both the theoretical and numerical point of view (see also ^[42] for an extension to the second order ARZ model). Furthermore, in ^[40] an optimal control problem using the maximal speed of the coordinated vehicle as control variable is stated. It is worth to point out the difference with the model we consider here: the model in ^[20,21], following the lines of ^[39], assumes that the moving bottleneck reduces the maximal density attainable in the road, together with its maximal flux (i.e. the capacity of the road). Moreover, it assumes that the moving bottleneck moves at constant speed for low car density. Instead, in the model we consider here the moving bottleneck does not reduce the maximal density attainable in the road, while it only reduces the maximal flux. Moreover, its velocity is always dependent on the car density. These features make the model fully consistent with a pure macroscopic point of view.

Let us also recall that the problem of modeling bottlenecks, moving or not, was addressed in ^[1,17,23] and the case of multiple moving bottlenecks was already considered in ^[33].

Finally, the terminology "moving bottleneck" is known also in transport engineering literature: see, among others, ^{[26,34,36,39]}. In these papers the problem of moving bottlenecks, both with a given path ^[26,36,39] and with a mutual interaction with the surrounding flow ^[34], is studied, but without a detailed description of the resulting schemes and of their convergence.

Goal. We propose two different numerical algorithms, the first allowing theoretical investigations, while the second allowing an actual implementation on a computer.

The first algorithm, called WFT-ODE, combines the Wave Front Tracking (WFT) method ^[6,31] to solve (1.1) with an exact method to solve (1.2). The moving bottleneck position is traced taking into consideration interactions with (shock or rarefaction) waves. Focusing on bounded variation data, we establish a convergence result for the approximate solution towards a solution to (1.1)–(1.2). We are also able to provide a convergence rate for $y$ in terms of total variation of initial datum of density, namely $TV(\bar\rho)$. To achieve the estimates, we measure the distance among successive approximations by means of tangent vectors. The latter represents distances among discontinuities (shocks or rarefaction shocks) for car density, and the difference of positions for moving bottleneck.

The second algorithm, called GOF, is a fractional step method combining the classical Godunov method ^[28] for (1.1) with an exact method for (1.2). There are two reasons for using the Godunov method: the easy implementation and the importance of this method in the transportation literature ^[35,36].

The GOF method is then extended to the case of multiple bottlenecks. For that, we consider two models. The first allows slow vehicles to pass each other, while the second inhibits such a possibility. The latter is more challenging from theoretical and numerical point of view, thus we focus on this case. Both methods are useful to deal, for instance, with the case of many buses on the same urban route ^[25].

Paper organization. In Section 2 we provide basic definitions and preliminary results. The theoretical scheme WFT-ODE is introduced in Section 3, while convergence results are given in Section 4. The second scheme GOF is described in Section 5. Then we deal with multiple bottlenecks in Section 6. Finally, numerical simulations are reported in Section 7.

2.   Basic definitions and preliminary results

In this section we provide the definition of solution to the system (1.1)–(1.2) and we detail the choice of the flux $f$ and the velocity $w$. We also give a description of Wave Front Tracking (WFT) algorithm for the case of a flux function depending on the space variable, as it occurs for (1.1). Finally, a result on possible interactions of the moving bottleneck with fronts in the WFT solution is proved. The latter will be extensively used in the rest of the paper.

We denote by $\mathcal{BV}(\mathbb{R})$ the space of functions with bounded variation. For every $T>0$, the definition of solution to the system (1.1)–(1.2) reads as follows.

Definition 1. A vector $(\rho, y)$ is a solution to (1.1)–(1.2) in $[0, T]$ if $\rho(t, \cdot)\in\mathcal{BV}(\mathbb{R})$ for a.e. $t\in [0, T]$ which solves the PDE in the sense of distributions, that is

for any differentiable function $\phi(t, x)$ with compact support in $t\geq 0$. Moreover, the position of the moving bottleneck $y(t)$ solves the ODE in $[0, T]$ in the sense of Filippov, namely $y(t)$ is an absolutely continuous function such that $y(0) = y_0$ and

for a.e. $t\in [0, T]$, where the set $\mathcal{I}[a, b]$ is defined as the smallest interval containing $a$ and $b$.

We assume hereafter that the maximal density of cars is $\rho_{max}$ = 1. We denote by $w(\rho)$ the velocity of the moving bottleneck and we assume that

for some constant $w_{max}>0$. We denote by $v(\rho, x-y)$ the velocity of cars and we assume that

where $\varphi\in C^1(\mathbb{R}; (0, +\infty))$ and it is strictly decreasing in $(-\beta, 0]$ and strictly increasing in $[0, \beta)$ for some $\beta>0$. We also assume that there exist two positive constants ${\underline v }$, ${\overline v }$ such that

and that ${\underline v }$ is reached for $\zeta = 0$ and ${\overline v }$ is reached for any $\zeta\in\mathbb{R}\backslash[-\beta, \beta]$, see Figure 1(a).

Note that, in particular $\varphi'(\zeta)$ is compactly supported and, for any fixed density $\rho$, the minimal velocity of cars ${\underline v }$ is taken when they overtake the moving bottleneck, i.e. for $x = y$. We finally assume

The last assumption is crucial for the model: it means that cars can overtake the bottleneck, see Figure 1(b).

2.1.   Preliminaries for Wave-Front Tracking algorithms

As it is well known, a solution to the Cauchy problem of a conservation law can be constructed using the WFT algorithm. Roughly speaking, first a $\mathcal{BV}$ initial datum is approximated by a piecewise constant function via sampling (thus with a smaller $\mathcal{BV}$ norm). For small times, a piecewise constant weak solution is obtained by piecing together solutions to Riemann problems, where rarefactions are replaced by a fan of rarefaction shocks. Then, when waves interact, a new Riemann problem is generated and solved again approximating rarefactions by rarefaction shocks' fans and so on. For details we refer the reader to ^[6] for the general theory and to ^[24] for the case of networks.

In our case the flux does not depend only on $\rho$, but also on $x$ and $t$ (via $y(t)$), so the WFT algorithm needs to be modified. Moreover, we will evolve, at the same time, the position of the bottleneck $y(t)$ according to equation (1.2). In order to perform the construction we first need a preliminary result comparing the speed of the bottleneck with that of neighboring waves of the WFT approximation of $\rho$.

Lemma 1. Assume (2.1)–(2.4) hold true. Assume that the discontinuity $(\rho_l, \rho_r)$, located at $x$, is the closest to the bottleneck position $y$ and there exists

such that $\rho_l, \rho_r \geq\rho_{min}$. Define

then

We remark that

Proof. Recall that the discontinuity travels with speed:

We divide the proof in two cases:

1. The wave is a shock wave.

2. The wave is a rarefaction-shock wave.

1. It must be ${\rho _r}>{\rho _\ell }$. Using (2.8), to prove inequalities in (2.7) we consider the following:

and

Let us first prove (2.10). We notice that

Thus, (2.10) can be written as

and then as

Since from (2.3) and (2.4) we have $\varphi(\zeta)\geq{\underline v }>w_{max}$, the worst case is when ${\rho _\ell } = \rho_{min}$ and ${\rho _r} = \rho_{min}+\varepsilon$ ($\varepsilon>0$ arbitrarily small):

then it is sufficient to prove

If we define $\psi(z) = \frac{z-w_{max}}{2 z-w_{max}}$, with $z = \varphi(\zeta) \in [{\underline v }, {\overline v }]$, we find that $\psi'(z) >0$, and consequently:

Therefore, in the worst case (taking also $\varepsilon = 0$)

which is true for $\mu$ defined as in (2.6). Since ${\rho _r}>{\rho _\ell }$, (2.10) implies (2.9) and then we conclude the proof.

2. In the case of rarefaction-shock we have $\rho_r = \rho_\ell - \delta_\nu$, where $0 < \delta_\nu < \frac1{\nu}$ is a parameter that will be defined by the WFT algorithm. We only prove the inequality in (2.9), then (2.10) follows because ${\rho _\ell }>{\rho _r}$. The inequality in (2.9) can be rewritten as

Reasoning as above, we have

Now the worst case occurs when ${\rho _r} = \rho_{min}$, ${\rho _\ell } = \rho_{min}+\delta_\nu$:

i.e.

Then we conclude as in step 1.

The WFT algorithm will need to be further split because of source effect due to the $x$ dependence of the flux. Indeed, consider a general equation of the type

where we assume

Therefore, formally deriving in space the flux function $g$ depending both on the density of cars $\rho(t, x)$ and on the position of cars $x$, one gets:

so that equation (2.13) can be written as a conservation law with source term

In our case we have

and the source term reads

(see Figure 2). Note that the assumptions (H1) are satisfied and that the $x$ dependence in the flux is indeed effective only for $\zeta\in(-\beta, \beta)$, that is when cars are close to the bottleneck's position.

Equation (2.13) can be solved applying an operator splitting method in which one alternates between solving the homogeneous conservation law

and the ordinary differential equation

3.   The semi-discrete WFT-ODE scheme and its properties

We are now ready to define the WFT-ODE scheme. In this section we describe a semi-discrete scheme for the coupled dynamics (1.1)–(1.2). The idea is to combine the WFT method for the conservation laws with exact solution to the ODE of the moving bottleneck. Indeed, since the approximate solution given by WFT is piecewise constant (in time and space), we can get a simple ODE with piecewise constant (in time) right-hand side for the moving bottleneck.

Let us now describe in detail the WFT-ODE scheme. Fix $T>0$, then for every $\nu\in\mathbb{N}$ and $\Delta t>0$, we define an approximate solution $(\rho_\nu, y_\nu)$ on the time interval $[0, T]$ via the following steps.

Step 1. Let $\bar\rho\in\mathcal{BV}(\mathbb{R})$ the initial datum. We define a piecewise constant function $\bar\rho_\nu$ with $N$ discontinuities $s_1 < \cdots < s_N$, by approximating $\bar\rho$ as follows (Figure 3):

Notice that $TV(\bar \rho_\nu)\le TV(\bar \rho)$.

Step 2. At time $t = 0$, we solve all Riemann problems at each discontinuity point of $\bar\rho_{\nu}$ for the homogeneous equation (2.16) and approximate every rarefaction wave with a rarefaction fan, formed by rarefaction shocks (non-entropic shocks) of strength less than $\delta_\nu$ travelling with speed given by the Rankine-Hugoniot condition.

More precisely: Assume, for simplicity, that $y_0\in]s_j, s_{j+1}[$ for some $j$ (i.e. $y_0\not = s_i$ for every $i = 1, \ldots, N$.) For small times we can get an approximate solution by solving the equations:

for the discontinuities $s_i$ of $\rho_\nu$ and the equation:

for the bottleneck, which gives $y(t) = y_0+tw(\rho(t, s_j+))$.

If two wave fronts $s_i$ and $s_{i+1}$ interacts, then we solve the new Riemann problem. Because of Lemma 1, $s_j(t) < y(t)$ for every $t\in[0, \Delta t]$. If $y(\bar t) = s_{j+1}(\bar t)$ for $\bar t\in[0, \Delta t[$, then the bottleneck interacts at $\bar t$ with the wave $s_{j+1}$ and we simply change the equation for the bottleneck to:

and so on for all the possible interactions with $\rho_\nu$ waves up to time $\Delta t$.

Finally, possibly slightly modifying the wave speeds, we can assume that $y_\nu(\Delta t) \neq s_i(\Delta t)$ for any $j$.

Step 3. To take into account the source term $\rho_\nu(\Delta t, x)$ is updated in order to include the effect of the source term. This is done by means of one step of the explicit Euler method approximating the ODE (2.17) by:

After this modification the function $\rho_\nu^{source}$ is no longer piecewise constant and we need to apply a specific sampling to ensure accurate estimates of continuous dependence.

Step 4. Let $\bar y = y(\Delta t)$. We update $\rho_\nu^{source}$ by sampling only on the interval $[\bar y-\beta, \bar y+\beta]$ (where $\rho^{source}$ is not piecewise constant.) Let $x_i$, $i = 1, \ldots, M$, be the point of discontinuity of $\rho_\nu^{source}$ on $[\bar y-\beta, \bar y+\beta]$. First we insert the points $y_j = \bar y-\beta+j\, \frac{2\beta}{2^\nu}$, $j = 1, \ldots, 2^\nu$. Moreover, if there exists $i$ such that no $y_j$ is in the interval $]x_i, x_{i+1}[$, then we insert the additional point $y_{i} = \frac{x_i+x_{i+1}}{2}$. Again possibly slightly modifying the wave speeds we can assume that $x_i\not = y_j$ for every $i$ and $j$. We define $\rho^{sampled}$ by sampling $\rho^{source}$ on subintervals generated by the partition: $P = \{x_i\}\cup\{y_j\}\cup \{y_i\}$ with the following rules:

-  On intervals $[x_i, y_j]$ we define $\rho_\nu^{updated} = \rho_\nu^{source}(x_i+)$.

-  On intervals $[y_j, x_{i+1}]$ we define $\rho_\nu^{updated} = \rho_\nu^{source}(x_{i+1}-)$.

-  On all other intervals $[a, b]$ we define $\rho_\nu^{updated} = \rho_\nu^{source}(a+)$.

The rationale behind this definition is the following. To obtain neat estimates of increase in jumps size of shocks and rarefaction shocks we need always to sample to the left and right of the discontinuities which were generated by the WFT step. In the sequel, we shall refer to these new generated (at positive time) waves as $s_{b}$ waves, being generated by the source term due to the bottleneck.

We can restart the procedure at step 2 in the interval $[\Delta t, 2\Delta t]$ and so on.

It is important to note that, while the sampling procedure does not increase the total variation, the updating step (3.1) can do it. For the first time step, the increment can be estimated easily as follows. Let us denote by $\omega(x; \rho_\nu, \Delta t)$ the contribution given by the source term at time $\Delta t$, i.e.

Recalling the assumption (H1), we have

Then, at every time step the additional total variation is bounded by $\max\limits_{\rho}\|\omega'(\cdot; \rho, \Delta t)\|_{L^1}$, and at final time $T$, namely after $\frac{T}{\Delta t}$ time steps, the total increment is bounded by $T \max\limits_{\rho}\|\partial_{xx}g(\rho, \cdot) \|_{L^1}$. More precisely we have the following:

Lemma 2. The increase of total variation in $\rho_\nu$ at Step 3 is bounded by $K_\varphi\ \Delta t$ where

Bounds on the $\mathcal{BV}$ norm and on number of waves and interactions are directly obtained as in the scalar case, when the flux is not depending on $x$, provided this dependence does not create resonance effects among waves. Indeed, we can rewrite a scalar conservation law with a $x$-dependent flux $g(\rho, x)$ as a $2\times2$ system as follows:

The above system is strictly hyperbolic, and therefore standard WFT analysis may be carried out, provided $\partial_\rho g(\rho, z)\neq 0$. In our case, $g(\rho, x; y(t)) = \rho(1-\rho)\varphi(x-y(t))$ and the above condition corresponds to the requirement that the bottleneck's position $y(t)$, where the $x$-dependence is effective as noticed before, travels with a speed uniformly different with the one of the closest $\rho$ waves. This is precisely the statement of Lemma 1.

The construction thus relies on Lemma 1, which in turn is applicable if $\rho_\nu>\frac{{\overline v }-w_{max}}{2 {\overline v }-w_{max}}$, see (2.5). Due to the source term $-\partial_x g$ this is not guaranteed for all times even if the initial data satisfies (2.5). However we have the following:

Lemma 3. Consider an initial datum $\bar{\rho}\in\mathcal{BV}(\mathbb{R})$ and assume $\eta = \min_x \frac{1}{\bar{\rho}(x)}\ \frac{{\overline v }-w_{max}}{2 {\overline v }-w_{max}} < 1$. Then for all positive times $T$ such that:

the condition (2.5) of Lemma 1 is verified for $\Delta t$ sufficiently small

Notice that for $\bar{v}\to w_{max}$ we have $\eta\to 0$, thus $T\to \infty$. In other words, if the maximal speed of cars is close the that of the moving bottleneck then the WFT-ODE scheme is defined for large times.

Proof. Notice that solutions to the ODE linked to the source term (2.17) satisfy $\partial_t \rho = \rho(1-\rho)\partial_x\varphi\geq -\|\varphi'\|_{\infty} \rho$, thus the solution satisfies $\rho(t)\geq e^{-t\|\varphi'\|_{\infty}} \rho(0)$. This gives the estimate (3.3). Since Step 3 is an approximation of the source term, the estimate remains valid for $\Delta t$ sufficiently small.

4.   Convergence of the WFT-ODE scheme

In this section we prove the convergence of the WFT-ODE scheme. The BV estimates on the WFT scheme allows to pass to the limit by standard arguments. On the other side, the convergence of the the bottleneck position requires additional work, in the same spirit of ^[11]. However, our case is more complex due to the complete coupling of the PDE for car density and the ODE for the bottleneck position. As in ^[11], to achieve this goal we use the technique of generalized tangent vectors ^[7,8], to the WFT approximate solutions $\rho_\nu$ and to the bottleneck position $y_\nu$.

More precisely, after introducing the tangent vector technique, we consider $\rho_{\nu+1}$ as obtained from $\rho_\nu$ by shifts of discontinuities, which generate generalized tangent vectors. Then we need to estimate:

1) The increase in the norm of tangent vector to $\rho_\nu$ due to the WFT algorithm.

2) The increase in the norm of tangent vector to $\rho_\nu$ due to the source term.

3) The increase in the norm of tangent vector to $y_\nu$ due to interactions with waves of $\rho_\nu$.

4.1.   Generalized tangent vectors

The technique is based on the idea of considering $L^1$ as a Finsler manifold with generalized tangent vectors with appropriate norm. We first consider the subspace of piecewise constant functions and "generalized tangent vectors" consisting of two components $((v^\theta, \xi^\theta), \eta^\theta)$, where $\xi^\theta$ describes the infinitesimal displacement of discontinuities and the scalar $\eta^\theta$ is the {infinitesimal} shift of the car trajectory. Let us take a family of piecewise constant functions $\{\rho ^{\theta }\}_{\theta \in [0, 1]}$, each of which has the same number of jumps, say at the points $s_{1}^{\theta} < ... < s_{N}^{\theta }$. Let us define the function

and also the quantities

and

Note that $v^\theta$ is not defined if $x = s^\theta_k$, $k = 1, \ldots, N$. Indeed, the contribution of the jumps to the tangent vector is given by $\{\xi^\theta_k\}_k$, which take into account the presence of the discontinuities. Moreover, $v^\theta$ solve the linearized equation along the trajectory, see ^[6] while the shifts change only at interactions times or, as detailed in next remark, due to Step 3 in the WFT-ODE scheme.

Remark 1. It is worth observing that, as a consequence of the shift $\eta ^{\theta}$ of the moving bottleneck's position, the waves $s_b$ generated in Step 4 of the scheme will be all shifted by the same quantity.

Then we say that the path $\gamma:\theta \rightarrow (\rho ^{\theta }, y^\theta)$ admits tangent vectors $((v^{\theta}, \xi ^{\theta }), \eta^{\theta })\in T_{\rho^\theta}\dot{ = } L^{1}(\mathbb{R}, \mathbb{R}) \times \mathbb{R}^{N}\times \mathbb{R}$. Note that in general such path is not differentiable w.r.t. the usual differential structure of $L^{1}$. One can compute the $L^{1}$-length of the path $\gamma$ in the following way:

According to (4.1), in order to compute the $L^{1}$-length of a path $\gamma$, we integrate the norm of its tangent vector which is defined as follows:

where $\Delta \rho^\theta_{k }$ is the jump across the discontinuity $s^\theta_{k }$. Notice that this is not the usual $L^1$ length of a path since there is the additional component $\eta$.

Let us introduce the following definition.

Definition 2. We say that a continuous map $\gamma :\theta \rightarrow (\rho^{\theta }, y^\theta)\dot{ = }\gamma (\theta)$ from $[0, 1]$ into $L_{loc}^{1}\times\mathbb{R}$ is a regular path if the following holds. All functions $\rho ^{\theta }$ are piecewise constant, with the same number of jumps, say at $s_{1}^{\theta } < ... < s_{N}^{\theta }$ and coincide outside some fixed interval $\left] -M, M\right[$. Moreover, $\gamma$ admits a generalized tangent vector $D\gamma (\theta) = ((v^{\theta }, \xi ^{\theta }), \eta^\theta)\in T_{\gamma (\theta)} = L^{1}(\mathbb{ R}; \mathbb{R})\times\mathbb{R}^{N}\times \mathbb{R}$, continuously depending on $\theta$.

4.2.   WFT approximations with shift evolution

Now, given two successive approximations of the initial data $\bar\rho$, namely $\bar\rho_\nu$ and $\bar\rho_{\nu+1}$, obtained by sampling, respectively, at lattice points of mesh $2^{-\nu}$ and $2^{-(\nu+1)}$, we can connect them by a regular path as follows (Figure 4):

Notice that such path satisfies $\rho^0 = \bar\rho_{\nu+1}$ and $\rho^1 = \bar\rho_{\nu}$. Moreover, $\rho^\theta$ is differentiable everywhere (possibly) except for $\theta = 0, 1$, where the number of jumps changes. Finally we have:

Lemma 4. The initial data $\bar{\rho}_\nu$ and $\bar{\rho}_{\nu+1}$ can be connected by a regular path with shifts whose norm is bounded by $2^{-(\nu+1)}$.

Since $\rho^\theta$ is piecewise constant, $v^\theta(x)$ is equal to 0 (where it is defined), then $\|v^\theta\|_{L^1} = 0$. At time $t = 0$, we also have $\eta^\theta = 0$ since $y_\nu(0) = y_{\nu+1}(0)$. By (4.1), we have

Now we consider the evolution in time of $\gamma$, denoted by $\gamma_t$, in the sense that we perform the wave front tracking for $\bar\rho_\nu$ and $\bar\rho_{\nu+1}$. It is easy to prove that $\gamma_t$ is still regular (some more point of not differentiability will arise because of interactions order). We aim at proving that

for an appropriate constant $K>0$. Then uniqueness and Lipschitz continuous dependence of solutions to Cauchy problems is straightforwardly achieved passing to the limit on the WFT-ODE approximate solutions.

4.3.   Estimates on tangent vectors to $\rho_\nu$ due to WFT

To estimate the shifts of waves at interactions times between waves, we can use the estimate [24,Lemma 2.7.2] (originally proved in ^[5]).

Lemma 5. Consider two waves with speeds $\lambda_1$ and $\lambda_2$ respectively, that interact together producing a wave with speed $\lambda_3$. If the first wave is shifted by (the tangent vector of the shift) $\xi_1$ and the second wave by (the tangent vector of the shift) $\xi_2$, then the (tangent vector of the) shift of the resulting wave at time of interaction is given by

Moreover we have

where $\Delta \rho_i$ are the signed strengths of the corresponding waves.

Proof. The interacting waves satisfy the ODE:

From the regularity of $\varphi$ we have that $x_i$ are smooth functions of time, thus we can linearize them around the interaction time and apply Lemma 2.7.2 of ^[24].

The wave shifts change also when the waves are located within the influence zone of the moving bottleneck, more precisely we have:

Lemma 6. Consider a wave $(\rho_l, \rho_r)$ with shift $\xi$ which is in the influence zone of the moving bottleneck on the time interval $[t_1, t_2]$ and does not interact with other waves of the bottleneck. Then we have:

Proof. Let $x(t; z)$, $t \in [t_1, t_2]$ be the unique solution to the ODE:

Let $\bar{x}$ be the position of the wave at time $t_1$ and consider the solution to shifted initial position $x(t; \bar{x}+\xi)$, $t \in [t_1, t_2]$. Note that, since $\lambda = \frac{\varphi(\zeta)\Delta (\rho(1-\rho))}{\Delta \rho}$, the solution $x(t; \cdot)$ satisfies

uniformly in $t \in [t_1, t_2]$, with $C = \|\varphi'\|_{\infty} = \max_{\mathbb{R}}|\varphi'|$ the uniform Lipschitz constant for $\lambda$, and then we also obtain the estimate

4.4.   Estimates on tangent vectors to $\rho_\nu$ due to the source

As specified above, our WFT algorithm includes a sampling procedure on the interval $[\bar y-\beta, \bar y+\beta]$ where $\rho^{source}$ is not piecewise constant due to the action of the source term; see (3.1). The so-defined piecewise constant approximation $\rho^{updated}$ has new tangent vectors which are estimated in the following lemma.

Lemma 7. Let $(\rho_l, \rho_r)$ be a discontinuity of $\rho^{source}$ with shift $\xi^-$ located inside the region of influence of the bottleneck, and let us denote with $|\Delta \rho^-|$ its strength: $|\Delta \rho^-| = | \rho_l - \rho_r|$. Then in $\rho^{updated}$ this discontinuity has unmodified shift and a new strength $|\Delta \rho^+|$ which verifies the following estimate

As a consequence, the tangent vector is estimated as follows:

where $\xi^+ = \xi^-$ is the (unmodified) shift of the wave in $\rho^{updated}$.

Proof. Due to Step 4 in our WFT algorithm (cfr. Section 3), a wave $(\rho_l, \rho_r)$ of $\rho^{source}$ located at a point in the region of influence of the bottleneck will be located in $\rho^{updated}$ at the same point, that is, $\xi^- = \xi^+$, indicating by ${}^\pm$ the values before and after the source action and the sampling. However, the jump will be changed by this procedure, namely $(\rho_l-\Delta t\ \partial_x\varphi\ \rho_l(1-\rho_l), \rho_r-\Delta t\ \partial_x\varphi\ \rho_r(1-\rho_r))$. Therefore

and in addition

4.5.   Estimates on tangent vectors to $y_\nu$

The tangent vector to $y_\nu$, i.e. $\eta$, changes only at interaction times of the moving bottleneck with waves. Both the moving bottleneck and the waves satisfy ODEs with smooth right-hand-side, thus by linearization we find the same interaction estimate as for the case treated in ^[11].

Proposition 1. Let $t^*$ be an interaction time between the moving bottleneck and a wave $({\rho _\ell }, {\rho _r})$ and indicate by $\eta_\pm$ the moving bottleneck shift after and before the interaction respectively. Then we have:

where $\mu$ is given by Lemma 1.

4.6.   Convergence of WFT-ODE

To prove convergence of the WFT-ODE scheme, we estimate the increase of the tangent vectors $\eta$ to the moving bottleneck position $y_\nu$. We fix $\Delta t$ and $\nu$ and start introducing some notation.

Definition 3. We denote by $(TV)_j^{p}$ the total variation of waves of $\rho_\nu(j\Delta t)$ which are to the right of the moving bottleneck (thus can potentially interact with it), and by $(TV)_j^{i}$ the total variation of waves of $\rho_\nu$ which interact with the moving bottleneck in the time interval $[j\Delta t, (j+1) \Delta t]$. Similarly we denote by $V_j^{p}$ the sum of the norms of tangent vectors to waves of $\rho_\nu(j\Delta t)$ which are to the right of the moving bottleneck (thus can potentially interact with it), and by $V_j^{i}$ the norms of tangent vectors to waves of $\rho_\nu$ which interact with the moving bottleneck in the time interval $[j\Delta t, (j+1) \Delta t]$. We denote by $\eta_j$ the tangent vector to $y_\nu$ at time $j\Delta t$.

With these notations we have the following estimates. From Lemma 2, we get:

From Lemmas 2, 5, 6 and 7, and Remark 1, we get:

More precisely, Lemma 2 guarantees that the total variation of new waves due to Step 3 is bounded by $K_\varphi \Delta t$, thus the new contribution to the tangent vector norm is bounded by $K_\varphi \Delta t\, \eta_{j+1}$ by Remark 1; Lemma 5 guarantees that the norm of tangent vector does not increase for waves interactions; Lemma 6 shows that the magnification due to permanence in the moving bottleneck influence zone is bounded by $e^{\|\varphi'\|_{\infty} \Delta t}$; Lemma 7 provides the bound $(1+\|\varphi'\|_{\infty} \Delta t)$ on magnification due to the source term. Finally, from Proposition 1, $\eta$ may increase by interacting with a wave $({\rho _\ell }, {\rho _r})$ having shift $\xi$ with multiplicative factor at most $(1+ \frac{w_{max}}{\mu} |{\rho _\ell }-{\rho _r}|)$ for any rarefaction wave and additive factor $\frac{w_{max}}{\mu} |\xi|\ |{\rho _\ell }-{\rho _r}|$ for any wave. The worst case scenario happens when first all shocks interact and then all rarefactions, this is bounded by the estimate:

We can rewrite the estimates (4.14)–(4.16) as:

where $K = \max \{K_\varphi, 2\|\varphi'\|_{\infty}, \frac{w_{max}}{\mu} \}$. We now have the following:

Lemma 8. Consider sequences $(TV)_j^p$, $(TV)_j^i$, $V_j^p$, $V_j^i$ and $\eta_j$ satisfying (4.17). Let $\mu_j$ be the solution to the system:

with initial datum $\mu_1 = K e^{KT}\, V_0$ and $V_1 = K\Delta t\, \mu_1$. Then we have $\eta_j\leq e^{K\, TV(\bar{\rho})+KT} \mu_j$.

Proof. First we have $\sum_i (TV)_j^i\leq TV(\bar{\rho})+KT$. Notice that $\eta_j$ is always increasing so the worst case for the multiplicative terms $e^{K (TV)_j^i}$ in (4.17) (third equation) is when the $(TV)_j^i = 0$ for $j = 1, \ldots, \frac{T}{\Delta t}-1$ and $(TV)_{\frac{T}{\Delta t}}^i = TV(\bar{\rho})+KT$. Also notice that the increase of $\eta_{j+1}$ due to the $V_j^i$ term is maximized when $V_j^i = V_j^p$. On the other side this may not achieve the maximal increase in $V^p_{j+1}$ because of the multiplicative term $(1+K\Delta t)$. However, such maximal increase of $V^p_{j+1}$ on the interval $[0, T]$ (thus after all time steps) is bounded by $e^{KT}$, which is the term appearing on the right-hand side of (4.18) (second equation). Therefore, with this corrected multiplicative term, the increase is maximized. Finally, $\eta_0 = 0$ thus the initial data for (4.18) are given by $\mu_1 = K e^{KT}\, V_0$ and $V_1 = K\Delta t\, \mu_1$.

We are now ready to estimate $\eta_j$. First notice that we have:

thus:

and finally:

Thus we obtain the following:

Proposition 2. Consider initial conditions $\bar\rho \in \mathcal{BV}(\mathbb{R})$, $y(0)\in\mathbb{R}$ and time horizon $T>0$. Using the notation of Lemma 3, assume that $\eta < 1$ and $T$ satisfy (3.3). Let $(\rho_\nu, y_\nu)$ be the approximate solutions computed by the WFT-ODE scheme, then:

where $K_1$ depends only on the total variation of $\bar\rho$, $T$, $w_{max}$, $\|\varphi\|_{\infty}$, $\|\varphi'\|_{\infty}$, $\mu$ defined in (2.6) and $K_\varphi$ defined in (3.2).

Proof. From (4.19), we have that $\eta_\nu$ can be estimated in terms of the total variation of $\bar{\rho}$, $T$, the constant $K$ and $V_0$. But $V_0$ is estimated by the total variation of $\bar{\rho}$ and the initial shifts, thus we conclude by Lemma 4.

5.   Godunov-ODE-FS scheme

Here we introduce a numerical scheme called Godunov-ODE-FS (GOF), which is based on fractional step method combining Godunov scheme for the PDE and exact solution for the ODE. The dynamics of (1.1) and (1.2) are thus solved separately at each iteration. The use of Godunov scheme is motivated by its easy implementation and its connection with the modelling of vehicular traffic problems, see ^[24,25].

5.1.   Godunov scheme for the PDE

Following the ideas described in Section 2.1 for the WFT algorithm, we report in the following the modified Godunov scheme for a conservation law with source term

(see (2.15)). We first introduce a numerical grid, denoting by $\Delta x$ the space mesh size and by $\Delta t$ the time mesh size. Moreover, we denote by $(t_l, x_m) = (l \Delta t, m \Delta x)$ the grid points for $l = 0, 1, \dots, L$, $m = 0, 1, \dots, M$, where $L$ and $M$ are, respectively, the number of time and space nodes of the grid, and by $C^l_m$ the discretization cell $(t_l, t_{l+1}) \times (x_{m-1}, x_m)$. For a function $u$ defined on the grid we write $u_m^l = u(t_l, x_m)$.

Let us denote by

the self-similar entropy solution of the unique Riemann problem defined on $C^l_m$ (with discontinuity point $x_{m-\frac{1}{2}}$) and let us define the numerical flux $F$ as $F(\rho_\ell, \rho_r, t, x) = g(W(0;\rho_\ell, \rho_r), t, x)$.

First, we replace the initial datum $\bar\rho(x)$ by a piecewise constant approximation,

Then, we alternate a single step of the classical Godunov scheme

with a single step of the Euler scheme

5.2.   GOF scheme

Before describing the GOF scheme, we point out that shock waves solutions to (1.1) have velocities depending on the moving bottleneck position $y$. The modified Godunov scheme can be thus used, provided that, in each discretization cell, the shock speed $\lambda({\rho _r}, {\rho _\ell }, \zeta)$ defined in (2.8) does not change sign as a function of $\zeta$. This is established by next Lemma.

Lemma 9. Let $\bar f (\rho) = \rho(1-\rho)$ and $\rho_\ell$ (resp., $\rho_r$) the left (resp., the right) state of a discontinuity. Set $\bar \lambda = \frac{\bar f(\rho_r) - \bar f(\rho_\ell)}{\rho_r - \rho_\ell}$. Then,

Proof. Since we have:

the result (5.2) is easily achieved noticing that $\varphi(\zeta) > 0$.

The situation described by Lemma (5.2) is depicted in Figure 5.

The discretization grid is defined as for the Godunov scheme and the value of density on grid points is called $\rho^l_m$. The moving bottleneck position will be defined for times $t_l$ and denoted by $y^l$.

Step 0. We approximate the initial datum $\bar\rho$ as for the Godunov scheme and set $y^0 = y(0)$.

Step 1. We assume that $y^l$ and $\rho^l_m$ were defined by previous step.

We keep fixed the position of the moving bottleneck, thus set $y(s) = y^l$ for $s \in [t_l, t_{l+1}]$. Define $W_m = W_m(t, x)$ to be the self similar solution to the Riemann problem defined in $x_{m-\frac{1}{2}}$ (for the bottleneck position fixed). This solution is formed by a wave, whose velocity varies in time but mantains its sign constant as proved by Lemma 9. We can use Gauss-Green formula as for Godunov scheme and define

Step 2. There exists $m$ such that $y^l\in [x_{m-1-\frac12}, x_{m-\frac12})$, see Figure 6. Using the computed density $\rho^{l+1}$, the initial velocity of the moving bottleneck is $w(\rho^{l+1}_{m-1})$ up to the interaction time with the line $x = x_{m-\frac12}$ then it travels with velocity $w(\rho^{l+1}_{m})$. More precisely, we compute the interaction time as

Now, if $\Delta t_{in}\leq\Delta t$ we set

otherwise we set

Remark 2. Notice that for the flux (2.2) the CFL condition reads

where $\zeta^l_m = x_m - y^l$ and $I(\rho^l_{m-1/2}, \rho^l_{m+1/2})$ is the interval with endpoints $\rho^l_{m-1/2}$ and $\rho^l_{m+1/2}$. In particular, CFL is always satisfied if $\Delta t \leq \frac{1}{2\bar v}\Delta x$.

6.   Multiple bottlenecks

In this section we introduce a model describing the evolution of car traffic flow along a road in presence of $P>1$ moving bottlenecks. In this case, we face the event that two slower vehicles can occupy the same position, which could lead to a non realistic dropping of flux capacity. To solve this problem, we introduce two modelling choices:

(A) Two moving bottlenecks can overtake each other, thus we re-define the flux function.

(B) Two moving bottlenecks cannot overtake each other. To achieve that, we re-define the velocity function of moving bottlenecks.

6.1.   Model A

Let us consider the following PDE-ODE coupled system

where $y_i = y_i(t)$ is the position of the $i$-th moving bottleneck, and $\rho = \rho(t, x)$ is the density on the road. The flux function $f$ is given by

where ${\bf{y}} = (y_1, \ldots, y_P)$ and the smooth function $v(x, \rho, {\bf{y}})$ is defined as

with the average velocity of cars $\hat v(\rho) = 1-\rho$. Each function $\varphi_i(\zeta)$ takes into account the flux capacity drop due to the presence of the $i$-th bottleneck and we assume:

Speeds of moving bottlenecks are defined as $w_i(\rho) = w_{max}^i(1-\rho)$ with ${\underline v }_{i}> w_{max}^i$.

Now we admit the presence of many slow-moving vehicles in the car traffic flow, with possibly different characteristics among each other. Therefore, in general, functions $w_{i}(\cdot)$ and $\varphi_{i}(\cdot)$ are different and $P$ can be arbitrarily large.

In view of our definition of the function $\Phi$ in the flux, when two slow vehicles are far enough, we go back to the case of a single moving bottleneck. On the other side, if two slow vehicles occupy the same position, the flux capacity drop is the largest among the two.

6.2.   Model B

We now introduce the alternative modelling choice in (B), which avoid overtaking between slow vehicles, again described by a coupled PDE-ODE system

with $y_{1, 0} < \ldots < y_{P, 0}$. The flux function is defined as in (6.1)–(6.2), but this time the function $\Phi$ is given by

We define differently also the velocities of moving bottlenecks. More precisely, to avoid overtaking among them, and indeed to avoid the superposition of the zones where functions $\varphi_i$ take values less than ${\overline v }$, we suitable define the velocity functions $w_i(\rho, y_i, y_{i+1})$. First let $\omega_i = \omega^i_{max}(1-\rho)$, with $0 < \omega^i_{max} < {\underline v }_i$, and then back recursively on $i$ set

and complete the definition for $y_{i+1} - y_{i}\in [\beta_{i+1} + \beta_{i}, 2(\beta_{i+1} + \beta_{i})]$ by smooth monotone interpolation. This choice introduces a follow-the-leader flavour in the microscopic dynamic and avoids the modelling problems described above. A bottleneck at position $y_{i}$ moves with velocity $\omega_{i}(\rho)$, depending only on the density of cars along the road, provided the distance from the vehicle ahead in position $y_{i+1}$ is sufficiently large, namely, larger than $2(\beta_{i+1} + \beta_{i})$. If such distance decreases, then $y_i$ starts to decelerate and eventually move with the same velocity as the vehicle ahead when $y_{i+1} - y_{i} = \beta_{i+1}+\beta_{i}$.

6.3.   GOF algorithm for multiple bottlenecks

We introduce a Godunov-ODE-FS algorithm for model (B), while model (A) can be treated in a manner completely similar to the case of a single bottleneck. Notice that for model (B) slow moving vehicles can influence each other, thus the numerics is more involved.

The discretization grid is defined as for the Godunov scheme and the value of density on grid points is called $\rho^l_m$. The moving bottlenecks' positions will be defined for times $t_l$ and denoted by $y_i^l$, $i = 1, \ldots, P$. For consistency we suppose

Notice that the $P$-th bottleneck plays the role of the leader and it is not influenced by positions of other slow vehicles. Thus its trajectory can be traced as in Section 5.2. Furthermore, because of definition of $w_i$, the $i$-th bottleneck's trajectory is only influenced by the $(i+1)$-th bottleneck trajectory. More precisely, in the scheme we let the $i$-th bottleneck proceed with velocity $\omega_i$ as long as the distance with the $(i+1)$-th bottleneck is larger than $\beta_{i+1}+\beta_{i}$, otherwise we let it proceed with the velocity of the $(i+1)$-th bottleneck.

Let us first note that, as for the case of a single bottleneck, the velocity of any wave inside each cell change in time, but its sign does not. Indeed, let $\lambda(\rho_r, \rho_\ell, x, {\bf{y}})$ be the speed of the wave, then, using notation of Lemma 9, we have

and, since $\Phi(x, {\bf{y}})>0$, it follows that

We are now ready to introduce the GOF algorithm.

Step 0. We approximate the initial datum $\bar\rho$ as for the Godunov scheme and set $y^0_i = y_i(0)$.

Step 1. We assume that $y_i^l$, for $i = 1, \ldots, P$, and $\rho^l_m$ were defined by previous step.

We keep fixed the positions of moving bottlenecks, thus set $y_i(s) = y_i^l$ for $s \in [t_l, t_{l+1}]$. Let $W_m = W_m(t, x)$ be the self similar solution to the Riemann problem defined in $x_{m-\frac{1}{2}}$ (for bottlenecks' positions fixed.) This solution is formed by waves, whose velocity sign is constant. Setting ${\bf{y}}^l = {\bf{y}}(t_l)$ and using the Gauss-Green formula, the scheme is expressed in the integral formulation as

Step 2. We compute $y^{l+1}_P$ using the density $\rho^{l+1}$ obtained at Step 1, as for the case of a single bottleneck. More precisely, there exists $m$ such that $y^l_P\in [x_{m-1-\frac12}, x_{m-\frac12})$; $y_P$ moves with velocity $\omega_P(\rho^{l+1}_{m-1})$ until the interaction time $\Delta t_{in}$ with the line $x = x_{m-\frac12}$, provided that $\Delta t_{in}\leq\Delta t$. After that time, $y_P$ moves with speed $\omega_P(\rho^{l+1}_{m})$. Finally, if $\Delta t\leq\Delta t_{in}$ then we set

otherwise

Step 3.(Figure 7) We compute $y^{l+1}_i$ by backward recursion on $i = P-1, \ldots, 1$. More precisely, we define a trajectory $y_i(t)$ for every $t\in [t_l, t_{l+1}]$ and set $y^{l+1}_i = y_i(t_{l+1})$.

First for $i = P$, if $\Delta t\leq\Delta t_{in}$ we set

otherwise

Now, fix $i$ and assume to have defined $y_j(t)$ for $j\geq i+1$ and $t\in [t_l, t_{l+1}]$. Let $m$ be such that $y^l_i\in [x_{m-1-\frac12}, x_{m-\frac12})$ and define $\Delta t_{in}$ the time at which the line $y^l_i+(t-t_l)\omega_i(\rho^{l+1}_{m-1})$ intersects the vertical line $x = x_{m-\frac12}$. If $\Delta t\leq\Delta t_{in}$ we set

otherwise

Now, if there exists a time $\Delta t_{b} < \Delta t$ such that $y_{i+1}(t_l+\Delta t_{b})-\tilde{y}_{i}(t_l+\Delta t_{b}) = \beta_{i+1}+\beta_i$, then we set

otherwise we simply set $y_i(t) = \tilde{y}_i(t)$ for every $t\in [t_l, t_{l+1}]$.

Remark 3. Notice that for the dynamics (6.3) the CFL condition reads

where $I(\rho^l_{m-1/2}, \rho^l_{m+1/2})$ is the interval with endpoints $\rho^l_{m-1/2}$ and $\rho^l_{m+1/2}$. In particular, CFL is always satisfied if $\Delta t \leq \frac{1}{2\max_i{\overline v }_i}\Delta x$.

7.   Numerical tests

In this section we show some simulations obtained by models (1.1)–(1.2) and (6.3). In particular, we focus on interactions of bottlenecks with different type of waves, namely shocks and rarefactions. For both models we implement the corresponding GOF scheme.

7.1.   Model (1.1)–(1.2)

To compute the approximate solution to (1.1)–(1.2) we use the numerical algorithm presented in Section 5.2. We set the space and the time mesh size of the grid to be $\Delta x = 0.02$ and $\Delta t = 0.01$. Setting parameters, we assume $w(\rho) = w_{max}(1-\rho)$ with $w_{max} = 0.4$ and

where ${\overline v } = 1$, ${\underline v } = 0.6$ and $\beta = 0.1$.

In Figure 8 a moving bottleneck starting from $y_0 = 0.5$ interacts with a rarefaction wave centred in $x = 1.4$, generated by the Riemann datum $\rho_\ell = 0.9$ for $x\leq1.4$ and $\rho_r = 0.45$ for $x>1.4$. The vehicle is initially moving with low speed ($w = 0.04$) due to the high density until it interacts with the first characteristic of the rarefaction; after that it accelerates since the density is decreasing.

In Figure 9 a bottleneck starting from $y_0 = 0.5$ is interacting with a backward moving shock starting from $x = 1.4$. The latter is generated by the Riemann datum $\rho_\ell = 0.3$ for $x\leq1.4$ and $\rho_r = 0.9$ for $x>1.4$, so that the speed of the shock, before it interacts with the bottleneck, is $\lambda = -0.2$. On a following time interval, during the interaction with the bottleneck (from around $t = 1$ until around $t = 2.5$ in Figure 9), the speed of the shock almost vanishes, because of the capacity drop caused by the bottleneck.

Figure 10 shows a shock and a rarefaction which interact with each other and with the bottleneck. The rarefaction wave is centred in $x = 0.6$ while the shock is starting from $x = 1.6$.

7.2.   Model (6.3)

Simulations for the system in (6.3) are provided by means of the numerical algorithm defined in Section 6.3. Interactions of three bottlenecks and different type of waves are considered. We set the space and time mesh size to be $\Delta x = 0.02$ and $\Delta t = 0.01$. Velocity functions for slow moving vehicles are $\omega_i(\rho) = \omega_{max}^i(1-\rho)$, with $\omega_{max}^1 = 0.49$ and $\omega_{max}^i = 0.4$ $i = 2, 3$. We set again

where, as before, ${\overline v } = 1$. In Figure 11, we set ${\underline v }_i = 0.5$ and $\beta_i = 0.25$ for $i = 1, 2, 3$. The initial datum for the density is piecewise constant with single discontinuity at $x = 2.5$ and values $\rho_\ell = 0.9$ and $\rho_r = 0.1$. This gives rise to a rarefaction wave. Each slow vehicle in the road is interacting with the rarefaction wave, but in different ways. The second vehicle, starting from $y_2(0) = 1.5$, is accelerating, but not as much as the first one, starting from $y_3(0) = 2$. This is due to differences in the rate of density decrease. The last one, starting from $y_1(0) = 1$, has maximal velocity $w_{max}^1$ higher than the others. However, it is not capable of strong acceleration, since its initial distance w.r.t. the second is equal to $\beta_1+\beta_2 = 0.5$, thus its speed is bounded by that of the second vehicle.

In Figure 12, we simulate an homogeneous case, i.e. all bottlenecks have the same characteristics. We set ${\underline v }_i = 0.5$ and $\beta_i = 0.25$ for $i = 1, 2, 3$. The initial datum is piecewise constant with $\rho_\ell = 0.85$ for any $x\leq3.5$ and $\rho_r = 0.95$ for any $x>3.5$, thus generating a shock with negative speed $\lambda = -0.8$. Initial data for bottlenecks are $y_{1, 0} = 1$, $y_{2, 0} = 2$ and $y_{3, 0} = 3$. Starting from $t = 0$, due to capacity dropping, the density behind each slow vehicle increases giving rise to a queue of cars. This phenomenon generates further waves travelling and interacting with bottlenecks as well as the main shock. When a bottleneck cross the shock, it enters a region with higher density, so that its velocity decreases. On the other hand, during the interaction with each slow vehicle, the speed of the shock decreases substantially, as it happens around $t = 1$, $t = 2.5$ and $t = 4$.

Conflict of interest

All authors declare no conflicts of interest in this paper.