A new study on the existence and stability to a system of coupled higher-order nonlinear BVP of hybrid FDEs under the $p$-Laplacian operator

1.   Introduction

Fractional calculus is a branch of mathematics that contains integrals and derivatives of arbitrary orders (including rational, real, complex numbers). Nowadays, the main focus of this field has been changed from pure mathematical formulations to applications on a variety of mathematical models such as modeling of pantograph systems ^[1], different kinds of Langevin equations ^[2,3], memristor-based chaotic circuits ^[4], etc. Indeed, a new approach for modeling physical phenomena within a simple and effective procedure is used via fractional differential equations (FDEs) with singular and nonsingular derivatives. Therefore, in recent decades, fractional operators have been appeared in science, engineering, physics, dynamics, biological models, as well as electrodynamics, and fluid mechanics. For more details, interested readers are referred to see the heat equation in the context of the Caputo-Fabrizio derivative ^[5], fractional graph representation of ethane ^[6], Hahn-integro-difference equations ^[7], investigation of some diseases such as anthrax in animals and mumps via non-singular derivatives ^[8,9], fuzzy fractional model of coronavirus ^[10], fractal-fractional modeling ^[11], uncertain fractional modeling ^[12] and the references therein. The reason beyond the success of fractional calculus in modeling natural phenomena is that its operators are nonlocal, which strongly makes them suitable and efficient to describe the long memory or nonlocal effects characterizing most physical phenomena, where the models of problems with uncertain parameters reflect the actual behavior of the systems and make the model more realistic. Obtaining an analytical solution for FDEs is more difficult in most cases. Thus, many researchers have proposed and improved several numerical techniques to get approximate solutions. In addition, finding approximate asymptotic solutions for these problems allows researchers in control and similar fields to find stable solutions quickly. In ^[13], asymptotic interval approximate solutions for FDEs are discussed with the Atangana-Baleanu derivative. On the other hand, studying the existence and uniqueness (EUS) of solutions is a hot topic in FDEs. In the sequel, we present some recent scientific contributions of researchers about EUS to FDEs. For example, in ^[14], the authors conducted research regarding the solution set of the following FDE:

under some simple restrictions on the functional coefficient $b(\gamma)$. Baleanu et al. ^[15] discussed two global solutions in relation to an initial value problem (IVP) involving a vast category of FDEs. Zhao et al. ^[16] developed an existence theorem of solution for FHDEs (fractional hybrid differential equations) involving Riemann-Liouville differential under mixed Lipschitz and Caratheodory conditions. Sitho et al. ^[17] used fixed point theorems to examine existence results for initial value problems for HFDEs. Khan et al. ^[18] investigated the EUS for a coupled system of FDEs with $p$-Laplacian operator. In the last few decades, an important class of FDEs has taken great attention in nonlinear differential equations, known as HDEs. For example, in ^[19], the authors studied the EUS to the ordinary HDEs with linear perturbation of the first type given by

where $\mathcal{Z}\in C\big([\varsigma_0, \varsigma_0+b]\times\mathbb{R}, \mathbb{R}-\{0\}\big)$, $b\in\mathbb{R}^{+}$, $[\varsigma_0, \varsigma_0+b]$ is bounded interval, $\mathcal{Z}(\varsigma, \mathcal{Q}(\varsigma))$ is continuous and $\mathcal{Y}(\varsigma, \mathcal{Q}(\varsigma))$ is a Caratheodory class of functions. Dhage et al. ^[20] studied the EUS of the ordinary HDEs under the linear perturbation of the second type given by

Some fundamental differential inequalities are used to prove the existence results for extremal solutions. Herzallah et al. ^[21] proved the EUS for two first and second types of HFDEs with the aid of fixed points, given as

where $\mathcal{D}^{\beta}$ is the Caputo fractional derivative of order $0 < \beta < 1$ and $\tau\in[0, T]$. For more recent hybrid models, we refer to ^{[22,23,24,25]}.

Recently, FDEs with $p$-Laplacian operators have been investigated by many researchers. For example, Khan et al. ^[26] established the EUS of a class of nonlinear HFDEs under p-Laplacian operator by a method based on topological degree, given by

For further information, we refer to read these papers ^[27]. To the best of the authors' knowledge, no publication exists that deals with the EUS and Hyers-Ulam stability (HU-stability) of general HFDEs under $p$-Laplacian operators equipped with mixed boundary conditions. Motivated by this reason, our main goal in this manuscript is to prove the EUS as well as the HU-stability for the following general system of nonlinear HFDEs under $p$-Laplacian operator which is formulated as

where $^c{\mathcal{D}}_0^{\gamma_i}$, $^c{\mathcal{D}}_0^{\beta_i}$, $i = 1, 2,$ are the Caputo fractional derivatives with $m-1 < \gamma_i, \, \beta_i\leq m,$ and $m$ is nonnegative integer number. For $k\in\mathbb{N}_1^4$, $\psi_k$ is a continuous function and belongs to $\mathcal{L}[0, 1]$, $\phi_p(y) = |y|^{p-2}y$ is a $p$-Laplacian operator, where $\frac{1}{p}+\frac{1}{q} = 1$ and $\phi_q = \phi_p^{-1}$.

The structure of this article is organized as follows: We present several fundamental theorems, definitions, and lemmas in Section 2 to be used in our study. In Section 3, by using the Green functions, we transform the coupled hybrid BVPs (1.1) into integral equations, Then by defining a set of operators, the integral equations are converted into a new equivalent fixed point problem. Thereafter, by using the topological degree method, the main results for existence and uniqueness are proved. In Section 4, the stability of the proposed system (1.1) is investigated via the Hyers-Ulam criterion. In Section 5, we give an example to show the validity and efficacy of the results, and finally, the conclusion is presented in Section 6.

2.   Auxiliary results

In this section, we recall some definitions and lemmas.

Definition 2.1. ^[28] The Caputo fractional derivative of order $\gamma$ of a real continuous function $\mathcal{Q}$ on $[0, \infty)$ is defined by

where $m-1 < \gamma < m$, $m$ is an integer number, if the right-hand side integral converges point-wise on the interval $(0, \infty)$.

Definition 2.2. ^[29] The Riemann-Liouville fractional integral of order $\gamma$ of a function $\mathcal{Q}$ (with above property) is defined by

Lemma 2.1. ^[30] Let $\gamma \in (m, m-1]$, $\mathcal{Q}\in C^{m-1}$. Then

for $b_i\in \mathbb{R}$, $i\in \mathbb{N}_1^m$.

The Banach space $\mathcal{I} = C\big((0, 1], \mathbb{R}\big)$ is a family of continuous functions with the functional norm $\|\mathfrak{f}\| = \sup\{|\mathfrak{f}(\varsigma)|:\varsigma\in[0, 1]\}.$ If we consider the product space $\mathcal{J} = \mathcal{I}\times\mathcal{I}$ via the norm $\|(\mathfrak{f}, \nu)\| = \|\mathfrak{f}\|+\|\nu\|$, then it is also a Banach space. In the following, some basic notations and results linked to the nonlinear coincidence degree theorem are recalled, which one can find them in ^[31,32,33].

Definition 2.3. ^[34] Let the category of all bounded sets $\mathcal{C}(\mathcal{J})$ be denoted by $\mathcal{S}.$ Then the mapping $\mathcal{B}:\mathcal{S}\rightarrow (0, \infty)$ entitled the Kuratowski measure of non-compactness is defined by

Definition 2.4. A mapping $\mathcal{A}:\mathcal{J}\rightarrow \mathcal{J}\, \,$is said to be a contraction if there is $1 > \xi > 0$ such that for each pairs of the points $(\mathfrak{f}, \nu)$ and $(\mathfrak{f}^*, \nu^*)$, the distance between the images of these points under $\mathcal{A}\, \,$ are closer than the distance between the points. Mathematically, we mean

Lemma 2.2. ^[34] For the measure $\mathcal{B}$, we have the following properties:

(1) For the relatively compact set $K$, $\mathcal{B}(K) = 0;$

(2) $\mathcal{B}$ is seminorm, i.e., $\mathcal{B}(\mu K) = |\mu|\mathcal{B}(K)$, $\mu\in\mathbb{R}$, and $\mathcal{B}(K_1+K_2)\leq\mathcal{B}(K_1)+\mathcal{B}(K_2)$;

(3) $K_1\subset K_2$ yields $\mathcal{B}(K_1) \leq \mathcal{B}(K_2)$;

(4) $\mathcal{B}({{conv}}K) = \mathcal{B}(K)$;

(5) $\mathcal{B}(\bar{K}) = \mathcal{B}(K)$.

Definition 2.5. ^[34] Suppose that $\mathcal{A}:\mathcal{S}\rightarrow \mathcal{I}$ is bounded and continuous such that $\mathcal{S}\subset\mathcal{I}.$ Then $\mathcal{A}$ is an $\mathcal{B}$-Lipschitz, if there is $\epsilon\geq0$ such that

In addition, $\mathcal{A}$ is called a strict $\mathcal{B}$-contraction under the condition $\epsilon < 1.$

Definition 2.6. ^[34] The function $\mathcal{A}$ is $\mathcal{B}$-contraction if

Therefore, $\mathcal{B}(\mathcal{A}(K)) > \mathcal{B}(K)$ yields $\mathcal{B}(K) = 0.$

Note that $\epsilon < 1$ implies that $\mathcal{A}$ is a strict contraction.

Assume that $\Omega_1$ and $\Omega_2$ are two bounded sets in $\mathcal{I}$ such that $0\in\Omega_1$, $\overline{\Omega_1}\subset\Omega_2$, and $\mathcal{A}:\mathcal{I}\cap(\overline{\Omega_2}\backslash\Omega_1)\rightarrow \mathcal{I}$ is an operator.

Definition 2.7. ^[35] An operator $\mathcal{A}:\mathcal{I}\cap(\overline{\Omega_2}\backslash\Omega_1)\rightarrow \mathcal{I}$ is

● ($\mathcal{U}_1$) Uniformly bounded if there is $\mathfrak{E} > 0$ such that $|\mathcal{A}(\varsigma)|\leq \mathfrak{E}$ for all $t\in \mathcal{I}\cap(\overline{\Omega_2}\backslash\Omega_1)$;

● ($\mathcal{U}_2$) Equicontinuous if for every $\epsilon > 0$, there is $\delta(\epsilon) > 0$ such that $|\mathcal{A}(\varsigma_1)-\mathcal{A}(\varsigma_2)| < \epsilon$ for all $\varsigma_1, \varsigma_2 \in\mathcal{I}\cap(\overline{\Omega_2}\backslash\Omega_1)$ with $|\varsigma_1-\varsigma_2| < \delta$.

Theorem 2.1. ^[36] Let $\mathcal{A}:\mathcal{I}\cap(\overline{\Omega_2}\backslash\Omega_1)\rightarrow \mathcal{I}$. Then $\mathcal{A}$ is compact if and only if it is uniformly bounded and equicontinuous.

Theorem 2.2. ^[37] Let $\mathcal{I}$ be a Banach space. Then a contraction $\mathcal{A}:\mathcal{I}\rightarrow \mathcal{I}$ has a unique solution $\varpi$ such that $\mathcal{A}(\varpi) = \varpi$ (the Banach contraction principle).

Proposition 2.1. ^[38] $\mathcal{A}$ is $\mathcal{B}$-Lipschitz with constant $\epsilon > 0$ if and only if $\mathcal{A}:\mathcal{S}\rightarrow \mathcal{I}$ is Lipschitz with constant $\epsilon > 0$.

Proposition 2.2. ^[38] The mapping $\mathcal{A}$ is $\mathcal{B}$-Lipschiz with constant $\epsilon = 0$ if and only if $\mathcal{A}:\mathcal{S}\rightarrow \mathcal{I}$ is said to be compact.

Theorem 2.3. ^[18] Let $\mathcal{A}:\mathcal{I}\rightarrow \mathcal{I}$ be a $\mathcal{B}$-contraction and

If $\mathfrak{D} \subset \mathcal{I}$is bounded, i.e., there is some $r > 0$ with $\mathfrak{D}\subset x_r(0)$ and

then $\mathcal{A}\, \,$ has at least one fixed point in $x_r(0).$

Lemma 2.3. ^[39,40] Consider the $p$-Laplacian operator $\phi_p$. Then

(1) If $1 < p\leq2$, $\mathcal{M}_1, \mathcal{M}_2 > 0$, and $|\mathcal{M}_1|, \mathcal{M}_2\geq\eta > o$, then

(2) If $p > 2$, and $|\mathcal{M}_1|, |\mathcal{M}_2|\leq\eta$, then

3.   Results regarding the EUS

In this section, based on the above auxiliary notions, we prove the following results on the EUS of the considered system (1.1).

Theorem 3.1. Assume that $\psi\in C[0, 1]\,$ is an integrable real function. Then for $\beta_{1}, \gamma_{1}\in(0, 1]$, the solution of

is

where $\mathcal{G}^{\gamma_1}(\varsigma, \eta)$ and $\mathcal{G}^{\beta_1}(\eta, \varsigma)\,$ are Green functions defined by

and

Proof. By applying the integral $\mathcal{I}_0^{\beta_1}$ on both sides of (3.1) and using Lemma 2.1, we get

By using the conditions $\big[\phi_p\big({^c{\mathcal{D}}_0^{\gamma_1}}\big(f(\varsigma)-\psi_2(\varsigma, \nu(\varsigma))\big)\big)\big]_{\varsigma = 0}^{(i)} = 0$, $i = 0, 2, 3, \ldots, m-1$, we get $c_1 =$ $c_3 = c_4 = \ldots = c_m = 0$.

The conditon $\big[\phi_p\big(^c{\mathcal{D}}_0^{\gamma_1}\big(f(\varsigma)-\psi_2(\varsigma, \nu(\varsigma))\big)\big)\big]_{\varsigma = \lambda}^{'} = 0$, $i\in\mathbb{N}_0^{m-1}-\{1\}$ gives

By inserting $c_i$, $i\in\mathbb{N}_1^m$ in (3.5), we get

We write (3.7) as

By applying $\phi_q = \phi_p^{-1}$ on both sides of (3.8), we have

By integrating on both sides of (3.9) by employing the operator $\mathcal{I}_0^{\gamma_1}$, we obtain

By using the condition ${f^{(i)}(\varsigma)}_{\varsigma = 0} = 0$, $i\in\mathbb{N}_1^{m-2}$, we get $k_2 = k_3 = \ldots = k_{m-1} = 0$.

By inserting the values of $k_i$, $i\in\mathbb{N}_2^{m-1}$ in (3.10), we obtain

Since $f^{(m-1)}(\varsigma)_{\varsigma = 1} = 0$, hence

Substituting $k_m$ in (3.11), we get

Also, the last condition $f(1)-\frac{1}{(m-1)!}f^{(m-1)}(0) = 0$ gives

Substituting $k_1$ in (3.13), we get

By using the Green function, we can write (3.15) as

This completes the proof.

By applying Theorem 3.1, the solution of (1.1) is corresponding to the following coupled integral equations:

where $\mathcal{G}^{\gamma_2}(\varsigma, \eta)$ and $\mathcal{G}^{\beta_2}(\varsigma, \eta)$ are the following Green functions:

and

Now, we define the operators $\mathcal{T}_j:\mathcal{I}\rightarrow \mathcal{I}$ for $j\in\mathbb{N}_1^4$ as follows:

Depending on $\mathcal{T}_j$, $j\in\mathbb{N}_1^4$, we define $\mathfrak{B}_1(f, \nu) = (\mathcal{T}_1, \mathcal{T}_2)$, $\mathfrak{B}_2(f, \nu) = (\mathcal{T}_3, \mathcal{T}_4)$, and then $\mathfrak{B}(f, \nu) = \mathfrak{B}_2(f, \nu)+\mathfrak{B}_1(f, \nu)$. So, by Theorem 3.1, the solutions of (3.17) and (3.18) are equivalent to the following fixed point problem:

In order to prove our next results, we introduce the following assumptions:

$(\mathfrak{R}_1)$ For positive constants $a_1, a_2, \mathbb{W}_{\psi_1}, \mathbb{W}_{\psi_3}$ and $d_1, d_2\in[0, 1]$, the two functions $\psi_1$, $\psi_3\,$ satisfy the following conditions:

$(\mathfrak{R}_2)$ For positive constants $a_3, a_4, \mathbb{W}_{\psi_2}, \mathbb{W}_{\psi_4}\,$ and $d_1, d_2\in[0, 1],$ the two functions $\psi_2, \; \psi_4\,$ satisfy the following conditions:

$(\mathfrak{R}_3)$ There are $\epsilon_{\psi_1}$ and $\epsilon_{\psi_2}\,$ such that for each $f, \nu, x, r\in\mathcal{B}$,

$(\mathfrak{R}_4)$ There are $\epsilon_{\psi_2}$ and $\epsilon_{\psi_4}\,$ such that for each $f, \nu, x, r\in\mathcal{B}$,

Remark 3.1. For the simplicity, we define the following notations:

In addition to that,

Theorem 3.2. Under the assumptions $(\mathfrak{R}_1)$–$(\mathfrak{R}_4)$, the operator $\mathfrak{B}:\mathcal{J}\rightarrow \mathcal{J}$ is continuous, and the growth condition holds under $\mathfrak{B}$ as follows:

for each $(f, \nu)\in \mathfrak{g}_r\subset\mathcal{J},$ and $\|(f, \nu)\|^d = \|f\|^d+\|\nu\|^d$ for $d = \max\{d_1, d_2\},$ where $\mathcal{K}$ and $\mathbb{M}^*$ are defined in Remark 3.1.

Proof. We consider a bounded set $\mathfrak{g}_r = \{(f, \nu)\in \mathcal{J}:\|(f, \nu)\|\leq r\}$ along with the sequence $(f_n, \nu_n)\,$ tending to $(f, \nu)$ in $\mathfrak{g}_r$. To prove $\|\mathfrak{B}(f_n, \nu_n)-\mathfrak{B}(f, \nu)\|\rightarrow0$, as $n\rightarrow \infty$, we first show the continuity of $\mathfrak{B}_1$. So, we estimate

This implies

For $\mathcal{T}_2\nu_n$ and $\mathcal{T}_2\nu$, we use the same steps to get

Since $f$ and $\nu$ are continuous, therefore $|\nu_n(\xi)-\nu(\xi)|\rightarrow0$, and $|f_n(\xi)-f(\xi)|\rightarrow0$ as $n\rightarrow \infty$. By using (3.22) and (3.23), we conclude that $|\mathcal{T}_1f_n(\varsigma)-\mathcal{T}_1f(\varsigma)|\rightarrow0$ and $|\mathcal{T}_2\nu_n(\varsigma)-\mathcal{T}_2\nu(\varsigma)|\rightarrow0$ as $n\rightarrow \infty$. This implies that $\mathcal{T}_1$ and $\mathcal{T}_2$ are continuous. The continuity of $\mathcal{T}_1$ and $\mathcal{T}_2$ leads to that of $\mathfrak{B}_1 = (\mathcal{T}_1, \mathcal{T}_2)$.

Secondly, to show the continuity of $\mathfrak{B}_2$, we consider

and

Since $f$ and $\nu$ are continuous, therefore, $|\nu_n(\xi)-\nu(\xi)|\rightarrow0$ and $|f_n(\xi)-f(\xi)|\rightarrow0$ as $n\rightarrow \infty$. By using (3.24) and (3.25), we conclude that $|\mathcal{T}_3f_n(\varsigma)-\mathcal{T}_3f(\varsigma)|\rightarrow0$ and $|\mathcal{T}_4\nu_n(\varsigma)-\mathcal{T}_4\nu(\varsigma)|\rightarrow0$ as $n\rightarrow \infty$. This gives that $\mathcal{T}_3$ and $\mathcal{T}_4$ are continuous. The continuity of $\mathcal{T}_3$ and $\mathcal{T}_4$ gives that of the operator $\mathfrak{B}_2 = (\mathcal{T}_3, \mathcal{T}_4)$. Consequently, $\mathfrak{B}(f, \nu) = (\mathfrak{B}_1(f, \nu), \mathfrak{B}_2(f, \nu))$ is continuous.

In order to prove the inequality (3.21), we use $\mathcal{T}_1, \; \mathcal{T}_2$ and the assumption $(\mathfrak{R}_1)$ as follows:

This leads to

For $\mathcal{T}_2$, by using the same steps, we obtain the following inequality:

By (3.26) and (3.27), we have

This implies

Next, by using $\mathcal{T}_3$ and $\mathcal{T}_4$ and the assumption $(\mathfrak{R}_2)$, we get

By (3.29) and (3.30), we have

Thus,

With the aid of (3.28) and (3.31), we get

This gives us the following inequality:

The proof is completed.

Theorem 3.3. With assumption $(\mathfrak{R}_1)$, the operator $\mathfrak{B}_1:\mathcal{J}\rightarrow \mathcal{J}$ is compact. Moreover, $\mathfrak{B}_1$ is $\mathcal{B}$-Lipschitz with constant zero.

Proof. Based on Theorem 3.2, we conclude that the operator $\mathfrak{B}_1:\mathcal{J}\rightarrow \mathcal{J}$ is bounded. Next, with the help of the assumption $(\mathfrak{R}_1)$ and by choosing $\varsigma_1, \varsigma_2\in[0, 1]$ arbitrarily, we have

This implies

Using the same steps with $\mathcal{T}_2\nu$, we get

By using (3.32) and (3.33), we have

As $\varsigma_1\rightarrow \varsigma_2$, the right-hand side of (3.34) tends to zero. Thus, $\mathfrak{B}_1 = \mathcal{T}_1f(\varsigma)+\mathcal{T}_2\nu(\varsigma)$ is an equicontinuous operator on $\mathcal{J}$. By Theorem 2.1, the operator $\mathfrak{B}_1(\mathcal{J})$ is compact. Therefore, $\mathfrak{B}_1$ is $\mathcal{B}$-Lipschitz with constant zero (by Proposition 2.2).

Theorem 3.4. Under the assumptions $(\mathfrak{R}_1)$ and $(\mathfrak{R}_3)$, and with the condition $\mathcal{K} = (\mathcal{Y}_1^*+\mathcal{Y}_2^*) < 1$, the system of coupled $p$-Laplacian HFDEs $(1.1)$ has a solution, and the set of solution is bounded in $\mathcal{J}$.

Proof. In order to prove the existence of the solution to our problem (1.1), Theorem 2.3 is used. Let

We consider a contrary procedure to demonstrate that $\mathcal{X}$ is bounded. For some $(f, \nu)\in\mathcal{J}$, we have $\|(f, \nu)\| = \mathbb{S}\rightarrow \infty$. From Theorem 3.2, we estimate

Since $\|(f, \nu)\| = \mathbb{S}$, then (3.36) implies

Divide both sides of (3.37) by $\|(f, \nu)\|$. Thus,

Therefore,

This leads to that $1\leq 0$ as $\mathbb{S}\rightarrow \infty$, but this is a contradiction of our assumption. Ultimately $\|(f, \nu)\| < \infty$, and hence $\mathcal{X}$ is bounded set and by Theorem 2.3, $\mathfrak{B}$ involves at least one fixed point which is the solution to supposed system of coupled $p$-Laplacian HFDEs (1.1), and the set of such solutions is bounded in $\mathcal{J}.$

Theorem 3.5. Let assumptions $(\mathfrak{R}_3)$ and $(\mathfrak{R}_4)$ to be held. Then the system of coupled $p$-Laplacian HFDEs $(1.1)$ has a unique solution if $\Delta = \mathcal{Z}_1+\mathcal{Z}_2 < 1.$

Proof. Firstly, by using $\mathcal{T}_1, \; \mathcal{T}_2$ and assumptions $(\mathfrak{R}_3), \; (\mathfrak{R}_4)$, we have

This implies

Using the same steps with $\mathcal{T}_2$, we get

By (3.38) and (3.39), we get

This leads to

Secondly, using $\mathcal{T}_3, \; \mathcal{T}_4$ and assumptions $(\mathfrak{R}_3)$, $(\mathfrak{R}_4)$, we get

By (3.41) and (3.42), we have

This implies

By using (3.40) and (3.43), we obtain

Thus,

By using the Banach's theorem (Theorem 2.2), and by considering $\Delta = \mathcal{Z}_1+\mathcal{Z}_2 < 1$, the inequality (3.44) implies that the contraction $\mathfrak{B}$ has a unique fixed point. It means that the system of coupled $p$-Laplacian HFDEs (1.1) has a unique solution.

4.   Hyers-Ulam Stability

In this part, we try to investigate HU-stability for the system of coupled $p$-Laplacian HFDEs (1.1). We refer that stability can also be studied with Lyapunov's direct method, which is considered as an approach for analyzing a stable system without solving the FDEs explicitly, see ^[41]. We present the following definition of HU-stability in the light of the definition offered in ^[42,43,44].

Definition 4.1. A system of coupled Hammerstein-type integral equations $(3.17)$ and $(3.18)$ is HU-stable if there are $M_f, M_\nu > 0$ such that for each $\mu_1, \mu_2 > 0$, and $(f, \nu)\in \mathcal{J}$ satisfying the following conditions:

If

then, there exists a pair, say $(\bar{f}(\varsigma), \bar{\nu}(\varsigma))\in \mathcal{J}$ satisfying

such that $|f(\varsigma)-\bar{f}(\varsigma)|\leq M_f\mu_2,$ and $|\nu(\varsigma)-\bar{\nu}(\varsigma)|\leq M_\nu\mu_1.$

Theorem 4.1. With the assumptions $(\mathfrak{R}_3)$ and $(\mathfrak{R}_4),$ the solution of the system of coupled $p$-Laplacian HFDEs $(1.1)$ is HU-stable.

Proof. By Theorem 3.5 and Definition 4.1, let $(f(\varsigma), \nu(\varsigma))$ satisfies the system of coupled integral equations (3.17) and (3.18). Let the pair $(\bar{f}(\varsigma), \bar{\nu}(\varsigma))$be any other approximation satisfying (4.3) and (4.4). In this case, we estimate

This implies

Using the same processes with $\nu(\varsigma)$ and $\bar{\nu}(\varsigma)$, we get

Therefore, with the aid of (4.5), (4.6) and assuming $M_f = (\epsilon_{\psi_2}+\mathcal{Y}_3)$, $M_\nu = (\epsilon_{\psi_4}+\mathcal{Y}_4)$, the coupled system (3.17), (3.18) is HU-stable. According to that, the system of coupled $p$-Laplacian HFDEs (1.1) is HU-stable.

5.   Example

An applied example of our system of coupled $p$-Laplacian HFDEs (1.1) will be addressed in this section to examine our results.

Example 5.1. The following example is provided to highlight our theoretical results for higher-order HFDE system with $p$-Laplacian operator when $m = 3$ given by

where $\nu\in [0, 1]$, $a_1 = a_2 = d_1 = d_2 = 0.3$, $p = 4$, $q = \frac{4}{3}$, $\lambda = 0.5$, $\gamma_1 = \frac{5}{2}$, $\gamma_2 = \frac{8}{3}$, $\beta_1 = \frac{9}{4}$, $\beta_2 = \frac{11}{4}$, $\eta_1 = \eta_2 = 1$, $\psi_1 = \frac{5}{16}\varsigma +\sin(\nu(\varsigma))$, $\psi_2 = \sqrt{\varsigma}(\frac{17}{26}+5 \cos(\nu(\varsigma)))$, $\psi_3 = \frac{4}{13}\varsigma^2+\cos(f(\varsigma))$, $\psi_4 = \varsigma^{\frac{1}{3}}(\frac{5}{16}\varsigma+\sin(f(\varsigma)))$, $\epsilon_{\psi_1} = \epsilon_{\psi_3} = \frac{3}{8}$, and $\epsilon_{\psi_2} = \epsilon_{\psi_4} = \frac{1}{8}$. By simple calculations, we get $\mathcal{Y}_3 = 0.111259$, $\mathcal{Y}_4 = 0.0540563$. So, $\mathcal{Z}_1 = 0.1653153$, $\mathcal{Z}_2 = 0.25$. Then, we have $\Delta = \mathcal{Z}_1+\mathcal{Z}_2 = 0.4153153 < 1$. By Theorem 3.5, we deduce that the higher-order HFDE system with $p$-Laplacian operator (5.1) has a unique solution. With similar fashion, the satisfication of the conditions of Theorem 4.1 can be easily established and according to that fact, the system of coupled $p$-Laplacian HFDEs (5.1) is HU-stable.

6.   Conclusions

As a result of the Banach contraction principle, Arzela-Ascoli's theorem, and nonlinear functional analysis, our work has provided suitable conditions for the existence and uniqueness of solution to the higher-order nonlinear boundary value problem of HFDEs which is more general and complex than many nonlinear problems in the literature. Moreover, the existence and uniqueness results were proved by using the topological degree method. The stability of the proposed system has been studied in the sense of Hyers-Ulam criterion. To validate our results, we provided an illustrative example. We also suggest for the researchers that the problem (1.1) has the potential to be studied for further aims, including multiplicity results and generalizing it with nonsingular operators.

Acknowledgments

This research was funded by National Science, Research and Innovation Fund (NSRF), and Suan Dusit University (contract No. 65-FF-010). Also, the sixth and seventh authors would like to thank Azarbaijan Shahid Madani University.

Conflict of interest

The authors declare no conflicts of interest.