$$D^\alpha(\Psi(z)) \simeq P^\alpha \Psi(z),$$

*where Pα is the* (*m* + 1) × (*m* + 1) *operational matrix of the fractional derivative of order α in the Caputo sense and is defined as follows:*

$$P^{\mathfrak{a}} = \begin{pmatrix} 0 & 0 & \dots & 0 \\ \vdots & \vdots & \dots & \vdots \\ 0 & 0 & \dots & 0 \\ 0 & 0 & \dots & 0 \\ \sum\_{k=0}^{i-\lceil\frac{\alpha}{2}\rceil} \xi\_{i,0,k} & \sum\_{k=0}^{i-\lceil\frac{\alpha}{2}\rceil} \xi\_{i,1,k} & \dots & \sum\_{k=0}^{i-\lceil\frac{\alpha}{2}\rceil} \xi\_{i,m,k} \\ \vdots & \vdots & \dots & \vdots \\ \sum\_{m=\lceil\frac{\alpha}{2}\rceil}^{m-\lceil\frac{\alpha}{2}\rceil} \xi\_{m,0,k} & \sum\_{k=0}^{m-\lceil\frac{\alpha}{2}\rceil} \xi\_{m,1,k} & \dots & \sum\_{k=0}^{\lceil\frac{\alpha}{2}\rceil} \xi\_{m,m,k} \\ \end{pmatrix} \text{ and } \xi\_{i,j,k} \text{ is given by} $$

$$\tilde{\xi}\_{i,j,k} = \begin{cases} i \sum\_{k=0}^{i-\lceil\frac{\alpha}{2}\rceil} \left(-1\right)^{k} \frac{4^{i-k} \Gamma(2i-k) \Gamma(i-k+1) \Gamma(i-k-m+1/2)}{\sqrt{\pi} \Gamma(k+1) \Gamma(2i-2k+1) \Gamma(i-k-m+1)^2}, & j=0, \\ 2i \sum\_{k=0}^{i-\lceil\frac{\alpha}{2}\rceil} \frac{1}{\pi} \sum\_{r=0}^{j} \frac{(-1)^{k+r}}{\sqrt{\pi}} \frac{4^{i-k} \Gamma(2i-k+1) \Gamma(i-k-m+1)}{\Gamma(k+1) \Gamma(2i-2k+1) \Gamma(i-k-m+1)}, & j=1, 2, 3, \dots \end{cases} \tag{19}$$

**Proof.** Applying the Caputo derivative to (12), we have

$$D^{\mathfrak{a}}(\mathrm{VL}\_i^\*(z)) = D^{\mathfrak{a}}\left(\sum\_{k=0}^i (-1)^k \frac{4^{i-k} 2i\Gamma(2i-k)}{\Gamma(k+1)\Gamma(2i-2k+1)} z^{i-k}\right). \tag{20}$$

Applying the linearity of the Caputo derivative, and using (5), we have

$$D^a(\mathrm{VL}\_i^\*(z)) = \sum\_{k=0}^{i-\lceil a \rceil} (-1)^k \frac{4^{i-k} 2i\Gamma(2i-k)\Gamma(i-k+1)}{\Gamma(k+1)\Gamma(2i-2k+1)\Gamma(i-k+1-a)} z^{i-k-a}, \; i = \lceil a \rceil, \ldots, n,\tag{21}$$

and

$$D^a(\text{VL}\_i^\*(z)) = 0, \quad i = 0, 1, \dots, \lceil a \rceil - 1. \tag{22}$$

Now, approximate *zi*−*k*−*α* by (*m* + 1) terms of the series, as

$$z^{i-k-n} \simeq \sum\_{j=0}^{m} c\_{kj} \operatorname{VL}\_j^\*(z),\tag{23}$$

where

$$\varepsilon\_{kj} = \frac{1}{\delta\_j \pi} \int\_0^1 \frac{u(z) \,\mathrm{VL}\_j^\*(z)}{\sqrt{z - z^2}} \, dz,\tag{24}$$

and *u*(*z*) = *zi*−*k*−*α*.

Now, inserting the value of *u*(*z*)*and* VL∗ *j* (*z*) into Equation (24), we obtain

$$c\_{kj} = \begin{cases} \frac{1}{4\pi} \int\_0^1 \frac{2z^{j-k-x}}{\zeta(z-z^2)} dz, & j=0, \\\frac{1}{2\pi} \int\_0^1 \frac{z^{j-k-x}}{\zeta(z-z^2)} \sum\_{r=0}^j \frac{(-1)^r 2j^{j-r} \Gamma(2j-r)}{\Gamma(r+1)\Gamma(2j-2r+1)} z^{j-r} dx, & j=1,2,3,\ldots, \end{cases}$$

$$= \begin{cases} \frac{1}{2\pi} \int\_0^1 \frac{z^{j-k-x}}{\zeta(z-z^2)} dz, & j=0, \\\frac{1}{\pi} \int\_0^1 \sum\_{r=0}^j \frac{(-1)^r j 4^{j-r} \Gamma(2j-r)}{\Gamma(r+1)\Gamma(2j-2r+1)} \frac{z^{j+j-k-r-x}}{\zeta(z-z^2)} dz, & j=1,2,3,\ldots, \end{cases}$$

$$= \begin{cases} \frac{1}{2\sqrt{\pi}} \frac{\Gamma(i-k-x+1/2)}{\Gamma(i-k-x+1)}, & j=0, \\\frac{1}{\sqrt{\pi}} \sum\_{r=0}^j \frac{(-1)^r j i^{j-r} \Gamma(2j-r)}{\Gamma(r+1)\Gamma(2j-2r)\Gamma(i+j-k-r-x+1/2)}, & j=1,2,3,\ldots,\end{cases} \tag{25}$$

by inserting the value of *zi*−*k*−*α* into Equation (21), we obtain

$$D^a(\operatorname{VL}\_i^\*(z)) \simeq \sum\_{j=0}^m S\_v(i,j) \operatorname{VL}\_j^\*(z)),\tag{26}$$

where *Sv*(*i*, *j*) = ∑*i*−%*α*& *k*=0 *ξi*,*j*,*k*, and

$$\mathcal{I}\_{i,j,k}^{\mathcal{I}} = \begin{cases} (-1)^k \frac{4^{i-k} i \Gamma(2i-k) \Gamma(i-k+1) \Gamma(i-k-a+1/2)}{\sqrt{\pi} \Gamma(k+1) \Gamma(2i-2k+1) \Gamma(i-k-a+1)^2}, & j=0, \\\\ \sum\_{r=0}^j \frac{(-1)^{k+r}}{\sqrt{\pi}} \frac{4^{i-k} 2i \Gamma(2i-k) \Gamma(i-k+1)}{\Gamma(k+1) \Gamma(2i-2k+1) \Gamma(i-k-a+1)} \\ \qquad \qquad \qquad \qquad \times \frac{4^{i-r} j \Gamma(2j-r) \Gamma(i+j-k-r-a+1/2)}{\Gamma(r+1) \Gamma(2j-2r+1) \Gamma(i+j+k-r-a+1)^{\prime}}, & j=1,2,3,\dots \end{cases} \tag{27}$$

Rewriting Equation (26) in vector form, we obtain

$$D^{\mathfrak{a}}(\operatorname{VL}\_{i}^{\ast}(z)) \simeq \left(\sum\_{k=0}^{i-\lceil a \rceil} \mathfrak{z}\_{i,0,k} \sum\_{k=0}^{i-\lceil a \rceil} \mathfrak{z}\_{i,1,k}, \dots, \sum\_{k=0}^{i-\lceil a \rceil} \mathfrak{z}\_{i,m,k} \right) \Psi(z). \tag{28}$$

For simplicity, we can write Equation (28) as:

$$D^{\mathfrak{a}}(\mathrm{VL}\_i^\*(z)) = P^{(a)}\Psi(z),\tag{29}$$

where

$$P^{(a)} = \left(\sum\_{k=0}^{i-\lceil a \rceil} \xi\_{i,0,k}, \sum\_{k=0}^{i-\lceil a \rceil} \xi\_{i,1,k}, \dots, \sum\_{k=0}^{i-\lceil a \rceil} \xi\_{i,m,k}\right). \tag{30}$$

Equations (22) and (29) prove the required result.

### **4. Application of Operational Matrices Method**

In this section, we apply the Vieta–Lucas operational matrix method to find the analytical-approximate solution of linear and nonlinear FDEs.

*4.1. Linear FDEs*

Consider the linear FDE

$$D^{\mu}v(z) = b\_1 D^{\theta\_1} v(z) + b\_2 D^{\theta\_2} v(z) + \dots + b\_k D^{\theta\_k} v(z) + b\_{k+1} v(z) + b\_{k+2} \mathfrak{g}(z),\tag{31}$$

with initial conditions

$$w^{(i)}(0) = d\_{i\prime} \qquad \qquad i = 0, \ldots, n,\tag{32}$$

where *bl*, for *l* = 1, ... , *k* + 2, are real constant coefficients and also *n* < *α* ≤ *n* + 1, and 0 < *ϑ*1 < *ϑ*2 < ... < *ϑk* < *α*.

The unknown function *υ*(*z*) and the source term *g*(*z*) can be approximated as:

$$\upsilon(z) \simeq \sum\_{i=0}^{M} \mathfrak{c}\_i \, \mathrm{VL}\_i^\*(z) = \mathbb{C}^T \Psi(z),\tag{33}$$

$$\log(z) \simeq \sum\_{i=0}^{M} h\_i \,\mathrm{VL}\_i^\*(z) = H^T \Psi(z),\tag{34}$$

where *H* = [*h*0, ... , *hM*] *T* is known, and *C* = [*c*0, ... , *cM*] *T* is an unknown to be determined. Now, using Equations (29) and (33), we have

$$D^a \upsilon(z) \simeq \mathbb{C}^T D^a \Psi(z) \simeq \mathbb{C}^T P^a \Psi(z),\tag{35}$$

$$D^{\theta\_j}v(z) \simeq \mathbb{C}^T D^{\theta\_j} \Psi(z) \simeq \mathbb{C}^T P^{\theta\_j} \Psi(z), \qquad j = 1, \ldots, k. \tag{36}$$

Using Equations (33)–(36), the residual *R*(*x*) for Equation (31) can be written as

$$R\_M(z) \simeq (\mathbb{C}^T P^\mathfrak{a} - b\_1 \mathbb{C}^T P^{\theta\_1} - \dots - b\_k \mathbb{C}^T P^{\theta\_k} - b\_{k+1} \mathbb{C}^T - b\_{k+2} \mathbb{C}^T) \mathbb{Y}(z). \tag{37}$$

Using the spectral tau method [36], a system of linear equations is generated by applying

$$
\langle \mathcal{R}\_M(z), \Psi(z) \rangle = \int\_0^1 \mathcal{R}\_M(z) \Psi(z) dz, \qquad j = 0, 1, \dots, M - n - 1. \tag{38}
$$

Moreover, by substituting Equation (33) in the initial conditions given in Equation (32), we obtain

$$v^{(i)}(0) = \mathbb{C}^T P^{(i)} \Psi(0) = d\_{i\prime} \qquad i = 0, 1, \ldots, n. \tag{39}$$

Equations (38) and (39) generate the (*M* − *n*) and (*n* + 1) set of linear equations, respectively. This system of linear equations can then be solved easily for the unknown coefficients. Consequently, we can approximate *υ*(*z*) given in Equation (33).

### *4.2. Nonlinear FDEs*

Consider the nonlinear fractional-order differential equation

$$F(z, \upsilon(z), D^{\vartheta\_1}\upsilon(z), D^{\vartheta\_2}\upsilon(z), \dots, D^{\vartheta\_k}\upsilon(z)) = 0,\tag{40}$$

with boundary conditions

$$H\_{\dot{l}}(\upsilon(\varsigma\_{\dot{l}}), \upsilon^{'}(\varsigma\_{\dot{l}}), \dots, \upsilon^{(s)}(\varsigma\_{\dot{l}})) = d\_{\dot{l}\prime} \tag{41}$$

where 0 ≤ *s* < *max*{*ϑj*, *j* = 1, . . . , *k*} ≤ *s* + 1, *ςj* ∈ [0, 1], *j* = 0, ... ,*s* and *Hj* are linear combinations of *υ*(*ςj*), *υ* (*ςj*),..., *υ*(*s*)(*ςj*). Now, using Theorem 1 and Equation (33), the terms of Equation (40) can be approximated as

$$\begin{cases} \boldsymbol{\nu}(\boldsymbol{z}) = \mathbb{C}^{T}\boldsymbol{\Psi}(\boldsymbol{z}),\\ \boldsymbol{D}^{\theta\_{1}}\boldsymbol{v}(\boldsymbol{z}) = \mathbb{C}^{T}\boldsymbol{P}^{\theta\_{1}}\boldsymbol{\Psi}(\boldsymbol{z}),\\ \boldsymbol{D}^{\theta\_{2}}\boldsymbol{v}(\boldsymbol{z}) = \mathbb{C}^{T}\boldsymbol{P}^{\theta\_{2}}\boldsymbol{\Psi}(\boldsymbol{z}),\\ \vdots \\ \boldsymbol{D}^{\theta\_{k}}\boldsymbol{v}(\boldsymbol{z}) = \mathbb{C}^{T}\boldsymbol{P}^{\theta\_{k}}\boldsymbol{\Psi}(\boldsymbol{z}).\end{cases} \tag{42}$$

Similarly, the terms of Equation (41) can be approximated as

$$\begin{cases} \boldsymbol{\nu}(\boldsymbol{\xi}\_{\boldsymbol{j}}) = \mathbb{C}^{T}\boldsymbol{\Psi}(\boldsymbol{\xi}\_{\boldsymbol{j}}),\\ \boldsymbol{\nu}'(\boldsymbol{\xi}\_{\boldsymbol{j}}) = \mathbb{C}^{T}P^{(1)}\boldsymbol{\Psi}(\boldsymbol{\xi}\_{\boldsymbol{j}}),\\ \boldsymbol{\nu}''(\boldsymbol{\xi}\_{\boldsymbol{j}}) = \mathbb{C}^{T}P^{(2)}\boldsymbol{\Psi}(\boldsymbol{\xi}\_{\boldsymbol{j}}),\\ \vdots \\ \boldsymbol{\nu}^{(s)}(\boldsymbol{\xi}\_{\boldsymbol{j}}) = \mathbb{C}^{T}P^{(s)}\boldsymbol{\Psi}(\boldsymbol{\xi}\_{\boldsymbol{j}}).\end{cases} \tag{43}$$

In light of (42) and (43), we may write Equations (40) and (41), respectively, as

$$F(z, \mathbb{C}^T \Psi(z), \mathbb{C}^T P^{\theta\_1} \Psi(z), \dots, \mathbb{C}^T P^{\theta\_k} \Psi(z)) = 0,\tag{44}$$

$$H\_{\mathbf{j}}(\mathbb{C}^{T}\Psi(\mathbf{g}\_{\mathbf{j}}), \mathbb{C}^{T}P^{(1)}\Psi(\mathbf{g}\_{\mathbf{j}}), \dots, \mathbb{C}^{T}P^{(s)}\Psi(\mathbf{g}\_{\mathbf{j}})) = d\_{\mathbf{j}}.\tag{45}$$

Now, to find the solution *υ*(*z*), we first collocate Equation (44) at (*M* − *s*) points. These equations, along with Equation (45), generate a system of algebraic equations, which can be solved to find *ci*, *i* = 0, . . . , *M*. Consequently, the function *υ*(*z*) can be approximated.

### **5. Error Estimate**

**Lemma 1.** ([37]) *The following assumptions for the function g*(*z*)*, such that g*(*k*) = *bk, must hold true:*

*Then,*

$$P\_{\mathfrak{m}} \le \int\_{\mathfrak{m}}^{\infty} \mathfrak{g}(z) dz.$$

**Theorem 2.** *If v*(*z*) ∈ *L*2 *w*(Δ), Δ = [0, 1]*, v*(*z*) = ∑∞ *k*=0 *bk* VL∗ *k* (*z*)*, bk is introduced in Equation* (15)*, and v*(*z*) ≤ *N, then we have*

$$\|v(z) - v\_m(z)\|\_w \le \frac{N}{\sqrt{96m^3\pi}}.\tag{46}$$

**Proof.** It is evident that the shifted VLPs are orthogonal on the interval [0, 1] with respect to the weight function, *w*(*z*) = √ 1 *z*−*z*2 . Hence, these polynomials form a complete *L*2 *w*(Δ) orthogonal set, where *L*2 *w*(Δ) represents the space of functions defined as *v* : Δ → R. Thus, the error in space *L*2 *w*(Δ) is determined as

$$\|\|v(z) - v\_m(z)\|\|\_{w}^{2} = \left(\int\_{0}^{1} |v(z) - v\_m(z)|^2 w(z) dz\right)^{\frac{1}{2}}.\tag{47}$$

Using Equations (13) and (17) in Equation (47), we have

$$\|\|v(z) - v\_{\mathcal{W}}(z)\|\|\_{w}^{2} = \left(\int\_{0}^{1} \left|\sum\_{k=m+1}^{\infty} b\_{k} \operatorname{VL}\_{k}^{\*}(z)\right|^{2} w(z) dz\right)^{\frac{1}{2}}.\tag{48}$$

Now, by applying the orthogonality property of shifted VLPs to Equation (48), we have

$$\|v(z) - v\_m(z)\|\_{\text{uv}}^2 = \frac{1}{\delta\_k \pi} \sum\_{k=m+1}^{\infty} |b\_k|^2. \tag{49}$$

Now, by using the substitution 4*x* − 2 = 2 cos(*θ*) in Equation (15), the coefficients, *bk*, *k* = 0, 1, ··· , *m*, can be determined as

$$b\_{k} = \frac{1}{4\delta\_{k}\pi} \int\_{0}^{\pi} v^{\prime\prime} \left(\frac{1+\cos(\theta)}{2}\right) \sin(\theta) \left(\frac{\sin(k-1)\theta}{k-1} - \frac{\sin(k+1)\theta}{k+1}\right) d\theta. \tag{50}$$

Equation (50) can also be expressed as

$$|b\_{k}| = \left| \frac{1}{4\delta\_{k}\pi} \int\_{0}^{\pi} v^{\prime\prime} \left( \frac{1+\cos(\theta)}{2} \right) \sin(\theta) \left( \frac{\sin(k-1)\theta}{k-1} - \frac{\sin(k+1)\theta}{k+1} \right) \right| d\theta. \tag{51}$$

Using *v*(*z*) ≤ *N*, and the properties of trigonometric functions, we may express Equation (51) as

$$|b\_k| \le \frac{N}{4k(k-1)(k+1)}, \ k > 2. \tag{52}$$

Now, using Equation (52) in Equation (49), we have

$$\|v(z) - v\_m(z)\|\_w^2 \le \frac{N^2}{32\pi} \sum\_{k=m+1}^{\infty} \frac{1}{k^4}.\tag{53}$$

Now, by using the Lemma 1, we have

$$\|v(z) - v\_m(z)\|\_w^2 \quad \le \quad \frac{N^2}{32\pi} \int\_m^\infty z^{-4} dz \tag{54}$$

$$\lambda = \frac{N^2}{32\pi} \times \frac{1}{3m^3} = \frac{N^2}{96m^3\pi}.\tag{55}$$

Finally, we have

$$||v(z) - v\_{\rm II}(z)||\_{\rm w} \le \frac{N}{\sqrt{96m^3 \pi}}.\tag{56}$$

### **6. Illustrative Examples**

In this section, we give some numerical examples to show the accuracy of our proposed method.

**Example 1.** *Consider the following fractional Bagley–Torvik equation*

$$a\_1 \mathcal{D}^2 \upsilon(z) + a\_2 \mathcal{D}^{3/2} \upsilon(z) + a\_3 \upsilon(z) = \mathcal{g}(z), \ z \in [0, 1], \tag{57}$$

*subject to the initial conditions with integer order*

$$
\upsilon(0) = 1 = \upsilon(0)'.\tag{58}
$$

*The source term g*(*z*) *is as follows:*

$$g(z) = 1 + z.\tag{59}$$

*The exact solution of the problem in Example 1 is:*

$$
\upsilon(z) = 1 + z.\tag{60}
$$

*Now, we apply the technique that is described in Section 4.1 by choosing the first three terms of VLPs. We may write the approximation solution as*

$$\begin{aligned} \upsilon(z) &= \mathbb{C}^{\Gamma} \Psi(z) \Longleftrightarrow y(0) = \mathbb{C}^{\Gamma} \Psi(0) = d\_0 = 1, \\ \upsilon'(z) &= \mathbb{C}^{\Gamma} H^{v} \Psi(z) \Longleftrightarrow y'(0) = \mathbb{C}^{\Gamma} H^{v} \Psi(0) = d\_1 = 1. \end{aligned} \tag{61}$$

*Now, we have Vieta–Lucas polynomials*

$$\begin{split} \Psi(z) &= (\text{VL}\_0^\*(z), \text{VL}\_1^\*(z), \text{VL}\_2^\*(z)), \\ &= (2, \ 4z - 2, \ 16z^2 - 16z + 2), \end{split} \tag{62}$$

*and*

$$\mathbf{G}^T = \begin{pmatrix} 0.75 \\ 0.25 \\ 0 \end{pmatrix}. \tag{63}$$

*The Vieta–Lucas operational matrices can be expressed as*

$$\begin{array}{rcl} \mathbf{D}^2 &=& \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 16 & 0 & 0 \end{pmatrix}, \\\\ \mathbf{D}^{3/2} &=& \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 11.4936 & 7.6624 & -1.5325 \end{pmatrix}, \\\ \mathbf{D}^1 &=& \begin{pmatrix} 0 & 0 & 0 \\ 2 & 0 & 0 \\ 0 & 8 & 0 \end{pmatrix}. \end{array} \tag{64}$$

*The residuals can be evaluated as*

$$\mathcal{R}(z) \quad = \left(\mathbf{C}\mathbf{D}^2 + a\_2\mathbf{C}\mathbf{D}^{3/2} + a\_2\mathbf{C} - a\_3\mathbf{G}\right)\Psi(z),$$

*where C* = [*c*0, *c*1, *c*2]. *Now, using initial conditions, we have:*

$$\begin{aligned} 2c\_0 - 2c\_1 + 2c\_2 &= 1, \\ 4c\_1 - 16c\_2 &= 1. \end{aligned} \tag{65}$$

*Moreover, using the inner product of the residual with the Vieta–Lucas polynomials, we obtain a system of equations. If we take one equation from this system and two equations from Equation* (65)*, then, by simultaneously solving these equations, we obtain c*0 = 0.75, *c*1 = 0.25, *c*2 = 0, *hence*

$$w(z) = (\frac{3}{4}, \frac{1}{4}, 0) \begin{pmatrix} 2 \\ 4z - 2 \\ 16z^2 - 16z + 2 \end{pmatrix} = 1 + z,\tag{66}$$

*which is the exact solution.*

**Remark 1.** *The numerical results computed using our method are compared with the method of [33] by choosing various n. We observe that our proposed method produces efficient numerical results as compared to the numerical results obtained by using the method of [33] (see Tables 1 and 2 and Figure 1). Moreover, for a small value of n* = 2*, the exact solution and the approximate solution computed by using our method coincide (see Table 1 and Figure 1).*

**Table 1.** Comparison of approximate solution of Example 1.

**Table 2.** Comparison of absolute errors of Example 1.

**Figure 1.** Approximate solutions of Example 1 computed by our method at *n* = 2 is compared with the method of [33] computed at *n* = 10.

**Remark 2.** *The numerical results of Example 1 computed at n* = 2 *by using our method are compared with the results obtained by using the methods of [32,34] at n* = 6*. We observe that, for a small value of n* = 2*, the approximate solution obtained using our method coincides with the exact* *solution of Example 1 (see Tables 3 and 4). However, the exact solution and the approximate solution computed by using the methods of [32,34] coincide at n* = 6*. This shows that our proposed method is numerically more efficient.*

**Table 3.** Comparison of approximate solution of Example 1.

**Table 4.** Comparison of approximate solution of Example 1.

**Example 2.** *Consider the following linear initial value problem [38]:*

$$\begin{aligned} D^a \upsilon(z) + \upsilon(z) &= 0, \ 0 < a < 2, \\ \upsilon(0) &= 1, \ \upsilon'(0) = 0. \end{aligned} \tag{67}$$

*The exact solution of the problem is υ*(*z*) = ∑∞ *k*=0 (−*zα*) *k* Γ(*αk*+1) *[39].*

*To solve the problem, we use the technique described in Section 4.1. The absolute error for α* = 0.85 *and n* = 2,5 *and* 8 *is shown in Table 5. An error plot is also shown in Figure 2 for these values.*

*We can see in Table 5 that a good approximation has been achieved by using some initial terms of VLPs. Moreover, the numerical results for υ*(*z*) *when n* = 10 *and α* = 0.5, 0.65, 0.8, 0.95 *and* 1 *are plotted in Figure 3. The exact solution for α* = 1 *is υ*(*z*) = *exp*(−*z*)*. It can be noted that the numerical solution converges to the analytical solution when α approaches 1. We also analyze the nonlocal behavior of the fractional derivative by computing the results at various non-integer values of α, which highlights the advantage of using the fractional derivatives, as the next state of the system depends not only upon its current state by also upon all of its historical states.*

**Table 5.** Absolute error for *α* = 0.85 and *n* = 2, 5 and 8 in Example 2.

**Figure 2.** Error plots of Example 2 at different scale levels.

**Figure 3.** Exact and approximate solutions of Example 2 are compared at different scale levels. **Example 3.** *Consider the following initial value problem [10]:*

$$\prescript{1}{}{\mathcal{D}}\upsilon(z) = a\_{1\mathsf{C}}\mathcal{D}^{\frac{1}{4}}\upsilon(z) - \upsilon(z) + \mathsf{g}(z), \; z \in [0, 1], \tag{68}$$

*subject to the initial condition with integer order*

$$
\upsilon(0) = 0.\tag{69}
$$

*The source term is given as*

$$g(z) = \frac{5z^{\frac{3}{2}}}{2} + z^{\frac{5}{2}} + \frac{15\sqrt{\pi}z^{\frac{3}{4}}}{8\Gamma(\frac{13}{4})}.$$

*The exact solution at a*1 = −1 *is given below:*

$$v(z) = z^2 \sqrt{z}.$$

*We test the behavior of our proposed method by solving Example* (3) *at various values of n. In Table 6, we list the L*∞ *and L*2 *errors for different values of n. We compare the numerical results obtained by using our proposed method with the numerical results obtained in [10]. It can be observed that the errors computed by using our method are much smaller than those computed by using the method presented in [10]; see Table 6. This highlights the efficiency of our method for this problem. Note that the symbol* " − " *means that the result for n is unavailable for the method [10].*

**Table 6.** Approximate results of Example 3 at various values of *n*.

**Example 4.** *Consider the following nonlinear initial value problem [40]:*

$$\begin{split} D^a v(z) &= \frac{40320}{\Gamma(9-a)} z^{8-a} - 3 \frac{\Gamma(5+a/2)}{\Gamma(5-a/2)} z^{4-a/2} + \frac{9}{4} \Gamma(a+1) + \left( \frac{3}{2} z^{a/2} - z^4 \right)^3 - v(z)^{\frac{3}{2}}, \\ v(0) &= 0, \qquad z^{'}(0) = 0, \quad 0 < a < 2. \end{split} \tag{70}$$

*The exact solution of the problem is υ*(*z*) = *z*8 − 3*z*(4+*α*/2) + 9 4 *zα [39].*

*We have solved the problem using the technique described in Section 4.2. The absolute error for α* = 0.85 *and n* = 2, 5 *and* 8 *is shown in Table 5. An error plot is also shown in Figures 4 and 5 for these values. We can see in Table 7 that a good approximation has been achieved. Numerical results for υ*(*z*) *when n* = 6 *and α* = 0.6, 0.7, 0.8, 0.9 *and* 1 *are plotted in Figure 6, along with the exact solutions at the given values of α. It can be noted that, as α approaches 1, the solution of the FDEs approaches that of the integer-order differential equations. We also analyze the nonlocal behavior of the fractional derivative by computing the results at various non-integer values of α, which highlights the advantage of using the fractional derivatives, as the next state of the system depends not only upon its current state but also upon all of its historical states (Table 8).*

**Table 7.** Absolute error for *α* = 0.85 and *n* = 2, 5 and 8 in Example 4.

**Figure 4.** Error plots of Example 4 at different scale levels.

**Figure 5.** Exact and approximate solutions of Example 4 are compared at *n* = 20.

**Figure 6.** Exact and approximate solutions of Example 4 are compared at different scale levels. **Table 8.** Approximate results of Example 4 at various values of *α*.

**Example 5.** *Consider the following nonlinear initial value problem:*

$$\mathcal{D}^3 v(z) + \mathcal{D}^{5/2} v(z) + v^2(z) = z^4, \ v(0) = v'(0) = 0, \ v''(0) = 2. \tag{71}$$

We solved this problem by using the same technique as described in Section 4.2 with *n* = 3.

The exact solution of the problem is *υ*(*z*) = *z*2, whereas the source term is *g*(*z*) = *z*4. The operational matrices can be expressed as

$$\begin{array}{rcl} \mathbf{D}^3 &=& \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 192 & 0 & 0 & 0 \end{pmatrix}, \\\\ \mathbf{D}^{5/2} &=& \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 137.9229 & 91.9486 & -18.3897 & 7.8813 \end{pmatrix}, \\\\ \mathbf{D}^1 &=& \begin{pmatrix} 0 & 0 & 0 & 0 \\ 2 & 0 & 0 & 0 \\ 0 & 8 & 0 & 0 \\ 6 & 0 & 12 & 0 \end{pmatrix}, \\\\ \mathbf{D}^2 &=& \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 16 & 0 & 0 & 0 \\ 0 & 96 & 0 & 0 \end{pmatrix}, \end{array} \tag{72}$$

where *C* = [*c*0, *c*1, *c*2, *c*3].

Now, using the initial condition, we obtain three equations:

$$\begin{aligned} 2c\_0 - 2c\_1 + 2c\_2 - 2c\_3 &= 0, \\ 4c\_1 - 16c\_2 + 36c\_3 &= 0, \\ 32c\_2 - 192c\_3 &= 2. \end{aligned} \tag{73}$$

Meanwhile, using the technique in Section 4.2, we obtain the following equation:

$$\mathbf{C}^{\mathsf{T}}\mathbf{D}^{(3)}\boldsymbol{\Psi}(z) + \mathbf{C}^{\mathsf{T}}\mathbf{D}^{(\frac{5}{2})}\boldsymbol{\Psi}(z) + \left[\mathbf{C}^{\mathsf{T}}\boldsymbol{\Psi}(z)\right]^2 - z^4 = 0. \tag{74}$$

Now, we collocate Equation (74) at the first root of *P*4(*z*), and we obtain *z*0 = 0.06698. Now, solving Equations (73) and (74), we obtain

$$w(z) = (\frac{3}{16}, \frac{1}{4}, \frac{1}{16}, 0) \begin{pmatrix} 2 \\ 4z - 2 \\ 16z^2 - 16z + 2 \\ 64z^3 - 96z^2 + 36z - 2 \end{pmatrix} \approx z^2,\tag{75}$$

which is the exact solution.

**Example 6.** *Consider the following non-homogenous multi-order fractional problem:*

$$\begin{split} c\mathcal{D}^{\mathfrak{x}}\upsilon(z) &= a\_{\mathbb{C}}\mathcal{D}^{\mathfrak{f}\_{0}}\upsilon(z) + b\_{\mathbb{C}}\mathcal{D}^{\mathfrak{f}\_{1}}\upsilon(z) + c\_{\mathbb{C}}\mathcal{D}^{\mathfrak{f}\_{2}}\upsilon(z) \\ &+ d\_{\mathbb{C}}\mathcal{D}^{\mathfrak{f}\_{3}}\upsilon(z) + \mathfrak{g}(z), \; z \in [0, 1], \; 0 < \mathfrak{a} < 2, \end{split} \tag{76}$$

*subject to the following initial conditions*

$$v(0) = 0, \ v'(0) = 0.$$

The source term is as below:

$$g(z) = 4z - z^2 - \frac{6776}{4503}z^{\frac{3}{2}} + 42z^5 - 14z^6 + z^7 + \frac{1516}{5629}z^{\frac{43}{2}} - 2. \tag{77}$$

The exact solution corresponding to *α* = 2 ,*a* = *c* = −1, *b* = 2, *d* = 0, *β*0 = 0, *β*1 = 1, *β*2 = 1 2 is given below:

$$v(z) = z^{\mathcal{T}} - z^2.$$

We can observe in Example 6 that a good approximation of the function has been achieved while using *n* = 7 as a scale level. The absolute error Table 9 at different scale levels is given below. (see Figure 7).

**Table 9.** Absolute error for *n* = 2, 5 and 8 in Example 6.

**Figure 7.** Graphs of exact solution and approximate solution of Example 6 with *n* = 7.

Moreover, a comparison between the exact solution and approximate solution at different scale levels for the values of *z* is given in Table 10.

**Table 10.** Comparison of exact solution with approximate solution (AS) at *m* = 2, 5 and 8 in Example 6.

### **7. Conclusions**

In the present study, we introduce a new fractional-order derivative operational matrix of VLPs in Caputo sense. The newly derived operational matrix is used to develop a computer-oriented numerical algorithm to solve the linear and nonlinear FDEs that include the Caputo fractional-order derivative. The proposed numerical algorithm has the advantage of transforming the problems into a system of algebraic equations that are easy to solve using any computational software. To the best of our knowledge, this is the first result where the numerical algorithm is presented using the operational matrix of VLPs and the solution of the problems is approximated using its basis vectors.

We tested the accuracy and efficiency of the algorithm by solving various linear and nonlinear FDEs with initial conditions. We found that with an increase in the values of *n*, the approximate solutions were in good agreement with the exact solutions. We also demonstrated the high efficiency of the method by determining the amount of absolute error and observed that as we increased *n*, this amount was decreased significantly. In addition to this, the numerical efficiency was also demonstrated by comparing the results obtained by using our method with results obtained otherwise in the literature [10,32–34]. We observed that our method produced more efficient results.

**Author Contributions:** Conceptualization, I.T. and T.A.; Formal analysis, Z.A.N., I.T. and M.A.A.; Funding acquisition, T.A. and M.A.A.; Investigation, Z.A.N. and T.A.; Methodology, Z.A.N. and I.T.; Resources, T.A.; Software, Z.A.N. and I.T.; Supervision, I.T. and T.A.; Validation, I.T., T.A. and M.A.A.; Visualization, Z.A.N.; Writing—original draft, Z.A.N.; Writing—review and editing, I.T., T.A. and M.A.A. All authors contributed equally to this article. All authors read and approved the final manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** Prince Nourah bint Abdulrahman University Researchers supporting project number (PNURSP2022R14), Princess Nourah bint Abdulrahman University, Riyadh, Saudi Arabia. The author T. Abdeljawad would like to thank Prince Sultan University for paying the APC and for the support through the TAS research lab.

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

### **References**

### *Article* **Numerical Analysis of Time-Fractional Whitham-Broer-Kaup Equations with Exponential-Decay Kernel**

**Humaira Yasmin**

Department of Basic Sciences, Preparatory Year Deanship, King Faisal University, Al-Ahsa 31982, Saudi Arabia<EMAIL_ADDRESS>
**Abstract:** This paper presents the semi-analytical analysis of the fractional-order non-linear coupled system of Whitham-Broer-Kaup equations. An iterative process is designed to analyze analytical findings to the specified non-linear partial fractional derivatives scheme utilizing the Yang transformation coupled with the Adomian technique. The fractional derivative is considered in the sense of Caputo-Fabrizio. Two numerical problems show the suggested method. Moreover, the results of the suggested technique are compared with the solution of other well-known numerical techniques such as the Homotopy perturbation technique, Adomian decomposition technique, and the Variation iteration technique. Numerical simulation has been carried out to verify that the suggested methodologies are accurate and reliable, and the results are revealed using graphs and tables. Comparing the analytical and actual solutions demonstrates that the proposed approaches effectively solve complicated non-linear problems. Furthermore, the proposed methodologies control and manipulate the achieved numerical solutions in a vast acceptable region in an extreme manner. It will provide us with a simple process to control and adjust the convergence region of the series solution.

**Keywords:** Adomian decomposition method; system of Whitham-Broer-Kaup equations; Caputo-Fabrizio derivative; Yang transform

### **1. Introduction**

Fractional calculus (FC) was invented by Newton, but it has recently piqued the interest of many academics. Fascinating breakthroughs in science and engineering applications have been found within the framework of FC over the last 30 years. Due to the complications involved with a heterogeneity issue, the notion of the fractional derivative has been industrialized. The behaviour of complex media with a diffusion mechanism may be captured using non-integer order differential operators [1–4]. It has proven a handy tool, and differential equations of any order may demonstrate various situations more efficiently and precisely. Numerous scholars began to work on calculus and its generalization to express their viewpoints while investigating many complicated events due to the rapid development of mathematical approaches using computer software [5–8].

Differential equations featuring non-linearities are used in science, technology, and engineering to explain a variety of phenomena, ranging from gravity to dynamical systems [9–11]. Non-linear partial differential equations (PDEs) are significant techniques for modeling non-linear dynamical events in a variety of domains, including mathematical biology, fluid mechanics, material science, and fluid dynamics, as shown in [12]. A sufficient set of partial differential equations can represent the bulk of dynamical systems. PDEs are also well-known for being utilized to solve mathematical difficulties like the Poincare and the Calabi conjectures.

It has already been demonstrated that the non-linear development of shallow-water waves may be represented using the technique of the Whitham-Broer-Kaup equation in fluid mechanics (WBKEs) [13]. Whitham, Broer, and Kaup [14–16] developed the integrated framework of the equations as mentioned earlier. The mentioned equations can be written

**Citation:** Yasmin, H. Numerical Analysis of Time-Fractional Whitham-Broer-Kaup Equations with Exponential-Decay Kernel. *Fractal Fract.* **2022**, *6*, 142. https://doi.org/ 10.3390/fractalfract6030142

Academic Editors: Asifa Tassaddiq and Muhammad Yaseen

Received: 10 February 2022 Accepted: 27 February 2022 Published: 2 March 2022

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

**Copyright:** © 2022 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 (https:// creativecommons.org/licenses/by/ 4.0/).

as the shallow water acoustic waves with various diversity connections, as shown [17]. In the classical order, the governing equations for the phenomenon mentioned above are represented by

$$\begin{cases} \mathbb{U}\_3 + \mathbb{U}\mathbb{U}\_\ell + \mathbb{V}\_\ell + q\mathbb{U}\_{\ell\varepsilon} = 0\\ \mathbb{V}\_3 + \mathbb{V}\mathbb{U}\_\ell + \mathbb{U}\mathbb{V}\_\ell - q\mathbb{V}\_{\ell\varepsilon} + p\mathbb{U}\_{\ell\varepsilon\varepsilon} = 0, \end{cases} \tag{1}$$

where U = U(*ε*, ), V = V(*ε*, ) indicates the horizontal velocity and height of the fluids, respectively, which differ greatly from the equilibrium, and *q*, *p* are the constants that are composed of various diffusion powers. For the past few decades, investigating the results of non-linear PDEs has been a major focus of research. Several authors have devised numerous mathematical methods to examine approximate results of non-linear PDEs. Mohyud Din et al. [18] investigated the analysis of many integer order PDEs using homotopy perturbation techniques. To solve the coupled set of Burgers and Brusselator equations, Biazar and Aminikhah [19] used the perturbation technique. For the numerical result of many traditional order differential equations by applying other techniques, interested readers can refer to Refs. [20–25]. Numerous strategies have been used to study the solution to the given non-linear coupled scheme (1) of PDEs. To address the classical order coupled systems of the WBK problem, Mohyud-Din et al. [26] employed perturbation methods. As a result, researchers like Xie et al. 2002 (who studied the solution using the hyperbolic technique) have used several powerful and efficient methods to investigate the problem of the WBK coupled equation of classical order PDEs. In the same way, El-Sayed and Kaya used the Adomian decomposition approach to investigate the scheme (1). Moreover, Ahmad et al. [12] used the Adomian decomposition method and He's polynomial to solve the coupled system (1).

Adomian proposed the Adomian decomposition method in 1980, which is a helpful technique for obtaining an explicit and numerical solution to a system of differential equations that represents a physical problem [27–29]. The Laplace transform technique is a vital technique in technology and applied mathematics. Combining the Adomian decomposition method and Yang transformation leads to a well-known technique named the Yang decomposition method. In this study, we convert differential equations to algebraic equations using the Laplace transform, and the non-linear terms are decomposed using Adomain polynomials. This numerical approach is effective for both deterministic and stochastic differential equation systems. It can be applied to a classical and fractional-order ordinary and a PDEs system, both linear and non-linear. There is no need for perturbation or liberalization in this procedure. Furthermore, unlike RK4, it does not require a pre-defined step size. In addition, this technique does not depend upon a parameter, as required for homotopy analysis and homotopy perturbation methods. However, the solutions achieved via this technique are the same as gained by the Adomian decomposition method (for detail, see [30–33]). It must be mentioned that the Yang decomposition method is more effective than the basic Adomian decomposition method.

The rest of this article is organized as follows. In Section 2, we present some basic definitions and properties. Section 3 describes the Yang decomposition method for solving fractional partial differential equations. The conclusion is presented at the end of the article.

### **2. Preliminaries Concepts**

In this section, we provide the fundamental definitions that will be used throughout the article. For the purpose of simplification, we write the exponential decay kernel as, *K*(, ) = *e*[−℘(−/1−℘)].

**Definition 1.** *If the Caputo-Fabrizio derivative is given as follows [34]:*

$${}^{CF}D\_{\odot}^{\wp}[\mathbb{P}(\heartsuit)] = \frac{N(\wp)}{1-\wp} \int\_{0}^{\heartsuit} \mathbb{P}'(\varrho)K(\mathbb{S},\varrho)d\varrho, \ n-1 < \wp \le n \tag{2}$$

*N*(℘) *is the normalization function with N*(0) = *N*(1) = 1*.*

$${}^{CF}D\_{\odot}^{\phi}[\mathbb{P}(\mathbb{S})] = \frac{N(\wp)}{1-\wp} \int\_{0}^{\mathbb{S}} [\mathbb{P}(\mathbb{S}) - \mathbb{P}(\varrho)] K(\mathbb{S}, \varrho) d\varrho. \tag{3}$$

**Definition 2.** *The fractional integral Caputo-Fabrizio is given as [34]*

$${}^{CF}I\_3^{\wp}[\mathbb{P}(\mathfrak{D})] = \frac{1-\wp}{N(\wp)}\mathbb{P}(\mathfrak{D}) + \frac{\wp}{N(\wp)}\int\_0^{\mathfrak{D}} \mathbb{P}(\varrho)d\varrho,\ \ \mathfrak{D} \ge 0, \wp \in (0,1]. \tag{4}$$

**Definition 3.** *For N*(℘) = 1*, the following result shows the Caputo-Fabrizio derivative of Laplace transformation [34]:*

$$L\left[ {}^{CF}D\_{\odot}^{\diamond}[\mathbb{P}(\mathbb{S})] \right] = \frac{\upsilon L\left[ \mathbb{P}(\mathbb{S}) - \mathbb{P}(0) \right]}{\upsilon + \wp(1 - \upsilon)}.\tag{5}$$

**Definition 4.** *The Yang transformation of* P() *is expressed as [35]*

$$\mathbb{Y}[\mathbb{P}(\mathbb{S})] = \chi(\upsilon) = \int\_0^\infty \mathbb{P}(\mathbb{S}) e^{-\frac{\mathbb{S}}{\upsilon}} d\mathbb{S}. \ \mathbb{S} > 0,\tag{6}$$

**Remark 1.** *The Yang transformation of a few useful functions is defined as:*

$$\begin{aligned} \mathbb{Y}[1] &= \upsilon, \\ \mathbb{Y}[3] &= \upsilon^2, \\ \mathbb{Y}[3^i] &= \Gamma(i+1)\upsilon^{i+1}. \end{aligned} \tag{7}$$

**Lemma 1.** *Let the Laplace transformation of* P() *is F*(*v*)*, then χ*(*v*) = *F*(1/*v*) *[36].*

**Proof.** From Equation (6), we can achieve another type of the Yang transformation by putting /*v* = *ζ* as

$$L[\mathbb{P}(\mathbb{S})] = \chi(v) = v \int\_0^\infty \mathbb{P}(v\zeta) e^{\zeta} d\zeta. \ \zeta > 0,\tag{8}$$

Since *L*[P()] = *F*(*v*), this implies that

$$F(\upsilon) = L[\mathbb{P}(\mathbb{S})] = \int\_0^\infty \mathbb{P}(\mathbb{S}) e^{-\upsilon \cdot \mathbb{S}} d\mathbb{S}.\tag{9}$$

Put = *ζ*/*v* in (9), we have

$$F(v) = \frac{1}{v} \int\_0^\infty \mathbb{P}\left(\frac{\zeta}{v}\right) e^{\zeta} d\zeta. \tag{10}$$

Thus, from Equation (8), we achieve

$$F(v) = \chi\left(\frac{1}{v}\right). \tag{11}$$

Additionally, from Equations (6) and (9), we achieve

$$F\left(\frac{1}{\upsilon}\right) = \chi(\upsilon). \tag{12}$$

The connections Equations between (11) and (12) represent the duality link between the Laplace and Yang transformation.

**Lemma 2.** *Let* P() *be a continuous function; then, the Caputo-Fabrizio derivative Yang transformation of* P() *is defined by [36]*

$$\mathbb{Y}[\mathbb{P}(\mathbb{S})] = \frac{\mathbb{Y}[\mathbb{P}(\mathbb{S}) - \upsilon \mathbb{P}(0)]}{1 + \wp(\upsilon - 1)}.\tag{13}$$

**Proof.** The Caputo-Fabrizio fractional Laplace transformation is given by

$$L[\mathbb{P}(\mathfrak{P})] = \frac{L[\mathfrak{P}(\mathfrak{P}) - \mathbb{P}(0)]}{v + \wp(1 - v)},\tag{14}$$

In addition, we have that the connection among Laplace and Yang property, i.e., *χ*(*v*) = *F*(1/*v*). To achieve the necessary result, we substitute *v* by 1/*v* in Equation (14), and get

$$\begin{split} \mathbb{V}[\mathbb{P}(\mathcal{S})] &= \frac{\frac{1}{\upsilon} \mathbb{V}[\mathbb{P}(\mathcal{S}) - \mathbb{P}(0)]}{\frac{1}{\upsilon} + \wp(1 - \frac{1}{\upsilon})}, \\ \mathbb{V}[\mathbb{P}(\mathcal{S})] &= \frac{\mathbb{V}[\mathbb{P}(\mathcal{S}) - \upsilon \mathbb{P}(0)]}{1 + \wp(\upsilon - 1)}. \end{split} \tag{15}$$

The proof is completed.

### **3. The Producer of YDM**

In this portion, we discuses the YDM producer for fractional partial differential equations.

$$\begin{aligned} \, \, ^{CF}D\_{\, 3}^{\phi} \mathbf{U}(\varepsilon, \mathbb{S}) + \mathcal{G}\_{1}(\mathbf{U}, \mathbb{V}) + \mathcal{L}\_{1}(\mathbb{U}, \mathbb{V}) - \mathcal{P}\_{1}(\varepsilon, \mathbb{S}) &= 0, \\ \, ^{CF}D\_{\, 3}^{\phi} \mathbf{V}(\varepsilon, \mathbb{S}) + \mathcal{G}\_{2}(\mathbb{U}, \mathbb{V}) + \mathcal{L}\_{2}(\mathbb{U}, \mathbb{V}) - \mathcal{P}\_{2}(\varepsilon, \mathbb{S}) &= 0, \quad 0 < \wp \le 1, \end{aligned} \tag{16}$$

with initial condition

$$\mathbb{U}(\varepsilon,0) = \mathbb{g}\_1(\varepsilon), \quad \mathbb{V}(\varepsilon,0) = \mathbb{g}\_2(\varepsilon). \tag{17}$$

where *D*℘ = *∂*℘ *∂*℘ is the Caputo fractional derivative of order ℘, G1, G2 and L1, L2 are the linear and non-linear functions, respectively, and P1,P2 are the source functions. Using the Yang transformation to Equation (16),

$$\begin{aligned} \mathbb{V}[D\_3^{\diamond}\mathbb{U}(\varepsilon, \mathbb{V})] + \mathbb{V}[\mathcal{G}\_1(\mathbb{U}, \mathbb{V}) + \mathcal{L}\_1(\mathbb{U}, \mathbb{V}) - \mathcal{P}\_1(\varepsilon, \mathbb{V})] &= 0, \\ \mathbb{V}[D\_3^{\diamond}\mathbb{V}(\varepsilon, \mathbb{V})] + \mathbb{V}[\mathcal{G}\_2(\mathbb{U}, \mathbb{V}) + \mathcal{L}\_2(\mathbb{U}, \mathbb{V}) - \mathcal{P}\_2(\varepsilon, \mathbb{V})] &= 0. \end{aligned} \tag{18}$$

Using the Yang transformation differentiation property, we have

$$\mathbb{V}[\mathbb{U}(\varepsilon,\mathbb{S})] = \nu \mathbb{U}(\varepsilon,0) + (1+\wp(\upsilon-1))\mathbb{Y}[\mathcal{P}\_1(\varepsilon,\mathbb{S})] - (1+\wp(\upsilon-1))\mathbb{Y}[\mathcal{G}\_1(\mathbb{U},\mathbb{V}) + \mathcal{L}\_1(\mathbb{U},\mathbb{V})],$$

$$\mathbb{V}[\mathbb{V}(\varepsilon,\mathbb{S})] = \nu \mathbb{V}(\varepsilon,0) + (1+\wp(\upsilon-1))\mathbb{Y}[\mathcal{P}\_2(\varepsilon,\mathbb{S})] - (1+\wp(\upsilon-1))\mathbb{Y}\{\mathcal{G}\_2(\mathbb{U},\mathbb{V}) + \mathcal{L}\_2(\mathbb{U},\mathbb{V})\}. \tag{19}$$

YDM describes the solution of infinite series U(*ε*, ) and V(*ε*, ),

$$\mathbb{U}(\mathfrak{e}, \mathfrak{Z}) = \sum\_{m=0}^{\infty} \mathbb{U}\_{m}(\mathfrak{e}, \mathfrak{Z}), \quad \mathbb{V}(\mathfrak{e}, \mathfrak{Z}) = \sum\_{m=0}^{\infty} \mathbb{V}\_{m}(\mathfrak{e}, \mathfrak{Z}). \tag{20}$$

Adomian polynomials of non-linear terms of L1 and L2 are represented as

$$\mathcal{L}\_1(\mathbb{U}, \mathbb{V}) = \sum\_{m=0}^{\infty} \mathcal{A}\_{\mathfrak{m}\_{\prime}} \quad \mathcal{L}\_2(\mathbb{U}, \mathbb{V}) = \sum\_{m=0}^{\infty} \mathcal{B}\_{\mathfrak{m}}.\tag{21}$$

The expression for Adomian polynomials is

$$\begin{split} \mathcal{A}\_{\text{ill}} &= \frac{1}{m!} \left[ \frac{\partial^{\text{m}}}{\partial \lambda^{\text{m}}} \left\{ \sum\_{m=0}^{\infty} \lambda^{\text{m}} \mathbb{U}\_{\text{m}} \sum\_{m=0}^{\infty} \lambda^{\text{m}} \mathbb{V}\_{\text{m}} \right\} \right]\_{\lambda=0} \\ \mathcal{B}\_{\text{ill}} &= \frac{1}{m!} \left[ \frac{\partial^{\text{m}}}{\partial \lambda^{\text{m}}} \left\{ \sum\_{m=0}^{\infty} \lambda^{\text{m}} \mathbb{U}\_{\text{m}} \sum\_{m=0}^{\infty} \lambda^{\text{m}} \mathbb{V}\_{\text{m}} \right\} \right]\_{\lambda=0} . \end{split} \tag{22}$$

Putting Equations (20) and (22) into Equation (19),

$$\begin{split} \mathbb{Y}\left[\sum\_{m=1}^{\infty} \mathbb{U}\_{m}(\varepsilon, \mathbb{S})\right] &= \ \_{s} \mathbb{U}(\varepsilon, 0) + (1 + \wp(\upsilon - 1)) \mathbb{Y}\{\mathcal{P}\_{1}(\varepsilon, \mathbb{S})\} \\ &- (1 + \wp(\upsilon - 1)) \mathbb{Y}\left\{\mathcal{G}\_{1}(\sum\_{m=0}^{\infty} \mathbb{U}\_{m}, \sum\_{m=0}^{\infty} \mathbb{V}\_{m}) + \sum\_{m=0}^{\infty} \mathcal{A}\_{m} \right\}, \\ \mathbb{Y}[\sum\_{m=1}^{\infty} \mathbb{V}\_{m}(\varepsilon, \mathbb{S})] &= \ \_{\upsilon} \mathbb{V}(\varepsilon, 0) + (1 + \wp(\upsilon - 1)) \mathbb{Y}\{\mathcal{P}\_{2}(\varepsilon, \mathbb{S})\} \\ &- (1 + \wp(\upsilon - 1)) \mathbb{Y}\left\{\mathcal{G}\_{2}(\sum\_{m=0}^{\infty} \mathbb{U}\_{m}, \sum\_{m=0}^{\infty} \mathbb{V}\_{m}) + \sum\_{m=0}^{\infty} \mathcal{B}\_{m} \right\}. \end{split} \tag{23}$$

The inverse Yang transformation is implemented on Equation (23),

$$\begin{aligned} \sum\_{m=1}^{\infty} \mathsf{U}\_{m}(\mathfrak{e}, \mathfrak{I}) &= \mathbb{Y}^{-1} [\upsilon \mathsf{U}(\mathfrak{e}, \mathfrak{I}) + (1 + \wp(\mathfrak{v} - 1)) \mathsf{Y} \{\mathscr{P}\_{\mathfrak{I}}(\mathfrak{e}, \mathfrak{I})\}] \\ &- \mathsf{Y}^{-1} \left[ (1 + \wp(\mathfrak{v} - 1)) \mathsf{Y} \left\{ \mathscr{G}\_{1} (\sum\_{m=0}^{\infty} \mathsf{U}\_{m}, \sum\_{m=0}^{\infty} \mathsf{V}\_{m}) + \sum\_{m=0}^{\infty} \mathscr{A}\_{m} \right\} \right], \\ \sum\_{m=1}^{\infty} \mathsf{V}\_{m}(\mathfrak{e}, \mathfrak{I}) &= \mathsf{Y}^{-1} [\upsilon \mathsf{V}(\mathfrak{e}, \mathfrak{I}) + (1 + \wp(\mathfrak{v} - 1)) \mathsf{Y} \{\mathscr{P}\_{\mathfrak{Z}}(\mathfrak{e}, \mathfrak{I})\}] \\ &- \mathsf{Y}^{-1} \left[ (1 + \wp(\mathfrak{v} - 1)) \mathsf{Y} \left\{ \mathscr{G}\_{2} (\sum\_{m=0}^{\infty} \mathsf{U}\_{m}, \sum\_{m=0}^{\infty} \mathsf{V}\_{m}) + \sum\_{m=0}^{\infty} \mathscr{B}\_{m} \right\} \right]. \end{aligned} \tag{24}$$

Find the U0 and V0 using the initial conditions and sources functions. The following terms are expressed:

$$\begin{aligned} \mathbb{U}\_{0}(\varepsilon,\mathbb{S}) &= \mathbb{Y}^{-1}[\mathbb{v}\mathbb{U}(\varepsilon,0) + (1+\wp(\upsilon-1))\mathbb{Y}\{\mathcal{P}\_{1}(\varepsilon,\mathbb{S})\}],\\ \mathbb{V}\_{0}(\varepsilon,\mathbb{S}) &= \mathbb{Y}^{-1}[\mathbb{v}\mathbb{V}(\varepsilon,0) + (1+\wp(\upsilon-1))\mathbb{Y}\{\mathcal{P}\_{2}(\varepsilon,\mathbb{S})\}].\end{aligned} \tag{25}$$

For *m* = 1

$$\begin{aligned} \mathcal{U}\_1(\mathfrak{e}, \mathfrak{z}) &= -\mathbb{Y}^{-1} [(1 + \wp(\upsilon - 1)) \mathbb{Y} \{ \mathcal{G}\_1(\mathbb{U}\_0, \mathbb{V}\_0) + \mathcal{A}\_0 \} ], \\ \mathcal{V}\_1(\mathfrak{e}, \mathfrak{z}) &= -\mathbb{Y}^{-1} [(1 + \wp(\upsilon - 1)) \mathbb{Y} \{ \mathcal{G}\_2(\mathbb{U}\_0, \mathbb{V}\_0) + \mathcal{B}\_0 \} ], \end{aligned}$$

the general for *m* ≥ 1, is given by

$$\begin{aligned} \mathbb{U}\_{m+1}(\mathfrak{e}, \mathfrak{d}) &= -\mathbb{Y}^{-1}[(1 + \wp(\upsilon - 1))\mathbb{Y}\{\mathcal{G}\_1(\mathbb{U}\_m, \mathbb{V}\_m) + \mathcal{A}\_m\}], \\ \mathbb{V}\_{m+1}(\mathfrak{e}, \mathfrak{d}) &= -\mathbb{Y}^{-1}[(1 + \wp(\upsilon - 1))\mathbb{Y}\{\mathcal{G}\_2(\mathbb{U}\_m, \mathbb{V}\_m) + \mathcal{B}\_m\}], \end{aligned}$$

### **4. Numerical Results**

**Example 1.** *Consider the fractional-order system of WBKEs [11]*

$$\begin{split} & \mathbb{E}^{\rm CF} D\_{\odot}^{\rm c} \mathbb{U}(\varepsilon, \mathfrak{z}) + \mathbb{U}(\varepsilon, \mathfrak{z}) \frac{\partial \mathbb{U}(\varepsilon, \mathfrak{z})}{\partial \mathbb{U}(\varepsilon, \mathfrak{z})} + \mathbb{V}(\varepsilon, \mathfrak{z}) \frac{\partial \mathbb{U}(\varepsilon, \mathfrak{z})}{\partial \mathbb{U}(\varepsilon, \mathfrak{z})} + 3 \frac{\partial^{3} \mathbb{U}(\varepsilon, \mathfrak{z})}{\partial \varepsilon^{3}} - \frac{\partial^{2} \mathbb{U}(\varepsilon, \mathfrak{z})}{\partial \varepsilon^{2}} = 0, \\ & 0 < \wp \le 1, \quad -1 < 3 \le 1, \quad -10 \le \varepsilon \le 10, \end{split} \tag{26}$$

*with the initial conditions*

$$\begin{aligned} \mathcal{U}(\varepsilon, 0) &= \frac{1}{2} - 8 \tanh(-2\varepsilon), \\ \mathcal{V}(\varepsilon, 0) &= 16 - 16 \tanh^2(-2\varepsilon). \end{aligned} \tag{27}$$

*Applying the Yang transformation of Equation (26), we have*

$$\begin{split} \mathbb{V}\left[\frac{\partial\mathcal{J}(\mathbf{c},\mathcal{J})}{\partial\mathcal{J}}\right] &= \mathbb{V}\left[\mathbb{I}(\mathbf{c},\mathcal{J})\mathbb{I}\right] - \mathbb{V}\left[\mathbb{I}(\mathbf{c},\mathcal{J})\mathbb{I}\right] - \mathbb{V}\left[\mathbb{I}(\mathbf{c},\mathcal{J})\mathbb{I}\right] \\ \mathbb{V}\left[\frac{\partial\mathcal{J}^{3}}{\partial\mathcal{J}^{3}}\right] &= -\mathbb{V}\left[\mathbb{I}(\mathbf{c},\mathcal{J})\frac{\partial\mathbb{I}}{\partial\mathcal{J}} + \mathbb{V}(\mathbf{c},\mathcal{J})\frac{\partial\mathbb{I}}{\partial\mathcal{J}}(\mathcal{E},\mathcal{J})\right] - \mathbb{V}\left[\mathbb{I}(\mathbf{c},\mathcal{J})\mathbb{I}\right] \\ \mathbb{V}\left[\frac{\partial\mathcal{J}^{3}}{\partial\mathcal{J}}\mathbb{I}(\mathbf{c},\mathcal{J})\right] &= -\mathbb{V}\left[\mathbb{I}(\mathbf{c},\mathcal{J})\frac{\partial\mathbb{I}}{\partial\mathcal{J}}(\mathcal{E},\mathcal{J}) + \frac{\partial\mathbb{I}}{\partial\mathcal{J}}(\mathcal{E},\mathcal{J})\frac{\partial\mathbb{I}}{\partial\mathcal{J}} + 3\frac{\partial\mathbb{I}}{\partial\mathcal{J}}(\mathcal{E},\mathcal{J})\right], \\ \mathbb{V}\left[\frac{1}{\partial\mathcal{J}(\mathcal{E},\mathcal{J})}\mathbb{I}(\mathcal{J}(\mathcal{E},\mathcal{J})) - \mathbb{V}(\mathcal{E},\mathcal{J}) &= -\mathbb{V}\left[\mathbb{I}(\mathbf{c},\mathcal{J})\frac{\partial\mathbb{I}}{\partial\mathcal{J}}(\mathcal{E},\mathcal{J}) + \frac{\partial\mathbb{I}}{\partial\mathcal{J}}(\mathcal{E$$

*The above equation is simplified*

$$\begin{split} \mathbb{Y}\{\mathbb{U}(\varepsilon,\mathbb{O})\} &= \mathbb{z}\{\mathbb{U}(\varepsilon,\mathbb{O})\} - (1+\wp(\upsilon-1))\mathbb{Y}\left[\mathbb{U}(\varepsilon,\mathbb{O})\frac{\partial\mathbb{U}(\varepsilon,\mathbb{O})}{\partial\varepsilon} + \frac{\varrho\varepsilon}{\mathfrak{d}\mathbb{I}(\varepsilon,\mathbb{O})} + \frac{\mathfrak{g}\varepsilon}{\mathfrak{d}\mathbb{I}(\varepsilon,\mathbb{O})} + \frac{\mathfrak{g}\varepsilon}{\mathfrak{d}\mathbb{I}(\varepsilon,\mathbb{O})}\right], \\ \mathbb{Y}\{\mathbb{V}(\varepsilon,\mathbb{O})\} &= \mathbb{z}\{\mathbb{V}(\varepsilon,\mathbb{O})\} - (1+\wp(\upsilon-1))\mathbb{Y}\left[\mathbb{U}(\varepsilon,\mathbb{O})\frac{\partial\mathbb{U}(\varepsilon,\mathbb{O})}{\partial\varepsilon} + \mathbb{V}(\varepsilon,\mathbb{O})\frac{\partial\mathbb{U}(\varepsilon,\mathbb{O})}{\partial\varepsilon} + 3\frac{\mathfrak{g}^{2}\mathbb{U}(\varepsilon,\mathbb{O})}{\partial\varepsilon^{2}} - \frac{\mathfrak{g}^{2}\mathbb{U}(\varepsilon,\mathbb{O})}{\partial\varepsilon^{2}}\right]. \end{split} \tag{28}$$

### *Using inverse Yang transform, we have*

$$\begin{split} \mathcal{U}(\varepsilon,\mathcal{O}) &= \mathcal{U}(\varepsilon,0) - \mathcal{Y}^{-1} \Big[ (1+\wp(\upsilon-1))\mathcal{Y} \Big\{ \mathcal{U}(\varepsilon,\mathcal{O}) \frac{\partial \mathcal{U}(\varepsilon,\mathcal{O})}{\partial \varepsilon} + \frac{\partial \mathcal{U}(\varepsilon,\mathcal{O})}{\partial \mathcal{E}} + \frac{\partial \mathcal{U}(\varepsilon,\mathcal{O})}{\partial \mathcal{V}} \Big\} \Big], \\ \mathcal{V}(\varepsilon,\mathcal{O}) &= \mathcal{V}(\varepsilon,0) - \mathcal{Y}^{-1} \Big[ (1+\wp(\upsilon-1))\mathcal{Y} \Big\{ \mathcal{U}(\varepsilon,\mathcal{O}) \frac{\partial \mathcal{U}(\varepsilon,\mathcal{O})}{\partial \varepsilon} + \mathcal{Y}(\varepsilon,\mathcal{O}) \frac{\partial \mathcal{U}(\varepsilon,\mathcal{O})}{\partial \varepsilon} + 3\frac{\partial^{2} \mathcal{U}(\varepsilon,\mathcal{O})}{\partial \varepsilon^{3}} - \frac{\partial^{2} \mathcal{V}(\varepsilon,\mathcal{O})}{\partial \varepsilon^{2}} \Big\} \Big]. \end{split} \tag{29}$$

*Assume that the* U(*ε*, ) *and the* V(*ε*, ) *infinite series solution functions as follows:*

$$\mathbb{U}(\mathfrak{e}, \mathfrak{Z}) = \sum\_{m=0}^{\infty} \mathbb{U}\_{\mathfrak{m}}(\mathfrak{e}, \mathfrak{Z}) \quad \text{and} \quad \mathbb{V}(\mathfrak{e}, \mathfrak{Z}) = \sum\_{m=0}^{\infty} \mathbb{V}\_{\mathfrak{m}}(\mathfrak{e}, \mathfrak{Z}).$$

*Remember that the Adomian polynomials are given as* UU*ε* = ∑∞ *m*=0 A*m,* UV*ε* = ∑∞ *m*=0 B*m and* VU*ε* = ∑∞ *m*=0 C*m*

$$\begin{split} &\sum\_{m=0}^{\infty} \mathsf{U}\_{\mathsf{m}}(\varepsilon, \mathbb{S}) = \mathsf{U}(\varepsilon, 0) - \mathbb{Y}^{-1} \left[ (1 + \wp(\upsilon - 1)) \mathbb{Y} \left\{ \sum\_{m=0}^{\infty} \mathcal{A}\_{m} + \frac{\partial \mathsf{U}(\varepsilon, \mathbb{S})}{\partial \varepsilon} + \frac{\partial \mathsf{V}(\varepsilon, \mathbb{S})}{\partial \varepsilon} \right\} \right], \\ &\sum\_{m=0}^{\infty} \mathsf{V}\_{\mathsf{m}}(\varepsilon, \mathbb{S}) = \mathsf{V}(\varepsilon, 0) - \mathbb{Y}^{-1} \left[ (1 + \wp(\upsilon - 1)) \mathbb{Y} \left\{ \sum\_{m=0}^{\infty} \mathcal{B}\_{m} + \sum\_{m=0}^{\infty} \mathcal{C}\_{m} + 3 \frac{\partial^{3} \mathsf{U}(\varepsilon, \mathbb{S})}{\partial \varepsilon^{3}} - \frac{\partial^{2} \mathsf{V}(\varepsilon, \mathbb{S})}{\partial \varepsilon^{2}} \right\} \right], \end{split}$$

$$\begin{split} \sum\_{m=0}^{\infty} \mathcal{U}\_{m}(\varepsilon, \mathbb{Q}) &= \frac{1}{2} - 8 \tanh(-2\varepsilon) - \mathbb{Y}^{-1} \Big[ (1 + \wp(v-1)) \mathbb{Y} \left\{ \sum\_{m=0}^{\infty} \mathcal{A}\_{m} + \frac{\partial \mathcal{U}(\varepsilon, \mathbb{Q})}{\partial \varepsilon} + \frac{\partial \mathcal{V}(\varepsilon, \mathbb{Q})}{\partial \varepsilon} \right\} \Big], \\ \sum\_{m=0}^{\infty} \mathcal{V}\_{m}(\varepsilon, \mathbb{Q}) &= 16 - 16 \tanh^{2}(-2\varepsilon) - \mathbb{Y}^{-1} \Big[ (1 + \wp(v-1)) \mathbb{Y} \left\{ \sum\_{m=0}^{\infty} \mathcal{B}\_{m} + \sum\_{m=0}^{\infty} \mathcal{C}\_{m} + 3 \frac{\partial^{3} \mathcal{U}(\varepsilon, \mathbb{Q})}{\partial \varepsilon^{3}} - \frac{\partial^{2} \mathcal{V}(\varepsilon, \mathbb{Q})}{\partial \varepsilon^{2}} \right\} \Big]. \end{split} \tag{30}$$

*With the aid of the Adomian polynomial, according to Equation (22), all forms of non-linear may be stated as*

$$\begin{split} \mathcal{A}\_{0} &= \mathbb{U}\_{0} \frac{\partial \mathbb{U}\_{0}}{\partial \varepsilon}, \quad \mathcal{A}\_{1} = \mathbb{U}\_{0} \frac{\partial \mathbb{U}\_{1}}{\partial \varepsilon} + \mathbb{U}\_{1} \frac{\partial \mathbb{U}\_{0}}{\partial \varepsilon}, \quad \mathcal{B}\_{0} = \mathbb{U}\_{0} \frac{\partial \mathbb{V}\_{0}}{\partial \beta}, \quad \mathcal{B}\_{1} = \mathbb{U}\_{0} \frac{\partial \mathbb{V}\_{1}}{\partial \beta} + \mathbb{U}\_{1} \frac{\partial \mathbb{V}\_{0}}{\partial \beta}, \\ \mathcal{C}\_{0} &= \mathbb{V}\_{0} \frac{\partial \mathbb{U}\_{0}}{\partial \varepsilon}, \quad \mathcal{C}\_{1} = \mathbb{V}\_{0} \frac{\partial \mathbb{U}\_{1}}{\partial \varepsilon} + \mathbb{V}\_{1} \frac{\partial \mathbb{U}\_{0}}{\partial \varepsilon}. \end{split}$$

*Therefore we can easily obtain*

$$\mathbb{U}\_0(\varepsilon, \mathbb{S}) = \frac{1}{2} - 8 \tanh(-2\varepsilon), \quad \mathbb{V}\_0(\varepsilon, \mathbb{S}) = 16 - 16 \tanh^2(-2\varepsilon).$$

*For m* = 0

$$\mathbb{U}\_{1}(\varepsilon,\mathbb{S}) = -8\sec h^{2} (-2\varepsilon) \{ 1 + \wp \mathbb{S} - \wp \}, \quad \mathbb{V}\_{1}(\varepsilon,\mathbb{S}) = -32\sec h^{2} (-2\varepsilon) \tanh(-2\varepsilon) \{ 1 + \wp \mathbb{S} - \wp \}.$$

$$\text{For } m = 1$$

$$\begin{split} \mathrm{U}\_{2}(\varepsilon,\mathbb{S}) &= -16\sec h^{2} (-2\varepsilon) \Big\{ 4\sec h^{2} (-2\varepsilon) - 8\tanh^{2} (-2\varepsilon) + 3\tanh(-2\varepsilon) \Big\} \Big\{ (1-\wp)2\wp \mathbb{S} + (1-\wp)^{2} + \frac{\wp^{2}\mathfrak{S}^{2}}{2} \Big\}, \\ \mathrm{V}\_{2}(\varepsilon,\mathbb{S}) &= -32\sec h^{2} (-2\varepsilon) \Big\{ 40\sec h^{2} (-2\varepsilon) \tanh(-2\varepsilon) + 96\tanh(-2\varepsilon) - 2\tanh^{2} (-2\varepsilon) \\ &- 32\tanh^{3} (-2\varepsilon) - 25\sec h^{2} (-2\varepsilon) \Big\} \Big\{ (1-\wp)2\wp \mathbb{S} + (1-\wp)^{2} + \frac{\wp^{2}\mathfrak{S}^{2}}{2} \Big\}. \end{split}$$

*The remaining steps of the YDM outcomes may be conveniently gathered from* U*m and* V*m* (*m* ≥ 2) *using the same methods. Then, we assess the sequence of possibilities as follows:*

$$\mathbb{U}(\varepsilon,\mathbb{G}) = \sum\_{m=0}^{\infty} \mathbb{U}\_{m}(\varepsilon,\mathbb{G}) = \mathbb{U}\_{0}(\varepsilon,\mathbb{G}) + \mathbb{U}\_{1}(\varepsilon,\mathbb{G}) + \mathbb{U}\_{2}(\varepsilon,\mathbb{G}) + \mathbb{U}\_{3}(\varepsilon,\mathbb{G}) + \cdots \cdot \mathbb{I}\_{n}$$

$$\mathbb{V}(\varepsilon,\mathbb{G}) = \sum\_{m=0}^{\infty} \mathbb{V}\_{m}(\varepsilon,\mathbb{G}) = \mathbb{V}\_{0}(\varepsilon,\mathbb{G}) + \mathbb{V}\_{1}(\varepsilon,\mathbb{G}) + \mathbb{V}\_{2}(\varepsilon,\mathbb{G}) + \mathbb{V}\_{3}(\varepsilon,\mathbb{G}) + \cdots \cdot \mathbb{I}\_{n}$$

$$\mathbb{U}(\varepsilon,\mathbb{G}) = \frac{1}{2} - 8\tanh(-2\varepsilon) - 8\sec h^{2}(-2\varepsilon)\left\{1 + \wp\mathbb{3} - \wp\right\} - 16\sec h^{2}(-2\varepsilon)\left\{4\sec h^{2}(-2\varepsilon) + \cdots + \mathbb{I}\_{n}\right\}$$

$$-8\tanh^{2}(-2\varepsilon) + 3\tanh(-2\varepsilon)\left\{(1-\wp)2\wp\mathbb{3} + (1-\wp)^{2} + \frac{\wp^{2}\mathfrak{3}^{2}}{2}\right\} - \cdots \cdot \frac{1}{2}$$

$$\begin{split} \mathbb{V}(\varepsilon,\mathbb{S}) &= \mathbb{V}(\varepsilon,\mathbb{S}) + \mathbb{V}(\varepsilon,\mathbb{S}) \left\{ (\mathbb{V}(\varepsilon,\mathbb{S})\mathbb{I}^2)^2 + (\mathbb{V}(\varepsilon,\mathbb{S}))^2 \right\} \\ \mathbb{V}(\varepsilon,\mathbb{S}) &= 16 - 16 \tanh^2(-2\varepsilon) - 32 \sec h^2(-2\varepsilon) \tanh(-2\varepsilon) \left\{ 1 + \wp\mathbb{S} - \wp \right\} - 32 \sec h^2(-2\varepsilon) \\ &\quad \left\{ 40 \sec h^2(-2\varepsilon) \tanh(-2\varepsilon) + 96 \tanh(-2\varepsilon) - 2 \tanh^2(-2\varepsilon) - 32 \tanh^3(-2\varepsilon) \right\} \\ &\quad - 25 \sec h^2(-2\varepsilon) \left\{ (1-\wp)2\wp\mathbb{S} + (1-\wp)^2 + \frac{\wp^2 3^2}{2} \right\} - \dotsb . \end{split}$$

*At integer order* ℘ = 1*, the following series form of solution is achieved:*

$$\begin{split} \mathrm{U}(\varepsilon,\mathbb{S}) &= \frac{1}{2} - 8\tanh(-2\varepsilon) - 8\operatorname{sech}(-2\varepsilon)^{2}\mathbb{S} + 8\operatorname{sech}^{2}(-2\varepsilon) \\ &\times \left\{ 3\tanh(-2\varepsilon) + 8\tanh(-2\varepsilon)^{2} + 4\operatorname{sech}^{2}(-2\varepsilon) \right\} \mathbb{S}^{2} + \cdots \ . \\ \mathrm{V}(\varepsilon,\mathbb{S}) &= 16 - 16\tanh^{2}(-2\varepsilon) - 32\operatorname{sech}^{2}(-2\varepsilon)\tanh(-2\varepsilon) \mathbb{S} \\ &\quad - 16\sec h^{2}(-2\varepsilon) \left\{ 96\tanh(-2\varepsilon) - 32\tanh^{3}(-2\varepsilon) \right. \\ &\left. + 40\sec h^{2}(-2\varepsilon)\tanh(-2\varepsilon) - 2\tanh^{2}(-2\varepsilon) - 25\sec h^{2}(-2\varepsilon) \right\} \mathbb{S}^{2} + \cdots \ . \end{split}$$

*The exact solution of Equation (26) at* ℘ = 1*,*

$$\begin{split} \mathcal{U}(\varepsilon, \mathbb{S}) &= \frac{1}{2} - 8 \tanh \left\{ -2 \left( \varepsilon - \frac{\mathbb{S}}{2} \right) \right\}, \\ \mathcal{V}(\varepsilon, \mathbb{S}) &= 16 - 16 \tanh^2 \left\{ -2 \left( \varepsilon - \frac{\mathbb{S}}{2} \right) \right\}. \end{split} \tag{31}$$

*In Figures 1 and 2, the actual and Yang decomposition method solutions at an integerorder* ℘ = 1 *are represented for both* U(*ε*, ) *and* V(*ε*, ) *of Example 1. It is observed that Yang decomposition method results are in good contact with the actual result of the models. In Figures 3 and 4, various fractional-order solutions of Example 1, at different fractional-orders,* ℘ = 1, 0.8, 0.6, 0.4 *are plotted. It is investigated that for Example 1, the fractional-order solutions are convergent to an integer-order solution for both* U(*ε*, ) *and* V(*ε*, )*. In Tables 1 and 2 show that yang decomposition method of different fractional order* ℘ *of Example 1. In Tables 3 and 4 compassion of different analytical and numerical methods of Example 1.*

**Figure 1.** The actual and YDM solution of U(*ε*, ) at ℘ = 1 of Example 1.

**Figure 2.** The actual and YDM solution of V(*ε*, ) at ℘ = 1 of Example 1.

**Figure 3.** The fractional-order solutions of V(*ε*, ) at ℘ of Example 1.

**Figure 4.** The fractional-order solutions of V(*ε*, ) at ℘ of Example 1.

**Table 2.** YDM solution of V(*ε*, ) at various fractional-order ℘ of Example 1.

**Table 3.** Comparison of absolute error (AE) of U(*ε*, ) at ℘ = 1 obtained by various methods.

**Table 4.** Comparison of absolute error of V(*ε*, ) at ℘ = 1 obtained by various methods.

**Example 2.** *Consider the fractional-order system of WBKEs [11]*

$$\begin{split} & \; ^{CF}D\_{3}^{\beta} \mathbb{U}(\varepsilon, \mathfrak{I}) + \mathbb{U}(\varepsilon, \mathfrak{I}) \frac{\partial \mathbb{U}(\varepsilon, \mathfrak{I})}{\partial \varepsilon} + \frac{1}{2} \frac{\partial \mathbb{U}(\varepsilon, \mathfrak{I})}{\partial \varepsilon} + \frac{\partial \mathbb{U}(\varepsilon, \mathfrak{I})}{\partial \varepsilon} + \frac{\partial \mathbb{U}(\varepsilon, \mathfrak{I})}{\partial \varepsilon} = 0, \\ & \; ^{CF}D\_{3}^{\beta} \mathbb{V}(\varepsilon, \mathfrak{I}) + \mathbb{U}(\varepsilon, \mathfrak{I}) \frac{\partial \mathbb{U}(\varepsilon, \mathfrak{I})}{\partial \varepsilon} + \mathbb{U}(\varepsilon, \mathfrak{I}) \frac{\partial \mathbb{U}(\varepsilon, \mathfrak{I})}{\partial \varepsilon} - \frac{1}{2} \frac{\partial^{2} \mathbb{V}(\varepsilon, \mathfrak{I})}{\partial \varepsilon^{2}} = 0, \\ & \; 0 < \wp \le 1, \quad 0 < \mathfrak{I} \le 1, \quad -100 \le \varepsilon \le 100, \end{split} \tag{32}$$

*with the initial conditions*

$$\begin{aligned} \mathcal{U}(\varepsilon, 0) &= \xi - \kappa \coth[\kappa(\varepsilon + \theta)], \\ \mathcal{V}(\varepsilon, 0) &= -\kappa^2 \text{cosech}^2[\kappa(\varepsilon + \theta)]. \end{aligned} \tag{33}$$

*Applying the Yang transformation of Equation (32), we have*

Y *∂*℘U(*ε*, ) *∂*℘ = −Y U(*ε*, ) *∂*U(*ε*, ) *∂ε* + 1 2 *∂*U(*ε*, ) *∂ε* + *∂*V(*ε*, ) *∂ε* , Y *∂*℘V(*ε*, ) *∂*℘ = −Y U(*ε*, ) *∂*V(*ε*, ) *∂ε* + V(*ε*, ) *∂*U(*ε*, ) *∂ε* − 1 2 *∂*2V(*ε*, ) *∂ε*2 , 1 (1 + ℘(*v* − 1))Y{U(*ε*, )} − *v*U(*ε*, 0) = −Y U(*ε*, ) *∂*U(*ε*, ) *∂ε* + 1 2 *∂*U(*ε*, ) *∂ε* + *∂*V(*ε*, ) *∂ε* , 1 (1 + ℘(*v* − 1))Y{V(*ε*, )} − *v*V(*ε*, 0) = −Y U(*ε*, ) *∂*V(*ε*, ) *∂ε* + V(*ε*, ) *∂*U(*ε*, ) *∂ε* − 1 2 *∂*2V(*ε*, ) *∂ε*2 . *The above equation is simplified*

$$\begin{split} \mathbb{V}\{\mathbb{U}(\varepsilon,\mathbb{S})\} &= \mathbb{V}\{\mathbb{U}(\varepsilon,0)\} - (1+\wp(\upsilon-1))\mathbb{V}\left[\mathbb{U}(\varepsilon,\mathbb{S})\frac{\partial\mathbb{U}(\varepsilon,\mathbb{S})}{\partial\varepsilon} + \frac{1}{2}\frac{\partial\mathbb{U}(\varepsilon,\mathbb{S})}{\partial\varepsilon} + \frac{\partial\mathbb{V}(\varepsilon,\mathbb{S})}{\partial\varepsilon}\right], \\ \mathbb{V}\{\mathbb{V}(\varepsilon,\mathbb{S})\} &= \mathbb{v}\{\mathbb{V}(\varepsilon,0)\} - (1+\wp(\upsilon-1))\mathbb{V}\left[\mathbb{U}(\varepsilon,\mathbb{S})\frac{\partial\mathbb{V}(\varepsilon,\mathbb{S})}{\partial\varepsilon} + \mathbb{V}(\varepsilon,\mathbb{S})\frac{\partial\mathbb{U}(\varepsilon,\mathbb{S})}{\partial\varepsilon} - \frac{1}{2}\frac{\partial^{2}\mathbb{V}(\varepsilon,\mathbb{S})}{\partial\varepsilon^{2}}\right]. \end{split} \tag{34}$$

*Using inverse Yang transform, we have*

$$\begin{split} \mathcal{U}(\varepsilon,\mathbb{O}) &= \mathcal{U}(\varepsilon,0) - \mathbb{Y}^{-1} \Big[ (1+\wp(\upsilon-1)) \mathbb{Y} \Big\{ \mathcal{U}(\varepsilon,\mathbb{O}) \frac{\partial \mathcal{U}(\varepsilon,\mathbb{O})}{\partial \varepsilon} + \frac{1}{2} \frac{\partial \mathcal{U}(\varepsilon,\mathbb{O})}{\partial \varepsilon} + \frac{\partial \mathcal{U}(\varepsilon,\mathbb{O})}{\partial \varepsilon} \Big\} \Big], \\ \mathcal{V}(\varepsilon,\mathbb{O}) &= \mathcal{V}(\varepsilon,0) - \mathbb{Y}^{-1} \Big[ (1+\wp(\upsilon-1)) \mathbb{Y} \Big\{ \mathcal{U}(\varepsilon,\mathbb{O}) \frac{\partial \mathcal{V}(\varepsilon,\mathbb{O})}{\partial \varepsilon} + \mathbb{V}(\varepsilon,\mathbb{O}) \frac{\partial \mathcal{U}(\varepsilon,\mathbb{O})}{\partial \varepsilon} - \frac{1}{2} \frac{\partial^{2} \mathcal{V}(\varepsilon,\mathbb{O})}{\partial \varepsilon^{2}} \Big\} \Big]. \end{split} \tag{35}$$

*Assume that the infinite series solution functions* U(*ε*, ) *and* V(*ε*, ) *are as follows:*

$$\mathrm{U}(\mathfrak{e}, \mathfrak{Z}) = \sum\_{m=0}^{\infty} \mathrm{U}\_{\mathfrak{m}}(\mathfrak{e}, \mathfrak{Z}), \quad \text{and} \quad \mathrm{V}(\mathfrak{e}, \mathfrak{Z}) = \sum\_{m=0}^{\infty} \mathrm{V}\_{\mathfrak{m}}(\mathfrak{e}, \mathfrak{Z}).$$

*Remember that* UU*ε* = ∑∞ *m*=0 A*m,* UV*ε* = ∑∞ *m*=0 B*m and* VU*ε* = ∑∞ *m*=0 C*m are the Adomian polynomials*

$$\begin{split} &\sum\_{m=0}^{\infty} \mathsf{U}\_{m}(\varepsilon, \mathbb{S}) = \mathsf{U}(\varepsilon, 0) - \mathsf{Y}^{-1} \Big[ (1 + \wp(\upsilon - 1)) \mathbb{Y} \Big{{} \left\{ \sum\_{m=0}^{\infty} \mathcal{A}\_{m} + \frac{1}{2} \frac{\partial \mathrm{U}(\varepsilon, \mathbb{S})}{\partial \varepsilon} + \frac{\partial \mathrm{V}(\varepsilon, \mathbb{S})}{\partial \varepsilon} \right\} \Big], \\ &\sum\_{m=0}^{\infty} \mathsf{V}\_{m}(\varepsilon, \mathbb{S}) = \mathsf{V}(\varepsilon, 0) - \mathsf{Y}^{-1} \Big[ (1 + \wp(\upsilon - 1)) \mathbb{Y} \Big{{} \left\{ \sum\_{m=0}^{\infty} \mathcal{B}\_{m} + \sum\_{m=0}^{\infty} \mathcal{C}\_{m} - \frac{1}{2} \frac{\partial^{2} \mathrm{V}(\varepsilon, \mathbb{S})}{\partial \varepsilon^{2}} \right\} \Big], \end{split}$$

$$\begin{split} \sum\_{m=0}^{\infty} \mathcal{U}\_m(\varepsilon, \mathbb{S}) &= \mathbb{\mathcal{Y}} - \kappa \operatorname{coth}[\kappa(\varepsilon + \theta)] - \mathbb{Y}^{-1} \Big[ (1 + \wp(v - 1)) \mathbb{Y} \Big{{} \Bigg{{} \times} \Big{{}} \mathcal{A}\_m + \frac{1}{2} \frac{\partial \mathcal{U}(\varepsilon, \mathbb{S})}{\partial \varepsilon} + \frac{\partial \mathcal{V}(\varepsilon, \mathbb{S})}{\partial \varepsilon} \Big{{} \Bigg{{} \}} \Big{{} \Bigg{{} \times} \mathcal{A}\_m \\ \sum\_{m=0}^{\infty} \mathcal{V}\_m(\varepsilon, \mathbb{S}) &= -\kappa^2 \operatorname{cscech}^2[\kappa(\varepsilon + \theta)] - \mathbb{Y}^{-1} \Big{{} \Bigg{{} \times} \Big{[} (1 + \wp(v - 1)) \mathbb{Y} \Big{{} \Bigg{{} \times} \Big{[} \mathcal{B}\_m + \sum\_{m=0}^{\infty} \mathcal{C}\_m - \frac{1}{2} \frac{\partial^2 \mathcal{V}(\varepsilon, \mathbb{S})}{\partial \varepsilon^2} \Big{{} \Bigg{{} \-} \mathcal{B}\_m \end{split} \tag{36}$$

*With the aid of the Adomian polynomial according to Equation (22), all forms of non-linear may be stated as*

$$\begin{split} \mathcal{A}\_{0} &= \mathbb{U}\_{0} \frac{\partial \mathbb{U}\_{0}}{\partial \varepsilon}, \quad \mathcal{A}\_{1} = \mathbb{U}\_{0} \frac{\partial \mathbb{U}\_{1}}{\partial \varepsilon} + \mathbb{U}\_{1} \frac{\partial \mathbb{U}\_{0}}{\partial \varepsilon}, \quad \mathcal{B}\_{0} = \mathbb{U}\_{0} \frac{\partial \mathbb{V}\_{0}}{\partial \beta}, \quad \mathcal{B}\_{1} = \mathbb{U}\_{0} \frac{\partial \mathbb{V}\_{1}}{\partial \beta} + \mathbb{U}\_{1} \frac{\partial \mathbb{V}\_{0}}{\partial \beta}, \\ \mathcal{C}\_{0} &= \mathbb{V}\_{0} \frac{\partial \mathbb{U}\_{0}}{\partial \varepsilon}, \quad \mathcal{C}\_{1} = \mathbb{V}\_{0} \frac{\partial \mathbb{U}\_{1}}{\partial \varepsilon} + \mathbb{V}\_{1} \frac{\partial \mathbb{U}\_{0}}{\partial \varepsilon}. \end{split}$$

*Hence, one can easily obtain*

$$
\mathbb{U}\_0(\mathfrak{e}, \mathfrak{z}) = \mathfrak{z} - \kappa \coth[\kappa(\mathfrak{e} + \theta)], \quad \mathbb{V}\_0(\mathfrak{e}, \mathfrak{z}) = -\kappa^2 \mathrm{cosech}^2[\kappa(\mathfrak{e} + \theta)].
$$

*For m* = 0

$$\begin{aligned} \mathbb{U}\_{1}(\varepsilon,\mathbb{S}) &= -\xi\kappa^{2}\mathrm{cosech}^{2}[\kappa(\varepsilon+\theta)]\left\{1+\wp\mathfrak{S}-\wp\right\}, \\ \mathbb{V}\_{1}(\varepsilon,\mathbb{S}) &= -\xi\kappa^{2}\mathrm{cosech}^{2}[\kappa(\varepsilon+\theta)]\,\mathrm{coth}[\kappa(\varepsilon+\theta)]\left\{1+\wp\mathfrak{S}-\wp\right\}. \end{aligned}$$

$$\text{For } m = 1$$

$$\begin{split} \mathcal{U}\_{2}(\varepsilon,\mathbb{Q}) &= \xi\kappa^{4}\operatorname{cscech}^{2}[\kappa(\varepsilon+\theta)] \left\{ 2\xi\kappa \left\{ (1-\wp)^{2}3\wp\odot + (1-\wp)^{3} + \frac{3\wp^{2}(1-\wp)3^{2}}{2} + \frac{\wp^{3}3^{3}}{3!} \right\} \\ &- (3\operatorname{coth}^{2}([\kappa(\varepsilon+\theta)]-1)) \left\{ (1-\wp)2\wp\odot + (1-\wp)^{2} + \frac{\wp^{2}3^{2}}{2} \right\} \right\}, \\ \mathcal{V}\_{2}(\varepsilon,3) &= [2\xi\kappa^{5}\operatorname{cscech}^{2}[\kappa(\varepsilon+\theta)]] \left[ \xi\kappa\operatorname{cscech}^{2}(3\operatorname{coth}^{2}([\kappa(\varepsilon+\theta)]-1)) \left\{ (1-\wp)^{2}3\wp\odot + (1-\wp)^{3} + \frac{3\wp^{2}(1-\wp)3^{2}}{2} \right\} \right] \\ &+ \frac{\wp^{3}3^{3}}{3!} \Big{]} + \frac{2\xi\kappa\operatorname{cscech}^{2}\operatorname{coth}^{2}([\kappa(\varepsilon+\theta)])\Im^{3\wp}}{\Gamma(\wp+1)\Gamma(3\wp+1)} - 2\xi\coth(3\operatorname{cscech}^{2}([\kappa(\varepsilon+\theta)]-1)) \Big{(}1-\wp)2\wp\odot + (1-\wp)^{2} + \frac{\wp^{2}3^{2}}{2} \Big{]} \Big{]}. \end{split}$$

*The remaining steps of the YDM results may be conveniently gathered from* U*m and* V*m* (*m* ≥ 2) *using the same procedure. The alternative series can then be assessed as follows:*

$$\mathbb{U}(\mathfrak{e},\mathfrak{z}) = \sum\_{m=0}^{\infty} \mathbb{U}\_{\mathfrak{m}}(\mathfrak{e},\mathfrak{z}) = \mathbb{U}\_{0}(\mathfrak{e},\mathfrak{z}) + \mathbb{U}\_{1}(\mathfrak{e},\mathfrak{z}) + \mathbb{U}\_{2}(\mathfrak{e},\mathfrak{z}) + \mathbb{U}\_{3}(\mathfrak{e},\mathfrak{z}) + \dotsb \dots$$

$$\mathbb{U}\_{\mathfrak{m}}(\mathfrak{e},\mathfrak{z}) = \sum\_{m=0}^{\infty} \mathbb{U}\_{\mathfrak{m}}(\mathfrak{e},\mathfrak{z}) = \mathbb{U}\_{0}(\mathfrak{e},\mathfrak{z}) + \mathbb{U}\_{1}(\mathfrak{e},\mathfrak{z}) + \mathbb{U}\_{2}(\mathfrak{e},\mathfrak{z}) + \mathbb{U}\_{3}(\mathfrak{e},\mathfrak{z}) + \dotsb \dots$$

U(*ε*, ) = *ξ* − *κ* coth[*κ*(*ε* + *θ*)] − *ξκ*2cosech2[*κ*(*ε* + *θ*)] 1 + ℘ − ℘ + *ξκ*4cosech2[*κ*(*ε* + *θ*)] 2*ξκ*Γ(2℘ + 1) (1 − ℘)23℘ + (1 − ℘)3 + 3℘2(1 − ℘)2 2 + ℘33 3! − (3 coth2([*κ*(*ε* + *θ*)] − 1)) (1 − ℘)2℘ + (1 − ℘)2 + ℘22 2 −··· . V(*ε*, ) = −*κ*2cosech2[*κ*(*ε* + *θ*)] − *ξκ*2cosech2[*κ*(*ε* + *θ*)] coth[*κ*(*ε* + *θ*)] 1 + ℘ − ℘ + [2*ξκ*5cosech2[*κ*(*ε* + *θ*)]] *ξκ*cosech2(3 coth2([*κ*(*ε* + *θ*)] − 1)) (1 − ℘)23℘ + (1 − ℘)3 + 3℘2(1 − ℘)2 2 + ℘33 3! + 2*ξκ*cosech2 coth2([*κ*(*ε* + *θ*)])3℘ Γ(℘ + 1)Γ(3℘ + 1) − 2*ξ* coth(3cosech2([*κ*(*ε* + *θ*)] − 1)) (1 − ℘)2℘ + (1 − ℘)2 + ℘22 2 # −··· .

*We achieve the following series solution at integer order* ℘ = 1, *κ* = 0.1*, ζ* = 0.005 *, θ* = 10*, defined by*

$$\begin{split} \mathcal{V}(\varepsilon,\mathfrak{z}) &= 0.005 - 0.1\coth(0.1\varepsilon + 10) - 0.0005\mathrm{cech}^2(0.1\varepsilon + 10)\mathfrak{z} + 5 \times 10^{-7}\mathrm{cech}\mathrm{sech}^2(0.1\varepsilon + 10) \mathrm{d}0.0033^3 \\ &- 0.5\Big(3\,\mathrm{coth}^2(0.1\varepsilon + 10) - 1.\Big) \mathfrak{z}^2, \\ \mathcal{V}(\varepsilon,\mathfrak{z}) &= -0.01\mathrm{cosh}\mathrm{sech}^2(0.1\varepsilon + 10) - 0.000010\mathrm{cosh}\mathrm{sech}^2(0.1\varepsilon + 10) \times \coth(0.1\varepsilon + 10)\mathfrak{z} + 1.0 \times 10^{-7}\mathrm{cscech}^2(0.1\varepsilon + 10) \\ &\times \Big[ 8.3 \times 10^{-5}\mathrm{3}^3 \mathrm{cscech}^2(0.1\varepsilon + 10)(3\,\mathrm{coth}(0.1\varepsilon + 10) - 1) - 3^2 \mathrm{coth}(0.1\varepsilon + 10) \Big(3\,\mathrm{cscech}^2(0.1\varepsilon + 10) - 1\,\mathrm{coth}(0.1\varepsilon + 10)\mathrm{d}0.0033 \,\mathrm{cscech}^2(0.1\varepsilon + 10)\mathrm{d}0.0033 \,\mathrm{cscech}^2(0.1\varepsilon + 10)\mathrm{d}0.0033 \,\mathrm{cscech}^2(0.1\varepsilon + 10)\mathrm{d}0.0033 \,\mathrm{cscech}^2(0.1\varepsilon + 10)\mathrm{d}0.0033 \,\mathrm{cscech}^2(0.1\varepsilon + 10)\mathrm{d}0.0033 \,\mathrm{cscech}^2(0.1\varepsilon + 10)\mathrm{d}0.0033 \,\mathrm{cscech}^2(0.1\varepsilon + 10)\mathrm{d}0.0033 \,\mathrm{$$

*The exact result of Equation (32) at* ℘ = 1 *and taking ξ* = 0.005*, θ* = 10 *and κ* = 0.1*.*

$$\begin{aligned} \mathcal{U}(\varepsilon, \mathfrak{z}) &= \mathfrak{z} - \kappa \coth[\kappa(\varepsilon + \theta - \mathfrak{z}\mathbb{C})], \\ \mathcal{V}(\varepsilon, \mathfrak{z}) &= -\kappa^2 \mathrm{cosech}^2[\kappa(\varepsilon + \theta - \mathfrak{z}\mathbb{C})]. \end{aligned} \tag{37}$$

*In Figures 5 and 6, the actual and Yang decomposition method solutions at an integerorder* ℘ = 1 *are represented for both* U(*ε*, ) *and* V(*ε*, ) *of Example 1. It is observed that Yang decomposition method results are in good contact with the actual result of the models. In Figures 7 and 8, various fractional-order solutions of Example 2, at different fractional-orders,* ℘ = 1, 0.8, 0.6, 0.4 *are plotted. It is investigated that for Example 2, the fractional-order solutions are convergent to an integer-order solution for both* U(*ε*, ) *and* V(*ε*, )*. In Tables 5 and 6 show that yang decomposition method of different fractional order* ℘ *of Example 2.*

**Figure 5.** The actual and YDM solution of U(*ε*, ) at ℘ = 1 of Example 2.

**Figure 6.** The actual and YDM solution of V(*ε*, ) at ℘ = 1 of Example 2.

**Figure 7.** The fractional-order solutions of V(*ε*, ) at ℘ of Example 2.

**Figure 8.** The fractional-order solutions of V(*ε*, ) at ℘ of Example 2.

**Table 6.** YDM solution of V(*ε*, ) at different fractional-order ℘ of Example 2.

### **5. Conclusions**

This research applies the Yang decomposition method to a fractional-order non-linear Whitham-Broer-Kaup equations system. The suggested approach has been thoroughly researched for fractional-order systems of linear and non-linear differential equations. The numerical results show that the approach is accurate and effective in achieving numerical solutions for non-linear fractional partial differential equations. The proposed methodology is an effective and convenient tool for evaluating numerical solutions to non-linear coupled systems of fractional PDEs compared to previous analytical techniques. Furthermore, the proposed scheme is easy and intuitive, requiring less computing time to solve additional fractional-order partial differential equations.

**Funding:** This work is supported by the Deanship of Scientific Research, King Faisal University, Saudi Arabia (Nasher (third edition) Grant Number: NA000190).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** No data were used to support this study.

**Acknowledgments:** This work is supported by the Deanship of Scientific Research, King Faisal University, Saudi Arabia (Nasher (third edition) Grant Number: NA000190).

**Conflicts of Interest:** The author declare that I have no conflict of interest.

### **References**

### *Article* **Extremal Solutions of Generalized Caputo-Type Fractional-Order Boundary Value Problems Using Monotone Iterative Method**

**Choukri Derbazi 1, Zidane Baitiche 1, Mohammed S. Abdo 2, Kamal Shah 3,4, Bahaaeldin Abdalla 3 and Thabet Abdeljawad 3,5,\***

**Abstract:** The aim of this research work is to derive some appropriate results for extremal solutions to a class of generalized Caputo-type nonlinear fractional differential equations (FDEs) under nonlinear boundary conditions (NBCs). The aforesaid results are derived by using the monotone iterative method, which exercises the procedure of upper and lower solutions. Two sequences of extremal solutions are generated in which one converges to the upper and the other to the corresponding lower solution. The method does not need any prior discretization or collocation for generating the aforesaid two sequences for upper and lower solutions. Further, the aforesaid techniques produce a fruitful combination of upper and lower solutions. To demonstrate our results, we provide some pertinent examples.

**Keywords:** *ϑ*-Caputo derivative; extremal solutions; monotone iterative method; sequences

### **1. Introduction**

Over the last few decades, fractional calculus has attracted the attention of many researchers in the community of science and technology. This is because of its significant applications in different fields of science and engineering such as mathematics, physics, chemistry, biology, economics, finance, rheology, etc. (for more details, see [1–3]). Further, the most important applications of fractional calculus can be found in the description of memory and hereditary processes. The mentioned processes can be more nicely explained by the concept of fractional-order derivatives as compared to traditional ones. Keeping their importance in mind, researchers have given much attention to the use of fractional-order derivatives and integrals in the mathematical modeling of real-world processes instead of classical derivatives and integrals. In this regard, several monographs, and plenty of papers and books have been published, in which various kinds of important results and applications have been reported. Some of these can be found in [4–7]. Nevertheless, the application of the aforesaid area has been traced out just in the last two decades. This is due to the progress in the area of chaos that revealed refined relationships with the concepts of fractional calculus. In addition, in recent times, the application of the theory of fractional calculus to robotics has opened promising aspects for future developments, where in these robots, joint-level control is usually planted by using PID-like schemes with position

**Citation:** Derbazi, C.; Baitiche, Z.; Abdo, M.S.; Shah, K.; Abdalla, B.; Abdeljawad, T. Extremal Solutions of Generalized Caputo-Type Fractional-Order Boundary Value Problems Using Monotone Iterative Method. *Fractal Fract.* **2022**, *6*, 146. https:// doi.org/10.3390/fractalfract6030146

Academic Editors: Asifa Tassaddiq and Muhammad Yaseen

Received: 6 February 2022 Accepted: 3 March 2022 Published: 7 March 2022

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

**Copyright:** © 2022 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 (https:// creativecommons.org/licenses/by/ 4.0/).

feedback. For instance, one famous machine using the fractional *PDα* controller is known as the Hexapod robot (see [8]). Further, some important applications of fractional calculus can be found about the dynamics in the trajectory control of redundant manipulators, where fractional-order derivatives have a more precise appearance than classical ones (see [9]). Further, the fractional derivative is a global operator that preserves greater degrees of freedom as compared to integer-order derivative, because a classical derivative with integer-order is a local operator. It is estimated that the fractional-order derivative operation contains some types of boundary conditions that involve information on the function further out. Some researchers have proved that fractional-order derivatives play a significant role in electrochemical analysis to elucidate the mechanistic behavior of the concentration of a substrate at the electrode surface to the current. Some researchers have proved that fractional-order modelers of contaminant flow in heterogeneous porous media involving fractional derivatives are more powerful than classical ones (see [10]). It should be kept in mind that fractional-order derivative operation of a function produces a complete spectrum or accumulation, which preserves the corresponding integer-order counterpart as a special case.

It is interesting that fractional-order derivatives have not yet been uniquely defined. Various renowned mathematicians have given their own definitions. Among the said definitions, some of them have gained much more popularity and proper attention from researchers, such as Riemann–Liouville (RL), Caputo, Hilfer, Caputo–Hadamard, Caputo– Katugampola, etc. It is interesting that aside from the aforesaid operators various other variants that contain singular kernels have been introduced recently. Hence the fractional differential operators have been divided into two classes including singular and non-singular. Here we state that this partition is not bad but provides a great degree of freedom in the choice of operator for the description of a particular phenomenon. It is remarkable that nearly all mentioned operators preserve memory in their respective kernels. Further, the two partitions of singular and non-singular kernels have their own benefits and drawbacks.

We remark that Hilfer-type fractional calculus unifies the aforesaid definitions. It is important that *ϑ*-Hilfer operators constitute a wide class of fractional derivatives (FDs). In this respect, some frequent results involving *ϑ*-FDs have been reported in [11–13]. However, various strategies exist in the literature to handle such types of problems of FDEs for computation of their solutions. Numerous tools and theories have been established so far. Iterative techniques of various kinds have key importance to investigate the aforementioned area. Among them are the monotone iterative algorithm, along with the upper and lower solutions method [14,15], fixed-point technique [16–18], and coincidence degree theory [19,20]. In particular, the monotone iterative method together with the technique of upper and lower solutions is an advantageous and effective tool for the existence as well as the approximation of solutions for nonlinear problems. In this regard, very useful results have been published so far. Among the iterative techniques, those introduced by Ladde, Lakshmikantham, and Vatsala [21] in 1985 for nonlinear differential equations have gained proper attention. Therefore, monotone iterative techniques associated with upper and lower solutions have been extensively used for nonlinear partial differential equations in the last few decades. In this regard, plenty of work has been published to date; a few can be found in [22–25]. Proposals have been made for classical differential equations for the first time [14,15,26–29]. In addition, the aforesaid techniques have been widely used to deal with FDEs subject to initial and boundary conditions. Some significant results can be found in [25,30–37]. We demonstrate that the said method is well known because it not only produces constructive proof for existence theorems but it also yields various comparison results, which are powerful tools to investigate the qualitative properties of solutions. Further, the sequence of iterations has useful behavior in the computation of numerical solutions to various boundary value and initial boundary value problems of classical as well FDEs. In addition, the method of upper and lower solutions is very useful

for the construction of Lyapunov functions. Construction of such functions is increasingly used to derive various stability theories in dynamical systems.

Motivated by the work cited above, and particularly by the work in [31], we determine the existence criteria of extremal solutions for the following *ϑ*-Caputo-type FDE in a Caputo sense with NBCs

$$\begin{cases} \,^c \mathbb{D}\_{a^+}^{\nu; \theta} \phi(\emptyset) = \mathbb{F}(\emptyset, \phi(\emptyset)), \; \emptyset \in \mathcal{I} := [a, b],\\ \mathbb{G}(\phi(a), \phi(b)) = 0, \end{cases} \tag{1}$$

where 0 < *ν* ≤ 1, *c*D*ν*;*ϑ a*+ is the *ϑ*-fractional operator of order *ν* in the Caputo sense and this is investigated. Further, F ∈ *C*(J × R, R), G ∈ *C*(R2, R).

The considered problem (1) in the current article includes a wide range of nonlinear FDEs involving the standard Caputo operator (for *ϑ*(*ς*) = *ς*), Caputo–Hadamard (for *ϑ*(*ς*) = log *ς*), and Caputo–Katugampola (for *ϑ*(*ς*) = *ς*p, p > 0). Further, fractional operators have been listed in Almeida [11] for further applications. Further, results acquired in the current article include the results of Franco et al. [26] if *a* → 0, *ϑ*(*ς*) → *ς*, and *ν* = 1 as a special case.

In this regard, we also point out some recent and similar findings that used operators on many fractional problems; see [38–40]. To the best of our knowledge in this regard, no one has considered the monotone iterative procedure to obtain the existence of extremal solutions involving a *ϑ*-Caputo derivative subject to NBCs. Therefore, motivated by the aforesaid gap, we have conducted this study.

The rest of this article is organized as follows. In Section 2, we insert some basic definitions and important results. Section 3 is devoted to studying the existence of extremal solutions for (1). In Section 4, we give two appropriate examples to highlight the feasibility of our abstract results.

### **2. Basic Results**

Some fundamental results about the *ϑ*-Caputo derivative and integral that are needed throughout this work are given below.

The function *ϑ* ∈ *C*(J, R) is non-decreasing differentiable with argument 0 < *ϑ* (*ς*), at every point of J.

**Definition 1** ([6,11])**.** *The ϑ-RL fractional integral of order ν* > 0 *for an integrable function φ*: J −→ R *is given by*

$$\mathbb{I}\_{a^{+}}^{\upsilon;\theta}\phi(\boldsymbol{\varsigma}) = \frac{1}{\Gamma(\upsilon)} \int\_{a}^{\varsigma} \theta'(s)(\theta(\boldsymbol{\varsigma}) - \theta(s))^{\upsilon - 1} \phi(s) \mathrm{d}s, \; \boldsymbol{\varsigma} > a.$$

**Definition 2** ([11])**.** *Let ϑ*, *φ* ∈ *Cn*(J, R)*. The ϑ-RL derivative of fractional order of a function φ with n* − 1 < *ν* < *n is given by*

$$\begin{aligned} \, \_0\mathbb{D}\_{a^+}^{\upsilon;\theta} \phi(\varsigma) &= \, \_0\left( \frac{D\_t}{\theta'(\varsigma)} \right)^n \, \_{a^+}^{\upsilon - \upsilon;\theta} \phi(\varsigma) \\ &= \, \_0\frac{1}{\Gamma(n-\nu)} \left( \frac{D\_t}{\theta'(\varsigma)} \right)^n \int\_a^{\varsigma} \theta'(s) (\theta(\varsigma) - \theta(s))^{n-\nu-1} \phi(s) \, \mathrm{d}s \, \_0\mathbb{D}\_{a^+}^{\upsilon} \phi \, \mathrm{d}s \end{aligned}$$

*where n* = [*ν*] + 1 *(n* ∈ N*), and Dt* = *d dt* .

**Definition 3** ([11])**.** *Let ϑ*, *φ* ∈ *Cn*(J, R)*. The ϑ-Caputo derivative of fractional order of function φ with n* − 1 < *ν* < *n is defined by*

$$\prescript{\circ}{}{\mathbb{D}}\_{a^{+}}^{\nu;\theta}\phi(\emptyset) = \prescript{\imath}{}{\mathbb{I}}\_{a^{+}}^{\nu-\nu;\theta}\phi\_{\theta}^{[n]}(\emptyset),$$

*where φ*[*n*] *ϑ* (*ς*) = *Dς ϑ*(*ς*) *n φ*(*ς*), *n* = [*ν*] + 1 *for ν* ∈/ N*, and n* = *ν for ν* ∈ N*.* *One has*

$$\prescript{\circ}{}{\mathbb{D}}\_{a^{+}}^{\upsilon;\theta} \phi(\boldsymbol{\varsigma}) = \begin{cases} \int\_{a}^{\mathbb{G}} \frac{\theta'(s)(\theta(\boldsymbol{\varsigma}) - \theta(\boldsymbol{s}))^{n-\upsilon-1}}{\Gamma(n-\upsilon)} \phi\_{\boldsymbol{\theta}}^{[n]}(\boldsymbol{s}) \, \mathrm{d}\boldsymbol{s} & \text{if } \upsilon \notin \mathbb{N}\_{\boldsymbol{\theta}},\\ \phi\_{\boldsymbol{\theta}}^{[n]}(\boldsymbol{\varsigma}) & \text{, } \quad \text{if } \upsilon \in \mathbb{N}. \end{cases}$$

**Lemma 1** ([11])**.** *Let ν*, *μ* > 0, *and φ* ∈ *C*(J, R), *for every ς* ∈ J

*1. c*D*ν*;*ϑ a*+ I *ν*;*ϑ a*+ *φ*(*ς*) = *φ*(*ς*)*, 2.* I *ν*;*ϑ a*+ *c*D*ν*;*ϑ a*+ *φ*(*ς*) = *φ*(*ς*) − *φ*(*a*), 0 < *ν* ≤ 1*, 3.* I *ν*;*ϑ a*+ (*ϑ*(*ς*) − *ϑ*(*a*))*μ*−1 = Γ(*μ*) Γ(*μ*+*ν*)(*ϑ*(*ς*) − *ϑ*(*a*))*μ*+*ν*−1,

$$4. \quad \, ^c \mathbb{D}\_{a^+}^{\nu;\theta} (\theta(\zeta) - \theta(a))^{\mu - 1} = \frac{\Gamma(\mu)}{\Gamma(\mu - \nu)} (\theta(\zeta) - \theta(a))^{\mu - \nu - 1}.$$

$$\text{5.} \quad \, ^c \mathbb{D}\_{a^+}^{\nu; \theta} (\theta(\emptyset) - \theta(a))^k = 0, \,\, ^\psi k < n \in \mathbb{N}.$$

**Definition 4** ([7])**.** *One- and two-parameter Mittag–Leffler functions (MLFs) are recalled as*

$$\mathbb{E}\_{\nu}(z) = \sum\_{k=0}^{\infty} \frac{z^k}{\Gamma(\nu k + 1)}, \quad (z \in \mathbb{R}, \: \nu > 0),$$

*and*

$$\mathbb{E}\_{\nu,\mu}(z) = \sum\_{k=0}^{\infty} \frac{z^k}{\Gamma(\nu k + \mu)}, \ \nu, \mu > 0 \text{ and } z \in \mathbb{R}\_{\nu}$$

*respectively. Clearly* E1,1(*z*) = E1(*z*) = exp(*z*)*.*

Further properties of MLFs are given below.

**Lemma 2** ([41])**.** *Let ν* ∈ (0, 1) *and z* ∈ R, *one has*

For further analysis, we recall the following Lemma [31] as:

**Lemma 3** ([31] (Lemma 4))**.** *Let ν* ∈ (0, 1]*, λ* ∈ R *and h* ∈ *C*(J, R), *then, the linear version*

$$\begin{cases} ^c\mathbb{D}\_{a^+}^{\vee;\theta}\phi(\emptyset) + \lambda\phi(\emptyset) = h(\emptyset), \ \emptyset \in \mathcal{J},\\ \phi(a) = \phi\_{a\prime} \end{cases}$$

*has a unique solution that is described as*

$$\begin{split} \boldsymbol{\phi}(\boldsymbol{\xi}) &= \theta\_{d} \mathbb{E}\_{\boldsymbol{\nu},1} \Big( -\lambda \big( \boldsymbol{\theta}(\boldsymbol{\xi}) - \boldsymbol{\theta}(\boldsymbol{a}) \big)^{\boldsymbol{\nu}} \big) \\ &+ \int\_{\boldsymbol{a}}^{\boldsymbol{\xi}} \boldsymbol{\theta}'(\boldsymbol{s}) (\boldsymbol{\theta}(\boldsymbol{\xi}) - \boldsymbol{\theta}(\boldsymbol{s}))^{\boldsymbol{\nu}-1} \mathbb{E}\_{\boldsymbol{\nu},\boldsymbol{\nu}} \Big( -\lambda (\boldsymbol{\theta}(\boldsymbol{\xi}) - \boldsymbol{\theta}(\boldsymbol{s}))^{\boldsymbol{\nu}} \big) h(\boldsymbol{s}) \mathrm{d}\boldsymbol{s}, \end{split}$$

*where* E*ν*,*μ*(·) *is the two-parametric MLF defined earlier.*

The given comparison results comprise a central rule in the following analysis.

**Lemma 4** ([31] (Lemma 5))**.** *Let ν* ∈ (0, 1]*, and λ* ∈ R*, if γ* ∈ *C*(J, R) *obey the given relation*

$$\begin{cases} \,^c \mathbb{D}\_{a^+}^{\upsilon;\theta} \gamma(\emptyset) \ge -\lambda \gamma(\emptyset), & \emptyset \in (a,b],\\ \gamma(a) \ge 0, & \end{cases}$$

*then γ*(*ς*) ≥ 0, *for all ς* ∈ J*.*

### **3. Main Results**

Here, key findings are established about the *ϑ*-Caputo FDEs (1). We develop two monotone iterative sequences for upper and lower solutions, respectively.

**Definition 5.** *The function φ* ∈ *C*(J, R) *such that c*D*ν*;*ϑ a*+ *φ exists and is continuous on* J *and is known to be a solution of the problem* (1)*. Further, φ gives the statistics of the equation c*D*ν*;*ϑ a*+ *φ*(*ς*) = F(*ς*, *φ*(*ς*))*, for each ς* ∈ J *and the NBCs*

$$\mathbb{G}(\phi(a), \phi(b)) = 0.$$

Subsequently, we mention the definitions of extremal solutions of *ϑ*-Caputo FDEs (1).

**Definition 6.** *The mapping φ*0 ∈ *C*(J, R) *is known as lower solution* (1)*, if it satisfies*

$$\begin{cases} ^c\mathbb{D}\_{a^+}^{\nu;\theta}\phi\_0(\emptyset) \le \mathbb{F}(\emptyset, \phi\_0(\emptyset)), \quad t \in (a, b],\\ ^c\mathbb{G}(\phi\_0(a), \phi\_0(b)) \le 0. \end{cases}$$

*An upper solution* 0 ∈ *C*(J, R) *of the problem* (1) *can be defined in a similar way by reversing the above inequality.*

Now to move forward, we will introduce the following conditions:

**Hypothesis 1.** *There exist φ*0 *and* 0 *as lower and upper solutions of problem* (1) *in C*(J, R) *respectively, with φ*0(*ς*) ≤ 0(*ς*), *ς* ∈ J*.*

**Hypothesis 2.** *There exists a constant λ* > 0 *with*

$$\mathbb{F}(\mathfrak{g}, \mathbf{y}) - \mathbb{F}(\mathfrak{g}, \mathbf{x}) \ge -\lambda(\mathbf{y} - \mathbf{x}) \quad \text{for } \mathfrak{g}\_0(\mathfrak{g}) \le \mathbf{x} \le \mathbf{y} \le \mathfrak{a}\_0(\mathfrak{g}), \mathfrak{g} \in \mathcal{I}.$$

**Hypothesis 3.** *There exist constants c* > 0 *and d* ≥ 0 *with φ*0(*a*) ≤ *u*1 ≤ *u*2 ≤ 0(*a*)*, φ*0(*b*) ≤ *v*1 ≤ *v*2 ≤ 0(*b*), *such that*

$$
\mathbb{G}(\mathfrak{u}\_2, \mathfrak{v}\_2) - \mathbb{G}(\mathfrak{u}\_1, \mathfrak{v}\_1) \le c(\mathfrak{u}\_2 - \mathfrak{u}\_1) - d(\mathfrak{v}\_2 - \mathfrak{v}\_1)\_\prime
$$

Now, we shall apply the monotonous method to prove our key findings.

**Theorem 1.** *Let* F : J × R −→ R *be continuous. Assume that Hypotheses 1–3 hold. Then there exist two monotone iterative sequences* {*φn*} *and* {*n*}*, which are converging uniformly on* J *to the extremal solutions of* (1) *in the sector* [*φ*0, 0]*, where*

$$[\mathfrak{q}\_0, \mathfrak{q}\_0] = \{ z \in \mathbb{C}(\mathfrak{J}, \mathbb{R}) : \mathfrak{q}\_0(\mathfrak{g}) \le z(\mathfrak{g}) \le \mathfrak{q}\_0(\mathfrak{g}), \quad \mathfrak{g} \in \mathfrak{J} \}.$$

**Proof.** First, for any *φ*0(*ς*), 0(*ς*) ∈ *C*(J, R) and *λ* > 0, we consider the following FDEs

$$\begin{cases} ^c\mathbb{D}\_{a^+}^{\nu;\theta}\phi\_{n+1}(\emptyset) = \mathbb{F}(\emptyset, \phi\_n(\emptyset)) - \lambda \left(\phi\_{n+1}(\emptyset) - \phi\_n(\emptyset)\right), & \emptyset \in \mathcal{J},\\ \phi\_{n+1}(a) = \phi\_n(a) - \frac{1}{c}\mathbb{G}(\phi\_n(a), \phi\_n(b)), & \end{cases} \tag{2}$$

and

$$\begin{cases} \,^c \mathbb{D}\_{a^+}^{\nu, \theta} \mathcal{o}\_{n+1}(\boldsymbol{\xi}) = \mathbb{F}(\boldsymbol{\xi}, \boldsymbol{\mathcal{o}}\_n(\boldsymbol{\xi})) - \lambda (\boldsymbol{\mathcal{o}}\_{n+1}(\boldsymbol{\xi}) - \boldsymbol{\mathcal{o}}\_n(\boldsymbol{\xi})), & \boldsymbol{\xi} \in \mathcal{J}, \\ \boldsymbol{\mathcal{o}}\_{n+1}(a) = \boldsymbol{\mathcal{o}}\_n(a) - \frac{1}{c} \mathbb{G}(\boldsymbol{\mathcal{o}}\_n(a), \boldsymbol{\mathcal{o}}\_n(b)), & \end{cases} \tag{3}$$

According to Lemma 3, one can deduce that (2) and (3) preserve at most one solution in *C*(J, R). Thus we have

$$\begin{split} \phi\_{n+1}(\boldsymbol{\xi}) &= \left(\phi\_{n}(\boldsymbol{a}) - \frac{1}{\boldsymbol{c}} \mathbb{G}(\phi\_{n}(\boldsymbol{a}), \phi\_{n}(\boldsymbol{b}))\right) \mathbb{E}\_{\boldsymbol{v},1} \Big( -\lambda \big( \boldsymbol{\theta}(\boldsymbol{\zeta}) - \boldsymbol{\theta}(\boldsymbol{a}) \big)^{\boldsymbol{\nu}} \big) \\ &+ \int\_{\boldsymbol{a}}^{\boldsymbol{\xi}} \boldsymbol{\theta}'(\boldsymbol{s}) (\boldsymbol{\theta}(\boldsymbol{\zeta}) - \boldsymbol{\theta}(\boldsymbol{s}))^{\boldsymbol{\nu} - 1} \mathbb{E}\_{\boldsymbol{v},\boldsymbol{\nu}} \Big( -\lambda \big( \boldsymbol{\theta}(\boldsymbol{\zeta}) - \boldsymbol{\theta}(\boldsymbol{s}) \big)^{\boldsymbol{\nu}} \big) \big( \mathbb{F}(\boldsymbol{s}, \phi\_{n}(\boldsymbol{s}) \big) + \lambda \boldsymbol{\phi}\_{\boldsymbol{n}}(\boldsymbol{s}) \big) \mathrm{d}\boldsymbol{s}, \ \boldsymbol{\xi} \in \operatorname{J}, \\\\ \boldsymbol{\mathcal{O}}\boldsymbol{\sigma}\_{n+1}(\boldsymbol{\xi}) &= \left(\boldsymbol{\varpi}\_{n}(\boldsymbol{a}) - \frac{1}{\boldsymbol{c}} \mathbb{G}(\boldsymbol{\omega}\_{n}(\boldsymbol{a}), \boldsymbol{\varpi}\_{n}(\boldsymbol{b})) \right) \big[ \mathbb{E}\_{\boldsymbol{v},\boldsymbol{\nu}} \Big( -\lambda \big( \boldsymbol{\theta}(\boldsymbol{\zeta}) - \boldsymbol{\theta}(\boldsymbol{a}) \big)^{\boldsymbol{v}} \big) \\ &+ \int\_{\boldsymbol{a}}^{\boldsymbol{\xi}} \boldsymbol{\theta}'(\boldsymbol{s}) (\boldsymbol{\vartheta}(\boldsymbol{\xi}) - \boldsymbol{\theta}(\boldsymbol{s}))^{\boldsymbol{v} - 1} \mathbb{E}\_{\boldsymbol{v},\boldsymbol{\nu}} \big( -\lambda \big( \$$

For appropriateness, the proof will be divided into a number of steps.

**Step 1**: The sequences *φn*(*ς*), *n*(*ς*)(*n* ≥ 1) are lower and upper solutions of (1), correspondingly. Moreover, we assume

$$\phi\_0(\emptyset) \le \phi\_1(\emptyset) \le \dots \le \phi\_n(\emptyset) \le \dots \le \phi\_n(\emptyset) \le \dots \le \phi\_1(\emptyset) \le \phi\_0(\emptyset), \quad \emptyset \in \mathcal{I}. \tag{4}$$

Firstly, we show that

$$
\phi\_0(\emptyset) \le \phi\_1(\emptyset) \le \phi\_1(\emptyset) \le \phi\_0(\emptyset), \quad \emptyset \in \mathcal{I}.
$$

Set *γ*(*ς*) = *φ*1(*ς*) − *φ*0(*ς*). From (2) and Definition 6, we obtain

$$\begin{split} \, ^c \mathbb{D}\_{a^+}^{\nu; \theta} \gamma(\emptyset) &= \, ^c \mathbb{D}\_{a^+}^{\nu; \theta} \Phi\_1(\emptyset) - \, ^c \mathbb{D}\_{a^+}^{\nu; \theta} \Phi\_0(\emptyset) \\ &\geq \, ^c \mathbb{F}(\emptyset, \phi\_0(\emptyset)) - \lambda \, (\phi\_1(\emptyset) - \phi\_0(\emptyset)) - \mathbb{F}(\emptyset, \phi\_0(\emptyset)) \\ &= -\lambda \, \gamma(\emptyset). \end{split}$$

Again, since *γ*(*a*) = −1 *c*G(*φ*0(*a*), *φ*0(*b*)) ≥ 0, *γ*(*ς*) ≥ 0, for *ς* ∈ J due to Lemma 4. Thus, *φ*0(*ς*) ≤ *φ*1(*ς*).

Similarly, we can find that 1(*ς*) ≤ 0(*ς*), *ς* ∈ J.

Now, let *γ*(*ς*) = 1(*ς*) − *φ*1(*ς*). Using (2) and (3) together with Hypotheses 2 and 3, we obtain

$$\begin{split} \mathbb{E}^{\mathbb{V},\mathcal{O}}\mathcal{D}^{\nu,\mathfrak{p}}\_{\mathfrak{q}^{+}}\gamma(\mathfrak{q}) &= \mathbb{F}\{\mathfrak{z},\mathfrak{o}\_{0}(\mathfrak{q})\} - \mathbb{F}\{\mathfrak{z},\mathfrak{q}\_{0}(\mathfrak{q})\} - \lambda\left(\mathfrak{o}\_{1}(\mathfrak{z}) - \mathfrak{o}\_{0}(\mathfrak{q})\right) + \lambda\left(\mathfrak{o}\_{1}(\mathfrak{q}) - \mathfrak{q}\_{0}(\mathfrak{q})\right) \\ &\geq -\lambda\left(\mathfrak{o}\_{0}(\mathfrak{q}) - \mathfrak{q}\_{0}(\mathfrak{q})\right) - \lambda\left(\mathfrak{o}\_{1}(\mathfrak{q}) - \mathfrak{o}\_{0}(\mathfrak{q})\right) + \lambda\left(\mathfrak{o}\_{1}(\mathfrak{q}) - \mathfrak{q}\_{0}(\mathfrak{q})\right) \\ &= -\lambda\gamma(\mathfrak{q}). \end{split}$$

Since

$$\begin{aligned} \gamma(a) &= \left(\varpi\_0(a) - \phi\_0(a)\right) - \frac{1}{\mathcal{c}} \left(\mathbb{G}\left(\varpi\_0(a), \varpi\_0(b)\right) - \mathbb{G}\left(\phi\_0(a), \phi\_0(b)\right)\right) \\ &\ge \frac{d}{\mathcal{c}} \left(\varpi\_0(b) - \phi\_0(b)\right) \ge 0, \end{aligned}$$

we obtain *φ*1(*ς*) ≤ 1(*ς*), *ς* ∈ J due to Lemma 4.

Secondly, we show that *φ*1(*ς*), 1(*ς*) are extremum solutions of (1). Since *φ*0 and 0 are lower and upper solutions of (1), by Hypotheses 2 and 3, we obtain

$$\mathbb{F}^{\mathfrak{c}} \mathbb{D}\_{a^{+}}^{\nu; \mathfrak{b}} \phi\_1(\mathfrak{g}) = \mathbb{F}(\mathfrak{g}, \mathfrak{so}(\mathfrak{g})) - \lambda \left(\phi\_1(\mathfrak{g}) - \mathfrak{so}(\mathfrak{g})\right) \le \mathbb{F}(\mathfrak{g}, \phi\_1(\mathfrak{g})),$$

and

$$\begin{aligned} \mathbb{G}(\phi\_1(a), \phi\_1(b)) &\leq \mathbb{G}(\phi\_0(a), \phi\_0(b)) + c(\phi\_1(a) - \phi\_0(a)) - d(\phi\_1(b) - \phi\_0(b)) \\ &= -d(\phi\_1(b) - \phi\_0(b)) \\ &\leq 0. \end{aligned}$$

Therefore, *φ*1(*ς*) is a lower solution of (1). Analogously, it is obvious that 1(*ς*) is an upper solution of (1).

Through the above debates and induction, we can show that *φn*(*ς*), *n*(*ς*),(*n* ≥ 1) are lower and upper solutions of (1), respectively, and the assumption (4) is true.

**Step 2:** *φn* → *φ* and *n* → .

First, we prove that {*φn*} is uniformly bounded. By considering supposition Hypothesis 2, we may write

$$\mathbb{P}(\mathfrak{g}, \mathfrak{\phi}\_0(\mathfrak{g})) + \lambda \mathfrak{\phi}\_0(\mathfrak{g}) \le \mathbb{P}(\mathfrak{g}, \mathfrak{\phi}\_n(\mathfrak{g})) + \lambda \mathfrak{\phi}\_n(\mathfrak{g}) \le \mathbb{P}(\mathfrak{g}, \mathfrak{a}\_0(\mathfrak{g})) + \lambda \mathfrak{a}\_0(\mathfrak{g}), \quad \mathfrak{g} \in \mathfrak{f}$$

i.e.,

$$\begin{aligned} 0 &\leq \mathbb{F}(\mathfrak{c}, \mathfrak{\phi}\_{n}(\mathfrak{c})) - \mathbb{F}(\mathfrak{c}, \mathfrak{\phi}\_{0}(\mathfrak{c})) + \lambda (\mathfrak{\phi}\_{n}(\mathfrak{c}) - \mathfrak{\phi}\_{0}(\mathfrak{c})) \\ &\leq \mathbb{F}(\mathfrak{c}, \mathfrak{\phi}\_{0}(\mathfrak{c})) - \mathbb{F}(\mathfrak{c}, \mathfrak{\phi}\_{0}(\mathfrak{c})) + \lambda (\mathfrak{\phi}\_{0}(\mathfrak{c}) - \mathfrak{\phi}\_{0}(\mathfrak{c})). \end{aligned}$$

Hence, we obtain

$$\begin{split} |\mathbb{F}(\mathfrak{f},\mathfrak{q}\_{\mathfrak{n}}(\mathfrak{c})) - \mathbb{F}(\mathfrak{c},\mathfrak{q}\_{\mathfrak{0}}(\mathfrak{c})) + \lambda(\mathfrak{q}\_{\mathfrak{n}}(\mathfrak{c}) - \mathfrak{q}\_{\mathfrak{0}}(\mathfrak{c}))| &\leq |\mathbb{F}(\mathfrak{g},\mathfrak{a}\_{\mathfrak{0}}(\mathfrak{c})) - \mathbb{F}(\mathfrak{g},\mathfrak{q}\_{\mathfrak{0}}(\mathfrak{c}))| \\ &+ \lambda(\mathfrak{a}\mathfrak{o}(\mathfrak{c}) - \mathfrak{q}\mathfrak{n}(\mathfrak{c}))|. \end{split}$$

Thus

$$\begin{split} |\mathbb{F}(\xi,\phi\_{n}(\xi)) + \lambda\phi\_{n}(\xi)| &\leq |\mathbb{F}(\xi,\phi\_{n}(\xi)) - \mathbb{F}(\xi,\phi\_{0}(\xi)) + \lambda\left(\phi\_{n}(\xi) - \phi\_{0}(\xi)\right)| \\ &\quad + |\mathbb{F}(\xi,\phi\_{0}(\xi)) + \lambda\phi\_{0}(\xi)| \\ &\leq |\mathbb{F}(\xi,\phi\_{0}(\xi)) - \mathbb{F}(\xi,\phi\_{0}(\xi)) + \lambda\left(\phi\_{0}(\xi) - \phi\_{0}(\xi)\right)| \\ &\quad + |\mathbb{F}(\xi,\phi\_{0}(\xi)) + \lambda\phi\_{0}(\xi)| \\ &\leq 2|\mathbb{F}(\xi,\phi\_{0}(\xi)) + \lambda\phi\_{0}(\xi)| + |\mathbb{F}(\xi,\phi\_{0}(\xi)) + \lambda\phi\_{0}(\xi)|. \end{split}$$

Since *φ*0, F are continuous on J, we can see a constant M independent of *n* with

$$|\mathbb{F}(\mathfrak{g}\_{\prime}\phi\_{n}(\mathfrak{g})) + \lambda\phi\_{n}(\mathfrak{g})| \leq \mathbb{M}.\tag{5}$$

Furthermore, from Hypothesis 3, we can obtain

$$
\phi\_0(a) - \frac{1}{\mathcal{c}} \mathbb{G}(\phi\_0(a), \phi\_0(b)) \le \phi\_n(a) - \frac{1}{\mathcal{c}} \mathbb{G}(\phi\_n(a), \phi\_n(b)) \le \mathcal{a}\_0(a) - \frac{1}{\mathcal{c}} \mathbb{G}(\phi\_0(a), \phi\_0(b)),
$$

i.e.,

$$\begin{aligned} 0 &\leq \phi\_n(a) - \phi\_0(a) - \frac{1}{c} \left( \mathbb{G}(\phi\_n(a), \phi\_n(b)) - \mathbb{G}(\phi\_0(a), \phi\_0(b)) \right) \\ &\leq \phi\_0(a) - \phi\_0(a) - \frac{1}{c} \left( \mathbb{G}(\phi\_0(a), \phi\_0(b)) - \mathbb{G}(\phi\_0(a), \phi\_0(b)) \right) .\end{aligned}$$

Hence, we obtain

$$\begin{aligned} & \left| \phi\_{\mathfrak{n}}(a) - \phi\_{\mathfrak{0}}(a) - \frac{1}{\mathfrak{c}} \left( \mathbb{G}(\phi\_{\mathfrak{n}}(a), \phi\_{\mathfrak{n}}(b)) - \mathbb{G}(\phi\_{\mathfrak{0}}(a), \phi\_{\mathfrak{0}}(b)) \right) \right| \leq \\ & \left| \begin{aligned} & \left| \mathscr{O}\_{\mathfrak{0}}(a) - \phi\_{\mathfrak{0}}(a) - \frac{1}{\mathfrak{c}} \left( \mathbb{G}(\mathfrak{c}\_{\mathfrak{0}}(a), \mathscr{O}\_{\mathfrak{0}}(b)) - \mathbb{G}(\phi\_{\mathfrak{0}}(a), \phi\_{\mathfrak{0}}(b)) \right) \right| \\ & \leq \left| \phi\_{\mathfrak{0}}(a) - \frac{1}{\mathfrak{c}} \mathbb{G}(\phi\_{\mathfrak{0}}(a), \phi\_{\mathfrak{0}}(b)) \right| + \left| \mathscr{O}\_{\mathfrak{0}}(a) - \frac{1}{\mathfrak{c}} \mathbb{G}(\phi\_{\mathfrak{0}}(a), \phi\_{\mathfrak{0}}(b)) \right|. \end{aligned} \end{aligned}$$

Thus

$$\begin{split} \left| \phi\_{\boldsymbol{n}}(\boldsymbol{a}) - \frac{1}{\boldsymbol{c}} \mathbb{G} (\phi\_{\boldsymbol{n}}(\boldsymbol{a}), \phi\_{\boldsymbol{n}}(\boldsymbol{b})) \right| &\leq \left| \phi\_{\boldsymbol{n}}(\boldsymbol{a}) - \phi\_{0}(\boldsymbol{a}) - \frac{1}{\boldsymbol{c}} \Big( \mathbb{G} (\phi\_{\boldsymbol{n}}(\boldsymbol{a}), \phi\_{\boldsymbol{n}}(\boldsymbol{b})) - \mathbb{G} (\phi\_{0}(\boldsymbol{a}), \phi\_{0}(\boldsymbol{b})) \Big) \right| \\ &+ \left| \phi\_{0}(\boldsymbol{a}) - \frac{1}{\boldsymbol{c}} \mathbb{G} (\phi\_{0}(\boldsymbol{a}), \phi\_{0}(\boldsymbol{b})) \right| \\ &\leq 2 \left| \phi\_{0}(\boldsymbol{a}) - \frac{1}{\boldsymbol{c}} \mathbb{G} (\phi\_{0}(\boldsymbol{a}), \phi\_{0}(\boldsymbol{b})) \right| + \left| \phi\_{0}(\boldsymbol{a}) - \frac{1}{\boldsymbol{c}} \mathbb{G} (\phi\_{0}(\boldsymbol{a}), \phi\_{0}(\boldsymbol{b})) \Big| \right|. \end{split}$$

Since *φ*0, 0 and G are continuous functions, we can find a constant K independent of *n*, such that

$$\left|\phi\_n(a) - \frac{1}{c} \mathbb{G}(\phi\_n(a), \phi\_n(b))\right| \le \mathbb{K}.\tag{6}$$

Moreover, by (2) and (3) we have

$$\begin{split} \left| \phi\_{n+1}(\boldsymbol{\xi}) \right| &= \left| \phi\_{n}(\boldsymbol{a}) - \frac{1}{c} \mathbb{G} (\phi\_{n}(\boldsymbol{a}), \phi\_{n}(\boldsymbol{b})) \right| \mathbb{E}\_{\boldsymbol{\nu},1} \Big( -\lambda (\boldsymbol{\theta}(\boldsymbol{\xi}) - \boldsymbol{\theta}(\boldsymbol{a}))^{\boldsymbol{\nu}} \Big) \\ &+ \int\_{\mathfrak{s}}^{\boldsymbol{\xi}} \boldsymbol{\theta}'(\boldsymbol{s}) (\boldsymbol{\theta}(\boldsymbol{\xi}) - \boldsymbol{\theta}(\boldsymbol{s}))^{\boldsymbol{\nu}-1} \mathbb{E}\_{\boldsymbol{\nu},\boldsymbol{\nu}} \Big( -\lambda (\boldsymbol{\theta}(\boldsymbol{\xi}) - \boldsymbol{\theta}(\boldsymbol{s}))^{\boldsymbol{\nu}} \Big) \Big| \mathbb{F} (\boldsymbol{s}, \phi\_{n}(\boldsymbol{s})) + \lambda \phi\_{n}(\boldsymbol{s}) \Big| \mathrm{d}s, \; \boldsymbol{\xi} \in \mathcal{J}. \end{split}$$

Using Lemma 2 along with (5) and (6), we obtain

$$\begin{split} |\phi\_{n+1}(\boldsymbol{\xi})| &= \mathbb{K} + \frac{\mathbb{M}}{\Gamma(\boldsymbol{\nu})} \int\_{a}^{\boldsymbol{\xi}} \theta'(\boldsymbol{s}) (\theta(\boldsymbol{\xi}) - \theta(\boldsymbol{s}))^{\boldsymbol{\nu}-1} \mathrm{d}\boldsymbol{s} \\ &\leq \mathbb{K} + \frac{\mathbb{M} (\boldsymbol{\theta}(\boldsymbol{b}) - \boldsymbol{\theta}(\boldsymbol{a}))^{\boldsymbol{\nu}}}{\Gamma(\boldsymbol{\nu}+1)}. \end{split}$$

Hence, {*φn*} is uniformly bounded in *C*(J, R). With same argument one can deduce that {*n*} is uniformly bounded.

It is necessary to derive that the sequences {*φn*} and {*n*} are equi-continuous on J. To do this, choosing *ς*1, *ς*2 ∈ J, with *ς*1 ≤ *ς*2. By (5) and (6) and Lemma 2, we have

By the continuity of E*ν*,1 −*λ*(*ϑ*(*ς*) − *ϑ*(*a*))*ν* on J, the right-hand-side of the preceding inequality approaches zero, when *ς*1 → *ς*2. This implies that {*φn*(*ς*)} is equi-continuous on J. Likewise, we can demonstrate that {*n*(*ς*)} is equi-continuous. Hence, by using Ascoli-Arzelás theorem, the sequence *φnk* (*ς*) → *φ*∗(*ς*) and *nk* (*ς*) → ∗(*ς*) as *k* → ∞. Hence the aforesaid relation combined under the monotonicity of sequences {*φn*(*ς*)} and {*n*(*ς*)} yields

$$\lim\_{n \to \infty} \phi\_{\mathbb{H}}(\emptyset) = \phi^\*(\emptyset) \quad \text{and} \quad \lim\_{n \to \infty} \phi\_{\mathbb{H}}(t) = \phi^\*(\emptyset).$$

uniformly on *ς* ∈ J and the limit functions *φ*∗, ∗ satisfy (1).

**Step 3:** *φ*∗ and ∗ are maximal and minimal solutions of (1) in [*φ*0, 0].

Let *z* ∈ [*φ*0, 0] be any solution of (1). Suppose that

$$
\phi\_n(\emptyset) \le z(\emptyset) \le \phi\_n(\emptyset), \quad \emptyset \in \mathcal{J}, \text{ for some } n \in \mathbb{N}^\*. \tag{7}
$$

Setting *γ*(*ς*) = *z*(*ς*) − *φn*+1(*ς*). It follows that

$$\begin{aligned} \mathbb{E}^{\mathbb{V};\theta} \mathcal{D}\_{a^{+}}^{\mathbb{V};\theta} \gamma(\boldsymbol{\xi}) &= \mathbb{F} \{ \boldsymbol{\xi}, \boldsymbol{z}(\boldsymbol{\xi}) \} - \mathbb{F} \{ \boldsymbol{\xi}, \phi\_{\mathbb{V}}(\boldsymbol{\xi}) \} + \lambda \left( \phi\_{\mathbb{V}+1}(\boldsymbol{\xi}) - \phi\_{\mathbb{II}}(\boldsymbol{\xi}) \right) \\ &\geq -\lambda \left( \boldsymbol{z}(\boldsymbol{\xi}) - \phi\_{\mathbb{II}}(\boldsymbol{\xi}) \right) + \lambda \left( \phi\_{\mathbb{II}+1}(\boldsymbol{\xi}) - \phi\_{\mathbb{II}}(\boldsymbol{\xi}) \right) \\ &= -\lambda \gamma(\boldsymbol{\xi}) .\end{aligned}$$

Furthermore

$$\begin{aligned} \phi\_{n+1}(a) &= \phi\_n(a) - \frac{1}{c} \mathbb{G}(\phi\_n(a), \phi\_n(b)) \\ &= \phi\_n(a) + \frac{1}{c} \mathbb{G}(z(a), z(b)) - \frac{1}{c} \mathbb{G}(\phi\_n(a), \phi\_n(b)) \\ &\leq z(a) - \frac{d}{c} (z(b) - \phi\_n(b)) \\ &\leq z(a), \end{aligned}$$

that is

$$\gamma(a) \ge 0.$$

In view of Lemma 4, we obtain *γ*(*ς*) ≥ 0, *ς* ∈ J, which implies

$$
\phi\_{n+1}(\emptyset) \le z(\emptyset) \text{, } \emptyset \in \mathcal{J}.
$$

Utilizing the same procedure, we can derive that

$$z(\emptyset) \le \mathcal{O}\_{n+1}(\emptyset) \text{ . } \emptyset \in \mathcal{J}.$$

Hence,

$$
\phi\_{n+1}(\emptyset) \le z(\emptyset) \le \phi\_{n+1}(\emptyset), \emptyset \in \mathcal{J}.
$$

Hence (7) is satisfied on J for all *n* ∈ N∗. Employing *n* → ∞ on (7) from either side, one has

$$
\phi^\* \le z \le \phi^\*.
$$

This confirms that *φ*∗, ∗ are the extremal solutions of (1) in [*φ*0, 0].

### **4. Illustrative Problems**

This section includes some test problems for the illustration of our main results.

**Example 1.** *Consider the ϑ-Caputo FDE* (1) *with*

$$\nu = 0.5, \ a = 1, \ b = e, \ \theta(\emptyset) = \ln \mathfrak{g}. \tag{8}$$

*In order to justify that Theorem 1 is valid, we take*

$$\begin{cases} \mathbb{F}(\boldsymbol{\xi}, \boldsymbol{\phi}(\boldsymbol{\xi})) = 1 - \boldsymbol{\phi}^2(\boldsymbol{\xi}) + 2\sqrt{\ln \boldsymbol{\xi}}, & \text{for } \boldsymbol{\xi} \in [1, \boldsymbol{\varepsilon}], \\ \mathbb{G}(\boldsymbol{\phi}(1), \boldsymbol{\phi}(\boldsymbol{\varepsilon})) = \boldsymbol{\phi}(1) - 1. \end{cases} \tag{9}$$

*Without difficulty, we can infer that*

$$
\phi\_0(\emptyset) = 1, \quad \phi\_0(\emptyset) = 1 + \sqrt{\ln \varrho}.
$$

*are lower and upper solutions of* (1)*, respectively. It is apparent that φ*0(*ς*) ≤ 0(*ς*)*, for ς* ∈ [1,*e*]. *In addition, for φ*0(*ς*) ≤ x ≤ y ≤ 0(*ς*)*, we have*

$$\mathbb{F}(\mathfrak{g}, \mathbf{y}) - \mathbb{F}(\mathfrak{g}, \mathbf{x}) \ge -4(\mathbf{y} - \mathbf{x}), \quad \mathfrak{g} \in [1, e].$$

*Hence the Hypothesis 2 holds with λ* = 4*. Further, if φ*0(1) ≤ *u*1 ≤ *u*2 ≤ 0(1)*, φ*0(*e*) ≤ *v*1 ≤ *v*2 ≤ 0(*e*), *we have*

$$
\mathbb{G}(\mathfrak{u}\_2, \mathfrak{v}\_2) - \mathbb{G}(\mathfrak{u}\_1, \mathfrak{v}\_1) \le (\mathfrak{u}\_2 - \mathfrak{u}\_1).
$$

*Therefore, Hypothesis 3 holds with c* = 1 *and d* = 0*. An application of Theorem 1 shows that the problem* (1) *with the data* (8) *and* (9) *has extremal solutions in the region* [*φ*0, 0]. *Moreover, the monotone iterative sequences* {*φn*}*n*∈N, {*n*}*n*∈N *can be acquired by*

$$\begin{split} \phi\_{n+1}(\boldsymbol{\xi}) &= \mathbb{E}\_{0.5,1} \left( -4\sqrt{\ln\frac{\boldsymbol{\xi}}{\boldsymbol{s}}} \right) + \int\_{1}^{\boldsymbol{\xi}} \left( \ln\frac{\boldsymbol{\xi}}{\boldsymbol{s}} \right)^{-0.5} \mathbb{E}\_{0.5,0.5} \left( -4\sqrt{\ln\frac{\boldsymbol{\xi}}{\boldsymbol{s}}} \right) \\ &\times \left( 1 - \phi\_{n}^{2}(\boldsymbol{s}) + 2\sqrt{\ln\boldsymbol{s}} + 4\phi\_{n}(\boldsymbol{s}) \right) \frac{d\boldsymbol{s}}{\boldsymbol{s}} \;/\; n \geq 0, \end{split} \tag{10}$$

$$\begin{split} \boldsymbol{\omega}\_{n+1}(\boldsymbol{\varsigma}) &= \mathbb{E}\_{0.5,1} \left( -4\sqrt{\ln \frac{\boldsymbol{\xi}}{\boldsymbol{s}}} \right) + \int\_{1}^{\boldsymbol{\varsigma}} \left( \ln \frac{\boldsymbol{\xi}}{\boldsymbol{s}} \right)^{-0.5} \mathbb{E}\_{0.5,0.5} \left( -4\sqrt{\ln \frac{\boldsymbol{\xi}}{\boldsymbol{s}}} \right) \\ &\times \left( 1 - \boldsymbol{\alpha}\_{n}^{2}(\boldsymbol{s}) + 2\sqrt{\ln \boldsymbol{s}} + 4\boldsymbol{\alpha}\_{n}(\boldsymbol{s}) \right) \frac{\mathbf{d} \mathbf{s}}{\mathbf{s}}, \ n \geq 0. \end{split} \tag{11}$$

**Example 2.** *Consider the following Caputo FDE*

$$\begin{cases} ^c\mathbb{D}\_{0^+}^\nu\phi(\emptyset) = \sin(\phi(\emptyset)) - \phi(\emptyset), \ \emptyset \in [0, 1],\\ 0.5\phi(0) - 3\phi(0)\phi(1) = 0. \end{cases} \tag{12}$$

*Note that problem (12) is a particular case of problem* (1)*, where*

 $\nu = 0.5$ ,  $a = 0$ ,  $b = 1$ ,  $\theta(\emptyset) = \emptyset$ .

*and*

$$\begin{cases} \mathbb{F}(\mathfrak{g}, \mathfrak{\phi}(\mathfrak{\zeta})) = \sin(\mathfrak{\phi}(\mathfrak{\zeta})) - \mathfrak{\phi}(\mathfrak{\zeta}), & \mathfrak{\zeta} \in \mathfrak{J}, \\ \mathbb{G}(\mathfrak{\phi}(0), \mathfrak{\phi}(1)) = 0.5\mathfrak{\phi}(0) - 3\mathfrak{\phi}(0)\mathfrak{\phi}(1) \end{cases}$$

*Taking φ*0(*ς*) = 0 *and* 0(*ς*) = √*ς, it is easy to verify that φ*0, 0 *are lower and upper solutions of* (12)*, respectively, and φ*0(*ς*) ≤ 0(*ς*), *for ς* ∈ [0, 1]*. Therefore, Hypothesis 1 of Theorem 1 holds. However, if φ*0(*ς*) ≤ x ≤ y ≤ 0(*ς*) *we have*

$$\mathbb{F}(\mathfrak{g}, \mathbf{y}) - \mathbb{F}(\mathfrak{g}, \mathbf{x}) \ge -2(\mathbf{y} - \mathbf{x})\_{\prime} \quad \mathfrak{g} \in [0, 1].$$

*Hence Hypothesis 2 holds with λ* = 2*, and if φ*0(0) ≤ *u*1 ≤ *u*2 ≤ 0(0)*, φ*0(1) ≤ *v*1 ≤ *v*2 ≤ 0(1), *we have*

$$
\mathbb{G}(\mathfrak{u}\_2, \mathfrak{v}\_2) - \mathbb{G}(\mathfrak{u}\_1, \mathfrak{v}\_1) \le (\mathfrak{u}\_2 - \mathfrak{u}\_1).
$$

*Therefore, Hypothesis* 3 *holds with c* = 1 *and d* = 0*. According to Theorem 1, there exist monotone iterative sequences* {*φn*} *and* {*n*} *that are uniformly converging to φ*∗ *and* ∗*, respectively, and φ*∗, ∗ *are the extremal solutions in* [*φ*0, 0] *of problem* (12)*.*

**Example 3.** *Consider the following ϑ-Caputo FDE*

$$\begin{cases} \,^c \mathbb{D}\_{0^+}^{0,3\chi^2\varepsilon} \phi(\boldsymbol{\varsigma}) = \frac{2}{\sqrt{\pi}} \sqrt{\epsilon^{2\varepsilon} - 1} + \epsilon^{2\varepsilon} - 1 - \sin(\epsilon^{2\varepsilon} - 1) + \sin(\phi(\boldsymbol{\varsigma})) - \phi(\boldsymbol{\varsigma}), \; \boldsymbol{\varsigma} \in [0, 1],\\ \phi(0) = \frac{2\pi}{\epsilon^2 - 1 + 2\pi} \phi(1). \end{cases} \tag{13}$$

*Comparing the above problem with problem* (1)*, we obtain*

 $\nu = 0.5$ ,  $a = 0$ ,  $b = 1$ ,  $\theta(\emptyset) = e^{2\xi}$ ,

*and*

$$\begin{cases} \mathbb{F}(\boldsymbol{\xi}, \boldsymbol{\phi}(\boldsymbol{\xi})) = \frac{2}{\sqrt{\pi}} \sqrt{\boldsymbol{\varepsilon}^{2\underline{\boldsymbol{\varepsilon}}} - 1} + \boldsymbol{\varepsilon}^{2\underline{\boldsymbol{\varepsilon}}} - 1 - \sin(\boldsymbol{\varepsilon}^{2\underline{\boldsymbol{\varepsilon}}} - 1) + \sin(\boldsymbol{\phi}(\boldsymbol{\xi})) - \boldsymbol{\phi}(\boldsymbol{\xi}), & \boldsymbol{\xi} \in \mathcal{J}\_{\boldsymbol{\varepsilon}} \\ \mathbb{G}(\boldsymbol{\phi}(\boldsymbol{0}), \boldsymbol{\phi}(\boldsymbol{1})) = \boldsymbol{\phi}(\boldsymbol{0}) - \frac{2\pi}{\boldsymbol{\varepsilon}^{2} - 1 + 2\pi} \boldsymbol{\phi}(\boldsymbol{1}). & \end{cases}$$

*One can verify that φ*(*ς*) = *e*2*ς* − 1 *is an exact solution of the problem* (1)*. Moreover, taking φ*0(*ς*) = 0 *and* 0(*ς*) = *e*2*ς* − 1 + 2*π, it is easy to verify that φ*0, 0 *are lower and upper solutions of* (13)*, respectively, and φ*0(*ς*) ≤ 0(*ς*), *for ς* ∈ [0, 1]*. Therefore, Hypothesis 1 of Theorem 1 holds. However, if φ*0(*ς*) ≤ x ≤ y ≤ 0(*ς*) *we have*

$$\mathbb{F}(\mathfrak{g}, \mathbf{y}) - \mathbb{F}(\mathfrak{g}, \mathbf{x}) \ge -2(\mathbf{y} - \mathbf{x}), \quad \mathfrak{g} \in [0, 1].$$

*Hence Hypothesis 2 holds with λ* = 2*, and if φ*0(0) ≤ *u*1 ≤ *u*2 ≤ 0(0)*, φ*0(1) ≤ *v*1 ≤ *v*2 ≤ 0(1), *we have*

$$
\mathbb{G}(\mu\_2, v\_2) - \mathbb{G}(\mu\_1, v\_1) \le (\mu\_2 - \mu\_1) - \frac{2\pi}{\varepsilon^2 - 1 + 2\pi}(v\_2 - v\_1).
$$

*Therefore, Hypothesis 3 holds with c* = 1 *and d* = 2*π e*2−1+2*π . According to Theorem 1, there exist monotone iterative sequences* {*φn*} *and* {*n*} *that are uniformly converging to φ*∗ *and* ∗*, respectively, and φ*∗, ∗ *are the extremal solutions in* [*φ*0, 0] *of problem* (13)*.*

### **5. Conclusions**

We have established sufficient results by using monotone iterative techniques together with upper and lower solutions for a class of boundary value problem involving a generalized form of Caputo derivative of fractional order. By using the mentioned tool, we have established fruitful combinations between lower and upper solutions. Further, the said method has the ability to produce two sequences of upper and lower solutions, respectively. For the construction of the aforesaid sequences this method does not need any kind of discretization or collocation like other methods usually require. The two sequences we have generated are of a monotonic type with increasing and decreasing behaviors, respectively. Moreover, the sequence that is monotonically decreasing converges to its lower bound. However, the other one that is monotonically increasing is converging to its upper bound. The bounds for upper and lower solutions have also been investigated for their uniqueness using Banach contraction theorem. For the justification of our results, we have provided some examples. Overall we have concluded that the proposed procedure is a powerful and efficient tool to study various classes of FDEs for their extremal solutions. In future this technique can be applied to investigate those classes of FDEs involving non-singular-type derivatives under boundary conditions for upper and lower solution. In addition, the mentioned tool can be applied to investigate fractal-fractional-type problems corresponding to boundary conditions. Overall we have concluded that the monotone iterative technique of applied analysis is a powerful technique for dealing with various kinds of problems involving different types of fractional-order operators.

**Author Contributions:** Conceptualization, C.D., Z.B.; writing—original draft preparation, C.D., Z.B. and M.S.A.; methodology, C.D. and Z.B.; software, K.S. and B.A.; validation, C.D., Z.B., M.S.A. and T.A.; formal analysis, C.D., Z.B. and M.S.A.; writing—review and editing, C.D., Z.B., M.S.A., K.S. and T.A.; investigation, M.S.A., K.S., B.A. and T.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by the Prince Sultan University, Riyadh P.O. Box 11586, Saudi Arabia.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** Authors Kamal Shah, Bahaaeldin Abdalla and Thabet Abdeljawad would like to thank Prince Sultan University for paying the APC and support through TAS research lab.

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

### **References**

### *Article* **Some** q**-Fractional Estimates of Trapezoid like Inequalities Involving Raina's Function**

**Kamsing Nonlaopon 1, Muhammad Uzair Awan 2,\*, Muhammad Zakria Javed 2, Hüseyin Budak 3 and Muhammad Aslam Noor 4**

**Abstract:** In this paper, we derive two new identities involving q-Riemann-Liouville fractional integrals. Using these identities, as auxiliary results, we derive some new q-fractional estimates of trapezoidal-like inequalities, essentially using the class of generalized exponential convex functions.

**Keywords:** convex; exponential convex; fractional; quantum; inequalities

**MSC:** 05A30; 26A33; 26A51; 34A08; 26D07; 26D10; 26D15

### **1. Introduction**

A set C ⊆ R is said to be convex, if

$$(1 - \mathbf{v})\mathbf{x} + \mathbf{v}\mathbf{y} \in \mathbf{c}$$

for all x, y ∈ C and ν ∈ [0, 1]. A function F : C → R is said to be convex, if

$$
\beta^{\mathbb{P}}((1-\mathsf{v})\mathsf{x}+\mathsf{v}\mathsf{y}) \lesssim (1-\mathsf{v})\mathcal{F}(\mathsf{x}) + \mathsf{v}\mathcal{F}(\mathsf{y}).
$$

for all x, y ∈ C and ν ∈ [0, 1].

The classical concepts of convexity are simple but have many applications in different fields of pure and applied sciences. For example, they play a significant role in the theory of optimization, mathematical economics, operations research, etc. In recent years, the classical concepts of convexity have been extended and generalized in different directions using novel and innovative ideas. It has been observed that these new generalizations of classical convexity enjoy some nice properties which classical convexity has. Recently, Cortez et al. [1] presented a new generalization of convexity class as follows:

**Definition 1** ([1])**.** *Let* ρ, λ > 0 *and* σ = (σ(0), ... , σ(k), ...) *be a bounded sequence of positive real numbers. A non-empty set* I *is said to be generalized convex, if*

$$
\mathfrak{w}\_1 + \mathfrak{r} \mathcal{R}\_{\mathfrak{p}, \lambda, \sigma} (\mathfrak{w}\_2 - \mathfrak{w}\_1) \in \mathcal{I}
$$

*for all* 1, 2 ∈ I *and* τ ∈ [0, 1]*.*

**Citation:** Nonlaopon, K.; Awan, M.U.; Javed, M.Z.; Budak, H.; Noor, A.M. Some q-Fractional Estimates of Trapezoid like Inequalities Involving Raina's Function. *Fractal Fract.* **2022**, *6*, 185. https://doi.org/10.3390/ fractalfract6040185

Academic Editors: Asifa Tassaddiq and Muhammad Yaseen

Received: 24 February 2022 Accepted: 22 March 2022 Published: 25 March 2022

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

**Copyright:** © 2022 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 (https:// creativecommons.org/licenses/by/ 4.0/).

Here, Rρ,λ,σ(z) is the Raina's function and is defined as:

$$\mathcal{R}\_{\mathfrak{p},\lambda,\sigma}(z) = \mathcal{R}\_{\mathfrak{p},\lambda}^{\sigma(0),\sigma(1),\cdots}(z) = \sum\_{\mathbf{k}=0}^{\infty} \frac{\sigma(\mathbf{k})}{\Gamma(\mathfrak{p}\mathbf{k} + \lambda)} z^{\mathbf{k}},$$

where ρ, λ > 0, |z| < R, σ = {σ(0), σ(1), ... , σ(k), ...} is a bounded sequence of positive real numbers and Γ (η) = ∞ xη−1e−xdx is the gamma function. For details, see [2].

0 Cortez et al. [1] also defined the class of generalized convex functions as:

**Definition 2** ([1])**.** *Let* ρ, λ > 0 *and* σ = (σ(0), ... , σ(k), ...) *be a bounded sequence of positive real numbers. A function* F : I → R *is said to be generalized convex, if*

$$\mathcal{F}(\mathfrak{\boldsymbol{\sigma}}\_1 + \mathfrak{\boldsymbol{\pi}}\mathbb{R}\_{\mathfrak{\boldsymbol{\rho}},\lambda,\mathfrak{\boldsymbol{\sigma}}}(\mathfrak{\boldsymbol{\sigma}}\_2 - \mathfrak{\boldsymbol{\sigma}}\_1)) \lessapprox (1 - \mathfrak{\boldsymbol{\pi}})\mathcal{F}(\mathfrak{\boldsymbol{\sigma}}\_1) + \mathfrak{\boldsymbol{\pi}}\mathcal{F}(\mathfrak{\boldsymbol{\sigma}}\_2)$$

*for all* 1, 2 ∈ I *and* τ ∈ [0, 1]*.*

Awan et al. [3] introduced the class of exponential convex functions as:

**Definition 3** ([3])**.** *A function* F : C → R *is said to be exponentially convex, if*

$$\mathcal{F}((1-\mathbf{v})\mathbf{x}+\mathbf{v}\mathbf{y}) \lessapprox (1-\mathbf{v})\frac{\mathcal{F}(\mathbf{x})}{\exp(\alpha \mathbf{x})} + \mathbf{v}\frac{\mathcal{F}(\mathbf{y})}{\exp(\alpha \mathbf{y})},$$

*for all* x, y ∈ C *and* ν ∈ [0, 1]*.*

Besides its applications, the theory of convexity has also played a dynamic role in developing the theory of inequalities. A wide class of inequalities is just a direct consequence of the applications of the convexity property of the functions. Hermite–Hadamard's inequality, also known as trapezium-like inequality, is one of the most studied results. It reads as:

Let F : I ⊆ R → R be a convex function, then

$$\left(\mathcal{P}\left(\frac{\Box\_1 + \Box\_2}{2}\right) \leqslant \frac{1}{\Box\_2 - \Box\_1} \int\_{\overline{\Box}\_1}^{\overline{\Box}\_2} \mathcal{P}(\mathbf{x})d\mathbf{x} \leqslant \frac{\overline{\mathcal{P}}(\Box\_1) + \overline{\mathcal{P}}(\Box\_2)}{2}.$$

For some recent developments related to Hermite–Hadamard's inequality and its applications, see [4].

In recent years, several new techniques have been used to obtain new versions of Hermite–Hadamard's inequality. For instance, Sarikaya et al. [5] utilized the concepts of fractional calculus and obtained the fractional analogues of Hermite–Hadamard's inequality. Alp et al. [6] obtained quantum analogue of Hermite–Hadamard's inequality. Awan et al. [3] obtained a new refinement of Hermite–Hadamard's inequality using the class of exponentially convex functions. Cortez et al. [1] obtained Hermite–Hadamard's inequality using the class of generalized convex functions. Kunt and Aljasem [7] obtained fractional quantum versions of Hermite–Hadamard type of inequalities. Noor et al. [8] obtained some more quantum estimates for Hermite–Hadamard inequalities using the class of convex functions. Sudsutad [9] obtained various new quantum integral inequalities for convex functions. Zhang et al. [10] obtained a new generalized quantum-integral identity and obtained new q-integral inequalities via (α, m)-convexity property of the functions.