# Euler’s Numerical Method on Fractional DSEK Model under ABC Derivative

Authors: Fareeha Sami Khan; M. Khalid; Omar Bazighifan; A. El-Mesady
Journal: Complexity
 (2022)

---

## Abstract

In this paper, DSEK model with fractional derivatives of the Atangana-Baleanu Caputo (ABC) is proposed. This paper gives a brief overview of the ABC fractional derivative and its attributes. Fixed point theory has been used to establish the uniqueness and existence of solutions for the fractional DSEK model. According to this theory, we will define two operators based on Lipschitzian and prove that they are contraction mapping and relatively compact. Ulam-Hyers stability theorem is implemented to prove the fractional DSEK model’s stability in Banach space. Also, fractional Euler’s numerical method is derived for initial value problems with ABC fractional derivative and implemented on fractional DSEK model. The symmetric properties contribute to determining the appropriate method for finding the correct solution to fractional differential equations. The numerical solutions generated using fractional Euler’s method have been plotted for different values ofα where α∈0,1 and different step sizes h. Result discussion will be given, describing the changes that occur due to the step size h.


---

## Body

## 1. Introduction

Fractional calculus is as historic as integer calculus but not until 1819, it was properly introduced in the form of definitions and functions. Many scientists, researchers, and mathematicians played their role in developing its theory such as [1–5]. Since fractional calculus was developed theoretically at first and had no practical application at the time, it was not as well known as integer calculus among other areas of science. However, after the contribution of Professor Mandelbrot’s fractal theory, fractional calculus theory developed rapidly and soon became the hot topic among all researchers around the globe.The existing theory of nonlinear science is now seemed to be only focused on fractional order calculus theory and the theory of chaos and dissipative structure ([6, 7]). The fact that fractional calculus describes the heredity and memory of any physical phenomenon is fascinating [8, 9]. As a result, it is now used more than integer calculus in fluid dynamics, quantum mechanics, mathematical biology applications, chemistry, control and signal theory, economics, image processing, etc. Models developed in many areas of science and engineering are observed to be best explained by fractional differential equations. The symmetries can be found by solving a related set of partial fractional differential equations. Since integer-order models lack memory and heredity, they cannot adequately and sufficiently describe physical phenomena in many cases. These applications have also led to the rapid development of fractional calculus theory. The authors of [10–13] have a great deal of literature on the subject describing applications and types of fractional derivatives. It is critical to note that every one of those fractional derivative order definitions has its own advantages and disadvantages.Aside from the mathematical satisfactions of the fractional-order Atangana-Baleanu derivative, the new derivative is being studied due to the necessity of implementing a model depicting the behavior of orthodox viscoelastic materials, thermal medium, and other materials. The proposed mechanism can depict material heterogeneities as well as some structure or media at multiple scales.The new kernel’s nonlocality enables the full description of memory inside of structure and media with multiple scales, which cannot be represented by classical fractional derivatives or those of the Caputo-Fabrizio type. Furthermore, we believe that Atangana-Baleanu derivatives can play an important role in the study of the microstructural behavior of some materials, particularly those involving nonlocal exchanges, which are important in defining the material’s properties states [14]. Atangana-Baleanu derivatives are thus extremely useful in describing a wide range of scientific, engineering, and technological problems.Descemet’s stripping endothelial keratoplasty (DSEK) is the name given to eye surgery [15–17] in which a damaged corneal layer is replaced with a healthy corneal layer from a donor or synthetic cornea. Cornea is a clear layer of the eye that is very important in the anterior part of the eye; if it is scratched or damaged, it affects vision. It itself is made up of five layers, see Figure 1. Keeping this in mind, the authors of [18, 19] created the DSEK model, which predicts the behavior of ocular parameters posttreatment. Because this procedure has never been studied mathematically as an ordinary system of differential equations, this work is extremely important. Since fractional calculus has been said to be the generalization of integer calculus, the fractional DSEK model is developed and studied theoretically and numerically in this paper.Figure 1
Five layers of the cornea.In this paper, the first section gives an overview of the literature background of fractional calculus and the DSEK model. The second section gives the preliminary concepts of fractional calculus that will be used in this work. Section three is based on the explanation of the fractional DSEK model. Section four shows the existence of a solution by the fixed point theory of the fractional DSEK model. Section five describes the Ulam-Hyers stability analysis of fractional DSEK. Section6 describes the computation of fractional Euler’s method for ABC fractional derivative and the application of fractional Euler’s method to fractional DSEK. The last section is the discussion of the results obtained and the conclusion of this paper.


## 2. Preliminaries of Fractional Calculus

### 2.1. Atangana-Baleanu Caputo Fractional Derivative

Definition 1.
The authors of [20] introduced a new Caputo fractional derivative as(1)Dtαgt=Nα1−α∫ctg′ye−α/1−αt−ydy,where g∈H1b,c,c>b,α∈0,1,Nα is the normalization function that follows the condition N0=N1=1.Definition 2.
If the function does not follow the conditionN0=N1=1, then it takes the form as(2)Dtαgt=αNα1−α∫ctgt−gye−α/1−αt−ydy.
This equation can also take the form of the conditionN0=N∞=1(3)Dtρgt=Nρρ∫ctg′ye−t−y/ρdy,where ρ=1−α/α∈0,∞, also α=1/1+ρ∈0,1.
This derivative was defined by [20] to involve an exponential kernel in fractional derivatives to represent the results of dynamic systems memory effects more accurately. With the passage of time, it occurred that this definition has a flaw in that it does not give the original function when α=1. To overcome this problem, the authors of [21] presented the accurate kernel and modified this definition accordingly.Definition 3.
Let the new fractional derivative be defined as(4)DcABCtαgt=Nα1−α∫ctg′yEαα/α−1t−yαdy,where g∈H1b,c,c>b also α∈0,1 and Nα has the same properties defined in [20]. Here, Eα is the generalized Mittag-Leffler function defined as Eα=Eα−tα=∑k=0∞−tαk/Гαk+1.
For the above definition, the constant function has a fractional derivative of zero. The above description would be helpful when solving real-world issues and will also provide a great benefit in utilizing the Laplace transform to solve any initial state physical problem. Nevertheless, if alpha is 0, we will not recover the initial function except when the function vanishes at the origin. We suggest the following definition, in order to avoid this problem.Definition 4.
Let the new fractional derivative be defined as(5)DcABRtαgt=Nα1−αddt∫ctgyEα−α/1−αt−yαdy,where g∈H1b,c,c>b also α∈0,1 and Nα has the same properties defined in [20]. Here, Eα is the generalized Mittag-Leffler function defined as Eα=Eα−tα=∑k=0∞−tαk/Гαk+1.
Both definitions have a nonlocal kernel. For calculations in this paper, we will use definitions in (4) and (5).


### 2.2. Properties of Atangana-Baleanu Caputo Fractional Derivative

(i)
Laplace transformation on equation (4)is(6)ℒ0ABRDtαgts=Nα1−αsαℒgtssα+α/1−α.(ii)
Laplace transformation on equation (5)is(7)ℒ0ABCDtαgts=Nα1−αsαℒgts−sα−1g0sα+α/1−α.(iii)
Letg∈H1b,c,c>b,α∈0,1, then the following relation exists [20]:(8)D0ABRtαgt=D0ABRtαgt+Ht.(iv)
Ifg is a continuous function on some closed interval a,b. Then, the following inequality can be written on a,b(9)D0ABRtαgt<Nα1−αK,ht=maxb≤t≤cht.(v)
Lipschitz condition Atangana-Baleanu Caputo fractional derivative satisfies the Lipschitz condition in Riemann and Caputo sense, and the following inequality exists:(10)D0ABCtαgt−D0ABCtαft≤Hgt−ft.Similarly, for (5), the Lipschitz condition exists as(11)0ABRDtαgt−0ABRDtαft≤Hgt−ft.(vi)
AB fractional integral for α∈0,1 the AB fractional integral for gt∈H10,t and nonlocal kernel is given as(12)IbABtαgt=1−αNαgt+αNαΓα∫btt−yα−1gydy,T>0.Whenα=1, the ordinary integral is obtained, and for α=0, the initial function is obtained.For proof of these, see [20].Lemma 1.
[21] Suggests that the proposed problem for α∈0,1 has a solution; that is,(13)DABCtαgt=ηtg0=g0,its solution is given by gt=g0+1−α/Nαηt+α/NαΓα∫0tt−yα−1gydy.


## 2.1. Atangana-Baleanu Caputo Fractional Derivative

Definition 1.
The authors of [20] introduced a new Caputo fractional derivative as(1)Dtαgt=Nα1−α∫ctg′ye−α/1−αt−ydy,where g∈H1b,c,c>b,α∈0,1,Nα is the normalization function that follows the condition N0=N1=1.Definition 2.
If the function does not follow the conditionN0=N1=1, then it takes the form as(2)Dtαgt=αNα1−α∫ctgt−gye−α/1−αt−ydy.
This equation can also take the form of the conditionN0=N∞=1(3)Dtρgt=Nρρ∫ctg′ye−t−y/ρdy,where ρ=1−α/α∈0,∞, also α=1/1+ρ∈0,1.
This derivative was defined by [20] to involve an exponential kernel in fractional derivatives to represent the results of dynamic systems memory effects more accurately. With the passage of time, it occurred that this definition has a flaw in that it does not give the original function when α=1. To overcome this problem, the authors of [21] presented the accurate kernel and modified this definition accordingly.Definition 3.
Let the new fractional derivative be defined as(4)DcABCtαgt=Nα1−α∫ctg′yEαα/α−1t−yαdy,where g∈H1b,c,c>b also α∈0,1 and Nα has the same properties defined in [20]. Here, Eα is the generalized Mittag-Leffler function defined as Eα=Eα−tα=∑k=0∞−tαk/Гαk+1.
For the above definition, the constant function has a fractional derivative of zero. The above description would be helpful when solving real-world issues and will also provide a great benefit in utilizing the Laplace transform to solve any initial state physical problem. Nevertheless, if alpha is 0, we will not recover the initial function except when the function vanishes at the origin. We suggest the following definition, in order to avoid this problem.Definition 4.
Let the new fractional derivative be defined as(5)DcABRtαgt=Nα1−αddt∫ctgyEα−α/1−αt−yαdy,where g∈H1b,c,c>b also α∈0,1 and Nα has the same properties defined in [20]. Here, Eα is the generalized Mittag-Leffler function defined as Eα=Eα−tα=∑k=0∞−tαk/Гαk+1.
Both definitions have a nonlocal kernel. For calculations in this paper, we will use definitions in (4) and (5).


## 2.2. Properties of Atangana-Baleanu Caputo Fractional Derivative

(i)
Laplace transformation on equation (4)is(6)ℒ0ABRDtαgts=Nα1−αsαℒgtssα+α/1−α.(ii)
Laplace transformation on equation (5)is(7)ℒ0ABCDtαgts=Nα1−αsαℒgts−sα−1g0sα+α/1−α.(iii)
Letg∈H1b,c,c>b,α∈0,1, then the following relation exists [20]:(8)D0ABRtαgt=D0ABRtαgt+Ht.(iv)
Ifg is a continuous function on some closed interval a,b. Then, the following inequality can be written on a,b(9)D0ABRtαgt<Nα1−αK,ht=maxb≤t≤cht.(v)
Lipschitz condition Atangana-Baleanu Caputo fractional derivative satisfies the Lipschitz condition in Riemann and Caputo sense, and the following inequality exists:(10)D0ABCtαgt−D0ABCtαft≤Hgt−ft.Similarly, for (5), the Lipschitz condition exists as(11)0ABRDtαgt−0ABRDtαft≤Hgt−ft.(vi)
AB fractional integral for α∈0,1 the AB fractional integral for gt∈H10,t and nonlocal kernel is given as(12)IbABtαgt=1−αNαgt+αNαΓα∫btt−yα−1gydy,T>0.Whenα=1, the ordinary integral is obtained, and for α=0, the initial function is obtained.For proof of these, see [20].Lemma 1.
[21] Suggests that the proposed problem for α∈0,1 has a solution; that is,(13)DABCtαgt=ηtg0=g0,its solution is given by gt=g0+1−α/Nαηt+α/NαΓα∫0tt−yα−1gydy.


## 3. Fractional DSEK Model

As mentioned in Section2, the definitions in (4) and (5) have nonlocal kernels, and therefore, Atangana-Baleanu Caputo fractional derivative operator’s performance in modeling eye surgery is better than any other definition. It inspired the valuable applications of several fractional operators in dynamic mathematical models; therefore, we are researching the dynamics of eye surgery derived in [18] by a system of nonlinear differential equations by involving ABC fractional derivative.(14)DABCtαpt=artst+γqt+δ,DABCtαqt=−δ−βrtst−γqt−stqt,DABCtαrt=−βrtst−qt−stpt,DABCtαst=stpt+stqt+qt,with initial conditions(15)p0=p0≥0,q0=q0≥0,r0=r0≥0,s0=s0≥0,where ABCDtα is the Atangana-Baleanu Caputo fractional derivative of order α. DSEK model is based on the same conditions given by [18]. Also, the defined parameters have the same description as given by [18]. Such as pt is the refractive index, qt is the axial length, rt is the corneal curvature, and st is the central corneal thickness.


### 3.1. Preliminaries for Fractional DSEK Model

For fractional analysis of the DSEK model, let us defineξ=p,q,r,s. To define the Banach space, let us say we have B=0,t where 0≤T≤t<∞. Then, the field can be written as G=CB,R4 under the norm supremum as(16)ξ=supT∈BξT:ξ∈G,where ξ=p+q+r+s. Also, p,q,r,s,N∈C0,t.Definition 5.
LetB be a Banach space. Then, ψ defined as B⟶B will be a Lipschitzian if there exists a constant l>0 for which the inequality exists such that(17)ψξ1−ψξ2≤lξ1−ξ2,for all ξ1,ξ2∈B. Where l is the Lipschitz constant for ψ. If l<1, then ψ is a contraction.Theorem 1.
Let B be a Banach space andψ:B⟶B be a contraction mapping. Then, there must exist a unique fixed point of ψ.Theorem 2.
A subset of Banach space B is supposed to beN. Let N be convex, closed, and nonempty. Suppose that F and G map N into G, and the following relations exist:(i)
Fu+Gv∈N∀ξ1,ξ2∈N(ii)
F is continuous and compact(iii)
G is a contraction mapping
Then, there existsξ∈N such that Fξ+Gξ=ξ.


## 3.1. Preliminaries for Fractional DSEK Model

For fractional analysis of the DSEK model, let us defineξ=p,q,r,s. To define the Banach space, let us say we have B=0,t where 0≤T≤t<∞. Then, the field can be written as G=CB,R4 under the norm supremum as(16)ξ=supT∈BξT:ξ∈G,where ξ=p+q+r+s. Also, p,q,r,s,N∈C0,t.Definition 5.
LetB be a Banach space. Then, ψ defined as B⟶B will be a Lipschitzian if there exists a constant l>0 for which the inequality exists such that(17)ψξ1−ψξ2≤lξ1−ξ2,for all ξ1,ξ2∈B. Where l is the Lipschitz constant for ψ. If l<1, then ψ is a contraction.Theorem 1.
Let B be a Banach space andψ:B⟶B be a contraction mapping. Then, there must exist a unique fixed point of ψ.Theorem 2.
A subset of Banach space B is supposed to beN. Let N be convex, closed, and nonempty. Suppose that F and G map N into G, and the following relations exist:(i)
Fu+Gv∈N∀ξ1,ξ2∈N(ii)
F is continuous and compact(iii)
G is a contraction mapping
Then, there existsξ∈N such that Fξ+Gξ=ξ.


## 4. Existence of Solutions for Fractional DSEK Model

By using the fixed-point theory, let us prove the uniqueness and existence of the DSEK model. To prove its uniqueness and existence, let us reformulate the DSEK model of (14).(18)ABCDtαpt=G1t,p,q,r,s,ABCDtαqt=G2t,p,q,r,s,ABCDtαrt=G3t,p,q,r,s,ABCDtαst=G4t,p,q,r,s,where(19)G1t,p,q,r,s=artst+γqt+δ,G2t,p,q,r,s=−δ−βrtst−γqt−stqt,G3t,p,q,r,s=−βrtst−qt−stpt,G4t,p,q,r,s=stpt+stqt+qt.Let us consider system (14) as(20)DABCtαptξt=Gt,ξt,with an initial condition ξ0=ξ0≥0 where(21)ξt=p,q,r,sT,ξ0=p0,q0,r0,s0T,Gt,ξt=Gnt,p,q,r,sT,n=1,2,3,4.In (21), the superscript T represents the transpose. By using Lemma 1 and AB fractional integral, the (20) becomes the fractional integral equation as(22)ξt=ξ0+1−αNαGt,ξt+αNαΓα∫0tt−yα−1Gy,ξydy.Now, to prove the existence uniqueness, we consider two hypotheses based on Lipschitzian and some growth condition assumptions.Hypothesis 1.
For two constantsϕE,θE, the inequality exists; that is,(23)Gt,ξt≤ϕEξ+θE,t∈0,T.Hypothesis 2.
For a constantME>0 such that(24)Gt,ξ1−Gt,ξ2≤MEξ1−ξ2,for each ξ∈B and T∈0,t.
Let us define two operatorsψ1 and ψ2 as(25)ψ1ξt=ξ0+1−αNαGT,ξT,ψ2ξt=αNαΓα∫0TT−yα−1Gy,ξydy,where B=ψ1+ψ2.Theorem 3.
Consider a closed convex setFϵ=ξ∈B:ξ≤ϵ where ϵ=β2/1−β1 such that β1=1−α/Nα+tα/NαΓαϕE<1,β2=ξ0+1−α/Nα+tα/NαΓαθE and prove that(26)ψ1ξ1+ψ2ξ2∈Fϵ,for ξ1,ξ2∈Fϵ(27)ψ1ξ1+ψ2ξ2≤maxT∈0,tξ0+1−αNαGT,ξT+αNαΓα∫0TT−yα−1Gy,ξydy≤ξ0+1−αNαϕEξ+θE+αNαΓα∫0TT−yα−1ϕEξ+θEdy=ξ0+1−αNα+tαNαΓαθE+1−αNα+tαNαΓαϕEϵ=β1+β2ϵ≤ϵ.
This confirms thatψ1ξ1+ψ2ξ2∈Fϵ.Theorem 4.
Prove thatψ1 is a contraction.
To prove thatψ1 is a contraction suppose ξ,ξ∗∈Fϵ. Then, by using Hypothesis 2, we have(28)ψ1ξ−ψ2ξ∗maxT∈0,t1−αNα|GT,ξT−GT,ξ∗T≤1−αNαMEmaxT∈0,tξT−ξ∗T≤1−αNαMEξ−ξ∗.
As we know that1−α/NαME<1,ψ1 is a contraction mapping.Theorem 5.
Prove thatψ2 is relatively compact.
We can prove thatψ2 is relatively compact by showing that ψ2 is continuous, uniformly bounded, and also equicontinuous.
As we know thatξT is continuous, then ψ2ξT is also continuous.
Let us assume thatξ∈Fϵ, then(29)ψ2ξ≤maxT∈0,tαNαΓα∫0TT−yα−1Gy,ξydy≤αNαΓα∫0TT−yα−1ϕEmaxT∈0,tξ+θEdy≤αNαΓα∫0TT−yα−1ϕEξ+θEdy≤tαNαΓαϕEϵ+θEdy.
Hence, proved thatψ2 is uniformly bounded on Fϵ. Now, we have to show that ψ2 is equicontinuous. Assume ξ∈Fϵ and T1,T2∈0,t where T1<T2. Then, we have(30)ψ2ξT2−ψ2ξT1≤αNαΓα∫T1T2T2−yα−1Gy,ξydy+αNαΓα∫0T1T1−yα−1−T2−yα−1Gy,ξydy≤ϕEϵ+θENαΓαT2−T1α+T1α−T2α+T2−T1α=2ϕEϵ+θENαΓαT2−T1α=limT1⟶T22ϕEϵ+θENαΓαT2−T1α⟶0.
Now, the Arzelá-Ascoli theorem suggests thatψ2 is relatively compact, and hence, it is completely continuous.Theorem 6.
If Hypothesis1 and Hypothesis 2 hold, then the fractional integral equation that is equation (20) which is the solution of equation (12) has at least one solution only if 1−α/NαME where β1 is(31)β1=1−αNα+tαNαΓαϕE<1.
By using Theorems2–5, it is proved that the integral equation given in (22) has at least one solution, and consequently, the DSEK model (14) under consideration also has at least one solution.Theorem 7.
Prove that integral equation (20) has a unique solution if β3=1−α/Nα+tα/NαΓαME<1 under Hypothesis 2.
As we haveψ:B⟶B defined as(32)ψξt=ξ0+1−αNαGT,ξT+αNαΓα∫0TT−yα−1Gy,ξydy.
Letξ,ξ∗∈B and T∈0,t. Then, we have(33)ψξt−ψξ∗≤maxT∈0,T1−αNαGT,ξT−GT,ξ∗T+maxT∈0,TαNαΓα∫0TT−yα−1Gy,ξy−Gy,ξ∗ydy≤1−αNα+tαNαΓαMEξ−ξ∗.
Hence,β3 suggests ψ is a contraction. Hence, (22) has a unique solution which suggests that (14) also has a unique solution.


## 5. Ulam-Hyers Stability for DSEK Model

Stability analysis of nonlinear dynamical models is a must. So, in this work, we use Ulam-Hyers stability for DSEK model12 with some nonlinear functional analysis concepts. Ulam-Hyers stability was introduced in 1940 by [22, 23] as a stability study for functional equations. This acted as a motivator for various researchers, and then, this stability was discussed in many forms. Using the fixed-point technique, the authors in [24] investigated the Hyers-Ulam-Rassias and Hyers-Ulam stability of the fractional Volterra integral-differential equation. In a Banach space, some results on generalized Hyers-Ulam stability of the linear differential equation were introduced in [25]. In [26], the authors investigated the Hyers-Ulam stability of first-order linear differential equations and extended previous results using the integral factor approach. In [27], the Hyers-Ulam-Rassias stability of a certain fractional differential equation was discussed, as well as the Hyers-Ulam stability of a certain fractional differential equation. For a particular family of fractional integrodifferential equations, the stability of Ulam-Hyers, Ulam-Hyers-Rassias, and semi-Ulam-Hyers-Rassias on some intervals was studied in [28]. The Ulam-Hyers and generalized Ulam-Hyers-Rassias stabilities for the solution of a fractional-order pseudoparabolic partial differential equation were investigated using the Gronwall inequality [29]. The existence and uniqueness of solutions, as well as Ulam-Hyers-Rassias stability, of an impulsive certain fractional differential equation were investigated in [30]. Sometimes, it is the stability analysis of differential equation ordinary or partial, integral equations, functional equations, etc. Various types have been formed of Ulam-Hyers stability theory, namely, Ulam-Hyers-Rassias, semi-Ulam-Hyers-Rassias [28, 31], and Ulam stability [32].Definition 6.
For someλ>0,ξ˜∈B if(34)DABCtαξ˜T−GT,ξ˜T≤λ,there must exist ξ∈B that satisfies the DSEK model (14) having an initial condition(35)ξ0=ξ˜0,where ξ˜−ξ≤ϵλ such that(36)ξt=p˜,q˜,r˜,s˜T,ξ0=p0˜,q0˜,r0˜,s0˜T,Gt,ξ˜t=Gnt,p˜,q˜,r˜,s˜T,n=1,2,3,4,λ=mλ1,λ2,λ3,λ4T,ϵ=mϵ1,ϵ2,ϵ3,ϵ4T,with this property, if there exists an ϵ>0, then it is said that the DSEK model (12) is UlamHyers stable.Remark 1.
Letf be a small perturbation such that f∈C0,t where f0=0 has the properties given as follows:(i)
fT≤λ, where T∈0,t and λ1>0(ii)
ForT∈0,t the model becomes(37)DABCtαξ˜T=GT,ξ˜T+fT,where fT=f1T,f2T,f3T,f4TT the superscript T represents the transpose.Lemma 2.
Perturbed system (35)(38)DABCtαξ˜T=GT,ξ˜T+fT,ξ˜0=ξ˜0,has a solution that satisfies the inequality(39)ξ˜fT−ξ˜T≤lλ,ξ˜fT represents the solution of the system (37),(40)l=Γα1−α+tαNαΓα.
By using Remark1 and Lemma 2, the solution of system (37) is given as(41)ξ˜fT=ξ˜0+1−αNαGT,ξ˜T+αNαΓα∫0TT−yα−1Gy,ξ˜ydy+1−αNαfT+αNαΓα∫0TT−yα−1Gydy.
Also, we know that(42)ξ˜T=ξ˜0+1−αNαGT,ξ˜T+αNαΓα∫0TT−yα−1Gy,ξ˜ydy.
Now, Remark1 suggests that(43)ξ˜fT−ξ˜T≤1−αNαfT+αNαΓα∫0TT−yα−1fyd≤Γα1−α+tαNαΓαλ=lλ.Theorem 8.
By using Theorem4, it is proven that the DSEK system (35) is Ulam-Hyers stable in B. Let the DSEK system (12) with initial conditions(44)ξ0=ξ˜0,has a unique solution as ξ∈B and ξ˜∈B is the solution of inequality (34), and then,(45)ξT=ξ0+1−αNαGT,ξT+αNαΓα∫0TT−yα−1Gy,ξydy.
Sinceξ0=ξ˜0 as suggested by an initial condition, hence, (45) becomes(46)ξ˜T=ξ˜0+1−αNαGT,ξ˜T+αNαΓα∫0TT−yα−1Gy,ξ˜ydy.
Then, by Lemma2 and the hypothesis above, we have(47)ξ˜T−ξT≤ξ˜T−ξ˜fT+ξ˜fT−ξT≤lλ+1−αNαGT,ξ˜T−GT,ξT+αNαΓα∫0Tt−yα−1GT,ξ˜T−GT,ξTdy+lλ≤2lλ+1−αNα+tαNαΓαMEξ˜−ξ,which implies that(48)ξ˜−ξ≤2lλ1−β3.
Sinceβ3<1 hence from ϵ=2l/1−β3, we obtain ξ˜−ξ≤ϵλ. Hence, proved that the DSEK system (35) is Ulam-Hyers stable.


## 6. Numerical Approximation of Fractional DSEK Model

In this section, the fractional DSEK model will be solved numerically by using fractional Euler’s method. There are several numerical techniques to compute the numerical results of a fractional system of differential equations, but in this case, even Euler’s method can analyze its solution. In order to do that, we first derive the fractional Euler’s method for Atangana-Baleanu Caputo fractional derivative.


### 6.1. Fractional Euler’s Method for Atangana-Baleanu Caputo Fractional Derivative

The authors of [33] proved that the generalized Taylor’s formula of Atangana-Baleanu Caputo fractional derivative is given as(49)ft=∑m=0nDABCtαn+1∑k=0n+1xαkn+1!αk1−α−k+n+1k!−k+n+1!Bαn+1Γkα+1+DABCtαm∑k=0mxαkm!αk1−αm−kk!m−k!BαmΓkα+1t−αfa.Suppose that we have an initial value problem(50)ABCDtαft=Gt,ft,f0=f00<α≤1,t>0.Let0,a be the interval on which we need to obtain the solution of our problem. For generalization instead of 0,a, we consider ti,fti and use this for our approximation. Let the k subintervals of equal width be h=a/k by using nodes ti=ih for i=0,1,…,k. Consider that ft,ABCDtαft, etc. are continuous on ti,fti, then, by using (50), we expand ft about t=t0 as(51)ft1=ft0+ABCDtαft0αt1−aαΓα+1Bα+⋯.Upon neglecting higher-order terms because step sizeh is considered as a smallest positive number and taking h=t1,ABCDtαft0=Gt0,ft0, (51) becomes(52)ft1=ft0+Gt0,ft0αh−aαΓα+1Bα.Equation (52) becomes the iterative equation for repeatedly calculating the points of t that approximates the solution of ft. Hence, the general form of fractional Euler’s method for solving initial value problems with Atangana-Baleanu Caputo fractional derivative is(53)ti+1=ti+h.(54)fti+1=fti+Gti,ftiαh−aαΓα+1Bα.It can be observed easily that forα=1, this becomes the classical Euler method.Now, to solve the fractional DSEK model numerically, we use the parameter values and initial conditions given in [18]. According to that table, a=100,β=50,γ=3.32mm, δ=0.015mm,p0=7.50mm,q0=24.39mm,r0=5.63mm, and s0=0.52mm. Now, the fractional DSEK system in (14) with the iterative formula (54) becomes(55)pti+1=p0+αh−aαΓα+1Bα3.32qti+100rtisti+0.015,qti+1=q0+αh−aαΓα+1Bαqti−sti−3.32qti−50rtisti−0.015,rti+1=r0+αh−aαΓα+1Bα−ptisti−qti−50rtisti,sti+1=s0+αh−aαΓα+1Bαptisti+qtisti+qti.By solving (55) with the help of software, we obtain the numerical solution for different values of α in the form of Figures 2–5 and 6. Figure 2(a) represents the fractional solution of refractive index pt for different values of α between 0<α≤1 and h=0.01, whereas Figure 2(b) is obtained for the step size h=0.001.Figure 2
Numerical simulation of refractive index for different values ofα and h implementing the fractional behavior. (a)h=0.01,0<α≤1 and t∈1,10. (b)h=0.0001,0<α≤1 and t∈1,10.
(a)(b)Figure 3
Numerical solution by fractional Euler’s method of axial length at different values ofα and h implementing the fractional behavior. (a)h=0.01,0<α≤1 and t∈1,10. (b)h=0.0001,0<α≤1 and t∈1,10.
(a)(b)Figure 4
Numerical solution by fractional Euler’s method of corneal curvature at different values ofα and h implementing the fractional behavior. (a)h=0.01,0<α≤1 and t∈1,10. (b)h=0.0001,0<α≤1 and t∈1,10.
(a)(b)Figure 5
Numerical solution by fractional Euler’s method of corneal thickness at different values ofα and h implementing the fractional behavior. (a)h=0.01,0<α≤1 and t∈1,10. (b)h=0.0001,0<α≤1 and t∈1,10.
(a)(b)Figure 6
Numerical solution by fractional Euler’s method of fractional DSEK system at different values ofα and h implementing the fractional behavior. (a)h=0.0001,0<α≤1 and t∈1,10. (b)h=0.01,0<α≤1 and t∈1,10.
(a)(b)The only difference among the solutions presented in Figures2(a) and 2(b) is the different values of h. If we observe as h⟶0, the fractional behavior is clearer to understand and gives us the accurate approximation for α=1, then h⟶1.Figures3(a) and3(b) represent the numerical solution of axial length for different values of α where 0<α≤1 and h=0.01,h=0.0001, respectively. By observing closely, the solutions depicted in Figure 3(a) show that qt⟶0 as t⟶∞. Since this model represents a real-life case of eye surgery, hence, this result is unacceptable. As for h⟶0, the graphical results in Figure 3(b) are more accurate because it suggests the qt≈24.0 or lies closer to 24.0 as t⟶∞.Figures4(a) and4(b) are the graphical illustration of corneal curvature for different values of 0<α≤1 and step sizes as h=0.01 and h=0.0001, respectively. Similar to the refractive index and axial length, the corneal curvature also depicts more realistic behavior when h⟶0.Figures5(a) and 5(b) represent the numerical solution of corneal thickness for different values of α where 0<α≤1 and h=0.01,h=0.0001, respectively. By observing closely, the solutions depicted in Figure 5(a) show that qt⟶0 as t⟶∞. Since this model represents the real-life case of eye surgery, hence, this result is unacceptable.As forh⟶0, the graphical results in Figure 5(b) are more accurate because it suggests the st≈0.52mm or lies closer to 0.52mm as. t⟶∞Figures6(a) and 6(b) show the numerical solution of fractional DSEK by fractional Euler’s method. In Figure 6(a), results have been presented for h=0.0001 and different values of α between 0,1.Graphical results are shown in Figures2–5 and 6 described that as h⟶0, the more accurate results we obtain. This is why the variables remain in ℜ+4 for h=0.0001 instead of h=0.001.Also, the fractional DSEK results have the hysteresis phenomenon, which means this system is influenced by the previous derivatives and values as well as the current conditions. The noninteger derivative given by different values ofα introduces the memory effect in the fractional DSEK model. As we explained in definitions (4) and (5), the exponential kernel when applied to the fractional DSEK model calculates the memory effect. This is why we can see the smoothness in Figures 2–5 and 6 as compared to graphical results in [18]. Results of refractive index, axial length, corneal curvature, and central corneal thickness are shown graphically of the ordinary system of differential equation in [18] showed huge oscillation whereas, in real life after surgery, the effect on vision is not that blurry or oscillated. The fractional DSEK model shows more realistic results of ocular parameters after Descemet’s stripping endothelial keratoplasty. It gives the same normal values but due to its fractal phenomenon, the oscillation among results is removed, and graphs are smoother giving the same normal values as the DSEK model in [18]. For more background about the numerical solutions of fractional-order differential equations, see [34–37].


## 6.1. Fractional Euler’s Method for Atangana-Baleanu Caputo Fractional Derivative

The authors of [33] proved that the generalized Taylor’s formula of Atangana-Baleanu Caputo fractional derivative is given as(49)ft=∑m=0nDABCtαn+1∑k=0n+1xαkn+1!αk1−α−k+n+1k!−k+n+1!Bαn+1Γkα+1+DABCtαm∑k=0mxαkm!αk1−αm−kk!m−k!BαmΓkα+1t−αfa.Suppose that we have an initial value problem(50)ABCDtαft=Gt,ft,f0=f00<α≤1,t>0.Let0,a be the interval on which we need to obtain the solution of our problem. For generalization instead of 0,a, we consider ti,fti and use this for our approximation. Let the k subintervals of equal width be h=a/k by using nodes ti=ih for i=0,1,…,k. Consider that ft,ABCDtαft, etc. are continuous on ti,fti, then, by using (50), we expand ft about t=t0 as(51)ft1=ft0+ABCDtαft0αt1−aαΓα+1Bα+⋯.Upon neglecting higher-order terms because step sizeh is considered as a smallest positive number and taking h=t1,ABCDtαft0=Gt0,ft0, (51) becomes(52)ft1=ft0+Gt0,ft0αh−aαΓα+1Bα.Equation (52) becomes the iterative equation for repeatedly calculating the points of t that approximates the solution of ft. Hence, the general form of fractional Euler’s method for solving initial value problems with Atangana-Baleanu Caputo fractional derivative is(53)ti+1=ti+h.(54)fti+1=fti+Gti,ftiαh−aαΓα+1Bα.It can be observed easily that forα=1, this becomes the classical Euler method.Now, to solve the fractional DSEK model numerically, we use the parameter values and initial conditions given in [18]. According to that table, a=100,β=50,γ=3.32mm, δ=0.015mm,p0=7.50mm,q0=24.39mm,r0=5.63mm, and s0=0.52mm. Now, the fractional DSEK system in (14) with the iterative formula (54) becomes(55)pti+1=p0+αh−aαΓα+1Bα3.32qti+100rtisti+0.015,qti+1=q0+αh−aαΓα+1Bαqti−sti−3.32qti−50rtisti−0.015,rti+1=r0+αh−aαΓα+1Bα−ptisti−qti−50rtisti,sti+1=s0+αh−aαΓα+1Bαptisti+qtisti+qti.By solving (55) with the help of software, we obtain the numerical solution for different values of α in the form of Figures 2–5 and 6. Figure 2(a) represents the fractional solution of refractive index pt for different values of α between 0<α≤1 and h=0.01, whereas Figure 2(b) is obtained for the step size h=0.001.Figure 2
Numerical simulation of refractive index for different values ofα and h implementing the fractional behavior. (a)h=0.01,0<α≤1 and t∈1,10. (b)h=0.0001,0<α≤1 and t∈1,10.
(a)(b)Figure 3
Numerical solution by fractional Euler’s method of axial length at different values ofα and h implementing the fractional behavior. (a)h=0.01,0<α≤1 and t∈1,10. (b)h=0.0001,0<α≤1 and t∈1,10.
(a)(b)Figure 4
Numerical solution by fractional Euler’s method of corneal curvature at different values ofα and h implementing the fractional behavior. (a)h=0.01,0<α≤1 and t∈1,10. (b)h=0.0001,0<α≤1 and t∈1,10.
(a)(b)Figure 5
Numerical solution by fractional Euler’s method of corneal thickness at different values ofα and h implementing the fractional behavior. (a)h=0.01,0<α≤1 and t∈1,10. (b)h=0.0001,0<α≤1 and t∈1,10.
(a)(b)Figure 6
Numerical solution by fractional Euler’s method of fractional DSEK system at different values ofα and h implementing the fractional behavior. (a)h=0.0001,0<α≤1 and t∈1,10. (b)h=0.01,0<α≤1 and t∈1,10.
(a)(b)The only difference among the solutions presented in Figures2(a) and 2(b) is the different values of h. If we observe as h⟶0, the fractional behavior is clearer to understand and gives us the accurate approximation for α=1, then h⟶1.Figures3(a) and3(b) represent the numerical solution of axial length for different values of α where 0<α≤1 and h=0.01,h=0.0001, respectively. By observing closely, the solutions depicted in Figure 3(a) show that qt⟶0 as t⟶∞. Since this model represents a real-life case of eye surgery, hence, this result is unacceptable. As for h⟶0, the graphical results in Figure 3(b) are more accurate because it suggests the qt≈24.0 or lies closer to 24.0 as t⟶∞.Figures4(a) and4(b) are the graphical illustration of corneal curvature for different values of 0<α≤1 and step sizes as h=0.01 and h=0.0001, respectively. Similar to the refractive index and axial length, the corneal curvature also depicts more realistic behavior when h⟶0.Figures5(a) and 5(b) represent the numerical solution of corneal thickness for different values of α where 0<α≤1 and h=0.01,h=0.0001, respectively. By observing closely, the solutions depicted in Figure 5(a) show that qt⟶0 as t⟶∞. Since this model represents the real-life case of eye surgery, hence, this result is unacceptable.As forh⟶0, the graphical results in Figure 5(b) are more accurate because it suggests the st≈0.52mm or lies closer to 0.52mm as. t⟶∞Figures6(a) and 6(b) show the numerical solution of fractional DSEK by fractional Euler’s method. In Figure 6(a), results have been presented for h=0.0001 and different values of α between 0,1.Graphical results are shown in Figures2–5 and 6 described that as h⟶0, the more accurate results we obtain. This is why the variables remain in ℜ+4 for h=0.0001 instead of h=0.001.Also, the fractional DSEK results have the hysteresis phenomenon, which means this system is influenced by the previous derivatives and values as well as the current conditions. The noninteger derivative given by different values ofα introduces the memory effect in the fractional DSEK model. As we explained in definitions (4) and (5), the exponential kernel when applied to the fractional DSEK model calculates the memory effect. This is why we can see the smoothness in Figures 2–5 and 6 as compared to graphical results in [18]. Results of refractive index, axial length, corneal curvature, and central corneal thickness are shown graphically of the ordinary system of differential equation in [18] showed huge oscillation whereas, in real life after surgery, the effect on vision is not that blurry or oscillated. The fractional DSEK model shows more realistic results of ocular parameters after Descemet’s stripping endothelial keratoplasty. It gives the same normal values but due to its fractal phenomenon, the oscillation among results is removed, and graphs are smoother giving the same normal values as the DSEK model in [18]. For more background about the numerical solutions of fractional-order differential equations, see [34–37].


## 7. Conclusion

In this paper, we investigated the fractional DSEK model presented by fractional derivatives of the Atangana-Baleanu Caputo type. We proved the uniqueness and existence of its solutions by using fixed point theory. For this, we defined a hypothesis based on Lipschitzian and two operatorsψi,i=1,2. Then, we proved that ψ1 and ψ2 are contraction and relatively compact and hence proven the uniqueness and existence of those defined hypotheses. Furthermore, for the fractional DSEK model, proving its stability was a must so by UlamHyers stability in Banach space, we proved that fractional DSEK is Ulam-Hyers stable. Moreover, we have discussed the advantages of using the ABC fractional derivative instead of any other. In this paper, we presented and investigated the fractional behavior of the DSEK model and performed the numerical investigation using mathematical software. The numerical method “Euler” which is used to solve fractional DSEK is derived for initial value problems with ABC fractional derivatives, in this paper. Since eye surgery is a crucial process and with the passage of time, the results of surgery can be observed but with the help of the fractional DSEK model, a clearer picture of this surgery will be given.



---

*Source: 4475491-2022-05-26.xml*