An Operational Matrix Based on Legendre Polynomials for Solving Fuzzy Fractional-Order Differential Equations

Abstract

This paper deals with the numerical solutions of fuzzy fractional differential equations under Caputo-type fuzzy fractional derivatives of order 𝛼∈(0,1). We derived the shifted Legendre operational matrix (LOM) of fuzzy fractional derivatives for the numerical solutions of fuzzy fractional differential equations (FFDEs). Our main purpose is to generalize the Legendre operational matrix to the fuzzy fractional calculus. The main characteristic behind this approach is that it reduces such problems to the degree of solving a system of algebraic equations which greatly simplifies the problem. Several illustrative examples are included to demonstrate the validity and applicability of the presented technique.

1. Introduction

The subject of fractional calculus has gained considerable popularity and importance during the past three decades. Some of the most recent works on this topic, that is, the theory of derivatives and integrals of fractional (noninteger) order, are such as the book of Podlubny [1], Baleanu et al. [2], Diethelm [3], Baleanu et al. [4], and Sabatier et al. [5]. Only in the last few years, the various applications of fractional calculus have been extended in the area of physics and engineering such as the modeling of nonlinear oscillation of earthquake [6], the fluid-dynamic models [7], continuum and statistical mechanics [8], and solid mechanics [9]. In a notably enormous number of recent works, one can find the robustness upon the usefulness of fractional calculus to derive particular solutions of different kinds of classical differential equations like Bessel differential equation of general order [10, 11]. Also, the most significant advantage of applying FDEs is their nonlocal property, which interprets that the next state of a system relies not only pon its current phase but also pon all of its past records of phases [12]. For example, with the fractional differentiability, the fluid dynamic traffic model can get rid of the shortage arising from the hypothesis of continuum traffic flow [12, 13].

On the other hand, the modeling of natural phenomena is stated using mathematical tools (mathematical arithmetic, mathematical logics, etc.). However, obtaining a deterministic model of such problems is not easy, even does not occur and always has some errors and vagueness. So, investigating a popular way to interpret such vagueness is important. Since 1965 with Zadeh’s well-known paper on introducing fuzzy sets, applications of fuzzy concept to the structure of any modeling has appeared more and more, instead of deterministic case. So the topic of fuzzy differential equations (FDEs) of integer order has been rapidly growing in recent years [14–20]. Additionally, the application of various techniques has been expanded by means of the interpolations and polynomials for approximating the fuzzy solutions of fuzzy integral equations vastly, like Bernstein polynomials [21, 22], Lagrange interpolation [23, 24], Chebyshev interpolation [25], Legendre wavelets [26], and Galerkin-type technique [27].

Recently, Agarwal et al. [28] proposed the concept of solutions for the fractional differential equations with uncertainty. They have considered the Riemann-Liouville’s differentiability with a fuzzy initial condition to solve FFDEs. In [29, 30], the authors considered the generalization of H-differentiability for the fractional case. Discovering a suitable approximate or exact solution for FFDEs is a significant task which has been aroused simultaneously with the emerging of FFDES, except for a few number of these equations, and we have hardship in finding their analytical solutions. Consequently, there have been limited efforts to develop new methods for gaining approximate solutions which reasonably estimate the exact solutions. Salahshour et al. [31] considered fuzzy laplace transforms for solving FFDEs under Riemann-Liouville H-differentiability. Also Mazandarani and Kamyad [32] generalized the fractional Euler method for solving FFDEs under Caputo-type derivative.

From another point of view, several methods have been exploited to solve fractional differential equations, and fractional partial differential equations, fractional integrodifferential equations such as Adomian’s decomposition method [7], He’s variational iteration method [33], homotopy perturbation method [34], and spectral methods [35, 36]. In this way, orthogonal functions have received considerable attention in dealing with the various kinds of fractional differential equations. The main characteristic behind the approach using this technique is that it reduces these problems to those of solving a system of algebraic equations thus greatly simplifying the problem. Saadatmandi and Dehghan [36] presented the shifted Legendre operational matrix for fractional derivatives and applied it with tau method for the numerical solution of fractional differential equations subject to initial conditions. Also in [37–39], the authors derived new formulas using shifted Chebyshev polynomials and shifted Jacobi polynomials of any degree, respectively and applied them together with tau and collocation spectral methods for solving multiterm linear and nonlinear fractional differential equations.

The essential target of this paper is to recommend a suitable way to approximate FFDEs using a shifted Legendre tau approach. This strategy demands a formula for fuzzy fractional-order Caputo derivatives of shifted Legendre polynomials of any degree which is provided and applied together with the tau method for solving FFDEs with initial conditions. Up till now, and to the best of our knowledge, few methods corresponding to those mentioned previously have been devoted to solve FFDEs and are traceless in the literature for FFDEs under Caputo differentiability. This partially motivates our interest in the operational matrix of fuzzy fractional derivative of shifted Legendre polynomials. Also another motivation is based on the reality that only a few terms of expansion of the shifted Legendre function is needed to reach to a high accuracy, therefore, it does not need to implement the method frequently for finding the approximate results in each particular point.

For finding the fuzzy solution, the shifted Legendre operational matrix is generalized for the fuzzy fractional derivative (0<𝛼<1) which is based on the Legendre tau method for solving numerically FFDEs with the fuzzy initial conditions. It is worthy to note here that the method based on using the operational matrix of the Legendre orthogonal function for solving FFDEs is computer oriented.

The aim of this paper is to introduce the shifted Legendre operational matrix of fuzzy fractional derivative which is based on Legendre tau method for solving FFDEs under generalized differentiability. Also, we introduce a suitable way to estimate the nonlinear fuzzy fractional initial problems on the interval [0,1], by spectral shifted Legendre collocation method based on Legendre operational matrix, to find the approximate fuzzy solution. Finally, the accuracy of the proposed algorithms is demonstrated by several test problems. We note that the two shifted Legendre and shifted Jacobi operational matrices have been introduced by Saadatmandi and Dehghan [36] and Doha et al. [39], respectively, in the crisp concept. We, therefore, motivated our interest in the shifted Legendre operational matrix in the fuzzy settings.

This paper is organized as follows: In Section 2, we begin by introducing some necessary definitions and mathematical preliminaries of the fuzzy calculus and fractional calculus. Some basic concepts, properties and theorems of fuzzy fractional calculus are presented in Section 3. Section 4 is devoted to the fuzzy Legendre functions and their properties. The shifted Legendre operational matrix of fuzzy fractional derivative for solving fuzzy fractional differential equation is obtained in Section 5. Section 6 illustrates the effectiveness of the proposed method through solving several examples which some of them are modelled based on the real phenomena. Finally, a conclusion is given in the last section.

2. Preliminaries

We give some definitions and introduce the necessary notation which will be used throughout the paper, see, for example, [40, 41]. Also for some definitions related to generalized fuzzy difference, one can find more in [42, 43].

We denote the set of all real numbers by ℝ. A fuzzy number is a mapping 𝑢∶ℝ→[0,1] with the following properties:(a)𝑢 is upper semicontinuous,(b)𝑢 is fuzzy convex, that is, 𝑢(𝜆𝑥+(1−𝜆)𝑦≥min{𝑢(𝑥),𝑢(𝑦)} for all 𝑥,𝑦∈ℝ,𝜆∈[0,1],(c)𝑢 is normal, that is, ∃𝑥₀∈ℝ for which 𝑢(𝑥₀)=1,(d)supp𝑢={𝑥∈ℝ∣𝑢(𝑥)>0} is the support of the 𝑢, and its closure cl (supp 𝑢) is compact.

Let 𝔼 be the set of all fuzzy number on ℝ. The 𝛼-level set of a fuzzy number 𝑢∈𝔼,0≤𝑟≤1, denoted by [𝑢]_𝑟, is defined as [𝑢]_𝑟={𝑢(𝑥)}𝑢𝑥∈ℝ∣≥𝑟if0<𝑟≤1clsuppif𝑟=0. (1) It is clear that the 𝑟-level set of a fuzzy number is a closed and bounded interval [𝑢(𝑟),𝑢(𝑟)], where 𝑢(𝑟) denotes the left-hand endpoint of [𝑢]_𝑟 and 𝑢(𝑟) denotes the right-hand endpoint of [𝑢]_𝑟. Since each 𝑦∈ℝ can be regarded as a fuzzy number 𝑦 defined by 𝑦(𝑡)=1if𝑡=𝑦,0if𝑡≠𝑦, (2)ℝ can be embedded in 𝔼.

The addition and scalar multiplication of fuzzy number in 𝔼 are defined as follows: (1)𝑢⊕𝑣̃=𝑢+𝑣,̃𝑢+𝑣̃,(2𝑢)=𝜆𝑢)(𝜆⊙(𝑟),𝜆𝑢(𝑟)𝜆𝜆≥0,𝑢(𝑟),𝜆𝑢(𝑟)𝜆<0. (3)

The metric structure is given by the Hausdorff distance 𝐷∶𝔼×𝔼→ℝ₊⋃0, 𝐷(𝑢,𝑣)̃=sup_[]𝑟∈0,1max𝑢(𝑟)−𝑣̃(𝑟),𝑢(𝑟)−𝑣̃(𝑟). (4) It is easy to see that 𝐷 is a metric in 𝔼 and has the following properties (i)𝐷(𝑢⊕𝑤,𝑣̃⊕𝑤)=𝐷(𝑢,𝑣)̃, for all 𝑢,𝑣,̃𝑤∈𝔼,(ii)𝐷(𝑘⊙𝑢,𝑘⊙𝑣)̃=|𝑘|𝐷(𝑢,𝑣)̃, for all 𝑘∈ℝ,𝑢,𝑣̃∈𝔼,(iii)𝐷(𝑢⊕𝑣,̃𝑤⊕𝑒)̃≤𝐷(𝑢,𝑤)+𝐷(𝑣,̃𝑒)̃, for all 𝑢,𝑣,̃𝑤∈𝔼,(iv)𝐷(𝑢+𝑣,̃0)≤𝐷(𝑢,0)+𝐷(𝑣,0), for all 𝑢,𝑣∈𝔼,(v)(𝔼,𝐷) is a complete metric space.

Definition 1. The property (iv) in the properties of the above metric space suggests the definition of a function ‖⋅‖_𝔼∶𝑅→𝔼 that ‖𝑢‖_𝔼∶=𝐷(𝑢,0),forall𝑢∈𝔼 that has the properties of usual norms. In [44], the properties of this function are presented as follows:(i)‖𝑢‖≥0, for all 𝑢∈𝔼 and‖𝑢‖=0 if and only if 𝑢=0,(ii)‖𝜆⋅𝑢‖=|𝜆|⋅‖𝑢‖ and ‖𝑢+𝑣‖≥‖𝑢‖+‖𝑣‖, for all 𝑢,𝑣∈𝔼,forall𝜆∈ℝ.(iii)|‖𝑢‖−‖𝑣‖|≤𝐷(𝑢,𝑣) and 𝐷(𝑢,𝑣)≤‖𝑢‖+‖𝑣‖forall𝑢,𝑣∈𝔼.

Definition 2 (see [45]). Let 𝑓 and 𝑔 be the two fuzzy-number-valued functions on the interval [𝑎,𝑏], that is, 𝑓,𝑔∶[𝑎,𝑏]→𝔼. The uniform distance between fuzzy-number-valued functions is defined by 𝐷^∗𝑓,𝑔∶=sup_{[]𝑥∈𝑎,𝑏}𝐷𝑓(𝑥),𝑔(𝑥). (5)

Remark 3 (see [46]). Let 𝑓∶[𝑎,𝑏]→𝔼 be fuzzy continuous. Then from property (iv) of Hausdorff distance, we can define 𝐷𝑓(𝑥),0=sup_𝑟∈[0,1]max𝑓^𝑟(𝑥),𝑓^𝑟(𝑥),∀𝑥∈[𝑎,𝑏]. (6)

Definition 4 (see [42]). Let 𝐾^𝑛_𝑐 be the space of nonempty compact and convex sets of ℝ^𝑛. The generalized Hukuhara difference of two sets 𝐴,𝐵∈𝐾^𝑛_𝑐 (gH-difference for short) is defined as follows: 𝐴⊖_gH(a)(b)()𝐵=𝐶⟺𝐴=𝐵+𝐶or𝐵=𝐴+−1𝐶. (7) In case (a) of the above equation, the gH-difference is coincident with the H-difference. Thus the gH-difference is a generalization of the H-difference.

Definition 5 (see [47]). Let 𝑥,𝑦∈𝔼. If there exists 𝑧∈𝔼 such that 𝑥=𝑦⊕𝑧, and then 𝑧 is called the H-difference of 𝑥 and 𝑦, and it is denoted by 𝑥⊖𝑦.

In this paper, the sign “⊖” always stands for H-difference and note that 𝑥⊕𝑦≠𝑥+(−𝑦). Also throughout the paper is assumed that the Hukuhara difference and Hukuhara generalized differentiability existed.

Definition 6 (see [42]). The generalized difference (g-difference for short) of two fuzzy numbers 𝑢,𝑣∈𝔼 is given by its level sets as 𝑢⊖_g𝑣_𝛼=cl_𝛽≥𝛼[𝑢]_𝛽⊖_gH[𝑣]_𝛽∀𝛼∈[0,1], (8) where the gH-difference ⊖_gH is with interval operands [𝑢]_𝛽 and [𝑣]_𝛽.

Proposition 7. The g-difference in Definition 6 is given by the expression 𝑢⊖_g𝑣_𝛼=inf_𝛽≥𝛼min𝑢_𝛽−𝑣_𝛽,𝑢_𝛽−𝑣_𝛽,sup_𝛽≥𝛼max𝑢_𝛽−𝑣_𝛽,𝑢_𝛽−𝑣_𝛽. (9)

Proof. See  [42].

The next proposition gives simplified notation for 𝑢⊖_g𝑣 and 𝑣⊖_g𝑢.

Proposition 8. For any two fuzzy numbers 𝑢,𝑣∈𝔼 the two g-difference 𝑢⊖_g𝑣 and 𝑣⊖_g𝑢 exist and, for any 𝛼∈[0,1], one 𝑢⊖_g𝑣=−(𝑣⊖_g𝑢) with 𝑢⊖_g𝑣_𝛼=𝑑_𝛼,𝑑_𝛼,𝑢⊖_g𝑣_𝛼=−𝑑_𝛼,−𝑑_𝛼, (10) where 𝑑_𝛼=inf𝐷_𝛼,𝑑_𝛼=sup𝐷_𝛼, (11) and the sets 𝐷_𝛼 are 𝐷_𝛼=𝑢_𝛽−𝑣_𝛽∣𝛽≥𝛼∪𝑢_𝛽−𝑣_𝛽∣𝛽≥𝛼. (12)

Proof.  See [42].

The following proposition prove that the g-difference is well-defined.

Proposition 9 (see [14]). For any fuzzy numbers 𝑢,𝑣∈𝔼 the g-difference 𝑢⊖_g𝑣 exists and it is a fuzzy number.

Proof. See [42].

The following property holds for g-derivative.

Proposition 10. Let 𝑢,𝑣∈𝔼 be two fuzzy numbers, and then(i)𝑢⊖_g𝑣=𝑢⊖_gH𝑣 whenever the expressions on the right exist, in particular, 𝑢⊖_g𝑢=0, (ii)(𝑢+𝑣)⊖_g𝑣=𝑢, (iii)0⊖_g(𝑢⊖_g𝑣)=𝑣⊖_g𝑢, (iv)𝑢⊖_g𝑣=𝑣⊖_g𝑢=𝑤 if and only if 𝑤=−𝑤. Furthermore, 𝑤=0 if and only if 𝑢=𝑣.

In this paper, we consider the following definition which was introduced by Bede and Gal in [14].

Definition 11 (see [14]). Let 𝑓∶(𝑎,𝑏)→𝔼 and 𝑥₀∈(𝑎,𝑏). One says that 𝑓 is strongly generalized differentiable at 𝑥₀, if there exists an element 𝑓^′(𝑥)∈𝔼, such that(i)for all ℎ>0 sufficiently small,  ∃𝑓(𝑥₀+ℎ)⊖𝑓(𝑥₀),∃𝑓(𝑥₀)⊖𝑓(𝑥₀−ℎ), and the limits (in the metric 𝐷) lim_ℎ→0+𝑓𝑥₀+ℎ⊖𝑓𝑥₀ℎ=lim_ℎ→0+𝑓𝑥₀⊖𝑓𝑥₀−ℎℎ=𝑓^′𝑥₀, (13)(ii)for all ℎ>0 sufficiently small,  ∃𝑓(𝑥₀)⊖𝑓(𝑥₀+ℎ),∃𝑓(𝑥₀−ℎ)⊖𝑓(𝑥₀), and the limits (in the metric 𝐷) lim_ℎ→0+𝑓𝑥₀⊖𝑓𝑥₀+ℎ−ℎ=lim_ℎ→0+𝑓𝑥₀−ℎ⊖𝑓𝑥₀−ℎ=𝑓^′𝑥₀, (14)(iii)for all ℎ>0 sufficiently small,  ∃𝑓(𝑥₀+ℎ)⊖𝑓(𝑥₀),∃𝑓(𝑥₀−ℎ)⊖𝑓(𝑥₀), and the limits (in the metric 𝐷) lim_ℎ→0+𝑓𝑥₀+ℎ⊖𝑓𝑥₀ℎ=lim_ℎ→0+𝑓𝑥₀−ℎ⊖𝑓𝑥₀−ℎ=𝑓^′𝑥₀, (15)(iv)for all ℎ>0 sufficiently small,  ∃𝑓(𝑥₀)⊖𝑓(𝑥₀+ℎ),∃𝑓(𝑥₀)⊖𝑓(𝑥₀−ℎ), and the limits (in the metric 𝐷) lim_ℎ→0+𝑓𝑥₀⊖𝑓𝑥₀+ℎ−ℎ=lim_ℎ→0+𝑓𝑥₀⊖𝑓𝑥₀−ℎℎ=𝑓^′𝑥₀. (16)

Remark 12. Throughout this paper, we say that 𝑓 is (1)-differentiable on (𝑎,𝑏), if 𝑓 is differentiable in the sense (i) of Definition 11 and also 𝑓 is (2)-differentiable on (𝑎,𝑏), if 𝑓 is differentiable in the sense (ii) of Definition 11.

Theorem 13 (see [17]). Let 𝑓∶(𝑎,𝑏)→𝔼 be a function and denote [𝐹(𝑡)]^𝑟=[𝑓_𝑟(𝑡),𝑔_𝑟(𝑡)], for each 𝑟∈[0,1]. Then(1)if 𝑓 is (1)-differentiable, then 𝑓_𝑟(𝑡) and 𝑔_𝑟(𝑡) are differentiable functions and 𝐹^′(𝑡)^𝑟=𝑓^′_𝑟(𝑡),𝑔^′_𝑟(𝑡), (17)(2)if 𝑓 is (2)-differentiable, then 𝑓_𝑟(𝑡) and 𝑔_𝑟(𝑡) are differentiable functions and 𝐹^′(𝑡)^𝑟=𝑔^′_𝑟(𝑡),𝑓^′_𝑟(𝑡). (18)

Definition 14 (see [42]). Let 𝑓∶(𝑎,𝑏)→𝔼 and 𝑥₀∈(𝑎,𝑏). We say that 𝑓 is g-differentiable at 𝑥₀, if there exists an element 𝑓^′(𝑥₀)∈𝔼 such that 𝑓^′𝑥₀=lim_ℎ→0𝑓(𝑥+ℎ)⊖_g𝑓(𝑥)ℎ. (19)

Next we review one of the main results from Bede [15] for fuzzy initial value problem (FIVP) under (1)-differentiability which Nieto et al. [48] generalized this results for FIVP under (2)-differentiability (let ‖⋅‖ denote the usual Euclidean norm).

Theorem 15 (see [15], characterization theorem). Let one consider the fuzzy initial value problem 𝑦^′(𝑥),𝑦𝑡=𝑓𝑥,𝑦₀=𝑦₀, (20) where 𝑓∶[𝑥₀,𝑥₀+𝑎]×𝔼→𝔼 is such that(i)[𝑓(𝑥,𝑦)]^𝑟=[𝑓^𝑟(𝑥,𝑦,𝑦),𝑓^𝑟(𝑥,𝑦,𝑦)], (ii)𝑓^𝑟 and 𝑓^𝑟 are equicontinuous (i.e., for any 𝜖>0 there is a 𝛿>0 such that |𝑓^𝑟(𝑥,𝑦,𝑧)−𝑓^𝑟(𝑥₁,𝑦₁,𝑧₁)|<𝜖 and |𝑓^𝑟(𝑥,𝑦,𝑧)−𝑓^𝑟(𝑥₁,𝑦₁,𝑧₁)|<𝜖 for all 𝑟∈[0,1], whenever (𝑥,𝑦,𝑧),(𝑥₁,𝑦₁,𝑧₁)∈[𝑥₀,𝑥₀+𝑎]×ℝ² and ‖(𝑥,𝑦,𝑧)−(𝑥₁,𝑦₁,𝑧₁)‖<𝛿 and uniformly bounded on any bounded set,(iii)there exists an 𝐿>0 such that |𝑓^𝑟(𝑥₂,𝑦₂,𝑧₂)−𝑓^𝑟(𝑥₁,𝑦₁,𝑧₁)|≤𝐿max{|𝑦₂−𝑦₁|,|𝑧₂−𝑧₁|} for all 𝑟∈[0,1], |𝑓^𝑟(𝑥₂,𝑦₂,𝑧₂)−𝑓^𝑟(𝑥₁,𝑦₁,𝑧₁)|≤𝐿max{|𝑦₂−𝑦₁|,|𝑧₂−𝑧₁|} for all 𝑟∈[0,1].
Then the FIVP (20) and system of ODEs𝑦^𝑟(𝑥)^′=𝑓^𝑟𝑥,𝑦^𝑟,𝑦^𝑟,𝑦^𝑟(𝑥)^′=𝑓^𝑟𝑥,𝑦^𝑟,𝑦^𝑟,𝑦^𝑟𝑥₀=𝑦^𝑟₀,𝑦𝑥₀=𝑦^𝑟₀, (21)are equivalent.

Corollary 16 (see [48]). If we consider FIVP (20) under (2)-differentiability then the FIVP (20) and the following system of ODEs are equivalent: 𝑦^𝑟(𝑥)^′=𝑓^𝑟𝑥,𝑦^𝑟,𝑦^𝑟,𝑦^𝑟(𝑡)^′=𝑓^𝑟𝑥,𝑦^𝑟,𝑦^𝑟,𝑦^𝑟𝑥₀=𝑦^𝑟₀,𝑦𝑥₀=𝑦^𝑟₀. (22)

Theorem 17 (see [49]). Let 𝑓(𝑥) be a fuzzy-valued function on [𝑎,∞) and it is represented by (𝑓(𝑥;𝑟),𝑓(𝑥;𝑟)). For any fixed 𝑟∈[0,1], assume that (𝑓(𝑥;𝑟) and 𝑓(𝑥;𝑟)) are Riemann-integrable on [𝑎,𝑏] for every 𝑏≥𝑎, and assume that there are two positive 𝑀(𝑟) and 𝑀(𝑟) such that: ^𝑏_𝑎𝑓(𝑥;𝑟)𝑑𝑥≤𝑀(𝑟),^𝑏_𝑎𝑓(𝑥;𝑟)𝑑𝑥≤𝑀(𝑟)𝑓𝑜𝑟𝑒𝑣𝑒𝑟𝑦𝑏≥𝑎. (23) Then 𝑓(𝑥) is improper fuzzy Riemann-integrable on [𝑎,∞) and the improper fuzzy Riemann-integral is a fuzzy number. Further more, we have ^∞_𝑎𝑓(𝑥)𝑑𝑥=^∞_𝑎𝑓(𝑥;𝑟)𝑑𝑥,^∞_𝑎𝑓(𝑥;𝑟)𝑑𝑥. (24)

Definition 18 (see [45]). 𝑓(𝑥)∶[𝑎,𝑏]→𝔼. We say that 𝑓 Fuzzy-Riemann integrable to 𝐼∈𝔼, if for any 𝜖>0, there exists 𝛿>0 such that for any division 𝑃={[𝑢,𝑣];𝜉} of [𝑎,𝑏] with the norms Δ(𝑃)<𝛿, we have 𝐷_𝑃^∗(𝑢−𝑣)⊙𝑓(𝜁),𝐼<𝜖, (25) where ∑^∗ means addition with respect to ⊕ in 𝔼,𝐼∶=(𝐹𝑅)^𝑏_𝑎𝑓(𝑥)𝑑𝑥. (26) We also call an 𝑓 as above, (𝐹𝑅)-integrable.

Definition 19 (see [50]). Consider the 𝑛×𝑛 linear system of equations 𝑎₁₁𝑥₁+𝑎₁₂𝑥₂+⋯+𝑎_1𝑛𝑥_𝑛=𝑦₁,𝑎₂₁𝑥₁+𝑎₂₂𝑥₂+⋯+𝑎_2𝑛𝑥_𝑛=𝑦₂,⋮𝑎_𝑛1𝑥₁+𝑎_𝑛2𝑥₂+⋯+𝑎_𝑛𝑛𝑥_𝑛=𝑦_𝑛. (27) The matrix form of the above equations is 𝐴𝑋=𝑌, (28) where the coefficient matrix 𝐴=(𝑎_𝑖𝑗), 1≤𝑖, 𝑗≤𝑛 is a crisp 𝑛×𝑛 matrix and 𝑦_𝑖∈𝔼,1≤𝑖≤𝑛. This system is called a fuzzy linear system (FLS).

Definition 20 (see [50]). A fuzzy number vector (𝑥₁,𝑥₂,…,𝑥_𝑛)^𝑡 given by 𝑥_𝑖=(𝑥_𝑖(𝑟),𝑥_𝑖(𝑟)),1≤𝑖≤𝑛,0≤𝑟≤1 is called a solution of the fuzzy linear system (27) if ^𝑛_𝑗=1𝑎_𝑖𝑗𝑥_𝑗=^𝑛_𝑗=1𝑎_𝑖𝑗𝑥_𝑗=𝑦_𝑖,_𝑛𝑗=1𝑎_𝑖𝑗𝑥_𝑗=^𝑛_𝑗=1𝑎_𝑖𝑗𝑥_𝑗=𝑦_𝑖. (29) If for a particular 𝑘, 𝑎_𝑘𝑗>0,1≤𝑗≤𝑛, we simply get ^𝑛_𝑗=1𝑎_𝑘𝑗𝑥_𝑗=𝑦_𝑘,^𝑛_𝑗=1𝑎_𝑘𝑗𝑥_𝑗=𝑦_𝑘. (30) To solve fuzzy linear systems, one can refer to [51–53].

Now we define some notations which are used for the fuzzy fractional calculus throughout the paper.(i)𝐿^𝔼_𝑝(𝑎,𝑏),1≤𝑝<∞ is the set of all fuzzy-valued measurable functions 𝑓 on [𝑎,𝑏] where ||𝑓||_𝑝=(∫¹₀(𝑑(𝑓(𝑡),0))^𝑝𝑑𝑡)^1/𝑝.(ii)𝐶^𝔼[𝑎,𝑏] is a space of fuzzy-valued functions which are continuous on [𝑎,𝑏].(iii)𝐴𝐶^𝔼[𝑎,𝑏] denotes the set of all fuzzy-valued functions which are absolutely continuous.

One can easily find this definition in the crisp sense in [1, 54].

Definition 21 (see [54]). The Riemann-Liouville fractional integral operator of order 𝛼, 𝑛−1<𝛼≤𝑛, of a function 𝑓∈𝐶[𝑎,𝑏] is defined as 𝐼^𝛼_𝑎+𝑓(𝑥)=1Γ(𝛼)^𝑥_𝑎𝑓(𝑡)(𝑥−𝑡)^1−𝛼𝐼𝑑𝑡𝑥>𝑎,_𝛼𝑏−𝑓(𝑥)=1Γ(𝛼)^𝑏_𝑥𝑓(𝑡)(𝑡−𝑥)^1−𝛼𝑑𝑡𝑥<𝑏. (31)

Properties of the operator 𝐼^𝛼 can be found in [1, 54, 55]. we refer to only the following

For 𝑓∈𝐶[𝑎,𝑏],  𝛼,𝛽≥0 and 𝛾>−1, (1)𝐼^𝛼𝐼^𝛽𝑓(𝑥)=𝐼^𝛼+𝛽𝑓(𝑥),(2)𝐼^𝛼𝐼^𝛽𝑓(𝑥)=𝐼^𝛽𝐼^𝛼𝑓(𝑥),(3)𝐼^𝛼𝑥^𝛾=Γ𝛾+1Γ𝑥𝛼+𝛾+1^𝛼+𝛾. (32)

Definition 22 (see [54]). The Riemann-Liouville fractional derivatives of order 0<𝛼<1 of a function 𝑓∈𝐶[𝑎,𝑏] are expressed by 𝐷^𝛼_𝑎+𝑓(𝑥)=1Γ()𝑑1−𝛼𝑑𝑥^𝑥_𝑎𝑓(𝑡)𝑑𝑡(𝑥−𝑡)^𝛼(),𝐷𝑥>𝑎_𝛼𝑏−𝑓(𝑥)1=−Γ()𝑑1−𝛼1𝑑𝑥Γ(1−𝛼)^𝑥_𝑎𝑓(𝑡)𝑑𝑡(𝑡−𝑥)^𝛼(𝑥<𝑏). (33)

Definition 23 (see [55]). The fractional Caputo derivatives ^𝑐𝐷^𝛼_𝑎+𝑓(𝑥) and ^𝑐𝐷^𝛼_𝑏−𝑓(𝑥) on [𝑎,𝑏] for 0<𝛼<1 are defined via the above Riemann-Liouville fractional derivatives by ^𝑐𝐷^𝛼_𝑎+𝑓(𝑥)=𝐷^𝛼_𝑎+𝑓(𝑡)−𝑓(𝑎)(_𝑐𝑥),𝐷^𝛼_𝑏−𝑓(𝑥)=𝐷^𝛼_𝑏−𝑓(𝑡)−𝑓(𝑏)(𝑥), (34) which can be simplified as ^𝑐𝐷^𝛼_𝑎+𝑓(𝑥)=𝐷^𝛼_𝑎+𝑓(𝑥)−𝑓(𝑎)Γ()()1−𝑎𝑥−𝑎^−𝛼,_𝑐𝐷^𝛼_𝑏−𝑓(𝑥)=𝐷^𝛼_𝑏−𝑓(𝑥)−𝑓(𝑏)Γ()()1−𝑎𝑏−𝑥^−𝛼. (35)

Also, the fractional Caputo derivative can be defined in a sense of integral form described in Definition 24.

Definition 24 (see [56]). The Caputo definition of the fractional-order derivative is defined as ^𝑐𝐷^𝛼𝑓(𝑥)=1Γ(𝑛−𝛼)^𝑥₀𝑓^𝑛(𝑡)(𝑥−𝑡)^{𝛼+1−𝑛}𝑑𝑡,𝑛−1<𝛼≤𝑛,𝑛∈ℕ, (36) where 𝛼>0 is the order of the derivative and 𝑛 is the smallest integer greater than 𝛼. For the Caputo derivative, we have ^𝑐𝐷^𝛼𝐶=0,(𝐶isaconstant_𝑐),𝐷^𝛼𝑥^𝛽=0,for𝛽∈ℕ₀⌈𝛼⌉,Γ,𝛽<𝛽+1Γ𝑥𝛽+1−𝛼^𝛽−𝛼,for𝛽∈ℕ₀⌈𝛼⌉⌊𝛼⌋.,𝛽≥or𝛽∉ℕ,𝛽> (37)(38)

The ceiling function ⌈𝛼⌉ is used to denote the smallest integer greater than or equal to 𝛼, and the floor function ⌊𝛼⌋ to denote the largest integer less than or equal to 𝛼. Also ℕ={1,2,…} and ℕ₀={0,1,2,…}.

Definition 25 (see [1]). Similar to the differential equation of integer order, the Caputo’s fractional differentiation is a linear operation, that is, ^𝑐𝐷^𝛼𝜆𝑓(𝑥)+𝜇𝑔(𝑥)=𝜆^𝑐𝐷^𝛼𝑓(𝑥)+𝜇^𝑐𝐷^𝛼𝑔(𝑥), (39) where 𝜆 and 𝜇 are constants.

Lemma 26. Let 0<𝛼<1 and 𝑓∈𝐶[𝑎,𝑏]∩𝐿_𝑝[𝑎,𝑏]. Then the Caputo fractional derivatives are bounded for any 𝑥∈[𝑎,𝑏] and (1≤𝑝<∞) as ^𝑐𝐷^𝛼_𝑎+𝑓(𝑥≤𝑓)^′(𝑥)_𝑝|Γ(]()1−𝛼)|[1−𝛼𝑥−𝑎^1−𝛼=𝑀^𝛼_𝑎,_𝑐𝐷^𝛼_𝑏+𝑓𝑥≤𝑓()^′(𝑥)_𝑝|Γ(]()1−𝛼)|[1−𝛼𝑏−𝑥^1−𝛼=𝑀^𝛼_𝑏. (40)

Proof. See [54, 57].

3. Fuzzy Caputo Fractional Derivatives

In this section, some definitions and theorems related to the fuzzy Caputo fractional derivatives are presented which are an extension of the fractional derivative in the crisp sense. The generalized differentiability should be considered to expand the concept of Caputo fractional derivatives for the fuzzy space. For more details, see [14, 30].

Definition 27. Let 𝑓∶𝐿^𝔼∩𝐶^𝔼 be a fuzzy set-value function and Φ(𝑥)=(1/Γ(1−𝛼))∫^𝑥_𝑎(𝑓(𝑡)𝑑𝑡/(𝑥−𝑡)^𝛼), and then 𝑓 is said to be g-Caputo fuzzy fractional differentiable at 𝑥, when ^g𝐷^𝛼_𝑎+𝑓(𝑥)=lim_ℎ→0Φ(𝑥+ℎ)⊖_gΦ(𝑥)ℎ, (41) where ^g𝐷^𝛼_𝑎+𝑓()=1𝑥;𝑟Γ(1−𝛼)^𝑥_𝑎𝑓^′(𝑡;𝑟)𝑑𝑡(𝑥−𝑡)^𝛼_g,𝐷^𝛼_𝑎+𝑓()=1𝑥;𝑟Γ(1−𝛼)^𝑥_𝑎𝑓^′(𝑡;𝑟)𝑑𝑡(𝑥−𝑡)^𝛼. (42)

Remark 28. A fuzzy-valued function 𝑓 is ^𝐶[1−𝛼]-differentiable, if it is differentiable as in Definition 27, Case (i), and it is ^𝐶[2−𝛼]-differentiable, if it is differentiable as in Definition 27, case (ii).

Theorem 29. Let 0<𝛼<1 and 𝑓∈𝐴𝐶^𝔼[𝑎,𝑏], then Caputo fuzzy fractional derivative exists almost everywhere on (𝑎,𝑏) and for all 0≤𝑟≤1 we have ^𝐶𝐷^𝛼_𝑎+𝑓)=1(𝑥;𝑟Γ(1−𝛼)^𝑥_𝑎𝑓^′(𝑡)(𝑥−𝑡)^𝛼1𝑑𝑡,Γ(1−𝛼)^𝑥_𝑎𝑓^′(𝑡)(𝑥−𝑡)^𝛼𝑑𝑡 (43)

when 𝑓 is (1)-differentiable, and ^𝐶𝐷^𝛼_𝑎+𝑓)=1(𝑥;𝑟Γ(1−𝛼)^𝑥_𝑎𝑓^′(𝑡)(𝑥−𝑡)^𝛼1𝑑𝑡,Γ(1−𝛼)^𝑥_𝑎𝑓^′(𝑡)(𝑥−𝑡)^𝛼𝑑𝑡 (44)

when 𝑓 is (2)-differentiable, in which (𝐼^𝛼_𝑎+𝑓)(𝑥)=(1/Γ(𝛼))∫^𝑥_𝑎(𝑓(𝑡)𝑑𝑡/(𝑥−𝑡)^1−𝛼) for 𝑥>𝑎.

Proof. It is straightforward by applying Definitions 20 and 24.

Theorem 30. Let one assume that 𝑓∈𝐶^𝔼[𝑎,𝑏], and then one has the following: 𝐼^𝛼_𝑎+^𝐶𝐷^𝛼_𝑎+𝑓(𝑥)=𝑓(𝑥)⊖𝑓(𝑎),0<𝛼<1, (45) when 𝑓 is ^𝐶[1−𝛼]-differentiable and 𝐼^𝛼_𝑎+^𝐶𝐷^𝛼_𝑎+𝑓(𝑥)=−𝑓(𝑎)⊖−𝑓(𝑥),0<𝛼<1, (46) when 𝑓 is ^𝐶[2−𝛼]-differentiable.

Proof. See [30].

Lemma 31. Let 0<𝛼<1 and 𝑓∈𝐴𝐶^𝔼[𝑎,𝑏], then the fuzzy Caputo derivative can be expressed by means of the fuzzy Riemann-Liouville integral as follows: ^𝐶𝐷^𝛼_𝑎+𝑓)=𝐼(𝑥;𝑟_{𝑎1−𝛼+}()=𝐼𝐷𝑓𝑥;𝑟_{𝑎1−𝛼+}𝐷𝑓(𝑥;𝑟),𝐼_𝑎^1−𝛼₊𝐷𝑓(𝑥;𝑟), (47) when 𝑓 is (1)-differentiable, and ^𝐶𝐷^𝛼_𝑎+𝑓)=𝐼(𝑥;𝑟_{𝑎1−𝛼+}()=𝐼𝐷𝑓𝑥;𝑟_{𝑎1−𝛼+}𝐷𝑓(𝑥;𝑟),𝐼_𝑎^1−𝛼₊𝐷𝑓(𝑥;𝑟), (48) when 𝑓 is (2)-differentiable.

Now we consider the generalization of Taylor’s formula for the fuzzy Caputo fractional derivative which was introduced in [32, 45]. It should be mentioned that this theorem is the extension of Taylor’s formula for the Caputo fractional derivative in the crisp context [58].

Theorem 32. Let 𝑓(𝑥)∈𝐴𝐶^𝔼(0,𝑏] and suppose that ^𝑐𝐷^𝑘𝛼𝑓(𝑥)∈𝐶^𝔼(0,𝑏] for 𝑘=0,1,…,𝑛+1 where 0<𝛼<1, 0≤𝑥₀≤𝑥 and 𝑥∈(0,𝑏]. Then we have 𝑓(𝑥)^𝑟=𝑓^𝑟(𝑥),𝑓^𝑟(𝑥,𝑓_𝑟)(𝑥)=^𝑛_𝑖=0𝑥^𝑖𝛼Γ(𝑖𝛼+1)^𝑐𝐷^𝑖𝛼𝑓^𝑟0⁺+^𝑐𝐷^(𝑛+1)𝛼𝑓^𝛼𝑥₀Γ()𝑥𝑛𝛼+𝛼+1^(𝑛+1)𝛼,𝑓^𝑟(𝑥)=^𝑛_𝑖=0𝑥^𝑖𝛼Γ(𝑖𝛼+1)^𝑐𝐷^𝑖𝛼𝑓^𝑟0⁺+^𝑐𝐷^(𝑛+1)𝛼𝑓^𝛼𝑥₀Γ()𝑥𝑛𝛼+𝛼+1^(𝑛+1)𝛼, (49) where ^𝑐𝐷^𝛼𝑓^𝑟(0)=^𝑐𝐷^𝛼𝑓^𝑟(𝑥)|_𝑥=0,^𝑐𝐷^𝛼𝑓^𝑟(0)=^𝑐𝐷^𝛼𝑓^𝑟(𝑥)|_𝑥=0.

Proof. See [32].

Now, characterization theorem (Theorem 15), which was introduced by Bede in [15] and established by Pederson and Sambandham in [59] for hybrid fuzzy differential equations, is extended for fuzzy Caputo-type fractional differential equations. To this end, we first consider the FFDEs under Caputo’s H-differentiability for 0<𝛼<1 as follows: ^𝑐𝐷^𝛼_𝑎+𝑦(𝑥)=𝑓𝑥,𝑦(𝑥),𝑦(𝑎)∈𝔼. (50)

Theorem 33 (characterization theorem). Let one consider the fuzzy fractional differential equation under Caputo’s H-differentiability (50) where 𝑓∶[𝑥₀,𝑥₀+𝑎]×𝔼→𝔼 and such that:(i)[𝑓(𝑥,𝑦)]^𝑟=[𝑓^𝑟(𝑥,𝑦,𝑦),𝑓^𝑟(𝑥,𝑦,𝑦)], (ii)𝑓^𝑟 and 𝑓^𝑟 are equicontinuous (i.e., for any 𝜖>0 there is a 𝛿>0 such that |𝑓^𝑟(𝑥,𝑦,𝑧)−𝑓^𝑟(𝑥₁,𝑦₁,𝑧₁)|<𝜖 and |𝑓^𝑟(𝑥,𝑦,𝑧)−𝑓^𝑟(𝑥₁,𝑦₁,𝑧₁)|<𝜖 for all 𝑟∈[0,1], whenever (𝑥,𝑦,𝑧),(𝑥₁,𝑦₁,𝑧₁)∈[𝑥₀,𝑥₀+𝑎]×ℝ² and ‖(𝑥,𝑦,𝑧)−(𝑥₁,𝑦₁,𝑧₁)‖<𝛿 and uniformly bounded on any bounded set,(iii) there exists an 𝐿>0 such that |𝑓^𝑟(𝑥₂,𝑦₂,𝑧₂)−𝑓^𝑟(𝑥₁,𝑦₁,𝑧₁)|≤𝐿max{|𝑦₂−𝑦₁|,|𝑧₂−𝑧₁|}forall𝑟∈[0,1], |𝑓^𝑟(𝑥₂,𝑦₂,𝑧₂)−𝑓^𝑟(𝑥₁,𝑦₁,𝑧₁)|≤𝐿max{|𝑦₂−𝑦₁|,|𝑧₂−𝑧₁|}forall𝑟∈[0,1].
Then, (50) and the following system of FDEs are equivalent when 𝑦(𝑥) is ^𝐶[1−𝛼]-differentiable ^𝑐𝐷^𝛼_𝑥⁺₀𝑦(𝑥;𝑟)=𝑓^𝑟𝑥,𝑦^𝑟,𝑦^𝑟,_𝑐𝐷^𝛼_𝑥⁺₀𝑦(𝑥;𝑟)=𝑓^𝑟𝑥,𝑦^𝑟,𝑦^𝑟,𝑦^𝑟𝑥₀=^𝑐𝑦^𝑟₀,𝑦𝑥₀=^𝑐𝑦^𝑟₀, (51) also (50) and the following system of FDEs is equivalent when 𝑦(𝑥) is ^𝐶[2−𝛼]-differentiable ^𝑐𝐷^𝛼_𝑥⁺₀𝑦(𝑥;𝑟)=𝑓^𝑟𝑥,𝑦^𝑟,𝑦^𝑟,_𝑐𝐷^𝛼_𝑥⁺₀𝑦(𝑥;𝑟)=𝑓^𝑟𝑥,𝑦^𝑟,𝑦^𝑟,𝑦^𝑟𝑥₀=^𝑐𝑦^𝑟₀,𝑦𝑥₀=^𝑐𝑦^𝑟₀. (52)

Proof. In the papers [15, 59], the authors proved for fuzzy ordinary differential equations and hybrid fuzzy differential equations. The result for FFDEs is obtained analogously by using Theorem 2 in [15] and Theorem  3.1 in [59].

4. Proposed Method for Solving FFDEs

Saadatmandi and Dehghan [36] introduced the shifted Legendre operational matrix for derivative with fractional order using a spectral method which has been followed by Doha et al. [37–39]. They presented the shifted Chebyshev polynomials and Jacobi polynomials for solving fractional differential equations by tau method. In this section, we try to approximate fuzzy solution using shifted Legendre polynomials under H-differentiability as follows.

4.1. Properties of Shifted Legendre Polynomials

The Legendre polynomials, denoted by 𝑃_𝑛(𝑧), are orthogonal with Legendre weight function: 𝑤(𝑧)=1 over [−1,1], namely [60], ¹₋₁𝑃_𝑛(𝑧)𝑃_𝑚(𝑧)2𝑑𝑧=𝛿2𝑛+1^𝑛𝑚, (53) where 𝛿_𝑛𝑚 is the Kronecker function and can be specified with the help of following recurrence formula: 𝑃₀(𝑧)=1,𝑃₁(𝑧)𝑃=𝑧,_𝑖+1(𝑧)=2𝑖+1𝑖+1𝑧𝑃_𝑖(𝑧)−𝑖𝑃𝑖+1^𝑖−1(𝑧),𝑖=1,2,…. (54) In order to use these polynomials on the interval [0,1], we define the so-called shifted Legendre polynomials by introducing the change of variable 𝑧=2𝑥−1. Let the shifted Legendre polynomials 𝑃_𝑛(2𝑥−1) be denoted by 𝐿_𝑛(𝑥). The shifted Legendre polynomials are orthogonal with respect to the weight function 𝑤_𝑠(𝑥)=1 in the interval (0,1) with the following orthogonality property: ¹₀𝐿_𝑛(𝑥)𝐿_𝑚(𝑥)1𝑑𝑥=𝛿2𝑛+1^𝑛𝑚. (55) The shifted Legendre polynomials are generated from the following three-term recurrence relation: 𝐿_𝑖+1(𝑡)=(2𝑖+1)(2𝑥−1)𝐿𝑖+1^𝑖(𝑥)−𝑖𝐿𝑖+1^𝑖−1(𝑥),𝐿𝑖=1,2,…₀(𝑥)=1,𝐿₁(𝑥)=2𝑥−1. (56) The analytic form of the shifted Legendre polynomial 𝐿_𝑛(𝑥) of degree 𝑛 is given by 𝐿_𝑛(𝑥)=^𝑛_𝑖=0(−1)^𝑛+𝑖(𝑛+𝑖)!()!𝑥𝑛−𝑖(𝑖!^𝑖)²=^𝑛_𝑖=0𝑒_𝑖,𝑛𝑥^𝑖, (57) in which 𝑒_𝑖,𝑛=(−1)^𝑛+𝑖(𝑛+𝑖)!(𝑛−𝑖)!(𝑖!)², (58) where 𝐿_𝑛(0)=(−1)^𝑛,𝐿_𝑛(1)=1. (59) A function 𝑢(𝑥) of independent variable defined for 0≤𝑥≤1 may be expanded in terms of shifted Legendre polynomials as 𝑢(𝑥)=^𝑚_𝑖=0𝑏_𝑖𝐿_𝑖(𝑥)=Λ^𝑇Φ(𝑥), (60) where the shifted Legendre coefficients matrix 𝐵 and the shifted Legendre vector Φ(𝑥) are given by Λ^𝑇=𝜆₁,𝜆₂,…,𝜆_𝑚,Φ(𝑥)=𝐿₀(𝑥),𝐿₁(𝑥),…,𝐿_𝑚(𝑥)^𝑇. (61) Also, the derivative of Φ(𝑥) can be expressed by 𝑑Φ(𝑥)(𝑥),𝑑𝑥=𝐷Φ (62) where 𝐷 is the (𝑚+1)×(𝑚+1) operational matrix of derivative given by 𝐷=𝑑_𝑖𝑗=22𝑗+1,for𝑗=𝑖−𝑘,𝑘=1,3,…,𝑚,if𝑚odd,𝑘=1,3,…,𝑚−1,if𝑚even,0,otherwise. (63) The graph of some shifted Legendre polynomials (for 3≤𝑚≤8) shown in Figure 1 to depict their behaviors.

Figure 1 

The shifted Legendre functions for different 𝑚.

Now, we use the shifted Legender functions due to approximate a fuzzy function.

4.2. The Approximation of Fuzzy Function

In this section, we propose a shifted Legendre approximation for the fuzzy-valued functions. To this end, we use Legendre’s nodes and fuzzy shifted Legendre polynomials to calculate the fuzzy best approximation. For more details, see [61–67].

Definition 34. For 𝑦∈𝐿^𝔼_𝑝(0,1)∩𝐶^𝔼(0,1) and Legendre polynomial 𝐿_𝑛(𝑥) a real valued function over (0,1), the fuzzy function is approximated by 𝑦(𝑥)=^∞∗_𝑗=0𝑐_𝑗⊙𝐿_𝑗(𝑥),𝑥∈(0,1), (64) where the fuzzy coefficients 𝑐_𝑗 are obtained by 𝑐_𝑗=2𝑗+1⊙¹₀𝑦(𝑥)⊙𝐿_𝑗(𝑥)𝑑𝑥, (65) in which 𝐿_𝑗(𝑥) is the same as in (57), and ∑^∗ means addition with respect to ⊕ in 𝔼.

Remark 35. In actuality, only the first (𝑚+1)-terms shifted Legendre polynomials are considered. So we have 𝑦(𝑥)≃𝑦_𝑚(𝑥)=^𝑚∗_𝑗=0𝑐_𝑗⊙𝐿_𝑗(𝑥)=𝐶^𝑇_𝑚⊙Φ_𝑚(𝑥), (66) hence 𝑦^𝑟(𝑥)≃𝑦^𝑟_𝑚(𝑥)=^𝑚∗_𝑗=0𝑐^𝑟_𝑗⊙𝐿_𝑗(𝑥), (67) that the fuzzy shifted Legendre coefficient vector 𝐶^𝑇_𝑚+1 and shifted Legendre function vector Φ_𝑚+1(𝑥) are defined as 𝐶^𝑇_𝑚=𝑐₀,𝑐₂,…,𝑐_𝑚,Φ_𝑚(𝑥)=𝐿₀(𝑥),𝐿₁(𝑥),…,𝐿_𝑚(𝑥). (68)(69)

Definition 36 (see [68]). A fuzzy-valued polynomial 𝑝^∗∈∏_𝑁 is the best approximation to fuzzy function 𝑓 on 𝜒={𝑥₀,𝑥₁,𝑥₂,…,𝑥_𝑁}, if max_{𝑖=0,1,2,…,𝑁}𝐷𝑝^∗𝑥_𝑖,𝑓_𝑖𝑝=min_∈∏_𝑁max_{𝑖=0,1,2,…,𝑁}𝐷𝑝𝑥_𝑖,𝑓_𝑖, (70) in which ∏_𝑁 is the set of all fuzzy valued polynomials.

The problem is referred to as the best shifted Legendre approximation, as we use Legendre’s nodes.

Theorem 37. The best approximation of a fuzzy function based on the Legendre nodes exists and is unique.

Proof. The proof is an instantaneous outcome of Theorem  4.2.1 in [68].

Now, we want to show that the fuzzy approximation converges of Legendre functions to function 𝑦(𝑥).

Lemma 38. Suppose that 𝑦(𝑥)∈𝐴𝐶^𝔼(0,1]∩𝐿^𝔼_𝑝(0,1] and ^𝑐𝐷^𝛼𝑦(𝑥)∈𝐶^𝔼(0,1], 0<𝛼<1, 0≤𝑥₀≤𝑥, and 𝑥∈(0,1], and also assume that 𝑦(𝑥) is ^𝑐[1−𝛼]-differentiable. Therefore using Theorem 32, we have 𝑦^𝑟(𝑥)−^𝑚_𝑖=0𝑥^𝑖𝛼Γ()𝐷𝑖𝛼+1^𝑖𝛼𝑦^𝑟0⁺≤𝑀_𝛼𝑥^(𝑚+1)𝛼Γ(),𝑚𝛼+𝛼+1𝑦^𝑟(𝑥)−^𝑚_𝑖=0𝑥^𝑖𝛼Γ()𝐷𝑖𝛼+1^𝑖𝛼𝑦^𝑟0⁺≤𝑀_𝛼𝑥^(𝑚+1)𝛼Γ(),𝑚𝛼+𝛼+1 (71) in which |𝐷^(𝑚+1)𝛼𝑦^𝑟(𝑥₀)|<𝑀_𝛼 and |𝐷^(𝑚+1)𝛼𝑦^𝑟(𝑥₀)|<𝑀_𝛼.

Proof. From Theorem 32, we have 𝑦^𝑟(𝑥)=^𝑚_𝑖=0𝑥^𝑖𝛼Γ(𝑖𝛼+1)^𝑐𝐷^𝑖𝛼𝑦^𝑟0⁺+^𝑐𝐷^(𝑚+1)𝛼𝑦^𝑟𝑥₀Γ()𝑥𝑚𝛼+𝛼+1^(𝑚+1)𝛼,𝑦^𝑟(𝑥)=^𝑚_𝑖=0𝑥^𝑖𝛼Γ(𝑖𝛼+1)^𝑐𝐷^𝑖𝛼𝑦^𝑟0⁺+^𝑐𝐷^(𝑚+1)𝛼𝑦^𝑟𝑥₀Γ()𝑥𝑚𝛼+𝛼+1^(𝑚+1)𝛼, (72) and the following relation can be obtained: 𝑦^𝑟(𝑥)−^𝑚_𝑖=0𝑥^𝑖𝛼Γ()𝐷𝑖𝛼+1^𝑖𝛼𝑦^𝑟0⁺≤^𝑐𝐷^(𝑚+1)𝛼𝑦^𝑟𝑥₀Γ()𝑥𝑚𝛼+𝛼+1^(𝑚+1)𝛼≤𝑀^𝑟_𝛼𝑥^(𝑚+1)𝛼Γ(),𝑚𝛼+𝛼+1𝑦^𝑟(𝑥)−^𝑚_𝑖=0𝑥^𝑖𝛼Γ()𝐷𝑖𝛼+1^𝑖𝛼𝑦^𝑟0⁺≤^𝑐𝐷^(𝑚+1)𝛼𝑦^𝑟𝑥₀Γ()𝑥𝑚𝛼+𝛼+1^(𝑚+1)𝛼≤𝑀^𝑟_𝛼𝑥^(𝑚+1)𝛼Γ(),𝑚𝛼+𝛼+1 (73) that 𝑀^𝑟_𝛼=|^𝑐𝐷^(𝑚+1)𝛼𝑦^𝑟(𝑥₀)| and 𝑀^𝑟_𝛼=|^𝑐𝐷^(𝑚+1)𝛼𝑦^𝑟(𝑥₀)|.

Remark 39. If we consider Lemma 38 and define 𝑀_𝛼=max{𝑀^𝑟_𝛼,𝑀^𝑟_𝛼}, then (71) can be stated in the following form, regarding to Section 2: 𝐷𝑦(𝑥),^𝑚_𝑖=0𝑥^𝑖𝛼Γ()𝐷𝑖𝛼+1^𝑖𝛼𝑦0⁺=sup_[]𝑟∈0,1max𝑦^𝑟(𝑥)−^𝑚_𝑖=0𝑥^𝑖𝛼Γ()𝐷𝑖𝛼+1^𝑖𝛼𝑦^𝑟0⁺,𝑦^𝑟(𝑥)−^𝑚_𝑖=0𝑥^𝑖𝛼Γ()𝐷𝑖𝛼+1^𝑖𝛼𝑦^𝑟0⁺≤sup_[]𝑟∈0,1max𝑀^𝑟_𝛼𝑥^(𝑛+1)𝛼Γ(),𝑛𝛼+𝛼+1𝑀^𝑟_𝛼𝑥^(𝑛+1)𝛼Γ(𝑛𝛼+𝛼+1)≤𝑀_𝛼𝑥^(𝑛+1)𝛼Γ().𝑛𝛼+𝛼+1 (74)

Theorem 40. Let ^𝑐𝐷^𝛼𝑦(𝑥)∈𝐶^𝔼(0,1]∩𝐿^𝔼_𝑝(0,1], and 0<𝛼<1. Also consider a sequence of finite dimensional fuzzy space 𝑋_𝔼⊂𝑋, 𝑚≥1, in which 𝑋_𝔼 have dimension 𝑑_𝑚+1. Additionally, 𝑋_𝔼 have a basis {𝐿_𝑖(𝑥)}^𝑚_𝑖=0. If one assumes that 𝑦_𝑚(𝑥)=𝐶^𝑇Φ is the best fuzzy approximation for fuzzy function 𝑦(𝑥) from {𝐿_𝑖(𝑥)}^𝑚_𝑖=0, then the error estimation is as follows: lim_𝑚→∞𝐷𝑦(𝑥),𝑦_𝑚(𝑥)=0. (75)

Proof. Let 𝑓∶𝑅→𝔼 be a fuzzy valued function such that [𝑓(𝑥)]^𝑟=[𝑓^𝑟(𝑥),𝑓^𝑟(𝑥)]. Also consider the fuzzy Taylor’s formula in Theorem 32, 𝑓^𝑟(𝑥)=∑^𝑚_𝑖=0(𝑥^𝑖𝛼/Γ(𝑖𝛼+1))^𝑐𝐷^𝑖𝛼𝑓^𝑟(0⁺) and 𝑓^𝑟(𝑥)=∑^𝑚_𝑖=0(𝑥^𝑖𝛼/Γ(𝑖𝛼+1))^𝑐𝐷^𝑖𝛼𝑓^𝑟(0⁺). From Lemma 38 we have 𝑦^𝑟(𝑥)−𝑓^𝑟(𝑥)≤𝑀^𝑟_𝛼𝑥^(𝑚+1)𝛼Γ(),𝑚𝛼+𝛼+1𝑦^𝑟(𝑥)−𝑓^𝑟(𝑥)≤𝑀^𝑟_𝛼𝑥^(𝑚+1)𝛼Γ().𝑚𝛼+𝛼+1 (76) From the assumption, 𝐶^𝑇Φ is the best fuzzy approximation to 𝑦 from {𝐿_𝑖(𝑥)}^𝑚_𝑖=0, and 𝑓∈𝑋_𝔼, 𝑥∈(0,1]. So one has 𝑦^𝑟(𝑥)−𝑦^𝑟_𝑚(𝑥)≤𝑦^𝑟(𝑥)−𝑓^𝑟(𝑥)≤𝑀^𝑟_𝛼𝑥^(𝑚+1)𝛼Γ(),𝑚𝛼+𝛼+1𝑦^𝑟(𝑥)−𝑦^𝑟_𝑚(𝑥)≤𝑦^𝑟(𝑥)−𝑓^𝑟(𝑥)≤𝑀^𝑟_𝛼𝑥^(𝑚+1)𝛼Γ(),𝑚𝛼+𝛼+1 (77) and thus, taking into account Theorem  1 in [69] and above relations, we have 𝑦^𝑟(𝑥)−𝑦^𝑟_𝑚(𝑥)²≤𝑦^𝑟(𝑥)−𝑓^𝑟(𝑥)²≤𝑀_𝛼^𝑟,2Γ(𝑚𝛼+𝛼+1)²¹⁰𝑥^{2(𝑚+1)𝛼}≤𝑀𝑑𝑥_𝛼𝑟,2Γ(𝑚𝛼+𝛼+1)²(2(𝑚+1)𝛼+1)²,𝑦^𝑟(𝑥)−𝑦^𝑟_𝑚(𝑥)²≤𝑦^𝑟(𝑥)−𝑓^𝑟(𝑥)²≤𝑀_𝛼^𝑟,2Γ(𝑚𝛼+𝛼+1)²¹⁰𝑥^{2(𝑚+1)𝛼}≤𝑑𝑥𝑀_𝛼^𝑟,2Γ(𝑚𝛼+𝛼+1)²(2(𝑚+1)𝛼+1)². (78) From Remark 39 and Lemma 38, we have 𝑦^𝑟(𝑥)−𝑦^𝑟_𝑚(𝑥≤𝑀)^𝛼Γ(𝑚𝛼+𝛼+1)1(2()),𝑚+1𝛼+1𝑦^𝑟(𝑥)−𝑦^𝑟_𝑚(𝑥≤𝑀)^𝛼Γ(𝑚𝛼+𝛼+1)1(2()),𝑚+1𝛼+1 (79) and if 𝑚→∞, we get ‖𝑦^𝑟(𝑥)−𝑦_𝑚^𝑟(𝑥)‖→0, ‖𝑦^𝑟(𝑥)−𝑦^𝑟_𝑚(𝑥)‖→0. Therefore, from Remark 39 and the definition of Hausdorff distance in Section 2, it can be implied that lim_𝑚→∞𝐷𝑦(𝑥),𝑦_𝑚(𝑥)=0, (80) which completes the proof.

Remark 41. The same result can be obtained for 𝑦(𝑥) under ^𝑐[2−𝛼]-differentiability.

4.3. Operational Matrix of Caputo Fractional Derivative

Lemma 42. The fuzzy Caputo fractional derivative of order 0<𝛼<1 over the shifted Legendre functions can be gained in the form of ^𝑐𝐷^𝛼𝐿_𝑖(𝑥)=^𝑖_𝑘=0𝑒^′_𝑘,𝑖Γ(𝑘+1)Γ()𝑥𝑘−𝛼+1^𝑘−𝛼, (81) where 𝑒^′_𝑘,𝑖=0 when 𝑖<⌈𝛼⌉ and for 𝑖≥⌈𝛼⌉, and so one has 𝑒^′_𝑘,𝑖=𝑒_𝑘,𝑖.

Proof. Employing the analytic form of shifted Legendre polynomials explained in Section 4.1 and Definition 25, we have: 𝐷^𝛼𝐿_𝑖(𝑥)=^𝑖_𝑘=0𝑒_𝑘,𝑖𝐷^𝛼⊙𝑥^𝑘−𝛼. (82) Now, by exploiting Definition 24, the lemma can be proved.

Lemma 43. Let 0<𝛼<1, and the integral of the product of the fuzzy Caputo fractional derivative with order 𝛼 over the shifted legendre functions can be obtained by ¹₀𝐷^𝛼𝐿_𝑖(𝑥)𝐿_𝑗(𝑥)𝑑𝑥=^𝑖𝑗_𝑘=0𝑙=0𝑒^′_𝑘,𝑖𝑒_𝑙,𝑗()Γ()𝑘+𝑙+1−𝛼𝑘+1Γ().𝑘−𝛼+1 (83)

Proof. Using Lemma 42 and the analytic form of shifted Legendre polynomials explained in Section 4.1, we can acquire: ¹₀𝐷^𝛼𝐿_𝑖(𝑥)𝐿_𝑗(𝑥)=𝑑𝑥_{𝑖𝑗𝑘=0𝑙=0}𝑒^′_𝑘,𝑖𝑒_𝑙,𝑗Γ(𝑘+1)Γ(𝑘−𝛼+1)¹₀𝑥^{(𝑘+𝑙+1)−𝛼−1}𝑑𝑥. (84)
The operational matrix of different orthogonal functions for solving differential equations was introduced in the crisp concept [36, 37, 39]. Here, the Legendre operational matrix (LOM) in [36] is applied to the FFDEs using Caputo-type derivative.
The Caputo fractional derivatives operator of order 0<𝛼<1 of the vector Φ defined in (69) can be stated by 𝐷^𝛼Φ(𝑥)≃𝐷^(𝛼)Φ(𝑥), (85) where 𝐷^(𝛼) is the (𝑚+1)-square operational matrix of fractional Caputo’s derivative of Legendre functions. Regarding the following theorem, the LOM elements 𝐷^(𝛼)_𝑖,𝑗 are determined under Caputo fractional derivative. This theorem is generalizing the operational matrix of derivatives of shifted Legendre given in Section 4.1 to the fractional calculus.

Theorem 44. Let Φ be Legendre functions vector. 𝐷^(𝛼)_𝑖,𝑗 is the (𝑚+1)-square operational matrix of fractional Caputo’s derivative of order 0<𝛼<1. Then the elements of 𝐷^(𝛼)_𝑖,𝑗 are achieved as 𝐷^(𝛼)_𝑖,𝑗=^𝑖_{⌈𝛼⌉𝑘=}𝜃_{𝑖,𝑗,𝑘}, (86) in which 𝜃_{𝑖,𝑗,𝑘} are acquired by 𝜃_{𝑖,𝑗,𝑘}=×2𝑗+1_𝑗𝑙=0(−1)^{𝑖+𝑗+𝑘+𝑙}(𝑖+𝑘)!𝑙+𝑗!(𝑖−𝑘)!𝑘!Γ(𝑘−𝛼+1)𝑗−𝑙!(𝑙!)²().𝑘+𝑙−𝛼+1 (87) Consider that in 𝐷^(𝛼), the first ⌈𝛼⌉ rows, are all zero.

Proof. Employing the relation (85) and the orthogonal properties of shifted Legendre functions (56), we have 𝐷^(𝛼)=𝐷^𝛼Φ(𝑥),Φ(𝑥)^𝑇𝐸⁻¹, (88) in which ⟨𝐷^𝛼Φ(𝑥),Φ(𝑥)^𝑇⟩ and 𝐸⁻¹ are (𝑚+1)-square matrixes defined as 𝐷^𝛼Φ(𝑥),Φ(𝑥)^𝑇=¹₀𝐷^𝛼𝐿_𝑖(𝑥)𝐿_𝑗(𝑥)𝑑𝑥_𝑖,𝑗=0,𝐸⁻¹=diag{(2𝑖+1)}^𝑚−1_𝑖=0. (89) Hence, applying Lemma 43 and inserting the above matrixes in the product 𝐷^(𝛼), the theorem be proved.

Remark 45. If 𝛼=𝑛∈ℕ, then Theorem 44 gives the mentioned result as in Section 4.1.
The following property of the product of two shifted Legendre polynomials fuzzy vectors will also be applied which is the extension of the crisp case introduced in [36], ΦΦ^𝑇Λ≃^∼ΛΦ, (90) where ^∼Λ is an (𝑚+1)-square product operational matrix for the fuzzy vector Λ=[Λ,Λ]. Using the above equation and by the orthogonal property equation (56), the elements {^∼Λ_𝑖𝑗}^𝑚_𝑖,𝑗=0 can be computed from ^∼Λ_𝑖𝑗=2𝑗+1^𝑚∗_𝑘=0𝜆_𝑘⊙𝑔_𝑖𝑗𝑘, (91) where ∑^∗ denotes the fuzzy summation and ⊙ indicates fuzzy multiplication. Moreover, 𝑔_𝑖𝑗𝑘 are obtained by 𝑔_𝑖𝑗𝑘=¹₀𝐿_𝑖(𝑥)𝐿_𝑗(𝑥)𝐿_𝑘(𝑥)𝑑𝑥, (92) that in the simpler form, it is given by 𝑔_𝑖𝑗𝑘=𝑑_𝑗−𝑙𝑑_𝑙𝑑_𝑖−𝑙2𝑖+2𝑗−2𝑙+1𝑑_{𝑖+𝑗−𝑙},𝑘=𝑖+𝑗−2𝑙;𝑙=0,1,…,𝑗,0,𝑘≠𝑖+𝑗−2𝑙;𝑙=0,1,…,𝑗, (93) in which 𝑑_𝑙=(2𝑙)!/2^𝑙(𝑙!)².

Remark 46. As a matter of fact, from (93), it is obvious that the elements of the product of shifted Legendre operational matrix are independent from the fuzzy vector.

5. Application of the LOM for Solving FFDEs

In this section, Legendre operational matrix of fractional Caputo’s derivative is considered to present its significance for solving fuzzy fractional differential equation which is alike to the concept of the Caputo-type derivative in crisp case [36]. The existence and uniqueness of the solutions under fuzzy Caputo’s type derivative to this problem are discussed in [30].

5.1. General Linear FFDEs

Let us consider the general linear fuzzy fractional differential equation ^𝑐𝐷^𝛼₀⁺𝑦(𝑥)=𝑎₁⊙𝑦(𝑥)⊕𝑎₂(𝑥)𝑦(0)⊙𝑓,0<𝛼≤1,=𝑦₀∈𝔼, (94) in which 𝑎_𝑗 for 𝑗=1,2 are fuzzy constant coefficients, 𝑦(𝑥)∶𝐿^𝔼∩𝐶^𝔼 is a continuous fuzzy-valued function, and ^𝑐𝐷^𝛼₀⁺ represents the fuzzy Caputo fractional derivative of order 𝛼.

Remark 47. In this section, considering more simplicity, it is assumed that for the fuzzy set-valued function 𝑦(𝑥)∶𝐼⊂ℝ→𝐿^𝔼∩𝐶^𝔼, we present the functions 𝑦^𝑟,𝑦^𝑟∶𝑅→𝐼⊂𝑅,𝑟∈[0,1] by 𝑦^𝑟₋(𝑥)=(𝑦(𝑥))^𝑟₋,𝑦^𝑟₊(𝑥)=(𝑦(𝑥))^𝑟₊,forall𝑥∈𝐼,forall𝑟∈[0,1], which are in the same previous manner, called the left and right 𝑟-cut functions of 𝑦(𝑥).

Let 𝑟∈[0,1] and 0<𝛼<1. Then from (94), we have that 𝑎₁⊙𝑦(𝑥)^𝑟=𝑎₁⊙𝑦(𝑥)^𝑟₋,𝑎₁⊙𝑦(𝑥)^𝑟₊=𝑎₁^𝑟₋𝑦(𝑥)^𝑟₋,𝑎₁^𝑟₊𝑦(𝑥)^𝑟₊=𝑎₁^(𝑟)𝑦(𝑥)^(𝑟),𝑎₁𝑦(𝑥)^𝑟_±=𝑎₁^𝑟_±𝑦(𝑥)^𝑟_±. (95) Therefore, ^𝑐𝐷^𝛼₀⁺𝑦^𝑟_±(𝑥)=𝑎₁^𝑟_±𝑦(𝑥)^𝑟_±+𝑎₂^𝑟_±𝑓(𝑥)^𝑟_±, (96) and thus ^𝑐𝐷^𝛼₀⁺𝑦^𝑟₋(𝑥),^𝑐𝐷^𝛼₀⁺𝑦^𝑟₊(𝑥=𝑎_1𝑟−)𝑦(𝑥)^𝑟₋+𝑎₂^𝑟₋𝑓(𝑥)^𝑟₋,𝑎₁^𝑟₊𝑦(𝑥)^𝑟₊+𝑎₂^𝑟₊𝑓(𝑥)^𝑟₊. (97) Hence, we obtain ^𝑐𝐷^𝛼₀⁺𝑦(𝑥)^(𝑟)=^𝑐𝐷^𝛼₀⁺𝑦^𝑟₋(𝑥),^𝑐𝐷^𝛼₀⁺𝑦^𝑟₊(𝑥=𝑎_1𝑟−)𝑦(𝑥)^𝑟₋+𝑎₂^𝑟₋𝑓(𝑥)^𝑟₋,𝑎₁^𝑟₊𝑦(𝑥)^𝑟₊+𝑎₂^𝑟₊𝑓(𝑥)^𝑟₊=𝑎₁⊙𝑦(𝑥)^(𝑟)+𝑎₂⊙𝑓(𝑥)^(𝑟)[],=𝑎,𝑟∈0,1₁⊙𝑦(𝑥)⊕𝑎₂⊙𝑓(𝑥)^(𝑟),𝑟∈[0,1]. (98) We can rewrite (94) in the operator form 𝑎₁^(𝑟)𝐼⊖_g^𝑐𝐷^𝛼𝑦^(𝑟)=𝑎₂^(𝑟)𝑓^(𝑟), (99) in which the fuzzy operator ^𝑐𝐷^𝛼=𝐼^1−𝛼𝐷 is supposed to be compact on a fuzzy Banach space 𝑋_𝔼 to 𝑋_𝔼 [45, 70] and ⊖_g is the notation for g-difference.

For solving fuzzy fractional differential equation (94), we attempt to find a fuzzy function 𝑦_𝑚^(𝑟)∈𝑋_𝔼; therefore, let (^𝑐𝐷^𝛼₀⁺𝑦)(𝑥), 𝑦(𝑥) and 𝑓(𝑥) be approximated using Definition 34 as: 𝑦(𝑥)≃𝑦_𝑚(𝑥)=^𝑚∗_𝑗=0𝑐_𝑗⊙𝐿_𝑗(𝑥)=𝐶^𝑇_𝑚⊙Φ_𝑚, (100) that 𝑦_𝑚(𝑥)^(𝑟)=^𝑚∗_𝑗=0𝑐_𝑗^(𝑟)⊙𝐿_𝑗(𝑥𝑓(𝑥)≃𝑓)𝑥∈𝐼⊂𝑅,_𝑚(𝑥)=^𝑚∗_𝑗=0𝑓_𝑗⊙𝐿_𝑗(𝑥)=𝐹^𝑇_𝑚⊙Φ_𝑚, (101) where 𝐹_𝑚+1=[𝑓₀,𝑓₁,…,𝑓_𝑚]^𝑇 is obtained as 𝑓_𝑗=2𝑗+1⊙¹₀𝑓(𝑥)⊙𝐿_𝑗(𝑥)𝑑𝑥. (102) Also, using relations (85) and (100), we obtain ^𝑐𝐷^𝛼𝑦(𝑥)≃𝐶^𝑇⊙𝐷^𝛼Φ_𝑚(𝑥)≃𝐶^𝑇⊙𝐷^(𝛼)Φ_𝑚(𝑥). (103) Substituting (100)–(103) in problem (94), the coefficients {𝑐_𝑗^(𝑟)}^𝑚_𝑗=0 are specified by imposing the equation to be almost fuzzy exact in the Legendre operational matrix form. Now, we establish fuzzy residual for the approximation of (94), when [𝑦(𝑥)]^(𝑟)≈[𝑦_𝑚(𝑥)]^(𝑟).

Theorem 48. Let 𝑦^𝑟∈𝐶^𝔼[0,1] and 0<𝛼≤1, and then (^𝑐𝐷^𝛼𝑦_𝑚)(𝑥)^(𝑟)=𝑅_𝑚(𝑥)⊕𝑎₂⊙𝑓_𝑚(𝑥)⊕𝑎₁⊙𝑦_𝑚(𝑥)^(𝑟). (104)

Proof. Let 𝑟∈[0,1]. We have [(^𝑐𝐷^𝛼₀⁺𝑦_𝑚)(𝑥)]^(𝑟)=[(^𝑐𝐷^𝛼₀⁺𝑦_𝑚)^𝑟₋(𝑥),(^𝑐𝐷^𝛼₀⁺𝑦_𝑚)^𝑟₊(𝑥)], [𝑦_𝑚(𝑥)]^(𝑟)=[𝑦^(𝑟)_𝑚−(𝑥),𝑦^(𝑟)_𝑚+(𝑥)] and regarding to what was conversed previously, ^𝑐𝐷^𝛼𝑦^(𝑟)_𝑚±(𝑥)=𝑅^(𝑟)_𝑚±(𝑥)−𝑎₁^𝑟_±𝑦_𝑚(𝑥)^𝑟_±−𝑎₂^𝑟_±𝑓_𝑚(𝑥)^𝑟_±, (105) so ^𝑐𝐷^𝛼𝑦_𝑚(𝑥)^(𝑟)=^𝑐𝐷^𝛼𝑦^(𝑟)_𝑚−(𝑥),^𝑐𝐷^𝛼𝑦^(𝑟)_𝑚+(𝑥=𝑅_{(𝑟)𝑚−})(𝑥),𝑅^(𝑟)_𝑚+(𝑥−𝑎_1𝑟−)𝑦_𝑚(𝑥)^𝑟₋,𝑎₁^𝑟₊𝑦_𝑚(𝑥)^𝑟₊−𝑎₂^𝑟₋𝑓_𝑚(𝑥)^𝑟₋,𝑎₂^𝑟₊𝑓_𝑚(𝑥)^𝑟₊=𝑅_𝑚(𝑥)^(𝑟)⊕(𝑎₁)⊙𝑦_𝑚(𝑥)^(𝑟)⊕(𝑎₂)⊙𝑓_𝑚(𝑥)^(𝑟). (106) Using g-difference, we have 𝑅_𝑚(𝑥)^(𝑟)=(^𝑐𝐷^𝛼𝑦_𝑚)(𝑥)^(𝑟)⊖_g(𝑎₁)⊙𝑦_𝑚(𝑥)^(𝑟)⊖_g(𝑎₂)⊙𝑓_𝑚(𝑥)^(𝑟), (107) or in the sense of fuzzy operator, 𝑅_𝑚^(𝑟)=𝑎₁^(𝑟)𝐼⊖_g^𝑐𝐷^𝛼𝑦^(𝑟)⊖_g𝑎₂^(𝑟)𝑓^(𝑟). (108)

It is anticipated that the deriving fuzzy function [𝑦_𝑚(𝑥)]^(𝑟) will be a suitable approximation of the exact solution [𝑦(𝑥)]^(𝑟). To this end, let 𝑋_𝔼=𝐿²_𝔼([0,1]), and let ⟨⋅,⋅⟩_𝔼 indicate the fuzzy inner product for 𝑋_𝔼. It is demanded that 𝑅_𝑚^(𝑟) satisfy 𝑅_𝑚^(𝑟),𝐿_𝑖𝐸=0,𝑖=0,1,…,𝑚−1,𝑟∈[0,1], (109) in which ⟨𝑅_𝑚^(𝑟),𝐿_𝑖⟩_𝐸=[(𝐹𝑅)∫¹₀𝑅_𝑚(𝑥)⊙𝐿_𝑖(𝑥)𝑑𝑥]^(𝑟). The left side is the shifted legendre coefficients associated with 𝐿_𝑗. If {𝐿_𝑖}^𝑚_𝑖=0 are the main members of shifted Legendre family Φ={𝐿_𝑖}^∞_𝑖=0 which is complete in 𝑋_𝔼, then (109) needs the main terms to be zero in the Fourier extension of 𝑅_𝑚^(𝑟) with respect to Φ which is called tau method in the crisp context and it is in a similar manner with the meaning of fuzzy (for more details, see [71–73]).

To discover 𝑦_𝑛^(𝑟), implementing (109) to (108), regarding the (100)–(103), we generate 𝑚 fuzzy linear equations as 𝑅_𝑚(𝑥)^(𝑟),𝐿_𝑖(𝑥=(_𝑐)𝐷^𝛼𝑦_𝑚)(𝑥)^(𝑟)⊖_g(𝑎₁)⊙𝑦_𝑚(𝑥)^(𝑟)⊖_g(𝑎₂)⊙𝑓_𝑚(𝑥)^(𝑟),𝐿_𝑖(𝑥)=0, (110) for 𝑖=0,1,…,𝑚−1. Equation (110) generate (𝑚) set of fuzzy linear equations. These fuzzy linear equations can be acquired as the following system of equations: ^𝑚−1∗_𝑗=0𝑐_𝑗^(𝑟)⊙𝐷^(𝛼)𝐿_𝑗,𝐿_𝑖−𝑎₁^(𝑟)𝐿_𝑗,𝐿_𝑖=^𝑚−1∗_𝑗=0𝑓_𝑗^(𝑟)⊙𝑎₂^(𝑟)𝐿_𝑗,𝐿_𝑖,𝑖=0,1,…,𝑚−1 (111) or ^𝑚−1∗_𝑗=0𝑐_𝑗^(𝑟)⊙(𝐹𝑅)¹₀𝐷^(𝛼)𝐿_𝑗(𝑥)𝐿_𝑖(𝑥)−()𝑑𝑡𝐹𝑅¹₀𝑎₁^(𝑟)𝐿_𝑗(𝑥)𝐿_𝑖(𝑥)=𝑑𝑡_{𝑚−1∗𝑗=0}𝑓_𝑗^(𝑟)⊙(𝐹𝑅)¹₀𝑎₂^(𝑟)𝐿_𝑗(𝑥)𝐿_𝑖(𝑥)𝑑𝑡. (112) Afterwards, substitution of (100) in the initial condition of (94) yields 𝑦(0)=^𝑚∗_𝑗=0𝑐_𝑗^(𝑟)⊙𝐿_𝑗(0)=𝑦₀, (113) that this equation be coupled with the previous fuzzy linear equations and constructed (𝑚+1) fuzzy linear equations. Clearly, after solving this fuzzy system, the coefficients {𝑐_𝑗}^𝑚_𝑗=0 will be gained.

5.1.1. Error Analysis

The aim of this section is to acquire the error bound for the Legendre approximation using shifted Legendre polynomials. We consider the best shifted Legendre approximation of a smooth fuzzy function under Caputo derivative for 0<𝛼≤1 to reach the result. It should be considered that the results of this section are the extension of this concept in the crisp context (see more in [69, 74–76]). Initially, we state the following lemma which can offer an upper bound for approximating the error of Caputo fractional derivative. So we define the error vector 𝐸_𝛼 as 𝐸_𝛼=𝐷^𝛼Φ−𝐷^(𝛼)Φ=𝐸_0,𝛼,𝐸_1,𝛼,…,𝐸_𝑚,𝛼^𝑇, (114) where 𝐸_𝑖,𝛼=𝐷^𝛼𝐿_𝑖(𝑥)−^𝑚_𝑗=0𝐷^(𝛼)_𝑖𝑗𝐿_𝑗(𝑥),𝑖=0,1,…,𝑚. (115)

Now, we can propose the next lemma by using Theorem  1 in [69] to depict the error bound of Caputo fractional derivative operator for the shifted Legendre polynomials.

Lemma 49. Let the error function of Caputo fractional derivative operator for Legendre polynomials 𝐸_𝑖,𝛼∶[𝑥₀,1]→ℝ be 𝑚+1 times continuously differentiable for 0<𝑥₀≤𝑥, 𝑥∈(0,1]. Also 𝐸_𝑖,𝛼∈𝐶_𝑚+1[𝑥₀,1] and 𝛼<𝑚+1 and then the error bound is presented as follows: 𝐸_𝑖,𝛼≤|Γ(𝑖+1)||Γ(𝑀1−𝛼)|Γ(𝑚𝛼+^𝛼𝛼+1)1(2())𝑥𝑚+1𝛼+1^0−𝛼. (116)

Proof. Firstly, we acquire a bound for ^𝑐𝐷^𝛼𝑥^𝑖, regarding Definition 24, as ^𝑐𝐷^𝛼𝑥^𝑖=Γ(𝑖+1)Γ()𝑥𝑖+1−𝛼^𝑖−𝛼≤|Γ(𝑚+2)||Γ(𝑥1−𝛼)|^0−𝛼,𝑖=0,1,…,𝑚. (117) Then using Lemma  1 in [69], we have: 𝐷^𝛼𝐿_𝑖(𝑥)=^𝑖_𝑘=0𝑒_𝑘,𝑖𝐷^𝛼𝑥^𝑘≤|Γ(𝑖+1)||Γ(𝑥1−𝛼)|^0𝑖−𝛼_𝑘=0𝑒_𝑘,𝑖=|Γ(𝑖+1)||Γ(𝑥1−𝛼)|^0−𝛼𝐿_𝑖(1)=|Γ(𝑖+1)||Γ(𝑥1−𝛼)|^0−𝛼, (118) in which 𝑒_𝑘,𝑖=(−1)^𝑘+𝑖(𝑘+𝑖)!(𝑖−𝑘)!(𝑘!)². (119) Now, utilizing Theorem  1 in [69], the lemma can be proved.

Therefore, the maximum norm of error vector 𝐸_𝛼 is attained as 𝐸_𝛼∞≤|Γ(𝑚+1)||Γ(𝑀1−𝛼)|Γ(𝑚𝛼+^𝛼𝛼+1)1(2())𝑥𝑚+1𝛼+1^0−𝛼. (120)

Now we extend these results for the fuzzy case. Hence, the error of fuzzy Caputo fractional derivative is defined as 𝐷^∗𝐷^𝛼Φ(𝑥),𝐷^(𝛼)Φ(𝑥)=sup_[]𝑥∈0,1𝐷𝐷^𝛼Φ(𝑥),𝐷^(𝛼)Φ(𝑥)=sup_[]𝑥∈0,1sup_[]𝑟∈0,1max𝐷^𝛼Φ(𝑥)^𝑟₋(𝑥)−𝐷^(𝛼)Φ(𝑥)^𝑟₋(𝑥,𝐷_𝛼)Φ(𝑥)^𝑟₊(𝑥)−𝐷^(𝛼)Φ(𝑥)^𝑟₊(𝑥)=sup_[]𝑟∈0,1max𝐷^𝛼Φ(𝑥)^𝑟₋(𝑥)−𝐷^(𝛼)Φ(𝑥)^𝑟₋(𝑥)_∞,𝐷^𝛼Φ(𝑥)^𝑟₊(𝑥)−𝐷^(𝛼)Φ(𝑥)^𝑟₊(𝑥)_∞ (121) regarding Lemmas 38 and 49 and Theorem 40, and we have: sup_[]𝑟∈0,1max𝐷^𝛼Φ(𝑥)^𝑟₋−𝐷^(𝛼)Φ(𝑥)^𝑟_−∞,𝐷^𝛼Φ(𝑥)^𝑟₊(𝑥)−𝐷^(𝛼)Φ(𝑥)^𝑟₊(𝑥)_∞Lemma49≤sup_[]𝑟∈0,1|Γ(max𝑚+1)||Γ(𝑀1−𝛼)|Γ(^𝑟𝛼,−)×𝑚𝛼+𝛼+11(2())𝑥𝑚+1𝛼+1^0−𝛼,|Γ(𝑚+1)||Γ(𝑀1−𝛼)|Γ(^𝑟𝛼,+)×𝑚𝛼+𝛼+11(2())𝑥𝑚+1𝛼+1^0−𝛼, (122) where 0<𝑥₀≤𝑥,  𝑥∈(0,1].

Lemma 50. Consider the error function of the fuzzy Caputo fractional derivative operator for shifted Legendre polynomials which is continuously differentiable and 0<𝛼<1. Then the error bound is given as follows: 𝐷^∗𝐷^𝛼Φ(𝑥),𝐷^(𝛼)Φ(𝑥)≤|Γ(𝑚+1)||Γ(𝑀1−𝛼)|Γ(𝑚𝛼+^𝛼𝛼+1)1(2())𝑥𝑚+1𝛼+1^0−𝛼, (123) that 𝑀_𝛼=max{𝑀_𝛼,−,𝑀_𝛼,+}.

Proof. It is straightforward from the definition of (𝐸_𝛼), Definition 1, and Remark 3.

For an error assessment of the approximation solution of (94), we assume that 𝑦_𝑚(𝑥) and 𝑦(𝑥) reveal the fuzzy approximate and exact solutions of the fuzzy fractional differential equations, respectively. Then we rewrite (94) as ^𝑐𝐷^𝛼𝑦^(𝑟)_𝑚±(𝑥)−𝑎_1±^(𝑟)𝑦_𝑚(𝑥)_±^(𝑟)=𝐻^(𝑟)_𝑚±(𝑥)+𝑎_2±^(𝑟)𝑓_𝑚(𝑥)_±^(𝑟), (124) that 𝐻^(𝑟)_𝑚±(𝑥) is the fuzzy perturbation function that depends only on 𝑦_𝑚(𝑥)_±^(𝑟). By subtracting (94) from above equation and using Lemmas 38 and 50, Theorem 40, and Definition 2, one can obtain 𝐻^(𝑟)_𝑚±(𝑥)_𝔼≤sup_[]𝑟∈0,1max𝐷^𝛼Φ(𝑥)^𝑟_±−𝐷^𝛼Φ(𝑥)^𝑟_±∞+𝑎_1±^(𝑟)sup_[]𝑟∈0,1max𝑦^𝑟_±(𝑥)−𝑦^𝑟_𝑚+1,±(𝑥)_∞. (125) If we substitute the error bound from Lemmas 38 and 50 and the proof of Theorem 40, then we have 𝐻^(𝑟)_𝑚±(𝑥)_𝐸≤|Γ(𝑚+1)||Γ(𝑀1−𝛼)|Γ(𝑚𝛼+𝛼^{(𝑟)𝛼,±}+1)1(2())𝑥𝑚+1𝛼+1^0−𝛼+𝑎_1±^(𝑟)𝑀^(𝑟)_𝛼,±Γ(𝑚𝛼+𝛼+1)1(2()),𝑚+1𝛼+1 (126) or 𝐻^(𝑟)_𝑚±(𝑥)_𝔼≤𝐶1+𝑥^𝛼₀, (127) at which 𝐶=(|Γ(𝑚+1)|/|Γ(1−𝛼)|)(𝑀_𝛼^(𝑟)/Γ(𝑚𝛼+𝛼+1))1/(2(𝑚+1)𝛼+1), for 0<𝑥₀≤𝑥, 𝑥∈(0,1]. Therefore, 𝐻_𝑚^(𝑟)(𝑥) is bounded.

5.2. Nonlinear FFDEs

Consider the nonlinear FFDE ^𝑐𝐷^𝛼₀⁺𝑦(𝑥)=𝐹𝑥,𝑦(𝑥),𝑦(𝑥)𝑦(0)=𝑦,0<𝛼<1,₀,𝑦₀, (128) where 𝐹 can be nonlinear in general, 𝑦(𝑥)∶𝐿^𝔼∩𝐶^𝔼 is a continuous fuzzy-valued function, and ^𝑐𝐷^𝛼₀⁺ indicates the fuzzy Caputo fractional derivative of order 𝛼.

In order to use LOM for this problem, we first approximate 𝑦(𝑥) and (^𝑐𝐷^𝛼₀⁺𝑦)(𝑥) as (100) and (103), respectively. By replacing these equations in (128), we have 𝐶^𝑇⊙𝐷^(𝛼)Φ(𝑥)≃𝐹𝑥,𝐶^𝑇Φ(𝑥),𝐶^𝑇Φ(𝑥). (129) We intend to find the fuzzy coefficients 𝐶_𝑚. Also by substituting (100) in the initial condition of nonlinear FFDE (128), we have 𝑦(0)=𝐶^𝑇Φ(0)=𝑦₀,𝑦(0)=𝐶^𝑇Φ(0)=𝑦₀. (130)

To find the approximate fuzzy solution 𝑦_𝑚(𝑥), we first collocate (129) at (𝑚) points. For appropriate collocation points we use the first (𝑚) shifted Legendre polynomials roots. These equations together with (130) generate (𝑚+1) nonlinear fuzzy equations which can be solved using Newton’s iterative method presented in [77].

Remark 51. The solvability of system (129) and (130) is a complicated problem and we cannot prove the existence and uniqueness of such a fuzzy solution. But in our accomplishment, we have solved this system, regarding the method in [77], using MATLAB functions. In the assumed example, these functions have prospered to gain an accurate fuzzy approximate solution of the system, even starting with a zero initial guess.

6. Numerical Examples

In this section, to demonstrate the effectiveness of the proposed method in the present paper, two different test examples are carried out. Also, the obtained numerical solutions be compared with exact solutions.

Example 52. Let us consider the following FFDE: ^𝑐𝐷^𝛼₀⁺𝑦(𝑥)(𝑥)𝑦()=𝑦=𝜆⊙𝑦,0<𝛼≤1,0;𝑟₀(𝑟),𝑦₀(𝑟),0<𝑟≤1, (131)

in which 𝑦(𝑥)∶𝐿^𝔼∩𝐶^𝔼 is a continuous fuzzy-valued function and ^𝑐𝐷^𝛼₀⁺ denotes the fuzzy Caputo fractional derivative of order 𝛼. We solve this example according to two following cases for 𝜆=1,−1.

Case 1. Suppose that 𝜆=1, and then using ^𝑐[1−𝛼]-differentiability and Theorem 33, we have the following: ^𝑐𝐷^𝛼₀⁺𝑦(𝑥;𝑟)=𝑦(𝑥)𝑦,0<𝛼≤1,(0;𝑟)=𝑦₀(𝑟),0<𝑟≤1,_𝑐𝐷^𝛼₀⁺𝑦(𝑥;𝑟)=𝑦(𝑥),0<𝛼≤1,𝑦(0;𝑟)=𝑦₀(𝑟),0<𝑟≤1, (132) in which 𝑦(0;𝑟)=[0.5+0.5𝑟,1.5−0.5𝑟].

The exact solution of FFDE is as follows: 𝑌(𝑥;𝑟)=(0.5+0.5𝑟)𝐸_𝛼,1𝑥^𝛼,0<𝛼≤1,𝑌(𝑥;𝑟)=(1.5−0.5𝑟)𝐸_𝛼,1𝑥^𝛼,0<𝑟≤1, (133) where 𝐸_𝛼,𝛼 is the classical Mittag-Leffler function 𝐸_𝛼,𝛼(𝑧)=^∞_𝑘=0𝑧^𝑘Γ(𝛼(.𝑘+1)) (134) Note that for 𝛼=1 this representation is still valid since Γ(1)=1 and 𝐸_1,1(𝑧)=𝑒^𝑧.

Remark 53. If 𝛼=1, then the fuzzy fractional differential equations (132) are converted to fuzzy integer-order differential equations. So the exact solution of problem under (1)-differentiability using Theorem 15 is as follows: 𝑌(𝑥;𝑟)=(0.5+0.5𝑟)𝑒^𝑥,𝑌(𝑥;𝑟)=(1.5−0.5𝑟)𝑒^𝑥,0<𝑟≤1. (135)

By utilizing the technique explained in Section 5, (131) with 𝜆=1 becomes: ^𝑚_𝑗=0𝑐^𝑟_𝑗𝑑^(𝛼)_𝑖,𝑗−𝐼𝐿_𝑗(𝑥;𝑟)=0𝑖=0,…,𝑚,_𝑚𝑗=0𝑐^𝑟_𝑗𝑑^(𝛼)_𝑖,𝑗−𝐼𝐿_𝑗(𝑥;𝑟)=0𝑖=0,…,𝑚, (136) or we can rewrite it in the matrix form as 𝐶^𝑇_𝑚,−𝐷^(𝛼)Φ(𝑥)𝐶−𝐼=0,_𝑇𝑚,+𝐷^(𝛼)−𝐼Φ(𝑥)=0, (137) in which (𝑐^𝑟_𝑗)=[𝑐^𝑟_𝑗,−,𝑐^𝑟_𝑗,+] and for 𝑗=0,1,…,𝑚. As it was described in Section 5, we produce 𝑚 fuzzy algebraic equations multiplied above fuzzy residual system by 𝐿_𝑖(𝑥) for 𝑖=0,1,…,𝑚−1 using orthogonal property, and so we obtain 𝐶^𝑇_𝑚,−𝐷^(𝛼)𝐶−𝐼=0,_𝑇𝑚,+𝐷^(𝛼)−𝐼=0. (138) Now, using the initial condition (131), we have 𝑦(0;𝑟)≃𝐶^𝑇_𝑚,−Φ_𝑚=(0.5+0.5𝑟),𝑦(0;𝑟)≃𝐶^𝑇_𝑚,+Φ_𝑚=(1.5−0.5𝑟). (139) The above three equations generate a set of 𝑚+1 fuzzy linear algebraic equations. As a result, for 𝛼=0.85 and 𝑚=3, one can gain 𝐷^(0.85)=00001.86390.3901−0.17550.1095−0.39014.52670.8696−0.40781.6885−0.47946.78311.3797 (140) the unknown coefficients 𝑐_𝑗 can be gained by substituting matrix 𝐷^(0.85) in the previous mentioned systems.

Case 2. Suppose that 𝜆=−1, then using ^𝑐[2−𝛼]-differentiability and Theorem 33, the FDEs are obtained as same as (132). Also the exact solution is given by 𝑌(𝑥;𝑟)=(0.5+0.5𝑟)𝐸_𝛼,1−𝑥^𝛼,0<𝛼≤1,𝑌(𝑥;𝑟)=(1.5−0.5𝑟)𝐸_𝛼,1−𝑥^𝛼,0<𝑟≤1. (141) By applying the technique described in Sections 4 and 5, namely 𝑚=2, we may write the approximate solution for Cases 1 and 2 in the forms 𝑦₂(𝑥)=𝑐₀𝐿₀(𝑥)+𝑐₁𝐿₁(𝑥)+𝑐₂𝐿₂(𝑥),𝑦₂(𝑥)=𝑐₀𝐿₀(𝑥)+𝑐₁𝐿₁(𝑥)+𝑐₂𝐿₂(𝑥). (142) So regarding the (100)–(110) and (142), the unknown parameters, 𝑐_𝑗,𝑐_𝑗,𝑗=0,1,2, are achieved for both cases.

Remark 54. If 𝛼=1, then the fuzzy fractional differential equation (131) with 𝜆=−1 is converted to fuzzy integer-order differential equation. So the exact solution of problem under (2)-differentiability using Corollary 16 is as follows: 𝑌(𝑥;𝑟)=(0.5+0.5𝑟)𝑒^−𝑥,𝑌(𝑥;𝑟)=(1.5−0.5𝑟)𝑒^−𝑥,0<𝑟≤1. (143)

The approximate solution, exact solution and absolute errors are depicted for Case 1 of Example 52 in Table 1 for 𝑥=1 with 𝛼=.85. It can be seen that a few terms of the shifted Legendre functions are required to achieve a suitable approximation which demonstrate the applicability of the proposed method for this problem. Also the absolute errors for 𝑚=2,4,8 with 𝛼=0.85 are plotted in Figure 2 which show the decreasing of the error with the increasing of the number of Legendre functions. The fuzzy approximate solution is shown in Figure 3 for 𝛼=0.85,𝑚=4.

Table 1 

The result of the proposed method for Case 1 of Example 52 with 𝛼=0.85 and 𝑚=4.

Figure 2 

The absolute error for different 𝑚 of Example 52, Case 1. 𝛼=0.85.

Figure 3 

The Fuzzy approximate solution of Example 52, Case 1. 𝛼=0.85,𝑚=4.

Remark 55. Figure 2 depicts the absolute errors for 𝑦(𝑥;𝑟) of Example 52, Case 1. In the same way, if we consider 𝑦(𝑥;𝑟), then analogously to the demonstration of Figure 2, we can obtain the absolute errors.

In the Table 2, the fuzzy approximate solution for Case 2 of Example 52 is considered. The fractional Caputo derivative is 𝛼=0.75 and the number of Legendre functions are 𝑚=9. The result are computed for 𝑥=1 with different 𝑟-cut. Again we can see that Table 2 demonstrates the validity of the method for this kind of problems. Furthermore, the absolute errors of Example 52, Case 2, are plotted in Figure 4 with 𝑚=5,9,11 and the fuzzy approximate solution is shown in Figure 5 for 𝛼=0.75,𝑚=9.

Table 2 

The result of the proposed method for Case 2 of Example 52 with 𝛼=0.75 and 𝑚=9.

Figure 4 

The Absolute Error for different 𝑚 of Example 52, Case 2. 𝛼=0.75.

Figure 5 

The fuzzy approximate solution of Example 52, Case 2. 𝛼=0.75,𝑚=9.

Remark 56. Figure 4 shows the absolute errors for 𝑦(𝑥;𝑟) of Example 52, Case 2. In the same way, if we consider 𝑦(𝑥;𝑟), then analogously to the demonstration of Figure 4, we can obtain the absolute errors.

Example 57. Let us consider the following FFDE: ^𝑐𝐷^𝛼₀⁺𝑦(𝑥)(𝑥)⊕()𝑦()=𝑦=𝜆⊙𝑦𝑥+1,0<𝛼≤1,0;𝑟₀(𝑟),𝑦₀(𝑟)0<𝑟≤1. (144)

Here, suppose that 𝜆=−1, then using ^𝑐[2−𝛼]-differentiability and Theorem 33 we have the following: ^𝑐𝐷^𝛼₀⁺𝑦(𝑥;𝑟)=−𝑦(𝑥)𝑦+𝑥+1,0<𝛼≤1,(0;𝑟)=𝑦₀(𝑟),0<𝑟≤1,_𝑐𝐷^𝛼₀⁺𝑦(𝑥;𝑟)=−𝑦(𝑥)+𝑥+1,0<𝛼≤1,𝑦(0;𝑟)=𝑦₀(𝑟),0<𝑟≤1, (145) where 𝑦(0;𝑟)=[0.5+0.5𝑟,1.5−0.5𝑟]. Also, regarding (145), the exact solution of FFDE (144) be determined as 𝑌(𝑥;𝑟)=(0.5+0.5𝑟)𝐸_𝛼,1𝑥^𝛼+^𝑥₀(𝑥−𝑡)^𝛼−1𝐸_𝛼,𝛼−(𝑥−𝑡)^𝛼(𝑥+1)𝑑𝑡,0<𝛼≤1,𝑌(𝑥;𝑟)=(1.5−0.5𝑟)𝐸_𝛼,1𝑥^𝛼+^𝑥₀(𝑥−𝑡)^𝛼−1𝐸_𝛼,𝛼−(𝑥−𝑡)^𝛼(𝑥+1)𝑑𝑡,0<𝑟≤1. (146) Applying the shifted Legendre method with LOM technique explained in Sections 4 and 5, we have ^𝑚_𝑗=0𝑐^𝑟_𝑗𝑑^(𝛼)_𝑖,𝑗+𝐼𝐿_𝑗(𝑥;𝑟)=^𝑚_𝑗=0𝑓^𝑟_𝑗𝐿_𝑗(𝑥;𝑟)𝑖=0,…,𝑚,_𝑚𝑗=0𝑐^𝑟_𝑗𝑑^(𝛼)_𝑖,𝑗+𝐼𝐿_𝑗(𝑥;𝑟)=^𝑚_𝑗=0𝑓^𝑟_𝑗𝐿_𝑗(𝑥;𝑟)𝑖=0,…,𝑚, (147) or it can be written in the matrix form as 𝐶^𝑇_𝑚,−𝐷^(𝛼)+𝐼Φ(𝑥)=𝐹^𝑇_𝑚,−Φ(𝑥),𝐶^𝑇_𝑚,+𝐷^(𝛼)+𝐼Φ(𝑥)=𝐹^𝑇_𝑚,+Φ(𝑥), (148) in which (𝑐^𝑟_𝑗)=[𝑐^𝑟_𝑗,−,𝑐^𝑟_𝑗,+] and (𝑓^𝑟_𝑗)=[𝑓^𝑟_𝑗,−,𝑓^𝑟_𝑗,+] for 𝑗=0,1,…,𝑚. As was illustrated in Section 5, we make 𝑚 fuzzy algebraic equations by the result of the inner product of fuzzy residual with 𝐿_𝑖(𝑥),  𝑖=0,1,…,𝑚−1, so we gain 𝐶^𝑇_𝑚,−𝐷^(𝛼)+𝐼=𝐹^𝑇_𝑚,−,𝐶^𝑇_𝑚,+𝐷^(𝛼)+𝐼=𝐹^𝑇_𝑚,+. (149) Also for the initial condition (144), we have 𝑦(0;𝑟)≃𝐶^𝑇_𝑚,−Φ_𝑚=(0.5+0.5𝑟),𝑦(0;𝑟)≃𝐶^𝑇_𝑚,+Φ_𝑚=(1.5−0.5𝑟). (150) The above equations produce a set of 𝑚+1 fuzzy linear algebraic equations. As a result, namely, with 𝑚=2, the approximate fuzzy solution should be written in the form 𝑦(𝑥)=²_𝑗=0𝑐_𝑗𝐿_𝑗(𝑥)=𝑐₀𝑐₁𝑐₂^𝑇Φ(𝑥),𝑦(𝑥)=²_𝑗=0𝑐_𝑗𝐿_𝑗(𝑥)=𝑐₀𝑐₁𝑐₂^𝑇Φ(𝑥). (151) Finally, the corresponding fuzzy approximate solution from (151) can be acquired.

For case 𝛼=0.95 and 𝑚=6, 𝑥=1 with 𝑟=0.1, we obtain the fuzzy approximate solution in a series expansion as 𝑦(𝑥,0.1)=0.5500+0.5533𝑥−0.1915𝑥²+1.1965𝑥³−1.9554𝑥⁴+1.5046𝑥⁵−0.4478𝑥⁶,𝑦(𝑥,0.1)=1.4500−0.5496𝑥+1.3466𝑥²−1.9084𝑥³+2.4940𝑥⁴−1.8195𝑥⁵−1.8195𝑥⁶. (152)

Remark 58. If 𝛼=1, then the fuzzy fractional differential equation (144) with 𝜆=−1 is converted to fuzzy integer-order differential equation. Therefore, the exact solution of problem under (2)-differentiability using Corollary 16 is as follows: 𝑌(𝑥;𝑟)=𝑥+(0.5+0.5𝑟)𝑒^−𝑥,𝑌(𝑥;𝑟)=𝑥+(1.5−0.5𝑟)𝑒^−𝑥. (153)

In order to evaluate the advantages and the accuracy of using the presented method for the fuzzy fractional differential equations, Example 57 is considered. The results are illustrated in Table 3 are the approximate solutions are compared with exact solutions and the absolute errors are derived for 𝛼=0.95 with 𝑚=6 at 𝑥=1. It is obvious that the fuzzy approximate solutions are in high agreement with the fuzzy exact solutions. Moreover, different numbers of Legendre polynomials are applied to obtain the absolute errors for this problem which can be seen in Figure 6. This graph shows that the method has a good convergence rate. Also, the fuzzy approximate solution is shown for 𝛼=0.95 in Figure 7. The absolute error of different order of fractional differentiability can be considered for 𝑥=1 in Figure 8 and the approximate solutions are shown in Figure 9 that shows that this approach can apply to solve the problem effectively with different fuzzy fractional order of derivatives with suitable errors. It can bee seen that the 𝛼 approaches an integer order, and the error has a tendency to decrease, as expected.

Table 3 

The result of the proposed method for Example 57 with 𝛼=0.95 and 𝑚=6.

Figure 6 

The absolute error for different 𝑚 of Example 57, 𝛼=0.95.

Figure 7 

The fuzzy approximate solution of Example 57, 𝛼=0.95,𝑚=6.

Figure 8 

The absolute error for different 𝛼 of Example 57, 𝑚=8.

Figure 9 

The approximate solution for 𝛼=0.75(∗),0.85(∘),and0.95(+) of Example 57, 𝑚=8.

Remark 59. Figure 6 illustrates the absolute errors for 𝑦(𝑥;𝑟) of Example 57. In the same way, if we consider for 𝑦(𝑥;𝑟), then analogously to the demonstration of Figure 6, we can obtain the absolute errors.

Example 60. For our third example, consider the inhomogeneous linear equation in [78] with fuzzy initial values, and so we have: ^𝑐𝐷^𝛼₀₊𝑦(𝑥)+𝑦(𝑥)=2𝑥^2−𝛼Γ()−𝑥3−𝛼Γ(2−^1−𝛼𝛼)+𝑥²𝑦()=[]−𝑥,0;𝑟−1+𝑟,1−𝑟,0<𝛼≤1,0≤𝑥≤1, (154) in which 𝑦(𝑥)∶𝐿^𝔼[0,1]∩𝐶^𝔼[0,1] is a continuous fuzzy-valued function and ^𝑐𝐷^𝛼₀⁺ denotes the fuzzy Caputo fractional derivative of order 𝛼.

Now, using [1−𝛼]-differentiability and Theorem 33, we have the following: ^𝑐𝐷^𝛼₀₊𝑦(𝑥;𝑟)+𝑦(𝑥;𝑟)=2𝑥^2−𝛼Γ()−𝑥3−𝛼Γ(2−^1−𝛼𝛼)+𝑥²𝑦−𝑥,(0;𝑟)=−1+𝑟,0<𝛼≤1,0≤𝑥≤1,_𝑐𝐷^𝛼₀₊𝑦(𝑥;𝑟)+𝑦(𝑥;𝑟)=2𝑥^2−𝛼Γ()−𝑥3−𝛼Γ(2−^1−𝛼𝛼)+𝑥²−𝑥,𝑦(0;𝑟)=1−𝑟,0<𝛼≤1,0≤𝑥≤1. (155) Solving (155) leads to determining the solution of FFDE (154) as follows: 𝑦(𝑥;𝑟)=(−1+𝑟)𝐸_𝛼,1−𝑥^𝛼+^𝑥₀(𝑥−𝑡)^𝛼−1𝐸_𝛼,𝛼−(𝑥−𝑡)^𝛼×2𝑥^2−𝛼Γ()−𝑥3−𝛼Γ(2−^1−𝛼𝛼)+𝑥²−𝑥𝑑𝑡,0≤𝑟≤1,𝑦(𝑥;𝑟)=(1−𝑟)𝐸_𝛼,1−𝑥^𝛼+^𝑥₀(𝑥−𝑡)^𝛼−1𝐸_𝛼,𝛼−(𝑥−𝑡)^𝛼×2𝑥^2−𝛼Γ()−𝑥3−𝛼Γ(2−^1−𝛼𝛼)+𝑥²−𝑥𝑑𝑡,0≤𝑟≤1. (156)

By exploiting the method proposed in Section 5, the equations are acquired by ^𝑚_𝑗=0𝑐^𝑟_𝑗𝑑^(𝛼)_𝑖,𝑗+𝐼𝐿_𝑗(𝑥;𝑟)=^𝑚_𝑗=0𝑓^𝑟_𝑗𝐿_𝑗(𝑥;𝑟)𝑖=0,…,𝑚,_𝑚𝑗=0𝑐^𝑟_𝑗𝑑^(𝛼)_𝑖,𝑗+𝐼𝐿_𝑗(𝑥;𝑟)=^𝑚_𝑗=0𝑓^𝑟_𝑗𝐿_𝑗(𝑥;𝑟)𝑖=0,…,𝑚, (157) or we can rewrite it in the matrix form as 𝐶^𝑇_𝑚,−𝐷^(𝛼)+𝐼Φ(𝑥)=𝐹^𝑇_𝑚,−Φ(𝑥),𝐶^𝑇_𝑚,+𝐷^(𝛼)+𝐼Φ(𝑥)=𝐹^𝑇_𝑚,+Φ(𝑥), (158) in which (𝑐^𝑟_𝑗)=[𝑐^𝑟_𝑗,−,𝑐^𝑟_𝑗,+] and (𝑓^𝑟_𝑗)=[𝑓^𝑟_𝑗,−,𝑓^𝑟_𝑗,+] for 𝑗=0,1,…,𝑚. As it was described in Section 5, we create 𝑚 fuzzy algebraic equations multiplied in above system by 𝐿_𝑖(𝑥) for 𝑖=0,1,…,𝑚−1 and implemented in the inner product using orthogonal property, and so we gain: 𝐶^𝑇_𝑚,−𝐷^(𝛼)+𝐼=𝐹^𝑇_𝑚,−,𝐶^𝑇_𝑚,+𝐷^(𝛼)+𝐼=𝐹^𝑇_𝑚,+. (159) Also for the initial condition (154), we have 𝑦(0;𝑟)≃𝐶^𝑇_𝑚,−Φ_𝑚=−1+𝑟,𝑦(0;𝑟)≃𝐶^𝑇_𝑚,+Φ_𝑚=1−𝑟. (160) Equations (159) and (160) produce a set of 𝑚+1 fuzzy linear algebraic equations. As a result, for 𝛼=0.75𝑟=0.1 and 𝑚=3 the unknown coefficients 𝑐_𝑗 can be achieved as 𝐶^𝑇_3,−=[−0.6914,0.3147,0.1144,0.0083] and 𝐶^𝑇_3,+=[0.3725,−0.1418,0.2954,−0.0903]. Also with these assumptions, one has 𝐷^(0.75)=00001.76520.5884−0.22630.1305−0.58843.66621.2114−0.49791.5389−0.62305.10811.8265. (161)

The approximate and exact solution of FFDEs (154) is illustrated in Table 4 at 𝑥=1 with the absolute errors. In this table, three shifted Legendre polynomials are considered to derive the approximate solution. It is clear that with the lower number of shifted Legendre functions, using proposed method, one can reach a resealable approximate solution. Also Figure 10 depicts the absolute error with 𝑚=5 and various values of 𝛼. This figure exhibits the applicability and validity of the proposed technique for this example. Furthermore, the fuzzy approximate solution for 𝛼=0.75 with 𝑚=5 is shown in Figure 11 at 𝑥=1. Additionally, in Figure 12, the fuzzy approximate solution is compared with different values of 𝛼. It can be seen that as 𝛼 tends to 1, the solution of the fuzzy fractional differential equations tends to that of the fuzzy integer-order differential equations. Finally, Figure 13 displays the absolute error for 𝛼=0.75 with different values of 𝑚 which is obvious that with increasing the number of shifted Legendre functions, the absolute error of the problem has been decreased gently.

Table 4 

The result of the proposed method for Example 60 with 𝛼=0.95 and 𝑚=3.

Figure 10 

The absolute error for different 𝛼 of Example 60, 𝑚=5.

Figure 11 

The fuzzy approximate solution of Example 60, 𝛼=0.75,𝑚=5.

Figure 12 

The fuzzy approximate solution for different 𝛼: (1: star, 0.95: square, 0.85: ∘, 0.75: ×, 0.65: pentagram) of Example 60, 𝑚=5.

Figure 13 

The absolute error for different 𝑚 of Example 60, 𝛼=0.75.

Remark 61. Figures 10 and 13 illustrate the absolute errors for 𝑦(𝑥;𝑟) of Example 60. In the same way, if we consider 𝑦(𝑥;𝑟), then analogously to the demonstration of Figures 10 and 13, we can obtain the absolute errors.

Now, we consider two applicable examples which are new under uncertainty represented by fuzzy-valued functions. Firstly, a brief history is given in the deterministic case of Example 62. For more details, see [1, 79].

There are two kinds of electrical circuits which are related to the fractional calculus. Circuits of the first types are supposed to consist of capacitors and resistors, which are described by conventional (integer-order) models. A circuit expressing fractional order behavior is called fractance. Circuits of the second type may consist of resistors, capacitors, and fractances. The term fractance was suggested initially by Mehaute and Crepy [80]. Now, in order to scrutinize the mentioned problem in a real case, we use the fuzzy initial value 𝑦(0;𝑟) and the concept of Caputo’s H-differentiability for fractional derivative of 𝑦(𝑥). So we have the following fuzzy fractional oscillation differential equation.

Example 62. Consider an electrical circuit (LR circuit) with an AC source. The current equation of this circuit can be written as follows: ^𝑐𝐷^𝛼₀⁺𝑦(𝑥)𝑅=−𝐿𝑦(𝑥)(𝑥)𝑦()=[],+𝑣,0≤𝑥≤1,0;𝑟0.96+0.04𝑟,1.01+0.01𝑟 (162) in which 𝑅 is the circuit resistance, and 𝐿 is a coefficient, corresponding to the solenoid and 0≤𝑟≤1. Suppose that 𝑣(𝑥)=sin(𝑥), 𝑅=1ohm and 𝐿=1H, so (162) can be rewritten as ^𝑐𝐷^𝛼₀⁺𝑦(𝑥)(𝑥)(𝑥)𝑦()=[].=−𝑦+sin,0≤𝑥≤1,0;𝑟0.96+0.04𝑟,1.01+0.01𝑟 (163)

Now, using [2−𝛼]-differentiability and Theorem 33, we have the following: ^𝑐𝐷^𝛼₀₊𝑦(𝑥;𝑟)+𝑦()(𝑥),𝑦𝑥;𝑟=sin(0;𝑟)=(0.96+0.04𝑟),0<𝛼≤1,0≤𝑥≤1,_𝑐𝐷^𝛼₀₊𝑦(𝑥;𝑟)+𝑦(𝑥;𝑟)=sin(𝑥),𝑦(0;𝑟)=(1.01+0.01𝑟),0<𝛼≤1,0≤𝑥≤1. (164) The exact solution of (163) under (2)-differentiability for integer order is given by 𝑌₁()=1𝑡;𝑟2(sin(𝑡)−cos(𝑡))+𝑒^−𝑡(1.46+0.04r),𝑌₁()=1𝑡;𝑟2(sin(𝑡)−cos(𝑡))+𝑒^−𝑡(1.51−0.01r). (165)

By applying the technique described in Section 5, we approximate fuzzy solution as ^𝑐𝐷^𝛼𝑦(𝑥)≃𝐶^𝑇⊙𝐷^(𝛼)Φ_𝑚(𝑥),𝑦(𝑥)≃𝑦_𝑚(𝑥)=𝐶^𝑇_𝑚⊙Φ_𝑚,sin(𝑥)≃𝑓_𝑚(𝑥)=𝐹^𝑇_𝑚⊙Φ_𝑚, (166) where vector 𝐹^𝑇_𝑚 is obtained as (102). With 𝑚=3, 𝛼=0.75, and 𝑟=0.1, we have 𝐷^(0.75)=,𝐹00001.76520.5884−0.22630.1305−0.58843.66621.2114−0.49791.5389−0.62305.10811.8265_𝑇4=0.45970.4279−0.0392−0.0072, (167) and subsequently, we generate 𝑚+1 fuzzy linear equations using (110) and also from substituting (100) in the initial conditions of (164), we can obtain the unknown fuzzy coefficients {𝑐_𝑗}³_𝑗=0 as 𝑐₀=0.7222,𝑐₁𝑐=−0.0497,₂=0.1272,𝑐₃=−0.0650,𝑐₀=0.7488,𝑐₁=−0.0611,𝑐₂=0.1317,𝑐₃=−0.0674, (168) and finally, the fuzzy approximate solution is given by 𝑦₄(𝑥,0.1)=0.9640−1.6420𝑥+2.7118𝑥²−1.2990𝑥³,𝑦₄(𝑥,0.1)=1.0090−1.7215𝑥+2.8128𝑥²−1.3483𝑥³. (169) In Table 5, the numerical results of (162) are depicted at 𝑥=1 for 𝛼=0.85 and 𝑚=8. The absolute error confirms that the proposed method approximates the fuzzy solution with a suitable accuracy. Figure 14 shows the absolute error of the proposed method with a different value of 𝛼 with 𝑚=10 which is clear as 𝛼 tends to 1, and the error is decreasing dramatically, also with the increasing of the number of Legendre functions, the absolute error is dereasing that is, shown in Figure 16. Finally the fuzzy approximate solution of (162) is demonstrated in Figure 15 in the interval 0≤𝑥≤1.

Table 5 

The result of the proposed method for Example 62 with 𝛼=0.85 and 𝑚=8.

Figure 14 

The absolute error for different 𝛼 of Example 62, 𝑚=10.

Figure 15 

The fuzzy approximate solution of Example 62, 𝛼=0.75,𝑚=10.

Figure 16 

The absolute error for different 𝑚 of Example 62, 𝛼=0.85.

Now we consider the second application of FFDES in the real world. But before we consider the problem of fuzzy model, We glimpse the vision of a non-fuzzy case.

A pharmacodynamic model is usually separated in two parts: a link model and a transduction model. The link model depicts the distribution of drug from an observed compartment (e.g., plasma) into a biophase (the effect compartment). Drug in the biophase induces a pharmacodynamic response (the effect), which is represented by a transduction model. For more details, see [81, 82].

The model describing the transfer of drug to the effect compartment, from the mechanistic interpretation, it can be considered as the pharmacokinetics of the drug. However, practically, and because drug condensation in the effect compartment does not appears, it is not often clear if the link model is actually linked with a site of action, or if it exhibits some pharmacodynamic related delay, or a combination of the two. It is begun by representing drug concentration in the effect compartment by the (Caputo) fractional differential equation as follows: ^𝛼𝐷^𝛼₀₊𝐶𝑒(𝑡)=𝑘_eo−𝐶𝑒(𝑡)+𝑓(𝑡), (170) in which 𝑓(𝑡) describes the drug input into the central compartment and 𝑘_eo is the elimination rate constant from the effect compartment. Also the action of drug is described by a nonlinear memory-less function of 𝐶𝑒.

Now, we consider the action of drug delivery by a fuzzy function under the fuzzy Caputo’s H-differentiability which is modelled by fuzzy fractional differential equation.

Example 63. Consider the following fuzzy pharmacodynamic model ^𝛼𝐷^𝛼₀₊𝐶𝑒(𝑡)=𝑘_eo(𝑡)(𝑡()=(),−𝐶𝑒+𝑓),0<𝛼<1,0≤𝑡≤1,𝐶𝑒0;𝑟0.5+0.5𝑟,1.5−0.5𝑟 (171) in which 𝑓(𝑡)=(Dose𝑘_𝑎/𝑉(𝑘_el−𝑘_𝑎))(𝑒^{−𝑘el𝑡}−𝑒^{−𝑘𝑎𝑡}) represents drug concentration in the central compartment following, for example, an oral administration of a dose, in this example dose=1 and parameters values are 𝑉=1,𝑘_𝑎=5.0,𝑘_el=0.5, and 𝑘_eo=1.

The exact solution of (171) with 𝛼=1 under generalized differentiability is given by 𝑌(𝑡;𝑟)=𝑒^−𝑡−209𝑒^𝑡/2−5𝑒18^−4𝑡+3+0.5𝑟,𝑌(𝑡;𝑟)=𝑒^−𝑡−209𝑒^𝑡/2−5𝑒18^−4𝑡+4−0.5𝑟. (172) We apply the method presented in Section 5 and solve this problem for 𝑚=3,  𝛼=0.95, and 𝑟=0.2. Hence, we have 𝐷^(0.95)=,𝐹00001.95660.1432−0.07430.0501−0.14325.48491.3447−0.18231.8823−0.20158.82580.5752_𝑇4=−0.6537−0.18860.2729−0.1322, (173) and thereafter the fuzzy unknown coefficients for the fuzzy approximation of (171) can be acquired easily by replacing the above results in (110) and solving this algebraic fuzzy linear equations system. So we have 𝑐₀=0.1744,𝑐₁𝑐=−0.4368,₂=0.0573,𝑐₃=0.0185,𝑐₀=0.6417,𝑐₁=−0.6635,𝑐₂=0.1043,𝑐₃=0.0095, (174) and ultimately the fuzzy approximate solution is given by 𝑦₄(𝑥,0.1)=0.6−0.9392𝑥−2.2482𝑥²+0.3823𝑥³,𝑦₄(𝑥,0.1)=1.4−1.8391𝑥+0.3419𝑥²+0.1891𝑥³. (175) In Table 6, the approximate fuzzy solutions (171) are compared with the exact solution at 𝑥=1 for 𝛼=0.95 and 𝑚=12. The absolute error shows that the proposed method approximates the fuzzy solution with a high accuracy. Moreover, Figure 17 describes the absolute error of the proposed method with a different value of 𝛼 with 𝑚=8. It is obvious that the method reaches a good approximation as 𝛼 approaches 1; also again with the increasing of the number of Legendre functions, the absolute error is decreasing that is, shown in Figure 19. The fuzzy approximate solution of the problem is displayed in Figure 18 in the interval 0≤𝑥≤1.

Table 6 

The result of the proposed method for Example 63 with 𝛼=0.95 and 𝑚=12.

Figure 17 

The absolute error for different 𝛼 of Example 63, 𝑚=8.

Figure 18 

The fuzzy approximate solution of Example 63, 𝛼=0.95,𝑚=8.

Figure 19 

The absolute error for different 𝑚 of Example 63, 𝛼=0.95.

Example 64. Consider the following initial value problem of nonlinear FFDE: ^𝑐𝐷^𝛼₀⁺𝑦(𝑥)=3𝐴⊙𝑦²𝑦()=√,0<𝛼<1,0≤𝑥≤1,0;𝑟0.5𝑟,0.2√1−𝑟+0.5, (176) in which 𝑦(𝑥)∶𝐿^𝔼[0,1]∩𝐶^𝔼[0,1] is a continuous fuzzy-valued function and ^𝑐𝐷^𝛼₀⁺ denotes the fuzzy Caputo fractional derivative of order 𝛼. Also 𝐴=[1+𝑟,3−𝑟] is a constant fuzzy number.

For approximating the fuzzy solution of (162), using the method described in Sections 4 and 5, if we approximate the solution by 𝑚 shifted Legendre functions, then it needs to consider the first 𝑚 roots of the shifted Legendre polynomial 𝐿_𝑚+1(𝑥). Initially using (162), we have ^𝑚_𝑗=0𝑐^𝑟_𝑗𝑑^(𝛼)_𝑖,𝑗𝐿_𝑗(𝑥;𝑟)−3(3−𝑟)^𝑚_𝑗=0𝑐^𝑟_𝑗𝐿_𝑗(𝑥;𝑟)²=0𝑖=0,…,𝑚,_𝑚𝑗=0𝑐^𝑟_𝑗𝑑^(𝛼)_𝑖,𝑗𝐿_𝑗(𝑥;𝑟)−3(1+𝑟)^𝑚_𝑗=0𝑐^𝑟_𝑗𝐿_𝑗(𝑥;𝑟)²=0𝑖=0,…,𝑚, (177) or we can rewrite it in the matrix form as 𝐶^𝑇_𝑚,−𝐷^(𝛼)Φ(𝑥)−3(3−𝑟)𝐶^𝑇_𝑚,−Φ(𝑥)²𝐶=0,_𝑇𝑚,+𝐷^(𝛼)Φ(𝑥)−3(1+𝑟)𝐶^𝑇_𝑚,+Φ(𝑥)²=0, (178) in which (𝑐^𝑟_𝑗)=[𝑐^𝑟_𝑗,−,𝑐^𝑟_𝑗,+] and for 𝑗=0,1,…,𝑚. As it was described in Section 5, we produce 𝑚 fuzzy algebraic equations by collecting (162) in the first 𝑚 root of 𝐿_𝑚+1(𝑥). Also using the initial condition (162), we have 𝑦(0;𝑟)≃𝐶^𝑇_𝑚,−Φ_𝑚(0)=0.5√𝑟,𝑦(0;𝑟)≃𝐶^𝑇_𝑚,+Φ_𝑚(0)=0.2√1−𝑟+0.5. (179) These roots can be put in the (178) and deriving three nonlinear fuzzy equation, then these equations be coupled with (179). Finally with solving the fuzzy nonlinear equations system, the unknown coefficients 𝑐_𝑗̃ are obtained.

Now, as an explanation, assume that 𝑚=3, then the first three roots of 𝐿₄(𝑥) are as follows: 𝑥₀=0.0694,𝑥₁=0.9306,𝑥₂=0.3300, (180) which for 𝛼=0.95 and 𝑟=0.1 are 𝑐₀=0.2593,𝑐₁𝑐=0.1106,₂=−0.0032,𝑐₃=0.0054,𝑐₀=0.1171,𝑐₁=−0.0203,𝑐₂=0.0045,𝑐₃=0.0155. (181) Finally, using the above results and (100), the approximate solution is computed in each particular point easily.

The approximate solution for Example 64 with different fractional Caputo order is derived in Table 7. The solution has been estimated using 3 shifted Legendre functions. The results in this table are at 𝑥=1. It is clear that the method is applicable and valid for different fuzzy fractional order 𝛼. Furthermore, by changing the number of shifted Legendre functions, the approximate solution for 𝛼=0.75 at 𝑥=1 is revealed in Figure 20 which shows that with only a few number of shifted Legendre polynomials, desired results are available. Ultimately, the fuzzy approximate solution is described in Figure 21 with 𝛼=0.75 and 𝑚=3.

Table 7 

The approximate solution of the proposed method for Example 64 at 𝑥=1 with 𝑚=3.

Figure 20 

The approximate solution for different 𝑚 of Example 64, 𝛼=0.75.

Figure 21 

The fuzzy approximate solution of Example 64, 𝛼=0.75,𝑚=3.

Remark 65. It is consequential to note that, for a crisp differential equation of integer or fractional order, the problem ^𝑐𝐷^𝛼₀⁺𝑦(𝑥)+𝐹𝑥,𝑦(𝑥)=0,𝑥∈[0,𝑋], (182) is the same as ^𝑐𝐷^𝛼₀⁺𝑦(𝑥)=−𝐹𝑥,𝑦(𝑥),𝑥∈[0,𝑋], (183) in which is 𝐹 is nonlinear operator, and in the case of space 𝔼, both problems are not equivalent [83].

7. Conclusion

In this paper, a cluster of orthogonal functions named shifted Legendre function is used to solve fuzzy fractional differential equations under Caputo type. The benefit of the shifted Legendre operational matrices method over other existing orthogonal polynomials is its simplicity of execution besides some other advantages. A general formulation of the proposed method is provided in details. Also a complete error analysis is considered in Section 5. In this paper, we tried to answer some questions:(i)how to derive a shifted Legendre operational matrice method over the fuzzy fractional integration and differentiation(ii)how to solve the fuzzy fractional order differential equations via the shifted Legendre operational matrices of the fuzzy fractional derivative.Moreover, it should be mentioned that another advantage of this technique is that it decreases these problems to the degree of solving a system of fuzzy algebraic equations thus extremely making the problems easier. Several examples was carried out to depict the effectiveness and the absence of complexity of the proposed method. The achieved solutions have a satisfactory results obtained with only a small number of Legendre polynomials.

For future research, we will consider this method for solving FFDEs with order 1<𝛼<2. Also we will apply it under Riemann-Liouville’s H-differentiability. Apart from this, the other orthogonal functions like Jacobi polynomials will be extended for solving FFDEs.

Acknowledgments

The authors would like to thank Professor Dumitru Baleanu, the Leader Guest Editor, for his consideration about this paper and the anonymous reviewers for their helpful suggestions and comments.