Existence in the Large for Caputo Fractional Multi-Order Systems with Initial Conditions

: One of the key applications of the Caputo fractional derivative is that the fractional order of the derivative can be utilized as a parameter to improve the mathematical model by comparing it to real data. To do so, we must ﬁrst establish that the solution to the fractional dynamic equations exists and is unique on its interval of existence. The vast majority of existence and uniqueness results available in the literature, including Picard’s method, for ordinary and/or fractional dynamic equations will result in only local existence results. In this work, we generalize Picard’s method to obtain the existence and uniqueness of the solution of the nonlinear multi-order Caputo derivative system with initial conditions, on the interval where the solution is bounded. The challenge presented to establish our main result is in developing a generalized form of the Mittag–Lefﬂer function that will cooperate with all the different fractional derivative orders involved in the multi-order nonlinear Caputo fractional differential system. In our work, we have developed the generalized Mittag–Lefﬂer function that sufﬁces to establish the generalized Picard’s method for the nonlinear multi-order system. As a result, we have obtained the existence and uniqueness of the nonlinear multi-order Caputo derivative system with initial conditions in the large. In short, the solution exists and is unique on the interval where the norm of the solution is bounded. The generalized Picard’s method we have developed is both a theoretical and a computational method of computing the unique solution on the interval of its existence.


Introduction
Fractional differential equations have various applications in widespread fields of science, such as engineering [1], chemistry [2][3][4], biology [5,6], physics [7][8][9], numerical analysis [10,11], and others [12,13]. In addition, Caputo fractional differential equations have the same initial and boundary conditions as the corresponding integer order dynamic equations. To use the order of the fractional derivative as a parameter to enhance the model, we need to show that the solution to the Caputo fractional differential equation exists and is unique on its interval of existence. However, a majority of the theoretical and computational methods of establishing existence and uniqueness results of dynamic equations are local existence results only. This includes Picard's method as well, which is both a theoretical and a computational method of proving the existence and uniqueness results for dynamic equations with initial conditions. In this work, we generalize Picard's method to provide an existence and uniqueness result in the large for the solution of the nonlinear multi-order Caputo differential system with initial conditions. The incorporation of fractional derivatives of order 0 < q < 1 allows for generalizations of the standard nonlinear system to multi-order systems where each iteration of the system has its own derivative order. A system of the form would be a simple illustration of a multi-order system. However, the computation of the solution on its interval of existence is not guaranteed by Picard's method. In [14], the authors have established that the fractional order of the derivative can be chosen as a parameter to enhance the model. With the modeling point in view, in this paper, we develop an existence result, namely existence in the large, for nonlinear multi-order fractional systems incorporating the Caputo fractional derivative with respect to time. In [15][16][17], the authors developed results for Riemann-Liouville multi-order systems. In their work, they constructed iterative techniques for various cases of multi-order systems, and their existence results were developed via upper and lower solution techniques. In [18], existence in the large was proven for the Caputo fractional reaction diffusion equation by Picard's method. In [19,20], the authors studied two and three-order systems of Caputo fractional differential equations of the same fractional order derivative with initial conditions using the Laplace transform method.
In this paper, we establish that the solution to the nonlinear multi-order Caputo differential system exists and is unique on the entire time interval where the norm of the solution is bounded. This is what is referred to as existence in the large.
A common challenge when working with multi-order systems is that the generalized exponential function, also known as the Mittag-Leffler function, depends on the order of the derivatives in question. That is, D c q E q (λt q ) = λE q (λt q ) only if the order of the derivative and the parameter for the Mittag-Leffler function are equal. Overcoming this challenge requires the construction of a generalized form of the Mittag-Leffler function that will cooperate with all the fractional derivative orders in the system.
In Section 2, we recall known results of fractional differential equations. In Section 3, we establish a nice formula for the generalized Mittag-Leffler function which cooperates with all the different order derivatives of the system and proves its convergence is uniform, allowing term-by-term fractional calculus operations. In Section 3, we also prove the existence of a unique solution of the linear Caputo multi-order system and then develop a comparison theorem. This comparison theorem leads to a Gronwall-type inequality that we incorporate in the proof of our main result. Our main result is developed in Section 4, where we establish a generalized Picard's existence and uniqueness result for the solution to the nonlinear Caputo multi-order system. Finally, we extend this result from a finite domain [0, T] for any T, T < ∞ where the solution of the system is bounded. It should be noted that in biological models, such as population models, the interval of existence is [0, ∞). On the other hand, the solution to u = u 2 , u(0) = 1 is known to exist on the interval [0, 1). However, the unique solution of D c q = u 2 , u(0) = 1, cannot be computed easily. Using a comparison theorem, we can establish that the interval of existence for any q < 1 is at best [0, 1 − ) for small , which depends on the value of q. Here the value of q can be chosen in such a way that our solution matches realistic data. In Section 5, we develop illustrative numerical results for biological models, involving different-order Caputo fractional derivatives.

Preliminaries
In this paper, we examine existence and uniqueness results for multi-order systems of Caputo fractional differential equations of orders 0 < q < 1 over the interval J = [0, T]. Our results employ the Caputo form of the fractional derivative given as for x ∈ C 1 (J). The corresponding fractional integral is defined as In this section, we recall results regarding the scalar nonlinear Caputo fractional initial value problem with t ∈ J and x ∈ C 1 (J). As seen in [12,21,22], if x ∈ C 1 (J), then (1) is equivalent to the fractional integral equation An important function in the discourse of fractional calculus is the Mittag-Leffler function which, for parameters α, β ∈ R, is defined as and is entire for α, β > 0. The one-parameter Mittag-Leffler function defined as is a special case of the Mittag-Leffler function that is integral to the study of Caputo fractional derivatives, since D c q E q (λt q ) = λE q (λt q ). We refer the reader to [23] for more details. Our results focus on the nonlinear multi-order N systems of the form where 0 < q i < 1, x = (x 1 , x 2 , x 3 , . . . , x N ) and η i ∈ R for each i ∈ C = {1, 2, 3, . . . , N}. For simplicity, we will always assume i, or other indices, are in C unless stated otherwise. For convenience, we will be using the following notation D c = ( D c q 1 , D c q 2 , D c q 3 , . . . , D c q N ), and I = (I q 1 , I q 2 , I q 3 , . . . , I q N ).
Thus, system (2) could be written as where η = (η 1 , η 2 , η 3 , . . . , η N ). The notation of multi-order integrals can be cumbersome, so let K i (t, s) = 1 Γ(q i ) (t − s) q i −1 ; additionally, we utilize K i when writing the Riemann-Liouville integral. Therefore, we can write t 0 K i (t, s)φds to mean the q i integral of any function φ. It should be noted that the Riemann-Liouville integral and the Caputo integral of the same order are the same.
In the next section, we develop a generalized Mittag-Leffler function and prove its uniform convergence and integral properties. From there, we establish a uniqueness result for a linear multi-order system and develop a Gronwall-style inequality for use in the proof of uniqueness in our main result for nonlinear systems.

Auxiliary Results
Results for the multi-order system involve a generalized Mittag-Leffler function. This function is an adjustment of the multi-order generalized exponential as discussed in [16]. Specifically, for t ∈ J and a constant λ ∈ R, we consider In [16], the functions constructed were generalizations of the two-parameter Mittag-Leffler function E q i ,q i , and they proved that the series converged uniformly. Thus, with only minor alterations, we can establish that Z converges uniformly on J as well. We show this proof next. Theorem 1. Z(λ, t) converges uniformly on J, and if λ > 0, then for every i, Proof. To begin the proof, we define Z i to be the first i sums of Z(λ, t); for example, . Thus, with this notation, Z N = Z(λ, t) and Z 1 = E q 1 (λt q 1 ), implying that Z 1 converges uniformly.
We will now prove that Z 2 converges uniformly. To do so, we utilize the beta function along with its decreasing nature. B(x, y) decreases in x and y for x, y > 0, since if 0 < x 1 ≤ x 2 , then for s ∈ (0, 1) we have s x 1 −1 ≥ s x 2 −1 , implying B(x 1 , y) ≥ B(x 2 , y). The symmetry of B establishes its monotonicity in both variables. Now, if we split off the k 2 = 0 term from Z 2 , we are left with two series: one only in terms of k 1 and another in terms of both k 1 and k 2 with k 2 ≥ 1. The result is shown below: The first term above is Z 1 , which converges uniformly, so we will show that the double series above converges uniformly.
For t ∈ J and every k 1 ≥ 1, k 2 ≥ 0, we have Now we note that the second fraction in line (4) makes up the terms of E q 2 (λt q 2 ). Next, we consider the first fraction, noting that Applying these results and the Weierstrass M-test, we find that Z 2 converges uniformly on J.
The proof that the remaining Z i series converge uniformly follows by similar arguments. We show the details for Z 3 , exemplifying how the later cases can be developed. We again start by isolating the k 3 = 0 term, yielding As we just proved, Z 2 converges uniformly, so we focus on the triple series above. Thus, for t ∈ J, k 1 , k 2 ≥ 0, and k 3 ≥ 1, we have Similar to previous work, the first fraction makes up the terms of Z 2 , and the second fraction makes up the terms of λT q 3 E q 3 ,q 3 (λT q 3 ). Therefore, by the Weierstrass M-test, Z 3 converges uniformly on J. Continuing this process, we can prove that every Z i converges uniformly, including the final step, which will give us that Z N , that is Z(λ, t), converges uniformly on J. Since N is fixed and finite in this case, we can employ the above steps exhaustively; however, we can also inductively prove that Z N will converge uniformly for all N by utilizing the same steps.
Since convergence is uniform, we can take the fractional integral of the series term by term. To evaluate the fractional integral, we consider i = 1 for simplicity and note that the same result will hold for every i.
Note that the sum in the second line of (5) is missing the k 1 = 0 term and, provided λ > 0, this term will be positive. Further, since the calculation found in (5) will hold for all i, we can say that I q i Z(λ, t) < 1 λ Z(λ, t). Now we look at the homogeneous linear multi-order system with constant coefficients where A is a real valued N × N matrix. Denoting the terms of A = (a i,j ), we can equivalently write (6) as Many of our future results depend on the value of T from the interval J = [0, T]. To conveniently carry out these results, define τ as the constant This constant is used frequently in our upcoming results, since t q i ≤ τ for all i, and t ∈ J for any value of T.
In the following proof that the linear system has a unique solution, we utilize Caccioppoli's extended Banach Contraction Mapping Principle, which we state below. Theorem 2. Let (X, d) be a complete metric space, and let f : X → X be a mapping such that for each n ≥ 1, there exists a constant c n such that d( f n (x), f n (y)) ≤ c n d(x, y), for all x, y ∈ X, where ∑ ∞ n=1 c n < ∞. Then, f has a unique fixed point.
We direct the reader to [24,25] for more details, including the proof. Now we will use this contraction mapping principle to prove that system (6) has a unique solution. Unless stated otherwise, we use the sum norm in our results, that is, Theorem 3. The linear system (6) has a unique solution on J, and the successive approximations converge uniformly to the solution.
Proof. Let J (x) be the mapping For simplicity, let M = max i,j |a i,j |, then Similarly, Continuing the process inductively, we can show that for all n, and since E q i (MNτ) converges uniformly for every i, we obtain Therefore, the mapping J has a unique fixed point by Theorem 2, implying (6) has a unique solution.
The succeeding inequality properties help us prove uniqueness in our main result. Corollary 1. The solution to (6) satisfies the inequality Proof. From Theorem 3, the sequence of successive approximations converges uniformly, and since x n = x 0 + ∑ n k=1 (x k − x k−1 ) for all n, the series ∑ ∞ k=1 (x k − x k−1 ) also converges uniformly. Now, note from the previous proof we have Similarly, continuing this process inductively, we can show that for all n. Now let us consider x 1 − x 0 ; for each i, Γ(q i +1) < 2 for each i. Thus, we have leading to for all n; since convergence is uniform, we have This directly yields which completes the proof.

Remark 1.
We note that Corollary 1 implies that if η = 0, then x = 0 on all of J.
Next, we develop a comparison theorem that provides conditions for when inequalities are preserved from original integral inequalities. This theorem is utilized in the uniqueness proof of our generalized Picard's method.
Proof. To begin, suppose the inequalities in (10) are strict inequalities; we show in this case that z < y. To prove this, suppose to the contrary that the set is non-empty. Let λ = inf(Λ). Since z(0) < y(0) and both are continuous over J, we can conclude that z i (λ) = y i (λ) for some i and that λ > 0. Therefore, since λ is the infimum of Λ and z ≤ y on [0, λ], ∑ N k=1 a ik z k ≤ ∑ N k=1 a ik y k on [0, λ]. Then, which is a contradiction. Therefore, z < y on J. Now we will prove that z ≤ y when the inequalities are not strict. Define the function z ε for each i as z εi = z i − εϕ, where ϕ = Z(2MN, t) as defined in (3). Then note that, since Then for each i, The final inequality comes from Theorem 1, which implies I q i ϕ < 1 2MN ϕ. This then implies that Then by the previous case with strict inequalities, we have that z ε < y on J. Then letting ε → 0, we obtain z ≤ y on J.
In the next section, we employ these results to prove existence in the large for system (2).

Main Result
What follows is a generalized Picard's existence for system (2), which is then extended to existence in the large.
on S. Then (2) has a unique solution on [0, T], and the successive approximations exist on [0, T] and converge uniformly to a solution x of (2).
Proof. First, we note that for all n ≥ 1 we have Therefore we utilize the Weierstrass M-test to show that the series ∑ ∞ k=1 (x k − x k−1 ) converges absolutely and uniformly, which in turn shows that (x n ) converges uniformly. The successive approximations described in (12) for t ∈ [0, T] and for all i. This implies that . Thus, For simplicity, let σ i = q k i ; continuing the previous process inductively, we can show for all n ≥ 2, since t σ i ≤ τ for each i as defined in (7). Now, in the direction of finding a bound on the above sums, let ρ = min i∈C {q i } and note that 1 Γ(z) is decreasing for real numbers z ≥ 2. Further note that there exists a W > 0 such that 1 + ∑ n i=1 σ i ≥ 2 for all n ≥ W. Therefore, we have that , for all n ≥ W.
Applying these results to the sums in (13), for n ≥ W, we get The bounds above describe the terms of the convergent series M L E ρ (LNτ). Therefore, by the Weierstrass M-test, the series ∑ ∞ n=1 (x n − x n−1 ) converges absolutely and uniformly on [0, T] and, in turn, x n converges uniformly on J. Therefore, letting x = lim x n and noting the convergence is uniform we obtain implying this limit is a solution to (2). Now we wish to show this solution is unique, so let y be another solution to (2), and let z be defined as z i = |x i − y i | for each i. Now note that z(0) = 0 and Then by Theorem 4, z ≤x; thus, by Remark 1, and since z ≥ 0, we have that z ≤ x = 0, implying that the solution to (2) is unique.
This result extends naturally to an existence in the large result given below. We note that we only extend this result to the positive t-axis due to the nature of fractional differential equations.
Theorem 6. Assume f is continuous and satisfies the Lipschitz condition (11) on each rectangle Then system (2) has a unique solution existing for all t.
This follows directly from the application of Theorem 5, as its hypotheses will be satisfied on each S λ .

Applications
In this section we construct illustrative numerical examples of the generalized Picard's method developed in Theorem 5. All computations and graph constructions were done with Python 3.9.6 on a 64-bit Windows 10 machine with a 3.4 GHz 6-core processor and 16 GB of RAM. We investigate the epidemic model studied in [20], first showing that the method applies to the single-order system and then expanding to illustrate the multi-order case.
The SIR model considered in this work was D c q S(t) = −βS(t)I(t), where β is the contact rate, and γ is the recovery rate. S(t) is the susceptible population, I(t) is the infected population, and R(t) is the population recovered from the disease. The authors assumed that the birth and death rates were the same for the small period of the epidemic. It is further assumed that there is no immigration and that the recovered population is immune to the disease, yielding a constant population. The authors of [20] utilized β = 0.01, γ = 0.02, L 1 = 20, L 2 = 10, and L 3 = 5 and examined values of q ∈ {0.8, 0.85, 0.9, 0.95, 1}. For the sake of comparison, we utilized the same values in our approximations. The numerical method employed computes the successive approximations from (8) up to x 99 with x 0 = (20, 10,5). The graphs in Figure 1 show the numerical results of the generalized Picard's method, and we find them reasonably close to results found in [20]. The computation speeds for these approximations averaged 7.762 s. For the second example, we extend the above system to the following multi-order system D c q 1 S(t) = −βS(t)I(t), S(0) = L 1 , D c q 2 I(t) = βS(t)I(t) − γI(t), utilizing the same constants as before but now considering the following fractional derivative order groups Q 1 = (0.8, 0.7, 0.8) Q 2 = (0.8, 0.8, 0.7) Q 3 = (0.8, 0.9, 0.9) Q 4 = (0.7, 0.8, 0.9) Q 5 = (0.7, 0.9, 0.8).
These orders are chosen to illustrate the behavior of the multi-order case, which exemplifies how various orders can be chosen to better fit certain data. The approximations are computed in the same way as before, up to x 99 . The computation time for the second example averaged 7.819 s. The results of the second example are shown in Figure 2. Concluding Remarks: In this work, we developed Picard's method to prove the existence and uniqueness of the solution to multi-order Caputo systems with initial conditions. The advantage of the multi-order system is that the order of the derivatives can be used as a parameter to make the solution closer to the data available for any mathematical model compared with the integer order or Caputo fractional system of the same order. We note that if the goal is to find the best fit for the approximations, i.e., the ideal orders, then, unfortunately, this can only be done via trial and error at this time. Determining an analytic way to find the optimal order is still an open problem and something we wish to study in the future.
In [26,27], the authors studied sequential Caputo fractional differential equations and their applications. It should be noted that the solution to the sequential Caputo differential equations can be reduced to a Caputo fractional differential system. The Caputo fractional linear system with constant coefficients and initial conditions can be solved using the Laplace transform method. The explicit and numerical computation of the linear system is essential to solving the corresponding nonlinear system by any iterative method. We plan to explore the explicit and numerical computation of the Caputo multi-order linear system. This will pave the way to developing iterative methods to solve nonlinear systems combined with upper and lower solutions.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: SIR Susceptible, Infected, and Recovered