New Series Solution of the Caputo Fractional Ambartsumian Delay Differential Equationation by Mittag-Lefﬂer Functions

: The fractional generalization of the Ambartsumian delay equation with Caputo’s fractional derivative is considered. The Ambartsumian delay equation is very difﬁcult to be solved neither in the case of ordinary derivatives nor in the case of fractional derivatives. In this paper we combine the Laplace transform with the Adomian decomposition method to solve the studied equation. The exact solution is obtained as a series which terms are expressed by the Mittag-Lefﬂer functions. The advantage of the present approach over the known in the literature ones is discussed.


Introduction
One of the successfully applied differential equations as a model is the Ambartsumian equation. This equation was derived in its original version by Ambartsumian [1] to describe the absorption of light by the interstellar matter. Later, many authors studied various properties of the solution of this equation. It is worth to mention the works about the existence and uniqueness [2] and about the introduced power series solution [3]. Note that several methods are applied to solve the classical Ambartsumian equation, such as, the Adomian decomposition method (ADM) combined with the Laplace Transform (LT) [4], the Homotopy perturbation method (HPM) [5]. Recently, based on the main property of fractional derivatives, their memory property, the ordinary derivative is replaced by a fractional one to make the model of the absorption of light by the interstellar matter more adequate. Besides, the standard Ambartsumian equation has been deducted by Ambartsumian [1] based on basic concepts in matter physics. However, the formulation of the current fractional model may requires additional information about the nature of interstellar constellation. Moreover, the present paper is devoted to obtaining analytical approximations in closed form which posses many advantages over the existing approaches in the relevant literature.
Several fractional generalizations of Ambartsumian equation are defined and studied (for example, the conformable derivative has been applied in [6,7], the Caputo fractional derivative has been used in [8]). In paper [8], the authors applied the homotopy transform analysis method (HTAM) to Caputo fractional generalization of Ambartsumian equation to obtain the solution as a power series in terms of powers of the fractional order.
Mathematics 2021, 9, 157 2 of 18 In this paper we study the fractional generalization of the Ambartsumian delay equation with Caputo fractional derivative. We combine both well known the ADM and LT to solve the studied equation. We obtain a new type of the exact solution for the given fractional model. This solution is in the form of a series in terms of the Mittag-Leffler functions. The explicit formulas for all terms of the series solution are provided. It is one of the advantages of the obtained solution comparatively with the known in the literature (for example, compare with the recently published paper [23] in which the explicit formulas only of the first several terms are obtained). The other advantage of the obtained solution is the proved convergence in the whole domain of the current model. Some numerical simulations are provided to illustrate the convergence, the influence of both the delay parameter and the initial value and to compare the obtained series solution with some known in the literature.

Some Preliminary Results from Fractional Calculus
Based on engineering reasoning, we will consider the case when the fractional order α ∈ (0, 1).
In this section we will provide some well known results from fractional calculus which will be used in our further study.
In this paper we will use the following definitions for fractional derivatives and integrals for scalar functions y : [0, T] → R with T ≤ ∞ (note in the case T = ∞ the interval will be half open): - where Γ is the Gamma function. -

Remark 1.
In the case α → 1 the Caputo fractional derivative is reduced to ordinary derivative (for more details, see, for example [24]).
The Laplace transform (LT) of the Caputo fractional derivative is where Y(s) is the LT of y(t). In our study we will use Mittag-Leffler functions of two parameter defined by in particular E α,1 (z) = E α (z). It follows from the definition that The LT of the expression of Mittag-Leffler function is ( [24]): Mathematics 2021, 9, 157 3 of 18 As partial cases of (3) we obtain the following equalities which will be used later (for more details see [24]): Some useful properties of Mittag-Leffler function are given by ( [24])

Statement of the Problem
In this paper we will study the following Caputo fractional Ambartsumian equation where y ∈ R, C 0 D α t y(t) is the Caputo fractional derivative of y(t), and q > 1 is a constant. The model is subjected to the initial condition: where λ is a constant. We will provide a procedure for obtaining a series solution of the Caputo fractional Ambartsumian Equationations (11) and (12). We will combine both methods LT and ADM (we will call this procedure combined LT-ADM).

General Slgorithm of Combined LT-ADM
Apply the LT to Equationations (11) and (12), apply equality (1) and obtain where Y(s) is the LT of y(t) and Y(qs) is the LT of 1 q y t q (see [25] for details of the LT). In order to apply the ADM, we rewrite Equation (13) in the canonical form: It is well known that the ADM assumes the solution Y(s) as a series in the form where Y i (s), i = 1, 2, . . . are unknown functions. Now, we will use LT to obtain these functions. Substitute the series (15) in (13) and obtain Following the ideas of [1] we get From the recurrence scheme (17) step by step we could obtain all terms Y i (s), i = 1, 2, . . . in the series (15). or Using partial fractions, we have where A 1 and A 2 are given by Substitute (20) in (19), we obtain the first term in the series (15) Let i = 2. From recurrence equality (17) and Equation (18) we have Also, Similarly to the case of i = 1, by partial fractions, we can write where the unknown coefficient are given by Substitute the coefficients B i , i = 1, 2, 3, defined by (25), in equalities (24) and (23) and obtain the second term in the series (15) Let i = 3. From equalities (17), and (22) we get and applying (26) and some simplifications we obtain where the coefficients C 1 , C 2 , C 3 , and C 4 are given by Hence, from (29) and (28) we obtain the third term in the series (16) where coefficients C i , i = 1, 2, 3, 4 are given by (29). By induction from recurrence equalities (17) we get Proceeding as above, one can obtain higher-order components of Y i (s) in the form where coefficients K (i) j , j = 1, 2, . . . , i + 1 depend on q and α. Note the explicit expressions for the coefficients K (i) j are very complicated and we will ignore them at this stage. Later we will use only the first three components which obtained above. Now, using the terms Y i (s), i = 1, 2, . . . , of the series (15) we will obtain the series solution of the initial value problem (11), (12). Apply the inverse LT to the series (15) to obtain the solution y(t) of (11) and apply (6) on Y(s). We obtain where from (17), (21), (4) and (6) we get Similarly, y 2 (t) and y 3 (t) can be obtained as and Similarly, from (32), the equalities (4)-(6) for LT, we can obtain y i , i ≥ 4, in the form: where coefficients K (i) j , j = 1, 2, 3, . . . , i − 1, depend on q and α. Unfortunately, we could obtain their explicit formulas only for i = 1, 2, 3. Therefore, the series solution given by (33) with terms y i (t) defined by (37) practically is not applicable. We could use it only for the third partial sum since only these terms are given explicitly by (34)-(36).

Series Solution with One-Parameter Mittag-Leffler Functions
This section applying the ideas for obtaining the series solution (37) with coefficients (33) and some obtained formulas, we will get general explicit formulas for the terms y i , i = 1, 2, 3, . . . .
We will use the presentation (31). Define Then from (31) for any i ≥ 1 we get Now, apply the inverse LT to (39) and obtain the general term y i (t) of the series solution in the form where ( * ) refers to the convolution of f (t) and g(t) with where W i (s) = ∏ i k=1 s + q −kα and W i (s) is its derivative. Substitute (41) in (40) and obtain From integral formula (10) with β = 1, γ = α, a = −q −kα , and b = −1 we obtain From (42) and (43) we obtain From (33) and (44) we obtain the general form of the series solution of (11) in the form It is easy to verify that in formula (45) the components y 1 (t), y 2 (t), and y 3 (t) are the same as the ones obtained in Section 3.2 and explicitly given by (34)-(36).
Although, all terms of the series solution (45) are given explicitly, we will continue to simplify further.
From formula (8) with z = −t α and z = −q −kα t α , respectively, we have Substitute (46) and (47) in (45) and get According to the definition of function W i (s) we have W i (s) = ∑ i m=1 ∏ i j=1, j =m (s + q −jα ) and W i (−q −kα ) = ∏ i j=1, j =k (q −jα − q −kα ). Finally, we obtain the formula for the series solution of Equation (11) given by It is noticed that in formula (49) only one-parameter Mittag-Leffler function is used. Accordingly, the formula (49) has an advantage comparatively with (45) in which twoparameter Mittag-Leffler functions are used. In addition, it can be easily verified that the expression (49) satisfies the initial condition (12), i.e., the series (49) is a solution of the initial value problem (11) and (12).
In connection with our further considerations we will define the n-partial sum ρ n of the series (49) and we will call it n-th approximate solution of (11) and (12): Also, define the absolute residual error RE n (t) by Applying the equality C 0 D α t E α [at α ] = aE α [at α ], t ≥ 0, a ∈ R we get

Limit Case as α → 1 of the Series Solution
According to Remark 1 the Caputo fractional derivative is reduced to an ordinary one as α → 1. Thus, the studied Caputo fractional generalzation (11) of Ambartsumian differential equation will be reduced to the well known and studied in the literature model Now, take limit α → 1 in the obtained series solution (49), use the equality E 1 [u] = e u , u ∈ R and obtaiñ As partial cases of (56) we obtain y 0 (t) = λe −t , , The Equation (55) with initial condition (12) is studied by many authors and different types of series solutions are obtained. In [15] a series solutions with exponential functions is obtained. If we compare the formulas (10), (11) (with corrected typos) for the first three in terms, obtained in [15], one can see, that they coincide the formulas (57), obtained as a limit case of the new formula (49). Therefore, the new series solution (49) in fractional case α ∈ (0, 1) is a generalization for the ordinary case known in the literature.

Numerical Simulations and Discussions
We will discuss the series solution (49) in different point of view by computer simulations. We will apply CAS Wolfram Mathematica to simulate and graph the results. Case 1. Convergence. Consider the n-th partial sums ρ n (t) of the series solution (49), called n-th approximate solution of (11), (12) and the error function RE n (t) defined by (50) and (51), respectively. Fix λ = 1, q = 1.5, α = 0.5 and graph the partial sums (n-th approximate solution of (11), (12) ) for various values of n = 2, 4, 7, 10, 15 (see Figure 1) and the corresponding errors (see Figure 2). In addition, the numerical results of the residual error RE 10 (t) at various values of α are tabulated in Table 1 when q = 1.5 and Table 2 when q = 2. From both Figures and Tables it could be seen the convergence of the n-th approximate solution of (11), (12) with increasing of n. Moreover, the obtained numerical values for the residual error RE 10 (t) reflects the accuracy of the present analysis.
RE 15 (t) Figure 2. Graphs of the n-th errors RE n (t) with n = 2, 4, 7, 10, 15 and λ = 1, q = 1.5, and α = 0.5. Case 2. Impact of the parameter q. We fix λ = 1, α = 0.5 and graph the n-th approximate solution of (11), (12) at n = 2, 4, 7, 10, 15 for various values of q = 1.5 (see Figure 1), q = 1.8 (see Figure 3), q = 2 (see Figure 4) and q = 2.5 (see Figure 5). A rapid convergence is observed from these figures using only a few terms of the series solutions (49). We could see, that the rate of convergence is increased for higher values of q, for example q > 2, where at q = 2.5 (see Figure 5) the 4-term, 7-term, 10-term, and 15-term are nearly identical. However, a rapid decrease in the surface brightness has been remarked by increasing the delay parameter q. This latest notice reveals that the curves of the approximate solutions tend faster to zero at higher values of q.
Case 3. Impact of the initial value λ. To study the impacts of the initial value λ on the approximation, we fix q = 1.5, α = 0.5, n = 10, and graph the 10-th approximate solution ρ 10 (t) of (11), (12) for various values of λ = 1, 2, 3 (see Figure 6). It can be seen from Figure 6 that the surface brightness is increased by increasing the given initial value λ.     6. Impact of the initial value λ on the 10-term approximate solution of (11), (12) at q = 1.5, and α = 0.5.

Case 4.
Comparing the series solution (49) with known ones in the literature. The initial value problem (11), (12) is studied by some other authors and they obtained some different types of series solution. For example, in [19], the following series solution is obtained: with n-th partial sum (n-th approximate solution of (11), (12)) defined by We fix λ = 1, α = 0.5, q = 1.4 and graph both n-the approximate solution ρ n (t) and Ψ n (t) of (11), (12) for various values of n = 10 (see Figure 7), n = 15 (see Figure 8). The numerical values of ρ 15 (t) (present), Ψ 15 (t) (Ref. [19]), and the difference d n (t) = ρ n (t) − Ψ n (t) are tabulated in Table 3. It can be seen from this table that the difference d n (t) increases as t increases and hence d n (t) is significant. The main observation here is that the values of Ψ 15 (t) are negative at some points such as at t = 8 and t = 10 (see Table 3) which is not physically acceptable. By this, the published solution in the literature [19] is not physical in the whole domain. Therefore, our approach is applicable and convergent in the whole domain which are the main advantages of the present work over the corresponding one in the literature [19]. Table 3. Comparison between ρ 15 (t) (present), Ψ 15 (t) (Ref. [19]) when λ = 1, q = 2 and α = 1. d n (t) = ρ n (t) − Ψ n (t) t ρ 15 (t) (Present) Ψ 15 (t) [

Conclusions
A combined approach based on the LT and the ADM was developed in this paper to solve the fractional model of the Ambartsumian delay equation with Caputo fractional derivative. The first step of such approach was applying the LT on current fractional equation and then solving the transformed equation by the ADM. It was also shown that a general formula for the y i (t)-components (Adomian's components) was successfully obtained. Hence, the exact solution was provided to the first time for fractional model of the Ambartsumian equation. Both of the approximate and the exact solutions were obtained and expressed in terms of the Mittag-Leffler functions. Moreover, the advantage of the present approach over the previous one in the relevant literature was that our solution converges in the whole domain as properties of the Mittag-Leffler functions. While the solution in the relevant literature [7] was given as a power series in terms of t α which converges in certain domains, i.e., not valid in the whole domain of the present model. The results obtained by simulations of the studied equation reveal that the approximate solution is highly accurate. Moreover, the absolute residual error approaches zero even at higher values of the delay parameter. Finally, the authors believe that the current developed approach deserves further extensions to include and investigate other higherorder fractional delay equations.