**3. The Stability and Convergence of the Constructed Difference Schemes**

First, we will start with some technical lemmas, which will be helpful in the context of convergence and stability.

The nonlinear delay term ˜ *f*(*μ*, *ν*, *x*, *t*) is sufficiently smooth and it satisfies the Lipschitz condition

$$\left| \bar{f}(\mu\_1, \nu, \mathbf{x}, t) - \bar{f}(\mu\_2, \nu, \mathbf{x}, t) \right| \le c\_1 |\mu\_1 - \mu\_2|, \quad \forall \mu\_1, \ \mu\_2 \in [0, L] \times [0, T], \tag{33}$$

$$|\tilde{f}(\mu, \nu\_1, \mathbf{x}, t) - \tilde{f}(\mu, \nu\_2, \mathbf{x}, t)| \le c\_1 |\nu\_1 - \nu\_2|, \quad \forall \nu\_1, \ \nu\_2 \in [0, L] \times [-\mathfrak{r}, T],\tag{34}$$

where *c*<sup>1</sup> and *c*<sup>2</sup> are two positive constants.

**Lemma 3.1.** *For any* V*<sup>k</sup> <sup>i</sup>* , U*<sup>k</sup> <sup>i</sup>* <sup>∈</sup> <sup>W</sup><sup>ˆ</sup> *<sup>h</sup>*, *if we define <sup>ε</sup><sup>k</sup> <sup>i</sup>* = <sup>V</sup>*<sup>k</sup> <sup>i</sup>* <sup>−</sup> <sup>U</sup>*<sup>k</sup> <sup>i</sup>* , *and also*

$$\mathcal{F}\_{i}^{k+\sigma} = \bar{f}\left( (\sigma+1)\mathbf{U}\_{i}^{k} - \sigma \mathbf{U}\_{i}^{k-1}, \sigma \mathbf{U}\_{i}^{k+1-n\_{0}} + (1-\sigma)\mathbf{U}\_{i}^{k-n\_{0}}, \mathbf{x}\_{i}, t\_{k+\sigma} \right),\tag{35}$$

*the following estimate is satisfied*

$$\left\|\left\|\boldsymbol{\mathcal{F}}\_{i}^{k+\sigma}-\bar{\boldsymbol{\mathcal{F}}}\_{i}^{k+\sigma}\right\|\right\|^{2}\leq 4\left[2c\_{1}^{2}\left(\left\|\boldsymbol{\varepsilon}\_{i}^{k}\right\|^{2}+\left\|\boldsymbol{\varepsilon}\_{i}^{k-1}\right\|^{2}\right)+c\_{2}^{2}\left(\left\|\boldsymbol{\varepsilon}\_{i}^{k+1-n\_{0}}\right\|^{2}+\left\|\boldsymbol{\varepsilon}\_{i}^{k-n\_{0}}\right\|^{2}\right)\right].\tag{36}$$

**Proof.** Recalling *F*˜*k*+*<sup>σ</sup> <sup>i</sup>* from (24) and under the assumptions of the nonlinear delay term (33), we can deduce the following estimate

$$|\mathcal{F}\_{i}^{k+\sigma} - \tilde{\mathcal{F}}\_{i}^{k+\sigma}| \le c\_1 |(\sigma + 1)\varepsilon\_i^k - \sigma \varepsilon\_i^{k-1}| + c\_2 |\sigma \varepsilon\_i^{k+1-n\_0} + (1 - \sigma)\varepsilon\_i^{k-n\_0}|, \quad 1 \le i \le M - 1, \quad 1 \le i \le N$$

and so

$$\mathbb{E}\left\|\tilde{F}\_{i}^{k+\sigma} - \tilde{\mathcal{F}}\_{i}^{k+\sigma}\right\|^{2} \leq h \sum\_{i=1}^{M-1} \left(c\_{1}|(\sigma+1)\varepsilon\_{i}^{k} - \sigma\varepsilon\_{i}^{k-1}| + c\_{2}|\sigma\varepsilon\_{i}^{k+1-\eta\_{0}} + (1-\sigma)\varepsilon\_{i}^{k-\eta\_{0}}|\right)^{2} \tag{37}$$

$$\leq 2c\_1^2 h \sum\_{i=1}^{M-1} \left[ \left| (\sigma + 1)\varepsilon\_i^k - \sigma \varepsilon\_i^{k-1} \right| \right]^2 + 2c\_2^2 h \sum\_{i=1}^{M-1} \left[ \left| \sigma \varepsilon\_i^{k+1-n\_0} + (1-\sigma)\varepsilon\_i^{k-n\_0} \right| \right]^2 \tag{38}$$

$$\leq 4\left[2c\_1^2\left(\left\|\boldsymbol{\varepsilon}\_i^k\right\|^2 + \left\|\boldsymbol{\varepsilon}\_i^{k-1}\right\|^2\right) + c\_2^2\left(\left\|\boldsymbol{\varepsilon}\_i^{k+1-n\_0}\right\|^2 + \left\|\boldsymbol{\varepsilon}\_i^{k-n\_0}\right\|^2\right)\right].\tag{39}$$

For the convenience of our analysis, the following Grönwall inequality is recalled.

**Lemma 3.2.** *Suppose that* {*H<sup>k</sup>* <sup>|</sup> *<sup>k</sup>* <sup>≥</sup> <sup>0</sup>} *is a non negative sequence that satisfies*

$$H^{k+1} \le A + B\tau \sum\_{j=1}^k H^j, \quad k = 0, 1, \dots, \tau$$

*then*

$$H^{k+1} \le A \exp(Bk\tau), \quad k = 0, 1, \dots, \tau$$

*in which the positivity of the constants A and B must be taken into account.*

**Lemma 3.3.** *For any p<sup>k</sup> <sup>i</sup>* , *<sup>q</sup><sup>k</sup> <sup>i</sup>* <sup>∈</sup> <sup>W</sup><sup>ˆ</sup> *<sup>h</sup>*, *the following estimates are satisfied*

$$
\left\langle \delta\_l \left( \delta\_{\mathcal{X}} \left( \varkappa \varrho \delta\_{\mathcal{X}} p \right)\_i^{1/2} \right), -2\sigma p\_i^1 \right\rangle = \frac{2r}{\pi} \left[ \left| p\_i^1 \right|\_{1, \varkappa\_2}^2 - \left\langle p\_i^0, p\_i^1 \right\rangle\_{\varkappa\_2} \right]\_{\varkappa\_2} \tag{40}
$$

$$
\left\langle \delta\_{\mathbf{x}} \left( \varkappa \varrho \delta\_{\mathbf{x}} q \right)\_{i}^{1/2}, -2\sigma p\_{i}^{1} \right\rangle = \sigma \left[ \left\langle q\_{i}^{1}, p\_{i}^{1} \right\rangle\_{\varkappa\_{2}} + \left\langle q\_{i}^{0}, p\_{i}^{1} \right\rangle\_{\varkappa\_{2}} \right], \tag{41}
$$

$$\left\langle \partial\_{\bar{l}} \left( \delta\_{x} \left( \varkappa\_{2} \delta\_{x} p \right)\_{i}^{k+\sigma} \right), -p\_{i}^{k+\sigma} \right\rangle \geq \frac{1}{4\pi} \left( \tilde{\mathcal{E}}^{k+1} - \tilde{\mathcal{E}}^{k} \right),\tag{42}$$

$$\left\langle \delta\_{\times} (\varkappa\_2 \delta\_x q)\_i^{k+\sigma}, -p\_i^{k+\sigma} \right\rangle = \langle q\_i^{k+\sigma}, p\_i^{k+\sigma} \rangle\_{\varkappa\_2}.\tag{43}$$

*where*

$$\bar{\mathcal{E}}^{k+1} = (2\sigma + 1)|w^{k+1}|\_{1,\varkappa\_2}^2 - (2\sigma - 1)|w^k|\_{1,\varkappa\_2}^2 + (2\sigma^2 + \sigma - 1)|w^{k+1} - w^k|\_{1,\varkappa\_2}^2, \quad k \ge 0. \tag{44}$$

**Proof.** Starting from the l.h.s of (40),

$$
\left\langle \delta\_1 \left( \delta\_\chi \left( \varkappa\_2 \delta\_\chi p \right)\_i^{1/2} \right), -2\sigma p\_i^1 \right\rangle = \frac{-2\sigma}{\pi} \left\langle \delta\_\chi \left( \left( \varkappa\_2 \delta\_\chi p \right)\_i^1 \right) - \delta\_\chi \left( \left( \varkappa\_2 \delta\_\chi p \right)\_i^0 \right), p\_i^1 \right\rangle \tag{45}
$$

$$=\frac{2\sigma}{\pi}h\sum\_{l=0}^{M-1}\left[ (\varkappa\_2)\_{l+1/2}\delta\_x p\_{l+1/2}^1 \delta\_x p\_{l+1/2}^1 - (\varkappa\_2)\_{l+1/2}\delta\_x p\_{l+1/2}^0 \delta\_x p\_{l+1/2}^1 \right],\tag{46}$$

then the r.h.s of (40) is achieved directly. Additionally, starting from the l.h.s of (41),

$$
\left\langle \delta\_{\mathbf{x}} \left( \varkappa \mathfrak{L} \delta\_{\mathbf{x}} q \right)\_{i}^{1/2}, -2\pi p\_{i}^{1} \right\rangle = \frac{-2\sigma}{2} \left\langle \delta\_{\mathbf{x}} \left( \left( \varkappa \mathfrak{L} \delta\_{\mathbf{x}} q \right)\_{i}^{1} \right) + \delta\_{\mathbf{x}} \left( \left( \varkappa \mathfrak{L} \delta\_{\mathbf{x}} q \right)\_{i}^{0} \right), p\_{i}^{1} \right\rangle \tag{47}
$$

$$=\sigma h \sum\_{i=0}^{M-1} \left[ (\varkappa\_2)\_{i+1/2} \delta\_{\ge} q\_{i+1/2}^1 \delta\_{\ge} p\_{i+1/2}^1 + (\varkappa\_2)\_{i+1/2} \delta\_{\ge} q\_{i+1/2}^0 \delta\_{\ge} p\_{i+1/2}^1 \right], \tag{48}$$

then the r.h.s of (41) is immediately held. Invoking the previous estimates and Lemma 2.10, we deduce

$$\left\langle \partial\_{\boldsymbol{l}} \left( \delta\_{\boldsymbol{x}} \left( \varkappa\_{2} \delta\_{\boldsymbol{x}} p \right)\_{i}^{k+\sigma} \right), -p\_{i}^{k+\sigma} \right\rangle = \left\langle \partial\_{\boldsymbol{l}} p\_{i}^{k+\sigma}, p\_{i}^{k+\sigma} \right\rangle\_{\mathcal{H}} \tag{49}$$

$$\geq \frac{1}{4\pi} \left( \mathcal{E}^{k+1} - \mathcal{E}^k \right),\tag{50}$$

where <sup>E</sup>˜*k*+<sup>1</sup> is defined by (44) and so (42) is achieved. The estimate (43) is simply calculated.

*Mathematics* **2020**, *8*, 1696

Now, we are in a position to combine the Lemmas 3.1–3.3 to prove the convergence and stability of the proposed compact difference scheme. To that end, Let

$$
\rho\_i^k = \mathbf{V}\_i^k - \mathbf{V}\_{i\prime\prime}^k \quad \varepsilon\_i^k = \mathbf{V}\_i^k - \nu\_{i\prime}^k \quad 0 \le i \le M, \ 0 \le k \le N.
$$

Subtracting (31) from (26)–(28), (20c)–(20e), respectively, we obtain the error equations, as follows

$$\sum\_{n=0}^{k} \mathsf{C}\_{k-\iota}^{k+1} \left( \mathcal{A} \rho\_i^{n+1} - \mathcal{A} \rho\_i^n \right) = \delta\_x \left( \varkappa\_2 \delta\_x e \right)\_i^{k+\vartheta} + \mathcal{A} \left( \mathsf{F}\_i^{k+\vartheta} - \mathcal{F}\_i^{k+\vartheta} \right) + \mathcal{S}\_i^{k+\vartheta}, \quad 1 \le i \le M-1, \quad 0 \le k \le N-1,\tag{51a}$$

$$
\delta\_l \left( \delta\_x \left( \varkappa\_2 \delta\_t e \right)\_i^{1/2} \right) = \delta\_x \left( \varkappa\_2 \delta\_t \rho \right)\_i^{1/2} + s\_i^{1/2}, \quad 1 \le i \le M - 1,\tag{51b}
$$

$$\partial\_{\overline{i}}\left(\delta\_{\overline{x}}\left(\varkappa\mathfrak{z}\delta\_{\overline{x}}e\right)\_{i}^{k+\sigma}\right) = \delta\_{\overline{x}}\left(\varkappa\mathfrak{z}\delta\_{\overline{x}}\rho\right)\_{i}^{k+\sigma} + s\_{i}^{k+\sigma}, \quad 1 \le i \le M-1, \quad 1 \le k \le N-1,\tag{51c}$$

$$
\epsilon\_l^k = 0, \quad 1 \le l \le M - 1, \quad -n\_0 \le k \le 0, \quad \rho\_l^0 = 0, \quad 1 \le l \le M - 1,\tag{51d}
$$

$$\epsilon\_0^k = \epsilon\_M^k = 0, \quad 0 \le k \le N,\tag{51e}$$

$$
\rho\_0^k = \rho\_M^k = 0, \quad 0 \le k \le N. \tag{51}
$$

**Theorem 1.** *Assume that <sup>u</sup>*(*x*, *<sup>t</sup>*) <sup>∈</sup> *<sup>C</sup>*6,4 *<sup>x</sup>*,*t*([0, *<sup>L</sup>*] <sup>×</sup> [−*τ*, *<sup>T</sup>*]) *is the smooth solution of* (5) *and* {V*<sup>k</sup> <sup>i</sup>* , *<sup>ν</sup><sup>k</sup> <sup>i</sup>* |0 ≤ *i* ≤ *M*, −*n*<sup>0</sup> ≤ *k* ≤ *N*} *the numerical solution of the scheme* (31)*. Subsequently, there exist positive constants h*<sup>0</sup> *and τ*0*, independent of h and τ, such that, when h* ≤ *h*<sup>0</sup> *and τ* ≤ *τ*0*, we have the error estimate*

$$\left\| \epsilon^k \right\|\_{\infty} \le \mathbb{C}\_1 \left( \tau^2 + h^4 \right), \quad \tau \sum\_{n=1}^k \left\| \rho^n \right\| \le \mathbb{C}\_1 \left( \tau^2 + h^4 \right), \quad -n\_0 \le k \le N.$$

**Proof.** The proof will be preformed in two steps. Let us tackle the first one. **Step 1.** When *k* = 0, the system (51) is as follows

$$
\hat{\mathcal{L}}\_{0}^{1}\left(\mathcal{A}\rho\_{i}^{1} - \mathcal{A}\rho\_{i}^{0}\right) = \sigma\delta\_{x}\left(\varkappa\_{2}\delta\_{x}\varepsilon\right)\_{i}^{1} + (1-\sigma)\delta\_{x}\left(\varkappa\_{2}\delta\_{x}\varepsilon\right)\_{i}^{0} + \mathcal{S}\_{i}^{\sigma}, \quad 1 \le i \le M-1,\tag{52a}
$$

$$
\delta\_{\mathbb{I}}\left(\delta\_{\mathbb{X}}\left(\varkappa\_{2}\delta\_{\mathbb{X}}e\right)\_{i}^{1/2}\right) = \delta\_{\mathbb{X}}\left(\varkappa\_{2}\delta\_{\mathbb{X}}\rho\right)\_{i}^{1/2} + s\_{i}^{1/2}, \quad 1 \le i \le M-1,\tag{52b}
$$

$$
\sigma\_i^k = 0, \quad 1 \le i \le M - 1, \quad -n\_0 \le k \le 0, \quad \rho\_i^0 = 0, \quad 1 \le i \le M - 1,\tag{52c}
$$

$$
\epsilon\_0^0 = \epsilon\_M^0 = 0,\tag{52d}
$$

$$
\rho\_0^0 = \rho\_M^0 = 0.\tag{52e}
$$

Taking the inner product of (52a) with *ρ*1, we obtain

*e k*

$$\begin{split} \left. \hat{\mathsf{C}}\_{0}^{1} \right\| \rho^{1} \Big|\_{\mathcal{A}}^{1} &= \hat{\mathsf{C}}\_{0}^{1} \langle \mathcal{A}\rho^{0}, \rho^{1} \rangle + \langle \delta\_{\mathtt{x}} \left( \varkappa\_{2} \delta\_{\mathtt{x}} \varepsilon e^{\sigma} \right)^{\sigma}, \rho^{1} \rangle + \langle \mathcal{S}^{\sigma}, \rho^{1} \rangle \\ &= \hat{\mathsf{C}}\_{0}^{1} \langle \mathcal{A}\rho^{0}, \rho^{1} \rangle + \langle \sigma \delta\_{\mathtt{x}} \left( \varkappa\_{2} \delta\_{\mathtt{x}} e^{\parallel} \right)^{1} + (1 - \sigma) \delta\_{\mathtt{x}} \left( \varkappa\_{2} \delta\_{\mathtt{x}} e^{\parallel} \right)^{0}, \rho^{1} \rangle + \langle \mathcal{S}^{\sigma}, \rho^{1} \rangle \\ &= \hat{\mathsf{C}}\_{0}^{1} \langle \mathcal{A}\rho^{0}, \rho^{1} \rangle - \sigma \langle \varepsilon^{1}, \rho^{1} \rangle\_{\mathbb{A}\_{2}^{0}} + (1 - \sigma) \langle \delta\_{\mathtt{x}} \left( \varkappa\_{2} \delta\_{\mathtt{x}} e^{\parallel} \right)^{0}, \rho^{1} \rangle + \langle \mathcal{S}^{\sigma}, \rho^{1} \rangle, \end{split} \tag{53}$$

Taking the inner product of (52b) with <sup>−</sup>2*σe*<sup>1</sup> and invoking Lemma 3.3, we arrive at

$$\frac{2\sigma}{\pi} \left| \varepsilon^{1} \right|\_{1,\varkappa\_{2}}^{2} = \frac{2\sigma}{\pi} \left\langle \varepsilon^{0}, \varepsilon^{1} \right\rangle\_{\varkappa\_{2}} + \sigma \left( \left\langle \varepsilon^{1}, \rho^{1} \right\rangle\_{\varkappa\_{2}} + \left\langle \rho^{0}, \varepsilon^{1} \right\rangle\_{\varkappa\_{2}} \right) - 2\sigma \langle \mathbf{s}^{1/2}, \varepsilon^{1} \rangle. \tag{54}$$

Adding (53) with (54) and using Young inequality and Lemma 2.6, we obtain

2 3 *C*ˆ1 0 ) ) )*ρ*1 ) ) ) 2 + 2*σ τ e* 1 2 1,κ<sup>2</sup> <sup>≤</sup> *<sup>C</sup>*ˆ1 <sup>0</sup> A*ρ*0, *<sup>ρ</sup>*1 + (<sup>1</sup> <sup>−</sup> *<sup>σ</sup>*)*δ<sup>x</sup>* (κ2*δxe*) <sup>0</sup> , *<sup>ρ</sup>*1 <sup>+</sup> S*σ*, *<sup>ρ</sup>*1 (55) + 2*σ τ* . *e* 0,*e* 1 / κ2 + *σ* . *ρ*0,*e* 1 / κ2 − 2*σs* 1/2,*e* 1 ≤ 2 9 *C*ˆ1 0 ) ) )*ρ*1 ) ) ) 2 + 9 8 *C*ˆ1 0 ) ) )*ρ*0 ) ) ) 2 + 2*C*ˆ1 0 9 ) ) )*ρ*1 ) ) ) 2 + 9(1 − *σ*) 8*C*ˆ1 0 ) ) ) *δ<sup>x</sup>* (κ2*δxe*) 0 ) ) ) 2 + 2*C*ˆ1 0 9 ) ) )*ρ*1 ) ) ) 2 + 9 8*C*ˆ1 0 S*σ*<sup>2</sup> + 2*σ* 9*τ e* 1 2 1,κ<sup>2</sup> + 9*σ* 2*τ e* 0 2 1,κ<sup>2</sup> + 2*σ* 9*τ e* 1 2 1,κ<sup>2</sup> + 9*τσ* 8 *e* 0 2 1,κ<sup>2</sup> + 2*σ* 9*τ e* 1 2 1,κ<sup>2</sup> + 3*στ* <sup>2</sup> *<sup>L</sup>*<sup>2</sup> ) ) )*s* 1/2) ) ) 2 , (56)

after a simplification, we get

$$\begin{split} \left\| \left. \epsilon^{1} \right|\_{1,\varkappa\_{2}}^{2} \leq & \frac{2\mathsf{T}\tau}{32\sigma} \hat{\mathsf{C}}\_{0}^{1} \left\| \left. \rho^{0} \right| \right\|^{2} + \frac{2\mathsf{T}\tau(1-\sigma)}{32\sigma \hat{\mathsf{C}}\_{0}^{1}} \left\| \left. \delta\_{\mathsf{x}} \left( \varkappa\_{2} \delta\_{\mathsf{x}} \epsilon \right) \right\| \right\|^{2} + \frac{2\mathsf{T}\tau}{32\sigma \hat{\mathsf{C}}\_{0}^{1}} \left\| \left. \mathcal{S}^{\mathsf{r}} \right\| ^{2} + \frac{2\mathsf{T}}{8} \left| \epsilon^{0} \right|\_{1,\varkappa\_{2}}^{2} \\ & + \frac{2\mathsf{T}\tau^{2}}{32} \left| \epsilon^{0} \right|\_{1,\varkappa\_{2}}^{2} + \frac{9\tau^{2}}{8} L^{2} \left\| \left. \mathsf{s}^{1/2} \right\| \right\|^{2}, \end{split} \tag{57}$$

it follows from (52b) that

$$\left(\boldsymbol{\delta}\_{\boldsymbol{x}}\left(\varkappa\boldsymbol{\varrho}\boldsymbol{\delta}\_{\boldsymbol{x}}\boldsymbol{\epsilon}\right)\right)\_{\mathrm{i}}^{1} = \boldsymbol{\delta}\_{\boldsymbol{x}}\left(\varkappa\boldsymbol{\varrho}\boldsymbol{\delta}\_{\boldsymbol{x}}\boldsymbol{\epsilon}\right)\_{\mathrm{i}}^{0} + \tau\boldsymbol{\delta}\_{\boldsymbol{x}}\left(\varkappa\boldsymbol{\varrho}\boldsymbol{\delta}\_{\boldsymbol{x}}\boldsymbol{\rho}\right)\_{\mathrm{i}}^{1/2} + \tau\boldsymbol{\varepsilon}\_{\mathrm{i}}^{1/2}, \quad 1 \le \mathrm{i} \le M - 1. \tag{58}$$

Substituting (58) into (52a), we have

$$\hat{\mathcal{L}}\_{0}^{1}\left(\mathcal{A}\rho\_{i}^{1} - \mathcal{A}\rho\_{i}^{0}\right) = \left[\tau\sigma\delta\_{x}\left(\varkappa\_{2}\delta\_{x}\rho\right)\_{i}^{1/2} + \tau s\_{i}^{1/2}\right] + \delta\_{x}\left(\varkappa\_{2}\delta\_{x}e\right)\_{i}^{0} + \mathcal{S}\_{i}^{\sigma}, \quad 1 \le i \le M - 1. \tag{59}$$

Taking the inner product of (59) with *ρ*1/2, we obtain

$$
\langle \mathcal{L}\_0^1 \langle \left( \mathcal{A} \rho^1 - \mathcal{A} \rho^0 \right), \rho^{1/2} \rangle = -\pi \sigma \left| \rho^{1/2} \right|\_{1, \varkappa\_2}^2 + \langle \pi \sigma \mathbf{s}^{1/2}, \rho^{1/2} \rangle + \langle \delta\_x \left( \varkappa\_2 \delta\_x \mathbf{e} \right)^0, \rho^{1/2} \rangle \tag{60}
$$

$$+\langle \mathcal{S}^{\sigma}, \rho^{1/2} \rangle,\tag{61}$$

By summation by parts and the Young inequality *ab* <sup>≤</sup> <sup>1</sup> <sup>2</sup>*<sup>θ</sup> <sup>a</sup>*<sup>2</sup> <sup>+</sup> *<sup>θ</sup>* <sup>2</sup> *<sup>b</sup>*2, with *<sup>θ</sup>* <sup>=</sup> <sup>3</sup> *C*ˆ1 0 , this yields

*C*ˆ1 0 2 ) ) )*ρ*1 ) ) ) 2 <sup>A</sup> <sup>−</sup> ) ) )*ρ*0 ) ) ) 2 A <sup>=</sup> *δ<sup>x</sup>* (κ2*δxe*) <sup>0</sup> , *<sup>ρ</sup>*1/2 − *τσ ρ*1/2 2 1,κ<sup>2</sup> <sup>+</sup> S*σ*, *<sup>ρ</sup>*1/2 <sup>+</sup> *τσ<sup>s</sup>* 1/2, *<sup>ρ</sup>*1/2 ≤ 3 2*C*ˆ1 0 ) ) ) *δ<sup>x</sup>* (κ2*δxe*) 0 ) ) ) 2 + *C*ˆ1 0 6 ) ) )*ρ*1/2) ) ) 2 + 3 2*C*ˆ1 0 S*σ*<sup>2</sup> <sup>+</sup> *C*ˆ1 0 6 ) ) )*ρ*1/2) ) ) 2 + 3*τ*2*σ*<sup>2</sup> 2*C*ˆ1 0 ) ) )*s* 1/2) ) ) 2 + *C*ˆ1 0 6 ) ) )*ρ*1/2) ) ) 2 ≤ *C*ˆ1 0 4 ) ) )*ρ*0 ) ) ) 2 + ) ) )*ρ*1 ) ) ) 2 + 3 2*C*ˆ1 0 ) ) ) *δ<sup>x</sup>* (κ2*δxe*) 0 ) ) ) 2 + 3 2*C*ˆ1 0 S*σ*<sup>2</sup> <sup>+</sup> 3*τ*2*σ*<sup>2</sup> 2*C*ˆ1 0 ) ) )*s* 1/2) ) ) 2 . (62)

Subsequently, by invoking Lemma 2.6 and some simple manipulations, we get

$$\left\|\left\|\rho^{1}\right\|\right\|^{2}\leq9\left\|\left\|\rho^{0}\right\|\right\|^{2}+\frac{18}{\left(\hat{\mathsf{C}}\_{0}^{1}\right)^{2}}\left\|\delta\_{\mathbf{x}}\left(\varkappa\_{2}\delta\_{x}\mathbf{c}\right)^{0}\right\|^{2}+\frac{18}{\left(\hat{\mathsf{C}}\_{0}^{1}\right)^{2}}\left\|\mathcal{S}^{\mathsf{V}}\right\|^{2}+\frac{18\operatorname{r}^{2}\sigma^{2}}{\left(\hat{\mathsf{C}}\_{0}^{1}\right)^{2}}\left\|\mathbf{s}^{1/2}\right\|^{2}.\tag{63}$$

**Step 2**. When *<sup>k</sup>* <sup>≥</sup> 1, we take the inner product of (51a) with *<sup>ρ</sup>k*+*<sup>σ</sup>* and obtain

*Mathematics* **2020**, *8*, 1696

$$\begin{split} \left\langle \sum\_{n=0}^{k} \mathsf{C}\_{k-n}^{k+1} \left( \mathcal{A} \rho^{n+1} - \mathcal{A} \rho^{n} \right), \rho^{k+\sigma} \right\rangle &= \left\langle \delta\_{\mathbf{x}} \left( \varkappa \mathfrak{z} \delta\_{\mathbf{x}} \mathfrak{e} \right)^{k+\sigma}, \rho^{k+\sigma} \right\rangle \\ &+ \left\langle \mathcal{A} \left( \mathcal{F}^{k+\sigma} - \mathcal{F}^{k+\sigma} \right), \rho^{k+\sigma} \right\rangle + \left\langle \mathcal{S}^{k+\sigma}, \rho^{k+\sigma} \right\rangle, \quad 1 \le k \le N - 1. \end{split} \tag{64}$$

By Lemmas 2.4 and 2.10, we deduce for 1 ≤ *k* ≤ *N* − 1,

$$\begin{split} \left\langle \sum\_{n=0}^{k} \mathsf{C}\_{k-n}^{k+1} \left( \mathcal{A} \boldsymbol{\rho}^{n+1} - \mathcal{A} \boldsymbol{\rho}^{n} \right), \boldsymbol{\rho}^{k+\sigma} \right\rangle &\geq \frac{1}{2} \sum\_{n=0}^{k} \mathsf{C}\_{k-n}^{k+1} \left[ \left\| \boldsymbol{\rho}^{n+1} \right\|\_{\mathcal{A}}^{2} - \left\| \boldsymbol{\rho}^{n} \right\|\_{\mathcal{A}}^{2} \right] \\ &= \frac{1}{2} \left( \sum\_{n=1}^{k+1} \mathsf{C}\_{k-n+1}^{k+1} \left\| \boldsymbol{\rho}^{n} \right\|\_{\mathcal{A}}^{2} - \sum\_{n=1}^{k} \mathsf{C}\_{k-n}^{k} \left\| \boldsymbol{\rho}^{n} \right\|\_{\mathcal{A}}^{2} - \mathsf{E}\_{k} \left\| \boldsymbol{\rho}^{1} \right\|\_{\mathcal{A}}^{2} - \mathsf{C}\_{k}^{k+1} \left\| \boldsymbol{\rho}^{0} \right\|\_{\mathcal{A}}^{2} \right). \end{split} \tag{65}$$

Young inequality is used for any *θ* > 0 to yield

$$\left| \langle \mathcal{S}^{k+\sigma}, \rho^{k+\sigma} \rangle \right| \le \theta \left\| \rho^{k+\sigma} \right\|^2 + \frac{1}{4\theta} \left\| \mathcal{S}^{k+\sigma} \right\|^2. \tag{66}$$

Using Lemma 3.1 and also Young inequality, this gives

$$\begin{split} \mathcal{A}\left(\mathcal{A}\left(\mathcal{F}^{k+\sigma}-\mathcal{F}^{k+\sigma}\right),\rho^{k+\sigma}\right) &\leq \theta \left\|\rho^{k+\sigma}\right\|^{2} + \frac{1}{4\theta} \left\|\mathcal{F}^{k+\sigma}-\mathcal{F}^{k+\sigma}\right\|\_{\mathcal{A}}^{2} \\ &\leq \frac{\theta}{2} \left\|\rho^{k+\sigma}\right\|^{2} + \frac{1}{2\theta} \left[2c\_{1}^{2}\left(\left\|\rho^{k}\right\|^{2} + \left\|\rho^{k-1}\right\|^{2}\right) + c\_{2}^{2}\left(\left\|\rho^{k+1-\mu\_{0}}\right\|^{2} + \left\|\rho^{k-\mu\_{0}}\right\|^{2}\right)\right]. \end{split} \tag{67}$$

Substituting (65)–(67) in (64), we get, for all 1 ≤ *k* ≤ *N* − 1,

$$\begin{split} \frac{1}{2} \left( \sum\_{n=1}^{k+1} \hat{\mathcal{L}}\_{k-n+1}^{k+1} \left\| \boldsymbol{\rho}^{n} \right\|\_{\mathcal{A}}^{2} - \sum\_{n=1}^{k} \hat{\mathcal{L}}\_{k-n}^{k} \left\| \boldsymbol{\rho}^{n} \right\|\_{\mathcal{A}}^{2} - \hat{\mathcal{b}}\_{k} \left\| \boldsymbol{\rho}^{1} \right\|\_{\mathcal{A}}^{2} - \hat{\mathcal{C}}\_{k}^{k+1} \left\| \boldsymbol{\rho}^{0} \right\|\_{\mathcal{A}}^{2} \right) \\ \leq & \left( \delta\_{x} \left( \varkappa\_{2} \delta\_{x} \boldsymbol{\epsilon} \right)^{k+\sigma}, \boldsymbol{\rho}^{k+\sigma} \right) + \theta \left\| \boldsymbol{\rho}^{k+\sigma} \right\|^{2} + \frac{1}{4\theta} \left\| \boldsymbol{S}^{k+\sigma} \right\|^{2} \\ \quad + \frac{\theta}{2} \left\| \boldsymbol{\rho}^{k+\sigma} \right\|^{2} + \frac{1}{2\theta} \left[ 2c\_{1}^{2} \left( \left\| \boldsymbol{\rho}^{k} \right\|^{2} + \left\| \boldsymbol{\rho}^{k-1} \right\|^{2} \right) + c\_{2}^{2} \left( \left\| \boldsymbol{\rho}^{k+1-n\_{0}} \right\|^{2} + \left\| \boldsymbol{\rho}^{k-n\_{0}} \right\|^{2} \right) \right] . \end{split} \tag{68}$$

Taking the inner product of (51c) with <sup>−</sup>*ek*+*σ*, we obtain

$$- \left\langle \partial\_{\overline{t}} \left( \delta\_{\overline{x}} \left( \varkappa\_2 \delta\_{\overline{x}} \varrho \right)^{k+\sigma} \right), \varepsilon^{k+\sigma} \right\rangle = - \left\langle \delta\_{\overline{x}} \left( \varkappa\_2 \delta\_{\overline{x}} \varrho \right)^{k+\sigma}, \varepsilon^{k+\sigma} \right\rangle - \left\langle \mathbf{s}^{k+\sigma}, \varepsilon^{k+\sigma} \right\rangle, \quad 1 \le k \le N-1. \tag{69}$$

For the l.h.s of (69), after recalling Lemmas 2.10 and 3.3, we get

$$- \left\langle \partial\_{\mathbf{f}} \left( \delta\_{\mathbf{x}} \left( \varkappa\_{2} \delta\_{\mathbf{x}} e \right)^{k+\upsilon} \right), \varepsilon^{k+\upsilon} \right\rangle \geq \frac{1}{4\pi} \left( \mathcal{E}^{k+1} - \mathcal{E}^{k} \right), \quad 1 \leq k \leq N - 1,\tag{70}$$

such that

$$\mathcal{E}^{k+1} = (2\sigma + 1)|\varepsilon^{k+1}|^2\_{1,\varkappa\_2} - (2\sigma - 1)|\varepsilon^k|^2\_{1,\varkappa\_2} + (2\sigma^2 + \sigma - 1)|\varepsilon^{k+1} - \varepsilon^k|^2\_{1,\varkappa\_2\prime} \quad k \ge 0, \varkappa\_2$$

and additionally,

$$\tilde{\mathcal{E}}^{k+1} \ge \frac{1}{\sigma} |e^{k+1}|\_{1, \varkappa\_2 \prime}^2 \quad k \ge 0. \tag{71}$$

By the Cauchy–Schwarz inequality, we have

$$\left| - \left\langle s^{k+\sigma}, e^{k+\sigma} \right\rangle \right| \le \frac{1}{2} \left\| e^{k+\sigma} \right\|^2 + \frac{1}{2} \left\| s^{k+\sigma} \right\|^2, \quad 1 \le k \le N-1,\tag{72}$$

so plugging (70) and (72) into (69), this gives

$$\frac{1}{4\pi} \left( \vec{\mathcal{E}}^{k+1} - \vec{\mathcal{E}}^{k} \right) \le - \left\langle \delta\_{\mathcal{X}} \left( \varkappa\_{2} \delta\_{\mathcal{X}} \rho \right)^{k+\sigma}, e^{k+\sigma} \right\rangle + \frac{1}{2} \left\| \epsilon^{k+\sigma} \right\|^{2} + \frac{1}{2} \left\| s^{k+\sigma} \right\|^{2}, \quad 1 \le k \le N - 1. \tag{73}$$

Adding (68) and (73), we obtain

$$\begin{split} \frac{1}{2} \Big( \sum\_{n=1}^{k+1} \hat{\epsilon}\_{k-n+1}^{k+1} \left\lVert \boldsymbol{\rho}^{\boldsymbol{u}} \right\rVert\_{\mathcal{A}}^{2} - \sum\_{n=1}^{k} \hat{\epsilon}\_{k-n}^{k} \left\lVert \boldsymbol{\rho}^{\boldsymbol{u}} \right\rVert\_{\mathcal{A}}^{2} - \hat{b}\_{k} \left\lVert \boldsymbol{\rho}^{\boldsymbol{1}} \right\rVert\_{\mathcal{A}}^{2} - \hat{\epsilon}\_{k}^{k+1} \left\lVert \boldsymbol{\rho}^{\boldsymbol{0}} \right\rVert\_{\mathcal{A}}^{2} \Big) + \frac{1}{4\pi} \left( \bar{\mathcal{E}}^{k+1} - \bar{\mathcal{E}}^{k} \right) \\ \leq \frac{1}{2} \left\lVert \boldsymbol{\epsilon}^{k+\sigma} \right\rVert^{2} + \frac{1}{2} \left\lVert \boldsymbol{\rho}^{k+\sigma} \right\rVert^{2} + \theta \left\lVert \boldsymbol{\rho}^{k+\sigma} \right\rVert^{2} + \frac{1}{4\theta} \left\lVert \boldsymbol{\rho}^{k+\sigma} \right\rVert^{2} \\ \quad + \frac{\theta}{2} \left\lVert \boldsymbol{\rho}^{k+\sigma} \right\rVert^{2} + \frac{1}{2\theta} \left[ 2\boldsymbol{\epsilon}\_{1}^{2} \left( \left\lVert \boldsymbol{\rho}^{k} \right\rVert^{2} + \left\lVert \boldsymbol{\rho}^{k-1} \right\rVert^{2} \right) + \boldsymbol{\epsilon}\_{2}^{2} \left( \left\lVert \boldsymbol{\rho}^{k+1-n} \right\rVert \right)^{2} + \left\lVert \boldsymbol{\rho}^{k-n} \right\rVert \right]^{2} \right]. \end{split}$$

Now, multiply both sides of (74) by 4*τ*, use Lemma 2.6 and do some arrangements, we then have

$$\begin{split} 2\tau \sum\_{n=1}^{k+1} \xi\_{k-n+1}^{k+1} \left\lVert \boldsymbol{\rho}^{n} \right\rVert\_{\mathcal{A}}^{2} + \mathcal{E}^{k+1} \leq 2\tau \left( \sum\_{n=1}^{k} \xi\_{k-n}^{k} \left\lVert \boldsymbol{\rho}^{n} \right\rVert\_{\mathcal{A}}^{2} \right) + \mathcal{E}^{k} + 2\tau \left( \delta\_{\mathtt{k}} \left\lVert \boldsymbol{\rho}^{1} \right\rVert^{2} + \mathcal{E}\_{\mathtt{k}}^{k+1} \left\lVert \boldsymbol{\rho}^{0} \right\rVert^{2} \right) \\ &+ \frac{4\tau}{2} \left\lVert \boldsymbol{\rho}^{k+\sigma} \right\rVert^{2} + \frac{4\tau}{2} \left\lVert \boldsymbol{\rho}^{k+\sigma} \right\rVert^{2} + \frac{4\tau}{4\theta} \left\lVert \boldsymbol{\mathcal{S}}^{k+\sigma} \right\rVert^{2} + 2\tau (1+2\theta) \left\lVert \boldsymbol{\rho}^{k+\sigma} \right\rVert^{2} \\ &+ \frac{4\tau}{2} \left[ 2c\_{1}^{2} \left( \left\lVert \boldsymbol{\rho}^{k} \right\rVert^{2} + \left\lVert \boldsymbol{\rho}^{k-1} \right\rVert^{2} \right) + c\_{2}^{2} \left( \left\lVert \boldsymbol{\rho}^{k+1-n\_{0}} \right\rVert^{2} + \left\lVert \boldsymbol{\rho}^{k-n\_{0}} \right\rVert^{2} \right) \right]. \end{split} \tag{74}$$

Denote <sup>H</sup>*k*+<sup>1</sup> <sup>=</sup> <sup>2</sup>*<sup>τ</sup>* <sup>∑</sup>*k*+<sup>1</sup> *<sup>n</sup>*=<sup>1</sup> *<sup>C</sup>*<sup>ˆ</sup> *<sup>k</sup>*+<sup>1</sup> *<sup>k</sup>*−*n*+<sup>1</sup> *ρn*<sup>2</sup> <sup>A</sup> <sup>+</sup> <sup>E</sup>˜*k*<sup>+</sup>1,, we can write

<sup>H</sup>*k*+<sup>1</sup> ≤ H*<sup>k</sup>* <sup>+</sup> <sup>2</sup>*<sup>τ</sup>* ˆ *bk* ) ) )*ρ*1 ) ) ) 2 + *C*ˆ *<sup>k</sup>*+<sup>1</sup> *k* ) ) )*ρ*0 ) ) ) 2 + 4*τ* 2 ) ) )*e k*+*σ* ) ) ) 2 + 4*τ* 2 ) ) )*s k*+*σ* ) ) ) 2 + 4*τ* 4*θ* ) ) ) <sup>S</sup>*k*+*<sup>σ</sup>* ) ) ) 2 + 8*τθ* ) ) )*ρk*+*<sup>σ</sup>* ) ) ) 2 + 4*τ θ* 2*c*<sup>2</sup> 1 ) ) )*ρk* ) ) ) 2 + ) ) )*ρk*−<sup>1</sup> ) ) ) 2 + *c*<sup>2</sup> 2 ) ) )*ρk*+1−*n*<sup>0</sup> ) ) ) 2 + ) ) )*ρk*−*n*<sup>0</sup> ) ) ) 2 ≤ H<sup>1</sup> <sup>+</sup> <sup>2</sup>*<sup>τ</sup> k* ∑ *n*=1 ˆ *bn* ) ) )*ρ*1 ) ) ) 2 + *k* ∑ *n*=1 *<sup>C</sup>*<sup>ˆ</sup> *<sup>k</sup>*+<sup>1</sup> *<sup>n</sup>* ) ) )*ρ*0 ) ) ) 2 + 8*τ* 2 *k*+1 ∑ *n*=1 *e <sup>n</sup>*<sup>2</sup> <sup>+</sup> 4*τ* 2 *k* ∑ *n*=1 ) )*s n*+*σ*) )<sup>2</sup> (75) + 4*τ* 4*θ k* ∑ *n*=1 ) )S*n*+*σ*) )<sup>2</sup> <sup>+</sup> <sup>4</sup>*τ*(<sup>1</sup> <sup>+</sup> <sup>2</sup>*θ*) *k*+1 ∑ *n*=1 *ρn*<sup>2</sup> <sup>+</sup> 4*τ* 2 *k* ∑ *n*=1 2*c*<sup>2</sup> 1 *ρn*<sup>2</sup> <sup>+</sup> ) ) )*ρn*−<sup>1</sup> ) ) ) 2 + *c*<sup>2</sup> 2 ) ) )*ρn*+1−*n*<sup>0</sup> ) ) ) 2 + ) )*ρn*−*n*<sup>0</sup> ) )2 , 1 ≤ *k* ≤ *N* − 1.

From (71) and Lemma 2.4,

$$\mathcal{H}^{k+1} \ge \tau \sum\_{r=0}^{m} p\_r \frac{(1 - \gamma\_r)^{T-\gamma\_r}}{\Gamma(2 - \gamma\_r)} (k + \sigma)^{1 - \gamma\_r} \sum\_{n=1}^{k+1} \|\rho^n\|\_{\mathcal{A}}^2 + \frac{1}{\sigma} |\epsilon^{k+1}|\_{1, \varkappa\_2 \prime}^2 \quad 1 \le k \le N - 1,\tag{76}$$

$$\mathcal{H}^1 = 2\tau \hat{\mathbb{C}}\_0^1 \left\| \rho^1 \right\|\_{\mathcal{A}}^2 + \mathcal{E}^1 \le 2\tau \hat{\mathbb{C}}\_0^1 \left\| \rho^1 \right\|^2 + (2\sigma + 1)|\varepsilon^1|\_{1, \varkappa\_2}^2 + (2\sigma^2 + \sigma - 1)|\varepsilon^1|\_{1, \varkappa\_2}^2. \tag{77}$$

Inserting the above two inequalities into (76) and considering

$$\mathcal{L} := \max \{ c\_1, c\_2 \}, \quad \mathcal{A} := \frac{2}{3} \sum\_{r=0}^{m} p\_r \frac{(1 - \gamma\_r) \, T^{-\gamma\_r}}{\Gamma(2 - \gamma\_r)} (k + \sigma)^{1 - \gamma\_r},$$

yields

$$\begin{split} \frac{1}{\tau} \|\boldsymbol{\varepsilon}^{k+1}\|\_{1,\varkappa\_{2}}^{2} + \tau \boldsymbol{A} \sum\_{n=1}^{k+1} \left\|\boldsymbol{\rho}^{n}\right\|^{2} &\leq 2\tau \mathsf{C}\_{0}^{1} \left\|\boldsymbol{\rho}^{1}\right\|^{2} + (2\boldsymbol{\sigma}^{2} + 3\boldsymbol{\tau}) \|\boldsymbol{\varepsilon}^{1}\|\_{1,\varkappa\_{2}}^{2} \\ &+ 2\tau \left(\sum\_{n=1}^{k} \mathsf{b}\_{n} \left\|\boldsymbol{\rho}^{1}\right\|^{2} + \sum\_{n=1}^{k} \mathsf{C}\_{n}^{k+1} \left\|\boldsymbol{\rho}^{0}\right\|^{2}\right) + 4\tau \sum\_{n=1}^{k+1} \left\|\boldsymbol{\varepsilon}^{n}\right\|^{2} + 2\tau \sum\_{n=1}^{k} \left\|\boldsymbol{s}^{n+\sigma}\right\|^{2} \\ &+ \frac{\tau}{\theta} \sum\_{n=1}^{k} \left\|\boldsymbol{\mathcal{S}}^{n+\sigma}\right\|^{2} + 4\tau (1 + 2\theta + 3\boldsymbol{\varepsilon}^{2}) \sum\_{n=1}^{k+1} \left\|\boldsymbol{\rho}^{n}\right\|^{2}. \end{split}$$

By choosing *<sup>θ</sup>* to achieve A <sup>≥</sup> <sup>4</sup>(<sup>1</sup> <sup>+</sup> <sup>2</sup>*<sup>θ</sup>* <sup>+</sup> <sup>3</sup>*c*2), invoking (57), (63), and denoting

G*k*+<sup>1</sup> := 2*τ C*ˆ1 <sup>0</sup> + *k* ∑ *n*=1 ˆ *bn* 9 ) ) )*ρ*0 ) ) ) 2 + 18 *C*ˆ1 0 2 ) ) ) *δ<sup>x</sup>* (κ2*δxe*) 0 ) ) ) 2 + 18 *C*ˆ1 0 <sup>2</sup> S*σ*<sup>2</sup> <sup>+</sup> 18*τ*2*σ*<sup>2</sup> *C*ˆ1 0 2 ) ) )*s* 1/2) ) ) 2 + (2*σ*<sup>2</sup> + 3*σ*) 27*τ* 32*σC*ˆ1 0 ) ) )*ρ*0 ) ) ) 2 + 27*τ*(1 − *σ*) 32*σC*ˆ1 0 ) ) ) *δ<sup>x</sup>* (κ2*δxe*) 0 ) ) ) 2 + 27*τ* 32*σC*ˆ1 0 S*σ*<sup>2</sup> <sup>+</sup> 27 8 *e* 0 2 1,κ<sup>2</sup> (78) + 27*τ*<sup>2</sup> 32 *e* 0 2 1,κ<sup>2</sup> + 9*τ*<sup>2</sup> <sup>8</sup> *<sup>L</sup>*<sup>2</sup> ) ) )*s* 1/2) ) ) 2 + 2*τ k* ∑ *n*=1 *<sup>C</sup>*<sup>ˆ</sup> *<sup>k</sup>*+<sup>1</sup> *<sup>n</sup>* ) ) )*ρ*0 ) ) ) 2 + 2*τ k* ∑ *n*=1 ) )*s n*+*σ*) )<sup>2</sup> <sup>+</sup> *<sup>τ</sup> θ k* ∑ *n*=1 ) )S*n*+*σ*) )2 ,

we obtain directly after following the assumptions (13),

$$|e^{k+1}|\_1^2 \le \frac{1}{b\_2} |e^{k+1}|\_{1, \varkappa\_2}^2 \le \frac{4\pi\sigma}{b\_2} \sum\_{n=1}^{k+1} ||e^n||^2 + \frac{\sigma}{b\_2} \mathcal{G}\_{k+1}, \quad 1 \le k \le N-1,\tag{79}$$

invoking Lemma 2.6 gives

$$\|e^{k+1}\|\_{1}^{2} \le \frac{2\tau L^{2}\sigma}{3b\_{2}}\sum\_{n=1}^{k+1} |e^{n}|\_{1}^{2} + \frac{\sigma}{b\_{2}}\mathcal{G}\_{k+1\prime} \quad 1 \le k \le N-1,\tag{80}$$

applying Grönwall Lemma 3.2 yields

$$|e^{k+1}|\_1^2 \le \frac{\sigma}{b\_2} \exp\left(\frac{4L^2\sigma}{3b\_2}\right) \mathcal{G}\_{k+1\prime} \quad 1 \le k \le N-1. \tag{81}$$

Accordingly, the proof is completed.

## **4. Almost Unconditional Stability**

To discuss the stability of the compact difference scheme (31), we also use the discrete energy method. Let {V¯ *<sup>k</sup> <sup>i</sup>* , *<sup>ν</sup>*¯*<sup>k</sup> <sup>i</sup>* |0 ≤ *i* ≤ *M*, −*n*<sup>0</sup> ≤ *k* ≤ *N*} be the solution of

*Mathematics* **2020**, *8*, 1696

$$\sum\_{n=0}^{k} \mathcal{L}\_{k-n}^{k+1} \left( \mathcal{A} \mathcal{V}\_{l}^{n+1} - \mathcal{A} \mathcal{V}\_{l}^{n} \right) = \delta\_{\mathbf{x}} \left( \varkappa\_{2} \delta\_{\mathbf{x}} \mathcal{O} \right)\_{l}^{k+\sigma} + \mathcal{A} \mathbb{P}\_{l}^{k+\sigma}, \quad 1 \le i \le M-1, \quad 0 \le k \le N-1,\tag{82a}$$

$$
\delta\_l \left( \delta\_\chi \left( \varkappa\_2 \delta\_\mathfrak{L} \mathfrak{I} \right)\_l^{1/2} \right) = \delta\_\chi \left( \varkappa\_2 \delta\_\mathfrak{L} \mathfrak{I} \right)\_l^{1/2}, \quad 1 \le i \le M - 1,\tag{82b}
$$

$$\left(\delta\_{\mathbb{X}}\left(\boldsymbol{\delta}\_{\mathbb{X}}\left(\varkappa\_{2}\delta\_{\mathbb{X}}\boldsymbol{\vartheta}\right)\right)^{\underline{k}+\sigma}\right) = \delta\_{\mathbb{X}}\left(\varkappa\_{2}\delta\_{\mathbb{X}}\boldsymbol{\vartheta}\right)^{\underline{k}+\sigma}, \quad 1 \le i \le M-1, \quad 1 \le k \le N-1,\tag{82c}$$

$$\begin{aligned} \mathfrak{d}\_{l}^{k} &= \mathfrak{k}\_{l}^{k} + \varrho\_{l}^{k}, \quad 1 \le i \le M - 1, \quad -n\_{0} \le k \le 0, \quad \mathfrak{d}\_{l}^{0} = \hat{\mathfrak{y}}(\mathbf{x}\_{l}), \quad 1 \le i \le M - 1,\tag{82d} \\\mathfrak{d}\_{0}^{k} &= \mathfrak{d}\_{M}^{k} = 0, \quad 0 \le k \le N. \end{aligned} \tag{82e}$$

$$
\vartheta\_0^\* = \vartheta\_M^\* = 0, \quad 0 \le k \le N,\tag{82e}
$$

$$
\mathcal{V}\_0^k = \mathcal{V}\_M^k = 0, \quad 0 \le k \le N,\tag{82f}
$$

where

$$\mathbb{P}\_{i}^{k+\sigma} = \tilde{f}\left( (\sigma+1)\,\boldsymbol{\upsilon}\_{i}^{k} - \sigma \boldsymbol{\upsilon}\_{i}^{k-1}, \sigma \boldsymbol{\upsilon}\_{i}^{k+1-n\_{0}} + (1-\sigma)\boldsymbol{\upsilon}\_{i}^{k-n\_{0}}, \boldsymbol{x}\_{i}, t\_{k+\sigma} \right). \tag{83}$$

where *<sup>k</sup> <sup>i</sup>* denotes an initial perturbation term that is very small.

**Theorem 2.** *Let ρ*¯*<sup>k</sup> <sup>i</sup>* <sup>=</sup> <sup>V</sup>¯ *<sup>k</sup> <sup>i</sup>* − V*<sup>k</sup> <sup>i</sup>* , *e*¯ *k <sup>i</sup>* = *<sup>ν</sup>*¯*<sup>k</sup> <sup>i</sup>* <sup>−</sup> *<sup>ν</sup><sup>k</sup> <sup>i</sup> , for* 0 ≤ *i* ≤ *M*, −*n*<sup>0</sup> ≤ *k* ≤ *N*. *Subsequently, there exist constants c*4, *c*5, *h*0, *τ*<sup>0</sup> *that fulfill*

$$\left\| \left. \ell^k \right\| \_{\infty} \le c\_4 \sum\_{k=-n}^0 \left\| \left. \ell^k \right\| \_{\ell'} \right\| \_{\ell'} \quad 0 \le k \le N\_{\prime}$$

*conditioned by*

$$h \le h\_{0\prime} \quad \tau \le \tau\_{0\prime} \quad \max\_{\substack{-n \le k \le 0 \\ 0 \le i \le M}} \left| \vartheta\_i^k \right| \le c\_5.$$

**Proof.** The perturbation equations in terms of *ρ*¯*<sup>k</sup> <sup>i</sup>* and *e*¯ *k <sup>i</sup>* come by subtracting (82) from (31) and similar to the proof of Theorem 1, the conclusion of stability holds immediately.

#### **5. Generalized Scheme for the Distributed Order Case**

We are now in a position to consider the distributed order form of dmfCDWEs, which means that

$$\int\_{1}^{2} \omega(u) \frac{\partial^{\mu} u(\mathbf{x}, t)}{\partial t^{\alpha}} du = \frac{\partial}{\partial \mathbf{x}} \left( q\_{1}(\mathbf{x}) \frac{\partial u}{\partial \mathbf{x}} \right) + q\_{2}(\mathbf{x}) \frac{\partial u}{\partial \mathbf{x}} + f(u(\mathbf{x}, t), u(\mathbf{x}, t-s), \mathbf{x}, t), \ 0 < t \le T, \ 0 \le \mathbf{x} \le L,\tag{84a}$$

with the following initial and boundary conditions

$$u(\mathbf{x},t) = d(\mathbf{x},t), \quad 0 \le \mathbf{x} \le L, \quad t \in [-s,0), \quad \frac{\partial u(\mathbf{x},0)}{\partial t} = \psi(\mathbf{x}) = \lim\_{t \to -0} \frac{\partial d(\mathbf{x},t)}{\partial t}, \tag{84b}$$

$$u(0,t) = \phi \circ (t), \quad u(L,t) = \phi\_L(t), \quad 0 < t \le T,\tag{84c}$$

Following the same manipulations illustrated before; starting from an exponential transformation technique, then a transformation to zero Dirichlet boundary conditions, we obtain the following system

$$\int\_{1}^{2} \omega(\mathbf{a}) \frac{\partial^{\mathbf{a}} \nu(\mathbf{x}, t)}{\partial t^{\mathbf{a}}} d\mathbf{a} = \frac{\partial}{\partial \mathbf{x}} \left( q\_{1}(\mathbf{x}) \frac{\partial \nu}{\partial \mathbf{x}} \right) + \tilde{f}(\nu(\mathbf{x}, t), \nu(\mathbf{x}, t - \mathbf{s}), \mathbf{x}, t), \ 0 < t \le T, \ 0 \le \mathbf{x} \le L,\tag{85a}$$

with the initial and boundary conditions as

$$v(\mathbf{x},t) = \mathbf{f}(\mathbf{x},t), \quad 0 \le \mathbf{x} \le L, \quad t \in [-s,0), \quad \frac{\partial v(\mathbf{x},0)}{\partial t} = \hat{\psi}(\mathbf{x}) = \lim\_{t \to -0} \frac{\partial r(\mathbf{x},t)}{\partial t}, \tag{85b}$$

$$
\upsilon(0,t) = \upsilon(L,t) = 0, \quad t > 0. \tag{85c}
$$

*Mathematics* **2020**, *8*, 1696

A numerical quadrature rule can be adapted to transform the system (85) to dmfCDWEs. We recall Simpson's rule (also known as the three-point Newton–Cotes quadrature rule), a proof of which can be found in any descent textbook.

**Lemma 5.1.** *Consider an equidistant partition of the interval* [1, 2] *into* 2*J subintervals, let* Δ*α* = <sup>1</sup> <sup>2</sup>*<sup>J</sup> and denote α<sup>l</sup>* = 1 + *l* Δ*α,* 0 ≤ *l* ≤ 2*J. Subsequently, the composite Simpson's rule reads*

$$\int\_{1}^{2} f(a) da = \Delta a \sum\_{l=0}^{2l} \gamma\_{l} f(a\_{l}) - \frac{(\Delta a)^{4}}{180} f^{(4)}(\zeta), \quad \zeta \in [1, 2]. \tag{86}$$

*where*

$$\gamma\_l = \begin{cases} \frac{1}{3}, & l = 0, 2J, \\ \frac{2}{3}, & l = 2, 4, \dots, 2J - 4, 2J - 2J, \\ \frac{4}{3}, & l = 1, 3, \dots, 2J - 3, 2J - 1. \end{cases}$$

Define the function *G*(· ; *xi*, *tj*) : *α* → *ω*(*α*) *∂αv*(*xi*,*tj*) *<sup>∂</sup>t<sup>α</sup>* . Suppose that *<sup>G</sup>*(*α*) <sup>∈</sup> C4([1, 2]), then by using Lemma 5.1, we approximate the distributed derivative as

$$\begin{split} \int\_{1}^{2} \omega(a) \frac{\partial^{a} \boldsymbol{v}(\mathbf{x}\_{i}, t\_{k+\sigma})}{\partial t^{a}} da &= \Delta a \sum\_{l=0}^{2l} \gamma\_{l} \omega(a\_{l}) \left. \frac{c}{0} D\_{t}^{a\_{l}} \boldsymbol{v}(\mathbf{x}\_{i}, t\_{k-1/2}) - \frac{(\Delta a)^{4}}{180} G^{(4)}(a; \mathbf{x}\_{i}, t\_{k+\sigma}) \right|\_{a = \tilde{\mathbf{x}}\_{i}^{k}},\\ &= \Delta a \sum\_{l=0}^{2l} \gamma\_{l} \omega(a\_{l}) \left. \frac{c}{0} D\_{t}^{a\_{l}} \boldsymbol{v}(\mathbf{x}\_{i}, t\_{k+\sigma}) + O \left( \Delta a \right)^{4}, \end{split} \tag{87}$$

for a *ζ<sup>k</sup> <sup>i</sup>* ∈ [1, 2]. Define

$$\mathring{\mathsf{C}}\_{n}^{(k+1)} = \Delta a \sum\_{l=0}^{2J} \gamma\_{l} \omega(a\_{l}) \frac{\mathsf{r}^{-\gamma\_{l}}}{\Gamma(2-\gamma\_{l})} \mathsf{C}\_{n}^{(k+1,\gamma\_{l})}, \quad n = 0, 1, \ldots, k, 1$$

then the constructed compact difference scheme has the following form

$$\sum\_{n=0}^{k} \mathcal{L}\_{k-n}^{k+1} \left( \mathcal{A} \mathcal{V}\_{l}^{n+1} - \mathcal{A} \mathcal{V}\_{l}^{n} \right) = \delta\_{\mathbf{x}} \left( \varkappa\_{2} \delta\_{\mathbf{x}} v \right)\_{l}^{k+\sigma} + \mathcal{A} \mathcal{F}\_{l}^{k+\sigma}, \quad 1 \le i \le M - 1, \quad 0 \le k \le N - 1,\tag{88a}$$

$$
\delta\_l \left( \delta\_\mathbf{x} \left( \varkappa\_2 \delta\_\mathbf{x} \psi \right)^{1/2}\_l \right) = \delta\_\mathbf{x} \left( \varkappa\_2 \delta\_\mathbf{x} \mathcal{V} \right)^{1/2}\_l, \quad 1 \le l \le M - 1,\tag{88b}
$$

$$
\partial\_l \left( \delta\_{\mathbf{x}} \left( \varkappa\_2 \delta\_{\mathbf{x}} \nu \right)\_l^{k+\sigma} \right) = \delta\_{\mathbf{x}} \left( \varkappa\_2 \delta\_{\mathbf{x}} \mathcal{V} \right)\_l^{k+\sigma}, \quad 1 \le i \le M-1, \quad 1 \le k \le N-1,\tag{88c}
$$

$$\nu\_l^k = \mathbb{M}\_l^k, \quad 1 \le l \le M - 1, \quad -u\_0 \le k \le 0, \quad \mathcal{V}\_l^0 = \hat{\mathfrak{p}}(\mathbf{x}\_l), \quad 1 \le l \le M - 1,\tag{88d}$$

$$
\omega\_0^k = \nu\_M^k = 0, \quad 0 \le k \le N,\tag{88e}
$$

$$\mathcal{V}\_0^k = \mathcal{V}\_M^k = 0, \quad 0 \le k \le N. \tag{886}$$

**Remark 1.** The local truncation error of the compact difference scheme (88) for the distributed order system (85) is of order O *τ*<sup>2</sup> + *h*<sup>4</sup> + (Δ*α*) 4 . The convergence and stability estimates can be derived in the same manner as in Theorem 1 and Theorem 2.
