On implicit coupled systems of fuzzy fractional delay di ﬀ erential equations with triangular fuzzy functions

: In this paper, we introduce and study an implicit coupled system of fuzzy fractional delay di ﬀ erential equations involving fuzzy initial values and fuzzy source functions of triangular type. We assume that these initial values and source functions are triangular fuzzy functions and deﬁne the solution of the implicit coupled system as a triangular fuzzy function matrix consisting of real functional matrices. The method of triangular fuzzy function, fractional steps and fuzzy terms separation are used to solve the implicit coupled systems. Further, we prove existence and uniqueness of solution for the considered systems, and also construct a solution algorithm. Finally, an example is given to illustrate our main results and some further work are presented.


Introduction
In an ecosystem, there exist universally competitions between two species, which is very important. And they tend to interact with each other so that there are some relationships such as predation or symbiosis among them (i.e. coupled systems) (see [9,39]). In addition, the delay and ambiguity may appear because of different reasons such as physical properties of equipment used in the system, signal transmission or measurement of system variables (see [14]), so we need to take them into account. Thus, in this paper, motivated by the work of Dong et al. [13] and Fatullayev et al. [14], we come up with the following problem of fuzzy fractional implicit delay differential coupled systems (FFIDDCs): and explore existence and uniqueness result of solution for the FFIDDC (1.1), where c D α f (t) denotes α-order Caputo fractional derivative of function f (t) ∈ C 1 α for 0 < α < 1 (see [22]), τ is the value of time delay, Φ 1 (t) and Φ 2 (t) are triangular fuzzy functions (TFFs) defined on [−τ, 0], F 1 (t, y(t), c D α x(t)) and F 2 (t, x(t), c D β y(t)) are also TFFs defined on (0, ∞)×R×R, and m 1 (t) and m 2 (t) are continuous crisp functions. The concept of TFFs used in this paper was first proposed by Gasilov et al. [17], then it has been widely employed [7,16,27]. Unlike in the past, Gasilov et al. [17] thought that a fuzzy function is a fuzzy bunch of real functions rather than a fuzzy number-valued function. Each of these real functions has a certain membership degree. The new concept of TFFs broadened the way (see [14]), and so we deal with fuzzy functions and fuzzy problems. Remark 1.1. In this paper, we consider the problem with Caputo fractional derivative, which is based on Riemann-Liouville fractional derivative [22]. We note that for temporal variable, the Caputo fractional derivative has been widely recognized in real application. For example, the Caputo fractional derivative of a constant is 0, but in the case of finite number at the lower limit of the interval, the Riemann-Liouville fractional derivative of a constant is not equal to 0. And the initial conditions given by Riemann-Liouville method cannot be explained physically. For more related work, see [28,34]. However, can one discuss the existence and uniqueness results for the solutions of the FFIDDC (1.1) with Riemann-Liouville fractional derivative or Hadamard fractional derivative? These are usually presented as important and significative problems in future research.
Some special cases of the FFIDDC (1.1) are listed as follows: (i) If τ = 0, 0 ≤ t ≤ 1 and α j = 1, where α j is the membership degrees (MDs) of F j for j = 1, 2, the FFIDDC (1.1) can be rewritten as the problem investigated by Dong et al. [13]: Making full use of fixed point theorem and the vector form of Gronwall inequality, Dong et al. [13] gained existence and uniqueness of solution to the initial value problem (IVP), and discussed the estimation of the solution.
(ii) If α = β = 1, the FFIDDC (1.1) reduces to the following fuzzy implicit delay differential coupled systems: which is brand new and is worth studying in the future. We note that the system (1.2) is an extension of the linear inhomogeneous fuzzy delay differential equations (FDDEs) considered by Fatullayev et al. [14]: The authors presented a method to represent the solution of FDDE as a fuzzy set of real functions, and proved existence and uniqueness of solution for FDDE involving TFFs. Furthermore, Fatullayev et al. [14] showed clearly that the proposed method can be extended to the system of FDDEs by using the research results of Gasilov and Amrahov [15]. The purpose of this paper is to discuss existence and uniqueness of solution for the FFIDDC (1.1).
Comparing with the simpler integer-order differential equations, fractional differential equations (FDEs) have more wide applications and play an important role in engineering, physics, finance and other fields. Nowadays, FDEs have been used to more accurately describe the dynamics of many systems when several complex phenomena in numerous seemingly diverse and widespread fields of science and engineering is modeled and investigated. See, for example, [3,6,22,28,34] and references therein. In addition, to study ecosystem problems in biology, one often discovers that an ecosystem does not have only a single species. That is to say, there must be multiple species, so there exists competition among them and coupling is taken into account (see [9,39]). Coupling relationship refers to the interaction and mutual influence between two or more objects. In fact, the coupled systems of FDEs have been studied by multitudinous researchers. See, for example, [2,38] and the references cited therein. Recently, by applying standard fixed point theorems for multivalued maps, existence of solutions for coupled systems of fractional differential inclusions with coupled boundary conditions was obtained by Ahmad et al. [5].
In science and engineering and other practical applications, time delay usually occurs on account of manual measurement, signal transmission, aging of equipments and so on (see [30]). Thereby, in recent years, the scholars have been interested in solving the fractional differential equations with time delay (DFDEs) (see, for instance, [30,35] and the references therein). Brzdek and Eghbali [10] considered the Ulam's stability of DFDEs and also proved that under some appropriate assumptions, each approximate solution of DFDEs is close to its only exact solution. In order to acquire existence of solutions for a time dependent delay differential equation with constant delay, Tabassum et al. [33] extended the previous approaches about fuzzy weakly contraction mapping principle. Of course, there are not just one type of FDEs, whether they are linear or nonlinear. And the common ones are Caputo type, Riemann-Liouville type and Hadamard type (see [11,21,36]). Moreover, the delay is not necessarily finite, it can be also infinite. Some researchers dedicated themselves to the study of FDEs with infinite delay. Banach fixed point theorem and the nonlinear alternative of Leray-Schauder type initiated by Granas and Dugundji [18], had been utilized by Benchohra et al. [8] to investigate existence results of FDEs with infinite delay.
Certainly, in actual application, there are not only time delay, but also exists the possibility of ambiguity. Thus, one needs to take the fuzzy into account, and then the fuzzy fractional differential equations (FFDEs) are proposed. Agarwal et al. [1] put forward a notion of fuzzy Riemann-Liouville differentiability based on Hukuhara differentiability, which can be employed to solve the IVPs of FFDEs. After that, many investigators have extended the concepts and expression forms of FFDEs, and bent themselves to obtain existence and uniqueness of solution for FFDEs. These theories are also well applied in real life. See, for example, [24,25,29] and the references therein. As Fatullayev et al. [14] pointed out, "the types of fuzzy items are more important than their numbers when a fuzzy problem needs to be solved". In fact, we can easily deal with n-th order homogeneous linear differential equations with fuzzy initial values (i.e. the problem with n fuzzy inputs) by using Zadeh's extension principle. But if a fuzzy forcing function appears in the problem, even for the delay differential equation of first order (i.e. with 2 fuzzy inputs), it will be indistinct how to use the principle.
It is well know that explicit equations are special cases of implicit equations, and the implicit equations are more general (see [12]). Cubiotti and Yao [12] introduced and studied a class of implicit second-order ordinary differential equations with known endpoint values. They also proved existence results under two conditions. However, as Hoa and Vu [21] stated briefly, "to the best of our knowledge, there exists no literature devoted to the uncertain fractional implicit differential equation in the fuzzy setting with the concepts of Riemann-Liouville, Caputo, Hadamard fractional derivatives". And then, Hoa and Vu [21] proposed some existence results for solutions of the fuzzy fractional implicit differential equations associated with the fuzzy Caputo-type, fuzzy Riemann-Liouville-type and fuzzy Hadamard-type concepts of fractional derivative, respectively. Very recently, Son and Dong [31] proposed global existence and some properties of solutions for the nonlocal problem of implicit fuzzy fractional differential systems. For more work on various forms of implicit FDE models, one can refer to [19,20,26] and the references therein.
The major methods of this paper are as follows: (i) We generalize TFF (Definition 2.1 of [17]) to a group of TFFs (see Example 2.1).
(ii) To solve the FFIDDC (1.1), we use the method of steps proposed by Fatullayev et al. [14]. Firstly, we figure out the solution of the FFIDDC (1.1) on interval [0, τ]. Afterwards, we similarly find the solution on interval [τ, 2τ]. And as we go along, eventually, the solution on interval [0, +∞) can be obtained.
(iii) Because of the fuzzy coupled fractional initial value problem with delays (FCFDIVP) (3.1), which is the identical transformation of the FFIDDC (1.1), is linear, the superposition principle can be used to work out it. As a consequence, we can decompose the FCFDIVP (3.1) into three subproblems (3.4)-(3.6) and solve them separately.
The remainder of this paper is organized as follows. The necessary preliminaries about the fuzzy theory that we are going to use are listed in Section 2. In Section 3, we introduce the concept of solution for the FFIDDC (1.1) and describe how to get it. Ultimately, a solution algorithm is proposed. In Section 4, we present an example to show that the hypotheses in all theorems can be met. Finally, in Section 5, we summarize the research results of this paper, and provide the content that one can study in the future.

Preliminaries
Based on the work of Zadeh [37], we define a fuzzy set A as a pair of the universal set U and the membership function (MF) µ : U → [0, 1]. The MF of a fuzzy set A can be denoted as µ A . For each After that, let U be the set of real numbers R, and a, c and b be real numbers which meet a ≤ c ≤ b. Then the set u with MF is called a triangular fuzzy number (TFN) and we denote it as u = (a, c, b). On the grounds of the geometric interpretation, the number c is called the vertex of u, we denote u = a and u = b to represent the left and the right end-points of u, respectively. Frequently, we express u = (a, c, b) as u = u cr + u un .
Here, u cr = c is the crisp part and u un = (a − c, 0, b − c) is the uncertain part of u. It is also useful to represent the fuzzy sets though their α-cuts. For each α ∈ (0, 1], the crisp set For the TFN u = (a, c, b), the α-cuts are intervals u α = [u α , u α ], where u α = a + α(c − a) and . These formulas can be rewritten as . From here we can see that an α-cut is homothetic to [a, b] (which is the 0-cut) with center c and with ratio (1 − α).
There are different notions about the fuzzy functions. In this study, we use the concept of the fuzzy function which was brought forward by Gasilov et al. [17], namely, fuzzy function is a bunch of fuzzy real functions. As a value F(t) of a fuzzy bunch F at time t, we understand the fuzzy set, which elements are the values of the real functions at t, with the higher MD of the corresponding functions. Mathematically, where "∧" and "→" are the logical conjunction and implication symbols, respectively. Definition 2.1. ( [17]) Let U be set of continuous functions defined on an interval I, and F a (·), F c (·), F b (·) ∈ U. We call the fuzzy subset F of U, determined by the MF as follows: as TFF and denote it as F = F a , F c , F b .
According to this definition, a TFF is a fuzzy set (or, fuzzy bunch) of real functions. Among them only two functions have the MD α: the functions y 1 = F a + α(F c − F a ) and Referring to the example of TFF given by Gasilov et al. [16], we give the following example of a group of TFFs.
Example 2.1. In Figure 1, we depict a group of TFFs as The value of a TFF at a time t ∈ I can be expressed by the following formula: We can easily find out that this value is a TFN, and a TFF F = F a , F b , F c is not a fuzzy numbervalued function. Actually, it is a fuzzy subset of the universe of continuous functions. Each element of this fuzzy subset is a real function with a certain MD.
. Further, when a TFF F is not regular, we call it as non-regular TFF. It means that for a non-regular TFF, in general, . And the graphs of functions F a , F c and F b can interchanged as t goes. In Figure 1, the TFFs F 1 and F 2 are non-regular and regular on [0, 2], respectively. Without loss of generality, no matter regular TFFs or non-regular TFFs, main results and the algorithm in this paper are all valid. In the sequel, we assume that TTFs are regular.

Main results and an algorithm
Now, we will consider the FFIDDC (1.1) in this section. With respect to the FFIDDC (1.1), we can solve it by using the method of steps introduced by Fatullayev et al. [14].
First of all, we deal with the FFIDDC (1. , the FFIDDC (1.1) is equivalent to FCFDIVP as follows: Thus, we transform the FCFDIVP (3.1) in matrix form: Here, is called to be a solution of the problem (3.2), which is also a solution to the FCFDIVP (3.1), where φ(0) = Z(0).
According to Definition 3.1, the solution Z is a fuzzy bunch of real functional matrices, which consists of functional matrices such as Z(t). If a functional matrix Z(t) satisfies for some functional matrices g ∈ supp( G) and φ(0) ∈ supp( Φ(0)), then it has a positive MD.
Since the FCFDIVP (3.1) is linear, we can solve it using the superposition principle. This follows that in order to find the solution of the FCFDIVP (3.1), we can separately consider the following subproblems: (1) The associated crisp problem: (2) The problem with initial TFFs: 3) The problem with fuzzy source functions and zero initial functions: x(0) = 0, y(0) = 0.
(3.6) By solving the above three subproblems (3.4)-(3.6), the solution of the FFIDDC (1.1) on interval [0, τ] can be obtained. Similarly, at the next step, we can find the solution on interval [τ, 2τ]. Hence, we can conclude that the solution of the FFIDDC (1.1) exists for t ∈ [0, ∞) using the method of steps [14].
In the sequel, we make clear how to solve each of these three subproblems (3.4)-(3.6).

The associated crisp problem
To solve the problem (3.4), we provide the following theorem.
Theorem 3.1. For the problem (3.4), we can get a unique solution to it, which is denoted as Z cr (t) = x cr (t) y cr (t) , when the following conditions fulfill: (H 2 ) There exist constants K i , L i ∈ (0, 1) for i = 1, 2 such that and for each x 1 , x 2 , y 1 , y 2 ∈ R, Proof. The proof is similar to Theorem 3.1 of Dong et al. [13] and it is omitted.

The problem with initial TFFs
In this subsection, we will give solvability of the problem (3.5) as follows.
, respectively, then the problem (3.5) has a unique solution Z φ , which is a TFF matrix given by Proof. On the one hand, according to Definition 3.1, it is easy for us to know that each Z(t) with non-zero MD from the bunch Z φ is a solution of the problem (3.7) for some On the other hand, in terms of Definition 2.1, the bunch φ 1 = φ 1a , 0, φ 1b consists of functions kφ 1a and kφ 1b ([0, 1] k = 1 − α) owing to φ 1c = 0. And the bunch φ 2 = φ 2a , 0, φ 2b is similar to the former.
From the above reasoning we can get the conclusion that the bunch Z φ consists of kZ a and kZ b . Therefore, the bunch Z φ is a TFF matrix determined to be Z φ = Z a , 0, Z b . Actually, we can get the solutions Z a (t) = φ 1a (0) φ 2a (0) and Z b (t) = φ 1b (0) φ 2b (0) .
We can express the value of the TFF matrix Z φ (3.8) at a time t by the formula Z φ (t) = (min{Z a (t), 0, Z b (t)}, 0, max{Z a (t), 0, Z b (t)}). Note that this value is a matrix of TFNs.

The problem with fuzzy source functions and zero initial functions
As for the solution of the problem (3.6), we give the theorem as follows.
are solutions of the following problem: x(0) = 0, in several, then the problem (3.6) has a unique solution Z g , which is a TFF matrix given by Proof. Since g 1a and g 2a satisfy (H 1 ) and (H 2 ), it follows from Theorem 3.1 that there will be a unique solution when g 1 g 2 = g 1a g 2a in (3.9). By the same token, from the assumption for g 1 g 2 = g 1b g 2b in the problem (3.9), only one solution will exist in this case. Thus, the proof can be done by the same way as in Theorem 3.2.
Thus, we can also express the value of the TFF matrix Z g (3.10) at a time t by Z g (t) = (min{Z u (t), 0, Z v (t)}, 0, max{Z u (t), 0, Z v (t)}), which is a matrix of TFNs.

A solution algorithm
By integrating Theorems 3.1-3.3 above, we can obtain the following solution algorithm for solving the FFIDDC (1.1). Step 2. Represent the initial values and source functions as Step 3. Find the solution Z cr (t) of the problem (3.4).
Step 4. Seek the solutions Z a (t) and Z b (t) of the problem (3.7) in regard to , respectively, and define Step 5. Solve the problem (3.9) and denote the solutions by Z u (t) and Z v (t), corresponding to g 1 g 2 = g 1a g 2a and , respectively, and let Step 6. Construct the unique solution of the FFIDDC (1.1) on interval [0, τ] as follows: Similarly, we can find the unique solution on interval [τ, 2τ], · · · . In consequence, we obtain a unique solution of the FFIDDC (1.1) which exists at t ≥ 0.

An example
In this section, referring to Example 1 in Fatullayev et al. [14], we take an example to clarify that the hypotheses in Theorems 3.1-3.3 can be satisfied, which can be also employed to verify the main results presented in Dong et al. [13].
where In the sequel, we only deal with the problem (4.1) on interval 0, π 2 . The unique solution on interval [τ, 2τ], · · · can be likewise found.
Next, we respectively solve the three subproblems of the problem (4.2) on 0, π 2 one by one. (i) With respect to g 1cr and g 2cr , we can rewrite them in another forms:  We can easily see that g 1cr , g 2cr : [0, π 2 ] × R × R → R are continuous and satisfy (H 1 ).
To sum up, all the hypotheses (H 1 ) and (H 2 ) in Theorem 3.1 are satisfied and we can obtain a unique solution Z cr (t), which is illustrated in Figure 2 to the following associated crisp subproblem via using the "L1 method" due to Li and Zeng [23]: which is the first subproblem of the problem (4.2). From Figure 2, one can easily see that the order of magnitude of the ordinate value in Figure 2 is 10 8 , and x cr and y cr established by relying on the "L1 method" discretization [23] of the Caputo fractional derivatives in (4.3), are growing very fast, and after iterating 43 times (t ≈ 0.66), the values of x cr and y cr approach 1 and 27.88, respectively. Since the ordinate value of Figure 2 is too large to observe between t = 0.6 and t = 0.7, we give the local graph Figure 3  (ii) Solve the following crisp problem: is the TFF matrix Z φ = Z a , 0, Z b , which is graphed in Figure 4, and one knows that Hence, Theorem 3.2 holds and there is one and only solution to the problem (4.4). Indeed, since the solutions Z a (t) and Z b (t) are all constant matrices in Figure 4, it is easy to see that the values of x φ 1 and y φ 2 are all constants when MD α ∈ [0, 1] is determined. That is, the values of x φ 1 and y φ 2 remain constant as t increases. For example, x φ 1a = 0.25, x φ 1b = −0.25, y φ 2a = 0.15 and y φ 2b = −0.15 when α = 0. (iii) Consider the following crisp problem: x(0) = 0, , and find the solutions respectively. That is, one solves the following crisp problems: and  As for the problem (4.5), g 1a and g 2a can be written in following forms:  It is clear to find that g 1a , g 2a : [0, π 2 ] × R × R → R are continuous. On the other hand, let K 1 = 0.4, K 2 = 0.1, L 1 = 0.6 and L 2 = 0.2. Then and for each u 1 , u 2 , v 1 , v 2 , m 1 , m 2 , n 1 , n 2 ∈ R, According to the above, all the hypotheses (H 1 ) and (H 2 ) hold. And that on the basis of Theorem 3.1, we can obtain the unique solution Z u (t) to the last subproblem of the problem (4.2) involving fuzzy source functions and zero initial functions as follows: Similarly, we can conclude that the problem (4.6) also has a unique solution Z v (t). It follows that the solution to the problem (4.7) is the TFF matrix Z g = Z u , 0, Z v , which is depicted in Figure 5, and we have Z g (t) = (min{Z u (t), 0, Z v (t)}, 0, max{Z u (t), 0, Z v (t)}) .
As a consequence, Theorem 3.3 holds and there is a unique solution to the problem (4.7). Moreover, it can be seen that, in Figure 5, the images of x g 1 and y g 2 are respectively symmetric when the corresponding MDs are chosen. For MD α ∈ [0, 1), as t increases, the upper branch of x g 1 firstly increases and then decreases (t ≈ 0.69 is the turning point), and the upper branch of y g 2 gets bigger and bigger. Especially, the values of x g 1 and y g 2 are constant 0 when the MD α = 1. However, the ordinate value of Figure 5 is of order of magnitude 10 −1 and very small compared with that in Figure 2, so it has little effect on the final solution Z(t).
In conclusion, by solving the problems (4.3), (4.4) and (4.7), it follows from Algorithm 3.1 that we can find the unique solution of the problem (4.1) on interval 0, π 2 as Z(t) = Z cr (t) + Z φ (t) + Z g (t), which can be found from Figure 6. It is easy to note that the image of fuzzy solution Z(t) in Figure 6 has a similar trend to that of the crisp solution Z cr (t) in Figure 2. That is because when MD is fixed, the unique solution Z(t) of the problem (4.1) is composed of three parts and the value of the crisp solution Z cr (t) for the problem (4.3) is too large (order of magnitude 10 8 of the ordinate value for Figure 2). Further, since the ordinate value of Figure 6 is very big (order of magnitude 10 8 ) and the curves of x, y can not be distinguished obviously, a local graph Figure 7 of Figure 6 is given for the fuzzy solutions of the problem (4.1) when the ordinate value interval is restricted to [0,30]. Choosing specially MD α = 0, we have x φ 1a = 0.93, x φ 1b = 1.07, y φ 2a = 27.68 and y φ 2b = 28.08. Before the end of this section, we remark that the unique solution of the problem (4.1) on interval π 2 , π , · · · can be found in the same way. Therefore, one can find the unique solution of the problem (4.1) for t > 0 via Algorithm (3.1).

Concluding remarks
In this paper, inspired by Fatullayev et al. [14] and Dong et al. [13], by using the methods of steps and separating fuzzy items, we solved the implicit coupled systems of fuzzy fractional delay differential equations with fuzzy initial values and source functions (FFIDDC) (1.1). Under the conditions (H 1 ) and (H 2 ), and based on the concept of TFFs and examples in [16,17], we obtained existence and uniqueness of solution for the FFIDDC (1.1). We note that the solution is a triangular fuzzy function matrix which is composed of real functional matrices. Furthermore, an algorithm for solution procedure and an example for graphical confirmation of explicit solution were given.
As Remark 1.1 mentioned, except for considering the problem (1.1) with Caputo fractional derivative, we will explore the problem of form (1.1) with fractional derivatives of Riemann-Liouville or Hadamard in the future. Recently, there are more and more researchers to explore the relevant (fractional) differentiable set-value problems and impulsive problems, such as Amrahov et al. [4] and Sun et al. [32]. Therefore, for further study, we can introduce the set-valued delay and the impulse into the FFIDDC (1.1), and extend the methods proposed in this paper to solve the new problems associated with set-values or impulsive elements.