On the Numerical Solution of Ordinary, Interval and Fuzzy Differential Equations by Use of F-Transform

: An interesting property of the inverse F-transform ˆ f of a continuous function f on a given interval [ a , b ] says that the integrals of ˆ f and f on [ a , b ] coincide. Furthermore, the same property can be established for the restrictions of the functions to all subintervals [ a , p k ] of the fuzzy partition of [ a , b ] used to deﬁne the F-transform. Based on this fact, we propose a new method for the numerical solution of ordinary differential equations (initial-value ordinary differential equation (ODE)) obtained by approximating the derivative · x ( t ) via F-transform, then computing (an approximation of) the solution x ( t ) by exact integration. For an ODE, a global second-order approximation is obtained. A similar construction is then applied to interval-valued and (level-wise) fuzzy differential equations in the setting of generalized differentiability (gH-derivative). Properties of the new method are analyzed and a computational section illustrates the performance of the obtained procedures, in comparison with well-known efﬁcient algorithms.


Introduction
The fuzzy transform (F-transform) of a continuous function f : [a, b] −→ R was introduced by Perfilieva in [1,2]. This special fuzzy method is particularly appealing and useful to handle many real-world problems and an extensive research activity has both analyzed its properties and its fields of applications; for the literature related to this paper we refer to, e.g., [3][4][5][6][7][8][9][10][11] and the references therein.
In recent research, attention has been paid to the numerically approximated solutions of ordinary differential equations (ODEs) · x(t) = F(t, x) of various types. In particular, it is shown that by using the inverse F-transform, it is possible to obtain good approximations of the solution x(t). The methods that use the F-transform are (computationally) superior with respect to other ones such as the second-order Runge-Kutta algorithm or basic multi-step algorithms (see [12][13][14][15][16]). In the final section of this paper, we will present some comments and a preliminary comparative valuation of the proposed methods.
In this paper, we propose a numerical method where the inverse F-transform is used to approximate the derivative · x(t); the solution x(t) is then obtained by exact integration of the approximated derivative: this is allowed by an interesting property which says that the integral of the inverse F-transform of f on [a, b] coincides (exactly) with the integral of the function f itself; this idea was presented in a preliminary form in [17].
For an ordinary differential equation, including the case of a system, a global second-order approximation is obtained. Properties of the new method are analyzed and a computational section illustrates the performance of the obtained procedures, in comparison with well-known efficient algorithms.
A similar construction is then applied to interval-valued (IDE) and (level-wise) fuzzy differential equations (FDE) in the setting of generalized differentiability (gH-derivative, see [18][19][20]). IDEs and FDEs are designed to model uncertainty and its propagation in a dynamical setting and it is well known that the fuzzy case can be expressed in terms of a family of IDEs by adopting the level-wise representation of fuzzy numbers and fuzzy-valued functions. For a modern introduction to FDEs under Hukuhara and generalized derivative, we refer to chapter 9 of Bede's book [21]. The interested reader is also referred to the recent book [22], in particular to chapter 4, and the references therein.
This paper is organized into five sections. In Section 2 we recall some of the basic definitions and properties of the F-transform as contained in [1,2,23]. Then, in Section 3, we describe our approach to the numerical solution of ordinary (and systems of) differential equations with initial conditions, usually referred as Cauchy problem. The F-transform is used to approximate the derivative of the solution to be founded and the unknown function is determined by point-wise by exact integration of the derivative; this section also contains the main properties of the method and a proof of its (global) convergence. Numerical examples and a computational comparison with other well-performing algorithms is presented in Section 4. Section 5 considers the case of interval differential equations and extends the use of F-transform to the numerical solution of IDEs and FDEs under (level-wise) generalized Hukuhara differentiability; in particular, the switching phenomenon is analyzed and rule to manage the switching is proposed and implemented in the proposed procedure. Several computational examples of interval and fuzzy differential equation are presented and discussed in Section 6. Section 7 presents some conclusions and present some ideas for further work.

Preliminaries
We briefly recall the basic definitions and properties of the F-transform setting. For all the details, we refer to the papers [1,2,23].
Let us define the following integrals and I k = I − k + I + k for k = 1, ..., n.
The direct F-transform of f with respect to (P, A) is the n-tuple of real numbers F (P,A) = (F 1 , F 2 , ..., F n ) defined as p k−1 f (t)A k (t)dt, k = 2, ..., n − 1 and, obtained from the direct fuzzy transform F (P,A) , the inverse F-transform (iF-transform) of f is the function f (P,A) : [a, b] −→ R given by The following properties (see [1]) are the fundamentals of the F-transform setting. 1. for any positive real ε, there exists a fuzzy partition (P ε , A ε ) such that the corresponding iF-transform 2. for all k = 1, ..., n, where h = max{p k+1 − p k ; k = 1, ..., n − 1}. 3. the direct and inverse F-transforms are linear, i.e., for any λ∈R and for any two continuous functions f , g : [a, b] −→ R, with direct F-transforms F (P,A) and G (P,A) with respect to the same partition (P, A); then 3.1. the direct F-transforms of λ f and f + g are, respectively, λF (P,A) and F (P,A) + G (P,A) , 3.2. the inverse F-transform of λ f and f + g are, respectively,λ f (P,A) and f (P,A) + g (P,A) .
In [10], the following property has been established: It is interesting to observe that a fuzzy partition has a "nesting" property (its proof is immediate): Proposition 3. Let (P, A) be a fuzzy partition of interval [a, b] with P = {a = p 1 < p 2 < ... < p n = b} and consider any subinterval [a, p k ] with k = 2, ..., n; let (P k , A k ) be the fuzzy partition of [a, p k ] defined by where the last basic function A k is given by the restriction of the basic function A k to the subinterval [p k−1 , p k ]. The fuzzy partition (P k , A k ), k = 1, 2, ..., n is called the k-th nested partition associated with (P, A) and clearly (P n , A n ) = (P, A).
A consequences of the properties above is the following proposition: Proposition 4. Let f : [a, b] −→ R be a continuous function and let F (P,A) = (F 1 , F 2 , ..., F n ) T be its F-transform with respect to (P, A). Then, for any k = 2, ..., n − 1, the F-transform of the restriction of f to the subinterval [a, p k ] (with respect to the k-th nested partition (P k , A k )), is given by where only the last component F k is changed with respect to the components in F (P,A) and is given by We have, for all k = 2, ..., n − 1 (the central summation is assumed to be zero if k = 2), With the notation in (1), equality (10) is written (for k > 1) as and we have: with h = b−a n−1 and let F (P,A) = (F 1 , F 2 , ..., F n ) T be the F-transform of f ; let also F (k) = (F 1 , ..., F k−1 , F k ) T be as in (9) for any k = 2, ..., n − 1. Then dt and, by the trapezoidal integration, . For the second equality (consider that A k (p k+1 ) = 0) we have

Numerical Solution of Initial-Value ODE by F-Transform
Let us consider the following initial-value ordinary differential equation (ODE): We assume the usual requirements on function f (t, x) that ensure existence and unicity of the solution x(t), t ∈ [t 0 , t 1 ].
We are interested to find an approximation of the final value x(t 1 ) of the solution x(t). Let (P, A) be a fixed (but arbitrary) uniform fuzzy partition of [t 0 , In the rest of the paper, we will make use of the following functions and notation: for any basic function A k ∈ A, k = 1, 2, ..., n, we will denote by B k the integral function defined by The components of the direct F-transform of · x(t) will be denoted by F 1 , F 2 , ..., F n and the inverse The following proposition follows immediately from (10).

Proposition 6.
If the exact values (F k ) k=1,...,n of the direct F-transform components of the function t −→ f (t, x (t)) on [t 0 , t 1 ] with respect to the fuzzy partition (P, A) are known, then the final value x (t 1 ) of the solution of (12) is exactly given by Furthermore, at the intermediate points p k , k = 2, ..., n − 1 of the decomposition P, we have that the solution x (p k ) is exactly given by Proof. By definition, we have and, from property (ii) in Proposition 1, also we then obtain By Proposition 2 applied to the derivative function · x (t) = f (t, x (t)), with (P, A) on [t 0 , t 1 ], and using (1)-(2), we have B j (p k ) = I j and B k (p k ) = I − k and the conclusion follows.
From the properties of F-transform as in Proposition 1(point 3.) we finally have the following Proposition 7. If the exact values (F k ) k=1,...,n of the direct F-transform components of the function t −→ f (t, x (t)) on [t 0 , t 1 ] with respect to a fuzzy partition (P, A) are known, then the solution x (t) of (12), for all k = 1, ..., n − 1, satisfies Proof. Apply the identity (14).

An F-Transform Algorithm for ODE
In view of Propositions (6) and (7), we then need a way to compute or to approximate the direct F-transform components (F k ) k=1,...,n of · x (t). By definition, we have that (here p 0 = p 1 and p n+1 = p n ), for k = 1, 2, ..., n, and, from Proposition (5), As a final step, substitute x(p k ) = x 0 + k−1 ∑ j=1 F j I j + F k I − k and obtain, for h −→ 0, Approximated values G k and G k of F k and F k , respectively, can be computed by solving the following equations (we will assume that the sums below will be zero if k = 1) To solve Equations (18) and (19) let us write them for different values of k = 1, ..., n. The first component G 1 can be determined from the initial condition When computing G 2 and G 2 , the value of G 1 is known, and Equations (18) and (19) become From the first equation we determine G 2 and, by substituting into the second, we compute For a general k, the values G 1 , ..., G k−1 are known and we need to solve Equation (19) only for G k then we set G k = 2I − k G k /I k . Summarizing, we determine the values G k and G k iteratively from (23) starting with (20). Each Equation (23) has the form of a fixed-point problem for G k and can be solved by any zero finder routine or, considering that the system is in (nonlinear) triangular form, by any iterative procedure.
We have the following approximation property for the solution of (12).
The last theorem allows development of the following algorithm for the numerical solution of ODE based on F-transform.
Step 2. For k = 1, 2, ..., n, compute the solutions G 1 , ..., G k−1 , G k of Equation (23) and define Step 3. The final value x (P,A) (t 1 ) (corresponding to p n = t 1 ) is the desired approximation of Theorem 1 ensures that the algorithm ODE-FT is (globally) convergent. Clearly, the convergence of the algorithm ODE-FT assumes that the solutions of Equation (23) are solved with high precision, independent on the number n of points p k of the fuzzy partition; in practice, if the exact solutions G 1 , ..., G k−1 , G k are approximated and substituted in the algorithm by quantities G * 1 , ..., G * k−1 , G * k , it is required that they are such that G j − G * j < tol j for small positive tolerances tol j < ε and fixed small ε > 0; in this case, taking into account that I j ≤ 2h, I − j ≤ h and consequently In the computations reported in this paper, we have solved the fixed-point problems either exactly (in the case of linear differential equations) or, in the nonlinear cases, with a tolerance tol ≤ 10 −12 in the absolute difference between two successive iterates of the used equation solver.

Extension to Systems of ODEs
The extension of the described procedure to solve systems of ordinary differential equations with initial conditions in the form We can find an approximation of the final value of each function x i (t 1 ), i = 1, 2, ..., d, in terms of a fixed fuzzy partition (P, A) of [t 0 , t 1 ], e.g., p 1 = t 0 , p k = p k−1 + h, k = 2, ..., n and h = t 1 −t 0 n−1 and basic functions A = {A 1 , .., A n }. Let · x i,(P,A) denote the iF-transform of · x i (t) with direct F-transform components (F i,k ) k=1,...,n , then, the final value x i (t 1 ) can be obtained, in terms of the components (F i,k ) k=1,...,n and the initial condition x i,0 : with respect to a fuzzy partition (P, A) are given, in this case, by d simultaneous systems of equalities (here p 0 = p 1 and p n+1 = p n ), for k = 1, 2, ..., n, and, for all i = 1, ..., d and k = 1, ..., n, we obtain, by substituting for h −→ 0, Approximated values G i,k and G i,k of F i,k and F i,k , respectively, are computed by solving the systems of d equations The first components G i,1 , from the initial condition and, for k ≥ 2, with the known values G i,1 , ..., G i,k−1 , we need to solve the system of d equations for G i,k where we have denoted Each system of Equation (29) has the form of a fixed-point problem for G i,k and can be solved by any zero finder routine or by an iterative procedure similar to the one already described, adapted for the simultaneous case: If G i,k , i = 1, ..., d and k = 1, ..., n, are the solutions of (28)-(29) and we define then, for the solution x i (p k ) of the system of ODEs at the points p k , k = 1, 2, ..., n, we have

The Case of Linear ODE
, then the algorithm ODE-FT can be simplified considerably, because Equation (23) for G k is and one obtains G 1 from I 1 G 1 = h (a(p 1 )x 0 + b(p 1 )) and, for k > 1, remark that the denominator can be controlled to be non-zero by choosing an appropriate (small) value of h; in particular when a(t) has positive values, in order to have 2 h − a(p k ) > 0 it is sufficient to choose h (or, equivalently, n) such that a(p k ) < 2 h for all k = 1, 2, ..., n.

Computational Results for ODEs
The algorithm ODE-FT described in Section 2 has been implemented in MATLAB and used to solve a series of d-dimensional ODEs with d = 1, 2, 3 and 4. For some of them the exact solution x j (t), j = 1, ..., d, is available and the performance of our algorithm is evaluated by comparing the found solutions x j,(P,A) (t) and x j (t) at a number M of points t ∈ {t 1 , t 2 , ..., t M }; for the problems without exact solutions, the comparison is performed with the solution x j,RK (t) found using the well-known Runge-Kutta-Fehlberg (RK) combination of fourth-and fifth-order approximation algorithm available in MATLAB function ode45, one of the possibly best implementations of the explicit RK method with variable step-size adaptation.
Our solution x i at points t i , that approximates x(t i ), is found by applying algorithm ODE-FT successively on each subinterval [t i , t i+1 ], i = 1, ..., 100, starting with the first interval [t 1 , x(t 2 )) and continuing by solving (32) on subsequent subintervals [t i , t i+1 ], i = 2, ..., 100, with the initial condition x(t i ) = x i . On each subinterval, algorithm ODE-FT is applied with a uniform fuzzy partition (P, A), with n = 101 nodes and triangular basic functions A k . We remark that with the decomposition of interval [0, 1] into 101 points and a fuzzy partition of each subinterval [t i , t i+1 ], i = 1, ..., 100, with n = 101 nodes, the resulting (fixed) step size for solving (32) To compare the found solution with the values obtained using routine erf() in MATLAB, we found the average absolute error 1 To test the validity of the ODE-FT algorithm, we have performed two numerical experiments on two series of ODEs.
For the first experiment, we have chosen five non-trivial problems with known exact solution; this allows comparison of the one obtained by ODE-FT with the exact one.
In the second experiment we solve other five typical well-known ODEs and we compare our solutions with the ones obtained by well-established routines in MATLAB, like ode45 (based on explicit single step Runge-Kutta method of orders 4 and 5 with step-size control), ode113 (a variable-step and variable-order method based on Adams-Bashforth-Moulton discretizations of orders 1 to 13) or ode15i (an implicit variable-step and variable-order solver of orders 1 to 5, particularly efficient for stiff problems).

ODEs with Known Exact Solutions
The five ODEs for the first experiment are denoted Ex1 (single equation, d = 1) and Ex2, Ex3, Ex4, Ex5 (systems of two equations, d = 2).

ODEs without Known Closed-Form Solution
The five ODEs for the second experiment are nonlinear d-dimensional systems denoted No1 These ODEs are generally considered to be numerically difficult to solve and are frequently used to test ODE solvers, in comparison with existing well-performing procedures. We have executed algorithm ODE-FT using MATLAB in parallel with well-established routines in MATLAB, in particular with routine ode45 (based on explicit single step Runge-Kutta method of orders 4 and 5 with step-size control) and, in one case, with routines ode113 (a variable-step and variable-order method based on Adams-Bashforth-Moulton discretization of orders 1 to 13) and ode15i (an implicit variable-step and variable-order solver of orders 1 to 5, particularly efficient for stiff problems). Now, as we have explained in Section 2, the performance of algorithm ODE-FT depends essentially on the step size h resulting from the two parameters MFT and N; on the other hand, our simple implementation of ODE-FT does not consider variable step size, designed to control adaptively the local discretization and/or approximation errors. For this reason, we have first executed routine ode45 to compute the solution at an appropriate number M of points t i , i = 1, 2, ..., M and, from the output of ode45, we have computed the number N of nodes for the fuzzy partition in such a way that h FT is of the same order as the optimal step size h ode45 determined by ode45. In this way, the two routines produce solutions by using a step size of the same magnitude and the comparison does not depend on this element; for example, suppose that in solving an ODE on time domain [a, b] at M uniform points t i , i = 1, 2, ..., M, routine ode45 returns, e.g., an optimal step size h ode45 , we then execute algorithm ODE-FT with MFT = M and with N∼ 1 .
We remark that to determine the step size h FT for ODE-FT, we can freely play with the two parameters MFT and N; in our experiment, in order to balance the number of Equation (29) to solve and the memory requirements, we have limited values of N between 11 and 501, so determining M to obtain the required step size, eventually by increasing M to take N in the range above.
Problem No1: Solution interval is t ∈ [0, 5π]; Routine ode45 solves this ODE with optimal step size h ode45  Figure 1 for the graphical representation of the found solutions. For this example, Figure 2 pictures the differences FT(t) − RK(t) between the FT solution and the RK solution for the three components x j (t), j = 1, 2, 3 at the MFT points. It is interesting to observe that the differences oscillate in sign, showing that in absolute value the differences are not systematically diverging.  Problem No2: (Van der Pol equations) Solution interval is t ∈ [0, 300]; The parameters are µ = 50, a = 3, ω = π 5 . This system is considered to be a hard ODE and requires very small step size to determine points where the solution changes suddenly (see Figure 3). To capture this phenomenon, we use M = 30,001. The solution by ode45 has optimal step size h ode45 = 2.196 × 10 −5 and final solution vector The parameters are α = 0.2, β = 0.2, γ = 5.0. The third equation is nonlinear, but this problem is considered to be not numerically hard. The solution by ode45 has optimal step size h ode45 = 0.
The parameters are a = 10, b = 28.0, c = 8 3 . It is well known that this nonlinear system is hard to solve numerically as it exhibits chaotic trajectories. Routine ode45 solves the Lorenz system by optimal step size h ode45 = 9.37 × 10  Figure 5 for the trajectories of components x j (t), j = 1, 2, 3 and Figure 6 for the 3D representations of the solutions.  We have also solved this system by a run of ODE-FT with MFT = 30,001, N = 101 and by running the two MATLAB routines ode113 and ode15i; the resulting solutions are visualized in Figure 7 and we see that all the routines tend to generate trajectories with very different behavior for large values of time t. Figure 7. Problem No4 (Lorenz system): The three components x j (t), j = 1, 2, 3 are obtained also by routines ode15i (cyan color) and ode113 (green color); FT solution is red-colored, RK solution is blue-colored. Remark that differences between the solutions of this system are essentially due to very small differences in the solutions found by the different methods.

Interval (Fuzzy) Differential Equations and F-Transform
In this section, we consider interval and fuzzy differential equations in the setting of generalized Hukuhara differentiability as described in [18,19,25]. A preliminary version of this section has been presented as a conference paper in [17].
Compact intervals of real numbers will be denoted by the usual endpoints notation A = [a − , a + ], B = [b − , b + ] or by the midpoint notation A = ( a; a), B = ( b; b) where a = 1 2 (a − + a + ) is the midpoint and a = 1 2 (a + − a − ) is the radius (half-length). The set of all compact real intervals will be denoted by K C .
The use of midpoint representation for intervals has been recently adopted to study several topics in the analysis of interval-valued functions, the single variable case (see [26,27] and the references therein), to which the interested reader is referred for a complete description of interval-valued generalized Hukuhara derivative (gH-derivative for short).
Given two intervals A, B ∈ K C , the gH-difference is the interval C ∈ K C (it always exists and is unique) such that Using midpoint notation, we have Please note that if A i gH B i are gH-differences of the same type for all i = 1, ..., n, i.e., all satisfy either (i) or (ii) above, then (see [25]) An interval-valued function F : [a, b] → K C will be denoted by F(t) = [ f − (t), f + (t)] or, in equivalent midpoint notation, by F(t) = ( f (t); f (t)), for t ∈ [a, b].
We will denote by R F the set of fuzzy numbers, i.e. normal, fuzzy convex, upper semi-continuous and compactly supported fuzzy sets defined over the real line R. The α-level set of u (or simply its α-cut) is defined by [u] α = {x|x ∈ R, u(x) ≥ α} and, for α = 0, it is the closure of the support A well-known result (see [21]) allows us to represent a fuzzy number as a pair u = (u − , u + ) of functions u − , u + : [0, 1] −→ R, defining the endpoints of the α-cuts as [u] α = [u − α , u + α ] = ( u α ; u α ) We refer to functions u − (.) and u + (.) as the lower and upper branches of u, respectively; the two functions u (.) , u (.) are the (level-wise) midpoint and radius functions.
Given fuzzy numbers u, v ∈ R F , the level-wise generalized Hukuhara difference (LgH-difference for short) is the family of intervals The following definition of generalized derivative can be applied both to interval-valued function or to the level-cuts of a fuzzy-valued function.

Definition 1 ([18]). Let t 0 ∈]a, b[ and h be such that t
If F : [a, b] → K C is an interval-valued function, its gH-derivative at t 0 is defined to be the limit, if it exists, -If F :]a, b[→ R F is a fuzzy-valued function, its level-wise gH-derivative (LgH-derivative for short) at t 0 is defined to be the family of the gH-derivatives of [F] α , if they exist for all α ∈ [0, 1], i.e., For an interval-valued function F : ], when f − (t) and f + (t) are both differentiable, we can distinguish two cases, corresponding to (i) and (ii) of (33) (see [18]) If F :]a, b[→ R F is a fuzzy-valued function, we define analogously (i)-LgH and (ii)-LgH differentiability of F, with the additional requirement that (38) (or (39), respectively) are valid for all [F] α .

As in [18], we say that a point t 0 ∈]a, b[ is an l-critical point of F if it is a critical point for the length function len([F(t)])
[ is a switching point for the gH-differentiability of F, if in any neighborhood V of t 0 there exist points t 1 < t 0 < t 2 such that type-I switch point): at t 1 (38) holds while (39) does not hold and at t 2 (39) holds and (38) does not hold, or type-II switch point): at t 1 (39) holds while (38) does not hold and at t 2 (38) holds and (39) does not hold. Analogous definitions can be given level-wise for a fuzzy-valued function.

Numerical Interval ODE by F-Transform
An interval differential equation (IDE) with initial condition can be written in the form We will approximate · x gH (t) by the iF-transforms of the two functions · x − gH (t) and · x + gH (t) on the same fuzzy partition (P,A), i.e., where for j = 1, 2, ..., n, using the same notation as in Section 2 with [a, b] = [t 0 , t 1 ], F − j and F + j are the direct F-transforms of ( · x − gH ) and ( · x + gH ). From the monotonicity properties of F-transform (see [4,10] ), we know that x + gH (t) for all t) and we can define the intervals F j = [F − j , F + j ]. Consequently, in terms of interval arithmetic operations, we also have the following approximation of the (interval-valued) On the other hand, the following function is well defined: and H(t 0 ) = 0 (here 0 stands for the interval [0, 0]). The interval-valued function H(t) will play a central role in our method to solve the IDE (40).
For simplicity, define the following interval-valued function Clearly, ϕ is a continuous interval-valued function and, from Theorem 43(i) in [18], the integral The following property is proved in [17].

Proposition 8. Consider the function H(t) in (43) and define the two interval-valued functions Φ(t) and Ψ(t)
Consequently, if function ϕ : [t 0 , t 1 ] −→ K C is generated by a fixed fuzzy partition (P,A) as in (44) and H(t) = ∑ n j=1 F j B j (t), then we have the following two cases: (1) If the intervals F j , j = 1, ..., n are such that (46) is a solution of (40).
(2) If the intervals F j , j = 1, ..., n are such that (45) is a solution of (40), provided that all the gH-differences H (t) gH (−x 0 ) are of the same type for all t ∈ [t 0 , t 1 ].

Remark 2.
From the properties of gH-difference, we have that It is interesting to observe that function Ψ (t) is (i)-gH-differentiable at all points, while function has always increasing length, but the same is not true for solution Finally, let G − j and G + j be the O(h 3 ) approximations, respectively, of F − j and F + j similar to (18)- (19); we can see that the interval-valued iF-transform approximation ( · x gH ) (P,A) (t) = ∑ n j=1 G j A j (t) of · x gH is useful in determining conditions for a switching point.
Indeed, if p k−1 < p k < p k+1 are three adjacent points of P (2 ≤ k ≤ n − 1) and we suppose that no switching point exists internally to the two subintervals [p k−1 , p k ] and [p k , p k+1 ] (i.e., possibly, the switching is exactly at p k ) then the solution x(t) must satisfy On the other hand, we have . Suppose now that p k is a switching point for the gH-differentiability of x(t). We have two cases: [. In terms of midpoint notation for intervals, we can summarize the discussion above as follows: Types of switching points: Let (P,A) be a fuzzy partition of [t 0 , t 1 ] and p k−1 < p k < p k+1 (2 ≤ k ≤ n − 1); let x(t) = ( x(t); x(t)) be a solution to (40) and H L k = ( H L k , H L k ) and H R k = ( H R k , H R k ) be given as in (48). Then, the midpoint values satisfy Assuming that only p k is eventually a switching point, we have the following four cases (a) if p k is a (i)-to-(ii) switch, then There are no general rules to locate a switching point. Denote by x(t) = [x − (t), x + (t)] = ( x(t); x(t)) a solution of the IDE (40); if x(t) is (i)-gH-differentiable, its length x(t) will not decrease, while x(t) will not increase if x(t) is (ii)-gH-differentiable. Some authors have noticed that possibly, the sequence of switching points can be pre-defined a priori, by positioning them in the time domain of the interval differential equation; this is true, at least in principle, provided that the found solution is guaranteed to have exactly them and not other switching points, but such purely exogenous proposal is not fully convincing.
An endogenous way may be preferred, similar to control strategies, to connect the type of gH-differentiability to the evolution of the trajectory. For example, it seems reasonable to locate the switching points depending on how the solution is evolving, by fixing a lower l(t) and an upper threshold u(t), say 0 ≤ L ≤ l(t) ≤ u(t) ≤ U and requiring that l(t) ≤ x(t) ≤ u(t) for all t; then, (a) a (i)-to-(ii) switch is decided at t = p k if the following condition is reached In some sense, the rule above will control the increasing and decreasing of "uncertainty" in an endogenous way, without any reference to the interval initial condition x 0 or to the interval-valued function F(t, x).
Furthermore, we are essentially free to decide the type of differentiability at the initial point; if the initial length x 0 is such that L < x 0 < U, we can start either with a type c) or a type d) and a unique solution is then found by applying the decided switching rule.
A different purely endogenous rule can be obtained by following the increase or decrease of x(t) and trying to intercept points t * where the function x(t) has a local maximum (for a (i)-to-(ii) switch) or a local minimum (for a (ii)-to-(i) switch). A necessary condition can be obtained according to the following simple result: Endogenous criteria for a switching point: Assume that x(t) = [x − (t), x + (t)] = ( x(t); x(t)) is such that x − (t) and x + (t) are differentiable so that its gH-derivative · x gH (t) can be expressed in terms of the derivatives (x − ) (t) and (x + ) (t). Let (P, A) be a fuzzy partition of [t 0 , t 1 ] and let ( . On the other hand, from the differentiability of x − (t) and x + (t), we have 0 = d dt x(t)) for t = p k , i.e., ( We then suggest the following (purely endogenous) switching rule.
Switching rule for IDE: Let (P, A) be a fuzzy partition of [t 0 , t 1 ] and suppose we have computed the approximate solution x(p i ) of IDE (40) at the points t 0 = p 1 < p 2 < ... < p k with k < n. Assume that h = t 1 −t 0 n−1 is sufficiently small and p i = p i−1 + h. In any case, the midpoint value is computed as x(p k ) = x(p k−1 ) , then p k is assumed to be a (i)-to-(ii) switch and the next iteration is performed according to (ii)-gH-differentiability.
, then p k is assumed to be a (ii)-to-(i) switch and the next iteration is performed according to (i)-gH-differentiability.

Numerical Procedure for IDE
We are now ready to summarize the proposed procedure to solve the IDE (40) on [t 0 , t 1 ], based on F-transform of the gH-derivative · x gH (t). Chose a uniform fuzzy partition (P,A) with n > 2 points P= {t 0 = p 1 < ... < p n = t 1 } and basic functions A= {A 1 , ..., A n }.
From the properties of the F-transforms on (P k ,A k ), we have for all k = 2, ..., n, where With the introduced notation, we solve (40) by the following procedure, assuming that no switching point exists on [p k−1 , p k ] and x 0 is the given initial condition.
Step 4. (Test for a switching at p k ) If F k < ε, then p k is declared to be a switch of gH-differentiability: if so, exchange the value of Meth from 1 to 2 or from 2 to 1; set H = f (p k , x(p k ))I + k , y = x(p k ). If instead F k > ε, then p k is not a switch: in this case, update H to H + F k I − k , without changing Meth.
Step 5. Continue from step 3 with next k, until the final point t 1 is reached.

Application to Fuzzy Differential Equations
If the fuzzy-valued function ϕ : [a, b] −→ R F is gH-differentiable with respect to t, or if it is level-wise gH-differentiable, we can approximate its gH-derivative ϕ gH or its LgH-derivative ϕ LgH by a family of interval-valued iF-transforms, parametrized with α ∈ [0, 1], in terms of its α-cuts (consider that A j (x) ≥ 0 for all x and all j) The F-transform of ϕ gH is then obtained by the family of F-transforms of the lower and upper and, considering that each A j is a continuous non-negative function, we can obtain an approximation of ϕ(x) by Consider now an FDE in gH-derivative form the approximation, we obtain Integrating with given initial condition ϕ(x 1 ) = ϕ 1 we obtain and, at the points x j ∈ P, defining ϕ j = ϕ(x j ), j = 2, ..., q, and, for j = 3, ..., n These are (generally nonlinear) non-differential fuzzy equations of the form where when solving for ϕ j , the fuzzy numbers ϕ 1 , ϕ 2 , ..., ϕ j−1 are determined in previous steps and the fuzzy quantity is known. If we search for a (i)-gH-differentiable solution, the fuzzy equation for ϕ j is if we search for a (ii)-gH-differentiable solution, the fuzzy equation for ϕ j is In terms of α-cuts, they can be written (for α ∈ [0, 1] or α ∈ {α 1 , α 2 , ..., α p }) respectively.

Computational Results for Interval and Fuzzy Differential Equations
For all problems we use M = 10, 001 as the number of points t i where the solution is computed and n = 11 as the number of nodes in the partition (P, A); it results that the step size of the method to The number of α-cuts is N α = 11 so that the LgH-differentiable solutions are computed for all The fuzzy initial condition is given in terms of the α-cuts of x 0 , [x 0 ] α = [1 + 0.5α, 3 − α], α ∈ [0, 1], with support [1,3] and core [1.5, 2].
In this case, the interval-valued solutions for the different values of α form a fuzzy-valued function, i.e., the family of IDEs (parametrized with α ∈ [0, 1]) solve the corresponding FDE.
For the two runs with Meth = 1 and Meth = 2, the step size is h = 2.000 ×10 −5 . The solution intervals (initial and final) and the solution intervals at the internal switching points are inserted in the following tables, for three values of α = 0, 0.5, 1.
Starting with Meth = 1, i.e., by an increasing-length solution (see Figure 9), the algorithm finds seven switching points at tw ∈ Sw = {.25, .5, .75, 1, 1.25, 1.5, 1.75} (see Table IDE1(a)).  Starting with Meth = 2, i.e., by a decreasing-length solution (see Figure 10), we find the same seven switching points (but with different interval values) as reported in Table IDE1 (b). In both cases, the fuzzy solution exists and is gH-differentiable. It seems remarkable that in this example, the two solutions are different, but share the same value at some of the switching points; indeed, denoting x 1 and x 2 the solution found with Meth = 1 and Meth = 2, we have that x 1 (.5) = x 2 (.5), x 1 (1) = x 2 (1), x 1 (1.5) = x 2 (1.5) while the solutions are not equal at all other points.   x where Ax(t) and c(t)B are obtained by standard interval operations and gH is gH-difference; the interval initial condition in terms of α-cuts is [ In this case, the step size is h = 4.712 ×10 −5 . The computations are performed with 11 α-cuts; in the next two tables we insert the results obtained with Meth = 1 and Meth = 2 for a subset of the computed α-cuts.
Remark the important fact that in this example, the switching points change in position for different values of α (see Figures 11 and 12 and Tables IDE2 (a)-(b)). Viewing the solution function as fuzzy-valued, it is not gH-differentiable and there exist three switching regions.   .
The α-cuts of the initial condition x 0 and the final fuzzy solution for Meth = 1 and Meth = 2 are inserted in Table FDE1.   The fuzzy solution is periodic with period T = 2π. For both Meth = 1 and Meth = 2, there are three internal switching points at tw ∈ Sw = {π, 2π, 3π}, where the length of .
x gH is zero and the length of the solution x(t) is maximal (at tw = π, 3π) or minimal (at tw = 2π.) At points t = 2π and t = 4π the solution coincides with the initial condition (see Figures 15 and 16).

Concluding Comments and Further Work
In this paper, we see that the F-transform approximation setting allows good numerical procedures to solve ordinary differential equations (ODEs) and to approach the numerical solution of interval and fuzzy differential equations. The computational comparison of the proposed ODE-FT method with other well-known and well behaving numerical routines available in MATLAB, such as ode45, ode15i or ode113, positions F-transform among the most promising mathematical tools for the approximation of functions.
It seems that in general, the different algorithms behave similarly, at least for the chosen step size h = 0.1 and h = 0.01 (consider that such h is a big one and values around h = 0.00001 or less are more adequate for a comparison, as suggested, e.g., by the values used in routine ode45 that chooses h dynamically). Possibly, more efficient and elaborated implementations of the proposed algorithms will require some additional analysis of F-transform properties to allow variable-order and/or step-size control.
As a tool for numerical solution of interval (IDEs) and fuzzy (FDEs) differential equations in terms of gH-derivative, the F-transform allows an immediate approach to handle the switching phenomenon, a still open problem in this area. This approach is obtained by the application of Equation (43) and proposition 8, which is possible because the interval-valued function H(t) offers an approximation of the interval solution x(t) at all points t ∈ [t 0 , t 1 ] and not only at the discretized points p k , as usual in the (explicit) single-or multi-step ODE solvers.
Further research can be planned in the design and experimentation of efficient numerical procedures to solve real-world applications. In this direction, a possible improvement in the approximation can be obtained by higher-order F d -transform (see [11] for recent results of its properties), by introducing local polynomials or, more generally, local parametric functions F k (t; θ) in place of constant direct components F k (coefficients of the polynomials or parameters θ ∈ R d are then estimated by least squares). In these cases, the inverse F-transform of a function f (t) on [a, b] has the form f d (P,A) (t) = n ∑ k=1 F k (t; θ (k) )A k (t), with estimated parameters θ (k) for the k-th direct component.
It is worth to remark that the same integral property used in this paper for the standard F-transform f (P,A) (t) is also valid for f d (P,A) (t), i.e., It should be interesting to see if higher-order F-transform approximations will be able to generate high orders O(h q ), q > 2 of convergence, and to obtain possibly increasing orders q by increasing d (two numeric schemes of order q = 2 based on the F 2 -transform are obtained in [16]).
Similar results can be obtained by considering the discrete F-transform f (P,A) (t j ) on a data set of points S = (t j , f j ), j = 1, 2, ..., m ; the integrals are substituted by summations and we have (see [10]) Finally, it is worth mentioning the possibility of applying the ideas presented in this paper to the numerical solution of other kinds of differential and integral equations, such as delay differential equations (e.g., [28]), differential equations on time scales, fractional differential equations ( [29]) and implicit differential algebraic equations (DAE, see, e.g., [30,31]); a general DAE with additional constraints, on an interval [t 0 , t 1 ] is expressed in the form In this cases, the discretization of [a, b] by a fuzzy partition (P, A) and the substitution of · x (t) and x(t) with the functions · x (P,A) (t) and x (P,A,x 0 ) (t) at points p k ∈ P will transform the DAE into a standard feasibility problem, consisting of finding feasible solutions for the transformed system at points p k , k = 1, ..., n.
Author Contributions: All authors contributed equally to the final version of this paper. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.