$$\int\_0^1 \lambda\_1 u\_0(t) \sin(\pi t) dt \ge (1+\varepsilon) \int\_0^1 \lambda\_1 u\_0(t) \sin(\pi t) dt - \frac{2\mathcal{C}}{\pi} \lambda\_1$$

which yields

$$\int\_0^1 \lambda\_1 u\_0(t) \sin(\pi t) dt \le \frac{2C}{\pi \varepsilon \lambda\_1}.$$

On the other hand, in view of the definition of cone *P*, one can easily obtain that

$$\|\|u\_0\|\int\_0^1 \sigma(t)\sin(\pi t)dt \le \int\_0^1 u\_0(t)\sin(\pi t)dt \le \frac{2\mathcal{C}}{\pi \varepsilon \lambda\_1}.$$

which means

4

4

$$\|\|\mu\_0\|\| \le \frac{2\mathcal{C}}{\pi \varepsilon \lambda\_1 \int\_0^1 \sigma(t) \sin(\pi t) dt} := R\_1^\*. \tag{10}$$

Therefore, if *R*1 > *R* ∗ 1 , immediately, one can get *µT*(*u*, *v*) 6= (*u*, *v*) for *µ* ≥ 1 and (*u*, *v*) ∈ *∂PR*1 . In addition, if *R*1 > *m* min *t*∈[ 1 4 , 3 4 ] *σ*(*t*) := *m σ* ∗ , then by the definition of cone *P*, one can get that for any *t* ∈ [ 1 , 3 ] and (*u*, *v*) ∈ *∂PR*1 ,

$$
\mu(t) \ge \min\_{t \in \left[\frac{1}{4}, \frac{3}{4}\right]} \mu(t) \ge \sigma^\* R\_1 > m. \tag{11}
$$

So, by (7), (11), and Lemma 1, one can get that for all (*u*, *v*) ∈ *∂PR*1 ,

k*T*(*u*, *v*)k ≥ *T*1(*u*, *v*)( 1 2 ) = Z 1 0 Z 1 0 *G*1,1( 1 2 , *τ*)*G*1,2(*τ*,*s*)*f*1(*s*, *u*(*s*), *v*(*s*))*dsdτ* ≥ *δ*1,1*δ*1,2*G*1,2( 1 2 , 1 2 )*C*1 Z 1 0 *G*1,2(*s*,*s*)*f*1(*s*, *u*(*s*), *v*(*s*))*ds* ≥ *δ*1,1*δ*1,2*G*1,2( 1 2 , 1 2 )*C*1(1 + *ε*)*λ*1 Z 3 4 1 4 *G*1,2(*s*,*s*)*u*(*s*)*ds* ≥ *δ*1,1*δ*1,2*G*1,2( 1 2 , 1 2 )*C*1(1 + *ε*)*λ*1*R*1*σ* ∗ Z 3 4 1 4 *G*1,2(*s*,*s*)*ds* ≥ *δ*1,1*δ*1,2*G*1,2( 1 2 , 1 2 )*C*1(1 + *ε*)*λ*1*m* Z 3 4 1 4 *G*1,2(*s*,*s*)*ds* > 0.

That is, inf (*u*,*v*)∈*∂PR*1 k*T*(*u*, *v*)k > 0. So, we can ultimately choose *R*1 > max{*R* ∗ 1 ,*r*1, *m σ* ∗ } such that

(9) holds.

Based on (5), (9), Lemma 2, and Lemma 3, we have

$$i(T\_\prime P\_{\overline{R}\_1} \backslash \overline{P}\_{r\_1 \prime} P) = i(T\_\prime P\_{\overline{R}\_1 \prime} P) - i(T\_\prime P\_{r\_1 \prime} P) = 0 - 1 = -1.$$

As a result, the conclusion of this theorem follows.

**Theorem 2.** *Assume that (H1) and (H3) hold. Then the BVP (1) has at least one nontrivial solution.*

**Proof of Theorem 2.** In the following, we divide the proof of Theorem 2 into three steps. *Step* 1. From condition (H3), there exist *ε* > 0 and *r*2 > 0 such that

$$f\_1(t, u, v) \ge (1 + \varepsilon)\lambda\_1 u\_\prime \,\,\forall t \in [0, 1], \,\,(u, v) \in \Lambda\_{[0, r\_2]}.\tag{12}$$

Subsequently, we claim that

$$\inf\_{(u,v)\in\partial P\_{r\_2}} \|T(u,v)\| > 0 \text{ and } \mu T(u,v) \neq (u,v), \; \forall \mu \ge 1, \; (u,v) \in \partial P\_{r\_2}.\tag{13}$$

In fact, if there exist *µ*0 ≥ 1 and (*u*0, *v*0) ∈ *∂Pr*2 such that *µ*0*T*(*u*0, *v*0) = (*u*0, *v*0), then by (6) and (12), one can obtain immediately

$$
\mu\_0^{(4)}(t) + \beta\_1 \mu\_0^{\prime\prime}(t) - \mu\_1 \mu\_0(t) \ge f\_1(t, \mu\_0(t), \nu\_0(t)) \ge (1 + \varepsilon)\lambda\_1 \mu\_0(t).
$$

Hence,

$$\int\_0^1 \lambda\_1 u\_0(t) \sin(\pi t) dt \ge (1+\varepsilon) \int\_0^1 \lambda\_1 u\_0(t) \sin(\pi t) dt.$$

Noticing that Z 1 0 *λ*1*u*0(*t*) sin(*πt*)*dt* > 0, we get a contradiction. In addition, it follows from Lemma 1 and (12) that for (*u*, *v*) ∈ *∂Pr*2 ,

$$\begin{split} \|T(u,v)\| &\geq T\_{1}(u,v)(\frac{1}{2}) \\ &= \int\_{0}^{1} \int\_{0}^{1} G\_{1,1}(\frac{1}{2},\tau) G\_{1,2}(\tau,s) f\_{1}(s,u(s),v(s)) ds d\tau \\ &\geq \delta\_{1,1} \delta\_{1,2} G\_{1,2}(\frac{1}{2},\frac{1}{2}) \mathbb{C}\_{1} \int\_{0}^{1} G\_{1,2}(s,s) f\_{1}(s,u(s),v(s)) ds \\ &\geq \delta\_{1,1} \delta\_{1,2} G\_{1,2}(\frac{1}{2},\frac{1}{2}) \mathbb{C}\_{1} (1+\varepsilon) \lambda\_{1} \int\_{0}^{1} G\_{1,2}(s,s) u(s) ds \\ &\geq \delta\_{1,1} \delta\_{1,2} G\_{1,2}(\frac{1}{2},\frac{1}{2}) \mathbb{C}\_{1} (1+\varepsilon) \lambda\_{1} r\_{2} \int\_{0}^{1} G\_{1,2}(s,s) \sigma(s) ds > 0, \end{split}$$

which yields inf (*u*,*v*)∈*∂Pr*2 k*T*(*u*, *v*)k > 0.

*Step* 2. The assumption (H3) implies that there exist *ε* ∈ (0, 1) and *m* > 0 such that

$$f\_1(t, \mu, v) \le (1 - \varepsilon)\lambda\_1 \mu\_\prime \,\,\forall t \in [0, 1], \,\,(\mu, v) \in \Lambda\_{[m, +\infty)}.\tag{14}$$

Moreover, by the continuity of *f*1 and *f*2, there exists *C* ∗ > 0 such that

$$f\_1(t, \mu, v) \le (1 - \varepsilon)\lambda\_1 \mu + \mathbb{C}^\* \text{ } \forall t \in [0, 1] \text{ } (\mu, v) \in \Lambda\_{\mathbb{R}^+}.\tag{15}$$

We claim that there exists a large enough *R*2 > *r*2 such that

$$
\mu T(\mu, \upsilon) \neq (\mu, \upsilon), \ \forall \mu \in (0, 1], \ (\mu, \upsilon) \in \partial P\_{\mathbb{R}\_2}.\tag{16}
$$

Suppose, on the contrary, there exist *µ*0 ∈ (0, 1] and (*u*0, *v*0) ∈ *∂PR*2 such that *µ*0*T*(*u*0, *v*0) = (*u*0, *v*0). Then (6) together with (15) guarantees

$$
\mu\_0^{(4)}(t) + \beta\_1 \mu\_0''(t) - \mathfrak{a}\_1 \mathfrak{u}\_0(t) \le f\_1(t, \mathfrak{u}\_0(t), \mathfrak{v}\_0(t)) \le (1 - \varepsilon)\lambda\_1 \mathfrak{u}\_0(t) + \mathcal{C}^\*.
$$

Consequently,

$$\int\_0^1 \lambda\_1 u\_0(t) \sin(\pi t) dt \le (1 - \varepsilon) \int\_0^1 \lambda\_1 u\_0(t) \sin(\pi t) dt + \frac{2\mathcal{C}^\*}{\pi} \lambda$$

namely,

$$\int\_0^1 u\_0(t) \sin(\pi t) dt \le \frac{2C^\*}{\pi e \lambda\_1}.$$

Moreover, based on the definition of cone *P*, we can immediately get

$$\|\|u\_0\|\int\_0^1 \sigma(t)\sin(\pi t)dt \le \int\_0^1 u\_0(t)\sin(\pi t)dt \le \frac{2C^\*}{\pi \varepsilon \lambda\_1}.$$

*Symmetry* **2020**, *12*, 1989

which means

$$\|\mathfrak{u}\_0\| \le \frac{2\mathcal{C}^\*}{\pi \varepsilon \lambda\_1 \int\_0^1 \sigma(t) \sin(\pi t) dt} \coloneqq R\_2^\*. \tag{17}$$

So, one can choose *R*2 > max{*R* ∗ 2 ,*r*2} such that (16) holds. *Step* 3. From (13), (16), Lemma 2, and Lemma 3, we deduce that

$$i(T\_\prime P\_{R\_2} \backslash \tilde{P}\_{r\_{2\prime}}P) = i(T\_\prime P\_{R\_{2\prime}}P) - i(T\_\prime P\_{r\_{2\prime}}P) = 1 - 0 = 1.1$$

As a result, BVP(1) has at least one nontrivial solution.

Up to now, some existence results of BVP(1) have been obtained by applying the fixed point index theory. In the following, the multiple solutions will be considered for BVP (1).

**Theorem 3.** *Assume that (H1) holds. In addition, suppose that*

*(1)* lim |*v*|≤*Nu u*→0+ sup max *t*∈[0,1] *f*1 (*t*,*u*,*v*) *u* < *λ*1*,* lim |*v*|≤*Nu u*→+∞ sup max *t*∈[0,1] *f*1 (*t*,*u*,*v*) *u* < *λ*1*;*

*(2) There exists r* > 0 *and a continuous nonnegative function* Φ*r such that*

$$f\_1(t, \mu, v) \ge \Phi\_r(t), \forall (t, \mu, v) \in [0, 1] \times (\sigma(t)r, r) \times [-Nr, Nr]^\perp$$

*and*

$$\max\_{t \in [0,1]} \int\_0^1 \int\_0^1 G\_{1,1}(t,\tau) G\_{1,2}(\tau,s) \Phi\_r(s) ds d\tau > r.c.$$

*Then the BVP (1) has at least two nontrivial solutions.*

**Proof of Theorem 3.** In order to obtain this conclusion, we firstly claim that

$$\inf\_{(u,v)\in\partial P\_r} \|T(u,v)\| > 0 \text{ and } \mu T(u,v) \neq (u,v), \forall \,\mu \ge 1, \ (u,v) \in \partial P\_r. \tag{18}$$

Suppose, on the contrary, there exist *µ*0 ≥ 1 and (*u*0, *v*0) ∈ *∂Pr* such that *µ*0*T*(*u*0, *v*0) = (*u*0, *v*0). Then,

$$\begin{split} \|u\_{0}\| &\quad \geq \|T(u\_{0},v\_{0})\| \geq T\_{1}(u\_{0},v\_{0})(t) \\ &= \int\_{0}^{1} \int\_{0}^{1} \mathbf{G}\_{1,1}(t,\tau) \mathbf{G}\_{1,2}(\tau,s) f\_{1}(s,u(s),v(s)) ds d\tau \\ &\geq \int\_{0}^{1} \int\_{0}^{1} \mathbf{G}\_{1,1}(t,\tau) \mathbf{G}\_{1,2}(\tau,s) \Phi\_{\tau}(s) ds d\tau. \end{split} \tag{19}$$

Taking the maximum for both sides of the above inequality in *t* ∈ [0, 1], we get that

$$\|\mu\_0\| \ge \max\_{t \in [0,1]} \int\_0^1 \int\_0^1 \mathbb{G}\_{1,1}(t,\tau) \mathbb{G}\_{1,2}(\tau,s) \Phi\_r(s) ds d\tau > r. \tag{20}$$

This means (*u*0, *v*0)∈ *∂Pr* , which is a contradiction. Moreover, one can easily see that inf (*u*,*v*)∈*∂Pr* k*T*(*u*, *v*)k > 0 holds from (19) and (20).

Next, similar to the process of proving (5) and (16), there exist *r*1 ∈ (0,*r*) and *R*2 ≥ max{*R* ∗ 2 ,*r*2,*r*} such that

$$
\mu T(u, v) \neq (u, v), \forall \mu \in (0, 1], \; \forall (u, v) \in \partial P\_{r\_1 \nu} \tag{21}
$$

$$
\mu T(u, v) \neq (u, v), \; \forall \mu \in (0, 1], \; \forall (u, v) \in \partial P\_{R\_2}.\tag{22}
$$

Thus, by (18), (21), (22), Lemma 2, and Lemma 3, one can immediately obtain that

$$i(T\_\prime P\_{R\_2} \backslash \tilde{P}\_{r\prime} P) = i(T\_\prime P\_{R\_2\prime} P) - i(T\_\prime P\_{r\prime} P) = 1 - 0 = 1\,,$$

$$i(T\_\prime P\_{r\_\perp} \backslash P\_{r\_{1'}}P) = i(T\_\prime P\_{r\_\prime}P) - i(T\_\prime P\_{r\_{1'}}P) = 0 - 1 = -1.$$

Namely, there exist (*u*1, *v*1) ∈ *Pr* \ *P*¯ *r*1 and (*u*2, *v*2) ∈ *PR*2 \ *P*¯ *r* satisfying *T*(*ui* , *vi*) = (*ui* , *vi*)(*i* = 1, 2), that is, (*ui* , *vi*)(*i* = 1, 2) is the solution of BVP(1).

Finally, we show (*u*1, *v*1) 6= (*u*2, *v*2). To see this we need only to prove BVP(1) has no solution on *∂Pr* . Suppose on the contrary, there exists (*u* ∗ , *v* ∗ ) ∈ *∂Pr* being a solution of BVP(1). Then *T*(*u* ∗ , *v* ∗ ) = (*u* ∗ , *v* ∗ ). By a similar process of obtaining (20), one can get k*u* ∗k = k*T*1(*u* ∗ , *v* ∗ )k > *r*, which is a contradiction. To sum up, Theorem 3 is proved.

From a process similar to the above, the following conclusion can be obtained.

**Theorem 4.** *Suppose that (H1) holds. In addition, suppose that*

*(1)* lim |*v*|≤*Nu u*→0+ inf min *t*∈[0,1] *f*1 (*t*,*u*,*v*) *u* > *λ*1*,* lim |*v*|≤*Nu u*→+∞ inf min *t*∈[0,1] *f*1 (*t*,*u*,*v*) *λ*1*u* > *λ*1*;*

*(2) There exists R* > 0*, and a continuous nonnegative function* Ψ*R such that*

$$f\_1(t, \iota, v) \le \Psi\_R(t), \forall (t, \iota, v) \in [0, 1] \times [\sigma(t)\mathbb{R}, \mathbb{R}] \times [-N\mathbb{R}, N\mathbb{R}]$$

and

$$\max\_{t \in [0,1]} \int\_0^1 \int\_0^1 G\_{1,1}(t,\tau) G\_{1,2}(\tau,s) \Psi\_R(s) ds d\tau < R.$$

*Then the BVP (1) has at least two nontrivial solutions.*

#### **Proof of Theorem 4.** We firstly prove that

$$
\mu T(\mu, v) \neq (\mu, v), \,\forall \mu \in (0, 1], \,(\mu, v) \in \partial P\_{\mathbb{R}}.\tag{23}
$$

To this end, suppose on the contrary that there exist *µ*0 ∈ (0, 1] and (*u*0, *v*0) ∈ *∂PR* such that *µ*0*T*(*u*0, *v*0) = (*u*0, *v*0). Hence, we get *u*0 = *µ*0*T*1(*u*0, *v*0), that is

$$u\_0(t) \le T\_1(u\_0, v\_0)(t) \quad \le \int\_0^1 \int\_0^1 G\_{1,1}(t, \tau) G\_{1,2}(\tau, s) \Psi\_R(s) ds d\tau < \mathcal{R}.\tag{24}$$

Noticing that (*u*0, *v*0) ∈ *∂PR*, this is a contradiction.

Next, from a process similar to (9) and (13), there exist *R*1 > max{*R*, *R* ∗ 1 ,*r*1, *m σ* ∗ } and *r*2 ∈ (0, *R*) such that

$$\inf\_{(u,v)\in\partial P\_{\mathbb{R}\_1}} \|T(u,v)\| > 0 \text{ and } \mu T(u,v) \neq (u,v), \,\forall \mu \ge 1, \,(u,v) \in \partial P\_{\mathbb{R}\_1}.\tag{25}$$

$$\inf\_{(u,v)\in\partial P\_{\mathcal{I}\_2}} \|T(u,v)\| > 0 \text{ and } \mu T(u,v) \neq (u,v), \; \forall \mu \ge 1, \; (u,v) \in \partial P\_{\mathcal{I}\_2}.\tag{26}$$

So, by (23)–(26), Lemma 2, and Lemma 3, one can get

$$i(T\_\prime P\_{R\_1} \backslash \tilde{P}\_{R\prime} P) = i(T\_\prime P\_{R\_1\prime} P) - i(T\_\prime P\_{R\prime} P) = 0 - 1 = -1,$$

$$i(T\_\prime P\_R \backslash \tilde{P}\_{r\_{2\prime}} P) = i(T\_\prime P\_{R\prime} P) - i(T\_\prime P\_{r\_{2\prime}} P) = 1 - 0 = 1.$$

Finally, from a process similar to the end of proof of Theorem 3, BVP(1) has at least two nontrivial solutions. As a result, the conclusion of this theorem follows.

#### **4. Examples**

In this section, two illustrative examples are worked out to show the effectiveness of the obtained results.

**Example 1.** *Consider the following BVP of fourth-order ordinary differential systems*

$$\begin{cases} u^{(4)}(t) + u''(t) - \pi^2 u(t) = f\_1(t, u, v), \ 0 < t < 1; \\ v^{(4)}(t) + \frac{1}{2} v''(t) - \frac{\pi^2}{2} v(t) = f\_2(t, u, v), \ 0 < t < 1; \\ u(0) = u(1) = u''(0) = u''(1) = 0; \\ v(0) = v(1) = v''(0) = v''(1) = 0, \end{cases} \tag{27}$$

where

$$f\_1(t,u,v) = \begin{cases} \begin{array}{llll} (\pi^4 - 2\pi^2)(1+\sin(\pi t))uv^{\frac{1}{4}} & \text{if } 0 < u < 1, \ |v| < u; \\ (\pi^4 - 2\pi^2)(1+\sin(\pi t))v^{\frac{1}{4}} & \text{if } u = 1, \ |v|  1, \ |v| < u; \end{cases} \\\\ f\_2(t,u,v) = \begin{cases} (\pi^4 - 2\pi^2)(1+\cos(\pi t))uv^{\frac{1}{4}} & \text{if } 0 < u < 1, \ |v| < u; \\ (\pi^4 - 2\pi^2)(1+\cos(\pi t))v^{\frac{1}{4}} & \text{if } u = 1, \ |v|  1, \ |v| < u. \end{cases}$$

Then, BVP (27) has at least two nontrivial solutions.

**Proof of Example 1 .** BVP (27) can be regarded as a BVP of the form (1). Choosing *α*1 = *π* 2 , *β*1 = 1, and *λ*1 = *π* 4 − 2*π* 2 > 0, then we have

$$\xi\_{1,1}^{\star} = \frac{-\beta\_1 + \sqrt{\beta\_1^2 + 4a\_1}}{2} = \frac{-1 + \sqrt{1 + 4\pi^2}}{2}, \\ \xi\_{1,2}^{\star} = \frac{-\beta\_1 - \sqrt{\beta\_1^2 + 4a\_1}}{2} = \frac{-1 - \sqrt{1 + 4\pi^2}}{2}.$$

Clearly, *α*1 and *β*1 satisfy the condition (2). Moreover, by careful calculation and Lemma 2.1 in Reference [32], one can obtain that

$$G\_{1,1}(t,s) = \begin{cases} \frac{\sinh \, w\_{1,1}t \sinh \, w\_{1,1}(1-s)}{w\_{1,1} \sinh \, w\_{1,1}} & 0 \le t \le s \le 1;\\ \frac{\sinh \, w\_{1,1}s \sinh \, w\_{1,1}(1-t)}{w\_{1,1} \sinh \, w\_{1,1}} & 0 \le s \le t \le 1; \end{cases}$$

$$G\_{1,2}(t,s) = \begin{cases} \frac{\sin \, w\_{1,2}t \sin \, w\_{1,2}(1-s)}{w\_{1,2} \sin \, w\_{1,2}} & 0 \le t \le s \le 1;\\ \frac{\sin \, w\_{1,2}s \sin \, w\_{1,2}(1-t)}{w\_{1,2} \sin \, w\_{1,2}} & 0 \le s \le t \le 1; \end{cases}$$

where *w*1,*i* = p |*ξ*1,*i* |(*i* = 1, 2).

Now, |*v*| ≤ 2*u*, | *f*2(*t*, *u*, *v*)| ≤ 2| *f*1(*t*, *u*, *v*)|, and *N* = *N*1*N*2*N*3. Thus, one can easily get that (H1) holds by choosing *N*3 ≥ max{2, 2 *N*1*N*2 }, where *Nj* = sup 0<*t*,*s*<1 *G*2,*j* (*t*,*s*) *G*1,*j* (*t*,*s*) , *j* = 1, 2.

In addition, by calculation, we get that

$$\lim\_{\substack{|v| \le Nu \\ u \to 0^+}} \sup\_{t \in [0,1]} \frac{f\_1(t, u, v)}{u} = \lim\_{\substack{|v| \le Nu \\ u \to 0^+}} \sup\_{t \in [0,1]} \frac{(\pi^4 - 2\pi^2)(1 + \sin(\pi t)) u v^{\frac{1}{4}}}{u} = 0 < \lambda\_1 \lambda\_2$$

$$\lim\_{\substack{|v| \le \mathbb{N}u \\ u \to +\infty}} \sup\_{t \in [0,1]} \frac{f\_1(t, u, v)}{u} = \lim\_{\substack{|v| \le \mathbb{N}u \\ u \to +\infty}} \sup\_{t \in [0,1]} \frac{(\pi^4 - 2\pi^2)(1 + \sin(\pi t))u^{\frac{1}{4}}v^{\frac{1}{4}}}{u} = 0 < \lambda\_1.$$

Choose

$$\tau = \min\{1, \left[\delta\_{1,1}\delta\_{1,2}\int\_0^1 \sqrt{\sigma(t)}\sin(\pi t)dt \min\_{t \in \left[\frac{1}{2}, \frac{3}{4}\right]} \left(G\_{1,1}(t,t)G\_{1,2}(t,t)\right)^2\right] > 0\}$$

and

$$
\Phi\_r(t) = \sqrt{r\sigma(t)}\sin(\pi t).
$$

Then, it is not difficult to obtain that the condition (2) in Theorem 3 holds. Hence, our conclusion follows from Theorem 3.

**Proof of Example 2 .** BVP (28) can be regarded as a BVP of the form (1). Using a similar process of the proof of Example 1, one can easily obtain that

$$\lim\_{\substack{|v| \le Nu \\ u \to 0^{+}}} \inf \min\_{t \in [0,1]} \frac{f\_1(t, u, v)}{u} = \lim\_{\substack{|v| \le Nu \\ u \to 0^{+}}} \inf \min\_{t \in [0,1]} \frac{(\pi^4 - 2\pi^2)(2+t)u^{\frac{1}{2}}v^{\frac{1}{2}}}{u} = +\infty > \pi^4 - 2\pi^2 = \lambda\_1 > 0$$

$$\lim\_{\substack{|v| \le Nu \\ u \to +\infty}} \inf \min\_{t \in [0,1]} \frac{f\_1(t, u, v)}{u} = \lim\_{\substack{|v| \le Nu \\ u \to +\infty}} \inf \min\_{t \in [0,1]} \frac{(\pi^4 - 2\pi^2)(2+t)u^{2}v^{\frac{1}{2}}}{u} + \infty > \pi^4 - 2\pi^2 = \lambda\_1 > 0$$

In addition, it is obvious that (H1) holds by choosing *N*3 = 2. In the following, set

$$R = \max\{1, \frac{2}{5\pi^2 \mathcal{C}\_{1,1} \mathcal{C}\_{1,2}} \max\_{t \in [0,1]} [\mathcal{G}\_{1,1}(t,t)\mathcal{G}\_{1,2}(t,t)]\} > 0$$

and

$$
\Psi\_R(t) = \pi^4 R^2 (2+t).
$$

Then, it is trivial to verify that assumption (2) of Theorem 3 is true. As a result, by Theorem 4, system (28) has at least two nontrivial solutions.

#### **5. Conclusions**

In this paper, we have obtained some appropriate results corresponding to multiple solutions for a class of nonlinear fourth-order boundary value problems with parameters. The multiple solutions for the considered systems are obtained under some suitable assumptions via fixed point index theory. The whole theoretical results has been demonstrated by providing two interesting examples. Hence, we claim that fixed point index theory can be used as a strong technique to study nonlinear fourth-order boundary value problems with parameters.

**Author Contributions:** Conceptualization and Visualization, Y.L.; Formal analysis and Investigation, D.Z.; Writing original draft and Investigation, L.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by NNSF of P.R.China (62073202 and 11671237), and a project of Shandong Province Higher Educational Science and Technology Program of China under the grant J18KA233.

**Acknowledgments:** The authors are thankful to the editor and anonymous referees for their valuable comments and suggestions.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Existence of Three Solutions for a Nonlinear Discrete Boundary Value Problem with** *φc***-Laplacian**

**Yanshan Chen 1,2 and Zhan Zhou 1,2,**∗

Received: 14 October 2020; Accepted: 3 November 2020; Published: 6 November 2020

**Abstract:** In this paper, based on critical point theory, we mainly focus on the multiplicity of nontrivial solutions for a nonlinear discrete Dirichlet boundary value problem involving the mean curvature operator. Without imposing the symmetry or oscillating behavior at infinity on the nonlinear term *f* , we respectively obtain the sufficient conditions for the existence of at least three non-trivial solutions and the existence of at least two non-trivial solutions under different assumptions on *f* . In addition, by using the maximum principle, we also deduce the existence of at least three positive solutions from our conclusion. As far as we know, our results are supplements to some well-known ones.

**Keywords:** *φc*-Laplacian; boundary value problem; critical point theory; three solutions

#### **1. Introduction**

Let Z and R denote all integers and real numbers, respectively. Let *N* be a fixed positive integer. Define Z(*a*, *b*) = {*a*, *a* + 1, · · · , *b*} with *a* ≤ *b* for any *a*, *b* ∈ Z.

Difference equations are widely used in various research fields, such as computer science, discrete optimization, economics and biological neural networks [1–4]. On the existence and multiplicity of solutions for the boundary value problems of difference equations, many authors have come to important conclusions by exploiting various methods, including the method of upper and lower solutions, Brouwer degree and invariant sets of descending flow [5–7]. Critical point theory was used largely to explore differential equations much earlier on in history. In 2003, Guo and Yu in [8] used critical point theory for the first time to obtain sufficient conditions for the existence of periodic solutions and subharmonic solutions of difference equations. This crucial breakthrough inspired many scholars to use critical point theory to study the dynamics of difference equations and many meaningful and interesting results have been obtained, especially in periodic solutions [9–11], homoclinic solutions [12–16] and boundary value problems [17–23].

In [24], Agarwal, Perera and O'Regan employed the critical point theory to establish the existence of at least two positive solutions of the following second order discrete boundary value problem

$$\begin{cases} -\Delta^2 \mathfrak{x}(k-1) = f(k, \mathfrak{x}(k)), \; k \in \mathbb{Z}(1, N),\\ \mathfrak{x}(0) = \mathfrak{x}(N+1) = 0, \end{cases}$$

where ∆*x*(*k*) = *x*(*k* + 1) − *x*(*k*), ∆ 2*x*(*k*) = ∆(∆*x*(*k*)) and *f* : Z(1, *N*) × R → R is a continuous function. In [18], by using the three critical point theorems proposed by Bonanno [25], Jiang and Zhou obtained

sufficient conditions for the existence of at least three solutions of the following Dirichlet boundary value problem with *φp*-Laplacian

$$\begin{cases} -\Delta(\phi\_p(\Delta u(k-1))) = \lambda f(k, u(k)), \; k \in \mathbb{Z}(1, N), \\ u(0) = u(N+1) = 0, \end{cases}$$

where *φp* is the *p*-Laplacian defined by *φp*(*u*) = |*u*| *p*−2*u* (*p* > 1) with *u* ∈ R and *λ* > 0 is a positive parameter. Different from the conclusion of [18], Bonanno in [26] obtained the existence of three positive solutions without the asymptotic condition of the nonlinear function *f* . In particular, Bonanno obtained the sufficient conditions for the existence of at least four nontrivial solutions when *f* satisfies the growth condition at zero and infinity [26], which improved the result in [18].

In [27], by using critical point theory, Nastasi and Vetro obtained the existence of at least two positive solutions to the following Dirichlet boundary value problem with (*p*, *q*)-Laplacian

$$\begin{cases} -\Delta(\phi\_p(\Delta u(z-1))) - \Delta(\phi\_\emptyset(\Delta u(z-1))) + a(z)\phi\_p(u(z)) + \beta(z)\phi\_\emptyset(u(z)) = \lambda g(z, u(z)), \; z \in \mathbb{Z}(1, N), \\ u(0) = u(N+1) = 0, \end{cases}$$

where 1 < *q* < *p* < +∞, *α*, *β* : Z(1, *N*) → R and *g* : [1, *N* + 1] × R → R are a continuous function with *g*(*N* + 1, *t*) = 0 for all *t* ∈ R.

It is well known that the differential equation with *φc*-Laplacian (*φc* is the mean curvature operator defined by *φc*(*s*) = *s*/ √ 1 + *s* 2 for *s* ∈ R) was studied by many scholars in the past decades [28–31]. It is usually regarded as a variant of the Liouville–Bratu–Gelfand problem, which is used to study the dynamic model of combustible gases. As mentioned above, we find that the research on the difference equation largely focuses on the case with *φp*-Laplacian. However, there are only a few results on the boundary value problems involving *φc*-Laplacian [32]. Recently, Zhou and Ling in [33] considered the existence of multiple solutions of the following discrete Dirichlet boundary value problem with *φc*-Laplacian.

Problem 1

$$\begin{cases} -\Delta \left( \frac{\Delta u(k-1)}{\sqrt{1+(\Delta u(k-1))^2}} \right) = \lambda f(k, u(k)), \; k \in \mathbb{Z}(1, N), \\ u(0) = u(N+1) = 0. \end{cases}$$

The authors in [33] found that the properties of the nonlinear term *F*(*k*, *u*) = R *u* 0 *f*(*k*, *t*)*dt* plays an important role on the existence of multiple solutions. When *F*(*k*, *u*) has oscillating behavior at infinity, there are infinite solutions to the boundary value (Problem 1) [33]. Naturally, we would like to ask: What will happen if *F*(*k*, *u*) does not oscillate at infinity?

To address this problem, in this paper, we will study the existence of solutions for the boundary value (Problem 1) without oscillating nonlinear terms. In fact, based on the theorems of G. Bonanno (Theorem 4.1 in [22] and Theorem 2.1 in [23]), we will give the conditions of the existence of at least three nontrivial solutions for (Problem 1), when *F*(*k*, *u*) does not have oscillation property at infinity. In addition, when *f*(*k*, 0) > 0, we obtain the sufficient conditions for the existence of at least three positive solutions of (Problem 1). Moveover, we give two examples to illustrate our main results.

For convenience, we end this section by recalling some classical definitions and two well-known lemmas, which are the main tools of this paper.

Let (*X*, k · k) be a real Banach space. We say that *I* : *X* → R is coercive on *X* if limk*u*k→+∞ *I*(*u*) = +∞. If *I* is a continuously Gâteaux differentiable functional, we say that *I* satisfies the Palais–Smale condition ((PS)-condition in short), if any sequence *un* ⊂ *X* such that {*I*(*un*)} is bounded and {*I* 0 (*un*)} is convergent to 0 in *X* ∗ , has a convergent subsequence in *X*.

(*H*) *Let* (*X*, k · k) *be a real finite dimensional Banach space and let* Φ*,* Ψ : *X* → R *be two continuously Gâteaux differentiable functionals with* Φ *coercive and such that*

$$\inf\_{X} \Phi = \Phi(0) = \Psi(0) = 0.$$

**Lemma 1** ([22])**.** *Assume that (H) holds and there exist r* > 0 *and x* ∈ *X, with* 0 < *r* < Φ(*x*)*, such that:*

$$\begin{aligned} (a\_1) \quad &\frac{\sup\_{\Phi(\mathbf{x}) \le r} \Psi(\mathbf{x})}{r} < \frac{\Psi(\overline{\mathbf{x}})}{\Phi(\overline{\mathbf{x}})};\\ (a\_2) \quad &\text{for each } \lambda \in \Lambda\_r := \left(\frac{\Phi(\overline{\mathbf{x}})}{\Psi(\overline{\mathbf{x}})}, \frac{r}{\sup\_{\Phi(\mathbf{x}) \le r} \Psi(\mathbf{x})}\right), \text{the functional } \Phi - \lambda \Psi \text{ is coercive.} \end{aligned}$$

*Then, for each λ* ∈ Λ*r , the functional Iλ* = Φ − *λ*Ψ *has at least three distinct critical points in X.*

**Lemma 2** ([23])**.** *Let X be a real Banach space and let* Φ*,* Ψ *: X* → R *be two continuously Gâteaux differentiable functionals such that* inf *X* Φ = Φ(0) = Ψ(0) = 0. *Assume that there are r* > 0 *and u*˜ ∈ *X, with* 0 < Φ(*u*˜) < *r, such that*

$$\frac{\sup\_{\Phi(u)\le r} \Psi(u)}{r} < \frac{\Psi(\mathfrak{u})}{\Phi(\tilde{u})} \tag{1}$$

*for each λ* ∈ Φ(*u*˜) Ψ(*u*˜) , *r* sup Φ(*u*)≤*r* Ψ(*u*) *, and the functional Iλ* = Φ − *λ*Ψ *satisfies (PS)-condition and it is unbounded from below. Then, for each λ* ∈ Φ(*u*˜) Ψ(*u*˜) , *r* sup Φ(*u*)≤*r* Ψ(*u*) *, the functional Iλ admits at least two non-zero critical points uλ*,1, *uλ*,2 *such that Iλ*(*uλ*,1) < 0 < *Iλ*(*uλ*,2).

#### **2. Preliminaries**

In this section, we recall some definitions, notations and properties. Consider the *N*-dimensional Banach space

$$S = \{ \mu : [0, N+1] \to \mathbb{R} : \mu(0) = \mathfrak{u}(N+1) = 0 \}$$

endowed with the norm

$$\|u\| = \left(\sum\_{k=1}^{N+1} |\Delta u(k-1)|^2\right)^{\frac{1}{2}}.$$

We define the functional Φ, Ψ : *S* → R in the following way

$$\Phi(u) = \sum\_{k=1}^{N+1} \left( \sqrt{1 + (\Delta u(k-1))^2} - 1 \right), \quad \Psi(u) = \sum\_{k=1}^{N} F(k, u(k)) \tag{2}$$

for each *u* ∈ *S*, where

$$F(k, \xi) := \int\_0^{\xi} f(k, t)dt, \quad \xi \in \mathbb{R}.$$

It is easy to check that Φ, Ψ ∈ *C* 1 (*S*, *R*) and we have

$$\Phi'(u)(v) = -\sum\_{k=1}^{N} \Delta(\phi\_c \left(\Delta u(k-1)\right)) v(k)$$

and

$$\Psi'(\boldsymbol{\mu})(v) = -\sum\_{k=1}^{N} f(k, \boldsymbol{\mu}(k)) v(k).$$

By a standard argument, it can be shown that the critical points of the functional

$$I\_{\lambda} = \Phi - \lambda \Psi$$

are the solutions of (Problem 1).

Let

$$\|\mu\|\_{\infty} = \max\left\{ |\mu(k)| : k \in \mathbb{Z}(1, N) \right\}.$$

We see that k · k∞ is another norm in *S*. From Lemma 2.2 of [18], we have

**Lemma 3** ([18])**.** *For any u* ∈ *S, the following relation holds*

$$\|u\|\_{\infty} \le \frac{\sqrt{N+1}}{2} \|u\|.$$

From (2.1) and (2.3) in [22], we have

**Lemma 4** ([22])**.** *For any u* ∈ *S, one has*

$$\frac{1}{\sqrt{N\lambda\_N}}||u|| \le ||u||\_{\infty} \le \frac{1}{\sqrt{\lambda\_1}}||u||\_{\prime\prime}$$

*where λ*1 = 4 sin2 *π* 2(*N*+1) *and λN* = 4 sin2 *Nπ* 2(*N*+1) *.*

Finally, in order to obtain the positive solutions of (Problem 1), we need the following strong maximum principle, which can be found in Theorem 2.1 of [33].

**Lemma 5** ([33])**.** *Assume u* ∈ *S such that either*

$$
\mu(k) > 0 \quad \text{or} \quad \Delta(\phi\_\mathfrak{c}(\Delta \mu(k-1))) \le 0,
$$

*for all k* ∈ Z(1, *N*). *Then, either u* > 0 *in* Z(1, *N*) *or u* ≡ 0*.*

#### **3. Main Results**

For convenience, set

$$F\_t := \sum\_{k=1}^{N} F(k, t)\_t \text{ for all } t > 0.$$

Our first result is the following theorem.

**Theorem 1.** *Assume that there exist two positive constants c and d with*

$$2\sqrt{1+d^2} > 1 + \sqrt{1 + \frac{4c^2}{N+1}}\tag{3}$$

*such that*

$$(i) \quad f(k, \xi) > 0 \text{ for each } k \in \mathbb{Z} \\ (1, N) \text{ and } \xi \in [-c, c]\varphi$$

$$\left(\left(\mathrm{ii}\right)\_{-}\right)\_{\frac{F\_d}{2\left(\sqrt{1+d^2}-1\right)}} > \frac{F\_c}{-1+\sqrt{1+\frac{4c^2}{N+1}}};$$

*(iii)* lim sup |*ξ*|→+∞ *F*(*k*,*ξ*) |*ξ*| < 2*Fc N*( √ 4*c* 2+*N*+1− √ *N*+1) *.*

$$\text{Then, for every } \lambda \in \Lambda := \left(\frac{2(\sqrt{1+d^2}-1)}{\mathbb{F}\_d}, \frac{-1+\sqrt{1+\frac{4r^2}{N+1}}}{\mathbb{F}\_d}\right), \text{(Problem 1) has at least three nontrivial solutions.}$$

**Proof.** We take *X* = *S*, Φ and Ψ as in (2). Clearly, Φ and Ψ are two continuously Gâteaux differentiable functionals. Now, we prove the coercivity of Φ. In fact, one has

$$\begin{aligned} \Phi(u) &= \sum\_{k=1}^{N+1} \left( \sqrt{1 + (\Delta u(k-1))^2} - 1 \right) \\ &\ge \left( \sum\_{k=1}^{N+1} \left[ 1 + (\Delta u(k-1))^2 \right] \right)^{\frac{1}{2}} - \sum\_{k=1}^{N+1} 1 \\ &\ge ||u|| - N - 1. \end{aligned}$$

This means lim k*u*k→+∞ Φ(*u*) = +∞ and verifies the coercivity of Φ. Moreover, by the definition of Φ and Ψ, we can obtain

$$\inf\_{X} \Phi = \Phi(0) = \Psi(0) = 0.$$

To summarize, condition (*H*) in Lemma 1 holds.

According to Lemma 1, it is clear that Theorem 1 holds if we can verify (*a*1) and (*a*2) of Lemma 1. Put

$$r = -1 + \sqrt{1 + \frac{4c^2}{N+1}}.$$

If Φ(*u*) ≤ *r*, let

$$v(k) = \sqrt{1 + (\Delta u(k))^2} - 1 \text{ for } k \in \mathbb{Z}(0, N).$$

Then, we have *N* ∑ *k*=0 *v*(*k*) ≤ *r* and

$$\sum\_{k=0}^{N} (\Delta u(k))^2 = \sum\_{k=0}^{N} \left( v(k)^2 + 2v(k) \right) \le \left( \sum\_{k=0}^{N} v(k) \right)^2 + 2\sum\_{k=0}^{N} v(k) \le r^2 + 2r = \frac{4c^2}{N+1}.$$

Thus, k*u*k∞ ≤ √ *N*+1 2 k*u*k ≤ *c* by Lemma 3.

By (i), we see that *F*(*k*, *ξ*) is increasing in *ξ* ∈ [−*c*, *c*]. Thus

$$\sup\_{\substack{\Phi(u)\le r\\r}} \frac{\Psi(u)}{r} \le \frac{\sup\_{\|u\|\_{\infty}\le c} \sum\_{k=1}^{N} F(k, u(k))}{-1 + \sqrt{1 + \frac{4c^2}{N+1}}} \le \frac{\sum\_{k=1}^{N} \max\_{|\tilde{\xi}| \le c} F(k, \tilde{\xi})}{-1 + \sqrt{1 + \frac{4c^2}{N+1}}} = \frac{F\_c}{-1 + \sqrt{1 + \frac{4c^2}{N+1}}}.\tag{4}$$

Then, it is easy to get

$$\sup\_{\Phi(u)\le r} \Psi(u) \\ \qquad \frac{1}{\lambda} < \frac{1}{\lambda}. \tag{5}$$

Now, let *u* ∈ *S* be defined by

$$
\overline{\mu}(k) = \begin{cases} d, & \text{if } k \in \mathbb{Z}(1, N), \\ 0, & \text{if } k = 0 \text{ or } k = N + 1. \end{cases}
$$

Then, we see from (3) that

$$\Phi(\overline{\pi}) = 2(\sqrt{1+d^2}-1) > \sqrt{1+\frac{4c^2}{N+1}} - 1 = r.s.$$

*Symmetry* **2020**, *12*, 1839

Moreover, it holds that

$$\frac{\Psi(\overline{u})}{\Phi(\overline{u})} = \frac{\sum\_{k=1}^{N} F(k, \overline{u}(k))}{2(\sqrt{1+d^2}-1)} = \frac{F\_d}{2(\sqrt{1+d^2}-1)}.\tag{6}$$

Therefore, we have

$$\frac{\Psi(\overline{u})}{\Phi(\overline{u})} > \frac{1}{\lambda}. \tag{7}$$

Hence, condition (*a*1) of Lemma 1 follows by combining (5) with (7). Next, we prove the coercivity of the functional *Iλ*. From (*iii*), there is an *ε* > 0 such that

$$\limsup\_{|\xi| \to +\infty} \frac{F(k, \xi)}{|\xi|} < \varepsilon < \frac{2F\_{\varepsilon}}{N(\sqrt{4c^2 + N + 1} - \sqrt{N + 1})}.$$

Then, there is a positive constant *h* such that

$$F(k, \xi) \le \varepsilon |\xi| + h$$

for each *ξ* ∈ R and *k* ∈ Z(1, *N*). By Lemma 3 and *λ* < −1+ q 1+ 4*c* 2 *N*+1 *Fc* , one has

$$\begin{aligned} \lambda \sum\_{k=1}^{N} F(k, u(k)) &\leq \lambda \sum\_{k=1}^{N} \left[ \varepsilon |u(k)| + h \right] \\ &\leq \lambda \frac{\varepsilon N \sqrt{N+1}}{2} ||u|| + \lambda Nh \\ &\leq \frac{\varepsilon N \left( \sqrt{N+1+4c^2} - \sqrt{N+1} \right)}{2F\_c} ||u|| + \frac{\left( -1 + \sqrt{1 + \frac{4c^2}{N+1}} \right) Nh}{F\_c} \end{aligned}$$

for each *u* ∈ *S*. Therefore, it is true that

$$\begin{split} I\_{\lambda}(u) &= \Phi(u) - \lambda \Psi(u) \\ &= \sum\_{k=1}^{N+1} \left[ \sqrt{1 + (\Delta u(k-1))^2} - 1 \right] - \lambda \sum\_{k=1}^{N} F(k, u(k)) \\ &\geq \left[ \sum\_{k=1}^{N+1} \left[ 1 + (\Delta u(k-1))^2 \right] \right]^{\frac{1}{2}} - N - 1 - \lambda \sum\_{k=1}^{N} F(k, u(k)) \\ &\geq \left[ 1 - \frac{\varepsilon N \left( \sqrt{N+1+4c^2} - \sqrt{N+1} \right)}{2F\_{\mathbb{C}}} \right] \|u\| - N - 1 - \frac{\left( -1 + \sqrt{1 + \frac{4c^2}{N+1}} \right) Nh}{F\_{\mathbb{C}}}. \end{split}$$

Thus, we get that *Iλ* is coercive, and condition (*a*2) of Lemma 1 is verified. In summary, all assumptions of Lemma 1 are proved, and so the functional *Iλ* has at least three distinct critical points in *X* for each *λ* ∈ Λ. Since *u* ≡ 0 is not a solution to the (Problem 1), the proof is completed.

**Remark 1.** *It is obvious that the mean curvature operator φc is odd symmetric (φc*(−*s*) = −*φc*(*s*)*, s* ∈ R*). When the nonlinear term f is also odd symmetric (f*(·, −*s*) = −*f*(·,*s*)*), the variational functional Iλ is even symmetric (Iλ*(−*u*) = *Iλ*(*u*)*, u* ∈ *S). In this case, it is easy to obtain multiple solutions to (Problem 1) by using the critical point theory with symmetries. However, in this paper, we obtain multiple solutions to (Problem 1) without the symmetry on f .*

Now, let

$$F^+(k, \mathfrak{f}) = \int\_0^\mathfrak{f} f(k, t^+) dt, \quad (k, \mathfrak{f}) \in \mathbb{Z}(1, N) \times \mathbb{R}\_+$$

where *t* + = max{0, *t*} and define *I* + *λ* = Φ − *λ*Ψ+, with Φ defined in (2) and

$$\Psi^+(\mathfrak{u}) := \sum\_{k=1}^N F^+(k, \mathfrak{u}(k)).$$

It is well known that *I* + *λ* ∈ *C* 1 (*S*, R) and the critical points of *I* + *λ* are precisely the solutions of the following problem.

Problem 2

$$\begin{cases} -\Delta(\phi\_c\left(\Delta u(k-1)\right)) = \lambda f(k, u^+(k)), \; k \in \mathbb{Z}(1, N),\\ u(0) = u(N+1) = 0. \end{cases}$$

We have the following corollary.

**Corollary 1.** *Assume that there exist two positive constants c and d with*

$$2\sqrt{1+d^2} > 1 + \sqrt{1 + \frac{4c^2}{N+1}}$$

*such that*

$$f(i) \quad f(k, \xi) > 0 \text{ for each } k \in \mathbb{Z}(1, N) \text{ and } \xi \in [0, c];$$

$$\left(\left(\mathrm{ii}\right)\_{-\frac{F\_d}{2\left(\sqrt{1+d^2}-1\right)}} > \frac{F\_c}{-1 + \sqrt{1 + \frac{4c^2}{N+1}}}\right;$$

$$(iii)\ \limsup\_{\xi \to +\infty} \frac{F(k,\xi)}{\xi} < \frac{2F\_0}{N(\sqrt{4c^2 + N + 1} - \sqrt{N + 1})}.$$

*Then, for every λ* ∈ Λ*, (Problem 1) has at least three positive solutions.*

**Proof.** For each *k* ∈ Z(1, *N*), consider (Problem 2) with

$$f^+(k,\xi) = \begin{cases} f(k,\xi), & \text{if } \xi > 0, \\ f(k,0), & \text{if } \xi \le 0. \end{cases}$$

One has that condition (*i*) of Theorem 1 holds. Besides, we have

$$\limsup\_{\mathfrak{f}\to-\infty} \frac{F^+(k,\mathfrak{f})}{|\mathfrak{f}|} = \limsup\_{\mathfrak{f}\to-\infty} \frac{\mathfrak{f}f(k,0)}{|\mathfrak{f}|} = -f(k,0) < \frac{2\mathcal{F}\_{\mathfrak{c}}}{N(\sqrt{4c^2+N+1}-\sqrt{N+1})} \cdot \frac{1}{2}$$

So, all conditions of Theorem 1 are true. Moreover, since *u* ≡ 0 is not a solution of problem (*P* + *λ*, *f* ), we can get that (Problem 2) has at least three nontrivial solutions. Assume *u* = {*u*(*k*)} is one of the nontrivial solution, for any *k* ∈ Z(1, *N*), one has either *u*(*k*) > 0 or

$$-\Delta(\phi\_{\varepsilon}\left(\Delta u(k-1)\right)) = \lambda f(k, u^+(k)) = \lambda f(k, 0) > 0.$$

Then, we have *u* > 0 for all *k* ∈ Z(1, *N*) by Lemma 5, i.e., *u* is a positive solution. In addition, if *u* is a positive solution of (Problem 2), then *u* is a positive solution of (Problem 1) obviously and Corollary 1 is proved.

Next, we will use Lemma 2 to obtain another conclusion of this paper.

**Theorem 2.** *Assume that there exist two positive constants c and d with*

$$2\sqrt{1+d^2} < 1 + \sqrt{1 + \frac{4c^2}{N+1}}\tag{8}$$

*such that*

$$\begin{aligned} (T\_1) \quad f(k, \xi) &> 0 \text{ for each } k \in \mathbb{Z}(1, N) \text{ and } \xi \in [-c, c];\\ (T\_2) \quad \max \left\{ \frac{2(\sqrt{1+d^2}-1)}{F\_d}, \frac{2\sqrt{N(N+1)}\sin\frac{N\pi}{2(N+1)}}{\beta} \right\} &< \frac{-1 + \sqrt{1 + \frac{4c^2}{N+1}}}{F\_c}; \end{aligned}$$

*(T*3*) there exist a positive constant β such that*

$$\liminf\_{|t| \to +\infty} \frac{F(k, t)}{|t|} > \beta$$

*for each k* ∈ Z(1, *N*).

$$\begin{aligned} \text{Then, for every } \lambda \in \left( \max\left\{ \frac{2(\sqrt{1+d^2}-1)}{\mathbb{F}\_d}, \frac{2\sqrt{N(N+1)}\sin\frac{N\pi}{2(N+1)}}{\mathbb{F}} \right\}, \frac{-1+\sqrt{1+\frac{4t^2}{N+1}}}{\mathbb{F}\_\mathbb{C}} \right), \text{ (Problem 1) has at } \lambda \in \left( \mathbb{F}\_\lambda(u\_{\lambda,1}), \begin{cases} \frac{2}{\sqrt{1+d^2}-1} \le \frac{2}{\sqrt{N(N+1)}}, \frac{2}{\sqrt{N(N+1)}} \le \frac{2}{\sqrt{N(N+1)}}, \\ 1-\frac{2}{\sqrt{1+d^2}-1} \le \frac{2}{\sqrt{N(N+1)}}, \end{cases} \right). \end{aligned}$$

**Proof.** Clearly, Φ, Ψ are two continuously Gâteaux differentiable functionals and

$$\inf\_{X} \Phi = \Phi(0) = \Psi(0) = 0.$$

Let *u* and *r* be the same as the ones defined in the proof of Theorem 1. Then we have

$$\frac{\sup\_{\Psi(u)\le r} \Psi(u)}{r} \le \frac{F\_c}{-1 + \sqrt{1 + \frac{4c^2}{N+1}}} < \frac{1}{\lambda}$$

and

$$\frac{\Psi(\overline{\mu})}{\Phi(\overline{\mu})} = \frac{F\_d}{2(\sqrt{1+d^2}-1)} > \frac{1}{\lambda}.$$

So inequality (1) in Lemma 2 holds. Besides, form (*T*3) and Lemma 4, there is a constant *h* such that, for any *u* ∈ *S* and *k* ∈ Z(1, *N*),

$$\begin{aligned} \lambda \sum\_{k=1}^N F(k, u(k)) &\geq \lambda \sum\_{k=1}^N (\beta |u(k)| - h) \\ &\geq \lambda \beta ||u||\_{\infty} - \lambda Nh \\ &\geq \frac{\lambda \beta}{\sqrt{N\lambda \lambda\_N}} ||u|| - \lambda Nh. \end{aligned}$$

#### Therefore, we have

$$\begin{split} \Lambda\_{\lambda}(u) &= \Phi(u) - \Psi(u) \\ &= \sum\_{k=1}^{N+1} \left[ \sqrt{1 + (\Delta u(k-1))^2} - 1 \right] - \lambda \sum\_{k=1}^{N} F(k, u(k)) \\ &\leq \sum\_{k=1}^{N+1} \left[ 1 + \sqrt{(\Delta u(k-1))^2} - 1 \right] - \frac{\lambda \beta}{\sqrt{N\lambda\lambda\_N}} ||u|| + \lambda Nh \\ &\leq \left( \sum\_{k=1}^{N+1} 1^2 \right)^{\frac{1}{2}} \left( \sum\_{k=1}^{N+1} (\Delta u(k-1))^2 \right)^{\frac{1}{2}} - \frac{\lambda \beta}{\sqrt{N\lambda\_N}} ||u|| + \lambda Nh \\ &\leq \left( \sqrt{N+1} - \frac{\lambda \beta}{\sqrt{N\lambda\_N}} \right) ||u|| + \lambda Nh \\ &= \left( \sqrt{N+1} - \frac{\lambda \beta}{2\sqrt{N}\sin\frac{N\pi}{2(N+1)}} \right) ||u|| + \lambda Nh. \end{split}$$

Thus, we see that lim k*u*k→+∞ *Iλ*(*u*) = −∞, which means that the functional *Iλ* is unbounded from below. Moreover, we can get that −*Iλ* is coercive. Therefore, *Iλ* satisfies the (PS)-condition and the proof is completed.

#### **4. Examples**

**Example 1.** *Fix N* ∈ Z(1, 20) *and consider the boundary value (Problem 1) with*

$$f(k,u) = f(u) = \begin{cases} 3u^2 + \frac{1}{4}\cos u, & \text{if } u \le \frac{3\pi}{2},\\ \frac{27\pi^2}{4}\cos(u + \frac{\pi}{2}), & \text{if } u > \frac{3\pi}{2}. \end{cases}$$

*for k* ∈ Z(1, *N*)*. Then, we have*

$$F(k,\xi) = F(\xi) = \begin{cases} \xi^3 + \frac{1}{4}\sin\xi, & \text{if } \xi \le \frac{3\pi}{2},\\ \frac{27\pi^3}{8} - \frac{1}{4} + \frac{27\pi^2}{4}\sin(\xi + \frac{\pi}{2}), & \text{if } \xi > \frac{3\pi}{2}. \end{cases}$$

*Let c* = 1 2 *and d* = 4*. Then one has*

$$2\sqrt{1+d^2} = 2\sqrt{17} > 1 + \sqrt{1+\frac{1}{N+1}} = 1 + \sqrt{1+\frac{4c^2}{N+1}}.$$

*In addition, as shown in Figure 1, f*(*k*, *ξ*) > 0 *for each ξ* ∈ [0, 1 2 ]*. Thus condition* (*i*) *of Corollary 1 follows. Moreover, we have*

$$\frac{F\_c}{-1 + \sqrt{1 + \frac{4e^2}{N+1}}} = \frac{N(1 + 2\sin\frac{1}{2})}{8\left(-1 + \sqrt{1 + \frac{1}{N+1}}\right)} \le \frac{N(1 + 2\sin\frac{1}{2})}{8\left(-1 + \sqrt{\frac{21}{20}}\right)} \approx 9.915N\tag{9}$$

*and*

$$\frac{F\_d}{2(\sqrt{1+d^2}-1)} = \frac{N(256+\sin 4)}{8(\sqrt{17}-1)} \approx 10.216N. \tag{10}$$

*Since*

$$\frac{256 + \sin 4}{8(\sqrt{17} - 1)} - \frac{1 + 2\sin\frac{1}{2}}{8\left(-1 + \sqrt{\frac{21}{20}}\right)} > 0,$$

*by combining (9) and (10), condition* (*ii*) *of Corollary 1 holds.*

**Figure 1.** The images of *f*(*u*) and *F*(*u*) in Example 1.

*We can further verify condition* (*iii*) *of Corollary 1, since*

$$\limsup\_{\substack{\xi \to +\infty \\ \xi \to +\infty}} \frac{F(k,\xi)}{\xi} = \limsup\_{\substack{\xi \to +\infty \\ \xi \to +\infty}} \frac{\frac{27\pi^3}{8} - \frac{1}{4} + \frac{27\pi^2}{4}\sin(\xi + \frac{\pi}{2})}{\xi} = 0 < \frac{2F\_{\mathbb{C}}}{N(\sqrt{4c^2 + N + 1} - \sqrt{N + 1})},$$

*To sum up, all the conditions of Corollary 1 are satisfied. Hence, for*

$$
\lambda \in \Lambda = \left( \frac{8(\sqrt{17}-1)}{N(256+\sin 4)}, \frac{8\left(-1+\sqrt{1+\frac{1}{N+1}}\right)}{N(1+2\sin\frac{1}{2})} \right) \sqrt{2}
$$

*the boundary value problem admits at least three positive solutions.*

**Example 2.** *Let N* = 3 *and β* = 10*, consider the boundary value (Problem 1) with*

$$f(k,u) = f(u) = \begin{cases} \frac{1}{6} \left[ (3 - 2u) \cos \frac{u}{2} + e^{u-2} \right], & \text{if } u \le 2\pi, \\\frac{1}{6} \left( 4\pi - 3 + e^{2\pi - 2} \right), & \text{if } u > 2\pi, \end{cases}$$

*for each k* ∈ Z(1, *N*). *Then, we have*

$$F(k,\xi) = F(\xi) = \begin{cases} \frac{1}{3}\left[\left(3 - 2\xi\right)\sin\frac{\xi}{2} + 4\left(1 - \cos\frac{\xi}{2}\right)\right] + \frac{\varepsilon^{-2}}{6}\left(e^{\xi} - 1\right), & \text{if } \xi \le 2\pi, \\\frac{1}{6}\left(4\pi - 3 + e^{2\pi - 2}\right)\xi - \frac{\pi}{3}\left(4\pi - 3 + e^{2\pi - 2}\right) + \frac{\varepsilon^{-2}}{6}\left(e^{2\pi - 1}\right) + \frac{8}{3}, & \text{if } \xi > 2\pi. \end{cases}$$

*Moreover, the images of f*(*u*) *and F*(*u*) *are shown in Figure 2. Letting c* = √ *N* + 1 = 2 *and d* = 1 2 *, we have*

$$12\sqrt{1+d^2} = \sqrt{5} < 1 + \sqrt{5} = 1 + \sqrt{1 + \frac{4c^2}{N+1}}.$$

*Now, we prove f*(*u*) > 0 *for each u* ∈ [−2, 2]*. Let*

$$f\_1(\mu) = (3 - 2\mu)\cos\frac{\mu}{2} \quad \text{and} \quad f\_2(\mu) = e^{\mu - 2}.$$

**Figure 2.** The images of *f*(*u*) and *F*(*u*) in Example 2.

*Then f*(*u*) = 1 6 (*f*1(*u*) + *f*2(*u*)) *for u* ≤ 2*π. Besides, it is obvious that f*2(*u*) > 0 *on* [−2, 2] *and f*1(*u*) ≥ 0 *on* - −2, 3 2 *. When u* ∈ - 3 2 , 2 *, we have*

$$f\_1'(u) = \frac{1}{2} \left[ (2u - 3) \sin \frac{u}{2} - 4 \cos \frac{u}{2} \right] \le \frac{1}{2} (\sin 1 - 4 \cos 1) < 0.$$

*Thus f*1(*u*) *is decreasing on* [ 3 2 , 2]*. Since f*2(*u*) *is an increasing function and f*1(2) + *f*2( 3 2 ) = − cos 1 + *e* − 1 2 > 0*, we can get f*(*u*) > 0 *on* [−2, 2]*. In addition, we have*

$$F\_d = \sum\_{k=1}^{N} F(k, d) = N \int\_0^d \frac{1}{6} \left[ (3 - 2u) \cos \frac{u}{2} + e^{u - 2} \right] du \approx 0.6634$$

*and*

$$F\_{\mathcal{C}} = \sum\_{k=1}^{N} F(k, \mathcal{c}) = N \int\_{0}^{\mathcal{C}} \frac{1}{6} \left[ (3 - 2u) \cos \frac{u}{2} + e^{u - 2} \right] du \approx 1.430...$$

*It follows that*

$$\frac{2(\sqrt{1+d^2}-1)}{F\_d} = \frac{\sqrt{5}-2}{F\_d} \approx 0.356\sqrt{5}$$

*and*

$$\frac{-1 + \sqrt{1 + \frac{4c^2}{N+1}}}{F\_c} = \frac{\sqrt{5} - 1}{F\_c} \approx 0.865.$$

*The above together with*

$$\frac{2\sqrt{N(N+1)}\sin\frac{N\pi}{2(N+1)}}{\beta} = \frac{4\sqrt{3}\sin\frac{3\pi}{8}}{10} \approx 0.640$$

*confirms condition* (*T*2) *of Theorem 2.*

*Besides,*

$$\begin{split} \liminf\_{|\tilde{\xi}| \to +\infty} \frac{F(k, \tilde{\xi})}{\tilde{\xi}} &= \liminf\_{\tilde{\xi} \to +\infty} \frac{\frac{1}{6} \left( 4\pi - 3 + e^{2\pi - 2} \right) \xi - \frac{\pi}{3} (4\pi - 3 + e^{2\pi - 2}) + \frac{e^{-2}}{6} (e^{2\pi - 1}) + \frac{8}{3} \xi} \\ &= \frac{1}{6} \left( 4\pi - 3 + e^{2\pi - 2} \right) \\ &\approx 13.673 > \beta \end{split}$$

*for each k* ∈ Z(1, *N*)*. We can see from Theorem 2 that the boundary value problem has at least two non-zero solutions for each λ* ∈ 4 √ 3 sin 3*π* 8 10 , √ 5−1 3 R 2 0 1 6 [(3−2*u*) cos *u* 2 +*e u*−2 ]*du* ≈ (0.640, 0.865)*.*

#### **5. Conclusions**

A discrete Dirichlet boundary value problem involving the mean curvature operator is studied in this paper. Unlike the existing result in [33], we obtained different sufficient conditions of the existence of multiple solutions without assuming that the nonlinear term oscillates at infinity, as shown in Theorems 1 and 2. First, according to the research results of Bonanno in [22], we obtain at least three non-trivial solutions in Theorem 1. In addition, as a supplement to Theorem 1, we prove the existence of at least three positive solutions through the maximum principle. Note that inequality (3) plays an important role in the proof of Theorem 1. For the situation that inequality (3) is not satisfied, under another suitable assumption on the nonlinear term, we still can obtain the existence of at least two non-trivial solutions based on Theorem 2.1 in [23]. It seems that the method used in this paper can be adapted to discuss the existence of homoclinic solutions or periodic solutions of difference equations with *φc*-Laplacian. This will be left as our future work.

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

**Funding:** This work is supported by the National Natural Science Foundation of China (Grant No. 11971126) and the Program for Changjiang Scholars and Innovative Research Team in University (Grant No. IRT 16R16).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

## **On the Absolute Stable Difference Scheme for Third Order Delay Partial Differential Equations**

#### **Allaberen Ashyralyev 1,2,3,**∗ **, Evren Hincal 1 and Suleiman Ibrahim 1**

Received: 16 April 2020; Accepted: 15 June 2020; Published: 19 June 2020

**Abstract:** The initial value problem for the third order delay differential equation in a Hilbert space with an unbounded operator is investigated. The absolute stable three-step difference scheme of a first order of accuracy is constructed and analyzed. This difference scheme is built on the Taylor's decomposition method on three and two points. The theorem on the stability of the presented difference scheme is proven. In practice, stability estimates for the solutions of three-step difference schemes for different types of delay partial differential equations are obtained. Finally, in order to ensure the coincidence between experimental and theoretical results and to clarify how efficient the proposed scheme is, some numerical experiments are tested.

**Keywords:** time delay; third order differential equations; difference scheme; stability

#### **1. Introduction**

Various problems in elasticity theory such as the problems of the longitudinal oscillations of a non-uniform viscoelastic rod, the problem of the longitudinal impact of a perfectly rigid body on a non-uniform finite-length viscoelastic rod with a variable cross-section, problems of wave propagation in a visco-elastic body, etc., lead to third order differential equations without the time delay term ([1–3]). Over the years, nonlocal and local boundary value problems have been of great interest due to their importance in the fields of engineering and science, especially in applied mathematics. Such problems have formed various research fields. Several nonlocal and local boundary value problems for differential equations have been investigated extensively in various works (for example, see [4–12] and the references given therein).

Differential equations having a delay term are used to model sociological, biological, as well as physical processes. They are used to model naturally occurring oscillation systems. A typical example of the occurrence of time delay can be seen in a sampled data control in control theory (see, for example, [13–17]). The presence of delay term in differential equations usually leads to difficulties in analyzing the differential equation. The boundedness and stability and the oscillation property of solutions for a third order delay ordinary differential and difference problems were widely studied (for example, see [18–25] and the references given therein).

Delay partial differential equations (DPDEs) arise in many applications such as control theory, climate models, medicine, biology, and much more (for example, see [26] and the references therein). The independent variables of partial differential equations having delay terms are time *t* together with one or more dimensional variable *x*, representing the position in space. It can also stand for the size of cells, relative DNA content, their level of mutation, as well as other parameters. The

solutions of partial differential equations having delay terms may stand for voltage, temperature, or densities or concentrations of various particles, for instance chemicals, cells, animals, bacteria, and so on. Numerical methods for partial differential equations with delay terms usually lead to specific difficulties, which are usually not present in equations without delay terms. The theory and applications of parabolic and hyperbolic partial differential equations having a time delay term were studied by numerous authors (for example, see [13,27–34] and the references given therein). Recent publications on third order DPDEs are not many.

Several physical models lead to initial-boundary value problems for third order DPDEs (see, e.g., [1,3,8]). It is known that such types of problems can be replaced with the initial value problem for a third order delay differential equation:

$$\begin{cases} \left. u\_{ttt}(t) + Au\_l(t) = bAu(t - w) + f(t), \, 0 < t < \infty, \\\left. u(t) = g(t), \, -w \le t \le 0 \end{cases} \right. \tag{1}$$

in a Hilbert space *H* with unbounded operator *A*. Here, *b* ∈ *R* 1 . Assume that *f*(*t*) is a continuous function on [0, ∞) and *f*(*t*) ∈ *D*(*A* 1/2), *g*(*t*) is a twice continuously differentiable function on [−*w*, 0] and *g* (*k*) (*t*) ∈ *D*(*A* (3−*k*)/2) for *k* = 0, 1, 2.

Let us give the main theorem of paper [35].

**Theorem 1.** *The solution of Problem (1) satisfies the stability estimates:*

$$a\_1 \le da\_0 + \int\_0^w \left\| f(s) \right\|\_{D(A^{\frac{1}{2}})} ds,\\ d = 2 + |b| \, w, \tag{2}$$

$$a\_{n+1} \le da\_n + \sum\_{j=1}^{n+1} \int\_{(j-1)\omega}^{jw} \left\Vert f(s) \right\Vert\_{D(A^{\frac{1}{2}})} ds, n = 1, 2, \cdots,\tag{3}$$

$$a\_0 = \max\left\{ \max\_{-w \le t \le 0} \left\| g\_{tt}(t) \right\|\_{D(A^{\frac{1}{2}})'} \max\_{-w \le t \le 0} \left\| g\_t(t) \right\|\_{D(A)'} \max\_{-w \le t \le 0} \left\| g(t) \right\|\_{D(A^{\frac{2}{2}})} \right\},$$

$$a\_{\mathcal{U}} = \max\left\{ \max\_{(n-1)w \le t \le nw} \left\| u\_{\mathcal{U}}(t) \right\|\_{D(A^{\frac{1}{2}})'} \max\_{(n-1)w \le t \le nw} \left\| u\_t(t) \right\|\_{D(A)'} \max\_{(n-1)w \le t \le nw} \left\| u(t) \right\|\_{D(A^{\frac{2}{2}})} \right\}.$$

In practice, stability estimates for the solution of several problems for third order DPDEs were obtained.

Moreover, publications on the theory and applications of difference schemes (DSs) for third order DPDEs are not available. Thus, the construction and investigation of stable DSs for the approximate solutions of third order DPDEs is of great importance. Our aim in this paper is to construct the absolute stable three-step DS of the first order of accuracy of the third order DPDE for the approximate solution of the problem (1). We consider the uniform set of grid points:

$$(-w, \infty)\_{\tau} = \{ t\_k : t\_k = k\tau, -N \le k < \infty, N\tau = \omega \}$$

with step *τ* > 0. Applying Taylor's decomposition method on three and two points (see [36,37]), we present the DS of the first order of accuracy:

$$\begin{cases} \frac{u\_{k+2} - 3u\_{k+1} + 3u\_{k} - u\_{k-1}}{\tau^3} + A \frac{u\_{k+2} - u\_{k+1}}{\tau} = bAu\_{k-N} + f(t\_k), \; k \ge 1, \\\ u\_{k} = g(t\_k), \; -N + 1 \le k \le 0, \\\ (I + \tau^2 A) \frac{u\_{1} - u\_{0}}{\tau} = g'(0), (I + \tau^2 A) \frac{u\_{2} - 2u\_{1} + u\_{0}}{\tau^2} = g''(0), \\\ (I + \tau^2 A) \frac{u\_{mN+1} - u\_{mN}}{\tau} = \frac{u\_{mN} - u\_{mN-1}}{\tau}, \\\ (I + \tau^2 A) \frac{u\_{mN+2} - 2u\_{mN+1} + u\_{mN}}{\tau^2} = \frac{u\_{mN} - 2u\_{mN-1} + u\_{mN-2}}{\tau^2}, m = 1, 2, \dots \end{cases} (4)$$

for the approximate solution of Problem (1).

The organization of this paper is as follows. In Section 2, the main theorem on the stability of DS (4) is established. In Section 3, stability estimates of DSs for the approximate solution of three problems for third order DPDEs are obtained. Numerical results are provided for one= and two=dimensional third order DPDEs in Section 4. Finally, Section 5 gives the conclusion and our future plans.

#### **2. Stability of DS**

All over the present paper, assume that *H* is a Hilbert space and *A* is a self-adjoint positive definite operator *A* ≥ *δI* in *H* and *R* = (*I* − *iτA* 1 2 ) −1 , *R*e = (*I* + *iτA* 1 2 ) −1 .

Note that three-step DS (4) can obviously be rewritten as the system of single-step and two-step delay DS:

$$\begin{cases} \frac{\underline{u}\_{k+1} - \underline{u}\_{k}}{\tau} = v\_{k}, k \ge 0, \\\\ \frac{v\_{k+1} - 2\underline{v}\_{k} + v\_{k-1}}{\tau^{2}} + Av\_{k+1} = p\_{k}, p\_{k} = bu\_{k-N} + f(t\_{k}), \; k \ge 1, \\\\ u\_{k} = g(t\_{k}), \; -N \le k \le 0, \\\\ (I + \tau^{2}A)v\_{0} = g'(0), (I + \tau^{2}A)\frac{v\_{1} - v\_{0}}{\tau} = g''(0), \\\\ (I + \tau^{2}A)v\_{mN} = \frac{u\_{mN} - u\_{mN-1}}{\tau}, \\\\ (I + \tau^{2}A)\frac{v\_{mN+1} - v\_{mN}}{\tau} = \frac{u\_{mN} - 2u\_{mN-1} + u\_{mN-2}}{\tau^{2}}, m = 1, 2, \dots \end{cases} (5)$$

for the solution of DS (4). Applying DC(5), we can obtain the formula for the solution of DS (4). For this, we will consider two cases 1 ≤ *k* ≤ *N* and *mN* + 1 ≤ *k* ≤ (*m* + 1)*N*, *m* = 1, 2, · · · , separately.

Let 1 ≤ *k* ≤ *N*. Applying (5), we get the following DS:

$$\begin{cases} \frac{\underline{u\_{k}} - \underline{u\_{k-1}}}{\tau} = v\_{k-1}, 1 \le k \le N, \\\frac{\underline{v\_{k+1}} - 2\underline{v\_{k}} + \underline{v\_{k-1}}}{\tau^2} + Av\_{k+1} = p\_k, p\_k = bu\_{k-N} + f(t\_k), 1 \le k \le N - 1, \\\ u\_k = g(t\_k), \ -N \le k \le 0, \\\ (I + \tau^2 A)v\_0 = g'(0), (I + \tau^2 A)\frac{v\_1 - v\_0}{\tau} = g''(0). \end{cases} \tag{6}$$

Therefore, we have that (see [38]):

$$
\mu\_0 = \mathcal{g}(0), \mu\_k = \mathcal{g}(0) + \sum\_{j=0}^{k-1} \tau v\_{j\prime} \ 1 \le k \le N\prime \tag{7}
$$

$$\begin{split} v\_{0} &= R\widetilde{R}g'(0), \ v\_{1} = R\widetilde{R}g'(0) + \tau R\widetilde{R}g''(0), \\ v\_{k} &= \frac{1}{2}[R^{k-1} + \widetilde{R}^{k-1}]R\widetilde{R}g'(0) + \frac{1}{2!}A^{-\frac{1}{2}}R(R^{k} - \widetilde{R}^{k})g''(0) \\ &\quad - \sum\_{s=1}^{k-1} \frac{\tau}{2}A^{-\frac{1}{2}}[R^{k-s} - \widetilde{R}^{k-s}]p\_{s} \\ &= \frac{1}{2}[R^{k-1} + \widetilde{R}^{k-1}]R\widetilde{R}g'(0) + \frac{1}{2!}A^{-\frac{1}{2}}R(R^{k} - \widetilde{R}^{k})g''(0) \\ &+ A^{-1}\left\{\frac{1}{2}\sum\_{s=2}^{k-1}[R^{k-s} + \widetilde{R}^{k-s}](p\_{s-1} - p\_{s}) + 2p\_{k-1} - [R^{k-1} + \widetilde{R}^{k-1}]p\_{1}\right\}, 2 \le k \le N, \\ p\_{k} &= bA g(t\_{k-N}) + f(t\_{k}), \ 1 \le k \le N - 1. \end{split} \tag{8}$$

Applying Formulas (6)–(8), we obtain:

$$u\_{k} = \begin{cases} g(0) + \tau R \tilde{R} g'(0), k = 0, \\ g(0) + \tau R \tilde{R} g'(0), k = 1, \\ g(0) + 2\tau R \tilde{R} g'(0) + \tau^{2} R \tilde{R} g''(0), k = 2, \\ g(0) + \frac{1}{2i} A^{-\frac{1}{2}} \left( R^{k-2} - \tilde{R}^{k-2} \right) R \tilde{R} g'(0) - \frac{1}{2} A^{-1} R \left( 2R \tilde{R} - \left( R^{k-1} + \tilde{R}^{k-1} \right) \right) g''(0), \\ \quad + \sum\_{j=2}^{k-1} \tau \sum\_{s=1}^{j} \frac{1}{2i} A^{-\frac{1}{2}} [R^{j-s} - \tilde{R}^{j-s}] p\_{s\prime}, 3 \le k \le N. \end{cases}$$

By an interchange of the order of summation, we get:

$$u\_{k} = \begin{cases} g(0) + \tau R \tilde{\mathcal{R}} g'(0), k = 0, \\ g(0) + \tau R \tilde{\mathcal{R}} g'(0), k = 1, \\ g(0) + 2\tau R \tilde{\mathcal{R}} g'(0) + \tau^2 R \tilde{\mathcal{R}} g''(0), k = 2, \\ g(0) + \frac{1}{2!} A^{-\frac{1}{2}} \left( R^{k-2} - \tilde{R}^{k-2} \right) R \tilde{\mathcal{R}} g'(0) - \frac{1}{2} A^{-1} R \left( 2R \tilde{\mathcal{R}} - \left( R^{k-1} + \tilde{R}^{k-1} \right) \right) g''(0) \\ + A^{-1} \sum\_{s=1}^{k-2} \frac{\tau}{2} [2I - \left( R^{k-1-s} + \tilde{R}^{k-1-s} \right)] p\_{s\prime}, 3 \le k \le N \end{cases} \tag{9}$$

for the solution of DS (4).

Let 1 + *mN* ≤ *k* ≤ (*m* + 1)*N*, *m* = 1, 2, · · · . Applying (5), we can get the DS:

$$\begin{cases} \frac{u\_k - u\_{k-1}}{\tau} = v\_{k-1} mN + 1 \le k \le (m+1)N, \\\frac{v\_{k+1} - 2v\_k + v\_{k-1}}{\tau^2} + Av\_{k+1} = p\_k, p\_k = bu\_{k-N} + f(t\_k), \\\ mN + 1 \le k \le (m+1)N - 1, \\\ u\_{mN} \text{ is given}, (I + \tau^2 A)v\_{mN} = \frac{u\_{mN} - u\_{mN-1}}{\tau}, \\\ (I + \tau^2 A)\frac{v\_{mN+1} - v\_{mN}}{\tau} = \frac{u\_{mN} - 2u\_{mN-1} + u\_{mN-2}}{\tau^2}. \end{cases} \tag{10}$$

Therefore, we have that (see [38]):

$$u\_k = u\_{mN} + \sum\_{j=mN+2}^{k-1} \tau v\_{j\prime} mN + 1 \le k \le (m+1)N,\tag{11}$$

$$\begin{split} v\_{mN} &= R\widetilde{\mathcal{R}}\frac{u\_{mN} - u\_{mN-1}}{\tau}, v\_{mN+1} = R\widetilde{\mathcal{R}}\frac{u\_{mN} - u\_{mN-1}}{\tau} + \tau R\widetilde{\mathcal{R}}\frac{u\_{mN} - 2u\_{mN-1} + u\_{mN-2}}{\tau^2} \\ &\quad v\_{k} = \frac{1}{2}[R^{k-mN-1} + \widetilde{\mathcal{R}}^{k-mN-1}]R\widetilde{\mathcal{R}}\frac{u\_{mN} - u\_{mN-1}}{\tau} \\ &\quad + \frac{1}{2}A^{-\frac{1}{2}}R(R^{k-mN} - \widetilde{\mathcal{R}}^{k-mN})\frac{u\_{mN} - 2u\_{mN-1} + u\_{mN-2}}{\tau^2} \\ &\quad - \sum\_{s=mN+1}^{k-1} \frac{\tau}{2}A^{-\frac{1}{2}}[R^{k-s} - \widetilde{\mathcal{R}}^{k-s}]p\_{s} \\ &= \frac{1}{2}[R^{k-1} + \widetilde{R}^{k-1}]R\widetilde{\mathcal{R}}\frac{u\_{mN} - u\_{mN-1}}{\tau} + \frac{1}{2\tau}A^{-\frac{1}{2}}R(R^{k} - \widetilde{R}^{k})g^{\prime\prime}(0) \\ &\quad + A^{-1}\left\{\frac{1}{2}\sum\_{s=mN+2}^{k-1}[R^{k-s} + \widetilde{R}^{k-s}](p\_{s-1} - p\_{s}) \\ &\quad + 2p\_{k-mN-1} - [R^{k-mN-1} + \widetilde{R}^{k-mN-1}]p\_{mN+1}\right\}, mN+2 \le k \le (m+1)N. \end{split} \tag{12}$$

Applying Formulas (10)–(12), we can obtain:

$$u\_{k} = \begin{cases} u\_{mN} + \tau R \tilde{\mathcal{R}} \frac{u\_{mN} - u\_{mN-1}}{\tau}, k = mN + 1, \\ u\_{mN} + 2\tau R \tilde{\mathcal{R}} \frac{u\_{mN} - u\_{mN-1}}{\tau} + \tau^{2} R \tilde{\mathcal{R}} \frac{u\_{mN} - 2u\_{mN-1} + u\_{mN-2}}{\tau^{2}}, k = mN + 2, \\ u\_{mN} + \frac{1}{2i} A^{-\frac{1}{2}} \left( R^{k-mN-2} - \widetilde{\mathcal{R}}^{k-mN-2} \right) R \tilde{\mathcal{R}} \frac{u\_{mN} - u\_{mN-1}}{\tau} \\ \quad - \frac{1}{2} A^{-1} R \left( 2R \widetilde{\mathcal{R}} - \left( R^{k-mN-1} + \widetilde{\mathcal{R}}^{k-mN-1} \right) \right) \frac{u\_{mN} - 2u\_{mN-1} + u\_{mN-2}}{\tau^{2}} \\ \quad + \sum\_{j=mN+2}^{k-1} \tau \sum\_{s=mN+1}^{\top} \frac{\tau}{2i} A^{-\frac{1}{2}} [R^{j-s} - \widetilde{\mathcal{R}}^{j-s}] p\_{s\prime}, mN + 3 \le k \le (m+1)N. \end{cases}$$

By an interchange of the order of summation, we get the solution of DS (4):

$$u\_{k} = \begin{cases} u\_{mN} + \tau R \tilde{\mathcal{R}} \frac{u\_{mN} - u\_{mN-1}}{\tau}, k = mN + 1, \\ u\_{mN} + 2\tau R \tilde{\mathcal{R}} \frac{u\_{mN} - u\_{mN-1}}{\tau} + \tau^{2} R \tilde{\mathcal{R}} \frac{u\_{mN} - 2u\_{mN-1} + u\_{mN-2}}{\tau^{2}}, k = mN + 2, \\ u\_{mN} + \frac{1}{2i} A^{-\frac{1}{2}} \left( R^{k-mN-2} - \tilde{\mathcal{R}}^{k-mN-2} \right) R \tilde{\mathcal{R}} \frac{u\_{mN} - u\_{mN-1}}{\tau} \\ \quad - \frac{1}{2} A^{-1} R \left( 2R \tilde{\mathcal{R}} - \left( R^{k-mN-1} + \tilde{\mathcal{R}}^{k-mN-1} \right) \right) \frac{u\_{mN} - 2u\_{mN-1} + u\_{mN-2}}{\tau^{2}} \\ \quad + A^{-1} \sum\_{s=mN+1}^{k-2} \frac{\tau}{2} [2I - \left( R^{k-1-s} + \tilde{R}^{k-1-s} \right)] p\_{s\prime}, mN + 3 \le k \le (m+1)N. \end{cases} \tag{13}$$

The following lemma will be needed in the sequel.

**Lemma 1.** *The following estimates are fulfilled:*

$$\|\|R\|\|\_{H\to H\prime} \parallel \widetilde{R} \|\|\_{H\to H} \leq 1\,\tag{14}$$

$$\|\|\boldsymbol{R}\widetilde{\boldsymbol{R}}^{-1}\|\|\_{\boldsymbol{H}\to\boldsymbol{H}\prime} \; \|\|\widetilde{\boldsymbol{R}}\boldsymbol{R}^{-1}\|\|\_{\boldsymbol{H}\to\boldsymbol{H}} \leq 1,\tag{15}$$

$$\|\operatorname{\tau}A^{\frac{1}{2}}\mathcal{R}\|\_{H\to H'}\|\operatorname{\tau}A^{\frac{1}{2}}\tilde{\mathcal{R}}\|\_{H\to H}\le 1.\tag{16}$$

The proof of the estimates (14)–(16) is based on the spectral theory of a self-adjoint operator in a Hilbert space [39].

Now, let us study the stability of DS (4).

**Theorem 2.** *The solution of DS (4) satisfies the following stability estimates:*

$$b\_1 \le (2 + \tau |b|(N-2))b\_0 + \tau \sum\_{s=1}^{N-2} \|A^{\frac{1}{2}} f(t\_s)\|\_{H\nu} \tag{17}$$

$$b\_{m+1} \le (2+\tau|b|(N-2))b\_m + \tau \sum\_{s=mN+1}^{(m+1)N} \|A^{\frac{1}{2}}f(t\_s)\|\_{H^s} \tag{18}$$

$$b\_0 = \max\left\{ \max\_{-N \le k \le 0} \|A^{\frac{1}{2}} g\_{lt}(t\_k)\|\_{H^\prime} \max\_{-N \le k \le 0} \|Ag\_{lt}(t\_k)\|\_{H^\prime} \max\_{-N \le k \le 0} \|A^{\frac{2}{2}} g(t\_k)\|\_{H} \right\},$$

$$b\_m = \max\left\{ \max\_{(m-1)N \le k \le mN-2} \|A^{\frac{1}{2}} \frac{u\_{k+2} - 2u\_{k+1} + u\_k}{\tau^2} \|\_{H^\prime} \right\}$$

$$\max\_{(m-1)N+1 \le k \le mN} \|A \frac{u\_k - u\_{k-1}}{\tau} \|\_{H^\prime} \frac{1}{2} \max\_{(m-1)N \le k \le mN} \|A^{\frac{2}{2}} u\_k \|\_{H} \right\}, m = 1, 2, \dots, \dots$$

**Proof.** Let us estimate *b*1.

From that and *u*0 = *g*(0), it follows that:

$$\begin{aligned} \frac{1}{2} \max\_{0 \le k \le N} \|A^{\frac{3}{2}} u\_k\|\_H \\\\ \le (2 + \tau |b|(N-2)) b\_0 + \tau \sum\_{s=1}^{N-2} \|A^{\frac{1}{2}} f(t\_s)\|\_H \end{aligned} \tag{19}$$

for the solution of DS (4). Applying Formulas (6)–(8), we can write:

$$\frac{u\_1 - u\_0}{\tau} = v\_0 = R\widetilde{R}g\_t(0),$$

$$\frac{u\_2 - u\_1}{\tau} = v\_1 = R\widetilde{R}g\_t(0) + \tau R\widetilde{R}g\_{tt}(0),$$

$$\frac{u\_k - u\_{k-1}}{\tau} = v\_{k-1} = \frac{1}{2}[R^{k-2} + \widetilde{R}^{k-2}]R\widetilde{R}g\_t(0) + \frac{1}{2i}R(R^{k-1} - \widetilde{R}^{k-1})A^{-\frac{1}{2}}g\_{tt}(0)$$

$$-\frac{\tau}{2i}\sum\_{s=1}^{k-2}[R^{k-1-s} - \widetilde{R}^{k-1-s}]bA^{\frac{1}{2}}g(t\_{s-N}) - \frac{\tau}{2i}\sum\_{s=1}^{k-2}[R^{k-1-s} - \widetilde{R}^{k-1-s}]A^{-\frac{1}{2}}f(t\_s), \, 3 \le k \le N.$$

Using this formula and the estimates (14)–(16), we obtain:

$$\|A\frac{u\_1 - u\_0}{\tau}\|\_{H} \le \|R\tilde{R}\|\|\_{H \to H} \|g\_t(0)\|\_{H} \le b\_0.$$

$$\|A\frac{u\_2 - u\_1}{\tau}\|\_{H} \le \|R\tilde{R}\|\_{H \to H} \|Ag\_t(0)\|\_{H} + \|\tau A^{\frac{1}{2}}\mathbb{R}\|\_{H \to H} \|\tilde{R}\|\_{H \to H} \|A^{\frac{1}{2}}g\_t(0)\|\_{H} \le 2b\_0.$$

$$\|A\frac{u\_k - u\_{k-1}}{\tau}\|\_{H} \le \frac{1}{2} \|\|R\|\_{H \to H}^{k-2} + \|\tilde{R}\|\_{H \to H}^{k-2} \|\|R\tilde{R}\|\_{H \to H}\|\|g\_t(0)\|\_{H}$$

$$+ \frac{1}{2} [\|R\|\_{H \to H}^{k-1} + \|\tilde{R}\|\_{H \to H}^{k-1}] \|R\|\_{H \to H} \|A^{\frac{1}{2}}g\_t(0)\|\_{H} + \frac{\tau}{2} \|b\|\sum\_{s=1}^{N-2} (\|R\|\_{H \to H}^{k-1-s} + \|\tilde{R}\|\_{H \to H}^{k-1-s}) \|A^{\frac{3}{2}}g(s\_{s-N})\|\_{H}$$

$$+ \frac{\tau}{2} \sum\_{s=1}^{N-2} [\|R\|\_{H \to \tilde{H}}^{k-1-s} + \|\tilde{R}\|\_{H \to H}^{k-1-s}] \|A^{\frac{1}{2}}f(s\_s)\|\_{H}$$

*Symmetry* **2020**, *12*, 1033

$$\leq (2+\tau|b|(N-2))b\_0 + \tau \sum\_{s=1}^{N-2} ||A^{\frac{1}{2}}f(t\_s)||\_{H^s} \\ \mathcal{S} \leq k \leq N.$$

Combining these estimates, we obtain:

$$\max\_{1 \le k \le N} \left\| A \frac{u\_k - u\_{k-1}}{\tau} \right\|\_H \le (2 + \tau |b|(N - 2))b\_0 + \tau \sum\_{s=1}^{N-2} \|A^{\frac{1}{2}} f(t\_s)\|\_H \tag{20}$$

for the solution of DS (4). Applying Formulas (6)–(8), we can write:

$$\frac{\mu\_2 - 2\mu\_1 + \mu\_0}{\tau^2} = \tau^{-1} \left(v\_1 - v\_0 \right) = \tilde{R}R\underline{g}\_1(0),$$

$$\frac{\mu\_3 - 2\mu\_2 + \mu\_1}{\tau^2} = \tau^{-1} \left(v\_2 - v\_1\right)$$

$$= -\tau A(R\bar{R})^2 g\_{tt}(0) + \left(I - \tau^2 A\right) \left(\tilde{R}\right)^2 g\_{tt}(0) + \tau R\tilde{R}p\_t$$

$$= -\tau A(R\bar{R})^2 g\_{tt}(0) + \left(I - \tau^2 A\right) \left(\tilde{R}R\right)^2 g\_{tt}(0) + \tau R\tilde{R}[bA g(t\_{k-N}) + f(t\_1)],$$

$$\frac{u\_{k+2} - 2u\_{k+1} + u\_k}{\tau^2} = \tau^{-1} \left(v\_{k+1} - v\_k\right)$$

$$= \tau^{-1} \left\{\frac{1}{2} [R^k + \tilde{R}^k] R\tilde{R}g\_t(0) + \frac{1}{2} A^{-\frac{1}{2}} [R^{k+1} - \tilde{R}^{k+1}] Rg\_{tt}(0)\right\}$$

$$- \sum\_{k=1}^k \frac{\tau}{2} A^{-\frac{1}{2}} [R^{k+1-s} - \tilde{R}^{k+1-s}] p\_s = \frac{1}{2} [R^{k+1} + \tilde{R}^{k+1}] R\tilde{R}g\_t(0)$$

$$- \frac{1}{2i} A^{-\frac{1}{2}} [R^k - \tilde{R}^k] Rg\_{tt}(0) + \sum\_{s=1}^k \frac{\tau}{2} A^{-\frac{1}{2}} [R^{k-s} - \tilde{R}^{k-s}] p\_s\right\}$$

$$= \tau^{-1} \left\{\frac{1}{2} [R^k + \tilde{R}^k - R$$

Using this formula and the estimates (14)–(16), we obtain:

 *A* 1 2 *u*2 − 2*u*1 + *u*0 *τ* 2 *H* ≤ k*RR*ek*H*→*H*k*A* 1 2 *gtt*(0)k*H* ≤ *b*0, *A* 1 2 *u*3 − 2*u*2 + *u*1 *τ* 2 *H* ≤ k*τA* 1 2 *R*k*H*→*H*k*R*k*H*→*H*k*R*ek 2 *H*→*H*k*Agt*(0)k*H* + k*RR*ek*H*→*H*k*A* 1 2 *gtt*(0)k*H* +*τ*|*b*|k*RR*ek*H*→*H*k*A* 3 2 *g*(*t*1−*N*)k*H* + *τ*k*RR*ek*H*→*H*k*A* 3 2 *f*(*t*1)k*H* ≤ (2 + *τ*|*b*|)*b*0 + *τ*k*A* 1 2 *f*(*t*1)k*H*, *A* 1 2 *uk*+2 − 2*uk*+1 + *uk τ* 2 *H* ≤ 1 2 [k*R*k *k H*→*H* + k*R*ek *k H*→*H*]k*RR*ek*H*→*H*k*Agt*(0)k*H*

$$\begin{split} \frac{1}{2} + \frac{1}{2} \| \| R \| \_{H \to H} \| \| R \| \_{H \to H}^{k+1} + \| \tilde{R} \| \_{H \to H}^{k+1} \| \| A^{\frac{1}{2}} g\_H(0) \| \_{H} + \frac{\tau}{2} \| b \| \sum\_{s=1}^{k-1} \| \| R \| \_{H \to H}^{k-s+1} + \| \tilde{R} \| \_{H \to H}^{k-s+1} \| \| \| A^{\frac{2}{2}} g(t\_{s-N}) \| \_{H} \\ + \frac{\tau}{2} \sum\_{s=1}^{k-1} \| \| R \| \_{H \to H}^{k-s+1} + \| \tilde{R} \| \_{H \to H}^{k-s+1} \| \| A^{\frac{1}{2}} f(t\_{s}) \| \_{H} + \tau \| b \| \| R \tilde{R} \| \_{H \to H} \| \| A^{\frac{2}{2}} g(t\_{k-N}) \| \_{H} \\ + \tau \| R \tilde{R} \| \_{H \to H} \| \| A^{\frac{1}{2}} f(t\_{k}) \| \_{H} \\ \leq (2 + \tau \| b \| (N - 2) ) b\_{0} + \tau \sum\_{s=1}^{N-2} \| \| A^{\frac{1}{2}} f(t\_{s}) \| \_{H}, 2 \leq k \leq N - 2. \end{split}$$

Combining these estimates, we can get:

$$\max\_{0 \le k \le N-2} \left\| A^{\frac{1}{2}} \frac{u\_{k+2} - 2u\_{k+1} + u\_k}{\tau^2} \right\|\_H$$

$$\le (2 + \tau |b|(N-2))b\_0 + \tau \sum\_{s=1}^{N-2} \|A^{\frac{1}{2}} f(t\_s)\|\_H \tag{21}$$

for the solution of DS (4). Estimate (17) follows from (19)–(21).

Now, let us estimate *bm*+1.

Combining these estimates, we obtain:

$$\frac{1}{2} \max\_{mN+1 \le k \le (m+1)N} ||A^{\frac{3}{2}} u\_k||\_H$$

$$\le (2 + \tau |b|(N-2))b\_m + \tau \sum\_{s=mN+1}^{(m+1)N-2} ||A^{\frac{1}{2}} f(t\_s)||\_H \tag{22}$$

for the solution of DS (4). Applying Formulas (10)–(12), we can write:

$$\frac{\mathfrak{u}\_{mN+1} - \mathfrak{u}\_{mN}}{\tau} = v\_{mN} = R\widetilde{R}\frac{\mathfrak{u}\_{mN} - \mathfrak{u}\_{mN-1}}{\tau}.$$

$$\begin{split} \frac{u\_{mN+2} - u\_{mN+1}}{\tau} &= v\_{mN+1} = R\widetilde{\bar{R}} \frac{u\_{mN} - u\_{mN-1}}{\tau} + \tau R\widetilde{\bar{R}} \frac{u\_{mN} - 2u\_{mN-1} + u\_{mN-2}}{\tau^2} \\ &\frac{u\_k - u\_{k-1}}{\tau} = v\_{k-1} = \frac{1}{2} [R^{k-mN-2} + \widetilde{\bar{R}}^{k-mN-2}] R\widetilde{\bar{R}} \frac{u\_{mN} - u\_{mN-1}}{\tau} \\ &+ \frac{1}{2i} R (R^{k-mN-1} - \widetilde{R}^{k-mN-1}) A^{-\frac{1}{2}} \frac{u\_{mN} - 2u\_{mN-1} + u\_{mN-2}}{\tau^2} \\ &- \frac{\tau}{2i} \sum\_{s=mN+1}^{k-2} [R^{k-1-s} - \widetilde{R}^{k-1-s}] b A^{\frac{1}{2}} u (t\_{s-N}) \\ &- \frac{\tau}{2i} \sum\_{s=mN+1}^{k-2} [R^{k-1-s} - \widetilde{R}^{k-1-s}] A^{-\frac{1}{2}} f(t\_{\bar{s}}), \; mN+3 \le k \le (m+1)N. \end{split}$$

Using this formula and the estimates (14)–(16), we can obtain:

k*A umN*+1 − *umN τ* k*H* ≤ k*RR*ek*H*→*H*k*A umN* − *umN*−1 *τ* k*H* ≤ *bm*, k*A umN*+2 − *umN*+1 *τ* k*H* ≤ k*RR*ek*H*→*H*k*A umN* − *umN*−1 *τ* k*H* +k*τA* 1 2 *R*k*H*→*H*k*R*ek*H*→*H*k*A* 1 2 *umN* − 2*umN*−1 + *umN*−2 *τ* 2 k*H* ≤ 2*bm*, k*A uk* − *uk*−1 *τ* k*H* ≤ 1 2 [k*R*k *k*−*mN*−2 *H*→*H* + k*R*ek *k*−*mN*−2 *H*→*H* ]k*RR*ek*H*→*H*k*A umN* − *umN*−1 *τ* k*H* + 1 2 [k*R*k *k*−1 *H*→*H* + k*R*ek *k*−1 *H*→*H* ]k*R*k*H*→*H*k*A* 1 2 *umN* − 2*umN*−1 + *umN*−2 *τ* 2 k*H* + *τ* 2 |*b*| (*m*+1)*N*−2 ∑ *s*=*mN*+1 [k*R*k *k*−1−*s H*→*H* + k*R*ek *k*−1−*s H*→*H* ]k*A* 3 2 *u*(*ts*−*N*)k*H* + *τ* 2 (*m*+1)*N*−2 ∑ *s*=*mN*+1 [k*R*k *k*−1−*s H*→*H* + k*R*ek *k*−1−*s H*→*H* ]k*A* 1 2 *f*(*ts*)k*H* ≤ (2 + *τ*|*b*|(*N* − 2))*bm* + *τ* (*m*+1)*N*−2 ∑ *s*=*mN*+1 k*A* 1 2 *f*(*ts*)k*H*, *mN* + 3 ≤ *k* ≤ (*m* + 1)*N*.

Combining these estimates, we obtain:

$$\max\_{m:N+1 \le k \le (m+1)N} \left\| A \frac{u\_k - u\_{k-1}}{\tau} \right\|\_H \le (2 + \tau |b|(N-2))b\_m + \tau \sum\_{s=mN+1}^{(m+1)N-2} \|A^{\frac{1}{2}} f(t\_s)\|\_H \tag{23}$$

for the solution of DS (4). Applying Formulas (10)–(12), we can write:

$$\frac{u\_{mN+2} - 2u\_{mN+1} + u\_{mN}}{\tau^2} = \tau^{-1} \left(v\_{mN+1} - v\_{mN}\right) = \tilde{R}R\frac{u\_{mN} - 2u\_{mN-1} + u\_{mN-2}}{\tau^2}$$

$$\frac{u\_{mN+3} - 2u\_{mN+2} + u\_{mN+1}}{\tau^2} = \tau^{-1} \left(v\_{mN+2} - v\_{mN+1}\right)$$

$$= -\tau A \left(R\tilde{R}\right)^2 \frac{u\_{mN} - u\_{mN-1}}{\tau} + \left(I - \tau^2 A\right) \left(\tilde{R}R\right)^2 \frac{u\_{mN} - 2u\_{mN-1} + u\_{mN-2}}{\tau^2} + \tau R \tilde{R} p\_{mN+1}$$

$$= -\tau A (R\tilde{R})^2 \frac{u\_{mN} - u\_{mN-1}}{\tau} + \left(I - \tau^2 A\right) \left(\tilde{R}R\right)^2 \frac{u\_{mN} - 2u\_{mN-1} + u\_{mN-2}}{\tau^2}$$

$$\frac{u\_{k+2} - 2u\_{k+1} + u\_k}{\tau^2} = \tau^{-1} \left(v\_{k+1} - v\_k\right)$$

+*τRR*e[*bAu*(*t*1−*N*) + *f*(*tmN*+1)], = *τ* −1 1 2 [*R k*−*mN* + *R*e*k*−*mN*]*RR*e *umN* − *umN*−1 *τ* + 1 2*i A* − 1 2 [*R k*−*Mn*+1 − *R*e*k*−*Mn*+1 ]*R umN* − 2*umN*−1 + *umN*−2 *τ* 2 − *k* ∑ *s*=*Mn*+1 *τ* 2*i A* − 1 2 [*R k*+1−*s* − *R*e*k*+1−*s* ]*ps* − 1 2 [*R k*−*Mn*−1 + *R*e*k*−*Mn*−1 ]*RR*e *umN* − *umN*−1 *τ* − 1 2*i A* − 1 2 [*R k*−*mN* − *R*e*k*−*mN*]*R umN* − 2*umN*−1 + *umN*−2 *τ* 2 + *k* ∑ *s*=*mN*+1 *τ* 2*i A* − 1 2 [*R k*−*s* − *R*e*k*−*s* ]*ps* = *τ* −1 1 2 [*R k*−*mN* + *R*e*k*−*mN* − *R k*−*mN*−1 − *R*e*k*−*mN*−1 ]*RR*e *umN* − *umN*−1 *τ* + 1 2*i A* −1 [*R k*+1 − *R*e*k*+1 − *R k* + *R*e*k* ]*R umN* − 2*umN*−1 + *umN*−2 *τ* 2 − *τ* 2*i A* − 1 2 [*R* − *R*e]*pk*−*mN* + *k*−1 ∑ *s*=*mN*+1 *τ* 2*i A* − 1 2 [*R k*−*s* − *R*e*k*−*s* − *R k*+1−*s* + *R*e*k*+1−*s* ]*ps* = − *i* 2 [*R k*−*mN* − *R*e*k*−*mN*]*RRA*e 1 2 *umN* − *umN*−1 *τ* + 1 2 *R*[*R k*+1 + *R*e*k*+1 ] *umN* − 2*umN*−1 + *umN*−2 *τ* 2 −*τRRAbu* e (*tk*−*N*) − *τRR f* e (*tk*−*mN*) − *τ* 2 *k*−1 ∑ *s*=*mN*+1 [*R k*−*s*+1 + *R*e*k*−*s*+1 ]*Abu*(*ts*−*N*) − *τ* 2 *k*−1 ∑ *s*=*mN*+1 [*R k*−*s*+1 + *R*e*k*−*s*+1 ] *f*(*ts*), *mN* + 2 ≤ *k* ≤ (*m* + 1)*N* − 2.

Using this formula and the estimates (14)–(16), we obtain that:

 *A* 1 2 *umN*+2 − 2*umN*+1 + *umN τ* 2 *H* ≤ k*RR*ek*H*→*H*k*A* 1 2 *umN* − 2*umN*−1 + *umN*−2 *τ* 2 k*H* ≤ *bm*, *A* 1 2 *umN*+3 − 2*umN*+2 + *umN*+1 *τ* 2 *H* ≤ k*τA* 1 2 *R*k*H*→*H*k*R*k*H*→*H*k*R*ek 2 *H*→*H*k*A umN* − *umN*−1 *τ* k*H* +k*RR*ek*H*→*H*k*A* 1 2 *umN* − 2*umN*−1 + *umN*−2 *τ* 2 k*H* +*τ*|*b*|k*RR*ek*H*→*H*k*A* 3 2 *u*(*t*1−*N*)k*H* + *τ*k*RR*ek*H*→*H*k*A* 3 2 *f*(*tmN*+1)k*H* ≤ (2 + *τ*|*b*|)*bm* + *τ*k*A* 1 2 *f*(*tmN*+1)k*H*, *A* 1 2 *uk*+2 − 2*uk*+1 + *uk τ* 2 *H* ≤ 1 2 [k*R*k *k*−*mN H*→*H* + k*R*ek *k*−*mN H*→*H* ]k*RR*ek*H*→*H*k*A umN* − *umN*−1 *τ* k*H*

$$\begin{split} &+\frac{1}{2}||\boldsymbol{R}||\_{H\to H} \| \| \boldsymbol{R} \|\_{H\to H}^{k+1} + \| \tilde{\boldsymbol{R}} \|\_{H\to H}^{k+1} \| \| \boldsymbol{A}^{\frac{1}{2}} \frac{u\_{mN} - 2u\_{mN-1} + u\_{mN-2}}{\tau^2} \| \boldsymbol{h} \\ &+\frac{\tau}{2} |\boldsymbol{b}| \sum\_{s=mN+1}^{k-1} \left( \| \boldsymbol{R} \|\_{H\to H}^{k-s+1} + \| \tilde{\boldsymbol{R}} \|\_{H\to H}^{k-s+1} \| \| \boldsymbol{A}^{\frac{2}{2}} \boldsymbol{u}(t\_{s-N}) \| \right) \|\_{H} \\ &+\frac{\tau}{2} \sum\_{s=mN+1}^{k-1} \left( \| \boldsymbol{R} \|\_{H\to H}^{k-s+1} + \| \tilde{\boldsymbol{R}} \|\_{H\to H}^{k-s+1} \| \| \boldsymbol{A}^{\frac{1}{2}} f(t\_{s}) \| \|\_{H} \\ &+\tau |\boldsymbol{b}| \| \boldsymbol{R} \tilde{\boldsymbol{R}} \|\_{H\to H} \| \boldsymbol{A}^{\frac{1}{2}} \boldsymbol{u}(t\_{k-(m+1)N}) \| \|\_{H} + \tau \| \boldsymbol{R} \tilde{\boldsymbol{R}} \|\_{H\to H} \| \boldsymbol{A}^{\frac{1}{2}} f(t\_{k-mN}) \| \|\_{H} \\ &\leq (2+\tau|\boldsymbol{b}|(N-2)) b\_{m} + \tau \sum\_{s=mN+1}^{(m+1)N-2} \| \boldsymbol{A}^{\frac{1}{2}} f(t\_{s}) \|\_{H'} \, mN + 2 \leq k \leq (m+1)N-2. \end{split}$$

Combining these estimates, we obtain:

$$\max\_{mN \le k \le (m+1)N-2} \left\| A^{\frac{1}{2}} \frac{u\_{k+2} - 2u\_{k+1} + u\_k}{\tau^2} \right\|\_{H}$$

$$\le (2 + \tau |b|(N-2))b\_m + \tau \sum\_{s=mN+1}^{(m+1)N-2} ||A^{\frac{1}{2}} f(t\_s)||\_H \tag{24}$$

for the solution of DS (4). Estimate (18) follows from (22)–(24). Theorem 2 is proven.

Note that applying Theorem 2, we can obtain the stability estimate:

$$\max\_{mN \le k \le (m+1)N - 2} \left\| A^{\frac{1}{2}} \frac{u\_{k+2} - 2u\_{k+1} + u\_k}{\tau^2} \right\|\_{H} + \max\_{mN + 1 \le k \le (m+1)N} \left\| A \frac{u\_k - u\_{k-1}}{\tau} \right\|\_{H} \tag{25}$$

$$+ \frac{1}{2} \max\_{mN + 1 \le k \le (m+1)N} \| A^{\frac{3}{2}} u\_k \|\_{H} \le (2 + \tau |b|(N-2))^m b\_0$$

$$+ \sum\_{j=1}^m (2 + \tau |b|(N-2))^{m-j} \tau \sum\_{s=(j-1)N+1}^{jN} \| A^{\frac{1}{2}} f(t\_s) \|\_{H} m = 0, 1, \dots$$

for the solution of DS (4).

#### **3. Applications**

Note that the generality of this approach permits studying of a general class of DPDEs. We consider the applications of Theorem 2 for three types of problems. First, the mixed problem for the one-dimensional DPDE with nonlocal conditions:

$$\begin{cases} \frac{\partial^3 u(t, \mathbf{x})}{\partial t^3} - (a(\mathbf{x}) u\_{\mathbf{x}\mathbf{x}}(t, \mathbf{x}))\_{\mathbf{x}} + \delta u\_{\mathbf{l}}(t, \mathbf{x}) \\ = b \left( - (a(\mathbf{x}) u\_{\mathbf{x}}(t - w, \mathbf{x}))\_{\mathbf{x}} + \delta u(t - w, \mathbf{x}) \right) + f(t, \mathbf{x})\_{\mathbf{x}} \\ 0 < t < \infty, \ 0 < \mathbf{x} < l, \\ u(t, 0) = u(t, l), \ u\_{\mathbf{x}}(t, 0) = u\_{\mathbf{x}}(t, l), \ 0 \le t < \infty, \\ u(t, \mathbf{x}) = g(t, \mathbf{x}), -w \le t \le 0, 0 \le \mathbf{x} \le l \end{cases} \tag{26}$$

is studied. Under compatibility conditions, Problem (26) has a unique solution *u*(*t*, *x*) for the given smooth functions *a*(*x*) ≥ *a* > 0, *x* ∈ (0, *l*), *δ* > 0, *a*(*l*) = *a*(0), *g*(*t*, *x*), −*w* ≤ *t* ≤ 0, 0 ≤ *x* ≤ *l*, *f*(*t*, *x*), 0 < *t* < ∞, 0 < *x* < *l*, and *b* ∈ *R* 1 .

The construction of full discretization to Problem (26) is completed in two stages. In the first stage, we consider the uniform grid space:

$$[0,l]\_{\mathfrak{h}} = \{ \mathfrak{x} = \mathfrak{x}\_{\mathfrak{n}} : \mathfrak{x}\_{\mathfrak{n}} = \mathfrak{n}h, \ 0 \le \mathfrak{n} \le M, \ \mathrm{Mh} = l \}$$

with step *h* > 0. Let *L*2*h* = *L*2([0, *l*]*h* ) be a Hilbert space of the grid functions *ϕ h* (*x*) = {*ϕn*} *M* 0 defined on [0, *l*]*h* , equipped with the norm:

$$||\varrho^h||\_{L\_{2h}} = \left(\sum\_{\mathbf{x}\in[0,l]\_h} |\varrho(\mathbf{x})|^2 h\right)^{1/2}.$$

Let *A x h* be the second order difference operator defined by:

$$A\_h^x \boldsymbol{\varrho}^h(\mathbf{x}) = \left\{-(a(\mathbf{x})\boldsymbol{\varrho}\_{\overline{\mathbf{x}}})\_{\mathbf{x},\mathcal{U}} + \delta \boldsymbol{\varrho}\_n\right\}\_1^{M-1} \tag{27}$$

acting in the space of grid functions *ϕ h* (*x*) = {*ϕn*} *M* 0 satisfying the conditions *ϕ*0 = *ϕM*, *ϕ*1 − *ϕ*0 = *ϕM* − *ϕM*−1. It is well known that *A x h* is a self-adjoint positive definite operator in *L*2*h* . Applying *A x h* in (26), we can obtain the initial value problem for an infinite system of third order differential equations:

$$\begin{cases} \begin{aligned} u\_{ttt}^h(t, \mathbf{x}) + A\_h^\mathbf{x} u\_t^h(t, \mathbf{x}) &= b A\_h^\mathbf{x} u(t - w, \mathbf{x}) + f^h(t, \mathbf{x}), \\ 0 < t < \infty, &\mathbf{x} \in [0, l]\_h \\ u^h(t, \mathbf{x}) = g^h(t, \mathbf{x}), -w \le t \le 0, \mathbf{x} \in [0, l]\_h. \end{aligned} \end{cases} \tag{28}$$

In the second stage, we use DS (4) for (28):

 *u h k*+2 (*x*)−3*u h k*+1 (*x*)+3*u h k* (*x*)−*u h k*−1 (*x*) *τ* 3 + *A x h u h k*+2 (*x*)−*u h k*+1 (*x*) *τ* = *bAx h u h k*−*N* (*x*) + *f h k* (*x*), *f h k* (*x*) = *f h* (*tk* , *x*), *k* ≥ 1, *x* ∈ [0, *l*]*h* , *u h k* (*x*) = *g h* (*tk* , *x*), −*N* ≤ *k* ≤ 0, (*Ih* + *τ* 2*A x h* ) *u h* 1 (*x*)−*u h* 0 (*x*) *τ* = *g h t* (0, *x*), (*Ih* + *τ* 2*A x h* ) *u h* 2 (*x*)−2*u h* 1 (*x*)+*u h* 0 (*x*) *τ* 2 = *g h tt*(0, *x*), *x* ∈ [0, *l*]*h* , (*Ih* + *τ* 2*A x h* ) *u h mN*+1 (*x*)−*u h mN*(*x*) *τ* = *u h mN*(*x*)−*u h mN*−1 (*x*) *τ* , (*Ih* + *τ* 2*A x h* ) *u h mN*+2 (*x*)−2*u h mN*+1 (*x*)+*u h mN*(*x*) *τ* 2 = *u h mN*(*x*)−2*u h mN*−1 (*x*)+*u h mN*−2 (*x*) *τ* 2 , *m* = 1, 2, ... (29)

**Theorem 3.** *The solutions of DS (29) obey the stability estimates:*

$$\begin{split} \max\_{mN \le k \le (m+1)N-2} \left\| \frac{\boldsymbol{u}\_{k+2}^{h} - 2\boldsymbol{u}\_{k+1}^{h} + \boldsymbol{u}\_{k}^{h}}{\tau^{2}} \right\|\_{W\_{2h}^{2}} + \max\_{mN+1 \le k \le (m+1)N} \left\| \frac{\boldsymbol{u}\_{k}^{h} - \boldsymbol{u}\_{k-1}^{h}}{\tau} \right\|\_{W\_{2h}^{2}} \\ + \frac{1}{2} \max\_{mN+1 \le k \le (m+1)N} \|\boldsymbol{u}\_{k}^{h}\|\_{W\_{2h}^{3}} \le C\_{1} \left[ (2 + \tau |b|(N-2))^{m} \boldsymbol{b}\_{0}^{h} \right] \\ + \sum\_{j=1}^{m} (2 + \tau |b|(N-2))^{m-j} \tau \sum\_{s=(j-1)N+1}^{jN} \|\boldsymbol{A}^{\frac{1}{2}} f(t\_{s})\|\_{W\_{2h}^{1}} \right], m = 0, 1, \dots \\ \boldsymbol{b}\_{0}^{h} = \max \left\{ \max\_{-N \le k \le 0} \|\boldsymbol{A}^{\frac{1}{2}} \boldsymbol{g}\_{th}^{h}(t\_{k})\|\_{W\_{2h}^{1}} \max\_{-N \le k \le 0} \|\boldsymbol{g}\_{t}^{h}(t\_{k})\|\_{W\_{2h}^{2}} \cdot \max\_{-N \le k \le 0} \|\boldsymbol{g}^{h}(t\_{k})\|\_{W\_{2h}^{3}} \right\}. \\ C\_{1} \text{ does not depend on  $\tau$  h  $ph$  and  $f$ } (k) \end{split}$$

*hold, where C*1 *does not depend on τ*, *h*, *g* (*tk* ), *and f h k* (*x*). **Proof.** DS (29) can be written in abstract form:

$$\begin{cases} \frac{u\_{k+2}^h - 3u\_{k+1}^h + 3u\_k^h - u\_k^h}{\tau^3} + A\_h \frac{u\_{k+2}^h - u\_{k+1}^h}{\tau} = bA\_h u\_{k-N}^h + f\_k^h, k \ge 1, \\\\ u\_k^h = g\_{k'}^h - N \le k \le 0, \\\\ (I\_h + \tau^2 A\_h) \frac{u\_1^h - u\_0^h}{\tau} = g\_1^h(0), (I\_h + \tau^2 A\_h) \frac{u\_2^h - 2u\_1^h + u\_0^h}{\tau^2} = g\_{tt}^h(0), \\\\ (I\_h + \tau^2 A\_h) \frac{u\_{mN+1}^h - u\_{mN}^h}{\tau} = \frac{u\_{mN}^h - u\_{mN-1}^h}{\tau}, \\\\ (I\_h + \tau^2 A\_h) \frac{u\_{mN+2}^h - 2u\_{mN+1}^h + u\_{mN}^h}{\tau^2} = \frac{u\_{mN}^h - 2u\_{mN-1}^h + u\_{mN-2}^h}{\tau^2}, m = 1, 2, \dots \end{cases} (30)$$

in a Hilbert space *L*2*h* with a self-adjoint positive definite operator *Ah* = *A x h* . Here, *g h k* = *g h k* (*x*), *f h k* = *f h k* (*x*), and *u h k* = *u h k* (*x*) are known and unknown abstract mesh functions defined on [0, *l*]*h* . Therefore, the estimate of Theorem 3 follows from the estimate (25). Theorem 3 is proven.

Second, let Ω be the unit open cube in the *m*-dimensional Euclidean space.

R*n* (*x* = (*x*1, · · · , *xm*) : 0 < *xk* < 1, *k* = 1, · · · , *n*) with boundary *S*, Ω = Ω ∪ *S*. In [0, ∞) × Ω, the mixed problem for the DPDE with the Dirichlet condition:

$$\begin{cases} u\_{\text{lift}}(t, \mathbf{x}) - \sum\_{r=1}^{n} \left( a\_{r}(\mathbf{x}) u\_{\text{ltx}}(t, \mathbf{x}) \right)\_{\mathbf{x}\_{r}} = -b \sum\_{r=1}^{n} \left( a\_{r}(\mathbf{x}) u\_{\text{ltx}}(t - w, \mathbf{x}) \right)\_{\mathbf{x}\_{r}}, \\ 0 < t < \infty, \mathbf{x} \in \Omega, \\ u(t, \mathbf{x}) = \mathbf{0}, \ \mathbf{x} \in \mathcal{S}, \ \mathbf{0} \le t < \infty, \\ u(t, \mathbf{x}) = \mathbf{g}(t, \mathbf{x}), -w \le t \le \mathbf{0}, \mathbf{x} \in \overline{\Omega} \end{cases} \tag{31}$$

is investigated. Under compatibility conditions, Problem (31) has a unique solution *u*(*t*, *x*) for the given smooth functions *ar*(*x*) ≥ *a* > 0, (*x* ∈ Ω), *g*(*t*, *x*), −*w* ≤ *t* ≤ 0, *x* ∈ Ω, *f*(*t*, *x*), 0 < *t* < ∞, *x* ∈ Ω, and *b* ∈ *R* 1 .

The construction of full discretization to Problem (31) is completed in two stages. In the first stage, we consider the uniform grid space:

$$\begin{aligned} \Omega\_{\mathbb{H}} &= \{ \mathbf{x} = \mathbf{x}\_{r} = (h\_{1}j\_{1}, \dots, h\_{n}j\_{n}), j = (j\_{1}, \dots, j\_{n}), \ 0 \le j\_{r} \le \mathbf{N}\_{r} \\ \mathbf{N}\_{\mathbf{r}}h\_{\mathbf{r}} &= 1, \ \mathbf{r} = 1, \dots, n \}, \ \Omega\_{\mathbb{H}} = \overline{\Omega}\_{\mathbb{H}} \cap \Omega, \mathcal{S}\_{\mathbb{H}} = \overline{\Omega}\_{\mathbb{H}} \cap \mathcal{S} \end{aligned}$$

and introduce the Hilbert space *L*2*h* = *L*2(Ω*h* ) of the grid functions *ϕ h* (*x*) = {*ϕ* (*h*1 *j*1, · · · , *hn jn*)} defined on Ω*h* equipped with the norm:

$$\left\| \left| \boldsymbol{\varrho}^{h} \right| \right\|\_{L\_{2h}} = \left( \sum\_{\boldsymbol{\chi} \in \Omega\_{h}} \left| \boldsymbol{\varrho}^{h} \left( \boldsymbol{\chi} \right) \right|^{2} h\_{1} \cdot \cdots \cdot h\_{n} \right)^{\frac{1}{2}}.$$

We consider the difference operator *A x h* defined by the formula:

$$A\_{\boldsymbol{h}}^{\boldsymbol{x}}\boldsymbol{u}^{\boldsymbol{h}} = -\sum\_{\boldsymbol{r}=1}^{n} \left(a\_{\boldsymbol{r}}(\boldsymbol{x})\boldsymbol{u}\_{\boldsymbol{x}\_{\boldsymbol{r}}}^{\boldsymbol{h}}\right)\_{\boldsymbol{x}\_{\boldsymbol{r}},\boldsymbol{j}\_{\boldsymbol{r}}} \tag{32}$$

acting in the space of grid functions *u h* (*x*), which satisfy the conditions *u h* (*x*) = 0 for all *x* ∈ *Sh* . It is well known that *A x h* is a self-adjoint positive definite operator in *L*2*h* . Applying *A x h* in (31), we can obtain that:

$$\begin{cases} \boldsymbol{u}\_{ttt}^{h}(t, \mathbf{x}) + \boldsymbol{A}\_{\boldsymbol{h}}^{x} \boldsymbol{u}\_{t}^{h}(t, \mathbf{x}) = \boldsymbol{b} \boldsymbol{A}\_{\boldsymbol{h}}^{x} \boldsymbol{u}(t - \boldsymbol{w}, \mathbf{x}) + \boldsymbol{f}^{h}(t, \mathbf{x}),\\ \boldsymbol{0} < t < \infty, \boldsymbol{x} \in \Omega\_{\mathbb{H}} \\\ \boldsymbol{u}^{h}(t, \mathbf{x}) = \boldsymbol{g}^{h}(t, \mathbf{x}), -\boldsymbol{w} \le t \le \mathbf{0}, \mathbf{x} \in \overline{\Omega}\_{\mathbb{H}}. \end{cases} \tag{33}$$

In the second stage, we also get the difference scheme as the one-dimensional problem case:

 *u h k*+2 (*x*)−3*u h k*+1 (*x*)+3*u h k* (*x*)−*u h k*−1 (*x*) *τ* 3 + *A x h u h k*+2 (*x*)−*u h k*+1 (*x*) *τ* = *bAx h u h k*−*N* (*x*) + *f h k* (*x*), *f h k* (*x*) = *f h* (*tk* , *x*), *k* ≥ 1, *x* ∈ Ω*h* , *u h k* (*x*) = *g h* (*tk* , *x*), −*N* ≤ *k* ≤ 0, (*Ih* + *τ* 2*A x h* ) *u h* 1 (*x*)−*u h* 0 (*x*) *τ* = *g h t* (0, *x*), (*Ih* + *τ* 2*A x h* ) *u h* 2 (*x*)−2*u h* 1 (*x*)+*u h* 0 (*x*) *τ* 2 = *g h tt*(0, *x*), *x* ∈ Ω*h* , (*Ih* + *τ* 2*A x h* ) *u h mN*+1 (*x*)−*u h mN*(*x*) *τ* = *u h mN*(*x*)−*u h mN*−1 (*x*) *τ* , *x* ∈ Ω*h* (*Ih* + *τ* 2*A x h* ) *u h mN*+2 (*x*)−2*u h mN*+1 (*x*)+*u h mN*(*x*) *τ* 2 = *u h mN*(*x*)−2*u h mN*−1 (*x*)+*u h mN*−2 (*x*) *τ* 2 , *x* ∈ Ω*h* , *m* = 1, 2, ... (34)

**Theorem 4.** *The solution of DS (34) obeys the following stability estimates:*

$$\begin{split} \max\_{mN \le k \le (m+1)N-2} \left\| \frac{u\_{k+2}^{h} - 2u\_{k+1}^{h} + u\_{k}^{h}}{\tau^{2}} \right\|\_{W\_{2h}^{1}} + \max\_{mN+1 \le k \le (m+1)N} \left\| \frac{u\_{k}^{h} - u\_{k-1}^{h}}{\tau} \right\|\_{W\_{2h}^{2}} \\ + \frac{1}{2} \max\_{mN+1 \le k \le (m+1)N} \|u\_{k}^{h}\|\_{W\_{2h}^{3}} \le C\_{2} \left[ (2 + \tau |b|(N-2))^{m} b\_{0}^{h} \right. \\ \left. \quad + \sum\_{j=1}^{m} (2 + \tau |b|(N-2))^{m-j} \tau \sum\_{s=(j-1)N+1}^{jN} \|A^{\frac{1}{2}} f(t\_{k})\|\_{W\_{2h}^{1}} \right], m = 0, 1, \dots, \\ b\_{0}^{h} = \max \left\{ \max\_{-N \le k \le 0} \|A^{\frac{1}{2}} g\_{t}^{h}(t\_{k})\|\_{W\_{2h'}^{1}} \max\_{-N \le k \le 0} \|g\_{t}^{h}(t\_{k})\|\_{W\_{2h'}^{2}} \max\_{-N \le k \le 0} \|g^{h}(t\_{k})\|\_{W\_{2h}^{3}} \right\}. \end{split}$$

*k*

**Proof.** DS (34) can be written in abstract form:

$$\begin{cases} \frac{u\_{k+2}^h - 3u\_{k+1}^h + 3u\_k^h - u\_{k-1}^h}{\tau^3} + A\_h \frac{u\_{k+2}^h - u\_{k+1}^h}{\tau} = bA\_h u\_{k-N}^h + f\_k^h, k \ge 1, \\\ u\_k^h = g\_{k'}^h - N \le k \le 0, \\\ (I\_h + \tau^2 A\_h) \frac{u\_1^h - u\_0^h}{\tau} = g\_t^h(0), (I\_h + \tau^2 A\_h) \frac{u\_2^h - 2u\_1^h + u\_0^h}{\tau^2} = g\_{tt}^h(0), \\\ (I\_h + \tau^2 A\_h) \frac{u\_{mN+1}^h - u\_{mN}^h}{\tau} = \frac{u\_{mN}^h - u\_{mN-1}^h}{\tau}, \\\ (I\_h + \tau^2 A\_h) \frac{u\_{mN+2}^h - 2u\_{mN+1}^h + u\_{mN}^h}{\tau^2} = \frac{u\_{mN}^h - 2u\_{mN-1}^h + u\_{mN-2}^h}{\tau^2}, m = 1, 2, \dots \end{cases}$$

in a Hilbert space *L*2*h* = *L*2(Ω*h* ) with self-adjoint positive definite operator *Ah* = *A x h* by Formula (32). Here, *g h k* = *g h k* (*x*), *f h k* = *f h k* (*x*), and *u h k* = *u h k* (*x*) are known and unknown abstract mesh functions defined on Ω*h* with the values in *L*2*h* . Then, the estimate of Theorem 4 follows from Estimate (25) and the following theorem.

**Theorem 5.** *The solution of the difference elliptic problem: [40]*

$$A\_h^\mathbf{x} u^h(\mathbf{x}) = \omega^h(\mathbf{x}), \; \mathbf{x} \in \Omega\_h; u^h(\mathbf{x}) = 0, \; \mathbf{x} \in \mathcal{S}\_h$$

*Symmetry* **2020**, *12*, 1033

*obeys the estimate:*

$$\sum\_{r=1}^n \left\| \left. u^h \right|\_{\mathcal{X}\_r \mathcal{X}\_7^r} \right\|\_{L\_{2h}} \le \mathsf{C}\_3 ||\omega^h||\_{L\_{2h}}.$$

*where C*3 *does not depend on h and ωh* .

Third, in [0, ∞) × Ω, the mixed problem for DPDE with the Neumann boundary condition:

$$\begin{cases} u\_{\rm{ff}}(t, \mathbf{x}) - \sum\_{r=1}^{n} \left( a\_{r}(\mathbf{x}) u\_{\rm{tr}\_{r}}(t, \mathbf{x}) \right)\_{\mathbf{x}\_{r}} + \delta u\_{\rm{f}}(t, \mathbf{x}) \\ = b \left( - \sum\_{r=1}^{n} \left( a\_{r}(\mathbf{x}) u\_{\rm{tr}\_{r}}(t - w, \mathbf{x}) \right)\_{\mathbf{x}\_{r}} + \delta u(t - w, \mathbf{x}) \right), \\ 0 < t < \infty, \mathbf{x} \in \Omega, \\ \frac{\partial u(t, \mathbf{x})}{\partial p} = \mathbf{0}, \ \mathbf{x} \in \mathcal{S}, \ 0 \le t < \infty \\ u(t, \mathbf{x}) = g(t, \mathbf{x}), -w \le t \le 0, \mathbf{x} \in \overline{\Omega} \end{cases} \tag{35}$$

is investigated. Here, −→*p* is the normal vector to *S*. Under compatibility conditions, Problem (31) has a unique solution *u*(*t*, *x*) for the given smooth functions *ar*(*x*) ≥ *a* > 0, (*x* ∈ Ω), *g*(*t*, *x*), −*w* ≤ *t* ≤ 0, *x* ∈ Ω, *f*(*t*, *x*), 0 < *t* < ∞, *x* ∈ Ω and *b* ∈ *R* 1 .

The construction of full discretization to Problem (35) is completed in two stages. In the first stage, we introduce the second order difference operator *A x h* defined by:

$$A\_h^\mathbf{x} u^h = -\sum\_{r=1}^n \left(\alpha\_r(\mathbf{x}) u\_{\mathbf{x}\_r}^h\right)\_{\mathbf{x}\_r, j\_r} + \delta u^h \,\tag{36}$$

acting in the space of grid functions *u h* (*x*) that satisfy the conditions *Dhu h* (*x*) = 0 for all *x* ∈ *Sh* . Here, *Dh* is the approximation of operator *∂ ∂* −→*p* . It is known that *A x h* is the self-adjoint positive definite operator in *L*2*h* . Using the difference operator *A x h* , we get the initial value problem (33). Therefore, in the second stage, we use DS (4) for Problem (33):

 *u h k*+2 (*x*)−3*u h k*+1 (*x*)+3*u h k* (*x*)−*u h k*−1 (*x*) *τ* 3 + *A x h u h k*+2 (*x*)−*u h k*+1 (*x*) *τ* = *bAx h u h k*−*N* (*x*) + *f h k* (*x*), *f h k* (*x*) = *f h* (*tk* , *x*), *k* ≥ 1, *x* ∈ Ω*h* , *u h k* (*x*) = *g h* (*tk* , *x*), −*N* ≤ *k* ≤ 0, (*Ih* + *τ* 2*A x h* ) *u h* 1 (*x*)−*u h* 0 (*x*) *τ* = *g h t* (0, *x*), (*Ih* + *τ* 2*A x h* ) *u h* 2 (*x*)−2*u h* 1 (*x*)+*u h* 0 (*x*) *τ* 2 = *g h tt*(0, *x*), *x* ∈ Ω*h* , (*Ih* + *τ* 2*A x h* ) *u h mN*+1 (*x*)−*u h mN*(*x*) *τ* = *u h mN*(*x*)−*u h mN*−1 (*x*) *τ* , *x* ∈ Ω*h* (*Ih* + *τ* 2*A x h* ) *u h mN*+2 (*x*)−2*u h mN*+1 (*x*)+*u h mN*(*x*) *τ* 2 = *u h mN*(*x*)−2*u h mN*−1 (*x*)+*u h mN*−2 (*x*) *τ* 2 , *x* ∈ Ω*h* , *m* = 1, 2, ... (37)

**Theorem 6.** *The solution of the difference scheme (37) obeys the stability estimates in Theorem 4.*

**Proof.** DS (37) can be written in abstract form:

$$\begin{cases} \frac{u\_{k+2}^h - 3u\_{k+1}^h + 3u\_k^h - u\_{k-1}^h}{\tau^3} + A\_h \frac{u\_{k+2}^h - u\_{k+1}^h}{\tau} = bA\_h u\_{k-N}^h + f\_k^h, k \ge 1, \\\ u\_k^h = \mathcal{g}\_{k'}^h - \mathcal{N} \le k \le 0, \\\ (I\_h + \tau^2 A\_h) \frac{u\_1^h - u\_0^h}{\tau} = \mathcal{g}\_t^h(0), (I\_h + \tau^2 A\_h) \frac{u\_2^h - 2u\_1^h + u\_0^h}{\tau^2} = \mathcal{g}\_{tt}^h(0), \\\ (I\_h + \tau^2 A\_h) \frac{u\_{mN+1}^h - u\_{mN}^h}{\tau} = \frac{u\_{mN}^h - u\_{mN-1}^h}{\tau}, \\\ (I\_h + \tau^2 A\_h) \frac{u\_{mN+2}^h - 2u\_{mN+1}^h + u\_{mN}^h}{\tau^2} = \frac{u\_{mN}^h - 2u\_{mN-1}^h + u\_{mN-2}^h}{\tau^2}, m = 1, 2, \dots \end{cases}$$

in a Hilbert space *L*2*h* = *L*2(Ω*h* ) with self-adjoint positive definite operator *Ah* = *A x h* by Formula (36). Here, *g h k* = *g h k* (*x*), *f h k* = *f h k* (*x*), and *u h k* = *u h k* (*x*) are known and unknown abstract mesh functions defined on Ω*h* with the values in *L*2*h* . Therefore, the estimate of Theorem 6 follows from the estimate (25) and the following theorem.

**Theorem 7.** *The solution of the elliptic difference problem: [40]*

$$A\_h^\mathbf{x} u^h(\mathfrak{x}) = \omega^h(\mathfrak{x}), \ \mathfrak{x} \in \Omega\_h; D^h u^h(\mathfrak{x}) = 0, \ \mathfrak{x} \in \mathcal{S}\_h$$

*satisfies the estimate:*

$$\sum\_{r=1}^n \left\| \left. \boldsymbol{\mu}^h \right|\_{\boldsymbol{\mathcal{X}}\_r \boldsymbol{\mathcal{X}}\_r^\*} \right\|\_{L\_{2h}} \le \mathsf{C}\_4 ||\boldsymbol{\omega}^h||\_{L\_{2h'}} $$

*where C*4 *is independent of h and ωh* .

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

It is well known that when the analytical methods fail to work properly, the numerical methods for getting partial differential equations' approximate solutions play a vital role in applied mathematics. In the operator approach, constants in theorems can be large; therefore, in this case a nice stability result must be supported numerically. For this reason, it is important to see that for such a type of theoretical result, we need numerical applications when one cannot know concrete values of constants in stability estimates. Therefore, the first order of accuracy DSs for the solution of one- and two-dimensional DPDEs are presented. To solve this problem, a procedure of modified Gauss elimination is applied. The result of the numerical experiment supports the theoretical statements for the solution of these DSs.

#### *4.1. One-Dimensional Problem*

First, we consider the mixed problem, with the exact solution *u* (*t*, *x*) = *e* −*t* cos *x*,

$$\begin{cases} u\_{ttt}(t, \mathbf{x}) - u\_{t\mathbf{x}\mathbf{x}}(t, \mathbf{x}) = -0.1u\_{\mathbf{x}\mathbf{x}}(t - 1, \mathbf{x}) \\ -2e^{-t}\cos\mathbf{x} - 0.1e^{-(t - 1)}\cos\mathbf{x}, t > 0, \; 0 < \mathbf{x} < \pi, \\\ u\_{\mathbf{x}}(t, 0) = u\_{\mathbf{x}}(t, \pi) = 0, \; 0 \le t < \infty, \\\ u(t, \mathbf{x}) = e^{-t}\cos\mathbf{x}, \; -1 \le t \le 0, \; 0 \le \mathbf{x} \le \pi \end{cases} \tag{38}$$

for the one-dimensional DPDE.

Applying DS (4), we get the following DS:

$$\begin{cases} \begin{aligned} &u\_{n}^{k+2} - 3u\_{n}^{k+1} + 3u\_{n}^{k} - u\_{n}^{k-1} \\ &- \sum\_{n\_{n}+1}^{\tau 3}^{\tau 3} - u\_{n+1}^{k+1} - 2\left(u\_{n}^{k+2} - u\_{n}^{k+1}\right) + u\_{n-1}^{k+2} - u\_{n-1}^{k+1} \\ &- 0.1e\_{n+1} - 2u\_{n}^{k-N} - u\_{n-1}^{k} - 2e^{-t\_{k}}\cos x\_{n} - 0.1e^{-(t\_{k-N})}\cos x\_{n}, \\ &t\_{k} = k\tau, \; IN + 1 \le k \le (l+1)N - 2, \; l = 0, 1, \dots \; 1 \le k \le N-1, \\ &N\tau = 1, \; x\_{n} = n\Lambda, \; 1 \le n \le M-1, \; Mk = \tau, \\ &u\_{n}^{k} = e^{-t\_{k}}\cos x\_{n}, \; -N \le k \le 0, \; 0 \le n \le M, \\ &u\_{n}^{k+1} - u\_{n}^{k} = -e^{-t\_{k}}\cos x\_{n}, \; -N \le k \le 0, \; 0 \le n \le M, \\ &u\_{n}^{k+2} - 2u\_{n}^{k+1} + u\_{n}^{k} = e^{-t\_{k}}\cos x\_{n}, \; -N \le k \le 0, \; 0 \le n \le M, \\ &u\_{1}^{k} - u\_{0}^{k} = u\_{M}^{k} - u\_{M-1}^{k} = 0, \; 0 \le k \le N. \end{aligned} \tag{39}$$

It can be written as the second order difference problem with matrix coefficients:

$$A\ u\_{n+1} + Bu\_n + Cu\_{n-1} = D\varphi\_{\nu\_\nu} \ 1 \le n \le M-1; \newline u\_0 = u\_1, \ \newline u\_M = u\_{M-1}.\tag{40}$$

Here and in the future, we put:

$$u\_s = \begin{bmatrix} u\_s^{IN} \\ \vdots \\ u\_s^{(l+1)N} \end{bmatrix}\_{(N+1)\times 1}, \quad s = n, n \pm 1.$$

Here,

$$A = \mathbb{C} = \begin{bmatrix} 0 & 0 & 0 & 0 & \cdot & 0 & 0 & 0 \\ 0 & 0 & a & -a & \cdot & 0 & 0 & 0 \\ 0 & 0 & 0 & a & \cdot & 0 & 0 & 0 \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 0 & 0 & 0 & 0 & \cdot & 0 & a & -a \\ 0 & 0 & 0 & 0 & \cdot & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \cdot & 0 & 0 & 0 \end{bmatrix}\_{(N+1)\times(N+1)}$$

$$B = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & \cdot & 0 & 0 & 0 & 0 \\ b & -3b & 3b - c & b - c & 0 & \cdot & 0 & 0 & 0 & 0 \\ 0 & b & -3b & 3b - c & c - b & \cdot & 0 & 0 & 0 & 0 \\ 0 & 0 & b & -3b & 3b - c & \cdot & 0 & 0 & 0 & 0 \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ 0 & 0 & 0 & 0 & 0 & \cdot & 3b - c & c - b & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \cdot & -3b & 3b - c & c - b & 0 \\ 0 & 0 & 0 & 0 & 0 & \cdot & b & -3b & 3b - c & c - b \\ -\frac{1}{\tau} & \frac{1}{\tau} & 0 & 0 & 0 & \cdot & 0 & 0 & 0 & 0 \\ \frac{1}{\tau^2} & -\frac{2}{\tau^2} & \frac{1}{\tau^2} & 0 & 0 & \cdot & 0 & 0 & 0 & 0 \\ \end{bmatrix}\_{(N+1)\times(N+1)}$$

where:

$$a = \frac{1}{\pi h^2}, \; b = -\frac{1}{\pi^3}, \; c = \frac{2}{\pi h^2}.$$

$$\boldsymbol{\varrho}\_{n} = \left[ \begin{array}{c} \boldsymbol{\varrho}\_{n}^{IN} \\ \boldsymbol{\varrho}\_{n}^{IN} \\ \vdots \\ \boldsymbol{\varrho}\_{n}^{(l+1)N} \end{array} \right]\_{(N+1)\times 1}, \quad \left\{ \begin{array}{c} \boldsymbol{\varrho}\_{n}^{IN} = \cos \mathbf{x}\_{n} \ -M \le k \le 0, \ 0 \le n \le M, \\\\ \boldsymbol{\varrho}\_{n}^{k} = \boldsymbol{f}(t\_{k}, \mathbf{x}\_{n}) = -0.1 \frac{\boldsymbol{u}\_{n+1}^{k-N} - 2\boldsymbol{u}\_{n}^{k-N} - \boldsymbol{u}\_{n-1}^{k-N}}{h^{2}} \\ -2e^{-t\_{k}} \cos \mathbf{x}\_{n} - 0.1e^{-(t\_{k-N})} \cos \mathbf{x}\_{n}, \\\\ \boldsymbol{t}\_{k} = k\pi, lN+1 \le k \le (l+1)N-2, \\\\ \boldsymbol{l} = 0, 1, \cdots, \; 1 \le n \le M-1, \\\\ \boldsymbol{\varrho}\_{n}^{(l+1)N-1} = -\cos \mathbf{x}\_{n} - M \le k \le 0, \ 0 \le n \le M, \\\\ \boldsymbol{\varrho}\_{n}^{(l+1)N} = \cos \mathbf{x}\_{n}, -M \le k \le 0, \ 0 \le n \le M, \end{array} \right\}$$

and *D* = *IN*+1 is the identity matrix.

To solve this second order difference problem, we use the following formula:

$$
\mu\_n = \mu\_{n+1}\mu\_{n+1} + \beta\_{n+1} \quad n = M-1, \dots, 1, 0,\tag{41}
$$

where *uM* = (*I* − *αM*) −1 *βM*, *αj* (*j* = 1, ..., *M* − 1) are (*N* + 1) × (*N* + 1) square matrices, *βj* (*j* = 1, ..., *M* − 1) are (*N* + 1) × 1 column matrices, *α*1 is the identity, and *β*1 is zero matrices, and:

$$\mathfrak{a}\_{n+1} = -\left(\mathcal{B} + \mathbb{C}\mathfrak{a}\_{\mathfrak{n}}\right)^{-1} A\_{\prime} \tag{42}$$

$$\mathfrak{f}\_{n+1} = \left(\mathcal{B} + \mathbb{C}\mathfrak{a}\_{\mathfrak{n}}\right)^{-1} \left(D\mathfrak{o}\_{\mathfrak{n}} - \mathbb{C}\mathfrak{f}\_{\mathfrak{n}}\right), n = 1, \dots, M-1.$$

The errors are computed by:

$$E\_M^N = \max\_{lN+1 \le k \le (l+1)N-1, 1 \le n \le M-1} \left| u(t\_k, \mathbf{x}\_n) - u\_n^k \right| \tag{43}$$

of the numerical solutions, where *u*(*tk* , *xn*) represents the exact solution and *u k n* represents the numerical solution at (*tk* , *xn*), and the results are given in Table 1.

**Table 1.** Errors of difference scheme (DS) (39).

#### *4.2. Two-Dimensional Problem*

Second, the mixed problem with the Dirichlet condition:

$$\begin{cases} \begin{aligned} &u\_{\rm{III}}(t,\mathbf{x},\mathbf{y}) - u\_{\rm{Ixx}}(t,\mathbf{x},\mathbf{y}) - u\_{\rm{Iyy}}(t,\mathbf{x},\mathbf{y}) = -0.1u\_{\rm{Xx}}(t-1,\mathbf{x}) - 0.1u\_{\rm{Jy}}(t-1,\mathbf{x},\mathbf{y}) \\ &-3e^{-t}\sin x \sin y - 0.2e^{-(t-1)}\sin x \sin y, \; t > 0, \; 0 < \mathbf{x}, \mathbf{y} < \pi, \\ &u(t,0,\mathbf{y}) = u(t,\pi,\mathbf{y}) = 0, \; 0 \le t < \infty, \; 0 \le \mathbf{y} \le \pi, \\ &u(t,\mathbf{x},0) = u(t,\mathbf{x},\pi) = 0, \; 0 \le t < \infty, \; 0 \le \mathbf{x} \le \pi \\ &u(t,\mathbf{x},\mathbf{y}) = e^{-t}\sin x \sin y, \; -1 \le t \le 0, \; 0 \le \mathbf{x}, \mathbf{y} \le \pi \end{aligned} \end{cases} \tag{44}$$

for the two-dimensional DPDE is considered. The exact solution of Problem (44) is *u* (*t*, *x*, *y*) = *e* −*t* sin *x* sin *y*.

Applying DS (4) to the problem (44), we get the following DS of the first order of accuracy in *t*:

 *u k*+2 *n*,*m* − 3*u k*+1 *n*,*m* + 3*u k n*,*m* − *u k*−1 *n*,*m τ* 3 − *u k*+2 *n*+1,*m* − *u k*+1 *n*+1,*m* − 2 *u k*+2 *n*,*m* − *u k*+1 *n*,*m* + *u k*+2 *n*−1,*m* − *u k*+1 *n*−1,*m τh* 2 − *u k*+2 *n*,*m*+1 − *u k*+1 *n*,*m*+1 − 2 *u k*+2 *n*,*m* − *u k*+1 *n*,*m* + *u k*+2 *n*,*m*−1 − *u k*+1 *n*,*m*−1 *τh* 2 = −0.1 *u k*−*N n*+1,*m* − 2*u k*−*N n*,*m* − *u k*−*N n*−1,*m h* 2 − *u k*−*N n*,*m*+1 − 2*u k*−*N n*,*m* − *u k*−*N n*,*m*−1 *h* 2 −3*e* −*tk* sin *xn* sin *ym* − 0.2*e* −*tk*−*N* sin *xn* sin *ym*, *tk* = *kτ*, *lN* + 1 ≤ *k* ≤ (*l* + 1)*N* − 2, *l* = 0, 1, ..., 1 ≤ *k* ≤ *N* − 1 , *Nτ* = 1, *xn* = *nh*, *ym* = *mh* 1 ≤ *n*, *m* ≤ *M* − 1, *Mh* = *π*, *u k n*,*m* = sin *xn* sin *ym*, −*N* ≤ *k* ≤ 0, 0 ≤ *n*, *m* ≤ *M*, *u k*+1 *n*,*m* − *u k n*,*m τ* = sin *xn* sin *ym*, −*N* ≤ *k* ≤ 0, 0 ≤ *n*, *m* ≤ *M*, *u k*+2 *n*,*m* − 2*u k*+1 *n*,*m* + *u k n*,*m τ* 2 = sin *xn* sin *ym*, −*N* ≤ *k* ≤ 0, 0 ≤ *n*, *m* ≤ *M*, *u k* 0,*m* = *u k M*,*m* = 0, 0 ≤ *k* ≤ *N*, 0 ≤ *m* ≤ *M*, *u k n*,0 = *u k n*,*M* = 0, 0 ≤ *k* ≤ *N*, 0 ≤ *n* ≤ *M*. (45)

It can be written as the second order difference problem with the matrix coefficients' form:

$$A\,u\_{n+1} + Bu\_n + Cu\_{n-1,\*} = D\varphi\_n, \quad 1 \le n \le M-1; \\ u\_0 = 0, \; u\_M = 0. \tag{46}$$

Here and in the future, we put:

$$\mathbf{u}\_{\mathbf{s}} = \left[ \mathbf{u}\_{0,\mathbf{s'}}^{lN} \cdots , \mathbf{u}\_{0,\mathbf{s}}^{(l+1)N}, \mathbf{u}\_{1,\mathbf{s'}}^{lN} \cdots , \mathbf{u}\_{1,\mathbf{s}}^{(l+1)N} \right.\\ \left. \text{ } \cdots \text{ } \right., \left. \mathbf{u}\_{M,\mathbf{s'}}^{lN} \cdots , \mathbf{u}\_{M,\mathbf{s}}^{(l+1)N} \right]\_{(N+1)(M+1)\times 1}^{T} \text{ s} = \mathbf{n}\_{\mathbf{s}} n \pm 1 \ldots n \pm 1.$$

*A*, *B*, *C*, *I* are (*N* + 1)(*M* + 1) × (*N* + 1)(*M* + 1) square matrices, and *I*, *R* are identity matrices.

Here,

$$a = \frac{1}{\tau h^2}, \ b = -\frac{1}{\tau^3}, \ c = \frac{2}{\tau h^2},$$

$$A = \mathbb{C} = \begin{bmatrix} 0 & 0 & \cdot & 0 & 0 \\ 0 & E & \cdot & 0 & 0 \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ 0 & 0 & \cdot & E & 0 \\ 0 & 0 & \cdot & 0 & 0 \end{bmatrix},$$

$$B = \begin{bmatrix} Q & O & O & O & \cdot & O & O & O & O \\ E & D & E & O & \cdot & O & O & O & O \\ O & E & D & E & \cdot & O & O & O & O \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ O & O & O & O & \cdot & 0 & E & D & E \\ O & O & O & O & O & O & O & O & O \end{bmatrix},$$

where:

$$E = \begin{bmatrix} 0 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 \\ 0 & 0 & a & -a & \cdots & 0 & 0 & 0 \\ 0 & 0 & a & -a & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & 0 & a & -a \\ 0 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 \end{bmatrix},$$

$$D = \begin{bmatrix} 1 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 \\ b & -3b & 3b - c & b - c & \cdots & 0 & 0 & 0 \\ 0 & b & -3b & 3b - c & \cdots & 0 & 0 & 0 \\ 0 & 0 & b & -3b & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & c - b & 0 & 0 \\ 0 & 0 & 0 & 0 & \cdots & -3b & 3b - c & c - b & 0 \\ 0 & 0 & 0 & 0 & \cdots & -3b & 3b - c & c - b \\ -\frac{1}{\tau} & \frac{1}{\tau} & 0 & 0 & \cdots & 0 & 0 & 0 \\ \frac{1}{\tau^2} & -\frac{2}{\tau^2} & \frac{1}{\tau^2} & 0 & \cdots & 0 & 0 & 0 \end{bmatrix}.$$

*Q* = *I* (*N*+1)×(*N*+1) , *O* = *O*(*N*+1)×(*N*+1) ,

*ϕn* = *ϕ lN* 0,*n* . . . *ϕ* (*l*+1)*N* 0,*n ϕ lN* 1,*n* . . . *ϕ* (*l*+1)*N* 1,*n* . . . *ϕ lN M*,*n* . . . *ϕ* (*l*+1)*N M*,*n* (*M*+1)(*N*+1)×1 , *ϕ lN m*,*n* = sin *xn* sin *ym*, −*M* ≤ *k* ≤ 0, 0 ≤ *n*, *m* ≤ *M*, *ϕ k m*,*n* = −0.1 *u k*−*N n*+1,*m* − 2*u k*−*N n*,*m* − *u k*−*N n*−1,*m h* 2 − *u k*−*N n*,*m*+1 − 2*u k*−*N n*,*m* − *u k*−*N n*,*m*−1 *h* 2 −3*e* −*tk* sin *xn* sin *ym* −0.2*e* −*tk*−*N* sin *xn* sin *ym*, *tk* = *kτ*, *lN* + 1 ≤ *k* ≤ (*l* + 1)*N* − 2, *l* = 0, 1, · · · , 1 ≤ *n*, *m* ≤ *M* − 1, *ϕ* (*l*+1)*N*−1 *m*,*n* = − sin *xn* sin *ym*, −*M* ≤ *k* ≤ 0, 0 ≤ *n*, *m* ≤ *M*, *ϕ* (*l*+1)*N m*,*n* = sin *xn* sin *ym*, −*M* ≤ *k* ≤ 0, 0 ≤ *n*, *m* ≤ *M*,

To solve this second order difference problem, we use the following formula:

$$
\mu\_n = \alpha\_{n+1}\mu\_{n+1} + \beta\_{n+1}, \quad n = M-1, \dots, 1,\\
\mu\_M = 0,\tag{47}
$$

where *αj* (*j* = 1, ..., *M* − 1) are (*N* + 1)(*M* + 1) × (*N* + 1)(*M* + 1) square matrices, *βj* (*j* = 1, ..., *M* − 1) are (*N* + 1)(*M* + 1) × 1 column matrices, and *α*1 and *β*1 are zero matrices and:

$$\mathfrak{a}\_{n+1} = -\left(\mathcal{B} + \mathbb{C}\mathfrak{a}\_{\mathfrak{n}}\right)^{-1} A\_{\mathfrak{n}\prime} \tag{48}$$

$$\mathfrak{F}\_{n+1} = \left(\mathcal{B} + \mathbb{C}\mathfrak{a}\_{\mathfrak{n}}\right)^{-1} \left(\mathcal{D}\mathfrak{q}\_{\mathfrak{n}} - \mathbb{C}\mathfrak{f}\_{\mathfrak{n}}\right), n = 1, \dots, \mathcal{M} - 1.$$

The errors are computed by:

$$lE\_M^N = \max\_{lN+1 \le k \le (l+1)N-1, 1 \le n, m \le M-1} \left| u(t\_k, \mathbf{x}\_{n\nu} y\_m) - u\_{n,m}^k \right| \tag{49}$$

of the numerical solutions, where *u*(*tk* , *xn*, *ym*) represents the exact solution and *u k n*,*m* represents the numerical solution at (*tk* , *xn*, *ym*), and the results are given in Table 2.

**Table 2.** Errors of the difference scheme (45).

As seen in Tables 1 and 2, we obtained some numerical results. If *M* and *N* are doubled, the values of the errors decrease by a factor of approximately 1/2 for DS (39) and (45), respectively.

#### **5. Conclusions**

$$\begin{cases} \ u\_{ttt}(t) + Au\_t(t) = f\left(t, u(t - w)\right), \ 0 < t < \infty, \\\ u(t) = g(t), -w \le t \le 0 \end{cases} \tag{50}$$

in a Hilbert space *H* with an unbounded operator *A*. Moreover, applying the method of [28], we can investigate the convergence of DSs for the numerical solution of Problem (50).

**Author Contributions:** Investigation, A.A., E.H. and S.I.. All authors read and approved the final version of this manuscript.

**Acknowledgments:** The authors thank the reviewers and A. L.Skubachevskii (RUDN) and Ch. Ashyralyyev (CU) for their helpful suggestions.

**Funding:** This research was funded by "Russian Foundation for Basic Research (RFBR) grant number 16–01–00450.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

**Sample Availability:** Samples of the compounds are available from the authors.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Symmetry* Editorial Office E-mail<EMAIL_ADDRESS>www.mdpi.com/journal/symmetry

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com ISBN 978-3-0365-3383-4