A Fuzzy Method for Solving Fuzzy Fractional Differential Equations Based on the Generalized Fuzzy Taylor Expansion

mathematics Article A Fuzzy Method for Solving Fuzzy Fractional Differential Equations Based on the Generalized Fuzzy Taylor Expansion Tofigh Allahviranloo 1,*, Zahra Noeiaghdam 2, Samad Noeiaghdam 3,4 and Juan J. Nieto 5 1Faculty of Engineering and Natural Sciences, Bahcesehir University, 34353 Istanbul, Turkey 2Department of Mathematics, Shahed University, Tehran 3319118651, Iran; [email protected] 3Department of Applied Mathematics and Programming, South Ural State University, Lenin Prospect 76, 454080 Chelyabinsk, Russia; [email protected] 4Baikal School of BRICS, Irkutsk National Research Technical University, 664074 Irkutsk, Russia 5Instituto de Matemáticas, Departamento de Estatística, Análise Matemática e Optimización, Universidade de Santiago de Compostela, 15782 Santiago de Compostela, Spain; [email protected] *Correspondence: [email protected] Received: 31 October 2020; Accepted: 30 November 2020; Published: 4 December 2020   Abstract: In this field of research, in order to solve fuzzy fractional differential equations, they are normally transformed to their corresponding crisp problems. This transformation is called the embedding method. The aim of this paper is to present a new direct method to solve the fuzzy fractional differential equations using fuzzy calculations and without this transformation. In this work, the fuzzy generalized Taylor expansion by using the sense of fuzzy Caputo fractional derivative for fuzzy-valued functions is presented. For solving fuzzy fractional differential equations, the fuzzy generalized Euler’s method is introduced and applied. In order to show the accuracy and efficiency of the presented method, the local and global truncation errors are determined. Moreover, the consistency , convergence, and stability of the generalized Euler’s method are proved in detail. Eventually, the numerical examples, especially in the switching point case, show the flexibility and the capability of the presented method. Keywords: fuzzy fractional differential equations; generalized fuzzy Taylor expansion; generalized fuzzy Euler’s method; global truncation error; local truncation error; convergence; stability 1. Introduction Fuzzy set theory is a powerful tool for modeling uncertain problems. Therefore, large varieties of natural phenomena have been modeled using fuzzy concepts. Particularly, the fuzzy fractional differential equation is a common model in different scientific fields, such as population models, evaluating weapon systems, civil engineering, and modeling electro-hydraulics. Therefore, the concept of the fractional derivative is a very important topic in fuzzy calculus. Therefore, fuzzy fractional differential equations have attracted much attention in mathematics and engineering fields. The first work devoted to the subject of fuzzy fractional differential equations is the paper by Agarwal et al. [1]. They have defined the Riemann–Liouville differentiability concept under the Hukuhara differentiability to solve fuzzy fractional differential equations. In recent years, fractional calculus has been introduced as an applicable topic to produce the accurate results of mathematical and engineering problems such as aerodynamics and control systems, signal processing, bio-mathematical problems, and others [1–5]. Furthermore, fractional differential equations in the fuzzy case [ 1 ] have been studied by many authors and they have been solved by various methods [ 6 – 9 ]. In [ 10 ], Hoa studied the fuzzy Mathematics 2020,8, 2166; doi:10.3390/math8122166 www.mdpi.com/journal/mathematics Mathematics 2020,8, 2166 2 of 24 fractional differential equations under Caputo gH-differentiability, and in [ 11 ] Agarwal et al. had a survey on mentioned problem to show the its relation with optimal control problems. Furthermore, Long et al. [12] illustrated the solvability of fuzzy fractional differential equations, and Salahshour et al. [13] applied the fuzzy Laplace transforms to solve this problem. There are many numerical methods to solve the fuzzy fractional differential equations by transforming to crisp problems [ 14 – 16 ]. In this paper, a new direct method is introduced to solve the mentioned problem without changing to the crisp form. The Taylor expansion method is one of the famous and applicable methods to solve the linear and nonlinear problems [ 17 – 19 ]. In this paper, the fuzzy generalized Taylor expansion based on the fuzzy Caputo fractional derivative is expanded. Then, Euler’s method is applied to solve the fuzzy fractional differential equations. Furthermore, the local and global truncation errors are considered, and finally the consistency, the convergence and the stability of the generalized fuzzy Euler’s method are demonstrated. Furthermore, some examples with the switching point are solved by using the presented method. The numerical results show the precision of the generalized Euler’s method to solve the fuzzy fractional differential equations. 2. Basic Concepts At first, a brief summary of the fuzzy details and some preliminaries is revisited [20–25]. Definition 1. Set RF={u:Rn→[0, 1]such that u satisfies in the conditions I to IV. } I. u is normal: there exists an x0∈Rnsuch that u(x0) = 1. II. u is fuzzy convex: for 0≤λ≤1, u(λx1+ (1−λx2)) ≥min{u(x1),u(x2)}. III. u is upper semi-continuous: for any x0∈Rn, it holds that u(x0)≥limx→x± 0u(x). IV. [u]0=supp(u) = cl{x∈Rn|u(x)>0}is a compact subset. is called the space of fuzzy numbers or the fuzzy numbers set. The r -level set is [u]r={x∈Rn|u(x)≥ r , 0 <r≤ 1 } . Then, from I to IV, it follows that, the r -level sets of u∈RF are nonempty, closed, and bounded intervals. Definition 2. A triangular fuzzy number is defined as a fuzzy set in RF , that is specified by an ordered triple u= (a , b , c)∈R3 with a≤b≤c such that u−(r) = a+ (b−a)r (or lower bound of u ) and u+(r) = c−(c−b)r (or upper bound of u) are the endpoints of r-level sets for all r ∈[0, 1]. A crisp number k is simply represented by u−(r) = u+(r) = k , 0 ≤r≤ 1 and called singleton. For arbitrary u , v∈RF and scalar k , we might summarize the addition and the scalar multiplication of two fuzzy numbers by /addition :[u⊕v]r= [u−(r) + v−(r),u+(r) + v+(r)], /scalar multiplication :     [ku]r= [ku−(r),ku+(r)],k≥0, [ku]r= [ku+(r),ku−(r)],k<0. The Hausdorff distance between fuzzy numbers is given by H:RF×RF→R+∪{0}as H(u,v) = sup 0≤r≤1 max{|u−(r)−v−(r)|,|u+(r)−v+(r)|}, where [u]r= [u−(r) , u+(r)] , [v]r= [v−(r) , v+(r)] . The metric space (RF , H) is complete, separable, and locally compact where the following conditions are valid for metric H: I. H(u⊕w,v⊕w) = H(u,v),∀u,v,w∈RF. II. H(λu,λv) =|λ| H(u,v),∀λ∈R,∀u,v∈RF. Mathematics 2020,8, 2166 3 of 24 III. H(u⊕v,w⊕z)≤ H(u,w) + H(v,z),∀u,v,w,z∈RF. Definition 3. Let u , v∈RF , if there exists w∈RF ,such that u=v+w ,then w is called the Hukuhara difference (H-difference) of u and v ,and it is denoted by uv . Furthermore, the generalized Hukuhara difference (gH-difference) of two fuzzy numbers u,v∈RFis defined as follows, ugH v=w⇔     (i)u=v⊕w, or (ii)v=u⊕(−1)w, it is easy to show that conditions (i) and (ii) are valid if and only if w is a crisp number. The conditions of the existence of ugH v∈RF are given in [ 22 ]. Through the whole of the paper, we suppose that the gH -difference exists. In this paper, the meaning of fuzzy-valued function is a function f:A→RF , A∈R where R is the set of all real numbers and [f(t)]r= [f−(t ; r) , f+(t ; r)] is the so-called r -cut or parametric form of the fuzzy-valued function f. Definition 4. A fuzzy-valued function f:[a , b]→RF is said to be continuous at t0∈[a , b] if for each ε> 0 there is δ> 0such that H(f(t) , f(t0)) <ε , whenever t∈[a , b] and |t−t0|<δ . We say that f is fuzzy continuous on [a,b]if f is continuous at each t0∈[a,b]. Throughout the rest of this paper, the notation Cf([a , b] , RF) is called the set of fuzzy-valued continuous functions which are defined on [a,b]. If f:[a , b]⊆R→RF is continuous by the metric H , then Rt af(s)ds is a continuous function in t∈[a,b]and the function fis integrable on [a,b]. Furthermore, it holds Zb af(s)dsr =Zb af−(s;r)ds,Zb af+(s;r)ds. Definition 5. Let f:[a , b]→RF , t0∈(a , b) with f−(t ; r) and f+(t ; r) both differentiable at t0 for all r∈[0, 1]and D fgH (gH-derivative) exists: I. The function f is F[(i)−gH]-differentiable at t0if [D fi.gH(t0)]r= [D f −(t0;r),D f +(t0;r)]. II. The function f is F[(ii)−gH]-differentiable at t0if [D fii.gH(t0)]r= [D f +(t0;r),D f −(t0;r)]. 3. Definitions and Properties of Fractional gH-Differentiability In this section, let us focus on some definitions and properties related to the fuzzy fractional generalized Hukuhara derivative which are useful in the sequel of this paper. Definition 6 ([ 26 ]) . Let f(t) be a fuzzy Lebesque integrable function. The fuzzy Riemann–Liouville fractional (for short (F.RL)-fractional) integral of order α>0is defined as follows, F.RL Iα [a,t]f(t) = 1 Γ(α)Zt a(t−s)α−1f(s)ds. Definition 7 ([ 26 ]) . Let f:[a , b]→RF . The fuzzy fractional derivative of f(t) in the Caputo sense is in the following form, FCDα ∗f(t) = F.RL Im−α [a,t](DmfgH)(t) = 1 Γ(m−α)Zt a(t−s)m−α−1DmfgH(s)ds, m−1<α<m,m∈N,t>a, Mathematics 2020,8, 2166 4 of 24 where ∀m∈N , DmfgH(s) ( gH -derivatives of f ) are integrable. In this paper, we consider fuzzy Caputo generalized Hukuhara derivative (for short FC[gH] -derivative) of order 0 <α≤ 1, for fuzzy-valued function f , so the FC[gH]-derivative will be expressed by FCDα ∗f(t) = F.RL I1−α [a,t](D fgH)(t) = 1 Γ(1−α)Zt a(t−s)−αD fgH(s)ds,t>a. (1) Lemma 1. Let f:[a , b]⊆R→RF be continuous. Then, F.RL Iα [a,t]f(t) for 0 <α≤ 1and t∈[a , b] is a continuous function. Proof. Under assumptions of the continuous functions, f(s) is a fuzzy Lebesque integrable function. On the other hand, since ∀ 0 <α≤ 1, (t−s)α−1≥ 0 is continuous, so Rt a(t−s)α−1f(s)ds is a continuous function and as a result F.RL Iα [a,t]f(t)is a continuous function in t∈[a,b]. Lemma 2. Let f∈ Cf(R , RF) , m∈N . Then, the fuzzy Riemann-Liouvil fractional integrals F.RL Iα [a,tm−1]f(tm−1),F.RL Iα [a,tm−2](F.RL Iα [a,tm−1]f)(tm−1)), ...,F.RL Iα [a,t](F.RL Iα [a,t1]...(F.RL Iα [a,tm−2] (F.RL Iα [a,tm−1]f)(tm−1)) ... ) for 0 <α≤ 1, are continuous functions in tm−1 , tm−2 , ..., t , respectively. Here, tm−1,tm−2, ..., t≥a, and they are real numbers. Proof. This lemma is a fairly straightforward generalization of Lemma 1. The proof will be done by introducing on m∈N . Assume that the lemma holds for (m) -times applying operator (F . RL) -fractional integrating for function f , we will prove it will correct for (m+ 1 ) -times applying operator (F . RL) -fractional integrating for function f . By Lemma 1, as f∈ Cf(R , RF) , thus F.RL Iα [a,tm−1]f(tm−1) is a continuous function in tm−1. Furthermore, under the hypothesis of induction, F.RL Iα [a,tm−2](F.RL Iα [a,tm−1]f)(tm−1), ..., (m)−times z }| { F.RL Iα [a,t](F.RL Iα [a,t1]...(F.RL Iα [a,tm−2](F.RL Iα [a,tm−1]f)(tm−1))...), are continuous functions in tm−1,tm−2,tm−3, ..., t, respectively. It follows easily that (m+1)−times z }| { F.RL Iα [a,tm+1](F.RL Iα [a,t](F.RL Iα [a,t1]...(F.RL Iα [a,tm−2](F.RL Iα [a,tm−1]f)(tm−1))...)), is a continuous function in tm+1, which is our claim. Definition 8 ([ 26 ]) . Let f:[a , b]→RF be the fuzzy Caputo generalized Hukuhara differentiable (for short FC[gH] -differentiable) at t0∈[a , b] .Thus, f is FC[(i)−gH] -differentiable at t0∈[a , b] if for 0≤r≤1 [FCDα ∗fi.gH(t0)]r= [CDα ∗f−(t0;r),CDα ∗f+(t0;r)], and that f is FC[(ii)−gH]-differentiable at t0if [FCDα ∗fii.gH(t0)]r= [CDα ∗f+(t0;r),CDα ∗f−(t0;r)], where CDα ∗f−(t0;r) = 1 Γ(1−α)Zt0 a(t0−s)−αD f −(s;r)ds, CDα ∗f+(t0;r) = 1 Γ(1−α)Zt0 a(t0−s)−αD f +(s;r)ds. Mathematics 2020,8, 2166 5 of 24 Definition 9 ([ 26 ]) . Let f:[a , b]→RF be a fuzzy-valued function on [a , b] . A point t0∈[a , b] is said to be a switching point for the FC[gH] -differentiability of f ,if in any neighborhood V of t0 there exist points t1<t0<t2such that (type I) f is FC[(i)−gH]-differentiable at t1while f is not FC[(ii)−gH]-differentiable at t1, and f is FC[(ii)−gH]-differentiable at t2while f is not FC[(i)−gH]-differentiable at t2, or (type II) f is FC[(ii)−gH]-differentiable at t1while f is not FC[(i)−gH]-differentiable at t1, and f is FC[(i)−gH]-differentiable at t2while f is not FC[(ii)−gH]-differentiable at t2. Theorem 1 ([ 1 ]) . If f:[a , b]→RF , [f(t)]r= [f−(t ; r) , f+(t ; r)] and f is integrable for 0 ≤r≤ 1, t∈[a , b] and α,β>0, then we have F.RL Iα [a,t](F.RL Iβ [a,t]f)(t) = F.RL Iα+β [a,t]f(t). Lemma 3 ([ 1 ]) . Suppose that f:[a , b]→RF be a fuzzy-valued function and D fgH is exist, then for 0<α≤1, F.RL Iα [a,t](FCDα ∗f)(t) = f(t)gH f(a), 0 ≤r≤1. The principal significance of this lemma is in the following theorem. Theorem 2 ([ 1 ]) . Let f:[a , b]→RF be the fractional gH -differentiable such that type of Caputo differentiability f in [a,b]does not change. Then, for a ≤t≤b and 0<α≤1, I. If f (s)is FC[(i)−gH]-differentiable then FCDα ∗fi.gH(t)is (F.RL)-integrable over [a,b]and f(t) = f(a)⊕F.RL Iα [a,t](FCDα ∗fi.gH)(t), II. If f (s)is FC[(ii)−gH]-differentiable then FCDα ∗fii.gH(t)is (F.RL)-integrable over [a,b]and f(t) = f(a)(−1)F.RL Iα [a,t](FCDα ∗fii.gH)(t). Lemma 4. Suppose that f:[a , b]→RF is the fractional gH -differentiable and FCDα ∗fgH(t)∈ Cf([a , b] , RF) then for 0<α≤1, F.RL Iα [t,a](FCDα ∗fi.gH)(t) = (−1)F.RL Iα [a,t](FCDα ∗fii.gH)(t), Proof. As FCDα ∗fi.gH(t) is continuous, it follows that FCDα ∗fi.gH(t) is the Riemann–Liouville integrable, and by using Lemma 3for 0 ≤r≤1 [F.RL Iα [t,a](FCDα ∗fi.gH)(t)]r= [F.RL Iα [t,a](FCDα ∗f−)(t;r),F.RL Iα [t,a](FCDα ∗f+)(t;r)] = [f−(a;r)−f−(t;r),f+(a;r)−f+(t;r)] = [f(a)f(t)]r, (2) moreover, [F.RL Iα [a,t](FCDα ∗fii.gH)(t)]r= [F.RL Iα [a,t](FCDα ∗f+)(t;r),F.RL Iα [a,t](FCDα ∗f−)(t;r)] = [f+(t;r)−f+(a;r),f−(t;r)−f−(a;r)] = [(−1)(f(a)f(t))]r, (3) by combining Equations (2) with (3) the lemma is proved. Mathematics 2020,8, 2166 6 of 24 Theorem 3. Let FCDα ∗f:[a,b]→RFand FCDnα ∗f∈ Cf([a,t],RF).For all t ∈[a,b]and 0<α≤1, I. Let FCDjα ∗f , j= 1, ..., n be the FC[(i)−gH] -differentiable, and they do not change in the type of differentiability on [a,b],then FCD(j−1)α ∗fi.gH(t) = FCD(j−1)α ∗fi.gH(a)⊕F.RL Iα [a,t](FCDjα ∗fi.gH)(t). II. If FCDjα ∗f , j= 1, ..., n are FC[(ii)−gH] -differentiable and the type of their differentiability does not change in the interval [a,b],then FCD(j−1)α ∗fii.gH(t) = FCD(j−1)α ∗fii.gH(a)⊕F.RL Iα [a,t](FCDjα ∗fii.gH)(t). III. Assume that FCDjα ∗f , j= 2 k− 1, k∈N are the FC[(i)−gH] -differentiable and they are FC[(ii)− gH]-differentiable, for j =2k,k∈Nthen FCD(j−1)α ∗fi.gH(t) = FCD(j−1)α ∗fi.gH(a)(−1)F.RL Iα [a,t](FCDjα ∗fii.gH)(t). IV. Suppose that FCDjα ∗f , j= 2 k− 1, k∈N are FC[(ii)−gH] -differentiable and they are FC[(i)− gH]-differentiable for j =2k,k∈N, so FCD(j−1)α ∗fii.gH(t) = FCD(j−1)α ∗fii.gH(a)(−1)F.RL Iα [a,t](FCDjα ∗fi.gH)(t). Proof. By assuming FCDjα ∗f∈ Cf([a , b] , RF) , j= 0, ..., n we give the proof only for parts II and III . Proving the other parts are similar. II. Our proof starts with the observation that FCDjα ∗f , j= 1, ..., n are FC[(ii)−gH] -differentiable. Therefore, using the properties of fuzzy Caputo derivative and Theorem 2, we have [F.RL Iα [a,t](FCDjα ∗fii.gH)(t)]r = [RL Iα [a,t](CDjα ∗f+)(t;r),RL Iα [a,t](CDjα ∗f−)(t;r)] = [RL Iα [a,t].CDα ∗(CD(j−1)α ∗f+)(t;r),RL Iα [a,t].CDα [a,t](CD(j−1)α ∗f−)(t;r)] = [CD(i−1)α ∗f+(t;r)−CD(j−1)α ∗f+(a;r),CD(j−1)α ∗f−(t;r)−CD(j−1)α ∗f−(a;r)] = [CD(j−1)α ∗f+(t;r),CD(j−1)α ∗f−(t;r)] −[CD(j−1)α ∗f+(a;r),CD(j−1)α ∗f−(a;r)] = [FCD(j−1)α ∗fii.gH(t)FCD(j−1)α ∗fii.gH(a)]r, thus, we obtain FCD(j−1)α ∗fii.gH(t) = FCD(j−1)α ∗fii.gH(a)⊕F.RL Iα [a,t](FCDjα ∗fii.gH)(t). Mathematics 2020,8, 2166 7 of 24 III. Under the conditions stated in the part III, FCDjα ∗f is FC[(i)−gH] -differentiable for j= 2 k− 1, k∈N and it is FC[(ii)−gH] -differentiable for j= 2 k , k∈N . In the sense of Section 2and by Theorem 2, we get [FCD(j−1)α ∗fi.gH(t)⊕(−1)F.RL Iα [a,t](FCDjα ∗fii.gH)(t)]r = [CD(j−1)α ∗f−(t;r),CD(j−1)α ∗f+(t;r)] + [−RL Iα [a,t](CDjα ∗f−)(t;r),−RL Iα [a,t](CDjα ∗f+)(t;r)] = [CD(j−1)α ∗f−(t;r),CD(j−1)α ∗f+(t;r)] + [CD(j−1)α ∗f−(a;r)−CD(j−1)α ∗f−(t;r),CD(j−1)α ∗f+(a;r)−CD(j−1)α ∗f+(t;r)] = [CD(j−1)α ∗f−(a;r),CD(j−1)α ∗f+(a;r)] = [FCD(j−1)α ∗fi.gH(a)]r, which completes the proof. 4. Fuzzy Generalized Taylor Theorem Theorem 4. Let T= [a , a+β]⊂R , with β> 0and FCDjα ∗f∈ Cf([a , b] , Rf) , j= 1, ..., n . For t ∈T, 0 <α≤1 I. If FCDjα ∗f , j= 0, 1, ..., n− 1are FC[(i)−gH] -differentiable, provided that type of fuzzy Caputo differentiability has no change. Then, f(t) = f(a)⊕FCDα ∗fi.gH(a)(t−a)α Γ(α+1)⊕FCD2α ∗fi.gH(a)(t−a)2α Γ(2α+1) ⊕... ⊕FCD(n−1)α ∗fi.gH(a)(t−a)(n−1)α Γ((n−1)α+1)⊕Rn(a,t), where Rn(a,t):=F.RL Iα [a,t](F.RL Iα [a,t1]...(F.RL Iα [a,tn−1](FCDnα ∗fi.gH)(tn))...). II. If FCDjα ∗f , j= 0, 1, ..., n− 1are FC[(ii)−gH] -differentiable, provided that type of fuzzy Caputo differentiability has no change, then f(t) = f(a)(−1)FCDα ∗fii.gH(a)(t−a)α Γ(α+1)(−1)FCD2α ∗fii.gH(a)(t−a)2α Γ(2α+1) (−1)... (−1)FCD(n−1)α ∗fii.gH(a)(t−a)(n−1)α Γ((n−1)α+1)(−1)Rn(a,t), where Rn(a,t):=F.RL Iα [a,t](F.RL Iα [a,t1]...(F.RL Iα [a,tn−1](FCDnα ∗fii.gH)(tn))...). III. If FCDjα ∗f,j=1, ..., n,exist and in each order the type of FC[gH]-differentiability changes on T f(t) = f(a)(−1)FCDα ∗fii.gH(a)(t−a)α Γ(α+1)⊕FCD2α ∗fi.gH(a)(t−a)2α Γ(2α+1) (−1)... (−1)FCD(j 2−1)α ∗fii.gH(a)(t−a)(j 2−1)α Γ(j 2α) ⊕FCD(j 2)α ∗fi.gH(a)(t−a)(j 2)α Γ(j 2α+1)(−1)... (−1)Rn(a,t), Mathematics 2020,8, 2166 8 of 24 where Rn(a,t):=F.RL Iα [a,t](F.RL Iα [a,t1]...(F.RL Iα [a,tn−1](FCDnα ∗fi.gH)(tn))...). IV. For FCDkα ∗f∈ Cf([a , b] , Rf) , k≥ 3, suppose that f on [a , ξ] is FC[(ii)−[gH] -differentiable and on [ξ , b] is FC[(i)−gH] -differentiable, in fact ξ is switching point (type II) for α -order derivative of f . Moreover, for t0∈[a , ξ] , let 2 α -order derivative of f in ξ1 of [t0 , ξ] have switching point (type I). On the other hand, the type of differentiability for FCDjα ∗f,j≤k on [ξ,b]does not change. Therefore, f(t) = f(t0)(−1)FCDα ∗fii.gH(t0)(ξ−t0)α Γ(α+1)⊕(−1)FCD2α ∗fi.gH(t0)(t0−ξ1)α Γ(α+1) (ξ−t0)α Γ(α+1)(−1)FCD2α ∗fii.gH(ξ1)(ξ−ξ1)2α Γ(2α+1)−(t0−ξ1)2α Γ(2α+1) ⊕FCDα ∗fi.gH(ξ)(t−ξ)α Γ(α+1)⊕FCD2α ∗fi.gH(ξ)(t−ξ)2α Γ(2α+1) ⊕F.RL Iα [t0,ξ].F.RL I[t0,ξ1].F.RL Iα [t0,t2](FCD3α ∗fi.gH)(t4) (−1)F.RL Iα [t0,ξ].F.RL I[ξ1,t1].F.RL Iα [ξ1,t3](FCD3α ∗fii.gH)(t5). ⊕F.RL I[ξ,t].F.RL Iα [ξ,s1].F.RL Iα [a,s2](FCD3α ∗fi.gH)(s3). Proof. Here, we prove parts II , III , and IV because the proof of part I is similar to the part II . Under the assumptions that FCDjα ∗f∈ Cf([a , b] , Rf) , j= 1, ..., n , we conclude that FCDjα ∗f are (F.RL)-fractional integrable on T, II. As fis a continuous function and FC[(ii)−gH]-differentiable, by Theorem 2, we get f(t) = f(a)(−1)F.RL Iα [a,t](FCDα ∗fii.gH)(t1). Under the hypotheses of Theorem, type of differentiability does not change, so by Theorem 3and by attention to (F.RL)-integrability of FCDα ∗fii.gH on T, we obtain FCDα ∗fii.gH(t1) = FCDα ∗fii.gH(a)⊕F.RL Iα [a,t1](FCD2α ∗fii.gH)(t2). Applying operator F.RL Iα [a,t]to FCDα ∗fii.gH(t1), gives F.RL Iα [a,t](FCDα ∗fii.gH)(t1) = FCDα ∗fii.gH(a)(t−a)α Γ(α+1)⊕F.RL Iα [a,t].F.RL Iα [a,t1](FCD2α ∗fii.gH)(t2). Lemma 2implies that the last double (F.RL)-fractional integral belongs to Rf. Therefore, f(t) = f(a)(−1)FCDα ∗fii.gH(a)(t−a)α Γ(α+1)(−1)F.RL Iα [a,t].F.RL Iα [a,t1](FCD2α ∗fii.gH)(t2), (4) by repeating the above argument, we get FCD2α ∗fii.gH(t2) = FCD2α ∗fii.gH(a)⊕F.RL Iα [a,t2](FCD3α ∗fii.gH)(t3). Therefore, we find that F.RL Iα [a,t1](FCD2α ∗fii.gH)(t2) = FCD2α ∗fii.gH(a)(t1−a)α Γ(α+1)⊕F.RL Iα [a,t1].F.RL Iα [a,t2](FCD3α ∗fii.gH)(t3), Mathematics 2020,8, 2166 9 of 24 moreover F.RL Iα [a,t].F.RL Iα [a,t1](FCD2α ∗fii.gH)(t2) = FCD2α ∗fii.gH(a)(t−a)2α Γ(2α+1) ⊕F.RL Iα [a,t].F.RL Iα [a,t1].F.RL Iα [a,t2](FCD3α ∗fii.gH)(t3). By Lemma 2, the last triple (F . RL) -fractional integral belongs to Rf . Therefore, substituting above equation into Equation (4), we find that f(t) = f(a)(−1)FCDα ∗fii.gH(a)(t−a)α Γ(α+1)(−1)FCD2α ∗fii.gH(a) (t−a)2α Γ(2α+1)(−1)F.RL Iα [a,t].F.RL Iα [a,t1].F.RL Iα [a,t2](FCD3α ∗fii.gH)(t3), the high order of the last formula by Lemma 2is a continuous function in terms of t so it belongs to Rf. With the same manner, we can demonstrate that part II is satisfied. III. Suppose that fis FC[(ii)−gH]-differentiable. Using Theorem 2, we have f(t) = f(a)(−1)F.RL Iα [a,t](FCDα ∗fii.gH)(t1). Under the hypothesis of theorem, as f is FC[(ii)−gH] -differentiable, FCDα ∗f is FC[(i)− gH]-differentiable. Therefore, by Theorem 3we get FCDα ∗fii.gH(t1) = FCDα ∗fii.gH(a)(−1)F.RL Iα [a,t1](FCD2α ∗fi.gH)(t2). Now, applying operator (F.RL)-integral to (FCDα ∗fii.gH)(t1)gives F.RL Iα [a,t](FCDα ∗fii.gH)(t1) =F.RL Iα [a,t](FCDα ∗fii.gH)(a)(−1)F.RL Iα [a,t].F.RL Iα [a,t1](FCD2α ∗fi.gH)(t2) =FCDα ∗fii.gH(a)(t−a)α Γ(α+1)(−1)F.RL Iα [a,t].F.RL Iα [a,t1](FCD2α ∗fi.gH)(t2). Lemma 2now leads to the last double (F.RL)-fractional integral belongs to Rf. Therefore, f(t) = f(a)(−1)FCDα ∗fii.gH(a)(t−a)α Γ(α+1)⊕F.RL Iα [a,t].F.RL Iα [a,t1](FCD2α ∗fi.gH)(t2). Similarly, as FCDα ∗f is FC[(i)−gH] -differentiable, FCD2α ∗f is FC[(ii)−gH] -differentiable and we get FCD2α ∗fi.gH(t2) = FCD2α ∗fi.gH)(a)(−1)F.RL Iα [a,t2](FCD3α ∗fii.gH)(t3), thus F.RL Iα [a,t1](FCD2α ∗fi.gH)(t2) = FCD2α ∗fi.gH(a)(t1−a)α Γ(α+1) (−1)F.RL Iα [a,t1].F.RL Iα [a,t2](FCD3α ∗fii.gH)(t3). Now, applying operator F.RL Iα [a,t]gives F.RL Iα [a,t].F.RL Iα [a,t1](FCD2α ∗fi.gH)(t2) = FCD2α ∗fi.gH(a)(t−a)2α Γ(2α+1) (−1)F.RL Iα [a,t].F.RL Iα [a,t1].F.RL Iα [a,t2](FCD3α ∗fii.gH)(t3), Mathematics 2020,8, 2166 16 of 24 as the inequality holds for all k, we get H(y(tk+1),yk+1)≤(1−hα Γ(α+1).`)(1−hα Γ(α+1).`)H(y(tk−1),yk−1)+r+r = (1−hα Γ(α+1).`)2.H(y(tk−1),yk−1)+r1+ (1−hα Γ(α+1).`). Repeated application of the above inequality enables us to write H(y(tk+1),yk+1) ≤(1−hα Γ(α+1).`)k+1.H(y(t0),y0)+r1+ (1−hα Γ(α+1).`) + ... + (1−hα Γ(α+1).`)k, obviously, this sum is a geometric series, so we have k ∑ i=0 (1−hα Γ(α+1).`)i=1−(1−hα Γ(α+1).`)k+1 hα Γ(α+1).`, that resulted to H(y(tk+1),yk+1)≤(1−hα Γ(α+1).`)k+1H(y(t0),y0)+rΓ(α+1) hα`1−(1−hα Γ(α+1).`)k+1, (23) with z=−hα Γ(α+1).`in Lemma 5concludes that (1−hα Γ(α+1).`)k+1≤e−hα Γ(α+1).`(k+1)≤e−`T Γ(α+1), where 0 ≤(k+1)hα≤Tfor (k+1)≤(N−1). Thus in Equation (23), we obtain H(y(tk+1),yk+1)≤e−`T Γ(α+1)H(y(t0),y0)+rΓ(α+1) hα`[1−e−`T Γ(α+1)]. Moreover, r=max 0≤k≤N−1H(rk, 0) = −h2α Γ(2α+1)max 0≤t≤TH(FCD2α ∗yii.gH(t), 0), and the accuracy of the initial value, concludes that H(y(t0),y0)=0, so H(y(tk+1),yk+1)≤ −hαΓ(α+1) `Γ(2α+1)[1−e−`T Γ(α+1)]max 0≤t≤TH(FCD2α ∗yii.gH(t), 0), now, letting h→0 then H(y(tk+1),yk+1)→ 0, which is the desired conclusion, and we can say that the fuzzy generalized Euler’s method is convergent in this step. Step II. To estimate the step II, consider y(t) is FC[(i)−gH] -differentiable, by using Equation (19) and let rk=h2α Γ(2α+1)FCD2α ∗yi.gH(tk) the proof of this step is similar to step I, so the fuzzy generalized Euler’s method is convergent. 6.3. Stability Now, the stability of the presented method is illustrated. For this aim, the following definition is presented. Definition 13. Assume that yk+1 , k+ 1 ≥ 0is the solution of fuzzy generalized Euler’s method where y0∈RF and also zk+1 is the solution of the same numerical method where z0=y0⊕δ0∈RF shows its Mathematics 2020,8, 2166 17 of 24 perturbed fuzzy initial condition. The fuzzy generalized Euler’s method is stable if there exists positive constant bh and Ksuch that ∀(k+1)hα≤T,k+1<N−1, h∈(0, bh)⇒ H(zk+1,yk+1)≤ Kδ whenever H(δ0, 0)≤δ. /Investigating the stability of the fuzzy generalized Euler’s method: The proof falls naturally into two steps: Step I. If y(t) is FC[(ii)−gH] -differentiable, by using Equation ( 18 ) the perturbed problem is in the following form, zk+1=zk(−1)hα Γ(α+1)f(tk,zk),z0=y0⊕δ0. (24) According to the Equations (19)and (24), we have H(zk+1,yk+1)≤ H(zk,yk)−hα Γ(α+1)H(f(tk,zk),f(tk,yk)), which we have been working under the assumption that specifications of the Hausdorff metric are satisfied. Using the Lipschitz condition, it can be concluded that H(zk+1,yk+1)≤(1−hα Γ(α+1)`)H(zk,yk), repeating with the inequality and applying Lemma 5lead us to the following inequality H(zk+1,yk+1)≤(1−hα Γ(α+1)`)k+1H(z0,y0) ≤e−hα Γ(α+1).`(k+1) Hz0gH y0, 0 ≤e− `T Γ(α+1)H(δ0, 0)≤ Kδ, where K=e− `T Γ(α+1) and for k+ 1 <N− 1 ⇒hα(k+ 1 )≤T . In this case, it is obvious the stability of the fuzzy generalized Euler’s method. Step II. For FC[(i)−gH] -differentiability of y(t) the same process can be used. In general, the above-mentioned analysis, points out that the fuzzy generalized Euler method is a stable approach. 7. Numerical Simulations In this section, several examples of the fractional differential equations are solved by using the full fuzzy generalized Euler method. Moreover, the numerical results are demonstrated on some tables for different values of hand t. Example 1. Let us consider the following initial value problem, FCDα ∗y(t) = (0, 1, 1.5)Γ(α+1), 0 <t≤1, Mathematics 2020,8, 2166 18 of 24 where y( 0 ) = 0, and y(t)=( 0, 1, 1.5 )tα is the exact FC[i−gH] -differentiable solution of problem. In order to find the numerical results we should construct the following iterative formula as yk+1=yk⊕hα Γ(α+1)[(0, 1, 1.5)Γ(α+1)],k=0, 1, ··· ,N−1, in Table 1, the numerical results for different values of t , α and h are demonstrated. In Figure 1, the exact solution and the Caputo gH-derivative for α=0.6 are demonstrated. (a) (b) Figure 1. ( a ) The exact solution ( b ) the Caputo gH-derivative of solution defined in Example 1for α=0.6. Table 1. Numerical results of Example 1for various t,αand h. α=0.3 α=0.6 α=0.9 t h =0.2 h=0.02 h=0.2 h=0.02 h=0.2 h=0.02 0.1 (0, 0.617034, 0.925551) (0, 0.309249, 0.463874) (0, 0.380731, 0.571096) (0, 0.0956352, 0.143453) (0, 0.234924, 0.352386) (0, 0.0295752, 0.0443627) 0.2 (0, 1.23407, 1.8511) (0, 0.618499, 0.927748) (0, 0.761462, 1.14219) (0, 0.19127, 0.286906) (0, 0.469848, 0.704771) (0, 0.0591503, 0.0887255) 0.3 (0, 1.8511, 2.77665) (0, 0.927748, 1.39162) (0, 1.14219, 1.71329) (0, 0.286906, 0.430359) (0, 0.704771, 1.05716) (0, 0.0887255, 0.133088) 0.4 (0, 2.46814, 3.7022) (0, 1.237, 1.8555) (0, 1.52292, 2.28438) (0, 0.382541, 0.573811) (0, 0.939695, 1.40954) (0, 0.118301, 0.177451) 0.5 (0, 3.08517, 4.62775) (0, 1.54625, 2.31937) (0, 1.90365, 2.85548) (0, 0.478176, 0.717264) (0, 1.17462, 1.76193) (0, 0.147876, 0.221814) 0.6 (0, 3.7022, 5.5533) (0, 1.8555, 2.78325) (0, 2.28438, 3.42658) (0, 0.573811, 0.860717) (0, 1.40954, 2.11431) (0, 0.177451, 0.266176) 0.7 (0, 4.31924, 6.47886) (0, 2.16475, 3.24712) (0, 2.66512, 3.99767) (0, 0.669447, 1.00417) (0, 1.64447, 2.4667) (0, 0.207026, 0.310539) 0.8 (0, 4.93627, 7.40441) (0, 2.474, 3.71099) (0, 3.04585, 4.56877) (0, 0.765082, 1.14762) (0, 1.87939, 2.81909) (0, 0.236601, 0.354902) 0.9 (0, 5.5533, 8.32996) (0, 2.78325, 4.17487) (0, 3.42658, 5.13987) (0, 0.860717, 1.29108) (0, 2.11431, 3.17147) (0, 0.266176, 0.399265) 1.0 (0, 6.17034, 9.25551) (0, 3.09249, 4.63874) (0, 3.80731, 5.71096) (0, 0.956352, 1.43453) (0, 2.34924, 3.52386) (0, 0.295752, 0.443627) Example 2. Consider the following problem, FCDα ∗y(t) = (−1)y(t), 0 ≤t≤1, where y( 0 )=( 0, 1, 2 ) and the exact FC[ii −gH] -differentiable solution of problem is in the form y(t) = (0, 1, 2)Eα(−tα). In order to solve the mentioned problem the following formula should be applied as y0= (0, 1, 2), yk+1=ykgH hα Γ(α+1)yk,k=0, 1, ··· ,N−1, the numerical results based on the presented method are obtained in Table 2for various t , α= 0.3, 0.6, 0.9 and h=0.2, 0.02. The figures of exact solution and the Caputo gH-derivative are shown in Figure 2. Mathematics 2020,8, 2166 19 of 24 (a) (b) Figure 2. ( a ) The exact solution ( b ) the Caputo gH-derivative of solution defined in Example 2for α=0.3. Table 2. Numerical results of Example 2for various t,αand h. α=0.3 α=0.6 α=0.9 th=0.2 h=0.02 h=0.2 h=0.02 h=0.2 h=0.02 0.1 (0, 0.312475, 0.624949) (0, 0.655421, 1.31084) (0, 0.573896, 1.14779) (0, 0.892967, 1.78593) (0, 0.755737, 1.51147) (0, 0.969249, 1.9385) 0.2 (0, 0.0976404, 0.195281) (0, 0.429577, 0.859154) (0, 0.329356, 0.658712) (0, 0.797391, 1.59478) (0, 0.571138, 1.14228) (0, 0.939444, 1.87889) 0.3 (0, 0.0305101, 0.0610203) (0, 0.281554, 0.563107) (0, 0.189016, 0.378032) (0, 0.712044, 1.42409) (0, 0.43163, 0.863261) (0, 0.910555, 1.82111) 0.4 (0, 0.00953365, 0.0190673) (0, 0.184536, 0.369072) (0, 0.108476, 0.216951) (0, 0.635832, 1.27166) (0, 0.326199, 0.652398) (0, 0.882555, 1.76511) 0.5 (0, 0.00297902, 0.00595805) (0, 0.120949, 0.241898) (0, 0.0622536, 0.124507) (0, 0.567777, 1.13555) (0, 0.246521, 0.493042) (0, 0.855415, 1.71083) 0.6 (0, 0.000930869, 0.00186174) (0, 0.0792725, 0.158545) (0, 0.0357271, 0.0714542) (0, 0.507007, 1.01401) (0, 0.186305, 0.37261) (0, 0.829111, 1.65822) 0.7 (0, 0.000290873, 0.000581746) (0, 0.0519568, 0.103914) (0, 0.0205036, 0.0410072) (0, 0.45274, 0.905481) (0, 0.140797, 0.281595) (0, 0.803615, 1.60723) 0.8 (0, 0.0000908904, 0.000181781) (0, 0.0340536, 0.0681072) (0, 0.0117669, 0.0235339) (0, 0.404282, 0.808565) (0, 0.106406, 0.212812) (0, 0.778903, 1.55781) 0.9 (0, 0.000028401, 0.0000568019) (0, 0.0223195, 0.0446389) (0, 0.00675299, 0.013506) (0, 0.361011, 0.722022) (0, 0.0804149, 0.16083) (0, 0.754951, 1.5099) 1.0 (0,8.87458 × 10−6,0.0000177492) (0, 0.0146286, 0.0292573) (0, 0.00387551, 0.00775103) (0, 0.322371, 0.644742) (0, 0.0607725, 0.121545) (0, 0.731735, 1.46347) Example 3. Let us to consider the following problem, FCDα ∗y(t) = −π2t2−αα2 (2−3α+α2)Γ(1−α)pFq1; 3 2−α 2, 2 −α 2;−1 4π2t2α20, 1 2, 1, 1 ≤t≤2, where y( 1 ) = 0, 1 2, 1cos(απ) and the exact solution is y(t) = 0, 1 2, 1cos(απt) . We know that this problem has the switching point at t= 1.40426. According to Equation (20), we should divide the interval [ 1, 2 ] to the N subinterval [tk , tk+1] , for k= 0, 1, ..., N− 1, and assuming the switching point belongs to [tj , tj+1] , then the following iterative formulas are applied as yk+1=yk⊕hα Γ(α+1)−π2t2−αα2 (2−3α+α2)Γ(1−α)pFqa;b;zt2 k0, 1 2, 1, k=0, 1, ··· ,j, yk+1=yk(−1)hα Γ(α+1)−π2t2−αα2 (2−3α+α2)Γ(1−α)pFqa;b;zt2 k0, 1 2, 1, k=j+1, ··· ,N−1, where a= 1, b=3 2−α 2, 2 −α 2 and z=−1 4π2α2 . Numerical results are demonstrated in Table 3for α= 0.8 and h= 0.2, 0.02, 0.002. In Figure 3, the graphs of the exact solution and the Caputo gH-derivative are presented for α=0.8. Mathematics 2020,8, 2166 20 of 24 (a) (b) Figure 3. ( a ) The exact solution ( b ) the Caputo gH-derivative of solution defined in Example 3for α=0.8. Table 3. Numerical results of Example 3for α=0.8, h=0.2, 0.02, 0.002 and various t. α=0.8 t h =0.2 h=0.02 h=0.002 1.1 (0, −0.452376, −0.918699) (0, −0.463549, −0.928996) (0, −0.464888, −0.929776) 1.2 (0, −0.489654, −0.984721) (0, −0.495423, −0.991934) (0, −0.496057, −0.992115) 1.3 (0, −0.489654, −0.984721) (0, −0.495423, −0.991934) (0, −0.496057, −0.992115) 1.4 (0, −0.452376, −0.918699) (0, −0.463549, −0.928996) (0, −0.464888, −0.929776) 1.5 (0, −0.400112, −0.800241) (0, −0.404397, −0.808832) (0, −0.404508, −0.809017) 1.6 (0, −0.307754, −0.628932) (0, −0.317689, −0.637143) (0, −0.318712, −0.637424) 1.7 (0, −0.20878, −0.417784) (0, −0.21233, −0.425584) (0, −0.21289, −0.425779) 1.8 (0, −0.0927856, −0.176232) (0, −0.0935876, −0.187196) (0, −0.0936907, −0.187381) 1.9 (0, 0.0304523, 0.0618529) (0, 0.0313271, 0.0627734) (0, 0.0313953, 0.0627905) 2.0 (0, 0.146488, 0.301241) (0, 0.15359, 0.308805) (0, 0.154508, 0.309017) Example 4. Consider the following nonlinear fractional differential equations under uncertainty, √ηFCDα ∗y(t) + y2(t) = g(x), 0 <α<1, t∈[0, 1], where g(x) = Γ(6) Γ(6−α)t5−α−3Γ(5) Γ(5−α)t3−α+Γ(5) Γ(4−α)t3−α+ (t5−3t4+2t3)2˜ η, and ˜ η(r) = ( 1, 2, 3 ) , y6= 0. Then, the exact solution of the problem is y(t) = √η(t5− 3 t4+ 2 t3) . By solving the problem under FC[(i)-gH]-differentiability using yk+1=yk+hα Γ(α+1)Γk,k=0, 1, ··· ,N−1, we obtain the numerical solution shown in Table 4, with different order of differentiability and step size. Table 4. Absolute error of Example 4at t=1. hα=0.1 α=0.3 α=0.5 α=0.7 α=0.9 1 10 7.20602 ×10−26.5418 ×10−25.8823 ×10−25.3707 ×10−25.0201 ×10−2 1 20 3.9603 ×10−23.3498 ×10−22.9368 ×10−22.6937 ×10−22.5668 ×10−2 1 40 2.0653 ×10−21.6611 ×10−21.4420 ×10−21.3384 ×10−21.2962 ×10−2 1 80 1.0448 ×10−48.1038 ×10−37.0482 ×10−36.6381 ×10−36.5091 ×10−3 In Figure 4, the graphs of the exact solution and Caputo gH-derivatives are presented for α=0.1, 0.3, 0.5, 0.7, 0.9, 1, and in Figure 5these Caputo gH-derivatives have been compared in r =0.5. Mathematics 2020,8, 2166 21 of 24 (a) (b) (c) (d) (e) (f) (g) (h) Figure 4. ( a ) The exact solution and ( b – h ) the Caputo gH-derivatives of the solution defined in Example 4for α=0, 0.1, 0.3, 0.5, 0.7, 0.9, 1, respectively. Mathematics 2020,8, 2166 22 of 24 0.2 0.4 0.6 0.8 1.0 t -0.3 -0.2 -0.1 0.1 0.2 0.3 Α=1 Α=0.9 Α=0.7 Α=0.5 Α=0.3 Α=0.1 Α=0.0 0.2 0.4 0.6 0.8 1.0 t -0.4 -0.2 0.2 0.4 Α=1 Α=0.9 Α=0.7 Α=0.5 Α=0.3 Α=0.1 Α=0.0 (a) (b) Figure 5. The Caputo gH-derivatives ( a ) y(t) and ( b ) y(t) in r= 0.5 of solution defined in Example 4 for different values of α. Remark 1. Although we have obtained the solution under FC[(i)−gH] -differentiability, it is easy to check that it is not FC[(i)−gH] -differentiable on ( 0, 1 ) . Actually, due to obtained results (see Table 5), we can consider the proper interval that the given exact solution and its approximation is FC[(i)−gH] -differentiable. Moreover, note that, we have computed the approximation of the solution of Example 4at point t= 1, which is clearly this point take place out of proper domain of FC[(i)−gH] -differentiability. In fact, the computed error at point t= 1, just obtained based on the lower-upper approximation of lower-upper of exact solution. For more clarification, we determined switching points regarding each order of differentiability. Table 5. Switching points for different values of α. α0.1 0.3 0.5 0.7 0.9 1 t0.9701 0.9109 0.8525 0.7949 0.7381 0.7101 Remark 2. Indeed, using the results of Table 5, we in fact deduce that by considering the problem of fractional order instead of integer order (here, first order), we obtain some wider interval than the first order case; on the other hand, when α= 1, the valid interval that the given exact solution verify the assumption FC[(i)− gH] -differentiability is ( 0, 0.7101 ) , while for α= 0.9 and α= 0.7, the valid interval are ( 0, 0.7381 ) and ( 0, 0.7949 ) , respectively. Actually, this is the first time in the literature that this new result, i.e., extending the length of valid interval that the type of differentiability remains unchanged, is investigated. 8. Conclusions Fractional differential equations are one of the important topics of fuzzy arithmetic which have many applications in sciences and engineering. Thus finding the numerical and analytical methods to solve these problems is very important. This paper was presented based on the two main topics. First, proving the generalized Taylor series expansion for fuzzy valued function based on the concept of generalized Hukuhara differentiability. Second, introducing the fuzzy generalized Euler’s method as an application of the generalized Taylor expansion and applying it to solve the fuzzy fractional differential equations. The capabilities and abilities of the presented method were shown by presenting several theorems about the consistency, the convergence, and the stability of the generalized Euler’s method. Furthermore, the accuracy and efficiency of the method were illustrated by considering the local and global truncation errors. The numerical results especially in the switching point case showed the precision of the generalized Euler’s method to solve the fuzzy fractional differential equations. Mathematics 2020,8, 2166 23 of 24 Author Contributions: Conceptualization, T.A. and J.J.N.; methodology, Z.N.; software, S.N.; validation, T.A., J.J.N. and S.N.; formal analysis, T.A.; investigation, Z.N.; resources, Z.N.; data curation, S.N.; writing—original draft preparation, Z.N.; writing—review and editing, J.J.N.; visualization, J.J.N.; supervision, T.A.; project administration, T.A.; funding acquisition, J.J.N. All authors have read and agreed to the published version of the manuscript. Funding: The work of J.J. Nieto has been partially supported by Agencia Estatal de Investigación (AEI) of Spain under grant MTM2016-75140-P, co-financed by the European Community fund FEDER, and XUNTA de Galicia under grants GRC2015-004 and ED431C 2019/02. Acknowledgments: The authors are grateful to the anonymous reviewers for their helpful, valuable comments and suggestions in the improvement of this manuscript. Conflicts of Interest: The authors declare no conflict of interest.. References 1. Agarwal, R.P.; Lakshmikantham, V.; Nieto, J.J. On the concept of solution for fractional differential equations with uncertainty. Nonlinear Anal. Theory Methods Appl. 2010,72, 2859–2862. [CrossRef] 2. Bhaskara, T.G.; Lakshmikanthama, V.; Leela, S. Fractional differential equations with a krasnoselskii-krein type condition. Nonlinear Anal. Hybrid Syst. 2009,3, 734–737. [CrossRef] 3. Kilbas, A.A.; Srivastava, M.H.; Trujillo, J.J. Theory and Application of Fractional Differential Equations; Elsevier: Amsterdam, The Netherlands, 2006; Volume 204. ISBN-13 978-0-444-51832-3. 4. Miller, K.S.; Ross, B. An Introduction to the Fractional Calculus and Fractional Diffrential Equations; John Wiley and Sons: New York, NY, USA, 1993. 5. Noeiaghdam, Z.; Allahviranloo, T.; Nieto, J.J. q -fractional differential equations with uncertainty. Soft Comput. 2019,23, 9507–9524. [CrossRef] 6. Ngo, V.H. Fuzzy fractional functional integral and differential equations. Fuzzy Sets Syst. 2015 ,280, 58–90. [CrossRef] 7. Allahviranloo, T. Uncertain Information and Linear Systems. In Studies in Systems, Decision and Control; Springer: Berlin/Heidelberg, Germany, 2020; Volume 254, pp. 109–119, ISBN 978-3-030-31323-4. 8. Allahviranloo, T.; Gouyandeh, Z.; Armand, A. A full fuzzy method for solving differential equation based on Taylor expansion. J. Intell. Fuzzy Syst. 2015,29, 1039–1055. [CrossRef] 9. Chehlabi, M.; Allahviranloo, T. Concreted solutions to fuzzy linear fractional differential equations. Appl. Soft Comput. 2016,44, 108–116. [CrossRef] 10. Hoa, N.V. Fuzzy fractional functional differential equations under Caputo gH-differentiability. Commun. Nonlinear Sci. And Numerical Simul. 2015,22, 1134–1157. [CrossRef] 11. Agarwal, R.P.; Baleanu, D.; Nieto, J.J.; Torres, D.F.M.; Zhou, Y. A survey on fuzzy fractional differential and optimal control nonlocal evolution equations. J. Comput. Appl. Math. 2018,339, 3–29. [CrossRef] 12. Long, H.V.; Son, N.T.K.; Tam, H.T.T. The solvability of fuzzy fractional partial differential equations under Caputo gH-differentiability. Fuzzy Sets Syst. 2017,309, 35–63. [CrossRef] 13. Salahshour, S.; Allahviranloo, T.; Abbasbandy, S. Solving fuzzy fractional differential equations by fuzzy Laplace transforms. Commun. Nonlinear Sci. Numer. Simul. 2012,17, 1372–1381. [CrossRef] 14. Khodadadi, E.; Celik, E. The Variational Iteration Method For Fuzzy Fractional Differential Equations With Uncertainty. Fixed Point Theory Appl. 2013,2013, 1–13. [CrossRef] 15. Mazandarani, M.; Kamyad, A.V. Modified fractional Euler method for solving Fuzzy Fractional Initial Value Problem. Commun. Nonlinear Sci. Numer. Simulat. 2013,18, 12–21. [CrossRef] 16. Ahmadian, A.; Suleiman, M.; Salahshour, S.; Baleanu, D. A Jacobi Operational Matrix for Solving A Fuzzy Linear Fractional Differential Equation. Adv. Differ. Equations 2013,2013, 1–29. [CrossRef] 17. Masjed-Jamei, M.; Moalemi, Z.; Area, I.; Nieto, J.J. A new type of Taylor series expansion. J. Inequalities Appl. 2018, 116. [CrossRef] 18. Noeiaghdam, S.; Sidorov, D.; Sizikov, V.; Sidorov, N. Control of accuracy on Taylor-collocation method to solve the weakly regular Volterra integral equations of the first kind by using the CESTAC method. Appl. Comput. Math. 2020,19, 87–105. 19. Zarei, E.; Noeiaghdam, S. Solving generalized Abel’s integral equations of the first and second kinds via Taylor-collocation method. arXiv 2018, arXiv:1804.08571. Mathematics 2020,8, 2166 24 of 24 20. Allahviranloo, T.; Chehlabi, M. Solving fuzzy differential equations based on the length function properties. Soft Comput. 2015,19, 307–320. [CrossRef] 21. Allahviranloo, T.; Abbasbandy, S.; Sedaghgatfar, O.; Darabi, P. A new method for solving fuzzy integro-differential equation under generalized differentiability. Neural Comput. Appl. 2012 ,21, 191–196. [CrossRef] 22. Bede, B.; Stefanini, L. Generalized differentiability of fuzzy-valued functions. Fuzzy Set Syst. 2013 ,230, 119–141. [CrossRef] 23. Goetschel, R.; Voxman, W. Elementary fuzzy calculus. Fuzzy Sets Syst. 1987,24, 31–43. [CrossRef] 24. Ma, M.; Friedman, M.; Kandel, A. A new fuzzy arithmetic. Fuzzy Set. Syst. 1999,108, 83–90. [CrossRef] 25. Stefanini, L.; Bede, B. Generalized Hukuhara differentiability of interval-valued functions and interval differential equations. Nonlinear Anal. 2009,71, 1311–1328. [CrossRef] 26. Chakraverty, S.; Tapaswini, S.; Behera, D. Fuzzy Arbitrary Order System: Fuzzy Fractional Differential Equations and Applications; Wiley: Hoboken, NJ, USA, 2016. [CrossRef] 27. Anastassiou, G.A. Fuzzy mathematics: Approximation theory. Stud. Fuzziness Soft Comput. 2010 ,251, 267–271. [CrossRef] 28. Epperson, J.F. An Introduction to Numerical Methods and Analysis; John Wiley and Sons: Hoboken, NJ, USA, 2007. [CrossRef] 29. Gear, C.W. Numerical Initial Value Problems in Ordinary Differential Equations; Prentice Hall: Englewood Clifs, NJ, USA, 1971. Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).