Numerical Solution of System of Fractional Differential Equations in Imprecise Environment

Abstract

Fractional calculus and fuzzy calculus theory, mutually, are highly applicable for showing different aspects of dynamics appearing in science. This chapter provides comprehensive discussion of system of fractional differential models in imprecise environment. In addition, presenting a new vast area to investigate numerical solutions of fuzzy fractional differential equations, numerical results of proposed system are carried out by the Grünwald‐Letnikov's fractional derivative. The stability along with truncation error of the Grünwald‐Letnikov’s fractional approach is also proved. Moreover, some numerical experiments are performed and effective remarks are concluded on the basis of efficient convergence of the approximated results towards the exact solutions and on the depictions of error bar plots.

1. Introduction

It is worthwhile mentioning, since last few decades, the theory of fractional calculus has gained significant importance in almost every branch of science, for having the capability to consider integrals and derivatives of any arbitrary order. The characteristic feature of generalizing the classic integer‐order differentiation and n‐fold integration to arbitrary fractional order have broadened its application in modeling several phenomena of physics, mathematics, and engineering. The differential models of fractional order, due to the nonlocal properties of fractional operator, are excellent instruments for providing information about the current as well as the historical state of the system. For these reasons, it is intensively developed and advanced, and existence of its solution is studied by well‐known authors, Euler, Laplace, Liouville, Riemann, Fourier, Abel, Caputo, etc., to further widen its scope in describing various real‐world problems of science, for instance see [1–6]. Another wide‐spreading exploration of mathematics is theory of fuzzy calculus, which has a lot of interesting applications in physics, engineering, mechanics, and many others. It is the theory of a particular type of interval‐valued functions, in which mapping is made in such a way that it takes all the possible values in $\left[0,1\right]$ and not only the crisp values as found in usual interval‐valued functions. After the inception of fuzzy set theory by Zadeh [7], its attributes have been extended and established to overcome impreciseness of parameters and structures in mathematical modeling, reasoning, and computing [8–12].

Advanced development of mathematical theories and techniques has gained very high standard. On the basis of classical theories, new theories are pioneered by undergoing its inadequacies and widening its scope in many disciplines. In a similar manner, the aforementioned theories have been brought together in modeling different aspects of applied sciences, to analyze the change in the respective system at each fractional step with the uncertain parameters. Agarwal et al. [13] initiatively incorporated uncertainty into dynamical system, modeled fractional differential equations with uncertainty, and studied its possible solutions. Ahmad et al. [14] described the situation of impreciseness of initial values of fractional differential equations and discussed its solutions by utilizing Zadeh’s extension principle. In [15, 16], authors considered the concept of Caputo and Riemann fractional derivative, respectively, together with the Hukuhara differentiability and demonstrated the fuzzy fractional differential equations and a lot of others [17–23].

In light of noteworthy applications of above‐mentioned theories, in this chapter, we demonstrate fractional order dynamical models in fuzzy environment to depict unequivocal fractional differential equations of dynamical system. Moreover, we investigate its numerical solutions using the well‐known Grünwald‐Letnikov's fractional definition. This definition is widely applicable as a numerical scheme to solve linear and nonlinear differential equations of fractional order [24–26]. It is considered as an extended form of the classical Euler method. Here it will be utilized, for the first time, to solve fractional differential equations of imprecise functions. Sequentially, this chapter features description of fuzzy theory and fuzzy‐valued functions for the explanation of impreciseness, modeling of system of nonlinear fractional order differential equations with imprecise functions, deliberation of Grünwald‐Letnikov’s fractional approach in conjunction with its truncation error for the proposed system, tabulated and pictorial investigations of some examples, and conclusive remarks of the undergone experiments and findings of the whole manuscript.

2. Basic descriptions

Fuzzy calculus theory is the branch of mathematical analysis that deals with the interval analysis of imprecise functions. This section comprises some rudiments of fuzzy calculus theory and acquaints the necessary notations that are prerequisite for the whole paper. All the below‐mentioned descriptions are widely elaborated and used in literature, for instance [13–23].

2.1. Fuzzy numbers

Let

be the set of subsets of the real axis

. If

and $\tau :\left[0,1\right]\to$

such that, $\tau$ is normal, fuzzy convex, upper semi‐continuous membership function and compactly supported on the real axis

, then

is said to be the space of fuzzy numbers τ. Any

can be represented in level sets explicitly, i.e. ${\left[\tau \right]}^{ƛ}=\left[\underset{\_}{\tau }\left(ƛ\right),\overline{\tau }\left(ƛ\right)\right]$

for $ƛ\in \left[0,1\right]$

, where $\underset{\_}{\tau }\left(ƛ\right)$

and $\overline{\tau }\left(ƛ\right)$

signify as the lower and upper branches of $\tau$, respectively, that satisfy the following conditions:

1. $\underset{\_}{\tau }\left(ƛ\right)$

   
   is bounded non‐decreasing lower function, left continuous on $\left(0,1\right]$ and right continuous at $ƛ=0$
   

2. $\underset{\_}{\tau }\left(ƛ\right)$

   
   is bounded non‐increasing upper function, left continuous on $\left(0,1\right]$ and right continuous at $ƛ=0$
   

3. $\underset{\_}{\tau }\left(ƛ\right)\le \overline{\tau }\left(ƛ\right)$

The sum and scalar product of any fuzzy number is the consequence of Zadeh’s extension principal. Let $\oplus$, $•$ and $\Theta$ be the symbols of addition, multiplication and subtraction, accordingly, for fuzzy numbers, which will be greatly used throughout the paper, then, for $ƛ\in \left[0,1\right]$:

1. ${\left[\tau \oplus \upsilon \right]}^{ƛ}={\left[\tau \right]}^{ƛ}\oplus {\left[\upsilon \right]}^{ƛ}=\left[\underset{\_}{\tau }\left(ƛ\right)+\underset{\_}{\upsilon }\left(ƛ\right),\overline{\tau }\left(ƛ\right)+\overline{\upsilon }\left(ƛ\right)\right]\tau ,\upsilon \in E$

2. For $a\in$${\left[a\tau \right]}^{ƛ}=a{\left[\tau \right]}^{ƛ}=\{\begin{array}{c}\left[a\underset{\_}{\tau }\left(ƛ\right),a\overline{\tau }\left(ƛ\right)\right]\text{if}a>0\\ \left\{0\right\}\text{if}a=0\\ \left[a\overline{\tau }\left(ƛ\right),a\underset{\_}{\tau }\left(ƛ\right)\right]\text{if}a<0\end{array}$

3. ${\left[\tau •\upsilon \right]}^{ƛ}=\left[\begin{array}{l}\text{min}\left\{\underset{\_}{\tau }\left(ƛ\right)\underset{\_}{\upsilon }\left(ƛ\right),\underset{\_}{\tau }\left(ƛ\right)\overline{\upsilon }\left(ƛ\right),\overline{\tau }\left(ƛ\right)\underset{\_}{\upsilon }\left(ƛ\right),\overline{\tau }\left(ƛ\right)\overline{\upsilon }\left(ƛ\right)\right\},\\ \text{max}\left\{\underset{\_}{\tau }\left(ƛ\right)\underset{\_}{\upsilon }\left(ƛ\right),\underset{\_}{\tau }\left(ƛ\right)\overline{\upsilon }\left(ƛ\right),\overline{\tau }\left(ƛ\right)\underset{\_}{\upsilon }\left(ƛ\right),\overline{\tau }\left(ƛ\right)\overline{\upsilon }\left(ƛ\right)\right\}\end{array}\right]$

4. $\tau \Theta \upsilon =\left[\text{min}\left\{\underset{\_}{\tau }\left(ƛ\right)-\underset{\_}{\upsilon }\left(ƛ\right),\overline{\tau }\left(ƛ\right)-\overline{\upsilon }\left(\wp \right)\right\},\text{max}\left\{\underset{\_}{\tau }\left(ƛ\right)-\underset{\_}{\upsilon }\left(ƛ\right),\overline{\tau }\left(ƛ\right)-\overline{\upsilon }\left(ƛ\right)\right\}\right]$

The distance between any two fuzzy numbers $\tau$ and $\upsilon$ is given by the Hausdorff metric

as:

$D\left(\tau ,\upsilon \right)=\underset{ƛ\in \left[\text{0,1}\right]}{\text{sup}}D\left({\left[\tau \right]}^{ƛ},{\left[\upsilon \right]}^{ƛ}\right)=\underset{ƛ\in \left[\text{0,1}\right]}{\text{sup}}\text{max}\left\{\left|\underset{\_}{\tau }\left(ƛ\right)-\underset{\_}{\upsilon }\left(ƛ\right)\right|,\left|\overline{\tau }\left(ƛ\right)-\overline{\upsilon }\left(ƛ\right)\right|\right\}$E1

Thus, (

,

) defines a complete metric space with the properties of Hausdorff metric for fuzzy numbers.

2.2. Fuzzy-valued Function and its fractional derivative

Any interval‐valued function

is said to be a fuzzy‐valued function if

is defined as

$\tilde{F}:R\to E$. Its

‐level set can be represented by real‐valued functions $\underset{\_}{F}\left(t;ƛ\right)$ and

$\overline{F}\left(t;ƛ\right)$ as its lower and upper branches, accordingly, i.e.

,

. Moreover, if

and

exist as finite fuzzy numbers, then

exists. Consequently, let

be the space of continuous fuzzy‐valued functions, then

if

and

are continuous. The arithmetic for any two fuzzy‐valued functions and

can be defined as previously mentioned in Section 2.1 for fuzzy numbers. Subsequent to existence of limit and continuity of

, the fuzzy‐valued function

(t) is said to be differentiable at each $t_0\in \left[a,b\right]$, if

exists, such that

${\tilde{F}}^{\prime }\left(t_0\right)=\underset{h\to 0}{\text{Lim}}\frac{\tilde{F}\left(t_0+h\right)\Theta \tilde{F}\left(t_0\right)}{h}$E2

where $h$ is taken in a way that $\left(t_0+h\right)\in \left(a,b\right)$. For

,

is said to be differentiable at $t\in \left[a,b\right]$ if its lower function

and upper function

are differentiable at $t\in \left[a,b\right]$, i.e. for all

$\in \left[0,1\right]$,

${\tilde{F}}^{\prime }\left(t\right)=\left[\text{min}\left\{\frac{d}{dt}\underset{\_}{F}\left(t;ƛ\right),\frac{d}{dt}\overline{F}\left(t;ƛ\right)\right\},\text{max}\left\{\frac{d}{dt}\underset{\_}{F}\left(t;ƛ\right),\frac{d}{dt}\overline{F}\left(t;ƛ\right)\right\}\right]$E3

In a similar manner, fractional order differential of $\tilde{F}\left(t\right)$ can be defined as, for all

$\in \left[0,1\right]$, if $\underset{\_}{F}\left(t;ƛ\right)$ and $\overline{F}\left(t;ƛ\right)$ are differentiable of order $\omega >0$, then $\tilde{F}\left(t\right)$ is differentiable of order $\omega >0$, i.e.

$D_t^{\omega }\tilde{F}\left(t\right)=\left[\text{min}\left\{D_t^{\omega }\underset{\_}{F}\left(t;ƛ\right),D_t^{\omega }\overline{F}\left(t;ƛ\right)\right\},\text{max}\left\{D_t^{\omega }\underset{\_}{F}\left(t;ƛ\right),D_t^{\omega }\overline{F}\left(t;ƛ\right)\right\}\right]$E4

where can be either fuzzy Riemann‐Liouville fractional differential operator or fuzzy Caputo-type fractional differential operator [15, 16, 19, 22, 23]. Here it is considered as fuzzy Caputo‐type fractional derivative that is approximated by Grünwald‐Letnikov's approach, illustrated in the next sequel.

2.3. System of fractional order fuzzy differential equations

In particular, modeling of differential equations of fractional order in imprecise characteristics is obtained by encompassing fuzzy‐valued functions. Let

, then fuzzy differential equation of fractional order $\omega \in \left(0,1\right]$, subjected to initial conditions, is structured as:

$D_t^{\omega }\tilde{X}\left(t\right)=\Psi \left(t,\tilde{X}\left(t\right)\right)$E5

$\tilde{X}\left(t_0\right)={\tilde{U}}_0$E6

where the unknown fuzzy‐valued function $\tilde{X}\left(t\right)$ can be written in form of

$$‐levels as, for all

$\in \left[0,1\right]$, $\tilde{X}\left(t\right)=\left[\underset{\_}{X}\left(t;ƛ\right),\overline{X}\left(t;ƛ\right)\right]$, where as

can be linear or nonlinear term in the form of fuzzy‐valued function and is the fuzzy number, which can also be expressed as

, for all

. Concisely, Eq. (5) is considered to have a unique and stable solution, for the reason that

is continuous and satisfies the Lipschitz condition, i.e. there exists $L>0$ such that for

$D\left(\Psi \left(t,\tilde{X}\right),\Psi \left(t,\tilde{W}\right)\right)\le L.D\left(\tilde{X},\tilde{W}\right)\forall \left(t,\tilde{X}\right),\left(t,\tilde{W}\right)\in R,\tilde{X},\tilde{W}\in E$E7

Many papers [14, 15, 22] comprise the theorems of stability and uniqueness of the solution of Eq. (5).

Here, we consider the system of fractional order fuzzy differential equations of the following form:

$\begin{array}{l}{D}_t^{{\omega }_1}{\tilde{X}}_1\left(t\right)=\Psi \left({\tilde{X}}_1\left(t\right),{\tilde{X}}_2\left(t\right),\cdots ,{\tilde{X}}_n\left(t\right)\right)\\ D_t^{{\omega }_2}{\tilde{X}}_2\left(t\right)=\Psi \left({\tilde{X}}_1\left(t\right),{\tilde{X}}_2\left(t\right),\cdots ,{\tilde{X}}_n\left(t\right)\right)\\ ⋮\\ D_t^{{\omega }_n}{\tilde{X}}_n\left(t\right)=\Psi \left({\tilde{X}}_1\left(t\right),{\tilde{X}}_2\left(t\right),\cdots ,{\tilde{X}}_n\left(t\right)\right)\end{array}$E8

with the initial conditions,

${\tilde{X}}_1\left(t_0\right)={\tilde{\nu }}_1,{\tilde{X}}_2\left(t_0\right)={\tilde{\nu }}_2,\dots ,{\tilde{X}}_n\left(t_0\right)={\tilde{\nu }}_n$E9

where ${\tilde{\nu }}_1,{\tilde{\nu }}_2,\dots ,{\tilde{\nu }}_n$ are the fuzzy numbers that can be written as, for all ${\tilde{\nu }}_n\left(ƛ\right)=\left[{\underset{\_}{\nu }}_n\left(ƛ\right),{\overline{\nu }}_n\left(ƛ\right)\right]$, n≥1, ω₁, ω₂, …, ω_n are the fractional orders such that ${\omega }_n\in \left(0,1\right]$ and the right hand side of Eq. (8) represent a system of fuzzy nonlinear equations with crisp coefficients $k_{ij}$, $i\ge 1,j\le n$, i.e.

$\Psi \left({\tilde{X}}_1\left(t\right),{\tilde{X}}_2\left(t\right),\cdots ,{\tilde{X}}_n\left(t\right)\right)={\displaystyle \sum _{j=1}^n{k}_{ij}{\tilde{X}}_j^m\left(t\right)},m\ge 1$E10

Therefore, Eq. (8) can be remodeled as:

$\begin{array}{l}{D}_t^{{\omega }_1}{\tilde{X}}_1\left(t\right)={\displaystyle \sum _{j=1}^n{k}_{1j}{\tilde{X}}_j^m\left(t\right)}=k_{11}{\tilde{X}}_1^m\left(t\right)\oplus k_{12}{\tilde{X}}_2^m\left(t\right)\oplus \cdots \oplus k_{1n}{\tilde{X}}_n^m\left(t\right),\\ D_t^{{\omega }_2}{\tilde{X}}_2\left(t\right)={\displaystyle \sum _{j=1}^n{k}_{2j}{\tilde{X}}_j^m\left(t\right)}=k_{21}{\tilde{X}}_1^m\left(t\right)\oplus k_{22}{\tilde{X}}_2^m\left(t\right)\oplus \cdots \oplus k_{2n}{\tilde{X}}_n^m\left(t\right),\\ ⋮\\ D_t^{{\omega }_2}{\tilde{X}}_2\left(t\right)={\displaystyle \sum _{j=1}^n{k}_{2j}{\tilde{X}}_j^m\left(t\right)}=k_{21}{\tilde{X}}_1^m\left(t\right)\oplus k_{22}{\tilde{X}}_2^m\left(t\right)\oplus \cdots \oplus k_{2n}{\tilde{X}}_n^m\left(t\right)\end{array}$E11

And as mentioned earlier,

are taken as the fuzzy Caputo‐type fractional differential operators and are numerically interpreted using Grünwald‐Letnikov’s fractional derivative definition.

3. Grünwald‐Letnikov’s fractional derivative

This section comprises the description of Grünwald‐Letnikov’s fractional derivative in conjunction with the algorithm to solve the system of Eq. (11) and undergoes some requisite theorem and lemma of the governing approach.

Consider a function

in finite interval [0, T], let the interval be divided into equidistant grids of step size h as:

$0={\eta }_0<{\eta }_1<\cdots <{\eta }_{\sigma }=t=\sigma h\text{with}{\eta }_{\sigma }-{\eta }_{\sigma -1}=h$E12

${}^{GL}D{}_t^{\omega }Y\left(t\right)=\underset{h\to 0}{\lim }\frac{1}{h^{\omega }}{\displaystyle \sum _{i=0}^{\left[\frac{t}{h}\right]}{\left(-1\right)}^i\left(\begin{array}{c}\omega \\ i\end{array}\right)}Y\left(t-ih\right)$E13

where $\left(\begin{array}{c}\omega \\ i\end{array}\right)$ are the binomial coefficients that are obtained by the formula:

$\left(\begin{array}{c}\omega \\ i\end{array}\right)=\frac{\Gamma \left(\omega +1\right)}{i!\Gamma \left(\omega -i+1\right)}$E14

and $\left[\frac{t}{h}\right]$ represents the integral part.

3.1. Lemma

Let

be a smooth function in $\left[0,T\right]$, such that it can be expressed as a power series for $\left[t\right]<T$, where $\left[t\right]$ is the integral part of $t$, then the Grünwald‐Letnikov’s approximation for each $0<t<T$, a series of step size $h$ and $t=\sigma h$ can be stated as:

${}^{GL}D{}_t^{\omega }Y\left(t\right)=\frac{1}{h^{\omega }}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}\omega \\ i\end{array}\right)}Y\left(t_{\sigma -i}\right)+O\left(h\right)\left(h\to 0\right)$E15

This definition is considered to be equivalent to the definition of Riemann‐Liouville fractional derivative and for equivalence to Caputo’s fractional definition the following term of initial value is added to the right hand side of Eq. (15), i.e.

${}^{GL}D{}_t^{\omega }Y\left(t\right)=\frac{1}{h^{\omega }}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}\omega \\ i\end{array}\right)}Y\left(t_{\sigma -i}\right)-\frac{t_{\sigma }^{-\omega }}{\Gamma \left(1-\omega \right)}Y\left(0\right)$E16

That becomes zero if initial values of Caputo‐type differential equations are homogeneous and again reduces to that of Riemann‐Liouville definition. Since here the fuzzy Caputo‐type fractional differential equations are considered with inhomogeneous initial values, the definition in Eq. (16) will be used for the approximation of Eq. (11).

Now let

be a fuzzy‐valued function such that

, then Grünwald‐Letnikov’s fractional derivative of

(t) is expressed as:

${}^{GL}D{}_t^{\omega }\tilde{F}\left(t\right)=\frac{1}{h^{\omega }}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}\omega \\ i\end{array}\right)}\tilde{F}\left(t_{\sigma -i}\right)\Theta \frac{t_{\sigma }^{-\omega }}{\Gamma \left(1-\omega \right)}\tilde{F}\left(0\right)$E17

and in

$$‐level sets it is sorted out as, for all

$\in \left[0,1\right]$,

${}^{GL}D{}_t^{\omega }\tilde{F}\left(t\right)=\left[\begin{array}{l}\frac{1}{h^{\omega }}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}\omega \\ i\end{array}\right)}\underset{\_}{F}\left(t_{\sigma -i},ƛ\right)-\frac{t_{\sigma }^{-\omega }}{\Gamma \left(1-\omega \right)}\underset{\_}{F}\left(0,ƛ\right),\\ \frac{1}{h^{\omega }}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}\omega \\ i\end{array}\right)}\overline{F}\left(t_{\sigma -i},ƛ\right)-\frac{t_{\sigma }^{-\omega }}{\Gamma \left(1-\omega \right)}\overline{F}\left(0,ƛ\right)\end{array}\right]$E18

Next consider the fractional system in Eq. (11), for the cases of inhomogeneous initial values. Assume the uniform grids $t_{\sigma }=\sigma h$, where $\sigma =1,\dots ,M$, such that $Mh=T,M\in$

. Applying Grünwald‐Letnikov’s fractional derivative on left hand sides of Eq. (11) we get,

$\begin{array}{l}\frac{1}{h^{{\omega }_1}}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}{\omega }_1\\ i\end{array}\right)}{\tilde{X}}_1\left(\left(\sigma -i\right)h\right)\Theta \frac{{\left(\sigma h\right)}^{-{\omega }_1}}{\Gamma \left(1-{\omega }_1\right)}{\tilde{X}}_1\left(0\right)={\displaystyle \sum _{j=1}^n{k}_{1j}{\tilde{X}}_j^m\left(\sigma h\right)},\\ \frac{1}{h^{{\omega }_2}}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}{\omega }_2\\ i\end{array}\right)}{\tilde{X}}_2\left(\left(\sigma -i\right)h\right)\Theta \frac{{\left(\sigma h\right)}^{-{\omega }_2}}{\Gamma \left(1-{\omega }_2\right)}{\tilde{X}}_2\left(0\right)={\displaystyle \sum _{j=1}^n{k}_{2j}{\tilde{X}}_j^m\left(\sigma h\right)},\\ ⋮\\ \frac{1}{h^{{\omega }_n}}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}{\omega }_n\\ i\end{array}\right)}{\tilde{X}}_n\left(\left(\sigma -i\right)h\right)\Theta \frac{{\left(\sigma h\right)}^{-{\omega }_n}}{\Gamma \left(1-{\omega }_n\right)}{\tilde{X}}_n\left(0\right)={\displaystyle \sum _{j=1}^n{k}_{nj}{\tilde{X}}_j^m\left(\sigma h\right)}\end{array}$E19

Solving above system fuzzy-valued functions of respective fuzzy functions are generated at different grid points.

3.2. Theorem: truncation error

Let fuzzy‐valued functions ${\tilde{X}}_1\left(t_{\sigma }\right),{\tilde{X}}_2\left(t_{\sigma }\right),\cdots ,{\tilde{X}}_n\left(t_{\sigma }\right)$ be the approximations to the true solutions ${\tilde{X}}_1\left(t_{\sigma }\right),{\tilde{X}}_2\left(t_{\sigma }\right),\cdots ,{\tilde{X}}_n\left(t_{\sigma }\right)$, respectively and consider $\Psi$ satisfies Lipchitz condition, then the local truncation error of the proposed numerical approach is $O\left(h^{1+{\omega }_n}\right)$, for $n\ge 1$, i.e.

$\begin{array}{l}{\tilde{X}}_1\left(t_{\sigma }\right)\Theta {\tilde{X}}_1\left(t_{\sigma }\right)=O\left(h^{1+{\omega }_1}\right),\\ {\tilde{X}}_2\left(t_{\sigma }\right)\Theta {\tilde{X}}_2\left(t_{\sigma }\right)=O\left(h^{1+{\omega }_2}\right),\\ ⋮\\ {\tilde{X}}_n\left(t_{\sigma }\right)\Theta {\tilde{X}}_n\left(t_{\sigma }\right)=O\left(h^{1+{\omega }_n}\right).\end{array}$E20

Proof:

Assume the nth equation of the system (19) and on applying Grünwald‐Letnikov’s fractional derivative we have,

$\frac{1}{h^{{\omega }_n}}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}{\omega }_n\\ i\end{array}\right)}{\tilde{X}}_n\left(t_{\sigma -i}\right)\Theta \frac{t_{\sigma }{}^{-{\omega }_n}}{\Gamma \left(1-{\omega }_n\right)}{\tilde{X}}_n\left(0\right)=\Psi \left({\tilde{X}}_1\left(t_{\sigma }\right),{\tilde{X}}_2\left(t_{\sigma }\right),\cdots ,{\tilde{X}}_n\left(t_{\sigma }\right)\right)$E21

for $n\ge 1$ and from Lemma 3.1 we can attain,

$\frac{1}{h^{{\omega }_n}}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}{\omega }_n\\ i\end{array}\right)}{\tilde{X}}_n\left(t_{\sigma -i}\right)\Theta \frac{t_{\sigma }{}^{-{\omega }_n}}{\Gamma \left(1-{\omega }_n\right)}{\tilde{X}}_n\left(0\right)+O\left(h\right)=\Psi \left({\tilde{X}}_1\left(t_{\sigma }\right),{\tilde{X}}_2\left(t_{\sigma }\right),\cdots ,{\tilde{X}}_n\left(t_{\sigma }\right)\right)$E22

Subtracting Eq. (22) from Eq. (21),

$\begin{array}{l}\frac{1}{h^{{\omega }_n}}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}{\omega }_n\\ i\end{array}\right)}{\tilde{X}}_n\left(t_{\sigma -i}\right)\Theta \frac{t_{\sigma }{}^{-{\omega }_n}}{\Gamma \left(1-{\omega }_n\right)}{\tilde{X}}_n\left(0\right)\Theta \frac{1}{h^{{\omega }_n}}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}{\omega }_n\\ i\end{array}\right)}{\tilde{X}}_n\left(t_{\sigma -i}\right)\\ \oplus \frac{t_{\sigma }{}^{-{\omega }_n}}{\Gamma \left(1-{\omega }_n\right)}{\tilde{X}}_n\left(0\right)+O\left(h\right)=\Psi \left({\tilde{X}}_1\left(t_{\sigma }\right),{\tilde{X}}_2\left(t_{\sigma }\right),\cdots ,{\tilde{X}}_n\left(t_{\sigma }\right)\right)\\ \Theta \Psi \left({\tilde{X}}_1\left(t_{\sigma }\right),{\tilde{X}}_2\left(t_{\sigma }\right),\cdots ,{\tilde{X}}_n\left(t_{\sigma }\right)\right)\end{array}$E23

Let, for $i=0,1,\dots \sigma -1$, ${\tilde{X}}_n\left(t_{\sigma }\right)={\tilde{X}}_n\left(t_{\sigma }\right)$, then on further manipulation we get,

$\begin{array}{l}\frac{1}{h^{{\omega }_n}}\left[{\tilde{X}}_n\left(t_{\sigma }\right)\Theta {\tilde{X}}_n\left(t_{\sigma }\right)\right]+O\left(h\right)=\Psi \left({\tilde{X}}_1\left(t_{\sigma }\right),{\tilde{X}}_2\left(t_{\sigma }\right),\cdots ,{\tilde{X}}_n\left(t_{\sigma }\right)\right)\\ \Theta \Psi \left({\tilde{X}}_1\left(t_{\sigma }\right),{\tilde{X}}_2\left(t_{\sigma }\right),\cdots ,{\tilde{X}}_n\left(t_{\sigma }\right)\right)\end{array}$E24

or it can be rearranged as:

$\begin{array}{l}\left[{\tilde{X}}_n\left(t_{\sigma }\right)\Theta {\tilde{X}}_n\left(t_{\sigma }\right)\right]=h^{{\omega }_n}D\left[\Psi \left({\tilde{X}}_1\left(t_{\sigma }\right),{\tilde{X}}_2\left(t_{\sigma }\right),\cdots ,{\tilde{X}}_n\left(t_{\sigma }\right)\right),\Psi \left({\tilde{X}}_1\left(t_{\sigma }\right),{\tilde{X}}_2\left(t_{\sigma }\right),\cdots ,{\tilde{X}}_n\left(t_{\sigma }\right)\right)\right]\\ +O\left(h^{1+{\omega }_n}\right)\end{array}$E25

where

defines Hausdroff distance. On using Lipschitz condition, i.e. Eq. (7), proof is completed by obtaining the following equation:

$\left(1-L_n{h}^{{\omega }_n}\right)\left[{\tilde{X}}_n\left(t_{\sigma }\right)\Theta {\tilde{X}}_n\left(t_{\sigma }\right)\right]\le O\left(h^{1+{\omega }_n}\right)\forall n\ge 1$E26

4. Numerical illustrations

Subsequent to the algorithm demonstrated in Section 3, here numerical experiments of some system of fuzzy fractional differential equations are presented. Results for fuzzy‐valued functions are depicted in tabular form in the finite interval $\left[0,1\right]$ at different values of $\omega \in \left(0,1\right]$. In addition, error bar pictorials are given for each respective example. All the exact values and calculations are carried out through Mathematica 10.

4.1. Example 1

Following nonlinear fractional system is solved in [27] using homotopy analysis method, here the system is restructured with imprecise functions

and

as:

$D_t^{{\omega }_1}{\tilde{X}}_1\left(t\right)=0.5{\tilde{X}}_1\left(t\right)$

$D_t^{{\omega }_2}{\tilde{X}}_2\left(t\right)={\tilde{X}}_2\left(t\right)\oplus {\tilde{X}}_1^2\left(t\right)$E27

with ${\omega }_1,{\omega }_2\in \left(0,1\right]$ and subjected to initial conditions

${\tilde{X}}_1\left(0\right)=\left[0.75+0.25ƛ,1.125-0.125ƛ\right],{\tilde{X}}_2\left(0\right)=\left[ƛ-1,1-ƛ\right]$E28

On applying Grünwald‐Letnikov’s fractional definition on left hand side of Eq. (27) and following the algorithm, the differential equations are reduced to nonlinear algebraic equations as:

$\frac{1}{h^{{\omega }_1}}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}{\omega }_1\\ i\end{array}\right)}{\tilde{X}}_1\left(\left(\sigma -i\right)h\right)\Theta \frac{{\left(\sigma h\right)}^{-{\omega }_1}}{\Gamma \left(1-{\omega }_1\right)}{\tilde{X}}_1\left(0\right)=0.5{\tilde{X}}_1\left(\sigma h\right),$

$\frac{1}{h^{{\omega }_2}}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}{\omega }_2\\ i\end{array}\right)}{\tilde{X}}_2\left(\left(\sigma -i\right)h\right)\Theta \frac{{\left(\sigma h\right)}^{-{\omega }_2}}{\Gamma \left(1-{\omega }_2\right)}{\tilde{X}}_2\left(0\right)={\tilde{X}}_2\left(\sigma h\right)\oplus {\tilde{X}}_1^2\left(\sigma h\right)$E29

which on expanding to

‐levels of ${\tilde{X}}_1\left(t\right)$ and ${\tilde{X}}_2\left(t\right)$ convert into system of four nonlinear equations, i.e. for all

$\in \left[0,1\right]$,

$\frac{1}{h^{{\omega }_1}}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}{\omega }_1\\ i\end{array}\right)}\underset{\_}{X_1}\left(\left(\sigma -i\right)h;ƛ\right)-\frac{{\left(\sigma h\right)}^{-{\omega }_1}}{\Gamma \left(1-{\omega }_1\right)}\underset{\_}{X_1}\left(0;ƛ\right)=0.5\underset{\_}{X_1}\left(\sigma h;ƛ\right),$E900

$\frac{1}{h^{{\omega }_1}}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}{\omega }_1\\ i\end{array}\right)}\overline{X_1}\left(\left(\sigma -i\right)h;ƛ\right)-\frac{{\left(\sigma h\right)}^{-{\omega }_1}}{\Gamma \left(1-{\omega }_1\right)}\overline{X_1}\left(0;ƛ\right)=0.5\overline{X_1}\left(\sigma h;ƛ\right),$

$\frac{1}{h^{{\omega }_2}}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}{\omega }_2\\ i\end{array}\right)}\underset{\_}{X_2}\left(\left(\sigma -i\right)h;ƛ\right)-\frac{{\left(\sigma h\right)}^{-{\omega }_2}}{\Gamma \left(1-{\omega }_2\right)}\underset{\_}{X_2}\left(0;ƛ\right)=\underset{\_}{X_2}\left(\sigma h;ƛ\right)+\underset{\_}{X_1^2}\left(\sigma h;ƛ\right),$

$\frac{1}{h^{{\omega }_2}}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}{\omega }_2\\ i\end{array}\right)}\overline{X_2}\left(\left(\sigma -i\right)h;ƛ\right)-\frac{{\left(\sigma h\right)}^{-{\omega }_2}}{\Gamma \left(1-{\omega }_2\right)}\overline{X_2}\left(0;ƛ\right)=\overline{X_2}\left(\sigma h;ƛ\right)+\overline{X_1^2}\left(\sigma h;ƛ\right)$E30

Table 1.

Numerical results and absolute errors of ${\tilde{X}}_1\left(t\right)$ for Example 1 at ${\omega }_1=1$, ${\omega }_2=1$, $h=0.001$ and $t=1$.

Solving this system, numerical approximations of Eq. (27) are obtained. Tables 1 and 2 represent absolute error of ${\tilde{X}}_1\left(t\right)$ and ${\tilde{X}}_2\left(t\right)$, respectively, for ${\omega }_1={\omega }_2=1$, $h=0.001$, $t=1$ and at different values of

, whereas Table 3 shows the approximations of ${\tilde{X}}_1\left(t\right)$ and ${\tilde{X}}_2\left(t\right)$ for ${\omega }_1=0.95$, ${\omega }_2=0.87,h=0.1$ and $t=1$, at different values of

. In Figures 1 and 2, the pointwise error variations of ${\tilde{X}}_1\left(t\right)$ and ${\tilde{X}}_2\left(t\right)$, accordingly, at each time within the given interval for ${\omega }_1={\omega }_2=1,h=0.1$ and

= 0.6, are plotted. In these graphs, each approximated point is plotted against the value of $\sigma$ in a discrete manner and each bar line on respective approximated point illustrates the measure of the absolute error at that point. Absolute error is obtained by taking the point‐to‐point difference between exact and the solutions calculated by Grünwald‐Letnikov’s fractional approach. Since these variations show small differences, this implies our results are in good agreement with the exact solutions.

Table 2.

Numerical results and absolute errors of ${\tilde{X}}_2\left(t\right)$ for Example 1 at ${\omega }_1=1$, ${\omega }_2=1$, $h=0.001$ and $t=1$.

Table 3.

Approximations of ${\tilde{X}}_1\left(t\right)$ and ${\tilde{X}}_2\left(t\right)$ of Example 1 for ${\omega }_1=0.95$, ${\omega }_2=0.87$, $h=0.1$ and $t=1$.

Figure 1.

Bar plot of $\sigma$ of ${\tilde{X}}_1\left(t\right)$ of Example 1 for $h=0.1,{\omega }_1={\omega }_2=1$ and = 0.6.

Figure 2.

Bar plot of approximate solutions and absolute error versus $\sigma$ of ${\tilde{X}}_2\left(t\right)$ of Example 1 for $h=0.1,{\omega }_1={\omega }_2=1$ and = 0.6.

4.2. Example 2

Consider the following nonlinear fractional system [27] with imprecise functions ${\tilde{X}}_1\left(t\right),{\tilde{X}}_2\left(t\right)$, and ${\tilde{X}}_3\left(t\right)$ as:

$D_t^{{\omega }_1}{\tilde{X}}_1\left(t\right)={\tilde{X}}_1\left(t\right),$

$D_t^{{\omega }_2}{\tilde{X}}_2\left(t\right)=2{\tilde{X}}_1^2\left(t\right)$

$D_t^{{\omega }_3}{\tilde{X}}_3\left(t\right)=3{\tilde{X}}_1\left(t\right)•{\tilde{X}}_2\left(t\right)$E31

with ${\omega }_1,{\omega }_2,{\omega }_3\in \left(0,1\right]$ and subjected to initial conditions

${\tilde{X}}_1\left(0\right)={\tilde{X}}_2\left(0\right)=\left[0.75+0.25ƛ,1.125-0.125ƛ\right],{\tilde{X}}_3\left(0\right)=\left[ƛ-1,1-ƛ\right]$E32

On employing Grünwald‐Letnikov’s approach, the differential equations are converted into nonlinear algebraic equations as:

$\frac{1}{h^{{\omega }_1}}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}{\omega }_1\\ i\end{array}\right)}{\tilde{X}}_1\left(\left(\sigma -i\right)h\right)\Theta \frac{{\left(\sigma h\right)}^{-{\omega }_1}}{\Gamma \left(1-{\omega }_1\right)}{\tilde{X}}_1\left(0\right)={\tilde{X}}_1\left(\sigma h\right),$

$\frac{1}{h^{{\omega }_2}}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}{\omega }_2\\ i\end{array}\right)}{\tilde{X}}_2\left(\left(\sigma -i\right)h\right)\Theta \frac{{\left(\sigma h\right)}^{-{\omega }_2}}{\Gamma \left(1-{\omega }_2\right)}{\tilde{X}}_2\left(0\right)=2{\tilde{X}}_1^2\left(\sigma h\right)$

$\frac{1}{h^{{\omega }_3}}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}{\omega }_3\\ i\end{array}\right)}{\tilde{X}}_3\left(\left(\sigma -i\right)h\right)\Theta \frac{{\left(\sigma h\right)}^{-{\omega }_3}}{\Gamma \left(1-{\omega }_3\right)}{\tilde{X}}_3\left(0\right)=3{\tilde{X}}_1\left(\sigma h\right)•{\tilde{X}}_2\left(\sigma h\right)$E33

and in

‐levels of ${\tilde{X}}_1\left(t\right),{\tilde{X}}_2\left(t\right)$,, and ${\tilde{X}}_3\left(t\right)$ the system above converts into six nonlinear equations, i.e. for all

,

$\frac{1}{h^{{\omega }_1}}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}{\omega }_1\\ i\end{array}\right)}\underset{\_}{X_1}\left(\left(\sigma -i\right)h;ƛ\right)-\frac{{\left(\sigma h\right)}^{-{\omega }_1}}{\Gamma \left(1-{\omega }_1\right)}\underset{\_}{X_1}\left(0;ƛ\right)=\underset{\_}{X_1}\left(\sigma h;ƛ\right),$

$\frac{1}{h^{{\omega }_1}}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}{\omega }_1\\ i\end{array}\right)}\overline{X_1}\left(\left(\sigma -i\right)h;ƛ\right)-\frac{{\left(\sigma h\right)}^{-{\omega }_1}}{\Gamma \left(1-{\omega }_1\right)}\overline{X_1}\left(0;ƛ\right)=\overline{X_1}\left(\sigma h;ƛ\right),$

$\frac{1}{h^{{\omega }_2}}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}{\omega }_2\\ i\end{array}\right)}\underset{\_}{X_2}\left(\left(\sigma -i\right)h;ƛ\right)-\frac{{\left(\sigma h\right)}^{-{\omega }_2}}{\Gamma \left(1-{\omega }_2\right)}\underset{\_}{X_2}\left(0;ƛ\right)=2\underset{\_}{X_1^2}\left(\sigma h;ƛ\right),$

$\frac{1}{h^{{\omega }_2}}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}{\omega }_2\\ i\end{array}\right)}\overline{X_2}\left(\left(\sigma -i\right)h;ƛ\right)-\frac{{\left(\sigma h\right)}^{-{\omega }_2}}{\Gamma \left(1-{\omega }_2\right)}\overline{X_2}\left(0;ƛ\right)=2\overline{X_1^2}\left(\sigma h;ƛ\right),$

$\frac{1}{h^{{\omega }_3}}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}{\omega }_3\\ i\end{array}\right)}\underset{\_}{X_3}\left(\left(\sigma -i\right)h;ƛ\right)-\frac{{\left(\sigma h\right)}^{-{\omega }_3}}{\Gamma \left(1-{\omega }_3\right)}\underset{\_}{X_3}\left(0;ƛ\right)=3\underset{\_}{X_1}\left(\sigma h;ƛ\right)\underset{\_}{X_2}\left(\sigma h;ƛ\right)$

$\frac{1}{h^{{\omega }_3}}{\displaystyle \sum _{i=0}^{\sigma }{\left(-1\right)}^i\left(\begin{array}{c}{\omega }_3\\ i\end{array}\right)}\overline{X_3}\left(\left(\sigma -i\right)h;ƛ\right)-\frac{{\left(\sigma h\right)}^{-{\omega }_3}}{\Gamma \left(1-{\omega }_3\right)}\overline{X_3}\left(0;ƛ\right)=3\overline{X_1}\left(\sigma h;ƛ\right)\overline{X_2}\left(\sigma h;ƛ\right)$E34

Thus, numerical results of Eq. (31) are obtained from the above system. Tables 4–6 present absolute error of ${\tilde{X}}_1\left(t\right),{\tilde{X}}_2\left(t\right)$, and ${\tilde{X}}_3\left(t\right)$, respectively, for ${\omega }_1={\omega }_2={\omega }_3=1$, $h=0.001$, $t=1$ and at different values of

. In Table 7, the approximations of ${\tilde{X}}_1\left(t\right),{\tilde{X}}_2\left(t\right)$, and ${\tilde{X}}_3\left(t\right)$ are rendered for $h=0.1,{\omega }_1=0.95$, ${\omega }_2=0.87$, ${\omega }_3=0.79$ and $t=1$, at different values of

. Additionally, the pointwise error variations between approximated and exact solutions of ${\tilde{X}}_1\left(t\right),{\tilde{X}}_2\left(t\right)$, and ${\tilde{X}}_3\left(t\right)$ at each time within the given interval for ${\omega }_1={\omega }_2={\omega }_3=1$ and

= 0.6 are plotted in Figures 3–5, respectively. It is to be noted that the small length of bar lines on each point is illustrating small differences between the exact and the result obtained by the proposed approach that shows the acceptable convergence of the solution towards the exact values.

Table 4.

Numerical results and absolute errors of

for Example 2 at ${\omega }_1={\omega }_2={\omega }_3=1$, $h=0.001$ and $t=1$.

Table 5.

Numerical results and absolute errors of

for Example 2 at ${\omega }_1={\omega }_2={\omega }_3=1$, $h=0.001$ and $t=1$.

Table 6.

Numerical results and absolute errors of ${\tilde{X}}_3\left(t\right)$ for Example 2 at ${\omega }_1={\omega }_2={\omega }_3=1$, $h=0.001$ and $t=1$.

Table 7.

Approximations of ${\tilde{X}}_1\left(t\right),{\tilde{X}}_2\left(t\right)$, and ${\tilde{X}}_3\left(t\right)$ of Example 2 for ${\omega }_1=0.95$, ${\omega }_2=0.87$ ${\omega }_3=0.79$, $h=0.1$ and $t=1$.

Figure 3.

Bar plot of approximate solutions and absolute error versus $\sigma$ of ${\tilde{X}}_1\left(t\right)$ of Example 2 for $h=0.1,{\omega }_1={\omega }_2={\omega }_3=1$ and = 0.6.

Figure 4.

Bar plot of approximate solutions and absolute error versus $\sigma$ of ${\tilde{X}}_2\left(t\right)$ of Example 2 for $h=0.1,{\omega }_1={\omega }_2={\omega }_3=1$ and = 0.6.

Figure 5.

Bar plot of approximate solutions and absolute error versus $\sigma$ of ${\tilde{X}}_3\left(t\right)$ of Example 2 for $h=0.1,{\omega }_1={\omega }_2={\omega }_3=1$ and = 0.6.

5. Conclusion

In this chapter, system of fractional differential equations with fuzzy‐valued functions was constructed to study the system in imprecise environment. We assessed numerical interpretations of the system using Grünwald‐Letnikov’s fractional derivative scheme, which has not been considered for fuzzy differential equations in literature hitherto. In addition, we illustrated the stability of the scheme for the system of fuzzy fractional differential equations. Furthermore, we conducted experiment on some nonlinear fuzzy fractional systems and successfully attained the approximated solutions. From the entire discussion and analysis, collectively, we come up with the following remarks:

• Scrutinizing differential models with arbitrary fractional order in combination with fuzzy theory is effectively advantageous to analyze the change in the system at each fractional step with imprecise parameters rather than crisp values.

• Grünwald‐Letnikov’s fractional definition is equivalent to either Riemann‐Liouville fractional definition or Caputo-type fractional definition in case of homogeneous and inhomogeneous initial values, respectively. Since Riemann‐Liouville fractional definition and Caputo‐type fractional definition are greatly applicable for defining fractional derivative of fuzzy‐valued functions, so is Grünwald‐Letnikov’s fractional definition found to be.

• Approximations of examples attained by undertaking Grünwald‐Letnikov’s fractional derivative approach are efficaciously convergent towards the exact solutions that prove the method to be appropriate for the solutions of fuzzy differential equations of fractional order to a great extent.

• Pointwise explanation of errors through bar graph is conspicuously helpful in locating the error between exact and calculated solutions at each time by simply measuring the length of the bar at the respective point.