# Numerical Solution of First-Order Linear Differential Equations in Fuzzy Environment by Runge-Kutta-Fehlberg Method and Its Application

Authors: Sankar Prasad Mondal; Susmita Roy; Biswajit Das
Journal: International Journal of Differential Equations
 (2016)

---

## Abstract

The numerical algorithm for solving “first-order linear differential equation in fuzzy environment” is discussed. A scheme, namely, “Runge-Kutta-Fehlberg method,” is described in detail for solving the said differential equation. The numerical solutions are compared with (i)-gH and (ii)-gH differential (exact solutions concepts) system. The method is also followed by complete error analysis. The method is illustrated by solving an example and an application.


---

## Body

## 1. Introduction

Fuzzy Differential Equation. In modeling of real natural phenomena, differential equations play an important role in many areas of discipline, exemplary in economics, biomathematics, science, and engineering. Many experts in such areas widely use differential equations in order to make some problems under study more comprehensible. In many cases, information about the physical phenomena related is always immanent with uncertainty.Today, the study of differential equations with uncertainty is instantaneously growing as a new area in fuzzy analysis. The terms such as “fuzzy differential equation” and “fuzzy differential inclusion” are used interchangeably in mention to differential equations with fuzzy initial values or fuzzy boundary values or even differential equations dealing with functions on the space of fuzzy numbers. In the year 1987, the term “fuzzy differential equation (FDE)” was introduced by Kandel and Byatt [1]. There are different approaches to discuss the FDEs: (i) the Hukuhara derivative of a fuzzy number valued function is used, (ii) Hüllermeier [2] and Diamond and Watson [3–5] suggested a different formulation for the fuzzy initial value problems (FIVP) based on a family of differential inclusions, (iii) in [6, 7], Bede et al. defined generalized differentiability of the fuzzy number valued functions and studied FDE, and (iv) applying a parametric representation of fuzzy numbers, Chen et al. [8] established a new definition for the differentiation of a fuzzy valued function and used it in FDE.Solution of Fuzzy Differential Equation by Numerical Techniques. Numerical methods are the methods by which we can find the solution of differential equation where the exact solution is critical to find. There exist various numerical methods for solving differential equation such as Setia et al. [9], Liu [10], and Setia et al. [11]. Our aim is to find the numerical techniques by which the solution of a linear or nonlinear first-order fuzzy differential equation comes easily and the solution is very close to the exact solution. There exist many techniques of numerical methods for finding the solution of fuzzy differential equation. Authors applied the method in certain types of fuzzy differential equation which shows that their techniques are best fit for that particular problem. The first paper on fuzzy differential equation and numerical analysis was published in 1999 by Ma et al. [12]. Allahviranloo et al. [13] apply the two-step method on fuzzy differential equations. Allahviranloo et al. [14] find the numerical solution by using predictor-corrector method. Allahviranloo et al. [15] find an algorithm for finding the solution Nth-order fuzzy linear differential equations using numerical techniques. Pederson and Sambandham in [16] use characterization theorem on hybrid fuzzy initial value problem. A soft computing technique, namely, artificial neural network, is implicated for solving FDE by Effati and Pakdaman [17]. Duraisamy and Usha [18] used modified Euler’s method. The extension principle method was compared by Euler’s method in Saberi Najafi et al. [19] article. Rostami et al. [20] find a numerical algorithm for solving nonlinear fuzzy differential equations. Moghadam and Dahaghin [21] apply two-step methods for numerical solution of FDE. Batiha in [22] finds an iterative solution of multispecies predator-prey model by variational iteration method. Ahmad and Hasan [23] proposed a new fuzzy version of Euler’s method for solving differential equations with fuzzy initial value. Nirmala and Chenthur Pandian [24] give an idea for improving the numerical result on FDE. Shafiee et al. [25] use predictor-corrector method for nonlinear fuzzy Volterra integral equations. Comparison results on some numerical techniques on first-order fuzzy differential equation are illustrated by Ghanbari [26]. The use of variational iteration method for solving Nth-order fuzzy differential equations is shown by Jafari et al. [27]. Tapaswini and Chakraverty [28] discuss a new approach to fuzzy initial value problem by Improved Euler method. The solution of FIVP is compared by Least Square method and Adomian Decomposition method by Ahmed and Fadhel [29]. Solution of differential equation by Euler’s method using fuzzy concept is developed by Saikia [30]. Ezzati et al. [31] find the numerical solution of Volterra-Fredholm integral equations with the help of inverse and direct discrete fuzzy transforms and collocation technique. The Adomian method is applied on second-order FDE by Wang and Guo [32]. Fard [33] uses iterative scheme to find the solution of generalized system of linear FDE, whereas Block method is used by Mehrkanoon et al. [34]. Asady and Alavi [35] apply a numerical method for solving Nth-order linear fuzzy differential equation.Solution of Fuzzy Differential Equation by Runge-Kutta Method. Runge-Kutta method is well known for finding the approximate or numerical solution. In the last decade Runge-Kutta method is applied in fuzzy differential equation for finding the numerical solution. The researchers are giving various types of view to apply these methods. Someone changes the order and someone applies different types on FDE, a comparison of another method to Runge-Kutta method. The details of published work done in Runge-Kutta method are summarized below.Numerical Solution of Fuzzy Differential Equations by Runge-Kutta method of order three is developed by Duraisamy and Usha [36]. Solution techniques for fourth-order Runge-Kutta method with higher order derivative approximations are developed by Nirmala and Chenthur Pandian [37]. Runge-Kutta method of order five is developed by Jayakumar et al. [38]. The techniques extended Runge-Kutta-like formulae of order four are developed by Ghazanfari and Shakerami [39]. Third-order Runge-Kutta method is developed by Kanagarajan and Sambath [40]. Runge-Kutta-Fehlberg method for hybrid fuzzy differential equation is solved by Jayakumar and Kanagarajan [41]. A different approach followed by Runge-Kutta method is applied by Akbarzadeh Ghanaie and Mohseni Moghadam [42]. “Numerical Solution of Fuzzy IVP with Trapezoidal and Triangular Fuzzy Numbers by Using Fifth-Order Runge-Kutta Method” is solved by Ghanbari [43]. “New Multi-Step Runge-Kutta Method for Solving Fuzzy Differential Equations” is solved by Nirmala and Chenthur Pandian [44]. “Numerical Solution of Fuzzy Hybrid Differential Equation by Third Order Runge-Kutta Nystrom Method” is solved by Saveetha and Chenthur Pandian [45]. A new approach to solve fuzzy differential equation by using third-order Runge-Kutta method is developed by Deshmukh [46]. Runge-Kutta method of order four is developed by Duraisamy and Usha [47] and order five is developed by Jayakumar and Kanagarajan [48].Application of Fuzzy Differential Equation. Fuzzy differential equations play a significant role in the fields of biology, engineering, and physics as well as among other fields of science, for example, in population models [49], civil engineering [50], bioinformatics and computational biology [51], quantum optics and gravity [52], modeling hydraulic [53], HIV model [54], decay model [55], predator-prey model [56], population dynamics model [57], friction model [58], growth model [59], bacteria culture model [60], bank account and drug concentration problem [61], barometric pressure problem [62], concentration problem [63], weight loss and oil production model [64], arm race model [65], vibration of mass [66], and fractional predator-prey equation [67].Novelties. Although some developments are done, some new interest and new work have been done by ourselves which are mentioned below:(i)
Differential equation is solved in fuzzy environment by numerical techniques; that is, coefficients and initial condition both are taken as fuzzy number of a differential equation and solved by numerical techniques.
(ii)
The numerical solution is compared with the exact solution ((i)-gH and (ii)-gH both cases).
(iii)
Runge-Kutta-Fehlberg method for solving fuzzy differential equation is used.
(iv)
For application purpose a mixture problem is considered.
(v)
The solutions are found using different step length for better accuracy of the result.
(vi)
The necessary algorithm for numerical solution is given.Structure of the Paper. The paper is organized as follows: in Preliminary Concepts, the preliminary concepts and basic concepts on fuzzy number and fuzzy derivative are given. The method for finding the exact solution is discussed in Exact Solution of Fuzzy Differential Equation. In Numerical Solution of Fuzzy Differential Equation we proposed Runge-Kutta-Fehlberg method in fuzzy environment. The convergence of the said method and algorithm for finding the numerical results are also discussed in this section. Numerical Example shows a numerical example. In Application an important application, namely, mixture problem, is illustrated in fuzzy environment. Finally conclusions and future research scope of this paper are drawn in last section, Conclusion.


## 2. Preliminary Concepts

Definition 1 (fuzzy set).
A fuzzy setA
~ is defined by A
~
=
x
,
μ
A
~
x
:
x
∈
A
,
μ
A
~
x
∈
[
0,1
]. In the pair x
,
μ
A
~
x the first element x belongs to the classical set A and the second element μ
A
~
x belongs to the interval [
0,1
], called membership function.Definition 2 (α-cut of a fuzzy set).
Theα-level set (or interval of confidence at level α or α-cut) of the fuzzy set A
~ of X is a crisp set A
α that contains all the elements of X that have membership values in A greater than or equal to α; that is, A
~
=
x
:
μ
A
~
x
≥
α
,
x
∈
X
,
0
<
α
≤
1.Definition 3 (fuzzy number).
The basic definition of fuzzy number is as follows [30]: if we denote the set of all real numbers by R and the set of all fuzzy numbers on R is indicated by R
F then a fuzzy number is mapping such that u
:
R
→
[
0,1
], which satisfies the following four properties:(i)
u is upper semicontinuous.
(ii)
u is a fuzzy convex; that is, u
λ
x
+
1
-
λ
y
≥
min
⁡
u
x
,
u
y
⁡ for all x
,
y
∈
R
,
λ
∈
[
0,1
].
(iii)
u is normal; that is, ∃
x
0
∈
R for which u
x
0
=
1.
(iv)
supp
⁡
u
=
{
x
∈
R
∣
u
x
>
0
} is support of u and the closure of (supp
⁡
u) is compact.Definition 4 (parametric form of fuzzy number [31]).
A fuzzy number is represented by an ordered pair of functions(
u
1
α
,
u
2
(
α
)
), 0
≤
α
≤
1, that satisfy the following condition:(1)
u
1
α is a bounded left continuous nondecreasing function for any α
∈
[
0,1
].
(2)
u
2
(
α
) is a bounded left continuous nonincreasing function for any α
∈
[
0,1
].
(3)
u
1
α
≤
u
2
(
α
) for any α
∈
[
0,1
].Note. If u
1
α
=
u
2
α
=
α, then α is a crisp number.Definition 5 (generalized Hukuhara difference [20]).
The generalized Hukuhara difference of two fuzzy numbersu
,
v
∈
R
F is defined as follows:(1)
u
⊖
g
H
v
=
w
⟺
i
u
=
v
⊕
w
,
or
ii
v
=
u
⊕
-
1
w
.Consider w
α
=
w
1
α
,
w
1
(
α
); then w
1
α
=
min
⁡
u
1
α
-
v
1
α
,
u
2
(
α
)
-
v
2
(
α
) and w
2
α
=
max
⁡
u
1
α
-
v
1
α
,
u
2
(
α
)
-
v
2
(
α
).Here the parametric representation of a fuzzy valued functionf
:
[
a
,
b
]
→
R
F is expressed by f
(
t
)
α
=
f
1
t
,
α
,
f
2
(
t
,
α
), t
∈
a
,
b
,
α
∈
[
0,1
].Definition 6 (generalized Hukuhara derivative for first order [20]).
The generalized Hukuhara derivative of a fuzzy valued functionf
:
(
a
,
b
)
→
R
F at t
0 is defined as (2)
f
′
t
0
=
lim
h
→
0
⁡
f
t
0
+
h
⊖
g
H
f
t
0
h
.If f
′
t
0
∈
R
F satisfying (2) exists, we say that f is generalized Hukuhara differentiable at t
0.Also we say thatf
(
t
) is (i)-gH differentiable at t
0 if (3)
f
′
t
0
α
=
f
1
′
t
0
,
α
,
f
2
′
t
0
,
α
,and f
(
t
) is (ii)-gH differentiable at t
0 if (4)
f
′
t
0
α
=
f
2
′
t
0
,
α
,
f
1
′
t
0
,
α
.Definition 7 (see [6]).
For arbitraryu
=
u
1
,
u
2 and v
=
v
1
,
v
2
∈
E
1, the quantity(5)
D
u
,
v
=
∫
0
1
u
1
-
v
1
2
+
∫
0
1
u
2
-
v
2
2
1
/
2is the distance between fuzzy numbers u and v.Definition 8 (triangular fuzzy number).
A triangular fuzzy number (TFN) denoted byA
~ is defined as a
,
b
,
c where the membership function(6)
μ
A
~
x
=
0
,
x
≤
a
,
x
-
a
b
-
a
,
a
≤
x
≤
b
,
1
,
x
=
b
,
c
-
x
c
-
b
,
b
≤
x
≤
c
,
0
,
x
≥
c
.Definition 9 (α-cut of a fuzzy set A
~).
Theα-cut of A
~
=
(
a
,
b
,
c
) is given by(7)
A
α
=
a
+
α
b
-
a
,
c
-
α
c
-
b
,
∀
α
∈
0,1
.Definition 10 (fuzzy ordinary differential equation (FODE)).
Consider a simple 1st-order linear ordinary differential equation as follows:(8)
d
x
d
t
=
k
x
+
x
0
with  initial  condition
x
t
0
=
γ
.The above ordinary differential equation is called fuzzy ordinary differential equation if any one of the following three cases holds:(i)
Onlyγ is a fuzzy number (Type-I).
(ii)
Onlyk is a fuzzy number (Type-II).
(iii)
Bothk and γ are fuzzy numbers (Type-III).


## 3. Exact Solution of Fuzzy Differential Equation

Consider the fuzzy initial value problem(9)
y
′
t
=
f
t
,
y
t
,
t
∈
I
=
0
,
T
with
y
0
=
y
0
,where f is a continuous mapping from R
+
×
R into R and y
0
∈
E with r-level sets(10)
y
0
r
=
y
1
0
;
α
,
y
2
0
;
α
,
α
∈
0,1
.We write f
t
,
y
=
[
f
1
t
,
y
,
f
2
(
t
,
y
)
] and f
1
t
,
y
=
F
[
t
,
y
1
,
y
2
], f
2
(
t
,
y
)
=
G
[
t
,
y
1
,
y
2
].Because ofy
′
t
=
f
t
,
y we have the following.Wheny
t
,
y is (i)-gH differentiable(11)
y
1
′
t
,
α
=
F
t
;
y
1
t
;
α
,
y
2
t
;
α
,
y
2
′
t
,
α
=
G
t
;
y
1
t
;
α
,
y
2
t
;
α
.When y
t
,
y is (ii)-gH differentiable(12)
y
2
′
t
,
α
=
F
t
;
y
1
t
;
α
,
y
2
t
;
α
,
y
1
′
t
,
α
=
G
t
;
y
1
t
;
α
,
y
2
t
;
α
,where, by using extension principle, we have the membership function (13)
f
t
;
y
t
s
=
Sup
⁡
y
t
τ
∖
s
=
f
t
,
τ
,
s
∈
R
.So fuzzy number is f
t
;
y
t. From this it follows that (14)
f
t
;
y
t
α
=
f
1
t
,
y
t
;
α
,
f
2
t
,
y
t
;
α
,
α
∈
0
;
1
,where(15)
f
1
t
,
y
t
;
α
=
F
t
;
y
1
t
;
α
,
y
2
t
;
α
=
min
⁡
f
t
,
u
∖
u
∈
y
1
t
;
α
,
y
2
t
;
α
,
f
2
t
,
y
t
;
α
=
G
t
;
y
1
t
;
α
,
y
2
t
;
α
=
max
⁡
f
t
,
u
∖
u
∈
y
1
t
;
α
,
y
2
t
;
α
.Note. (1) Both cases ((i)-gH and (ii)-gH) can be applied to a FDE for finding exact solution.(2) After takingα-cut of the given FDE, it transforms to system of ordinary differential equation.


## 4. Numerical Solution of Fuzzy Differential Equation

### 4.1. Runge-Kutta-Fehlberg Method for Ordinary (Crisp) Differential Equation

Consider the initial value problemy
′
t
=
f
t
,
y
t
;
y
t
0
=
y
0.The Runge-Kutta-Fehlberg method (denoted as RKF45) is one way to try to resolve this problem.The problem is to solve the initial value problem in above equation by means of Runge-Kutta methods of order 4 and order 5.First we need some definitions:(16)
k
1
=
h
f
t
i
,
y
i
,
k
2
=
h
f
t
i
+
1
4
h
,
y
i
+
1
4
k
1
,
k
3
=
h
f
t
i
+
3
8
h
,
y
i
+
3
32
k
1
+
9
32
k
2
,
k
4
=
h
f
t
i
+
12
13
h
,
y
i
+
1932
2197
k
1
-
7200
2197
k
2
+
7296
2197
k
3
,
k
5
=
h
f
t
i
+
h
,
y
i
+
439
216
k
1
-
8
k
2
+
3680
513
k
3
-
845
4104
k
4
,
k
6
=
h
f
t
i
+
h
,
y
i
-
8
27
k
1
+
2
k
2
-
3544
2565
k
3
+
1859
4104
k
4
-
11
40
k
5
. Then an approximation to the solution of initial value problem is made using Runge-Kutta method of order 4:(17)
y
i
+
1
=
y
i
+
25
216
k
1
+
1408
2565
k
3
+
2197
4101
k
4
-
1
5
k
5
.A better value for the solution is determined using a Runge-Kutta method of order 5:(18)
z
i
+
1
=
y
i
+
16
135
k
1
+
6656
12
,
825
k
3
+
28
,
561
56
,
430
k
4
-
9
50
k
5
+
2
55
k
6
.The optimal step size s
h can be determined by multiplying the scalar s times the step size h. The scalar s is(19)
s
=
ϵ
h
2
z
i
+
1
-
y
i
+
1
1
/
4
=
0.0840896
ϵ
h
z
i
+
1
-
y
i
+
1
1
/
4
,where ϵ is the specified error control tolerance.Note that RK4 requires 4 function evaluations and RK5 requires 6 evaluations, that is, 10 for RK4 and RK5. Fehlberg devised a method to get RK4 and RK5 results using only 6 function evaluations by using some ofK values in both methods.


### 4.2. Runge-Kutta-Fehlberg Method for Solving Fuzzy Differential Equations

LetY
=
Y
1
,
Y
2 be the exact solution and let y
=
y
1
,
y
2 be the approximated solution of the fuzzy initial value problem.Let[
Y
(
t
)
]
α
=
[
Y
1
(
t
,
α
)
,
Y
2
(
t
,
α
)
]
,
[
y
(
t
)
]
r
=
[
y
1
(
t
,
α
)
,
y
2
(
t
,
α
)
].Throughout this argument, the value ofr is fixed. Then the exact and approximated solution at t
n are, respectively, denoted by(20)
Y
t
n
α
=
Y
1
t
n
,
α
,
Y
2
t
n
,
α
,
y
t
n
α
=
y
1
t
n
,
α
,
y
2
t
n
,
α
.The grid points at which the solution is calculated are h
=
(
T
-
t
0
)
/
N
,
t
i
=
t
0
+
i
h
,
0
≤
i
≤
N.Then we obtained(21)
y
1
t
n
+
1
,
α
=
y
1
t
n
,
α
+
16
135
K
1
+
6656
12
,
825
K
3
+
28
,
561
56
,
430
K
4
-
9
50
K
5
+
2
55
K
6
,where(22)
K
1
=
h
F
t
n
,
y
1
t
n
,
α
,
y
2
t
n
,
α
,
K
2
=
h
F
t
n
+
1
4
h
,
y
1
t
n
,
α
+
1
4
K
1
,
y
2
t
n
,
α
+
1
4
K
1
,
K
3
=
h
F
t
n
+
3
8
h
,
y
1
t
n
,
α
+
3
32
K
1
+
9
32
K
2
,
y
2
t
n
,
α
+
3
32
K
1
+
9
32
K
2
,
K
4
=
h
F
t
n
+
12
13
h
,
y
1
t
n
,
α
+
1932
2197
K
1
-
7200
2197
K
2
+
7296
2197
K
3
,
y
2
t
n
,
α
+
1932
2197
K
1
-
7200
2197
K
2
+
7296
2197
K
3
,
K
5
=
h
F
t
n
+
h
,
y
1
t
n
,
α
+
439
216
K
1
-
8
K
2
+
3680
513
K
3
-
845
4104
K
4
,
y
2
t
n
,
α
+
439
216
K
1
-
8
K
2
+
3680
513
K
3
-
845
4104
K
4
,
K
6
=
h
F
t
n
+
h
,
y
1
t
n
,
α
-
8
27
K
1
+
2
K
2
-
3544
2565
K
3
+
1859
4104
K
4
-
11
40
K
5
,
y
2
t
n
,
α
-
8
27
K
1
+
2
K
2
-
3544
2565
K
3
+
1859
4104
K
4
-
11
40
K
5
,
(23)
y
2
t
n
+
1
,
α
=
y
2
t
n
,
α
+
16
135
L
1
+
6656
12
,
825
L
3
+
28
,
561
56
,
430
L
4
-
9
50
L
5
+
2
55
L
6
,where(24)
L
1
=
h
G
t
n
,
y
1
t
n
,
α
,
y
2
t
n
,
α
,
L
2
=
h
G
t
n
+
1
4
h
,
y
1
t
n
,
α
+
1
4
L
1
,
y
2
t
n
,
α
+
1
4
L
1
,
L
3
=
h
G
t
n
+
3
8
h
,
y
1
t
n
,
α
+
3
32
L
1
+
9
32
L
2
,
y
2
t
n
,
α
+
3
32
L
1
+
9
32
L
2
,
L
4
=
h
G
t
n
+
12
13
h
,
y
1
t
n
,
α
+
1932
2197
L
1
-
7200
2197
L
2
+
7296
2197
L
3
,
y
2
t
n
,
α
+
1932
2197
L
1
-
7200
2197
L
2
+
7296
2197
L
3
,
L
5
=
h
G
t
n
+
h
,
y
1
t
n
,
α
+
439
216
L
1
-
8
L
2
+
3680
513
L
3
-
845
4104
L
4
,
y
2
t
n
,
α
+
439
216
L
1
-
8
L
2
+
3680
513
L
3
-
845
4104
L
4
,
L
6
=
h
G
t
n
+
h
,
y
1
t
n
,
α
-
8
27
L
1
+
2
L
2
-
3544
2565
L
3
+
1859
4104
L
4
-
11
40
L
5
,
y
2
t
n
,
α
-
8
27
L
1
+
2
L
2
-
3544
2565
L
3
+
1859
4104
L
4
-
11
40
L
5
.


### 4.3. Convergence of Fuzzy Runge-Kutta-Fehlberg Method

The solution is calculated by grid points ata
=
t
0
≤
t
1
≤
⋯
≤
t
N
=
b and h
=
(
b
-
a
)
/
N
=
t
n
+
1
-
t
n.Therefore, we have(25)
Y
1
t
n
+
1
,
α
=
Y
1
t
n
,
α
+
F
t
n
,
Y
1
t
n
,
α
,
Y
2
t
n
,
α
,
Y
2
t
n
+
1
,
α
=
Y
2
t
n
,
α
+
G
t
n
,
Y
1
t
n
,
α
,
Y
2
t
n
,
α
,
y
1
t
n
+
1
,
α
=
y
1
t
n
,
α
+
F
t
n
,
y
1
t
n
,
α
,
y
2
t
n
,
α
,
y
2
t
n
+
1
,
α
=
y
2
t
n
,
α
+
G
t
n
,
y
1
t
n
,
α
,
y
2
t
n
,
α
.Clearly, y
1
(
t
,
α
) and y
2
(
t
,
α
) converge to Y
1
(
t
,
α
) and Y
2
(
t
,
α
), respectively, whenever h
→
0; that is, (26)
lim
h
→
0
⁡
y
1
t
,
α
=
Y
1
t
,
α
,
lim
h
→
0
⁡
y
2
t
,
α
=
Y
2
t
,
α
.Proof.
Before we go to the main proof we need to know some results.Lemma 11.
Let the sequence of numbersW
n
=
0
N satisfy(27)
W
n
+
1
≤
A
W
n
+
B
,
0
≤
n
≤
N
-
1
,for some given positive constants A and B. Then(28)
W
n
≤
A
n
W
0
+
B
A
n
-
1
A
-
1
,
0
≤
n
≤
N
.Lemma 12.
Let the sequence of numbersW
n
=
0
N and V
n
=
0
N satisfy(29)
W
n
+
1
≤
W
n
+
A
max
⁡
W
n
,
V
n
+
B
,
V
n
+
1
≤
V
n
+
A
max
⁡
W
n
,
V
n
+
B
,for some given positive constants A and B, and denote(30)
U
n
=
W
n
+
V
n
,
0
≤
n
≤
N
.Then (31)
U
n
≤
A
-
n
U
0
+
B
-
A
-
n
-
1
A
-
-
1
,
0
≤
n
≤
N
,where A
-
=
1
+
2
A and B
-
=
2
B.LetF
(
t
,
u
,
v
) and G
(
t
,
u
,
v
) be obtained by substituting y
1
t
,
α
,
y
2
t
,
α
=
[
u
,
v
] in (21) and (23); that is,(32)
F
t
,
u
,
v
=
16
135
K
1
t
,
u
,
v
+
6656
12
,
825
K
3
t
,
u
,
v
+
28
,
561
56
,
430
K
4
t
,
u
,
v
-
9
50
K
5
t
,
u
,
v
+
2
55
K
6
t
,
u
,
v
,
G
t
,
u
,
v
=
16
135
L
1
t
,
u
,
v
+
6656
12
,
825
L
3
t
,
u
,
v
+
28
,
561
56
,
430
L
4
t
,
u
,
v
-
9
50
L
5
t
,
u
,
v
+
2
55
L
6
t
,
u
,
v
.The domain where F and G are defined is as(33)
H
=
t
,
u
,
v
∣
0
≤
t
≤
T
,
-
∞
<
v
<
∞
,
-
∞
<
u
≤
v
.Theorem 13.
LetF
t
,
u
,
v and G
t
,
u
,
v belong to C
p
-
1
(
K
) and let the partial derivative of F and G be bounded over K. Then for arbitrary fixed 0
≤
α
≤
1, the numerical solution of (9),  [
y
1
(
t
,
α
)
,
y
2
(
t
,
α
)
] converges to the exact solution [
Y
1
(
t
,
α
)
,
Y
2
(
t
,
α
)
].Proof (see [46]).
By using Taylor’s theorem we get(34)
Y
1
t
n
+
1
,
α
=
Y
1
t
n
,
α
+
h
F
t
n
,
Y
1
t
n
,
α
,
Y
2
t
n
,
α
+
h
p
+
1
p
+
1
!
Y
1
p
+
1
ξ
n
,
1
,
Y
2
t
n
+
1
,
α
=
Y
2
t
n
,
α
+
h
G
t
n
,
Y
1
t
n
,
α
,
Y
2
t
n
,
α
+
h
p
+
1
p
+
1
!
Y
2
p
+
1
ξ
n
,
2
,where ξ
n
,
1
,
ξ
n
,
2
∈
(
t
n
,
t
n
+
1
).
Now if we denote(35)
W
n
=
Y
1
t
n
,
α
-
y
1
t
n
,
α
,
V
n
=
Y
2
t
n
,
α
-
y
2
t
n
,
α
,then the above two expressions converted to(36)
W
n
+
1
=
W
n
+
h
F
t
n
,
Y
1
t
n
,
α
,
Y
2
t
n
,
α
-
F
t
n
,
y
1
t
n
,
α
,
y
2
t
n
,
α
+
h
p
+
1
p
+
1
!
Y
1
p
+
1
ξ
n
,
1
,
V
n
+
1
=
V
n
+
h
G
t
n
,
Y
1
t
n
,
α
,
Y
2
t
n
,
α
-
G
t
n
,
y
1
t
n
,
α
,
y
2
t
n
,
α
+
h
p
+
1
p
+
1
!
Y
2
p
+
1
ξ
n
,
2
.Hence we can write(37)
W
n
+
1
≤
W
n
+
2
L
h
max
⁡
W
n
,
V
n
+
h
p
+
1
p
+
1
!
M
,
V
n
+
1
≤
V
n
+
2
L
h
max
⁡
W
n
,
V
n
+
h
p
+
1
p
+
1
!
M
,where M
=
max
⁡
max
⁡
Y
1
p
+
1
(
t
;
α
)
,
max
⁡
Y
2
p
+
1
(
t
;
α
) for t
∈
[
0
,
T
] and L
>
0 is a bound for the partial derivative of F and G.
Therefore we can write(38)
W
n
≤
1
+
4
L
h
n
U
0
+
2
h
p
+
1
p
+
1
!
M
1
+
4
L
h
n
-
1
4
L
h
,
V
n
≤
1
+
4
L
h
n
U
0
+
2
h
p
+
1
p
+
1
!
M
1
+
4
L
h
n
-
1
4
L
h
,where U
0
=
W
0
+
V
0.
In particular,(39)
W
N
≤
1
+
4
L
h
N
U
0
+
2
h
p
+
1
p
+
1
!
M
1
+
4
L
h
T
/
h
-
1
4
L
h
,
V
N
≤
1
+
4
L
h
N
U
0
+
2
h
p
+
1
p
+
1
!
M
1
+
4
L
h
T
/
h
-
1
4
L
h
.Since W
0
=
V
0
=
0, we have(40)
W
N
≤
M
e
4
L
T
-
1
2
L
p
+
1
!
h
p
,
V
N
≤
M
e
4
L
T
-
1
2
L
p
+
1
!
h
p
.Thus, if h
→
0, we get W
N
→
0 and V
N
→
0, which completes the proof.


### 4.4. Algorithm for Finding the Numerical Solution

Step 1.
F
(
t
,
y
1
,
y
2
)
← “Function to be supplied”⁡
G
(
t
,
y
1
,
y
2
)
← “Function to be supplied”Step 2.
Readt
(
0
)
,
y
1
(
0
)
,
y
2
(
0
)
,
h, limit.Step 3.
Fori
=
0
(
1
) limit⁡
K
1
←
h
F
(
t
i
,
y
1
t
i
,
r
,
y
2
(
t
i
,
r
)
)
⁡
K
2
←
h
·
F
(
t
i
+
(
1
/
4
)
h
,
y
1
t
i
,
r
+
(
1
/
4
)
K
1
,
y
2
(
t
i
,
r
)
+
(
1
/
4
)
K
1
)
⁡
K
3
←
h
·
F
(
t
i
+
(
3
/
8
)
h
,
y
1
(
t
i
,
r
)
+
(
3
/
32
)
K
1
+
(
9
/
32
)
K
2
,
y
2
(
t
i
,
r
)
+
(
3
/
32
)
K
1
+
(
9
/
32
)
K
2
)
⁡
K
4
←
h
·
F
(
t
i
+
(
12
/
13
)
h
,
y
1
(
t
i
,
r
)
+
(
1932
/
2197
)
K
1
-
(
7200
/
2197
)
K
2
+
(
7296
/
2197
)
K
3, y
2
(
t
i
,
r
)  +  (
1932
/
2197
)
K
1
-
(
7200
/
2197
)
K
2
+
(
7296
/
2197
)
K
3
)
⁡
K
5
←
h
F
(
t
i
+
h
,
y
1
(
t
i
,
r
)
+
(
439
/
216
)
K
1
-
8
K
2
+
(
3680
/
513
)
K
3
-
(
845
/
4104
)
K
4
,
y
2
(
t
i
,
r
)
+
(
439
/
216
)
K
1  −  8
K
2
+
(
3680
/
513
)
K
3
-
(
845
/
4104
)
K
4
)
⁡
K
6
←
h
F
(
t
i
+
h
,
y
1
t
i
,
r
-
(
8
/
27
)
K
1
+
2
K
2
-
(
3544
/
2565
)
K
3
+
(
1859
/
4104
)
K
4
-
(
11
/
40
)
K
5
,
y
2
t
i
,
r  −  (
8
/
27
)
K
1
+
2
K
2
-
(
3544
/
2565
)
K
3
+
(
1859
/
4104
)
K
4
-
(
11
/
40
)
K
5
)
⁡
L
1
←
h
·
G
t
i
,
y
1
t
i
,
r
,
y
2
t
i
,
r
⁡
L
2
←
h
G
t
i
+
(
1
/
4
)
h
,
y
1
t
i
,
r
+
(
1
/
4
)
L
1
,
y
2
t
i
,
r
+
(
1
/
4
)
L
1
⁡
L
3
←
h
G
t
i
+
(
3
/
8
)
h
,
y
1
t
i
,
r
+
(
3
/
32
)
L
1
+
(
9
/
32
)
L
2
,
y
2
t
i
,
r
+
(
3
/
32
)
L
1
+
(
9
/
32
)
L
2
⁡
L
4
←
h
G
(
t
i
+
(
12
/
13
)
h
,
y
1
t
i
,
r
+
(
1932
/
2197
)
L
1
-
(
7200
/
2197
)
L
2
+
(
7296
/
2197
)
L
3
,
y
2
t
i
,
r  + (
1932
/
2197
)
L
1
-
(
7200
/
2197
)
L
2
+
(
7296
/
2197
)
L
3
)
⁡
L
5
←
h
G
(
t
i
+
h
,
y
1
t
i
,
r
+
(
439
/
216
)
L
1
-
8
L
2
+
(
3680
/
513
)
L
3
-
(
845
/
4104
)
L
4
,
y
2
t
i
,
r
+
(
439
/
216
)
L
1  −  8
L
2
+
(
3680
/
513
)
L
3
-
(
845
/
4104
)
L
4
)
⁡
L
6
←
h
G
(
t
i
+
h
,
y
1
t
i
,
r
-
(
8
/
27
)
L
1
+
2
L
2
-
(
3544
/
2565
)
L
3
+
(
1859
/
4104
)
L
4
-
(
11
/
40
)
L
5
,
y
2
t
i
,
r  −  (
8
/
27
)
L
1
+
2
L
2
-
(
3544
/
2565
)
L
3
+
(
1859
/
4104
)
L
4
-
(
11
/
40
)
L
5
)Step 4.
⁡
y
1
t
i
+
1
,
r
=
y
1
t
i
,
r
+
(
16
/
135
)
K
1
+
(
6656
/
12
,
825
)
K
3
+
(
28
,
561
/
56
,
430
)
K
4
-
(
9
/
50
)
K
5
+
(
2
/
55
)
K
6Step 5.
⁡
y
2
(
t
i
+
1
,
r
)
=
y
2
(
t
i
,
r
)
+
(
16
/
135
)
L
1
+
(
6656
/
12
,
825
)
L
3
+
(
28
,
561
/
56
,
430
)
L
4
-
(
9
/
50
)
L
5
+
(
2
/
55
)
L
6Step 6.
t
i
+
1
=
t
i
+
h. Write y
1
t
i
+
1
,
r
,
y
2
(
t
i
+
1
,
r
), t
i
+
1.Step 7.
NextiStep 8.
End.


## 4.1. Runge-Kutta-Fehlberg Method for Ordinary (Crisp) Differential Equation

Consider the initial value problemy
′
t
=
f
t
,
y
t
;
y
t
0
=
y
0.The Runge-Kutta-Fehlberg method (denoted as RKF45) is one way to try to resolve this problem.The problem is to solve the initial value problem in above equation by means of Runge-Kutta methods of order 4 and order 5.First we need some definitions:(16)
k
1
=
h
f
t
i
,
y
i
,
k
2
=
h
f
t
i
+
1
4
h
,
y
i
+
1
4
k
1
,
k
3
=
h
f
t
i
+
3
8
h
,
y
i
+
3
32
k
1
+
9
32
k
2
,
k
4
=
h
f
t
i
+
12
13
h
,
y
i
+
1932
2197
k
1
-
7200
2197
k
2
+
7296
2197
k
3
,
k
5
=
h
f
t
i
+
h
,
y
i
+
439
216
k
1
-
8
k
2
+
3680
513
k
3
-
845
4104
k
4
,
k
6
=
h
f
t
i
+
h
,
y
i
-
8
27
k
1
+
2
k
2
-
3544
2565
k
3
+
1859
4104
k
4
-
11
40
k
5
. Then an approximation to the solution of initial value problem is made using Runge-Kutta method of order 4:(17)
y
i
+
1
=
y
i
+
25
216
k
1
+
1408
2565
k
3
+
2197
4101
k
4
-
1
5
k
5
.A better value for the solution is determined using a Runge-Kutta method of order 5:(18)
z
i
+
1
=
y
i
+
16
135
k
1
+
6656
12
,
825
k
3
+
28
,
561
56
,
430
k
4
-
9
50
k
5
+
2
55
k
6
.The optimal step size s
h can be determined by multiplying the scalar s times the step size h. The scalar s is(19)
s
=
ϵ
h
2
z
i
+
1
-
y
i
+
1
1
/
4
=
0.0840896
ϵ
h
z
i
+
1
-
y
i
+
1
1
/
4
,where ϵ is the specified error control tolerance.Note that RK4 requires 4 function evaluations and RK5 requires 6 evaluations, that is, 10 for RK4 and RK5. Fehlberg devised a method to get RK4 and RK5 results using only 6 function evaluations by using some ofK values in both methods.


## 4.2. Runge-Kutta-Fehlberg Method for Solving Fuzzy Differential Equations

LetY
=
Y
1
,
Y
2 be the exact solution and let y
=
y
1
,
y
2 be the approximated solution of the fuzzy initial value problem.Let[
Y
(
t
)
]
α
=
[
Y
1
(
t
,
α
)
,
Y
2
(
t
,
α
)
]
,
[
y
(
t
)
]
r
=
[
y
1
(
t
,
α
)
,
y
2
(
t
,
α
)
].Throughout this argument, the value ofr is fixed. Then the exact and approximated solution at t
n are, respectively, denoted by(20)
Y
t
n
α
=
Y
1
t
n
,
α
,
Y
2
t
n
,
α
,
y
t
n
α
=
y
1
t
n
,
α
,
y
2
t
n
,
α
.The grid points at which the solution is calculated are h
=
(
T
-
t
0
)
/
N
,
t
i
=
t
0
+
i
h
,
0
≤
i
≤
N.Then we obtained(21)
y
1
t
n
+
1
,
α
=
y
1
t
n
,
α
+
16
135
K
1
+
6656
12
,
825
K
3
+
28
,
561
56
,
430
K
4
-
9
50
K
5
+
2
55
K
6
,where(22)
K
1
=
h
F
t
n
,
y
1
t
n
,
α
,
y
2
t
n
,
α
,
K
2
=
h
F
t
n
+
1
4
h
,
y
1
t
n
,
α
+
1
4
K
1
,
y
2
t
n
,
α
+
1
4
K
1
,
K
3
=
h
F
t
n
+
3
8
h
,
y
1
t
n
,
α
+
3
32
K
1
+
9
32
K
2
,
y
2
t
n
,
α
+
3
32
K
1
+
9
32
K
2
,
K
4
=
h
F
t
n
+
12
13
h
,
y
1
t
n
,
α
+
1932
2197
K
1
-
7200
2197
K
2
+
7296
2197
K
3
,
y
2
t
n
,
α
+
1932
2197
K
1
-
7200
2197
K
2
+
7296
2197
K
3
,
K
5
=
h
F
t
n
+
h
,
y
1
t
n
,
α
+
439
216
K
1
-
8
K
2
+
3680
513
K
3
-
845
4104
K
4
,
y
2
t
n
,
α
+
439
216
K
1
-
8
K
2
+
3680
513
K
3
-
845
4104
K
4
,
K
6
=
h
F
t
n
+
h
,
y
1
t
n
,
α
-
8
27
K
1
+
2
K
2
-
3544
2565
K
3
+
1859
4104
K
4
-
11
40
K
5
,
y
2
t
n
,
α
-
8
27
K
1
+
2
K
2
-
3544
2565
K
3
+
1859
4104
K
4
-
11
40
K
5
,
(23)
y
2
t
n
+
1
,
α
=
y
2
t
n
,
α
+
16
135
L
1
+
6656
12
,
825
L
3
+
28
,
561
56
,
430
L
4
-
9
50
L
5
+
2
55
L
6
,where(24)
L
1
=
h
G
t
n
,
y
1
t
n
,
α
,
y
2
t
n
,
α
,
L
2
=
h
G
t
n
+
1
4
h
,
y
1
t
n
,
α
+
1
4
L
1
,
y
2
t
n
,
α
+
1
4
L
1
,
L
3
=
h
G
t
n
+
3
8
h
,
y
1
t
n
,
α
+
3
32
L
1
+
9
32
L
2
,
y
2
t
n
,
α
+
3
32
L
1
+
9
32
L
2
,
L
4
=
h
G
t
n
+
12
13
h
,
y
1
t
n
,
α
+
1932
2197
L
1
-
7200
2197
L
2
+
7296
2197
L
3
,
y
2
t
n
,
α
+
1932
2197
L
1
-
7200
2197
L
2
+
7296
2197
L
3
,
L
5
=
h
G
t
n
+
h
,
y
1
t
n
,
α
+
439
216
L
1
-
8
L
2
+
3680
513
L
3
-
845
4104
L
4
,
y
2
t
n
,
α
+
439
216
L
1
-
8
L
2
+
3680
513
L
3
-
845
4104
L
4
,
L
6
=
h
G
t
n
+
h
,
y
1
t
n
,
α
-
8
27
L
1
+
2
L
2
-
3544
2565
L
3
+
1859
4104
L
4
-
11
40
L
5
,
y
2
t
n
,
α
-
8
27
L
1
+
2
L
2
-
3544
2565
L
3
+
1859
4104
L
4
-
11
40
L
5
.


## 4.3. Convergence of Fuzzy Runge-Kutta-Fehlberg Method

The solution is calculated by grid points ata
=
t
0
≤
t
1
≤
⋯
≤
t
N
=
b and h
=
(
b
-
a
)
/
N
=
t
n
+
1
-
t
n.Therefore, we have(25)
Y
1
t
n
+
1
,
α
=
Y
1
t
n
,
α
+
F
t
n
,
Y
1
t
n
,
α
,
Y
2
t
n
,
α
,
Y
2
t
n
+
1
,
α
=
Y
2
t
n
,
α
+
G
t
n
,
Y
1
t
n
,
α
,
Y
2
t
n
,
α
,
y
1
t
n
+
1
,
α
=
y
1
t
n
,
α
+
F
t
n
,
y
1
t
n
,
α
,
y
2
t
n
,
α
,
y
2
t
n
+
1
,
α
=
y
2
t
n
,
α
+
G
t
n
,
y
1
t
n
,
α
,
y
2
t
n
,
α
.Clearly, y
1
(
t
,
α
) and y
2
(
t
,
α
) converge to Y
1
(
t
,
α
) and Y
2
(
t
,
α
), respectively, whenever h
→
0; that is, (26)
lim
h
→
0
⁡
y
1
t
,
α
=
Y
1
t
,
α
,
lim
h
→
0
⁡
y
2
t
,
α
=
Y
2
t
,
α
.Proof.
Before we go to the main proof we need to know some results.Lemma 11.
Let the sequence of numbersW
n
=
0
N satisfy(27)
W
n
+
1
≤
A
W
n
+
B
,
0
≤
n
≤
N
-
1
,for some given positive constants A and B. Then(28)
W
n
≤
A
n
W
0
+
B
A
n
-
1
A
-
1
,
0
≤
n
≤
N
.Lemma 12.
Let the sequence of numbersW
n
=
0
N and V
n
=
0
N satisfy(29)
W
n
+
1
≤
W
n
+
A
max
⁡
W
n
,
V
n
+
B
,
V
n
+
1
≤
V
n
+
A
max
⁡
W
n
,
V
n
+
B
,for some given positive constants A and B, and denote(30)
U
n
=
W
n
+
V
n
,
0
≤
n
≤
N
.Then (31)
U
n
≤
A
-
n
U
0
+
B
-
A
-
n
-
1
A
-
-
1
,
0
≤
n
≤
N
,where A
-
=
1
+
2
A and B
-
=
2
B.LetF
(
t
,
u
,
v
) and G
(
t
,
u
,
v
) be obtained by substituting y
1
t
,
α
,
y
2
t
,
α
=
[
u
,
v
] in (21) and (23); that is,(32)
F
t
,
u
,
v
=
16
135
K
1
t
,
u
,
v
+
6656
12
,
825
K
3
t
,
u
,
v
+
28
,
561
56
,
430
K
4
t
,
u
,
v
-
9
50
K
5
t
,
u
,
v
+
2
55
K
6
t
,
u
,
v
,
G
t
,
u
,
v
=
16
135
L
1
t
,
u
,
v
+
6656
12
,
825
L
3
t
,
u
,
v
+
28
,
561
56
,
430
L
4
t
,
u
,
v
-
9
50
L
5
t
,
u
,
v
+
2
55
L
6
t
,
u
,
v
.The domain where F and G are defined is as(33)
H
=
t
,
u
,
v
∣
0
≤
t
≤
T
,
-
∞
<
v
<
∞
,
-
∞
<
u
≤
v
.Theorem 13.
LetF
t
,
u
,
v and G
t
,
u
,
v belong to C
p
-
1
(
K
) and let the partial derivative of F and G be bounded over K. Then for arbitrary fixed 0
≤
α
≤
1, the numerical solution of (9),  [
y
1
(
t
,
α
)
,
y
2
(
t
,
α
)
] converges to the exact solution [
Y
1
(
t
,
α
)
,
Y
2
(
t
,
α
)
].Proof (see [46]).
By using Taylor’s theorem we get(34)
Y
1
t
n
+
1
,
α
=
Y
1
t
n
,
α
+
h
F
t
n
,
Y
1
t
n
,
α
,
Y
2
t
n
,
α
+
h
p
+
1
p
+
1
!
Y
1
p
+
1
ξ
n
,
1
,
Y
2
t
n
+
1
,
α
=
Y
2
t
n
,
α
+
h
G
t
n
,
Y
1
t
n
,
α
,
Y
2
t
n
,
α
+
h
p
+
1
p
+
1
!
Y
2
p
+
1
ξ
n
,
2
,where ξ
n
,
1
,
ξ
n
,
2
∈
(
t
n
,
t
n
+
1
).
Now if we denote(35)
W
n
=
Y
1
t
n
,
α
-
y
1
t
n
,
α
,
V
n
=
Y
2
t
n
,
α
-
y
2
t
n
,
α
,then the above two expressions converted to(36)
W
n
+
1
=
W
n
+
h
F
t
n
,
Y
1
t
n
,
α
,
Y
2
t
n
,
α
-
F
t
n
,
y
1
t
n
,
α
,
y
2
t
n
,
α
+
h
p
+
1
p
+
1
!
Y
1
p
+
1
ξ
n
,
1
,
V
n
+
1
=
V
n
+
h
G
t
n
,
Y
1
t
n
,
α
,
Y
2
t
n
,
α
-
G
t
n
,
y
1
t
n
,
α
,
y
2
t
n
,
α
+
h
p
+
1
p
+
1
!
Y
2
p
+
1
ξ
n
,
2
.Hence we can write(37)
W
n
+
1
≤
W
n
+
2
L
h
max
⁡
W
n
,
V
n
+
h
p
+
1
p
+
1
!
M
,
V
n
+
1
≤
V
n
+
2
L
h
max
⁡
W
n
,
V
n
+
h
p
+
1
p
+
1
!
M
,where M
=
max
⁡
max
⁡
Y
1
p
+
1
(
t
;
α
)
,
max
⁡
Y
2
p
+
1
(
t
;
α
) for t
∈
[
0
,
T
] and L
>
0 is a bound for the partial derivative of F and G.
Therefore we can write(38)
W
n
≤
1
+
4
L
h
n
U
0
+
2
h
p
+
1
p
+
1
!
M
1
+
4
L
h
n
-
1
4
L
h
,
V
n
≤
1
+
4
L
h
n
U
0
+
2
h
p
+
1
p
+
1
!
M
1
+
4
L
h
n
-
1
4
L
h
,where U
0
=
W
0
+
V
0.
In particular,(39)
W
N
≤
1
+
4
L
h
N
U
0
+
2
h
p
+
1
p
+
1
!
M
1
+
4
L
h
T
/
h
-
1
4
L
h
,
V
N
≤
1
+
4
L
h
N
U
0
+
2
h
p
+
1
p
+
1
!
M
1
+
4
L
h
T
/
h
-
1
4
L
h
.Since W
0
=
V
0
=
0, we have(40)
W
N
≤
M
e
4
L
T
-
1
2
L
p
+
1
!
h
p
,
V
N
≤
M
e
4
L
T
-
1
2
L
p
+
1
!
h
p
.Thus, if h
→
0, we get W
N
→
0 and V
N
→
0, which completes the proof.


## 4.4. Algorithm for Finding the Numerical Solution

Step 1.
F
(
t
,
y
1
,
y
2
)
← “Function to be supplied”⁡
G
(
t
,
y
1
,
y
2
)
← “Function to be supplied”Step 2.
Readt
(
0
)
,
y
1
(
0
)
,
y
2
(
0
)
,
h, limit.Step 3.
Fori
=
0
(
1
) limit⁡
K
1
←
h
F
(
t
i
,
y
1
t
i
,
r
,
y
2
(
t
i
,
r
)
)
⁡
K
2
←
h
·
F
(
t
i
+
(
1
/
4
)
h
,
y
1
t
i
,
r
+
(
1
/
4
)
K
1
,
y
2
(
t
i
,
r
)
+
(
1
/
4
)
K
1
)
⁡
K
3
←
h
·
F
(
t
i
+
(
3
/
8
)
h
,
y
1
(
t
i
,
r
)
+
(
3
/
32
)
K
1
+
(
9
/
32
)
K
2
,
y
2
(
t
i
,
r
)
+
(
3
/
32
)
K
1
+
(
9
/
32
)
K
2
)
⁡
K
4
←
h
·
F
(
t
i
+
(
12
/
13
)
h
,
y
1
(
t
i
,
r
)
+
(
1932
/
2197
)
K
1
-
(
7200
/
2197
)
K
2
+
(
7296
/
2197
)
K
3, y
2
(
t
i
,
r
)  +  (
1932
/
2197
)
K
1
-
(
7200
/
2197
)
K
2
+
(
7296
/
2197
)
K
3
)
⁡
K
5
←
h
F
(
t
i
+
h
,
y
1
(
t
i
,
r
)
+
(
439
/
216
)
K
1
-
8
K
2
+
(
3680
/
513
)
K
3
-
(
845
/
4104
)
K
4
,
y
2
(
t
i
,
r
)
+
(
439
/
216
)
K
1  −  8
K
2
+
(
3680
/
513
)
K
3
-
(
845
/
4104
)
K
4
)
⁡
K
6
←
h
F
(
t
i
+
h
,
y
1
t
i
,
r
-
(
8
/
27
)
K
1
+
2
K
2
-
(
3544
/
2565
)
K
3
+
(
1859
/
4104
)
K
4
-
(
11
/
40
)
K
5
,
y
2
t
i
,
r  −  (
8
/
27
)
K
1
+
2
K
2
-
(
3544
/
2565
)
K
3
+
(
1859
/
4104
)
K
4
-
(
11
/
40
)
K
5
)
⁡
L
1
←
h
·
G
t
i
,
y
1
t
i
,
r
,
y
2
t
i
,
r
⁡
L
2
←
h
G
t
i
+
(
1
/
4
)
h
,
y
1
t
i
,
r
+
(
1
/
4
)
L
1
,
y
2
t
i
,
r
+
(
1
/
4
)
L
1
⁡
L
3
←
h
G
t
i
+
(
3
/
8
)
h
,
y
1
t
i
,
r
+
(
3
/
32
)
L
1
+
(
9
/
32
)
L
2
,
y
2
t
i
,
r
+
(
3
/
32
)
L
1
+
(
9
/
32
)
L
2
⁡
L
4
←
h
G
(
t
i
+
(
12
/
13
)
h
,
y
1
t
i
,
r
+
(
1932
/
2197
)
L
1
-
(
7200
/
2197
)
L
2
+
(
7296
/
2197
)
L
3
,
y
2
t
i
,
r  + (
1932
/
2197
)
L
1
-
(
7200
/
2197
)
L
2
+
(
7296
/
2197
)
L
3
)
⁡
L
5
←
h
G
(
t
i
+
h
,
y
1
t
i
,
r
+
(
439
/
216
)
L
1
-
8
L
2
+
(
3680
/
513
)
L
3
-
(
845
/
4104
)
L
4
,
y
2
t
i
,
r
+
(
439
/
216
)
L
1  −  8
L
2
+
(
3680
/
513
)
L
3
-
(
845
/
4104
)
L
4
)
⁡
L
6
←
h
G
(
t
i
+
h
,
y
1
t
i
,
r
-
(
8
/
27
)
L
1
+
2
L
2
-
(
3544
/
2565
)
L
3
+
(
1859
/
4104
)
L
4
-
(
11
/
40
)
L
5
,
y
2
t
i
,
r  −  (
8
/
27
)
L
1
+
2
L
2
-
(
3544
/
2565
)
L
3
+
(
1859
/
4104
)
L
4
-
(
11
/
40
)
L
5
)Step 4.
⁡
y
1
t
i
+
1
,
r
=
y
1
t
i
,
r
+
(
16
/
135
)
K
1
+
(
6656
/
12
,
825
)
K
3
+
(
28
,
561
/
56
,
430
)
K
4
-
(
9
/
50
)
K
5
+
(
2
/
55
)
K
6Step 5.
⁡
y
2
(
t
i
+
1
,
r
)
=
y
2
(
t
i
,
r
)
+
(
16
/
135
)
L
1
+
(
6656
/
12
,
825
)
L
3
+
(
28
,
561
/
56
,
430
)
L
4
-
(
9
/
50
)
L
5
+
(
2
/
55
)
L
6Step 6.
t
i
+
1
=
t
i
+
h. Write y
1
t
i
+
1
,
r
,
y
2
(
t
i
+
1
,
r
), t
i
+
1.Step 7.
NextiStep 8.
End.


## 5. Numerical Example

Example. Solve y
′
=
-
y
+
t
+
1 with initial condition y
0
=
(
0.96,1
,
1.01
). Then find the solution at t
=
0.1.Solution. For (i)-gH differentiable case the exact solution is(41)
y
1
r
,
t
=
t
+
0.96
+
0.04
r
e
-
t
,
y
2
r
,
t
=
t
+
1.01
-
0.01
r
e
-
tand for (ii)-gH differentiable case the exact solution is(42)
y
1
r
,
t
=
1
+
t
+
-
0.04
+
0.04
r
e
t
,
y
2
r
,
t
=
1
+
t
+
0.01
-
0.01
r
e
t
.Remark 14.
From Figure1 and Table 1 we conclude that the lower exact solution ((i)-gH case) is approximately equal to the numerical solution when we take the step length h
=
0.01 (for h
=
0.001 is nearly equal), whereas the lower exact solution ((ii)-gH case) is approximately equal to the numerical solution when we take the step length h
=
0.1.Table 1
Comparison of the exact solutions and numerical solutions for different step lengths att
=
0.1.
r
Exact solution for (i)-gH differentiable case
Exact solution for (ii)-gH differentiable case
Numerical solution forh
=
0.1 by RKF method
⁢Numerical solution forh
=
0.01 by RKF method
⁢Numerical solution forh
=
0.001 by RKF method
Y
1
Y
2
Y
1
Y
2
y
1
y
2
y
1
y
2
y
1
y
2
0
0.9686
1.0139
1.0558
1.1111
1.0597
1.1149
0.9696
1.0201
0.9610
1.0110
0.1
0.9723
1.0130
1.0602
1.1099
1.0642
1.1138
0.9737
1.0191
0.9650
1.0100
0.2
0.9759
1.0121
1.0646
1.1088
1.0686
1.1127
0.9777
1.0181
0.9690
1.0090
0.3
0.9795
1.0112
1.0691
1.1077
1.0730
1.1116
0.9818
1.0171
0.9730
1.0080
0.4
0.9831
1.0103
1.0735
1.1066
1.0774
1.1105
0.9858
1.0161
0.9770
1.0070
0.5
0.9867
1.0094
1.0779
1.1055
1.0818
1.1094
0.9898
1.0151
0.9810
1.0060
0.6
0.9904
1.0085
1.0823
1.1044
1.0862
1.1083
0.9939
1.0141
0.9850
1.0050
0.7
0.9940
1.0076
1.0867
1.1033
1.0906
1.1072
0.9979
1.0131
0.9890
1.0040
0.8
0.9976
1.0066
1.0912
1.1022
1.0951
1.1061
1.0020
1.0121
0.9930
1.0030
0.9
1.0012
1.0057
1.0956
1.1011
1.0995
1.1050
1.0060
1.0110
0.9970
1.0020
1
1.0048
1.0048
1.1000
1.1000
1.1039
1.1039
1.0100
1.0100
1.0010
1.0010Figure 1
Comparison of lower exact solutions and numerical solution for different step lengths att
=
0.1.Remark 15.
From Figure2 and Table 1 we conclude that the upper exact solution ((i)-gH case) is approximately equal to the numerical solution when we take the step length h
=
0.001 (for h
=
0.01 is nearly equal), whereas the upper exact solution ((ii)-gH case) is approximately equal to the numerical solution when we take the step length h
=
0.1.Figure 2
Comparison of upper exact solutions and numerical solution for different step lengths att
=
0.1.


## 6. Application

Problem. A tank initially contains 300 gals of brine which has dissolved in c lbs of salt. Coming into the tank at 3 gals/min is brine with concentration k lbs salt/gals and the well stirred mixture leaves at the rate 3 gals/min. Let y
x lbs be the salt in the tank at any time t
≥
0. Then d
y
(
x
)
/
d
x
+
(
1
/
100
)
y
x
=
k, x
∈
[
0,0.5
] with y
0
=
c, if the initial condition is being modeled as fuzzy numbers c
=
(
1,2
,
3
) and k
=
(
1,2
,
4
). Find solution at x
=
0.4.Solution. For (i)-gH differentiable case the exact solution is(43)
y
1
x
,
α
=
1
+
α
100
1
+
99
e
-
1
/
100
x
,
y
2
x
,
α
=
2
-
α
50
+
148
-
49
α
50
e
-
1
/
100
x
.For (ii)-gH differentiable case the exact solution is(44)
y
1
x
,
α
=
149
-
149
α
e
1
/
100
x
+
-
248
+
50
α
e
-
1
/
100
x
+
100
+
100
α
,
y
2
x
,
α
=
-
149
-
149
α
e
1
/
100
x
+
-
248
+
50
α
e
-
1
/
100
x
+
400
-
200
α
.Remark 16.
From Figure3 and Table 2 we conclude that the lower exact solution ((i)-gH case) is approximately equal to the numerical solution when we take the step length h
=
0.001 (for h
=
0.01 is nearly equal), whereas the lower exact solution ((ii)-gH case) is not equal to any numerical solution.Table 2
Comparison of the exact solutions and numerical solutions for different step lengths atx
=
0.4.
α
Exact solution for (i)-gH differentiable case
Exact solution for (ii)-gH differentiable case
Numerical solution forh
=
0.1 by RKF method
⁢Numerical solution forh
=
0.01 by RKF method
Numerical solution forh
=
0.001 by RKF method
Y
1
Y
2
Y
1
Y
2
y
1
y
2
y
1
y
2
y
1
y
2
0
0.9960
2.9882
2.5872
3.3928
1.1039
3.3117
1.0100
3.0301
1.0010
3.0030
0.1
1.0957
2.8886
2.6075
3.3326
1.2143
3.2013
1.1110
2.9291
1.1011
2.9029
0.2
1.1953
2.7890
2.6279
3.2723
1.3247
3.0909
1.2120
2.8281
1.2012
2.8028
0.3
1.2949
2.6894
2.6482
3.2121
1.4351
2.9805
1.3130
2.7271
1.3013
2.7027
0.4
1.3945
2.5897
2.6685
3.1519
1.5455
2.8701
1.4141
2.6261
1.4014
2.6026
0.5
1.4941
2.4901
2.6888
3.0916
1.6558
2.7597
1.5151
2.5251
1.5015
2.5025
0.6
1.5937
2.3905
2.7091
3.0314
1.7662
2.6493
1.6161
2.4241
1.6016
2.4024
0.7
1.6933
2.2909
2.7295
2.9711
1.8766
2.5390
1.7171
2.3231
1.7017
2.3023
0.8
1.7929
2.1913
2.7498
2.9109
1.9870
2.4286
1.8181
2.2221
1.8018
2.2022
0.9
1.8925
2.0917
2.7701
2.8507
2.0974
2.3182
1.9191
2.1211
1.9019
2.1021
1
1.9921
1.9921
2.7904
2.7904
2.2078
2.2078
2.0201
2.0201
2.0020
2.0020Figure 3
Comparison of lower exact solutions and numerical solution for different step lengths att
=
0.4.Remark 17.
From Figure4 and Table 2 we conclude that the upper exact solution ((i)-gH case) is approximately equal to the numerical solution when we take the step length of h
=
0.01 and h
=
0.001. For h
=
0.1 it is nearly equal, whereas the upper exact solution ((ii)-gH case) is not equal to any numerical solution.Figure 4
Comparison of upper exact solutions and numerical solution for different step lengths att
=
0.4.


## 7. Conclusion

In this paper we applied Runge-Kutta-Fehlberg method for finding the numerical solution of first-order linear differential equation in fuzzy environment. The numerical solution is compared with the exact solution ((i)-gH and (ii)-gH both cases). The results presented in the contribution show that Runge-Kutta-Fehlberg method is a powerful mathematical tool for solving first-order linear differential equation in fuzzy environment. The convergence of Runge-Kutta-Fehlberg method has been discussed. The process method is applied to a mechanical problem in fuzzy environment which shows that it is a promising method to solve the said types of differential equation. In the future we can apply these methods for solving higher order linear and nonlinear differential equation in fuzzy environment.



---

*Source: 8150497-2016-05-12.xml*