On the numerical solution of ordinary, interval and fuzzy differential equations by use of F-transform

axioms Article On the Numerical Solution of Ordinary, Interval and Fuzzy Differential Equations by Use of F-Transform Davide Radi 1,2, Laerte Sorini 3and Luciano Stefanini 3,* 1Department of Economics and Management, University of Pisa, Via C. Ridolfi, 10, 56124 Pisa (PI), Italy; [email protected] or [email protected] 2Department of Finance, Faculty of Economics, VŠB—Technical University of Ostrava, Sokolská tr. 33, 70121 Ostrava, Czech Republic 3DESP, Department of Economics, Society, Politics, University of Urbino Carlo Bo, Via A. Saffi 42, 61029 Urbino, Italy; [email protected] *Correspondence: [email protected] Received: 24 November 2019; Accepted: 28 January 2020; Published: 5 February 2020   Abstract: 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 , pk] of the fuzzy partition of [a , b] used to define 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 efficient algorithms. Keywords: F-transform; initial-value ODE; numerical ODE solver; interval differential equations; gH-Derivative; fuzzy differential equations 1. 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–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 – 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]. Axioms 2020,9, 15; doi:10.3390/axioms9010015 www.mdpi.com/journal/axioms Axioms 2020,9, 15 2 of 37 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 – 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 2we 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 5considers 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. 2. 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]. A fuzzy partition (P , A) of the compact interval [a , b] is defined by a finite decomposition P= {a=p1<p2< ... <pn=b} of [a , b] with n points and by a family A={A1 , A2 ,..., An} of n continuous basic functions Ak:[a , b]−→ [ 0,1 ] with the following properties (the decomposition P is not required to be uniform): 1. Ak(t) = 0 for t/∈]pk−1 , pk+1[ ( k= 2,..., n− 1), A1(t) = 0 for t/∈[p1 , p2[ , An(t) = 0 for t/∈ ]pn−1,pn]; 2. Ak(pk) = 1 for all k=1, 2, ..., nand A1(t) + A2(t) + ... +An(t) = 1 for all t∈[a,b]; 3. for k= 2,..., n− 1, Ak is increasing on [pk−1 , pk] and decreasing on [pk , pk+1] , A1 is decreasing on [p1,p2],Anis increasing of [pn−1,pn]. Let us define the following integrals      I− 1=0 , I+ 1=Rp2 p1A1(t)dt I− k=Rpk pk−1Ak(t)dt ,I+ k=Rpk+1 pkAk(t)dt I− n=Rpn pn−1A1(t)dt ,I+ n=0 (1) and Ik=I− k+I+ kfor k=1, ..., n. (2) Axioms 2020,9, 15 3 of 37 The direct F-transform of f with respect to (P , A) is the n -tuple of real numbers F(P,A)= (F1,F2,..., Fn)defined as F1=1 I1 p2 R p1 f(t)Ak(t)dt,Fn=1 In pn R p1 f(t)Ak(t)dt, (3) Fk=1 Ik pk+1 R pk−1 f(t)Ak(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]−→ Rgiven by f(P,A)(t) = n ∑ k=1FkAk(t)for t∈[a,b]. (4) The following properties (see [1]) are the fundamentals of the F-transform setting. Proposition 1. (from [1]) If f :[a,b]−→ Ris a continuous function then; 1. for any positive real ε , there exists a fuzzy partition (Pε , Aε) such that the corresponding iF-transform f(Pε,Aε):[a,b]−→ Rsatisfy f(t)−f(Pε,Aε)(t)<εfor all t ∈[a,b]. 2. for all k =1, ..., n, Fk=f(P,A)(pk)(5) f(pk) = Fk+O(h2)as h −→ 0 (6) where h =max{pk+1−pk;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: Proposition 2. Let f , g:[a , b]−→ R be continuous functions and let (P , A) be a fuzzy partition of [a , b] . Then (i) the iF-transform satisfies b Ra f(P,A)(t)dt = b Ra f(t)dt (7) (ii) if f ang g have the same direct F-transform components F(P,A)=G(P,A), then Zb af(t)dt =Zb ag(t)dt. It is interesting to observe that a fuzzy partition has a “nesting” property (its proof is immediate): Axioms 2020,9, 15 4 of 37 Proposition 3. Let (P , A) be a fuzzy partition of interval [a , b] with P={a=p1<p2<... <pn=b} and consider any subinterval [a,pk]with k =2, ..., n; let (Pk,Ak)be the fuzzy partition of [a,pk]defined by Pk={a,p2,..., pk}(8) Ak={A1,..., e Ak} where the last basic function e Ak is given by the restriction of the basic function Ak to the subinterval [pk−1 , pk] . The fuzzy partition (Pk , Ak) , k= 1,2,..., n is called the k -th nested partition associated with (P , A) and clearly (Pn,An) = (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)= (F1 , F2 ,..., Fn)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 , pk] (with respect to the k -th nested partition (Pk , Ak) ), is given by F(k)= (F1 ,..., Fk−1 , e Fk)T where only the last component e Fkis changed with respect to the components in F(P,A)and is given by e Fk= pk Ra f(t)Ak(t)dt pk Ra Ak(t)dt .(9) We have, for all k =2, ..., n−1(the central summation is assumed to be zero if k =2), f(Pk,Ak)(t) = F1A1(t) + k−1 ∑ j=2 FjAj(t) + e Fke Ak(t) and pk Ra f(t)dt =F1 pk Ra A1(t)dt + k−1 ∑ j=2 Fj pj+1 Ra Aj(t)dt +e Fk pk Ra Ak(t)dt. (10) With the notation in (1), equality (10) is written (for k>1) as pk Ra f(t)dt =F1I+ 1+ k−1 ∑ j=2 FjIj+e FkI− k(11) and we have: Proposition 5. Let f:[a , b]−→ R be a continuous function, let (P , A) be a uniform fuzzy partition of [a , b] with h=b−a n−1 and let F(P,A)= (F1 , F2 ,..., Fn)T be the F-transform of f ; let also F(k)= (F1 ,..., Fk−1 , e Fk)T be as in (9) for any k =2, ..., n−1. Then I− ke Fk=h 2f(pk) + O(h3)as h −→ 0 IkFk=h f (pk) + O(h3)as h −→ 0. Proof. Consider 2 ≤k≤n−1. We have I− ke Fk= pk R pk−1 f(t)Ak(t)dt and, by the trapezoidal integration, I− ke Fk=h 2[f(pk−1)Ak(pk−1) + f(pk)Ak(pk)] + O(h3); Axioms 2020,9, 15 5 of 37 on the other hand, Ak(pk) = 1 and Ak(pk−1) = 0, so I− ke Fk=h 2f(pk) + O(h3) . For the second equality (consider that Ak(pk+1) = 0) we have IkFk=I− ke Fk+h 2f(pk) + O(h3) =h 2f(pk) + h 2f(pk) + O(h3) = h f (pk) + O(h3). 3. Numerical Solution of Initial-Value ODE by F-Transform Let us consider the following initial-value ordinary differential equation (ODE): (· x(t)=f(t,x(t)) ,t∈[t0,t1] x(t0)=x0(12) We assume the usual requirements on function f(t , x) that ensure existence and unicity of the solution x(t),t∈[t0,t1]. We are interested to find an approximation of the final value x(t1) of the solution x(t) . Let (P , A) be a fixed (but arbitrary) uniform fuzzy partition of [t0,t1] with p1=t0 , pk=pk−1+h , k= 2,..., n and h=t1−t0 n−1 ; let A={A1,.., An} be the basic functions. Let · x(P,A) denote the iF-transform of · x(t) with (exact) direct F-transform components (Fk)k=1,...,n. In the rest of the paper, we will make use of the following functions and notation: for any basic function Ak∈A,k=1, 2,..., n, we will denote by Bkthe integral function defined by Bk(t)=Zt t0 Ak(s)ds. The components of the direct F-transform of · x(t) will be denoted by F1 , F2 ,..., Fn and the inverse F-transform of · x(P,A)of is the function · x(P,A)(t) = n ∑ k=1FkAk(t)for t∈[a,b]. (13) The following proposition follows immediately from (10). Proposition 6. If the exact values (Fk)k=1,...,n of the direct F-transform components of the function t−→ f(t,x(t)) on [t0,t1] with respect to the fuzzy partition (P , A) are known, then the final value x(t1) of the solution of (12) is exactly given by x(t1)=x0+ n ∑ k=1 FkIk. Furthermore, at the intermediate points pk , k= 2,..., n− 1of the decomposition P , we have that the solution x (pk)is exactly given by x(pk)=x0+ k−1 ∑ j=1 FjIj+e FkI− k(14) Proof. By definition, we have · x(P,A)(t)= n ∑ k=1 FkAk(t) Axioms 2020,9, 15 6 of 37 and, from property (ii) in Proposition 1, also Zt1 t0 · x(P,A)(t)dt =Zt1 t0 · x(t)dt; we then obtain x(t1)=x0+Zt1 t0 · x(P,A)(s)ds =x0+ n ∑ k=1 FkBk(t1). By Proposition 2 applied to the derivative function · x(t)=f(t,x(t)) , with (P , A) on [t0,t1] , and using (1)–(2), we have Bj(pk)=Ijand Bk(pk)=I− kand the conclusion follows. Remark 1. Consider that in general we have x(t)6=x0+Rt t0· x(P,A)(s)ds, for t /∈ {p1,..., pn}. From the properties of F-transform as in Proposition 1 (point 3.) we finally have the following Proposition 7. If the exact values (Fk)k=1,...,n of the direct F-transform components of the function t−→ f(t,x(t)) on [t0,t1] with respect to a fuzzy partition (P , A) are known, then the solution x(t) of (12), for all k=1, ..., n−1, satisfies x(t)=x(P,A,x0)(pk) + Zt pk f(s,x(s)) ds, for all t ∈[pk,pk+1[, where x(P,A,x0)(t) = x0+Zt t0 · x(P,A)(s)ds (15) Proof. Apply the identity x(t)=x0+Rpk t0f(s,x(s)) ds +Rt pkf(s,x(s)) ds and (14). 3.1. An F-Transform Algorithm for ODE In view of Propositions 6and 7, we then need a way to compute or to approximate the direct F-transform components (Fk)k=1,...,nof · x(t). By definition, we have that (here p0=p1and pn+1=pn), for k=1, 2,..., n, IkFk=Zpk+1 pk−1 · x(t)Ak(t)dt =Zpk+1 pk−1 f(t,x(t))Ak(t)dt I− ke Fk=Zpk pk−1 · x(t)Ak(t)dt =Zpk pk−1 f(t,x(t))Ak(t)dt and, from Proposition 5, IkFk=h f (pk,x(pk)) + O(h3) I− ke Fk=h 2f(pk,x(pk)) + O(h3). As a final step, substitute x(pk) = x0+ k−1 ∑ j=1 FjIj+e FkI− kand obtain, for h−→ 0, IkFk=h f pk,x0+ k−1 ∑ j=1 FjIj+e FkI− k!+O(h3)(16) I− ke Fk=h 2f pk,x0+ k−1 ∑ j=1 FjIj+e FkI− k!+O(h3). (17) Axioms 2020,9, 15 7 of 37 Approximated values Gk and e Gk of Fk and e Fk , respectively, can be computed by solving the following equations (we will assume that the sums below will be zero if k=1) IkGk=h f pk,x0+ k−1 ∑ j=1 GjIj+e GkI− k!(18) I− ke Gk=h 2f pk,x0+ k−1 ∑ j=1 GjIj+e GkI− k!. (19) To solve Equations (18) and (19) let us write them for different values of k= 1,..., n . The first component G1 can be determined from the initial condition · x(t0)=f(t0 , x(t0)) = f(p1 , x0) , obtained from I1F1=h f (p1,x(p1)) + O(h3)by ignoring the term O(h3): G1=h f (t0,x0) I1. (20) When computing G2and e G2, the value of G1is known, and Equations (18) and (19) become I− 2e G2=h 2fp2,x0+G1I1+e G2I− 2(21) I2G2=h f p2,x0+G1I1+e G2I− 2. (22) From the first equation we determine e G2 and, by substituting into the second, we compute G2=2I− 2e G2/I2. For a general k , the values G1 ,..., Gk−1 are known and we need to solve Equation (19) only for e Gk I− ke Gk=h 2f pk,x0+ k−1 ∑ j=1 GjIj+e GkI− k!; (23) then we set Gk=2I− ke Gk/Ik. Summarizing, we determine the values e Gk and Gk iteratively from (23) starting with (20). Each Equation (23) has the form of a fixed-point problem for e Gk e Gk=h 2I− k f pk,x0+ k−1 ∑ j=1 GjIj+e GkI− 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). Theorem 1. Let Gk= (G1, ..., Gk−1,e Gk)be solutions of (20)–(23) and define xGk(pk)=x0+ k−1 ∑ j=1 GjIj+e GkI− kfor k =1, 2, ..., n. Then, for the solution x (pk)of (12) at the points pk, k =1, 2,..., n, we have x(pk)=xGk(pk)+O(kh3)as h −→ 0, x(t1) = xGn(t1)+ (t1−t0)O(h2)as h −→ 0. Axioms 2020,9, 15 8 of 37 Proof. From (14), (24) and (25) we have |x(pk)−xGk(pk)| ≤ k−1 ∑ j=1|Fj−Gj|Ij+|e Fk−e Gk|I− k = k ∑ j=1 O(h3) = kO(h3)as h−→ 0. From h=t1−t0 n−1 we get n= (t1−t0)O(1 h) and nO(h3)=(t1−t0)O(1 h)O(h3)=(t1−t0)O(h2) as h−→ 0 and the conclusion follows. The last theorem allows development of the following algorithm for the numerical solution of ODE based on F-transform. Algorithm ODE-FT : Find an approximated final value x(P,A)(t1) of the ODE · x(t)=f(t,x(t)) , t∈[t0,t1]with initial condition x(t0)=x0. Step 1. Choose a uniform fuzzy partition (P , A) of [t0,t1] with n points p1=t0 , pk=pk−1+h , k=2, ..., nand h=t1−t0 n−1; let I− k,I+ kand Ikas in (1)-(2). Step 2. For k=1, 2, ..., n, compute the solutions G1,..., Gk−1,e Gkof Equation (23) and define x(P,A)(pk)=x0+ k−1 ∑ j=1 GjIj+e GkI− kfor k=1, 2,..., n. Step 3. The final value x(P,A)(t1) (corresponding to pn=t1 ) is the desired approximation of x(t1) with x(t1)−x(P,A)(t1)= (t1−t0)O(h2)as h−→ 0. Theorem 1ensures 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 pk of the fuzzy partition; in practice, if the exact solutions G1 ,..., Gk−1 , e Gk are approximated and substituted in the algorithm by quantities G∗ 1 ,..., G∗ k−1 , e G∗ k , it is required that they are such that Gj−G∗ j<tolj for small positive tolerances tolj<ε and fixed small ε> 0; in this case, taking into account that Ij≤ 2 h , I− j≤h and consequently k−1 ∑ j=1|Gj−G∗ j|Ij+|e Gk− e G∗ k|I− k≤(2k−1)εh=2k−1 n−1(t1−t0)ε, the precision of ODE-FT is such that |x(pk)−xG∗ k(pk)| ≤ k−1 ∑ j=1|Fj−Gj|Ij+ k−1 ∑ j=1|Gj−G∗ j|Ij+|e Fk−e Gk|I− k+|e Gk−e G∗ k|I− k = k ∑ j=1 O(h3) + k−1 ∑ j=1|Gj−G∗ j|Ij+|e Gk−e G∗ k|I− kas h−→ 0 ≤kO(h3) + 2k−1 n−1(t1−t0)εas h−→ 0 and, for k=n , |x(pn)−xG∗ n(pn)| ≤ (t1−t0)O(h2+ 2 ε)as h−→ 0 and n−→ ∞ . 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. Axioms 2020,9, 15 9 of 37 3.2. Extension to Systems of ODEs The extension of the described procedure to solve systems of ordinary differential equations with initial conditions in the form      . xi(t) = fi(t,x1(t),..., xn(t)),i=1, 2,..., d xi(t0) = xi,0,i=1, 2,..., d t∈[t0,t1] We can find an approximation of the final value of each function xi(t1) , i= 1,2,..., d , in terms of a fixed fuzzy partition (P , A) of [t0,t1] , e.g., p1=t0 , pk=pk−1+h , k= 2,..., n and h=t1−t0 n−1 and basic functions A={A1,.., An} . Let · xi,(P,A) denote the iF-transform of · xi(t) with direct F-transform components (Fi,k)k=1,...,n , then, the final value xi(t1) can be obtained, in terms of the components (Fi,k)k=1,...,nand the initial condition xi,0: xi,(P,A)(t) = x0+Zt t0 · xi,(P,A)(s)ds =x0+ n ∑ k=1 Fi,kBk(t) where Bk(t)=Rt t0Ak(s)ds,k=1, ..., n. The direct F-transform components of the functions t−→ fi(t,x1(t),..., xd(t)) on [t0,t1] with respect to a fuzzy partition (P , A) are given, in this case, by d simultaneous systems of equalities (here p0=p1and pn+1=pn), for k=1, 2, ..., n, IkFi,k=Zpk+1 pk−1 · xi(t)Ak(t)dt =Zpk+1 pk−1 fi(t,x1(t),..., xd(t))Ak(t)dt I− ke Fi,k=Zpk pk−1 · xi(t)Ak(t)dt =Zpk pk−1 fi(t,x1(t),..., xd(t))Ak(t)dt and, for all i= 1,..., d and k= 1,..., n , we obtain, by substituting xi(pk) = xi,0 + k−1 ∑ j=1 Fi,jIj+e Fi,kI− k , for h−→ 0, IkFi,k=h f pk,x1,0 + k−1 ∑ j=1 F1,jIj+e F1,kI− k,..., xd,0 + k−1 ∑ j=1 Fd,jIj+e Fd,kI− k!+O(h3)(24) I− ke Fi,k=h 2f pk,x1,0 + k−1 ∑ j=1 F1,jIj+e F1,kI− k,..., xd,0 + k−1 ∑ j=1 Fd,jIj+e Fd,kI− k!+O(h3). (25) Approximated values Gi,k and e Gi,k of Fi,k and e Fi,k , respectively, are computed by solving the systems of dequations IkGi,k=h fi pk,x1,0 + k−1 ∑ j=1 G1,jIj+e G1,kI− k,..., xd,0 + k−1 ∑ j=1 Gd,jIj+e Gd,kI− k!(26) I− ke Gi,k=h 2fi pk,x1,0 + k−1 ∑ j=1 G1,jIj+e G1,kI− k,..., xd,0 + k−1 ∑ j=1 Gd,jIj+e Gd,kI− k!. (27) The first components Gi,1, from the initial condition · xi(t0)=fi(p1,x1,0,..., xd,0), are Gi,1 =h fi(t0,x1,0, ..., xd,0) I1(28) Axioms 2020,9, 15 16 of 37 Problem No2: (Van der Pol equations) Solution interval is t∈[0,300];          . x1=x2 . x2=µ(1−x2 1)x2−x1+asin(ωt) x1(0) = 2 x2(0) = 0. 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 hode45 = 2.196 × 10 −5 and final solution vector x(1) = 1.700147566 ×100,x(2) = −1.819270275 ×10−2. With MFT = 30,001, we choose N = 501 to have hstep = 2.0 × 10 −5 and ODE-FT finds the final value x( 1 ) = 1.700150689 × 10 0 , x( 2 ) = − 1.819263248 × 10 −2 . The comparison gives AveDiffFTRK = 1.398228278980211 ×10−4and MaxDiffFTRK = 0.397462807473403. Figure 3. Problem No2 (Van der Pol): The two components xj(t) , j= 1,2 are displayed in the order from top to bottom; FT solution is red-colored, RK solution is blue-colored. Remark the eight jumps of the solution from positive to negative values or vice versa. Problem No3: (Rössler’s equations) Solution interval is t∈[0, 20];      . x1=−x1−x3x1(0) = 1 . x2=x1+αx2x2(0) = 1 . x3=β+x3(x1−γ)x3(0) = 1 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 hode45 = 0.00144 and final solution vector x( 1 ) = −3.722813228 ×10−1,x(2) = 7.646204266 ×100,x(3) = 3.722813233 ×10−1. Axioms 2020,9, 15 17 of 37 With MFT = 101, we choose N = 41 to have hstep = 0.005 and ODE-FT finds the final value x( 1 ) = − 3.722813228 × 10 −1 , x( 2 ) = 7.646166172 × 10 0 , x( 3 ) = 3.722813233 × 10 −1 . The comparison (see Figure 4) gives AveDiffFTRK = 1.0 ×10−5∗(0.025127112474346,0.947056407667963,0.015362751349660), MaxDiffFTRK = 1.0 ×10−4∗(0.030880110357123, 0.380942339370804,0.039961137723310). Figure 4. Problem No3 (Rössler system):The three components xj(t) , j= 1,2, 3 are displayed in the order from top to bottom; FT solution is red-colored, RK solution is blue-colored. Problem No4: (Lorenz’s system) Solution interval is t∈[0,100];      . x1=a(x2−x1)x1(0) = 10 . x2=x1(b−x3)−x1x2(0) = 20 . x3=x1x2−cx3x3(0) = 10 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 hode45 = 9.37 × 10 −7 and computes the final value x( 1 ) = − 2.480991301578959 × 10 0 , x( 2 ) = 4.324784533758371 × 10 −1 , x( 3 ) = 2.470283374011438 × 10 1 . Taking M = 30,001 the step size hstep = 1.0 × 10 −6 is obtained with N = 1001; the found final value using ODE-FT is x( 1 ) = − 2.453228212514948 × 10 0 , x( 2 ) = 3.798179329138090 × 10 −1 , x( 3 ) = 2.457087763696305 × 10 1 with AveDiffFTRK = (0.004740929429043, 0.006537238822312, 0.008150973471168) and MaxDiffFTRK = (0.149405722393993, 0.285676182819206, 0.328764820957794). See Figure 5for the trajectories of components xj(t) , j= 1,2,3 and Figure 6for the 3D representations of the solutions. Axioms 2020,9, 15 18 of 37 Figure 5. Problem No4 (Lorenz system): The three components xj(t) , j= 1,2, 3 are displayed in the order from top to bottom; FT solution is red-colored, RK solution is blue-colored. Remark that this system, with the given values of parameters, exhibits a strong sensitivity to changes in initial conditions and the numerical solutions are a famous example of chaotic trajectory. Figure 6. Problem No4 (Lorenz system): 3D representations of the solutions obtained by ODE-FT ( left ) and ode45 (right); they appear to be coincident to graphical precision. Axioms 2020,9, 15 19 of 37 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 7and 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 xj(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. Problem No5 : (Periodic system with period T= 8, see [ 24 ], Section 6.8): Solution interval is t∈[0,8];                . x1=x3x1(0) = 1−k . x2=x4x2(0) = 0 . x3=−α2x1 (x2 1+x2 2)3 2x3(0) = 0 . x4=−α2x2 (x2 1+x2 2)3 2x4(0) = α1+k 1−k The parameters are k= 0.25, α=π 4q1+k 1−k ; we remark that from periodicity, the exact solution has xi(8) = xi(0),i=1,...,4. The initial condition is x(1) = 7.50 ×10−1,x(2) = 0.0, x(3) = 0.0, x(4) = 1.013944668993403. With MFT = 801, N = 201 and the step size hstep = 5.0 × 10 −5 , the final solution found by ODE-FT is x( 1 ) = 7.500000000000394 × 10 −1 , x( 2 ) = − 1.890855 × 10 −8 , x( 3 ) = 1.923939 × 10 −8 , x( 4 ) = 1.013944668993378, and the one computed by ode45 is x( 1 ) = 7.499999999999980 × 10 −1 , x( 2 ) = − 5.566185 × 10 −14 , x( 3 ) = 5.123376 × 10 −14 , x( 4 ) = 1.013944668993403; comparatively (see Figure 8), we get AveDiffFTRK = 1.0 ×10−8∗(0.4273611,0.3695634,0.3303905,0.3234898), MaxDiffFTRK = 1.0 ×10−7∗(0.1278607, 0.18908499,0.19239342,0.13127243). Axioms 2020,9, 15 20 of 37 Figure 8. Problem No5 (Periodic system): The four components xj(t) , j= 1,2,3,4 are displayed in the order from top to bottom; FT solution is red-colored, RK solution is blue-colored and the two coincide at graphical precision. 5. 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= (ba ; ea) , B= (bb ; eb) where ba=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 KC. 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∈ KC , the gH-difference is the interval C∈ KC (it always exists and is unique) such that AgH B=C⇐⇒ ((i) A =B+C or (ii) B =A−C. (33) Using midpoint notation, we have AgH B= (ba−bb ; |ea−eb|) and (i) is verified if ea≥eb and C= (ba−bb ; ea−eb) , (ii) is verified when ea≤eb and C= (ba−bb ; eb−ea) . If ea=eb then clearly AgH B= (ba−bb;0) = {ba−bb}is a singleton (real number). Please note that if AigH Bi are gH-differences of the same type for all i= 1,..., n , i.e., all satisfy either (i) or (ii) above, then (see [25]) n ∑ i=1AigH Bi=n ∑ i=1AigH n ∑ i=1Bi. (34) Axioms 2020,9, 15 21 of 37 An interval-valued function F:[a , b]→ KC will be denoted by F(t)=[f−(t) , f+(t)] or, in equivalent midpoint notation, by F(t) = ( b f(t);e f(t)), for t∈[a,b]. We will denote by RF 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 [u]0=cl{x|x∈R,u(x)>0}. 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+ α]=(b uα;e uα) We refer to functions u− (.) and u+ (.) as the lower and upper branches of u , respectively; the two functions b u(.),e u(.)are the (level-wise) midpoint and radius functions. Given fuzzy numbers u , v∈RF , the level-wise generalized Hukuhara difference (LgH-difference for short) is the family of intervals uLgH v=[u]αgH [v]α|α∈[0, 1]. A fuzzy-valued function F:[a , b]→RF will have α -cuts denoted by [F(t)]α= [ f− α(t) , f+ α(t)] or, equivalently, by (b fα(t) ; e fα(t)) , for t∈[a , b] . Remark that for each α∈[ 0,1 ] , the functions [F]α:[a , b]→ KC defined by [F]α(t) = [F(t)]α are properly interval-valued functions and their family (sometimes called their bunch) {[F]α|α∈[0, 1]}gives a unique and equivalent representation of F. 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 t0∈]a,b[and h be such that t0+h∈]a,b[. - If F:[a , b]→ KC is an interval-valued function, its gH-derivative at t0 is defined to be the limit, if it exists, F0 gH(t0) = lim h→0F(t0+h)gH F(t0) h.(35) - If F:]a , b[→RF is a fuzzy-valued function, its level-wise gH-derivative (LgH-derivative for short) at t0 is defined to be the family of the gH-derivatives of [F]α, if they exist for all α∈[0,1], i.e., F0 LgH(t0) = n([F]α)0gH(t0)|α∈[0,1]owhere (36) ([F]α)0gH(t0) = lim h→0 1 h[F(t0+h)]αgH [F(t0)]α.(37) For an interval-valued function F:[a,b]→ KC,F(t)=[f−(t),f+(t)], when f−(t)and f+(t)are both differentiable, we can distinguish two cases, corresponding to (i) and (ii) of (33) (see [18]) Definition 2. Let F :[a,b]−→ KCand t0∈]a,b[with fα(t)and fα(t)both differentiable at t0. We say that - F is (i)-gH-differentiable at t0if (i.) F0 gH(t0) = hf−0(t0),f+0(t0)i(38) - F is (ii)-gH-differentiable at t0if (ii.) F0 gH(t0) = [f+0(t0),f−0(t0)].(39) If F:]a , b[→RF 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]α. Axioms 2020,9, 15 22 of 37 As in [ 18 ], we say that a point t0∈]a , b[ is an l -critical point of F if it is a critical point for the length function len([F(t)]) = f+(t)−f−(t) . A point t0∈]a , b[ is a switching point for the gH-differentiability of F, if in any neighborhood Vof t0there exist points t1<t0<t2such that type-I switch point): at t1 (38) holds while (39) does not hold and at t2 (39) holds and (38) does not hold, or type-II switch point): at t1 (39) holds while (38) does not hold and at t2 (38) holds and (39) does not hold. Analogous definitions can be given level-wise for a fuzzy-valued function. 5.1. Numerical Interval ODE by F-Transform An interval differential equation (IDE) with initial condition can be written in the form (· xgH (t)=F(t,x(t)) ,t∈[t0,t1] x(t0)=x0(40) where F(t,x(t)) = [F−(t,x(t)) , F+(t,x(t))] and x(t)= [x−(t) , x+(t)] are intervals for all t∈[t0,t1] and x0= [x− 0,x+ 0]is an interval initial value. The gH-derivative of x(t)is denoted by · xgH (t)= [ · x− gH (t),· x+ gH (t)]. We will approximate · xgH (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., (· x− gH)(P,A)(t)= n ∑ j=1 F− jAj(t)(41) (· x+ gH)(P,A)(t)= n ∑ j=1 F+ jAj(t)(42) where for j= 1,2,..., n , using the same notation as in Section 2with [a , b]=[t0 , t1] , 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 F− j≤F+ j (because · x− gH (t)≤· x+ gH (t) for all t ) and we can define the intervals Fj= [F− j , F+ j] . Consequently, in terms of interval arithmetic operations, we also have the following approximation of the (interval-valued) inverse F-transform (· xgH)(P,A)(t)of · xgH: (· xgH)(P,A)(t)= n ∑ j=1 FjAj(t). On the other hand, the following function is well defined: H(t) = t Rt0 n ∑ j=1 FjAj(s)ds!= n ∑ j=1 FjBj(t),t∈[t0,t1](43) and H(t0) = 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). Axioms 2020,9, 15 23 of 37 For simplicity, define the following interval-valued function ϕ(t) = n ∑ j=1 FjAj(t),t∈[t0,t1]; (44) Clearly, ϕ is a continuous interval-valued function and, from Theorem 43(i) in [ 18 ], the integral function H(t) = t Rt0 ϕ(s)ds is gH-differentiable with H0 gH(t) = ϕ(t),H(t0) = 0. 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) on [t0,t1]by Φ(t) = H(t)gH (−x0)(45) Ψ(t) = H(t) + x0.(46) Then (1) Ψ(t)is gH-differentiable with Ψ0gH(t) = ∑n j=1FjAj(t)and Ψ(t0) = x0. (2) If all the gH-differences H(t)gH (−x0) are of the same type for all t , then Φ(t) is gH-differentiable with Φ0gH(t) = ∑n j=1FjAj(t)and Φ(t0) = x0. Consequently, if function ϕ:[t0,t1]−→ KC is generated by a fixed fuzzy partition (P,A) as in (44) and H(t) = ∑n j=1FjBj(t), then we have the following two cases: (1) If the intervals Fj,j=1, ..., nare such that ϕ(t)=F(t,H(t) + x0)for t∈[t0,t1] then Ψ(t)defined on [t0,t1]as in (46) is a solution of (40). (2) If the intervals Fj,j=1, ..., nare such that ϕ(t)=ft,H(t)gH (−x0)for t∈[t0,t1] then Φ(t) defined on [t0,t1] as in (45) is a solution of (40), provided that all the gH-differences H(t)gH (−x0)are of the same type for all t∈[t0,t1]. Remark 2. From the properties of gH-difference, we have that Ψ(t)gH x0=H(t) Φ(t)gH x0= (H(t)gH (−x0)) gH x0. It is interesting to observe that function Ψ(t) is (i)-gH-differentiable at all points, while function Φ(t) is (i)-gH-differentiable if the differences H(t)gH (−x0) are of type (i), i.e., if H(t) = Φ(t)−x0 , and is (ii)-gH-differentiable if the differences H(t)gH (−x0) are of type (ii), i.e., if x0=Φ(t)−H(t) . Consequently, solution Ψ(t)= [Ψ−(t),Ψ+(t)]=(b Ψ(t);e Ψ(t)) has always increasing length, but the same is not true for solution Φ(t)= [Φ−(t),Φ+(t)]=(b Φ(t);e Φ(t)) Axioms 2020,9, 15 24 of 37 in the case where Φ(t)is (ii)-gH-differentiable. Finally, let G− j and G+ j be the O(h3) approximations, respectively, of F− j and F+ j similar to (18)–(19); we can see that the interval-valued iF-transform approximation (· xgH)(P,A)(t)=∑n j=1GjAj(t) of · xgH is useful in determining conditions for a switching point. Indeed, if pk−1<pk<pk+1 are three adjacent points of P (2 ≤k≤n− 1) and we suppose that no switching point exists internally to the two subintervals [pk−1 , pk] and [pk , pk+1] (i.e., possibly, the switching is exactly at pk) then the solution x(t)must satisfy x(pk)gH x(pk−1) = HL kand x(pk+1)gH x(pk) = HR k(47) where HL k= pk R pk−1 · xgH (t)dt and HR k= pk+1 R pk · xgH (t)dt. (48) On the other hand, we have HL k=Fk−1I+ k−1+FkI− kand HR k=FkI+ k+Fk+1I− k+1. Suppose now that pkis a switching point for the gH-differentiability of x(t). We have two cases: (I): (i)-to-(ii) switch: x(t) is (i)-gH-differentiable on ]pk−1 , pk[ and is (ii)-gH-differentiable on ]pk , pk+1[ , i.e., x(pk) = x(pk−1) + HL k and x(pk) = x(pk+1)−HR k so that x(pk+1) = (x(pk−1) + HL k)gH (−HR k), where the difference is of type (i); (II): (ii)-to-(i) switch: x(t) is (ii)-gH-differentiable on ]pk−1 , pk[ and is (i)-gH-differentiable on ]pk , pk+1[ , i.e., x(pk−1) = x(pk)−HL k and x(pk+1) = x(pk) + HR k so that x(pk+1) = x(pk−1)gH (−HL k)+HR k, where the difference is of type (i). Instead, if pk is not a switching point and x(pk+1) = x(pk−1) + HL k+HR k we have a type (i) solution on ]pk−1 , pk+1[ ; or, if x(pk−1) = x(pk+1)−HL k+HR k we have a type (ii) solution on ]pk−1 , pk+1[ . 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 [t0,t1] and pk−1<pk<pk+1 (2 ≤k≤ n− 1); let x(t) = (b x(t) ; e x(t)) be a solution to (40) and HL k= ( b HL k , e HL k) and HR k= ( b HR k , e HR k) be given as in (48). Then, the midpoint values satisfy b x(pk) = b x(pk−1) + b HL k b x(pk+1) = b x(pk−1) + b HL k+b HR k. Assuming that only pkis eventually a switching point, we have the following four cases (a) if pkis a (i)-to-(ii) switch, then e x(pk) = e x(pk−1) + e HL k e x(pk+1) = e x(pk−1) + e HL k−e HR k≥0 (b) if pkis a (ii)-to-(i) switch, then e x(pk) = e x(pk−1)−e HL k≥0 e x(pk+1) = e x(pk−1)−e HL k+e HR k≥0 (c) if x(t)is (i)-gHdifferentiable on [pk−1,pk+1], then e x(pk) = e x(pk−1) + e HL k e x(pk+1) = e x(pk−1) + e HL k+e HR k Axioms 2020,9, 15 25 of 37 (d) if x(t)is (ii)-gHdifferentiable on [pk−1,pk+1], then e x(pk) = e x(pk−1)−e HL k≥0 e x(pk+1) = e x(pk−1)−(e HL k+e HR k)≥0. There are no general rules to locate a switching point. Denote by x(t) = [x−(t) , x+(t)] = (b x(t) ; e x(t)) a solution of the IDE (40); if x(t) is (i)-gH-differentiable, its length e x(t) will not decrease, while e 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)≤Uand requiring that l(t)≤e x(t)≤u(t)for all t; then, (a) a (i)-to-(ii) switch is decided at t=pkif the following condition is reached (e x(pk) = e x(pk−1) + e HL k≤u(pk) e x(pk+1) = e x(pk−1) + e HL k+e HR k>u(pk+1) (b) a (ii)-to-(i) switch is decided at t=pkif the following condition is reached (e x(pk) = e x(pk−1)−e HL k≥l(pk) e x(pk+1) = e x(pk−1)−e HL k−e HR k<l(pk+1). 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 x0 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 e x0 is such that L<e x0<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 e x(t) and trying to intercept points t∗ where the function e 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)] = (b x(t) ; e x(t)) is such that x−(t) and x+(t) are differentiable so that its gH-derivative · xgH(t) can be expressed in terms of the derivatives (x−)0(t) and (x+)0(t) . Let (P , A) be a fuzzy partition of [t0 , t1] and let (· xgH)(P,A)(t)= ∑n k=1FkAk(t) be the iF-transform of · xgH(t) with interval-valued components Fk= (b Fk ; e Fk) . If pk is a local minimum or maximum of e x(t), then e Fk=0+O(h2). Proof. From the properties of F-transform, we have F− k= ( · xgH)−(pk) + O(h2) , F+ k= ( · xgH)+(pk) + O(h2) and, in particular, e Fk=g · xgH(pk) + O(h2) . On the other hand, from the differentiability of x−(t) and x+(t) , we have 0 =d dt e x(t)) for t=pk , i.e., (x+)0(pk)−(x−)0(pk) = 0. It follows that g · xgH(pk) = (x+)0(pk)−(x−)0(pk) 2=0, i.e., e Fk=g · xgH(pk) + O(h2) = 0+O(h2). We then suggest the following (purely endogenous) switching rule. Axioms 2020,9, 15 32 of 37 Table 3. (Interval-valued problem IDE2-a): For five values of α∈{0,0.2,0.4,0.6,0.8} , the three switching points tw corresponding to Meth=1 are given in the first column; the second column contains the interval-valued solution X(tw) = [x− α(tw) , x+ α(tw)] . If α= 1, the solution is single-valued and no switching point exists. α= 0.00 Initial value = [−1.0, 1.0] tw = 1.7954×10−1X(tw) = [−1.0090673275×10−1, 1.0090673275×10−1] tw = 2.5918×100X(tw) = [−5.0989773431×10−1, 5.0989773431×10−1] tw = 4.3297×100X(tw) = [−8.9894395242×10−1, 8.9894395242×10−1] Final value: [−8.35571783×10−1, 8.35571783×10−1] α= 0.20 Initial value: [−8.00×10−1, 8.00×10−1] tw = 2.7944×10−1X(tw) = [−8.1796697089×10−1, 8.1796697089×10−1] tw = 2.5810×100X(tw) = [−4.6779933026×10−1, 4.6779933026×10−1] tw = 4.3542×100X(tw) = [−7.9775755945×10−1, 7.9775755945×10−1] Final value: [−7.52396593×10−1, 7.52396593×10−1] α= 0.40 Initial value: [−6.00×10−1, 6.00×10−1] tw = 3.7935×10−1X(tw) = [−6.2502008124×10−1, 6.2502008124×10−1] tw = 2.5579×100X(tw) = [−4.0293011585×10−1, 4.0293011585×10−1] tw = 4.3844×100X(tw) = [−6.7059314969×10−1, 6.7059314969×10−1] Final value: [−6.41450486×10−1, 6.41450486×10−1] α= 0.60 Initial value: [−4.00×10−1, 4.00×10−1] tw = 4.8538×10−1X(tw) = [−4.2690473012×10−1, 4.2690473012×10−1] tw = 2.5174×100X(tw) = [−3.0871431900×10−1, 3.0871431900×10−1] tw = 4.4226×100X(tw) = [−5.0673707018×10−1, 5.0673707018×10−1] Final value: [−4.91260393×10−1, 4.91260393×10−1] α= 0.80 Initial value: [−2.00×10−1, 2.00×10−1] tw = 6.0271×10−1X(tw) = [−2.2003436269×10−1, 2.2003436269×10−1] tw = 2.4500×100X(tw) = [−1.7731361249×10−1, 1.7731361249×10−1] tw = 4.4716×100X(tw) = [−2.9073091262×10−1, 2.9073091262×10−1] Final value: [−2.85325600×10−1, 2.85325600×10−1] α= 1.00 Initial value: [ 0.00, 0.00] Final value: [0.00000000×100, 0.00000000×100] Table 4. (Interval-valued problem IDE2-b): For five values of α∈{0,0.2,0.4,0.6,0.8} , the switching points tw corresponding to Meth=2 are given in the first column; the second column contains the interval-valued solution X(tw) = [x− α(tw) , x+ α(tw)] . The number of switching points changes with α and if α=1, the solution is single-valued and no switching point exists. α= 0.00 Initial value = [−1.0, 1.0] tw = 1.9321×10−1X(tw) = [−9.8998479170×10−1, 9.8998479170×10−1] Final value: [−3.50656260×100, 3.50656260×100] α= 0.20 Initial value: [−8.00×10−1, 8.00×10−1] tw = 3.0866×10−1X(tw) = [−7.7947575655×10−1, 7.7947575655×10−1] Final value: [−1.87166877×100, 1.87166877×100] α= 0.40 Initial value: [−6.00×10−1, 6.00×10−1] tw = 4.2600×10−1X(tw) = [−5.7091483288×10−1, 5.7091483288×10−1] tw = 3.0348×100X(tw) = [−9.5345951680×10−1, 9.5345951680×10−1] tw = 4.2887×100X(tw) = [−8.4078653821×10−1, 8.4078653821×10−1] Final value: [−8.90550965×10−1, 8.90550965×10−1] α= 0.60 Initial value: [−4.00×10−1, 4.00×10−1] tw = 5.4899×10−1X(tw) = [−3.6865624377×10−1, 3.6865624377×10−1] tw = 2.7506×100X(tw) = [−5.2242232939×10−1, 5.2242232939×10−1] tw = 4.5230×100X(tw) = [−3.3741623443×10−1, 3.3741623443×10−1] Final value: [−3.44451472×10−1, 3.44451472×10−1] α= 0.80 Initial value: [−2.00×10−1, 2.00×10−1] tw = 6.8094×10−1X(tw) = [−1.7696972113×10−1, 1.7696972113×10−1] tw = 2.5277×100X(tw) = [−2.1957335579×10−1, 2.1957335579×10−1] tw = 4.6445×100X(tw) = [−8.1718701678e−02, 8.1718701678e−02] Final value: [−8.21485320e−02, 8.21485320e−02] α= 1.00 Initial value: [0.00, 0.00] Final value: [0.00000000×100, 0.00000000×100] Axioms 2020,9, 15 33 of 37 Problem FDE1: Solution interval is t∈[0, 1 2]; (. xgH(t) = −1 2x(t) + 2sin(3t) x(0) = x0 where [x0]α= [−1+α,1 −α],α∈[0,1]. The step size in this case is h= 4.000 ×10−5. The solution found with both Meth = 1 and Meth = 2 are fuzzy-valued with lengths of the α -cuts all increasing (Meth = 1, Figure 13) or all decreasing (Meth = 2, Figure 14) and there are no switching points. The α -cuts of the initial condition x0 and the final fuzzy solution for Meth = 1 and Meth = 2 are inserted in Table 5. Table 5. (Fuzzy-valued problem FDE1): For eleven values of α∈ni−1 10 |i=1,...,10o , the table contains the interval level-wise initial condition (column 2), the final interval value with Meth = 1 (column 3) and the final interval value with Meth = 2. αInitial Condition Final solution (Meth = 1) Final solution (Meth = 2) 0.0 [−1.0000×100, 1.0000×100] [−7.9066×100, 6.8715×100] [−6.5292×10−1,−3.8225×10−1] 0.1 [−9.0000×10−1, 9.0000×10−1] [−7.1677×100, 6.1326×100] [−6.3939×10−1,−3.9579×10−1] 0.2 [−8.0000×10−1, 8.0000×10−1] [−6.4288×100, 5.3937×100] [−6.2586×10−1,−4.0932×10−1] 0.3 [−7.0000×10−1, 7.0000×10−1] [−5.6899×100, 4.6548×100] [−6.1232×10−1,−4.2285×10−1] 0.4 [−6.0000×10−1, 6.0000×10−1] [−4.9510×100, 3.9158×100] [−5.9879×10−1,−4.3639×10−1] 0.5 [−5.0000×10−1, 5.0000×10−1] [−4.2121×100, 3.1769×100] [−5.8526×10−1,−4.4992×10−1] 0.6 [−4.0000×10−1, 4.0000×10−1] [−3.4732×100, 2.4380×100] [−5.7172×10−1,−4.6345×10−1] 0.7 [−3.0000×10−1, 3.0000×10−1] [−2.7343×100, 1.6991×100] [−5.5819×10−1,−4.7699×10−1] 0.8 [−2.0000×10−1, 2.0000×10−1] [−1.9954×100, 9.6022×10−1] [−5.4465×10−1,−4.9052×10−1] 0.9 [−1.0000×10−1, 1.0000×10−1] [−1.2565×100, 2.2132×10−1] [−5.3112×10−1,−5.0405×10−1] 1.0 [0.0000×100, 0.0000×100] [−5.1759×10−1,−5.1759×10−1] [−5.1759×10−1,−5.1759×10−1] Figure 13. Problem FDE1, Meth = 1: Fuzzy-valued gH-differentiable solution ( left ) and its gHderivative (right). There are no switching points. Figure 14. Problem FDE1, Meth = 2: Fuzzy-valued gH-differentiable solution ( left ) and its gH-derivative (right). For all α-cuts, there are no switching points. Axioms 2020,9, 15 34 of 37 Problem FDE2: Solution interval is t∈[0, 4π]; (. xgH(t) = sin(t)x(t) x(0) = x0 with [x0]α= [−1+α,1 −α],α∈[0,1]. In this case, the step size is h= 1.257 ×10−4. 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 . xgH 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). Figure 15. Problem FDE2, Meth = 1: gH-differentiable solution ( left ) and its gH-derivative ( right ). There are three switching points, in the same position for all α-cuts. Figure 16. Problem FDE2, Meth = 2: gH-differentiable solution ( left ) and its gH-derivative ( right ). There are three switching points, in the same position for all α-cuts. 7. 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. One of our conclusions is then that developing numerical procedures based on F-transform is a promising area of research, anticipated by some successes in recent research such as [ 13 – 16 ]; a complete comparison of (and between) the various F-transform-based proposed methods and our approach was not a scope of our study, where we have chosen standard well-performing routines as benchmarks and we have evaluated algorithm ODE-FT with the same and sufficiently small step size h on typical (including hard) ODEs. Two of the examples in Section 4are also presented in [ 14 ], Axioms 2020,9, 15 35 of 37 where the quantity MSE (mean squared error of approximate solution and the exact one) is computed. Ex4 is Example 1 in [ 14 ] and Ex5 is Example 3 in [ 14 ]. The best MSE quantities obtained by [ 14 ] and by ODE-FT are reported in Table 6: Table 6. For the two ODE problems in example Ex4 and Ex5, using different step-sizes, the table contains the computed MSE for the two solution-variables x1(t) and x2(t) , obtained by ODE-FT and Scheme II in [14]. Example/(Algorithm, Step Size) MSE(x1)MSE(x2) Ex4/(Scheme II in [14], h = 0.01) 2.241 ×10−23.775 ×10−4 Ex4/(ODE-FT, h = 0.01) 1.865 ×10−62.027 ×10−6 Ex5/(Scheme II in [14], h = 0.1) 1.721 ×10−54.102 ×10−5 Ex5/(ODE-FT, h = 0.1) 1.242 ×10−61.685 ×10−5 Ex5/(ODE-FT, h = 0.01) 1.264 ×10−10 2.693 ×10−9 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∈[t0 , t1] and not only at the discretized points pk , as usual in the (explicit) singleor 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 Fd -transform (see [ 11 ] for recent results of its properties), by introducing local polynomials or, more generally, local parametric functions Fk(t;θ) in place of constant direct components Fk (coefficients of the polynomials or parameters θ∈Rd are then estimated by least squares). In these cases, the inverse F-transform of a function f(t) on [a , b] has the form fd (P,A)(t) = n ∑ k=1Fk(t ; θ(k))Ak(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 fd (P,A)(t), i.e., b Ra f(t)dt = b Ra fd (P,A)(t)dt =n ∑ k=1 b Ra Fk(t;θ(k))Ak(t)dt. (49) It should be interesting to see if higher-order F-transform approximations will be able to generate high orders O(hq) , q> 2 of convergence, and to obtain possibly increasing orders q by increasing d (two numeric schemes of order q=2 based on the F2-transform are obtained in [16]). Similar results can be obtained by considering the discrete F-transform f(P,A)(tj) on a data set of points S=(tj,fj),j=1, 2,..., m ; the integrals are substituted by summations and we have (see [ 10 ]) m ∑ j=1f(P,A)(tj) = m ∑ j=1fj. (50) 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 Axioms 2020,9, 15 36 of 37 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 [t0,t1]is expressed in the form          Ft,x(t),· x(t)=0 G(t,x(t))=0 H(t,x(t))≤0 x(t0)=x0. (51) 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,x0)(t) at points pk∈P will transform the DAE into a standard feasibility problem, consisting of finding feasible solutions for the transformed system at points pk,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. Acknowledgments: Davide Radi gratefully acknowledges the support of the Czech Science Foundation (GACR) under project [20-16701S] and the VŠB-TU Ostrava under the SGS project SP2020/11. Conflicts of Interest: The authors declare no conflict of interest. References 1. Perfilieva, I. Fuzzy Transforms: Theory and Applications. Fuzzy Sets Syst. 2006,157, 993–1023. [CrossRef] 2. Perfilieva, I. Fuzzy Transforms: A challenge to conventional transform. In Advances in Images and Electron Physics; Hawkes, P.W., Ed.; Elsevier Academic Press: Cambridge, CA, USA, 2007; Volume 147, pp. 137–196. 3. Coroianu, L.; Stefanini, L. General approximation of fuzzy numbers by F-transform. Fuzzy Sets Syst. 2016 , 288, 46–74. [CrossRef] 4. Coroianu, L.; Stefanini, L. Properties of fuzzy transform obtained from Lp minimization and a connection with Zadeh’s extension principle. Inf. Sci. 2019,478, 331–354. [CrossRef] 5. Guerra, M.L.; Stefanini, L. Quantile and expectile smoothing based on L1 -norm and L2 -norm fuzzy transforms. Int. J. Approx. Reason. 2019,107, 17–43. [CrossRef] 6. Perfilieva, I.; de Baets, B. Fuzzy transforms of monotone functions with applications to image compression. Inf. Sci. 2010,180, 3304–3315. [CrossRef] 7. Perfilieva, I.; Novak, V.; Dvorak, A. Fuzzy transform in the analysis of data. Int. J. Approx. Reason. 2008 ,48, 36–46. [CrossRef] 8. Stefanini, L. Fuzzy Transform with Parametric LU-Fuzzy Partitions. In Computational Intelligence in Decision and Control; Ruan et Al, D., Ed.; World Scientific: Singapore, 2008; pp. 399–404. 9. Stefanini, L. Fuzzy Transform and Smooth Functions. In Proceedings of the IFSA-EUSFLAT 2009 Conference, Lisbon, Portugal, 20–24 July 2009; pp. 579–584. 10. Stefanini, L. F-Transform with Parametric Generalized Fuzzy Partitions. Fuzzy Sets Syst. 2011 ,180, 98–120. [CrossRef] 11. Zeinali, M.; Alikhani, R.; Shahmorad, S.; Bahrani, F.; Perfilieva, I. On the structural properties of Fm -transform with applications. Fuzzy Sets Syst. 2018,342, 32–52. [CrossRef] 12. Ahmad, M.Z.; Hasan, M.K.; de Baets, B. Analytical and numerical solutions of fuzzy differential equations. Inf. Sci. 2013,236, 156–167. [CrossRef] 13. Kasasbeh, H.A.L.; Perfilieva, I.; Ahmad, M.Z.; Yahya, Z.R. New Fuzzy Numerical Methods for Solving Cauchy Problems. Appl. Syst. Innov. 2018,1, 15. [CrossRef] 14. Kasasbeh, H.A.L.; Perfilieva, I.; Ahmad, M.Z.; Yahya, Z.R. New Approximation Methods Based on Fuzzy Transform for solving SODEs: I. Appl. Syst. Innov. 2018,1, 29. [CrossRef] 15. Kasasbeh, H.A.L.; Perfilieva, I.; Ahmad, M.Z.; Yahya, Z.R. New Approximation Methods Based on Fuzzy Transform for solving SODEs: II. Appl. Syst. Innov. 2018,1, 30. [CrossRef] Axioms 2020,9, 15 37 of 37 16. Khastan, A.; Perfilieva, I.; Alijani, Z. A new fuzzy approximation method to Cauchy problems by fuzzy transform. Fuzzy Sets Syst. 2016,288, 75–95. [CrossRef] 17. Radi, D.; Stefanini, L. Fuzzy Differential Equations by F-Transform. In Proceedings of the NAFIPS 2015 Annual Meeting, Redmond, WA, USA, 17–19 August 2015; pp. 384–389. 18. Bede, B.; Stefanini, L. Generalized differentiability of fuzzy-valued functions. Fuzzy Sets Syst. 2013 ,230, 119–141; doi:10.1016/j.fss.2012.10.003. [CrossRef] 19. Stefanini, L.; Bede, B. Generalized Hukuhara differentiability of interval-valued functions and interval differential equations. Nonlinear Anal. 2009,71, 1311–1328. [CrossRef] 20. Stefanini, L.; Bede, B. Generalized fuzzy differentiability with LU-parametric representation. Fuzzy Sets Syst. 2014,257, 184–203. [CrossRef] 21. Bede, B. Mathematics of Fuzzy Sets and Fuzzy Logic; Springer: Berlin, Germany, 2013. 22. Gomes, L.T.; de Barros, L.C.; Bede, B. Fuzzy Differential Equations in Various Approaches; Springer: Berlin, Germany, 2015. 23. Stefanini, L. A generalization of Hukuhara difference and division for interval and fuzzy arithmetic. Fuzzy Sets Syst. 2010,161, 1564–1584. [CrossRef] 24. Forsythe, G.E.; Malcolm, M.A.; Moler, C.B. Computer Methods for Mathematical Computations; Prentice-Hall: Englewood Cliffs, NJ, USA, 1977. 25. Stefanini, L.; Arana-Jimenez, M. Karush-Kuhn-Tucker conditions for interval and fuzzy optimization in several variables under total and directional generalized differentiability. Fuzzy Sets Syst. 2019 ,362, 1–34. [CrossRef] 26. Stefanini, L.; Guerra, M.L.; Amicizia, B. Interval Analysis and Calculus for Interval-Valued Functions of a Single Variable. Part I: Partial Orders, gH-derivative, Momotonicity. Axioms 2019,8, 113. [CrossRef] 27. Stefanini, L.; Sorini, L.; Amicizia, B. Interval Analysis and Calculus for Interval-Valued Functions of a Single Variable. Part II: Extremal Points, Convexity, Periodicity. Axioms 2019,8, 114. [CrossRef] 28. Tomasiello, S. An alternative use of fuzzy transform with application to a class of delay differential equations. Int. J. Comput. Math. 2017,94, 1719–1726. [CrossRef] 29. Stamova, I.M.; Stamov, G.T. Functional and Impulsive Differential Equations of Fractional Order—Qualitative Analysis and Applications; CRC Press: Boca Raton, FL, USA, 2017. 30. Ascher, U.M.; Petzold, L.R. Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations; SIAM: Philadelphia, PA, USA, 1998. 31. Kunkel, P.; Mehrmann, V. Differential-Algebraic Equations—Analysis and Numerical Solution; European Mathematical Society: Zürich, Switzerland, 2006. c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).