**2. Numerical Method**

We begin with the description of the numerical scheme. We must employ the integration along the characteristic curves; therefore, we should remind the reader of the definition of *x*(*t*; *t* <sup>∗</sup>, *x*∗), *x*<sup>∗</sup> ∈ [*x*0, *xM*(*t*)], *t* ∗ > 0, as the solution of

$$\begin{cases} \frac{d}{dt}x(t;t^\*,\mathbf{x}\_\*) = \mathbf{g}(x(t;t^\*,\mathbf{x}\_\*),\mathbf{s}(t),t), \quad t \ge t^\*,\\ x(t^\*;t^\*,\mathbf{x}\_\*) = \mathbf{x}\_\*, \end{cases} \tag{5}$$

that represents the characteristic curve that at time *t* <sup>∗</sup> starts at *x*∗. Now, we consider *w*(*t*; *t* <sup>∗</sup>, *x*∗) = *u*(*x*(*t*; *t* <sup>∗</sup>, *x*∗), *t*), for each *x*<sup>∗</sup> ∈ [*x*0, *xM*(*t*)], *t* ∗ > 0. It satisfies the following initial value problem

$$\begin{cases} \quad \frac{d}{dt}w(t;t^\*,\mathbf{x}\_\*) = -\mu^\*(\mathbf{x}(t;t^\*,\mathbf{x}\_\*),\mathbf{s}(t),t) \, w(t;t^\*,\mathbf{x}\_\*), \quad t > t^\*,\\\quad w(t^\*;t^\*,\mathbf{x}\_\*) = u(\mathbf{x}\_\*,t^\*), \end{cases} \tag{6}$$

where *μ*∗(*x*,*s*(*t*), *t*) = *μ*(*x*,*s*(*t*), *t*) + *gx*(*x*,*s*(*t*), *t*), whose solution is given by the following formula

$$u(\mathbf{x}(t;t^\*,\mathbf{x}\_\*),t) = u(\mathbf{x}\_\*,t^\*) \exp\left\{-\int\_{t^\*}^{t} \mu^\*\left(\mathbf{x}\left(\alpha;t^\*,\mathbf{x}\_\*\right),s(\alpha),\alpha\right) d\alpha\right\}, t \ge t^\*. \tag{7}$$

For the numerical method, we discretize this expression, together with the boundary condition and the initial data in (1), and the initial value problem in (3).

The integration is carried out on a finite time interval [0, *<sup>T</sup>*]. Then, given *<sup>J</sup>*, *<sup>L</sup>* <sup>∈</sup> <sup>N</sup>, we define the discretization parameters *h* = (*x*<sup>0</sup> *<sup>M</sup>* − *x*0)/*J*, *k* = *τ*/*L*, and the number of discrete time levels *N* = [*T*/*k*], which are given by *t <sup>n</sup>* <sup>=</sup> *n k*, <sup>−</sup>*<sup>L</sup>* <sup>≤</sup> *<sup>n</sup>* <sup>≤</sup> *<sup>N</sup>*. We begin the description of the numerical method with the initial values, which, in this case, are the initial size discretization **X**0, *X*<sup>0</sup> *<sup>j</sup>* = *x*<sup>0</sup> + *j h*, <sup>0</sup> <sup>≤</sup> *<sup>j</sup>* <sup>≤</sup> *<sup>J</sup>*, the initial condition on the initial grid, **<sup>U</sup>**0, *<sup>U</sup>*<sup>0</sup> *<sup>j</sup>* = *<sup>u</sup>*0(*X*<sup>0</sup> *<sup>j</sup>* ), 0 ≤ *j* ≤ *J*, and the resource, that has to be initialized on the interval [−*τ*, 0], with the values *<sup>S</sup><sup>n</sup>* <sup>=</sup> *<sup>s</sup>*0(*<sup>t</sup> <sup>n</sup>*), <sup>−</sup>*<sup>L</sup>* <sup>≤</sup> *<sup>n</sup>* <sup>≤</sup> 0.

We compute, for 0 ≤ *n* ≤ *N* − 1, the approximation at the general level *t <sup>n</sup>*<sup>+</sup>1: {**X***n*<sup>+</sup>1, *<sup>S</sup>n*<sup>+</sup>1, **<sup>U</sup>***n*+1}, **<sup>X</sup>***n*+<sup>1</sup> = (*Xn*+<sup>1</sup> <sup>0</sup> , *<sup>X</sup>n*+<sup>1</sup> <sup>1</sup> , ... , *<sup>X</sup>n*+<sup>1</sup> *<sup>J</sup>*+*n*+1), **<sup>U</sup>***n*+<sup>1</sup> = (*Un*+<sup>1</sup> <sup>0</sup> , *<sup>U</sup>n*+<sup>1</sup> <sup>1</sup> , ... , *<sup>U</sup>n*+<sup>1</sup> *<sup>J</sup>*+*n*+1) when the values at the previous time steps are known: {**X***m*, *<sup>S</sup>m*, **<sup>U</sup>***m*}, **<sup>X</sup>***<sup>m</sup>* = (*X<sup>m</sup>* <sup>0</sup> , *<sup>X</sup><sup>m</sup>* <sup>1</sup> , ... , *<sup>X</sup><sup>m</sup> <sup>J</sup>*+*m*), **<sup>U</sup>***<sup>m</sup>* = (*U<sup>m</sup>* <sup>0</sup> , *<sup>U</sup><sup>m</sup>* <sup>1</sup> , ... , *<sup>U</sup><sup>m</sup> <sup>J</sup>*+*m*), 0 ≤ *m* ≤ *n*. We should point out that we design a second-order method based on the modified Euler scheme and on the trapezoidal quadrature rule. This means that we need an intermediate stage in which we

compute a first-order approximation: {**X***n*+1,∗, *<sup>S</sup>n*+1,∗, **<sup>U</sup>***n*+1,∗}, **<sup>X</sup>***n*+1,<sup>∗</sup> = (*Xn*+1,<sup>∗</sup> <sup>0</sup> , *<sup>X</sup>n*+1,<sup>∗</sup> <sup>1</sup> , ... , *<sup>X</sup>n*+1,<sup>∗</sup> *<sup>J</sup>*+*n*+1), **U***n*+1,<sup>∗</sup> = (*Un*+1,<sup>∗</sup> <sup>0</sup> , *<sup>U</sup>n*+1,<sup>∗</sup> <sup>1</sup> ,..., *<sup>U</sup>n*+1,<sup>∗</sup> *<sup>J</sup>*+*n*+1), given by the equations

$$\begin{cases} X\_0^{n+1,\*} = \mathbf{x}\_0, \quad X\_{j+1}^{n+1,\*} = X\_j^n + k \, \mathcal{g}(X\_j^n, S^n, t^n), \quad 0 \le j \le J+n, \\ S^{n+1,\*} = S^n + k \, f(S^n, S^{n-1}, \mathcal{Q}(\mathbf{X}^n, \boldsymbol{\gamma}^n(\mathbf{X}, S) \cdot \mathbf{U}^n), t^n), \\ \mathcal{U}\_{j+1}^{n+1,\*} = \mathcal{U}\_j^n \, \exp\left\{-k \, \boldsymbol{\mu}^\*(X\_j^n, S^n, t^n)\right\}, \quad 0 \le j \le J+n, \\ \mathcal{g}(\mathbf{X}\_0^{n+1,\*}, S^{n+1,\*}, t^{n+1}) \, \mathcal{U}\_0^{n+1,\*} = \mathcal{Q}(\mathbf{X}^{n+1,\*}, \boldsymbol{\mathfrak{a}}^{n+1,\*}(\mathbf{X}, S) \cdot \mathbf{U}^{n+1,\*}). \end{cases} \tag{8}$$

As soon as we have computed the intermedium quantities {**X***n*+1,∗, *<sup>S</sup>n*+1,∗, **<sup>U</sup>***n*+1,∗} we can obtain the approximations at the advanced time level with

⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ *<sup>X</sup>n*+<sup>1</sup> <sup>0</sup> <sup>=</sup> *<sup>x</sup>*0, *<sup>X</sup>n*+<sup>1</sup> *<sup>j</sup>*+<sup>1</sup> = *<sup>X</sup><sup>n</sup> <sup>j</sup>* + *k* 2 *g*(*X<sup>n</sup> <sup>j</sup>* , *<sup>S</sup>n*, *<sup>t</sup> <sup>n</sup>*) + *g*(*Xn*+1,<sup>∗</sup> *<sup>j</sup>*+<sup>1</sup> , *<sup>S</sup>n*+1,∗, *<sup>t</sup> <sup>n</sup>*+1) , 0 ≤ *j* ≤ *J* + *n*, *Sn*+<sup>1</sup> = *S<sup>n</sup>* + *k* 2 *<sup>f</sup>*(*Sn*, *<sup>S</sup>n*−*L*, <sup>Q</sup>(**X***n*, *<sup>γ</sup><sup>n</sup>* (**X**, *<sup>S</sup>*) · **<sup>U</sup>***n*), *<sup>t</sup> n*) <sup>+</sup> *<sup>f</sup>*(*Sn*+1,∗, *<sup>S</sup>n*+1−*L*, <sup>Q</sup>(**X***n*+1,∗, *<sup>γ</sup>n*+1,<sup>∗</sup> (**X**, *<sup>S</sup>*)·**U***n*+1,∗), *<sup>t</sup> <sup>n</sup>*+1) , *Un*+<sup>1</sup> *<sup>j</sup>*+<sup>1</sup> = *<sup>U</sup><sup>n</sup> <sup>j</sup>* exp<sup>0</sup> −*k* 2 *μ*∗(*X<sup>n</sup> <sup>j</sup>* , *<sup>S</sup>n*, *<sup>t</sup> <sup>n</sup>*) + *μ*∗(*Xn*+1,<sup>∗</sup> *<sup>j</sup>*+<sup>1</sup> , *<sup>S</sup>n*+1,∗, *<sup>t</sup> <sup>n</sup>*+1) 1 , 0 ≤ *j* ≤ *J* + *n*, *<sup>g</sup>*(*Xn*+<sup>1</sup> <sup>0</sup> , *<sup>S</sup>n*<sup>+</sup>1, *<sup>t</sup> <sup>n</sup>*+1) *<sup>U</sup>n*+<sup>1</sup> <sup>0</sup> <sup>=</sup> <sup>Q</sup>(**X***n*<sup>+</sup>1, *<sup>α</sup>n*+<sup>1</sup> (**X**, *<sup>S</sup>*) · **<sup>U</sup>***n*+1). (9)

In (8) and (9), at each time level, *γ<sup>n</sup>* (**X**, *S*) and *α<sup>n</sup>* (**X**, *S*) represent vectors with components *γn <sup>j</sup>* (**X**, *<sup>S</sup>*) = *<sup>γ</sup>*(*X<sup>n</sup> <sup>j</sup>* , *<sup>S</sup>n*, *<sup>t</sup> <sup>n</sup>*) and *α<sup>n</sup> <sup>j</sup>* (**X**, *<sup>S</sup>*) = *<sup>α</sup>*(*X<sup>n</sup> <sup>j</sup>* , *<sup>S</sup>n*, *<sup>t</sup> <sup>n</sup>*), 0 <sup>≤</sup> *<sup>j</sup>* <sup>≤</sup> *<sup>J</sup>* <sup>+</sup> *<sup>n</sup>*, respectively, and *<sup>γ</sup>n*,<sup>∗</sup> (**X**, *<sup>S</sup>*) and *αn*,<sup>∗</sup> (**X**, *S*) vectors with components *γn*,<sup>∗</sup> *<sup>j</sup>* (**X**, *<sup>S</sup>*) <sup>=</sup> *<sup>γ</sup>*(*Xn*,<sup>∗</sup> *<sup>j</sup>* , *<sup>S</sup>n*,∗, *<sup>t</sup> <sup>n</sup>*) and *αn*,<sup>∗</sup> *<sup>j</sup>* (**X**, *<sup>S</sup>*) <sup>=</sup> *<sup>α</sup>*(*Xn*,<sup>∗</sup> *<sup>j</sup>* , *<sup>S</sup>n*,∗, *<sup>t</sup> n*), <sup>0</sup> <sup>≤</sup> *<sup>j</sup>* <sup>≤</sup> *<sup>J</sup>* <sup>+</sup> *<sup>n</sup>* <sup>+</sup> 1, respectively. The notation *<sup>γ</sup><sup>n</sup>* (**X**, *<sup>S</sup>*) · **<sup>U</sup>***n*, *<sup>α</sup><sup>n</sup>* (**X**, *<sup>S</sup>*) · **<sup>U</sup>***n*, *<sup>γ</sup>n*,<sup>∗</sup> (**X**, *<sup>S</sup>*) · **<sup>U</sup>***n*,<sup>∗</sup> and *<sup>α</sup>n*,<sup>∗</sup> (**X**, *<sup>S</sup>*) · **<sup>U</sup>***n*,<sup>∗</sup> represent the componentwise products of the corresponding vectors, and <sup>Q</sup> is an appropriate second-order quadrature rule, as given below. We observed that, in this procedure, the total number of nodes in the grid increases in one at each time step due to the new node that fluxes from the left boundary into the domain. This feature led us to name this kind of scheme the aggregation grid nodes method (AGN) [12]. The algorithm is explicit and totally straightforward and can be resumed as in the pseudocode given in Box 1.

**Box 1.** Pseudocode of the algorithm AGN (and SGN).

```
% Choose L, J integers
k = τ/L; h = (x0
               M − x0)/J
N = [T/k]
% Initialization
do j=0,J
       X0
        j = x0 + j h % Initial grid points
      U0
        j = u0(X0
                j )
end
do n=-L,0
      Sn = s0(t
               n)
end
% Advancing the solution on time, step by step
do n=0, N-1
      Formulae in (8) % First stage to Xn+1,∗, Sn+1,∗, Un+1,∗
      Formulae in (9) % Second stage to Xn+1, Sn+1, Un+1
% Optional (difference between AGN and SGN methods)
      Formulae in (10) % Selection procedure in case of SGN method
end
```
In the numerical experimentation, we use a more efficient procedure in order to keep the number of nodes constant (and consequently, to not increase the computational cost) at each time step. Toward that goal, we eliminate a chosen node from the grid at each time step. The selection procedure we use in the method depends on the dynamics of the grid, as it was introduced in [13]. We refer to this kind of scheme as the selection grid nodes method (SGN) [12]. The node picked out is the first *Xn*+<sup>1</sup> *l* that satisfies

$$|X\_{l+1}^{n+1} - X\_{l-1}^{n+1}| = \min\_{1 \le j \le l} |X\_{j+1}^{n+1} - X\_{j-1}^{n+1}|.\tag{10}$$

Once it is elected, we rename the nodes to keep the index labels from *j* = 0 to *j* = *J*.

Finally, we pay attention to the nonlocal terms; they are discretized by means of a composite quadrature rule, based on the trapezoidal formula, on the grid points **<sup>X</sup>** = (*X*0, *<sup>X</sup>*1,..., *<sup>X</sup>*ℵ), ℵ ∈ <sup>N</sup>,

$$\mathcal{Q}(\mathbf{X}, \mathbf{V}) = \sum\_{l=0}^{\aleph - 1} \frac{X\_{l+1} - X\_l}{2} \left( V\_l + V\_{l+1} \right), \quad \text{where} \quad \mathbf{V} = (V\_0, V\_1, \dots, V\_{\aleph}).$$

#### **3. Convergence Analysis**

The convergence analysis here follows the discretization framework developed in [14]. In this framework, suitable discrete normed spaces and operators are introduced to formulate the equations of the numerical method. Then, appropriate properties of consistency and nonlinear stability are established. The numerical method is an adaptation of the one given in [9], and so is the analysis. However, for the sake of completeness, we provide the details of the numerical method formulation, and for the sake of simplicity, we present only those parts of the convergence analysis proofs associated with estimates derived from the discretization of the delay term in the differential equation. That means that the step by step recurrences for the errors at each time level increase the order of the difference equation as the time discretization parameter tends to zero. For the theoretical analysis, we consider the AGN method.

With this aim, we fix *T* > 0 and assume that the problem (1), (3) has a unique solution *u*(*x*, *t*), *x* ∈ [*x*0, *xM*(*t*)], and *s*(*t*), *t* ∈ [−*τ*, *T*], with the following regularity assumptions:

**Hypothesis 1** (H1)**.** *<sup>u</sup>* ∈ C2([*x*0, *xM*(*t*)] <sup>×</sup> [0, *<sup>T</sup>*]) *and is nonnegative.*

**Hypothesis 2** (H2)**.** *<sup>s</sup>* ∈ C2([−*τ*, *<sup>T</sup>*]) *and is nonnegative.*

Additionally, we assume that there exists a compact neighborhood *D* of {*s*(*t*), 0 ≤ *t* ≤ *T*}, such that

**Hypothesis 3** (H3)**.** *<sup>γ</sup>* ∈ C2([*x*0, *xM*(*t*)] <sup>×</sup> *<sup>D</sup>* <sup>×</sup> [0, *<sup>T</sup>*]) *and is nonnegative.*

**Hypothesis 4** (H4)**.** *<sup>μ</sup>* ∈ C2([*x*0, *xM*(*t*)] <sup>×</sup> *<sup>D</sup>* <sup>×</sup> [0, *<sup>T</sup>*]) *and is nonnegative.*

**Hypothesis 5** (H5)**.** *<sup>α</sup>* ∈ C2([*x*0, *xM*(*t*)] <sup>×</sup> *<sup>D</sup>* <sup>×</sup> [0, *<sup>T</sup>*]) *and is nonnegative.*

**Hypothesis 6** (H6)**.** *<sup>g</sup>* ∈ C3([*x*0, *xM*(*t*)] <sup>×</sup> *<sup>D</sup>* <sup>×</sup> [0, *<sup>T</sup>*]) *and there exists a positive constant <sup>C</sup> such that g*(*x*0,*s*, *t*) ≥ *C, s*, *t* ≥ 0*. In addition, the characteristic curves x*(*t*; *t* <sup>∗</sup>, *x*∗) *defined by (5), are continuous and differentiable with respect to the initial values* (*t* <sup>∗</sup>, *x*∗) ∈ [0, *T*] × [*x*0, *xM*(*t*)]*.*

Finally, we suppose that there are compact neighborhoods, *Df* of {*s*(*t*), −*τ* ≤ *t* ≤ *T*}, and *DI* of 0 *xM*(*t*) *x*0 *γ*(*x*,*s*(*t*), *t*) *u*(*x*, *t*) *dx*, 0 ≤ *t* ≤ *T* 1 , such that

**Hypothesis 7** (H7)**.** *<sup>f</sup>* ∈ C2(*Df* <sup>×</sup> *Df* <sup>×</sup> *DI* <sup>×</sup> [0, *<sup>T</sup>*]) *and is nonnegative.*

Now, we choose *<sup>L</sup>* <sup>∈</sup> <sup>N</sup> and assume that the time discretization parameter, *<sup>k</sup>*, takes values in the set

$$K = \{ k : k = \pi/(\nu L), \,\nu \in \mathbb{N} \}.$$

Then for *<sup>k</sup>* <sup>∈</sup> *<sup>K</sup>*, we set *<sup>N</sup>* = [*T*/*k*] and choose *<sup>J</sup>* <sup>∈</sup> <sup>N</sup> such that *<sup>h</sup>* = (*x*<sup>0</sup> *<sup>M</sup>* − *x*0)/*J*, and the ratio *r* := *k*/*h* is a positive constant fixed throughout the analysis. Then, *x*<sup>0</sup> *<sup>M</sup>* is always the last node in the size grid. For each *k* ∈ *K*, we define the space

$$\mathcal{A}\_k = \prod\_{n=0}^N \left( \mathbb{R}^{J+n} \times \mathbb{R}^{J+n+1} \right) \times \mathbb{R}^{L+N+1} \text{.}$$

where vectors **y**0, **V**0,..., **y***N*, **V***N*, **a** ∈ A*<sup>k</sup>* are used to describe the following approximations: to the inner and right-hand boundary grid nodes, given by **<sup>y</sup>***<sup>n</sup>* <sup>∈</sup> <sup>R</sup>*J*+*n*, and to the theoretical solution at the complete grid, **<sup>V</sup>***<sup>n</sup>* <sup>∈</sup> <sup>R</sup>*J*+*n*<sup>+</sup>1, for each time level *<sup>t</sup> <sup>n</sup>*, 0 <sup>≤</sup> *<sup>n</sup>* <sup>≤</sup> *<sup>N</sup>*; and to the theoretical solution to the delay differential equation at time levels *t <sup>n</sup>*, <sup>−</sup>*<sup>L</sup>* <sup>≤</sup> *<sup>n</sup>* <sup>≤</sup> *<sup>N</sup>*, provided by **<sup>a</sup>** <sup>∈</sup> <sup>R</sup>*L*+*N*<sup>+</sup>1. We also consider the space

$$\mathcal{B}\_k = \left(\mathbb{R}^I \times \mathbb{R}^{I+1} \times \mathbb{R}^{L+1}\right) \times \mathbb{R}^N \times \prod\_{n=1}^N \left(\mathbb{R}^{I+n} \times \mathbb{R}^{I+n}\right) \times \mathbb{R}^N,$$

where each component of vector **Y**0, **P**0, **A**0, **P**0, **Y**1, **P**1,..., **Y***N*, **P***N*, **A** ∈ B*<sup>k</sup>* describes the residuals of the discrete solution when it is substituted on the discrete equations that define the numerical method: we discriminate **Y**0, **P**0, **A**<sup>0</sup> = (*A*−*L*, *A*−*L*<sup>+</sup>1, ... , *A*0), given by the approximations to the initial conditions; **P**0, provided by the solution at the left boundary node at each time level; **Y***n*, **P***n*, <sup>0</sup> <sup>≤</sup> *<sup>n</sup>* <sup>≤</sup> *<sup>N</sup>*, given by the grid nodes and solution to them; and **<sup>A</sup>** = (*A*1, *<sup>A</sup>*2, ... , *<sup>A</sup>N*), due to the solution to the delay differential equation, as we show below. Both spaces, A*<sup>k</sup>* and B*k*, have the same dimension.

In order to measure the size of the errors, we define

$$\|\eta\|\|\_{\infty} = \max\_{1 \le j \le \aleph} |\eta\_j|\prime \quad \eta = (\eta\_1, \eta\_2, \dots, \eta\_{\aleph}) \in \mathbb{R}^{\aleph},$$

$$\|\mathbf{V}^n\|\|\_1 = \sum\_{j=0}^{J+n} h \, |V\_j^n\vert \, \prime \quad \mathbf{V}^n = (V\_0^n, V\_1^n, \dots, V\_{J+n}^n) \in \mathbb{R}^{J+n+1}.$$

Thus, we endow the spaces <sup>A</sup>*<sup>k</sup>* and <sup>B</sup>*<sup>k</sup>* with the following norms. If **y**0, **V**0,..., **y***N*, **V***N*, **a** ∈ A*k*,

$$\|\left(\mathbf{y}^{0},\mathbf{V}^{0},\dots,\mathbf{y}^{N},\mathbf{V}^{N},\mathbf{a}\right)\|\_{\mathcal{A}\_{k}} = \max\left(\max\_{0\le n\le N} \|\mathbf{y}^{n}\|\_{\infty}, \max\_{0\le n\le N} \|\mathbf{V}^{n}\|\_{\infty}, \|\mathbf{a}\|\_{\infty}\right).$$

On the other hand, if **Y**0, **P**0, **A**0, **P**0, **Y**1, **P**1,..., **Y***N*, **P***N*, **A** ∈ B*k*,

$$\begin{split} & \| \left( \mathbf{Y}^{0}, \mathbf{P}^{0}, \mathbf{A}^{0}, \mathbf{P}\_{0}, \mathbf{Y}^{1}, \mathbf{P}^{1}, \dots, \mathbf{Y}^{N}, \mathbf{P}^{N}, \mathbf{A} \right) \|\_{\mathcal{B}\_{k}} \\ &= \| \mathbf{Y}^{0} \|\_{\infty} + \| \mathbf{P}^{0} \|\_{\infty} + \| \mathbf{A}^{0} \|\_{\infty} + \| \mathbf{P}\_{0} \|\_{\infty} + \sum\_{n=1}^{N} k \, \| \mathbf{Y}^{n} \|\_{\infty} + \sum\_{n=1}^{N} k \, \| \mathbf{P}^{n} \|\_{\infty} + \sum\_{n=1}^{N} k \, \| A^{n} \|. \end{split} $$

Now, for each *k* ∈ *K*, we define

$$\mathbf{x}^{n} = (\mathbf{x}\_{1}^{n}, \dots, \mathbf{x}\_{j+n}^{n}) \in \mathbb{R}^{I+n}, \mathbf{x}\_{j}^{0} = \mathbf{x}\_{0} + jh, \ 1 \le j \le I; \ \mathbf{x}\_{j}^{n} = \mathbf{x}(t\_{n}; t\_{n-1}, \mathbf{x}\_{j-1}^{n-1}), \ 1 \le j \le I + n, \ 1 \le n \le N; \tag{11}$$

and we denote *x<sup>n</sup>* <sup>0</sup> = *x*0, *n* ≥ 0. Recall that *x*(*t*; *t* <sup>∗</sup>, *x*∗) represents the theoretical solution to problem (5), *t* <sup>∗</sup> ∈ [0, *T*], *x*<sup>∗</sup> ∈ [*x*0, *xM*(*t*)]. In addition, if *u* represents the theoretical solution to (1) we define

$$\mathbf{u}^{n} = (u\_0^n, u\_1^n, \dots, u\_{f+n}^n) \in \mathbb{R}^{f+n+1}, \quad u\_j^n = \mathfrak{u}(\mathbf{x}\_j^n, t\_n), \quad 0 \le j \le f+n, \quad 0 \le n \le N. \tag{12}$$

Finally, if *s* is the theoretical solution to (3) then we define

$$\mathbf{s}\_{k} = (\mathbf{s}^{-L}, \mathbf{s}^{-L+1}, \dots, \mathbf{s}^{N}), \quad \mathbf{s}^{\boldsymbol{n}} = \mathbf{s}(t^{\boldsymbol{n}}), \quad -L \le \boldsymbol{n} \le N. \tag{13}$$

Therefore, **<sup>u</sup>**˜ *<sup>k</sup>* = (**x**0, **<sup>u</sup>**0, **<sup>x</sup>**1, **<sup>u</sup>**1,..., **<sup>x</sup>***N*, **<sup>u</sup>***N*, **<sup>s</sup>***k*) ∈ A*k*.

In the following, we introduce the discretization operator. Let *R* be a positive constant and let we denote by *<sup>B</sup>*A*<sup>k</sup>* (**u**˜ *<sup>k</sup>*, *R kp*) ⊂ A*<sup>k</sup>* the open ball with center **<sup>u</sup>**˜ *<sup>k</sup>* and radius *R kp*, 1 <sup>&</sup>lt; *<sup>p</sup>* <sup>&</sup>lt; 2. Then, we define

$$\Phi\_k: B\_{\mathcal{A}\_k}(\mathfrak{u}\_k, R\, k^P) \to \mathcal{B}\_{k\prime}$$

$$\Phi\_k\left(\mathbf{y}^0, \mathbf{V}^0, \dots, \mathbf{y}^N, \mathbf{V}^N, \mathbf{a}\right) = \left(\mathbf{Y}^0, \mathbf{P}^0, \mathbf{A}^0, \mathbf{P}\_{0\prime}, \mathbf{Y}^1, \mathbf{P}^1, \dots, \mathbf{Y}^N, \mathbf{P}^N, \mathbf{A}\right),\tag{14}$$

given by the following equations, first

$$\mathbf{y}^{0} = \mathbf{y}^{0} - \mathbf{x}^{0} \in \mathbb{R}^{I} \, \tag{15}$$

$$\mathbf{P}^0 = \mathbf{V}^0 - \mathbf{U}^0 \in \mathbb{R}^{I+1},\tag{16}$$

$$\mathbf{A}^0 = \mathbf{a}^0 - \mathbf{S}^0 \in \mathbb{R}^{L+1},\tag{17}$$

where **a**<sup>0</sup> = (*a*−*L*, *a*−*L*<sup>+</sup>1, ... , *a*0). Vectors **X**<sup>0</sup> and **U**<sup>0</sup> represent approximations at *t* = 0, respectively, to the initial grid nodes and to the theoretical solution at these nodes. Vector **S**<sup>0</sup> represents an approximation to the initial resource in the interval [−*τ*, 0]. Second,

$$P\_0^{n+1} = V\_0^{n+1} - \frac{\mathcal{Q}\left(\mathbf{y}^{n+1}, \mathbf{a}^{n+1}(\mathbf{y}, a) \cdot \mathbf{V}^{n+1}\right)}{\mathcal{g}\left(\mathbf{x}\_0, a^{n+1}, t^{n+1}\right)},\tag{18}$$

$$\mathcal{Y}\_{j+1}^{n+1} = \frac{1}{k} \left\{ y\_{j+1}^{n+1} - y\_j^n - \frac{k}{2} \left( \mathbf{g}(y\_{j'}^n, a^n, t^n) + \mathbf{g}(y\_{j+1}^{n+1,\*}, a^{n+1,\*}, t^{n+1}) \right) \right\},\tag{19}$$

$$P\_{j+1}^{n+1} = \frac{1}{k} \left\{ V\_{j+1}^{n+1} - V\_j^n \exp\left( -\frac{k}{2} \left( \mu^\* \left( y\_j^n, a^n, t^n \right) + \mu^\* \left( y\_{j+1}^{n+1,\*}, a^{n+1,\*}, t^{n+1} \right) \right) \right) \right\},\tag{20}$$

0 ≤ *j* ≤ *J* + *n* − 1,

$$A^{n+1} = \frac{1}{k} \left\{ a^{n+1} - a^n - \frac{k}{2} \left( f\left( a^n, a^{n-1}, \mathcal{Q}(\mathbf{y}^n, \boldsymbol{\gamma}^n(\mathbf{y}, a) \cdot \mathbf{V}^n), t^n \right) \right. \right. \\ \left. + f^n \left( a^{n+1,\*} , \mathcal{Q}(\mathbf{y}^{n+1,\*}, \boldsymbol{\gamma}^{n+1,\*}(\mathbf{y}, a) \cdot \mathbf{V}^{n+1,\*}) \right) \right), \tag{21}$$
 
$$+ f\left( a^{n+1,\*}, a^{n-L+1}, \mathcal{Q}(\mathbf{y}^{n+1,\*}, \boldsymbol{\gamma}^{n+1,\*}(\mathbf{y}, a) \cdot \mathbf{V}^{n+1,\*}), t^{n+1} \right)) \; \right\}, \tag{22}$$

0 ≤ *n* ≤ *N* − 1. Where, with the notation introduced in Section 2,

$$(y\_{j+1}^{n+1,\*} = y\_j^n + k \, g(y\_j^n, a^n, t^n), \tag{22}$$

$$V\_{j+1}^{n+1,\*} = V\_j^n \exp\left(-k\,\mu^\* \left(y\_{j}^n, a^n, t^n\right)\right),\tag{23}$$

0 ≤ *j* ≤ *J* + *n* − 1,

$$\mathcal{V}\_0^{n+1,\*} = \frac{\mathcal{Q}(\mathbf{y}^{n+1,\*}, \mathbf{a}^{n+1,\*} (\mathbf{y}, a) \cdot \mathbf{V}^{n+1,\*})}{\mathcal{g}(\mathbf{x}\_0, a^{n+1,\*}, t\_{n+1})},\tag{24}$$

$$a^{n+1,\*} = a^n + k \left( a^n, a^{n-L}, \mathbb{Q}(\mathbf{y}^n, \gamma^n(\mathbf{y}, a) \cdot \mathbf{V}^n), t^n \right), \tag{25}$$

0 ≤ *n* ≤ *N* − 1. We rewrite the quadrature rule employed in (18)–(25) in the next general form Q(**X**, **V**) = *J*+*n* ∑ *l*=0 *ql*(**X**) *Vl*. Written in this way, we highlight that the number of nodes considered at each time level *t <sup>n</sup>* is *J* + *n* + 1, counting both the boundary node *X<sup>n</sup>* <sup>0</sup> = *<sup>x</sup><sup>n</sup>* <sup>0</sup> = *x*0, 0 ≤ *n* ≤ *N* and the interior grid nodes *X<sup>n</sup> <sup>j</sup>* , 1 ≤ *j* ≤ *J* + *n*, 0 ≤ *n* ≤ *N*. This notation is also valid even when we consider quadrature rules whose nodes are, at each time level, a subgrid of **<sup>X</sup>***n*, 0 <sup>≤</sup> *<sup>n</sup>* <sup>≤</sup> *<sup>N</sup>*, by asumming that the corresponding weights *ql*(**X**) of some of the nodes are zero .

Note that **Φ***<sup>k</sup>* takes into account all the grid nodes and their corresponding solution values at each time level, and it employs quadrature rules possibly based on a subgrid. If **U**˜ *<sup>k</sup>* = (**X**0, **<sup>U</sup>**0, **<sup>X</sup>**1, **<sup>U</sup>**1,..., **<sup>X</sup>***N*, **<sup>U</sup>***N*, **<sup>s</sup>***k*) <sup>∈</sup> *<sup>B</sup>*A*<sup>k</sup>* (**u**˜ *<sup>k</sup>*, *R kp*), satisfies

$$\Phi\_k(\mathbf{\tilde{U}}\_k) = \mathbf{0} \in \mathcal{B}\_{k\prime} \tag{26}$$

the nodes, **<sup>X</sup>***<sup>n</sup>* and 0 <sup>≤</sup> *<sup>n</sup>* <sup>≤</sup> *<sup>N</sup>* and the corresponding values of the solution, **<sup>U</sup>***n*, 0 <sup>≤</sup> *<sup>n</sup>* <sup>≤</sup> *<sup>N</sup>*, at such nodes are numerical solutions to the scheme defined by (9) when the composite trapezoidal quadrature rule is used on a subgrid at each time step. On the other hand, the numerical solution to the scheme defined by (9) satisfies (26).

Henceforth, *C* will denote a positive constant, independent of *k*, *h* (*k* = *r h*), *j* (0 ≤ *j* ≤ *J* + *n*) and *n* (−*L* ≤ *n* ≤ *N*), and *C* may have different values at different places.

As in [12], it is possible to establish the following property that we called (SG) for the selection procedure given by (10): if \$ *xn j n l* %*M*(*n*) *<sup>l</sup>*=<sup>0</sup> , denotes a subgrid of the grid \$ *xn j* %*J*+*n*+<sup>1</sup> *<sup>j</sup>*=<sup>0</sup> at the time level *<sup>t</sup> n*, 0 ≤ *n* ≤ *N*,

(SG) There exists a positive constant *C* such that, for *h* sufficiently small, *x<sup>n</sup> j n l*+1 <sup>−</sup> *<sup>x</sup><sup>n</sup> j n l* ≤ *C h*, <sup>0</sup> <sup>≤</sup> *<sup>l</sup>* <sup>≤</sup> *<sup>M</sup>*(*n*) <sup>−</sup> 1, *<sup>x</sup><sup>n</sup> j n* 0 = *x*0, *x<sup>n</sup> j n M*(*n*) = *x<sup>n</sup> <sup>J</sup>*+*n*, 0 ≤ *n* ≤ *N*.

The property (SG) condenses the essential information about the adaptive selection procedure to allow the proof, under the hypotheses (H1)–(H7), of the following general properties of the composite trapezoidal quadrature Q (**X**, **V**) based on the subgrids,

(P1) |*I*(*t <sup>n</sup>*) − Q (**x***n*, *<sup>γ</sup>n*(**x**,*s*) · **<sup>u</sup>***n*)<sup>|</sup> <sup>≤</sup> *C h*2, when *<sup>h</sup>* <sup>→</sup> 0, 0 <sup>≤</sup> *<sup>n</sup>* <sup>≤</sup> *<sup>N</sup>*.

(P2) *xM*(*t n*) *x*0 *α*(*x*,*s*(*t <sup>n</sup>*), *t <sup>n</sup>*) *u*(*x*, *t <sup>n</sup>*) *dx* − Q (**x***n*, *<sup>α</sup>n*(**x**,*s*) · **<sup>u</sup>***n*) <sup>≤</sup> *C h*2, when *<sup>h</sup>* <sup>→</sup> 0, 0 <sup>≤</sup> *<sup>n</sup>* <sup>≤</sup> *<sup>N</sup>*.


$$\left| \sum\_{i=0}^{J+n} \left( q\_i(\mathbf{y}^n) - q\_i(\mathbf{z}^n) \right) \gamma(z\_i^n, a^n, t^n) \, V\_i^n \right| \le C ||\mathbf{y}^n - \mathbf{z}^n||\_{\infty} $$

when *k* → 0, 0 ≤ *n* ≤ *N*.

(P6) Let *<sup>R</sup>* and *<sup>p</sup>* be positive constants with 1 <sup>&</sup>lt; *<sup>p</sup>* <sup>&</sup>lt; 2. If **<sup>y</sup>***n*, **<sup>z</sup>***<sup>n</sup>* <sup>∈</sup> *<sup>B</sup>*∞(**x***n*, *R kp*), **<sup>V</sup>***<sup>n</sup>* <sup>∈</sup> *<sup>B</sup>*∞(**u***n*, *R kp*) and *<sup>a</sup><sup>n</sup>* <sup>∈</sup> *<sup>B</sup>*∞(*sn*, *R kp*); then

$$\left| \sum\_{i=0}^{l+n} \left( q\_i(\mathbf{y}^n) - q\_i(\mathbf{z}^n) \right) \text{ a } (z\_i^n, a^n, t^n) \text{ } V\_i^n \right| \le \mathbb{C} ||\mathbf{y}^n - \mathbf{z}^n||\_{\infty} $$

when *h* → 0, 0 ≤ *n* ≤ *N*.

*Mathematics* **2020**, *8*, 1440

Then, the properties (P1)–(P6) allow us to establish that the nonlinear discrete operators describing **Φ***<sup>k</sup>* are well defined for *k* ∈ *K* small enough, as we formulate in the following theorem (we refer to [9] about details of the proof).

**Theorem 1.** *Assume that hypotheses (H1)–(H7) hold and that the quadrature rules used in (18)–(25) satisfy the properties (P1)–(P6). If*

$$\left(\mathbf{X}^0, \mathbf{V}^0, \dots, \mathbf{X}^N, \mathbf{V}^N, \mathbf{S}\right) \in B\_{\mathcal{A}\_k}(\bar{\mathbf{u}}\_{k'} \, R \, k^p)\_{'},$$

*where R is a fixed positive constant and* 1 < *p* < 2*, then, for k sufficiently small,*

$$\mathbb{Q}(\mathbf{X}^n, \gamma^n(\mathbf{X}, \mathbf{S}) \cdot \mathbf{V}^n) \in D\_{I\_1} \tag{27}$$

0 ≤ *n* ≤ *N. Furthermore, there exists a positive constant R , independent of k, such that, as k* → 0*,* **<sup>X</sup>***n*,<sup>∗</sup> <sup>∈</sup> *<sup>B</sup>*∞(**x***n*, *<sup>R</sup> <sup>k</sup>p*)*, Sn*,<sup>∗</sup> <sup>∈</sup> *<sup>B</sup>*∞(*sn*, *<sup>R</sup> <sup>k</sup>p*) *and* **<sup>V</sup>***n*,<sup>∗</sup> <sup>∈</sup> *<sup>B</sup>*∞(**u***n*, *<sup>R</sup> <sup>k</sup>p*)*, and*

$$\mathcal{Q}(\mathbf{X}^{n,\*}, \gamma^{n,\*}(\mathbf{X}, \mathbf{S}) \cdot \mathbf{V}^{n,\*}) \in D\_{I\prime} \tag{28}$$

1 ≤ *n* ≤ *N.*

Now, we define the local discretization error as

$$\boldsymbol{\Phi}\_k(\boldsymbol{\upmu}\_k) = (\mathbf{Z}^0, \mathbf{L}^0, \sigma^0, \mathbf{L}\_0, \mathbf{Z}^1, \mathbf{L}^1, \dots, \mathbf{Z}^N, \mathbf{L}^N, \sigma) \in \mathcal{B}\_{k,\sigma}$$

and we say that the discretization (14) is consistent if, as *k* → 0,

$$\lim \| \Phi\_k(\vec{\mathbf{u}}\_k) \|\_{\mathcal{B}\_k} = 0.$$

The following theorem establishes the consistency of the numerical scheme defined by Equations (26).

**Theorem 2** (Consistency)**.** *Assume that hypotheses (H1)–(H7) hold and that the considered quadrature rules satisfy properties (P1)–(P6). Then, as k* → 0*, the local discretization error satisfies*

$$\|\Phi\_k(\bar{\mathbf{u}}\_k)\|\_{\mathcal{B}\_k} = \|\mathbf{u}^0 - \mathbf{U}^0\|\_{\infty} + \|\mathbf{x}^0 - \mathbf{X}^0\|\_{\infty} + \max\_{-L \le n \le 0} |\mathbf{s}^n - S^n| + O(h^2 + k^2). \tag{29}$$

**Proof.** The only novelty with respect to the proof of the consistency in [9] is to establish bounds for the truncation errors in the numerical approximations to the solution of the delay differential equation that drives the dynamics of the resource. Therefore, the regularity hypotheses (H1)–(H7), properties (P1), (P2) of the quadrature rule and error bounds for the explicit Euler method allow us to bound

$$\begin{split} |s^{n} - s^{n,\*}| &\leq \left| s(t^{n}) - s^{n-1} - k \left( s^{n-1}, s^{n-1-L}, I(t^{n-1}), t^{n-1} \right) \right| \\ &+ k \left| f\left( s^{n-1}, s^{n-1-L}, I(t^{n-1}), t^{n-1} \right) - f\left( s^{n-1}, s^{n-1-L}, \mathcal{Q}(\mathbf{x}^{n-1}, \boldsymbol{\gamma}^{n-1}(\mathbf{x}, \mathbf{s}) \cdot \mathbf{u}^{n-1}, t^{n-1} \right) \right) \right| \\ &\leq \mathcal{C} \left( k^{2} + k \left| I(t^{n-1}) - \mathcal{Q}(\mathbf{x}^{n-1}, \boldsymbol{\gamma}^{n-1}(\mathbf{x}, \mathbf{s}) \cdot \mathbf{u}^{n-1}) \right| \right) \\ &\leq \mathcal{C} \left( k^{2} + k \, h^{2} \right) . \end{split} \tag{30}$$

As was made in [9], for *h* sufficiently small, we also can attain,

$$|\left|\mathcal{Q}\left(\mathbf{x}^{\mathrm{n}},\gamma^{\mathrm{n}}(\mathbf{x},\mathbf{s})\,\mathbf{u}^{\mathrm{n}}\right)-\mathcal{Q}\left(\mathbf{x}^{\mathrm{n},\ast},\gamma^{\mathrm{n},\ast}(\mathbf{x},\mathbf{s})\,\mathbf{u}^{\mathrm{n},\ast}\right)\right|\leq C\left(h^{2}+k^{2}\right),\tag{31}$$

and bounds on the truncation errors corresponding to the numerical approximations to the nodal grids, **<sup>Z</sup>***n*, 1 <sup>≤</sup> *<sup>n</sup>* <sup>≤</sup> *<sup>N</sup>*,

$$|Z\_j^n| \quad \le \quad \mathcal{C}(k^2 + h^2 k), \quad 1 \le j \le J + n, \quad 1 \le n \le N. \tag{32}$$

Next, (21), the regularity hypotheses (H1)–(H7), the property (P1), inequalities (30) and (31) and the error bound of the modified Euler scheme employed, allow us to achieve

<sup>|</sup>*σn*| ≤ <sup>1</sup> *k* 0 *s <sup>n</sup>* <sup>−</sup> *<sup>s</sup> <sup>n</sup>*−<sup>1</sup> <sup>−</sup> *<sup>k</sup>* 2 *f s <sup>n</sup>*−1,*s <sup>n</sup>*−1−*L*, *I*(*t <sup>n</sup>*−1), *t n*−1 + *f s <sup>n</sup>*,*s <sup>n</sup>*−*L*, *I*(*t <sup>n</sup>*), *t n* + *k* 2 *f s <sup>n</sup>*−1,*s <sup>n</sup>*−1−*L*, *I*(*t <sup>n</sup>*−1), *t n*−1 − *f s <sup>n</sup>*−1,*s <sup>n</sup>*−1−*L*, <sup>Q</sup>(**x***n*−1, *<sup>γ</sup>n*−1(**x**,*s*) · **<sup>u</sup>***n*−1), *<sup>t</sup> n*−1 + *k* 2 *f s <sup>n</sup>*,*s <sup>n</sup>*−*L*, *I*(*t <sup>n</sup>*), *t n* − *f s <sup>n</sup>*,∗,*s <sup>n</sup>*−*L*, *I*(*t <sup>n</sup>*), *t n* + *k* 2 *f s <sup>n</sup>*,∗,*s <sup>n</sup>*−*L*, *I*(*t <sup>n</sup>*), *t n* − *f s <sup>n</sup>*,∗,*s <sup>n</sup>*−*L*, <sup>Q</sup>(**x***n*,∗, *<sup>γ</sup>n*,∗(**x**,*s*) · **<sup>u</sup>***n*,∗), *<sup>t</sup> n* 1 ≤ *C* \$ *<sup>k</sup>*<sup>2</sup> <sup>+</sup> <sup>|</sup>*<sup>s</sup> <sup>n</sup>* <sup>−</sup> *<sup>s</sup> <sup>n</sup>*,∗| <sup>+</sup> *I*(*t <sup>n</sup>*−1) − Q(**x***n*−1, *<sup>γ</sup>n*−1(**x**,*s*) · **<sup>u</sup>***n*−1) + |*I*(*t <sup>n</sup>*) − Q(**x***n*, *<sup>γ</sup>n*(**x**,*s*) · **<sup>u</sup>***n*)<sup>|</sup> <sup>+</sup> |Q(**x***n*, *<sup>γ</sup>n*(**x**,*s*) · **<sup>u</sup>***n*) − Q(**x***n*,∗, *<sup>γ</sup>n*,∗(**x**,*s*) · **<sup>u</sup>***n*,∗)|} ≤ *C k*<sup>2</sup> + *h*<sup>2</sup> , (33)

<sup>1</sup> <sup>≤</sup> *<sup>n</sup>* <sup>≤</sup> *<sup>N</sup>*. Finally, the bounds for the truncation errors produced by *<sup>U</sup><sup>n</sup> <sup>j</sup>* , 0 ≤ *j* ≤ *J* + *n*, 1 ≤ *n* ≤ *N*, are derived as in [9],

$$|L\_j^n| \quad \le \quad \mathcal{C}(k^2 + h^2), \quad 0 \le j \le J + n, \quad 1 \le n \le N. \tag{34}$$

Therefore, (29) follows from (32)–(34).

Another piece of notion that plays an important role in the analysis of the numerical method is the stability with k-dependent thresholds. For *k* ∈ *K*, let *Rk* be a real number ( the stability threshold) with 0 < *Rk* < ∞. We say that the discretization (26) is stable for **u**˜ *<sup>k</sup>* restricted to the thresholds *Rk*, if there exist two positive constants *k*<sup>0</sup> ∈ *K* and S ( the stability constant) such that, for any *k* ∈ *K* with *<sup>k</sup>* <sup>≤</sup> *<sup>k</sup>*0, the open ball *<sup>B</sup>*A*<sup>k</sup>* (**u**˜ *<sup>k</sup>*, *Rk*) is contained in the domain of **<sup>Φ</sup>***k*" and for all **<sup>V</sup>**˜ *<sup>k</sup>*, **<sup>W</sup>**˜ *<sup>k</sup>* in that ball,

$$\|\|\mathbf{\bar{V}}\_{k} - \mathbf{\bar{W}}\_{k}\|\|\_{\mathcal{A}\_{k}} \leq \circledast \|\|\boldsymbol{\Phi}\_{k}(\mathbf{\bar{V}}\_{k}) - \boldsymbol{\Phi}\_{k}(\mathbf{\bar{W}}\_{k})\|\|\_{\mathcal{B}\_{k}}.$$

We introduce the following auxiliary result, proved in [9] where the same quadrature rule and the same procedure of selection of the subgrid at each time step were used.

**Proposition 1.** *Assume that hypotheses (H1)–(H7) hold and that the considered quadrature rules satisfy properties (P1)–(P6). Let be* **<sup>y</sup>***n*, **<sup>z</sup>***<sup>n</sup>* <sup>∈</sup> *<sup>B</sup>*∞(**x***n*, *R kp*)*,* **<sup>V</sup>***n*, **<sup>W</sup>***<sup>n</sup>* <sup>∈</sup> *<sup>B</sup>*∞(**u***n*, *R kp*) *and <sup>a</sup>n, <sup>b</sup><sup>n</sup>* <sup>∈</sup> *<sup>B</sup>*∞(*sn*, *R kp*)*, where R is a positive constant and* 1 < *p* < 2*. Then, as k* → 0*,*

$$\left| \left| \mathcal{Q}(\mathbf{y}^{n}, \gamma^{n}(\mathbf{y}, a) \cdot \mathbf{V}^{n}) - \mathcal{Q}(\mathbf{z}^{n}, \gamma^{n}(\mathbf{z}, b) \cdot \mathbf{W}^{n}) \right| \right| \leq \mathbb{C} \left( \left\| \mathbf{V}^{n} - \mathbf{W}^{n} \right\|\_{1} + \left\| \mathbf{y}^{n} - \mathbf{z}^{n} \right\|\_{\infty} + \left| a^{n} - b^{n} \right| \right), \tag{35}$$

$$|\mathcal{Q}(\mathbf{y}^{n,\*},\boldsymbol{\gamma}^{n,\*}(\mathbf{y},\boldsymbol{a})\cdot\mathbf{V}^{n,\*})-\mathcal{Q}(\mathbf{z}^{n,\*},\boldsymbol{\gamma}^{n,\*}(\mathbf{z},\boldsymbol{b})\cdot\mathbf{W}^{n,\*})| \leq \mathbb{C}\left(||\mathbf{V}^{n,\*}-\mathbf{W}^{n,\*}||\_{1}+||\mathbf{y}^{n,\*}-\mathbf{z}^{n,\*}||\_{\infty}+|\boldsymbol{a}^{n,\*}-\boldsymbol{b}^{n,\*}|\right),\tag{36}$$

1 ≤ *n* ≤ *N.*

Now, we introduce the theorem that establishes the *stability* of the discretization defined by Equations (26).

**Theorem 3** (Stability)**.** *Assume that hypotheses (H1)–(H7) hold and that the considered quadrature rules satisfy properties (P1)–(P6). Then, the discretization is stable for* **u**˜ *<sup>k</sup> with Rk* = *R kp,* 1 < *p* < 2*.*

*Mathematics* **2020**, *8*, 1440

**Proof.** We denote

for

$$
\Phi\_k\left(\mathbf{y}^0, \mathbf{V}^0, \mathbf{y}^1, \mathbf{V}^1, \dots, \mathbf{y}^N, \mathbf{V}^N, \mathbf{a}\right) = \left(\mathbf{Y}^0, \mathbf{P}^0, \mathbf{A}^0, \mathbf{P}\_0, \mathbf{Y}^1, \mathbf{P}^1, \dots, \mathbf{y}^N, \mathbf{P}^N, \mathbf{A}\right),
$$

$$
\Phi\_k\left(\mathbf{z}^0, \mathbf{W}^0, \mathbf{z}^1, \mathbf{W}^1, \dots, \mathbf{z}^N, \mathbf{W}^N, \mathbf{b}\right) = \left(\mathbf{Z}^0, \mathbf{R}^0, \mathbf{B}^0, \mathbf{R}\_0, \mathbf{Z}^1, \mathbf{R}^1, \dots, \mathbf{z}^N, \mathbf{R}^N, \mathbf{B}\right).
$$

$$
\left(\mathbf{y}^0, \mathbf{V}^0, \mathbf{y}^1, \mathbf{V}^1, \dots, \mathbf{y}^N, \mathbf{V}^N, \mathbf{a}\right), \left(\mathbf{z}^0, \mathbf{W}^0, \mathbf{z}^1, \mathbf{W}^1, \dots, \mathbf{z}^N, \mathbf{W}^N, \mathbf{b}\right) \in \mathcal{B}\_{\mathcal{A}\_k}(\bar{\mathbf{u}}\_k, \mathbb{R}\_k). \text{ Now, we set}
$$

$$\mathbf{E}^n = \mathbf{V}^n - \mathbf{W}^n \in \mathbb{R}^{l+n+1}, \mathbf{A}^n = \mathbf{y}^n - \mathbf{z}^n \in \mathbb{R}^{l+n}, \mathbf{0} \le n \le N; \sigma^n = b^n - a^n \in \mathbb{R}\_r - L \le n \le N\_r$$

$$\mathbb{E}^{\mathsf{u},\*} = \mathsf{V}^{\mathsf{u},\*} - \mathsf{W}^{\mathsf{u},\*} \in \mathbb{R}^{l+n+1}, \\ \Delta^{\mathsf{u},\*} = \mathsf{y}^{\mathsf{u},\*} - \mathsf{z}^{\mathsf{u},\*} \in \mathbb{R}^{l+n}, \\ \sigma^{\mathsf{u},\*} = b^{\mathsf{u},\*} - a^{\mathsf{u},\*} \in \mathbb{R}\_{\succ}$$

1 ≤ *n* ≤ *N*. From (22) and hypothesis (H6), we obtain

$$\|\Delta^{n,\*}\|\_{\infty} \le (1+\mathcal{C}k) \left\|\Delta^{n-1}\right\|\_{\infty} + \mathcal{C}k \left|\sigma^{n-1}\right|,\tag{37}$$

1 ≤ *n* ≤ *N*. Next, from (25), by means of hypothesis (H7) and inequality (35), we arrive at

$$\begin{split} |\boldsymbol{\sigma}^{n,\*}| &\leq \quad |\boldsymbol{\sigma}^{n-1}| + k \left| f\left(\boldsymbol{a}^{n-1}, \boldsymbol{a}^{n-1-L}, \boldsymbol{Q}\left(\mathbf{y}^{n-1}, \boldsymbol{\gamma}^{n-1}(\mathbf{y}, \boldsymbol{a}) \cdot \mathbf{V}^{n-1}\right), \boldsymbol{t}^{n-1}\right) \right. \\ &\quad - f\left(\boldsymbol{b}^{n-1}, \boldsymbol{b}^{n-1-L}, \boldsymbol{Q}\left(\mathbf{z}^{n-1}, \boldsymbol{\gamma}^{n-1}(\mathbf{z}, \boldsymbol{b}) \cdot \mathbf{W}^{n-1}\right), \boldsymbol{t}^{n-1}\right) \right| \\ &\leq \quad \left(1 + \mathbb{C}k\right) |\boldsymbol{\sigma}^{n-1}| + \mathbb{C}k \left| \boldsymbol{\sigma}^{n-1-L}| + \mathbb{C}k \left\{ \left\| \boldsymbol{\Delta}^{n-1} \right\|\_{\infty} + \left\| \mathbf{E}^{n-1} \right\|\_{1} \right\}, \end{split} \tag{38}$$

<sup>1</sup> <sup>≤</sup> *<sup>n</sup>* <sup>≤</sup> *<sup>N</sup>*. Now, from (23), hypotheses (H4), (H6) and the boundedness of **W***n*∞, we have, for *<sup>k</sup>* sufficiently small,

$$\left| \left| E\_{\rangle}^{\text{H},\*} \right| \le \left( 1 + \mathcal{C}k \right) \left\| \mathbf{E}^{n-1} \right\|\_{1} + \mathcal{C}k \left\{ \left| \sigma^{n-1} \right| + \left| \left| \Delta^{n-1} \right| \right|\_{\infty} \right\}\_{\omega} \tag{39}$$

1 ≤ *j* ≤ *J* + *n*. Additionally, (24), hypothesis (H6) and inequalities (36)–(38), allow us to obtain, for *k* sufficiently small,

$$\mathbb{P}\left|E\_0^{n,\*}\right| \le \mathbb{C}\left\{ \left| \sigma^{n-1} \right| + \left| \sigma^{n-1-L} \right| + \left\| \Delta^{n-1} \right\|\_{\infty} + \left\| \mathbf{E}^{n-1} \right\|\_1 + \left\| \mathbf{E}^{n,\*} \right\|\_1 \right\},\tag{40}$$

<sup>1</sup> <sup>≤</sup> *<sup>n</sup>* <sup>≤</sup> *<sup>N</sup>*. Then, we multiply <sup>|</sup>*En*,<sup>∗</sup> *<sup>j</sup>* <sup>|</sup> by *<sup>h</sup>* and sum in *<sup>j</sup>*, 0 <sup>≤</sup> *<sup>j</sup>* <sup>≤</sup> *<sup>J</sup>* <sup>+</sup> *<sup>n</sup>*, to achieve a bound of *En*,∗1, 1 ≤ *n* ≤ *N*. Therefore, from (39), (40) and that *k* = *r h*, we derive, for *h* sufficiently small,

$$\begin{split} \|\boldsymbol{E}^{n,\*}\|\_{1} &\leq \mathcal{C}h\left\{ \left|\boldsymbol{\sigma}^{n-1}\right| + \left|\boldsymbol{\sigma}^{n-1-L}\right| + \left\|\boldsymbol{\Delta}^{n-1}\right\|\_{\infty} + \left\|\boldsymbol{\mathbf{E}}^{n-1}\right\|\_{1} + \left\|\boldsymbol{\mathbf{E}}^{n,\*}\right\|\_{1} \right\} \\ &+ h\sum\_{j=1}^{l+n} \left\{ (1+\mathcal{C}k)\left\|\boldsymbol{\mathbf{E}}^{n-1}\right\|\_{1} + \mathcal{C}k\left( \left|\boldsymbol{\sigma}^{n-1}\right| + \left\|\boldsymbol{\Delta}^{n-1}\right\|\_{\infty} \right) \right\} \\ &\leq \mathcal{C}h\left\|\boldsymbol{\mathbf{E}}^{n,\*}\right\|\_{1} + \mathcal{C}\left\{ \left|\boldsymbol{\sigma}^{n-1}\right| + \left|\boldsymbol{\sigma}^{n-1-L}\right| + \left\|\boldsymbol{\Delta}^{n-1}\right\|\_{\infty} + \left\|\boldsymbol{\mathbf{E}}^{n-1}\right\|\_{1} \right\}, \end{split} \tag{41}$$

1 ≤ *n* ≤ *N*, and for *h* sufficiently small,

$$\|E^{n,\*}\|\_1 \le \mathcal{C} \left\{ \left| \sigma^{n-1} \right| + \left| \sigma^{n-1-L} \right| + \left\| \Delta^{n-1} \right\|\_\infty + \left\| \mathbf{E}^{n-1} \right\|\_1 \right\},\tag{42}$$

1 ≤ *n* ≤ *N*. Next, (19), hypotheses (H6)–(H7), (37) and (38) enable us to arrive at

$$\begin{split} \left| \left| \Delta\_{j}^{n} \right| \right| &\leq \left( 1 + \mathsf{C}k \right) \left| \Delta\_{j-1}^{n-1} \right| + \mathsf{C}k \left\{ \left| \sigma^{n-1} \right| + \left| \Delta\_{j}^{n,\*} \right| + \left| \sigma^{n,\*} \right| \right\} + k \left| Y\_{j}^{n} - Z\_{j}^{n} \right| \\ &\leq \left( 1 + \mathsf{C}k \right) \left| \Delta\_{j-1}^{n-1} \right| + \mathsf{C}k \left\{ \left| \left| \Delta^{n-1} \right| \right|\_{\infty} + \left| \left| \mathsf{E}^{n-1} \right| \right|\_{1} + \left| \sigma^{n-1} \right| + \left| \sigma^{n-1-L} \right| \right\} + k \left| Y\_{j}^{n} - Z\_{j}^{n} \right|, \end{split} \tag{43}$$

1 ≤ *j* ≤ *J* + *n*, 1 ≤ *n* ≤ *N*. Thus, from (43), when *N* ≥ *n* > *j* ≥ 1, we have

$$\begin{split} |\Lambda\_{j}^{\mathfrak{n}}| &\leq \quad \mathbb{C}k \sum\_{l=0}^{j-1} (1+\mathbb{C}k)^{l} \left( ||\mathbb{E}^{\mathfrak{n}-1-l}||\_{1} + ||\mathbb{A}^{\mathfrak{n}-1-l}||\_{\infty} + \left| \sigma^{\mathfrak{n}-1-l} \right| + \left| \sigma^{\mathfrak{n}-1-l-L} \right| \right) \\ &+ k \sum\_{l=0}^{j-1} (1+\mathbb{C}k)^{l} |Y\_{j-l}^{\mathfrak{n}-l} - Z\_{j-l}^{\mathfrak{n}-l}| . \end{split} \tag{44}$$

and when *J* + *n* ≥ *j* ≥ *n* ≥ 1, it follows that

$$\begin{split} |\Lambda\_{j}^{n}| \leq & \quad (1+\mathsf{C}k)^{n} \left| \Delta\_{j-n}^{0} \right| + \mathsf{C}k \sum\_{l=0}^{n-1} (1+\mathsf{C}k)^{l} \left( \left\| \mathsf{E}^{n-1-l} \right\|\_{1} + \left\| \mathsf{A}^{n-1-l} \right\|\_{\infty} + \left| \sigma^{n-1-l} \right| + \left| \sigma^{n-1-l-L} \right| \right) \\ & + k \sum\_{l=0}^{n-1} (1+\mathsf{C}k)^{l} \left| Y\_{j-l}^{n-l} - Z\_{j-l}^{n-l} \right|. \end{split} \tag{45}$$

Then, by means of (44) and (45), we can conclude that

$$\|\|\Delta^{\mathfrak{m}}\|\|\_{\infty} \leq \mathbb{C} \left\{ \|\Delta^{\mathfrak{0}}\|\|\_{\infty} + \sum\_{m=0}^{n-1} k \, \|\mathbf{E}^{\mathfrak{m}}\|\|\_{1} + \sum\_{m=0}^{n-1} k \, \|\Delta^{\mathfrak{m}}\|\|\_{\infty} + \sum\_{m=-L}^{n-1} k \, |\sigma^{\mathfrak{m}}| + \sum\_{m=1}^{n} k \, \|\mathbf{Y}^{\mathfrak{m}} - \mathbf{Z}^{\mathfrak{m}}\|\_{\infty} \right\}, \tag{46}$$

1 ≤ *n* ≤ *N*. On the other hand, from (21), (H7) and (35)–(38) and (42), for *k* small enough, we obtain that

<sup>|</sup>*σn*| ≤ (<sup>1</sup> <sup>+</sup> *C k*)|*σn*−1<sup>|</sup> <sup>+</sup> *C k* \$ <sup>|</sup>*σn*,∗| <sup>+</sup> *σn*−*<sup>L</sup>* <sup>+</sup> |Q(**y***n*,∗, *<sup>γ</sup>n*,∗(**y**, *<sup>a</sup>*) · **<sup>V</sup>***n*,∗) − Q(**z***n*,∗, *<sup>γ</sup>n*,∗(**z**, *<sup>b</sup>*) · **<sup>W</sup>***n*,∗)<sup>|</sup> + Q(**y***n*−1, *<sup>γ</sup>n*−1(**y**, *<sup>a</sup>*) · **<sup>V</sup>***n*−1) − Q(**z***n*−1, *<sup>γ</sup>n*−1(**z**, *<sup>b</sup>*) · **<sup>W</sup>***n*−1) + *σn*−*L*−<sup>1</sup> % <sup>+</sup> *<sup>k</sup>* <sup>|</sup>*A<sup>n</sup>* <sup>−</sup> *<sup>B</sup>n*<sup>|</sup> <sup>≤</sup> (<sup>1</sup> <sup>+</sup> *C k*)|*σn*−1<sup>|</sup> <sup>+</sup> *C k* \$) ) ) **Δ***n*−<sup>1</sup> ) ) ) <sup>∞</sup> <sup>+</sup> ) ) ) **E***n*−<sup>1</sup> ) ) ) 1 % <sup>+</sup> *C k* \$ **Δ***n*,∗<sup>∞</sup> <sup>+</sup> **E***n*,∗<sup>1</sup> <sup>+</sup> <sup>|</sup>*σn*,∗| <sup>+</sup> *σn*−*L*−<sup>1</sup> + *σn*−*<sup>L</sup>* % <sup>+</sup> *<sup>k</sup>* <sup>|</sup>*A<sup>n</sup>* <sup>−</sup> *<sup>B</sup>n*<sup>|</sup> <sup>≤</sup> (<sup>1</sup> <sup>+</sup> *C k*)|*σn*−1<sup>|</sup> <sup>+</sup> *C k* \$ *σn*−*L*−<sup>1</sup> + *σn*−*<sup>L</sup>* + ) ) ) **Δ***n*−<sup>1</sup> ) ) ) <sup>∞</sup> <sup>+</sup> ) ) ) **E***n*−<sup>1</sup> ) ) ) 1 % <sup>+</sup> *<sup>k</sup>* <sup>|</sup>*A<sup>n</sup>* <sup>−</sup> *<sup>B</sup>n*|, (47)

1 ≤ *n* ≤ *N*. Thus,

$$\begin{aligned} \left|\sigma^{n}\right| &\leq \left(1+\mathcal{C}k\right)^{n}\left|\sigma^{0}\right| + k\sum\_{l=0}^{n-1} (1+\mathcal{C}k)^{l}\left|A^{n-l} - B^{n-l}\right| \\ &+ \mathcal{C}k\sum\_{l=0}^{n-1} (1+\mathcal{C}k)^{l}\left\{ \left|\sigma^{n-l-L-1}\right| + \left|\sigma^{n-l-L}\right| + \left\|\Delta^{n-l-1}\right\|\_{\infty} + \left\|\mathbf{E}^{n-l-1}\right\|\_{1} \right\}, \end{aligned} \tag{48}$$

1 ≤ *n* ≤ *N*. Therefore,

$$\mathbb{P}\left|\sigma^{\mathfrak{n}}\right| \leq \mathbb{C}\left\{ \left|\sigma^{0}\right| + \sum\_{m=-L}^{n-L} k \left|\sigma^{m}\right| + \sum\_{m=1}^{n} k \left|A^{m} - B^{m}\right| + \sum\_{m=0}^{n-1} k \left\|\Delta^{\mathfrak{n}}\right\|\_{\infty} + \sum\_{m=0}^{n-1} k \left\|\mathbf{E}^{m}\right\|\_{1} \right\},\tag{49}$$

1 ≤ *n* ≤ *N*. Finally, when we study the residuals caused by the approximation to the solution, from (20), hypotheses (H4), (H6), (H7), the inequalities (35), (37), (38) and **W***n*−1<sup>∞</sup> <sup>≤</sup> *<sup>C</sup>*, we have, for *k* sufficiently small,

$$\begin{split} |E\_{j}^{n}| \leq & k \left| P\_{j}^{n} - R\_{j}^{n} \right| + (1 + \mathbb{C}k) \left| E\_{j-1}^{n-1} \right| + \mathbb{C}k \left\{ \left| \Delta\_{j-1}^{n-1} \right| + \left| \sigma^{n-1} \right| + \left| \Delta\_{j}^{n,\*} \right| + \left| \sigma^{n,\*} \right| \right\} \\ \leq & k \left| P\_{j}^{n} - R\_{j}^{n} \right| + (1 + \mathbb{C}k) \left| E\_{j-1}^{n-1} \right| + \mathbb{C}k \left\{ \left\| \Delta^{n-1} \right\|\_{\infty} + \left\| \mathbf{E}^{n-1} \right\|\_{1} + \left| \sigma^{n-1} \right| + \left| \sigma^{n-1-L} \right| \right\}, \end{split} \tag{50}$$

1 ≤ *j* ≤ *J* + *n*, 1 ≤ *n* ≤ *N*. Now, we derive the stability estimate for the boundary node from (18). We employ hypotheses (H5) and (H6), the properties (P3) and (P6), and that **W***n*<sup>∞</sup> <sup>≤</sup> *<sup>C</sup>*, to achieve

$$|E\_0^n| \le |P\_0^n - R\_0^n| + \mathcal{C} \left\{ |\sigma^n| + ||\mathbf{A}^n||\_\infty + ||\mathbf{E}^n||\_1 \right\},\tag{51}$$

1 ≤ *n* ≤ *N*. Thus, from (50), when *N* ≥ *n* > *j* ≥ 1, we obtain

$$\begin{split} |E\_j^n| &\le (1+\mathcal{C}k)^j \left| E\_0^{n-j} \right| + k \sum\_{l=0}^{j-1} (1+\mathcal{C}k)^l \left| P\_{j-l}^{n-l} - R\_{j-l}^{n-l} \right| \\ &+ \mathcal{C}k \sum\_{l=0}^{j-1} (1+\mathcal{C}k)^l \left\{ \left\| \mathbf{E}^{n-1-l} \right\|\_{1} + \left\| \mathbf{A}^{n-1-l} \right\|\_{\infty} + |\sigma^{n-1-l}| + |\sigma^{n-1-l-L}| \right\}. \end{split} \tag{52}$$

Additionally, when *J* + *n* ≥ *j* ≥ *n* ≥ 1, it follows that

$$\begin{split} |E\_j^n| &\le \left(1 + \mathbb{C}k\right)^n |E\_{j-n}^0| + k \sum\_{l=0}^{n-1} (1 + \mathbb{C}k)^l \left| P\_{j-l}^{n-l} - R\_{j-l}^{n-l} \right| \\ &+ \mathbb{C}k \sum\_{l=0}^{n-1} (1 + \mathbb{C}k)^l \left\{ \| \mathbf{E}^{n-1-l} \| \_1 + \| \mathbf{A}^{n-1-l} \| \_\infty + |\sigma^{n-1-l}| + |\sigma^{n-1-l-L}| \right\}. \end{split} \tag{53}$$

Thus, we can conclude that

$$\mathbb{E}\left|E\_{\boldsymbol{\beta}}^{\boldsymbol{m}}\right| \leq \mathbb{C}\left\{ \left\|\mathbf{E}^{0}\right\|\_{1} + \sum\_{m=0}^{n-1} k \left\|\mathbf{E}^{m}\right\|\_{1} + \sum\_{m=0}^{n-1} k \left\|\mathbf{A}^{m}\right\|\_{\infty} + \sum\_{m=-L}^{n-1} k \left|\sigma^{m}\right| + \sum\_{m=1}^{n} k \left\|\mathbf{P}^{m} - \mathbf{R}^{m}\right\|\_{\infty} \right\}.\tag{54}$$

As we are interested in a maximum norm bound of **E***n*, next, we proceed to achieve a bound for the term **E***n*1, 1 <sup>≤</sup> *<sup>n</sup>* <sup>≤</sup> *<sup>N</sup>*. Then, for *<sup>k</sup>* small enough,

$$\|\mathbb{E}^{\pi}\|\_{1} \leq \mathbb{C}\left\{\|\mathbb{E}^{0}\|\_{1} + \sum\_{m=1}^{\pi}k\,\|\mathbb{E}^{m}\|\_{1} + \sum\_{m=0}^{\pi}k\,\|\Delta^{m}\|\_{\infty} + \sum\_{m=-L}^{\pi}k\,|\sigma^{m}| + \sum\_{m=1}^{\pi}k\,\|\mathbb{P}^{m} - \mathbb{R}^{m}\|\_{\infty} + \|\mathbb{P}\_{0} - \mathbb{R}\_{0}\|\_{\infty}\right\},\tag{55}$$

1 ≤ *n* ≤ *N*. Thus, by means of the discrete Gronwall lemma,

$$\|\|\mathbf{E}^{\mathbf{n}}\|\|\_{1} \leq \mathbb{C} \left\{ \|\mathbf{E}^{0}\|\|\_{1} + \sum\_{m=0}^{n} k \, \|\Delta^{m}\|\|\_{\infty} + \sum\_{m=-L}^{n} k \, |\sigma^{m}| + \sum\_{m=1}^{n} k \, \|\mathbf{P}^{m} - \mathbf{R}^{m}\|\|\_{\infty} + \|\mathbf{P}\_{0} - \mathbf{R}\_{0}\|\|\_{\infty} \right\},\tag{56}$$

1 ≤ *n* ≤ *N*. Next, we substitute (56) into (46), and by means of the discrete Gronwall Lemma, for *k* sufficiently small, it follows that

$$\|\|\boldsymbol{\Delta}^{\boldsymbol{\mathsf{w}}}\|\|\_{\infty} \leq \mathbb{C} \left\{ \|\boldsymbol{\Delta}^{0}\|\|\_{\infty} + \|\mathbb{E}^{0}\|\boldsymbol{\varepsilon} + \sum\_{m=-L}^{n-1} k \|\boldsymbol{\varepsilon}^{m}\| + \|\mathbb{P}\_{0} - \mathbb{R}\_{\boldsymbol{0}}\|\_{\infty} + \sum\_{m=1}^{n-1} k \,\|\mathbb{P}^{m} - \mathbb{R}^{m}\|\|\_{\infty} + \sum\_{m=1}^{n} k \,\|\mathbb{Y}^{m} - \mathbb{Z}^{m}\|\_{\infty} \right\},\tag{57}$$

1 ≤ *n* ≤ *N*. Now, we substitute (56) and (57) in (49), and apply again the discrete Gronwall Lemma, to attain

$$\begin{split} |\sigma^{n}| \quad \leq \quad \mathbb{C} \left\{ \||\mathbf{A}^{0}\|\!\_{\infty} + \||\mathbf{E}^{0}\|\!\_{1} + \|\sigma^{0}\|\!\_{\infty} + \sum\_{m=1}^{n} k \left| A^{m} - B^{m} \right| + \||\mathbf{P}\_{0} - \mathbf{R}\_{0}\|\!\_{\infty} \right. \\ \left. + \sum\_{m=1}^{n-1} k \left\|\mathbf{P}^{m} - \mathbf{R}^{m}\right\|\!\_{\infty} + \sum\_{m=1}^{n-1} k \left\|\mathbf{Y}^{m} - \mathbf{Z}^{m}\right\|\!\_{\infty} \right\}, \end{split} \tag{58}$$

1 ≤ *n* ≤ *N*. Then, we substitute (58) in (56) and (57) to obtain

$$\begin{split} \|\mathbf{A}^{\mathrm{m}}\|\_{\infty} &\leq \quad \mathbb{C} \left\{ \|\mathbf{A}^{\mathrm{0}}\|\_{\infty} + \|\mathbf{E}^{\mathrm{0}}\|\_{1} + \|\sigma^{\mathrm{0}}\|\_{\infty} + \|\mathbf{P}\_{0} - \mathbf{R}\_{0}\|\_{\infty} + \sum\_{m=1}^{n-1} k \, |A^{\mathrm{m}} - B^{\mathrm{m}}| \\ &+ \sum\_{m=1}^{n-1} k \, \|\mathbf{P}^{\mathrm{m}} - \mathbf{R}^{\mathrm{m}}\|\_{\infty} + \sum\_{m=1}^{n} k \, \|\mathbf{Y}^{\mathrm{m}} - \mathbf{Z}^{\mathrm{m}}\|\_{\infty} \right\}, \end{split} \tag{59}$$

and

$$\begin{split} \|\mathbf{E}^{\mathbf{u}}\|\_{1} &\leq \quad \mathbb{C} \left\{ \|\mathbf{A}^{0}\|\|\_{\infty} + \|\mathbf{E}^{0}\|\_{1} + \|\sigma^{0}\|\_{\infty} + \|\mathbf{P}\_{0} - \mathbf{R}\_{0}\|\|\_{\infty} + \sum\_{m=1}^{n} k \, |A^{m} - B^{m}| \\ &+ \sum\_{m=1}^{n} k \, \|\mathbf{P}^{m} - \mathbf{R}^{m}\|\_{\infty} + \sum\_{m=1}^{n} k \, \|\mathbf{Y}^{m} - \mathbf{Z}^{m}\|\_{\infty} \right\} \,. \end{split} \tag{60}$$

1 ≤ *n* ≤ *N*. Finally, we substitute (58)–(60) in (51) and (54) to arrive at

$$\begin{split} \|\|\mathbf{E}^{\rm H}\|\|\_{\infty} &\leq \quad \mathbb{C} \left\{ \|\mathbf{A}^{0}\|\|\_{\infty} + \|\|\mathbf{E}^{0}\|\|\_{1} + \|\sigma^{0}\|\_{\infty} + \|\mathbf{P}\_{0} - \mathbf{R}\_{0}\|\|\_{\infty} + \sum\_{m=1}^{n} k \, |A^{m} - B^{\rm H}| \\ &+ \sum\_{m=1}^{n} k \, \|\mathbf{P}^{m} - \mathbf{R}^{m}\|\_{\infty} + \sum\_{m=1}^{n} k \, \|\mathbf{Y}^{m} - \mathbf{Z}^{m}\|\|\_{\infty} \right\} \, \end{split} \tag{61}$$

1 ≤ *n* ≤ *N*. Thus, due to (58), (59) and (61) we have

$$\begin{split} & \| \left( \Delta^{0}, \mathbf{E}^{0}, \dots, \Delta^{N}, \mathbf{E}^{N}, \sigma \right) \|\_{\mathcal{A}\_{k}} \\ & \leq \quad \mathbb{C} \, \| \left( \Delta^{0}, \mathbf{E}^{0}, \sigma^{0}, \mathbf{P}\_{0} - \mathbf{R}\_{0}, \mathbf{Y}^{1} - \mathbf{Z}^{1}, \mathbf{P}^{1} - \mathbf{R}^{1}, \dots, \mathbf{Y}^{N} - \mathbf{Z}^{N}, \mathbf{P}^{N} - \mathbf{R}^{N}, \mathbf{A} - \mathbf{B} \right) \|\_{\mathcal{B}\_{k}}. \end{split}$$

We finish the analysis of the numerical method with the convergence. The global discretization error is defined as

$$
\bar{\mathbf{e}}\_k = \mathbf{\ddot{u}}\_k - \mathbf{\ddot{U}}\_k \in \mathcal{A}\_k.
$$

We say that the discretization (26) is convergent if there exists *k*<sup>0</sup> ∈ *K* such that, for each *k* ∈ *K* with *<sup>k</sup>* <sup>≤</sup> *<sup>k</sup>*0, (26) has a solution **<sup>U</sup>**˜ *<sup>k</sup>* for which, as *<sup>k</sup>* <sup>→</sup> 0,

$$\lim \|\|\mathbf{\bar{u}}\_k - \mathbf{\bar{O}}\_k\|\|\_{\mathcal{A}\_k} = \lim \|\|\mathbf{\bar{e}}\_k\|\|\_{\mathcal{A}\_k} = 0.$$

In our analysis, we use the following result of the general discretization framework introduced in [14].

**Theorem 4.** *Assume that (26) is consistent and stable with thresholds Rk. If* Φ*<sup>k</sup> is continuous in B*(**u**˜ *<sup>k</sup>*, *Rk*) *and* **l***k*B*<sup>k</sup>* = *<sup>o</sup>*(*Rk*) *as k* → <sup>0</sup>*, then:*

*i) For k sufficiently small, the discrete Equations (26) possess a unique solution in B*(**u**˜ *<sup>k</sup>*, *Rk*)*.*

*ii) As k* → <sup>0</sup>*,* **e**˜ *<sup>k</sup>*A*<sup>k</sup>* = *<sup>O</sup>*(**l***k*B*<sup>k</sup>* )*.*

Finally, we propose the next theorem which establishes the *convergence* of the numerical method defined by Equations (26).

**Theorem 5.** *Assume that hypotheses (H1)–(H7) hold and that the considered quadrature rules satisfy properties (P1)–(P6). Then, for k sufficiently small, the numerical method defined by Equations (26) has a unique solution* **<sup>U</sup>**˜ *<sup>k</sup>* <sup>∈</sup> *<sup>B</sup>*(**u**˜ *<sup>k</sup>*, *Rk*) *and*

$$\|\|\mathbf{U}\_{k} - \bar{\mathbf{u}}\_{k}\|\|\_{\mathcal{A}\_{k}} \le \mathcal{C} \left( \|\mathbf{x}^{0} - \mathbf{X}^{0}\|\_{\infty} + \|\mathbf{u}^{0} - \mathbf{U}^{0}\|\_{\infty} + \|\mathbf{s}^{0} - \mathbf{S}^{0}\|\_{\infty} + O(h^{2} + k^{2}) \right). \tag{62}$$

The proof of Theorem 5 is immediately derived by means of consistency (Theorem 2), stability (Theorem 3) and Theorem 4. Specifically, if **X**<sup>0</sup> = **x**0, **U**<sup>0</sup> = **u**<sup>0</sup> and **S**<sup>0</sup> = **s**0, the proposed numerical scheme is second-order accurate.

Next, we also establish an error bound for the differences between the theoretical solution *u* evaluated at the numerical values of the grid nodes, and the approximation obtained with the numerical method.

**Theorem 6.** *Assume that hypotheses (H1)–(H7) hold and that the considered quadrature rules satisfy properties (P1)–(P6). Then, for k sufficiently small, let*

$$\mathbf{u}\_k^\dagger = \left(\mathbf{u}^{0,\dagger}, \mathbf{u}^{1,\dagger}, \dots, \mathbf{u}^{N,\dagger}\right) \in \prod\_{n=0}^N \mathbb{R}^{J+n}\text{'} $$

*be defined by* **u***n*,† = *u*(*X<sup>n</sup>* <sup>0</sup> , *t <sup>n</sup>*), *u*(*X<sup>n</sup>* <sup>1</sup> , *t <sup>n</sup>*),..., *u*(*X<sup>n</sup> <sup>J</sup>*+*n*, *t n*) <sup>∈</sup> <sup>R</sup>*J*<sup>+</sup>*n,* <sup>0</sup> <sup>≤</sup> *<sup>n</sup>* <sup>≤</sup> *N, where <sup>X</sup><sup>n</sup> <sup>j</sup> ,* 0 ≤ *j* ≤ *J* + *n,* <sup>0</sup> <sup>≤</sup> *<sup>n</sup>* <sup>≤</sup> *N, are the grid nodes given in (26). Then, the solution* **<sup>U</sup>**˜ *<sup>k</sup> satisfies*

$$\max\_{0 \le n \le N} \left\{ \max \left\{ \| \mathbf{u}^{n, \dagger} - \mathbf{U}^{0} \|\_{\infty} |s(t^{\mathbf{u}}) - \mathbf{S}^{\mathbf{u}}| \right\} \right\} \le C \left( \| \mathbf{x}^{0} - \mathbf{X}^{0} \|\_{\infty} + \| \mathbf{u}^{0} - \mathbf{U}^{0} \|\_{\infty} + \| \mathbf{s}^{0} - \mathbf{S}^{0} \|\_{\infty} + O(\hbar^{2} + k^{2}) \right). \tag{63}$$

**Proof.** We only have to consider the triangle inequality

$$|u(X\_{\mathbf{j}}^{n},t^{n}) - \mathcal{U}\_{\mathbf{j}}^{n}| \leq |u(X\_{\mathbf{j}}^{n},t^{n}) - u(\mathbf{x}\_{\mathbf{j}}^{n},t^{n})| + |u(\mathbf{x}\_{\mathbf{j}}^{n},t^{n}) - \mathcal{U}\_{\mathbf{j}}^{n}|.$$

0 ≤ *j* ≤ *J* + *n*, 0 ≤ *n* ≤ *N*. The smoothness hypothesis (H1) on *u* is enough to derive that the first term on the right hand side of the inequality is *<sup>O</sup>*(|*X<sup>n</sup> <sup>j</sup>* <sup>−</sup> *<sup>x</sup><sup>n</sup> <sup>j</sup>* <sup>|</sup>) = *<sup>O</sup>*(*h*<sup>2</sup> <sup>+</sup> *<sup>k</sup>*2), as proved in Theorem 5. Additionally, the second-order rate of convergence proved in such a Theorem shows that the second term is *O*(*h*<sup>2</sup> + *k*2).

The last convergence theorem remains true even if the method employs a cuadrature rule at each time-step in which only a selection of nodes are involved, whenever the subgrid of these cuadrature nodes satisfies the (SG) property. In the case of the selection procedure given by (10), the convergence can be demonstrated in two stages. First, as proven in [12], we can establish that the selection procedure (10) creates subgrids with the property (SG) if the procedure is applied to a grid whose nodes are in the neighborhood of the theoretical ones within a radius *R kp*. This is just the case, as the convergence Theorem 6 establishes, step by step, if the solutions at each fixed time level are obtained by the discrete operator (14).

#### **4. Numerical Results**

We carried out different simulations with the aim of exploring the behavior of the numerical method for the structured consumer population with a delayed resource model (1), (3). We tried to corroborate the convergence theoretical results with an academical test. It consisted of a test problem with meaningful nonlinearities both from mathematical and biological points of view in which the delay differential equation for the resource employed *τ* = 1. We used the following size-specific growth, fertility and mortality rates

$$g(x,z,t) = \frac{\lambda}{2} \frac{1+z}{z} \left( \left(\frac{z}{1+z}\right)^2 - x^2 \right) + \frac{x\rho}{1+z} \left(\frac{29}{30} - \frac{z}{\kappa}\right),$$

$$a(x,z,t) = \frac{3}{2}\lambda \frac{1 + \left(\frac{z}{\Lambda\left(\frac{29}{30} - \frac{z}{\kappa}\right)}\right)^{\frac{-30\lambda}{29\rho}}}{1 + 2\left(\frac{z}{\Lambda\left(\frac{29}{30} - \frac{z}{\kappa}\right)}\right)^{\frac{-30\lambda}{29\rho}}},$$

$$\mu(x,z,t) = \lambda \frac{1+z}{z} \left(\frac{z}{1+z} + 2x\right) - \frac{3\rho}{1+z} \left(\frac{29}{30} - \frac{z}{\kappa}\right).$$

The weight function in (4) was *γ*(*x*, *z*, *t*) = *x*2, and finally, the function that drove the dynamics of the resource (3) was chosen to be

$$f(z(t), z(t-1), \zeta, t) = \rho \, z(t) \left( 1 - \frac{z(t-1)}{\kappa} \right) - \rho \, \zeta \, \left( \kappa + 30 \left( z(t) - z(t-1) \right) \right) \frac{(1 + z(t))^5}{\kappa \left( z(t) \right)^4 (1 + 4\epsilon^{-\lambda})}.$$

With this choice of data functions, the problem (1), (3) has the following solution:

$$u(\mathbf{x},t) = \left(\frac{s(t)}{1+s(t)} - \mathbf{x}\right)^2 + e^{-\lambda t} \left(\left(\frac{s(t)}{1+s(t)}\right)^2 - \mathbf{x}^2\right),$$

$$s(t) = \frac{29}{30} \frac{\Lambda e^{29\rho t/30}}{1 + \Lambda e^{29\rho t/30}/\kappa}.$$

where Λ = 24; *κ* = 5; and *ρ* and *λ* are parameters that must be fixed. In the experiment, we employed *ρ* = 0.1 and *λ* = 0.3.

The knowledge of the exact solution to the problem allowed us to compare it with the numerical solution and to compute the error caused by the numerical approximation. Then, given *h* and *k* the discretization parameters in size and time, and the corresponding numerical solution with (**X**0, **U**0, **X**1, **U**1,..., **X***N*, **U***N*, **S**), we computed

$$\mathcal{A}\_{h,k}^{\dagger} = \max \left\{ \max\_{0 \le n \le N} \left\{ \max\_{0 \le j \le J} |u(X\_j^n, t^n) - \mathcal{U}\_j^n| \right\}\_{0 \le n \le N} \left\{ |s(t^n) - S^n| \right\} \right\}. \tag{64}$$

Note that the exact positions of the grid nodes, *x<sup>n</sup> <sup>j</sup>* , 0 ≤ *j* ≤ *J*, 0 ≤ *n* ≤ *N*, are unknown, so we compared them with the solution on the computed grid *X<sup>n</sup> <sup>j</sup>* , 0 ≤ *j* ≤ *J*, 0 ≤ *n* ≤ *N*. In Figure 1, we show the evolution of the computed error in both consumer and resource populations. We observe that the error is produced mainly at the birth (minimum size) and increases with time.

We can also obtain the experimental order of convergence by means of the following well-known formula:

$$
\widehat{order\_{2\,h,2\,k}} = \frac{\log\left(\mathfrak{E}\_{2\,h,2\,k}^{\mathfrak{t}}/\mathfrak{E}\_{h\,k}^{\mathfrak{t}}\right)}{\log 2}.
$$

We can show the convergence of our method and that it is second-order accurate by means of Table 1. In the test, the initial size interval was taken as [0, *x*<sup>0</sup> *<sup>M</sup>*], with *<sup>x</sup>*<sup>0</sup> *<sup>M</sup>* = 0.8, and the numerical integration was carried out on the time interval [0, *T*], with *T* = 10. Each column and each row of the table represent a computation with different values of the time and size discretization parameters, respectively. The upper number of each entry in columns two to six of said table represents the error

computed from (64) and the lower quantity is the experimental order of convergence. We also remark that the ratio *r* := *k*/*h* was kept fixed as soon as we considered approximations along each diagonal of the table. We realize that the results in Table 1 clearly confirm the expected (theoretically) second-order convergence for these approximations. Different values of parameters *ρ* and *λ* confirm the same behavior of the error.

**Figure 1.** *<sup>k</sup>* = 1.25 <sup>×</sup> <sup>10</sup>−2, *<sup>h</sup>* = 1.25 <sup>×</sup> <sup>10</sup>−1.Plot on the left: evolution of the error in the consumer population; plot on the right: evolution of the error in the resource population.

**Table 1.** T = 10. For each value of *h* and *k*, the corresponding upper number is *e*ˆ † *<sup>h</sup>*,*k*; the lower number is *order* † *<sup>h</sup>*,*k*.


## **5. Conclusions**

We proposed a numerical method specially adapted to integrate a model that describes the evolution of a size-structured consumer population feeding on a dynamical resource. The model couples a first-order hyperbolic partial differential equation with nonlocal and nonlinear boundary condition and a differential equation driving the dynamics of the resource. The novelty of the model is the appearance of a delay in the resource evolution rate. Numerical methods are the only feasible approach, except in exceptional cases, for the approximation of the solution of the problem and eventually studying the dynamical behavior of the model. The authors have in preparation [15], using this model, a numerical study on the effect of the delay in the resource evolution rate on the dynamics of a population of the *Daphnia magna* (water flea).

The analyzed numerical method integrated the problem along its characteristic curves. We proved optimal second-order of convergence to the solution; this was confirmed numerically by means of an academic test problem. This is the first convergence analysis to our knowledge for the discretization of this kind of problem.

**Author Contributions:** All authors contributed equally to this article. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded in part by project MTM2017-85476-C2-1-P of the Spanish Ministerio de Economía y Competitividad and European FEDER Funds. VA138G18 of the Consejería de Educación, Junta de Castilla y León, Spain.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.
