Diagonally Implicit Multistep Block Method of Order Four for Solving Fuzzy Differential Equations Using Seikkala Derivatives

: In this paper, the solution of fuzzy differential equations is approximated numerically using diagonally implicit multistep block method of order four. The multistep block method is well known as an efﬁcient and accurate method for solving ordinary differential equations, hence in this paper the method will be used to solve the fuzzy initial value problems where the initial value is a symmetric triangular fuzzy interval. The triangular fuzzy number is not necessarily symmetric, however by imposing symmetry the deﬁnition of a triangular fuzzy number can be simpliﬁed. The symmetric triangular fuzzy interval is a triangular fuzzy interval that has same left and right width of membership function from the center. Due to this, the parametric form of symmetric triangular fuzzy number is simple and the performing arithmetic operations become easier. In order to interpret the fuzzy problems, Seikkala’s derivative approach is implemented. Characterization theorem is then used to translate the problems into a system of ordinary differential equations. The convergence of the introduced method is also proved. Numerical examples are given to investigate the performance of the proposed method. It is clearly shown in the results that the proposed method is comparable and reliable in solving fuzzy differential equations.


Introduction
Modeling of real-world problems involves ordinary differential equations (ODEs) that are not always perfect.For example, the initial value may not be known exactly and the function may contain uncertain parameters.This inexactness leads to the necessity of fuzzy differential equations (FDEs) to overcome the situation.Examples of FDEs application in real life include models such as in solid waste management systems [1], hydraulic differential servo cylinders [2], population models [3], hyperchaotic systems [4], and service composition [5].
Chang and Zadeh in [6] first introduced the fuzzy derivative concept.The concept was then extended by Dubois and Prade [7] and the authors also addressed the problem of linearity.Goetschel and Voxman [8] and Puri and Ralescu [9] studied and proposed several important definitions and theorems in the fuzzy derivative.A significant contribution for fuzzy initial value problems (FIVPs) is done by Seikkala [10] and Kaleva [11,12].
The exact solutions for certain FDEs are inconvenient and tedious to determine because of the complexity of the FDE problems.Therefore, numerical methods are used for solving the FDEs.

Preliminaries
Some basic definitions and notations that are significant for this paper will be reviewed in this section.For further details, refer [15,[24][25][26] (to name a few).Definition 1.Let η : R → [0, 1] is a mapping of fuzzy number with R is set of all real numbers.The mapping has following properties: 1.
Consider E as set of all fuzzy numbers on R.
Definition 2. The r-level set of a fuzzy number η ∈ E, 0 ≤ r ≤ 1 and denoted by [η] r , is defined as and is closed and bounded interval by η(r), η(r) where η(r) and η(r) denotes the left-hand endpoint and right-hand endpoint of [η] r respectively.
Given the Hausdorff distance D : E × E → R + ∪ 0 as the distance between two fuzzy numbers and a metric in E as follows (E, D) is a complete metric space.
Let f : R → E be a fuzzy valued function.
Hukuhara differentiability (H-derivative) for fuzzy mappings is introduced by Puri and Ralescu in [9].The differentiability is based on H-difference sets.Definitions 4 and 5 explained the concept.Definition 4. Let x, y ∈ E. If there exist z ∈ E such that x = y ⊕ z, then z is called the H-difference of x and y, and it is denoted by x y.The sign " " stands for H-difference and note that x y = x + (−1)y.
exist and are equal to f (t 0 ).
Here, the limits are taken in the metric space (E, D) since we have defined

Fuzzy Initial Value Problems
In this paper, we will consider first order FIVP as follows where y is a fuzzy function of t, f (t, y) is a fuzzy function of crisp variable t and fuzzy variable y, and y is Hukuhara fuzzy derivative of y with the given initial value y(t 0 ) = y 0 .Let y(t) = y(t), y(t) , if y(t) is Hukuhara differentiable then y (t) = y (t), y (t) , and f (t, y(t)) = f t, y(t), y(t) , f t, y(t), y(t) .The parametric form of y is given by As referred to [21], the Seikkala derivative of a fuzzy function is defined by provided that the equation defines a fuzzy number.
For solving (1) by using numerical methods for ODEs, it must be translated to system of ODEs.Bede in [27] proposed characterization theorems for the solutions of FDEs by using ODEs.
there exist an L > 0 such that Then FIVP (1) and the system of ODEs (2) are equivalent.

Derivation of 2-Point Diagonally Implicit Multistep Block Method
Consider the initial value problems (IVPs) for ODEs of the form where a and b are finite.The interval [a, b] will be divided into a series of blocks.In this proposed method, each block contained two points, y n+1 and y n+2 which will be computed simultaneously in a block with the step size h at the points x n+1 and x n+2 respectively as in Figure 1.
Consider the initial value problems (IVPs) for ODEs of the form where and are finite.The interval , will be divided into a series of blocks.In this proposed method, each block contained two points, and which will be computed simultaneously in a block with the step size ℎ at the points and respectively as in Figure 1.The Lagrange interpolation polynomial is applied for the derivation of predictor and corrector formulas.To obtain the predictor formulas, the set of points ( , ), ( , ), ( , ) is interpolated and the order of the formulas is one less than the corrector formulas.The corrector formulas will use the set of interpolation points i.e., ( , ), ( , ), ( , ), ( , ) for derivation of and ( , ), ( , ), ( , ), ( , ) for derivation of .The The Lagrange interpolation polynomial is applied for the derivation of predictor and corrector formulas.To obtain the predictor formulas, the set of points {(x n−2 , f n−2 ), (x n−1 , f n−1 ), (x n , f n )} is interpolated and the order of the formulas is one less than the corrector formulas.The corrector formulas will use the set of interpolation points i.e., {( The formulas of y n+1 will used three back values (x n , x n−1 , x n−2 ), while the y n+2 will used two back values i.e., (x n , x n−1 ).
Formulas for y n+1 and y n+2 are derived by integrating (3) such that where the interval of integration for y n+1 and y n+2 are [x n , x n+1 ] and [x n , x n+2 ] respectively.These are equivalent to and The function f (x, y) in ( 5) and ( 6) is replaced by the Lagrange polynomial which interpolates the set of corresponding mentioned points.Evaluating the integral in (5) using MAPLE gives the formula for the first and second point of the block multistep method as follows: The first point, and the second point, As referred to [28,29], form the general matrix of linear multistep method where A k and B k are r by r matrices with elements a lm and b lm for l, m = 1, 2, ..., r.The associated difference operator L is given by where z(x) is the exact solution and assume to be sufficiently differentiable.
L is said to be order p and C p+1 is called the error constant.
Symmetry 2018, 10, 42 The general form of constant C q is defined as To determine the order of the method, ( 7) and ( 8) is written in matrix form as follows 1 0 0 1 From ( 12), , B 4 = and . Hence, by applying (11),

Implementation for Solving FIVPs
Consider the FIVP (1) and let the exact solution is denoted by Y(t; r), Y(t; r) while y(t; r), y(t; r) is the approximate solution.To integrate the system (2), from t 0 let T > t 0 , and replace [t 0 , T] by Y n (t; r) = Y n (t; r), Y n (t; r) denoted as the exact solution at t n , with the approximate solution y n (t; r) = y n (t; r), y n (t; r) .Given t n = t 0 + h where 0 ≤ n ≤ N and h = T−t 0 N .The general fuzzy 2PDO4 method by considering (7) and ( 8) is written as for the exact solutions and y(t n+1 ; r) = y(t n ; r) + h 24 9 f t n+1 , y n+1 (r), 16) for the approximate solutions.

Algorithm
The following algorithm is based on using 2PDO4 method.To approximate the solution of the following FIVP For N as an arbitrary positive integer,

Convergence
From ( 16) and ( 17) is the convergence condition for the approximates.The following lemma is considered to proof the convergence.
Lemma 1. [9] Let {w n } N n=0 (a sequence of numbers) satisfy: for some given positive constants A, B and C. Then

Results and Discussion
To implement the 2PDO4 method, PE(CE) m mode is used where P, C and E stands for predictor, corrector and evaluation of a function respectively.The CE implementation is repeated m times for each step and step size h = 0.1 is considered in the numerical results.
The following FIVPs were tested to show the performance of the proposed method: where triangular fuzzy number c = (1/2/3).
[c is a triangular fuzzy number where the basis of the membership function of c is in the interval [c 1 , c 3 ] with the known summit c 2 .Hence, the triangular fuzzy numbers are denoted by c = (c 1 /c 2 /c 3 )].
The numerical results of 2PDO4 compared to ERK4, RK4 and ANN are display in Tables 1-5 when solving Problem 1-5.Note that these results are at t = 1 with h = 0.1.
From Tables 1-3, it is obvious that 2PDO4 gives better results compared to ERK4, RK4 and ANN in terms of accuracy.It is also observed that the execution times of 2PDO4 for solving Problem 3 and 5 are faster than the RK4.Table 4 shows that 2PDO4 used less number of function evaluations to obtain comparable accuracy compared to ERK4 and RK4.
For Problem 5, the exact solutions cannot be calculated analytically.Hence, we compared the approximate solutions of 2PDO4 and RK4.The approximate solutions of Problem 5 at t = 1 with h = 0.1 are presented in Table 5.It can be observed that 2PDO4 is applicable to solve Problem 5 with an advantage of less total steps and faster execution times compared to RK4.
It is obvious in 1−5 that 2PDO4 required less number of total steps compared to other methods.This is expected since 2PDO4 approximates the solutions at two points simultaneously.It is clear in the figures that the solutions by the 2PDO4 method almost coincide with the exact solutions.While the plot of approximate solutions for Problem 5 by using 2PDO4 and RK4 is shown in Figure 6.It is clear in the figures that the solutions by the 2PDO4 method almost coincide with the exact solutions.While the plot of approximate solutions for Problem 5 by using 2PDO4 and RK4 is shown in Figure 6.

Conclusions
In this paper, the diagonally implicit block method of order four (2PDO4) has been implemented for solving linear and nonlinear FDEs.Numerical results have shown that the proposed method gave comparable results in terms of accuracy and fewer total steps.The total functions are comparable or less compared to the existing method.The 2PDO4 is able to obtain less execution time than RK4 when solving two tested problems.Hence, the 2PDO4 method can be named in the literature as an

Conclusions
In this paper, the diagonally implicit block method of order four (2PDO4) has been implemented for solving linear and nonlinear FDEs.Numerical results have shown that the proposed method gave comparable results in terms of accuracy and fewer total steps.The total functions are comparable or less compared to the existing method.The 2PDO4 is able to obtain less execution time than RK4 when

Conclusions
In this paper, the diagonally implicit block method of order four (2PDO4) has been implemented for solving linear and nonlinear FDEs.Numerical results have shown that the proposed method gave comparable results in terms of accuracy and fewer total steps.The total functions are comparable or less compared to the existing method.The 2PDO4 is able to obtain less execution time than RK4 when solving two tested problems.Hence, the 2PDO4 method can be named in the literature as an additional and reliable method for solving FDEs.

Problem 5 .
(Ma et al. [13]) Consider nonlinear FIVP:y (t) = e −y 2 (t) , y(0) = [0.75+ 0.25r, 1.5 − 0.5r], t, r ∈ [0, 1].The following notations are used in the tables: r: Fuzzy numbers with bounded r-level intervals Y, Y: Lower and upper bounded exact solution y, y: Lower and upper bounded approximated solution FCN: Total function calls TS: Total steps taken TIME: The execution time taken in seconds 2PDO4: 2-point diagonally implicit multistep method of order four ERK4: Extended Runge-Kutta-like formulae of order 4[26] results of 2PDO4 and RK4 have been obtained in the same platform.

Figures 2 -
Figures 2-5 is a plot of the exact and approximate solutions by using 2PDO4 for Problems 1-4.It is clear in the figures that the solutions by the 2PDO4 method almost coincide with the exact solutions.While the plot of approximate solutions for Problem 5 by using 2PDO4 and RK4 is shown in Figure6.

Figures 2 -
Figures 2-5 is a plot of the exact and approximate solutions by using 2PDO4 for Problems 1-4.It is clear in the figures that the solutions by the 2PDO4 method almost coincide with the exact solutions.While the plot of approximate solutions for Problem 5 by using 2PDO4 and RK4 is shown in Figure6.

Figure 2 .
Figure 2. Exact solutions and approximate solutions for Problem 1 at t = 1 with h = 0.1.

Figures 2 -
Figures 2-5 is a plot of the exact and approximate solutions by using 2PDO4 for Problems 1-4.It is clear in the figures that the solutions by the 2PDO4 method almost coincide with the exact solutions.While the plot of approximate solutions for Problem 5 by using 2PDO4 and RK4 is shown in Figure6.

Figure 3 .
Figure 3. Exact solutions and approximate solutions for Problem 2 at = 1 with ℎ = 0.1.Figure 3. Exact solutions and approximate solutions for Problem 2 at t = 1 with h = 0.1.

Figure 4 .
Figure 4. Exact solutions and approximate solutions for Problem 3 at = 1 with ℎ = 0.1.Figure 4. Exact solutions and approximate solutions for Problem 3 at t = 1 with h = 0.1.

Figure 4 .
Figure 4. Exact solutions and approximate solutions for Problem 3 at = 1 with ℎ = 0.1.Figure 4. Exact solutions and approximate solutions for Problem 3 at t = 1 with h = 0.1.

Figure 4 .
Exact solutions and approximate solutions for Problem 3 at = 1 with ℎ = 0.1.

Figure 5 .
Figure 5. Exact solutions and approximate solutions for Problem 4 at t = 1 with h = 0.1.

Figures 7 - 21 Figures 7 -Figure 7 .
Figures 7-10 is a plot of the absolute errors versus r for lower and upper bounds of Problems 1-4.It can be seen in Figures7-9that the errors of 2PDO4 are smaller than the comparison methods, but in Figure10the errors of 2PDO4 are smaller at certain points.

Figure 7 .
Figure 7. Graph of the numerical results when solving Problem 1: (a) Graph of absolute errors versus r for lower bounds; (b) Graph of absolute errors versus r for upper bounds.

Figure 7 .Figure 8 .
Figure 7. Graph of the numerical results when solving Problem 1: (a) Graph of absolute errors versus for lower bounds; (b) Graph of absolute errors versus for upper bounds.

Figure 8 .Figure 9 .Figure 10 .
Figure 8. Graph of the numerical results when solving Problem 2: (a) Graph of absolute errors versus r for lower bounds; (b) Graph of absolute errors versus r for upper bounds.Symmetry 2018, 10, x FOR PEER REVIEW 19 of 21

Figure 9 .Figure 9 .Figure 10 .
Figure 9. Graph of the numerical results when solving Problem 3: (a) Graph of absolute errors versus r for lower bounds; (b) Graph of absolute errors versus r for upper bounds.

Figure 10 .
Figure 10.Graph of the numerical results when solving Problem 4: (a) Graph of absolute errors versus r for lower bounds; (b) Graph of absolute errors versus for upper bounds.

Table 3 .
Numerical results of 2PDO4 and RK4 for solving Problem 3.

Table 5 .
Numerical results of 2PDO4 and RK4 for solving Problem 5.
Note: Simulation results of 2PDO4 and RK4 has been obtained in the same platform.