Nonlinear higher order fractional terminal value problems

: Terminal value problems for systems of fractional di ﬀ erential equations are studied with an especial focus on higher-order systems. Discretized piecewise polynomial collocation methods are used for approximating the exact solution. This leads to solving a system of nonlinear equations. For solving such a system an iterative method with a required tolerance is introduced and analyzed. The existence of a unique solution is guaranteed with the aid of the ﬁxed point theorem. Order of convergence for the given numerical method is obtained. Numerical experiments are given to support theoretical results.


Introduction
An important tool for describing the natural process is a system of fractional differential equations [1,2]. An ordinary system of fractional differential equations (FDEs) is described by C 0 D α t (q(t)) = p(t, q(t)), t ∈ [0, a] (1.1) where a ∈ R, q : [0, a] → R ν is a ν-dimensional vector function and p : [0, a] × R ν → R ν . Here, C 0 D α t is the Caputo fractional derivative of order α > 0. D α q(t) = p(t, q(t)), t ∈ [0, a], q(a) = q a . (1.2) where D α is fractional derivative and q a is a given ν-dimensional vector. At the first glance, these systems seem to be well-posed with regular source functions [15,16]. Surprisingly, Cong and Tuan [17] revealed some counterexamples. The papers [14,18] pointed out that the well-posedness depends on the terminal value a. For larger a we may not obtain a unique solution. Surprisingly, TVPs for fractional differential equations of higher-order derivatives α ≥ 1 are not well-studied. An order α ∈ (1, 2] problem on an infinite interval have been studied in [19]. These problems have the form C 0 D α t q(t) = p(t, q(t), q (t)), t ∈ [0, M], q(∞) = q ∞ , where q (∞) ∈ R is a given number and C 0 D α t stands for Caputo derivative. This paper's important aim is to study the existence results and regularity of the higher-order terminal value problem for systems of FDEs. The other important aim of this paper is to introduce an analyzed high-order numerical method for solving such problems. Thus, the piecewise polynomial collocation method (PPCM) as a numerical solver is introduced and analyzed in detail.
The PPCMs have some superior advantages in solving differential or integral equations. For example, if we use polynomials collocation methods (known as spectral methods) we may encounter Runge-phenomena. The piecewise characteristic of the PPCMs avoids such probable divergence.
In comparison with Runge-Kutta methods, the coefficient and parameters of the given methods are obtained constructively. It is remarkable to mention that for ordinary differential equations PPCMs, for some piecewise polynomial spaces are equivalent to Runge-Kutta methods. Conclusively, having more parameters (collocation parameters, degree of polynomial space, step size, meshing method) for controlling error, complexity and order of convergence makes the PPCMs competitively superior category of methods.
This paper is organized as follows: In Sections 2 and 3 we obtain the existence and uniqueness of a mild solution for higher-order FDEs. The regularity for the linear case is studied in Section 4. A numerical method is introduced in Section 5 and analysand in Section 6. Supportive numerical experiments are given in Section 7.

Higher order terminal value problem
where n − 1 < α < n and n ∈ N. The ν-dimensional vector function q = [q 1 , . . . , q ν ] is an unknown vector and q (i) a ∈ R ν (i = 0, . . . , n − 1) are known vectors. Further, we impose the Lipschitz condition on components of the function p : where L i (i = 0, . . . , ν) are constants and are not depend on t. Let z = q (n−1) , (z = q (n) ). Taking repeatedly integrals and using terminal values, we obtain By definition of Caputo derivative, we have Acting Riemann-Liouville integral 0 J α−n+1 t on both sides of Eq (1.1), we obtain Finally, putting t = a and using terminal conditions, we obtain (2.4) Remark 1. The solution of the coupled system of (2.2) and (2.4) can be regarded as a mild solution of the system (2.1).

Existence of a unique solution
Since we use vector functions, we need a norm combined with vector norm and function norm. This combinations brings us to face complexity in calculating induced norm. In this respect, the max norm defined by seems to be more easier. We establish well-posedness of the inverse problems (1.1) and (2.1). We define the operators P and Q : (C[0, a]) 2ν → (C[0, a]) ν by the right hand sides of systems (2.2) and (2.4), i. e., It is straightforward to show that for w 1 and w 2 ∈ (C[0, a]) 2ν . We can write and where L M = max i=1,...,ν L i . Now, the Eqs (3.4) and (3.5) give Conclusively, we can estate the following theorem.

Regularity for linear FDEs
The regularity of the solution is important for analyzing numerical methods, especially finite difference methods and methods based on projection [20]. The regularity speaks about the order of differentiability of solutions. Usually, the regularity is investigated in the space C m [a, b], (see Corollary 1). But, some good functions such as contains such functions. This space is not Banach or complete. Therefore, we need to introduce another Banach space such that In this paper, we study regularity on the following weighted space C m,α (0, a] introduced in [13].
Definition 4.1. [13] Let 0 ≤ α < 1 and m ∈ N. Then q ∈ C m,α (0, a] if there exist functions q i ∈ C[0, a] for i = 0, · · · , m + 1 such that q = t α q 0 + q 1 and is a Banach space. Proof. Let q n ∈ C m,α (0, a], n ∈ N be a Cauchy sequence with norm . α,m . Thus, by the definition of C m,α (0, a] there exists functions q n,i in C[0, a] for i = 0, · · · , m + 1 such that q n (t) = t α q n,0 (t) + q n,1 (t) and From the definition of the norm . α,m , the sequences q n,i for a fixed i are Cauchy sequence in the Banach space C[0, a] and thus have a unique limit say and hence for all t ∈ (a, b]. is a system of linear FDEs. Here A is a given ν × ν dimensional matrix function and b ∈ (C[0, a]) ν is a given source function. To study linear systems we introduce the operators W 1 and W 2 : and Let us also define the vector valued functions G 2 and G 1 : and Then, the Eqs (2.2) and (2.4) can be written as an inhomogeneous system Proof. Since T is a combination of weakly singular integral operators, it is a compact linear operator on (C[0, a]) 2ν and (C m,α (0, a]) 2ν . Also, it is clear that System (25) has a unique solution on (C[0, a]) 2ν by Theorem 3.1. Therefore, the homogeneous system w − T(w) = 0 has the trivial null space in (C[0, a]) 2ν and thus (C m,α (0, a]) 2ν . Conclusively, alternative Fredholm theorem asserts that the system (4.7) has a unique solution on (C m,α (0, a]) 2ν .
Remark 2. Regularity of the nonlinear case needs further investigations.

Numerical method
Considering the memory, one of the best methods for solving the coupled systems (2.2)-(2.4) is to use collocation methods on piecewise polynomial spaces. To implement such methods we partition the solution interval [0, a] into sub-intervals σ i = [t i , t i+1 ], i = 0, . . . , N − 1, with length h i := t i+1 − t i where 0 = t 0 < . . . < t N = a are nodes of a chosen mesh (uniform or graded mesh) and N ∈ N. A graded mesh with exponent r ≥ 1 is described by Let 0 < c 1 < . . . < c m ≤ 1 (m ∈ N) be collocation parameters, t n,i = t n + c i h n (n = 0, . . . , N − 1) be collocation points and letq N (t) andẑ N (t) be approximate solutions. The restriction of approximate functions to the σ k is fully determined by Lagrange polynomials interpolation formulâ where Q k, j =q N (t n, j ), Z k, j =ẑ N (t k, j ), t k, j = t k + hc j and L j ( j = 1, . . . , m) are Lagrange polynomials of degree m − 1. The integrals in the operators P and Q of systems (3.1) and (3.2) can be discretized as: Similarly, we discretize integrals of the operator Q. We have By interpolating p(t l + sh l , q(t l + sh l )) on c i we obtain and thus The discretized mappings P N and Q N : (C[0, a]) 2ν → (PC[0, a]) ν (PC stands for piecewise continuous space) related to systems (3.1) and (3.2) can be defined by and the dense solution can be evaluated in all points of the desired interval by Eqs (5.1) and (5.2). We note that the six parameters γ n,i,l, j , η n,l, j ,η n,k, j , n,l, j , ζ n,l, j andζ n,k, j fully determines the method for each m.
The equations of these six parameters for m = 1 and m = 2 are provided in Tables 1 and 2, respectively. Table 1. The coefficient of the collocation method of order m = 1. Table 2. The coefficient of the collocation method of order m = 2. Here, we used k = t k + vh k − t l and Θ = a − t l for simplifying the presentation of formulas.

Numerical analysis
For simplifying our analysis we recall some notations. Let S −1 m−1 (I h ) be the space of piecewise polynomials of degree less than m, (m ∈ N) on the mesh partitioning I h = {t i : i = 0, · · · , N}. The projection operator Π m−1,N : is uniquely determined by interpolation on collocation points such that Thus, Eq (5.10) can be written as A useful theorem regarding Π m−1,N that simplifies existence and convergence analysis is the following theorem.
Proof. In each sub-interval σ n of [0, a] by Lagrange interpolation formula, we have The rest of the proof is straightforward by taking the max norm.

Existence results for the proposed numerical method
One of the most fundamental questions is whether the system (6.1) has a unique solution? The answer is affirmative.
the operators P N of (5.8) and P of (3.1) are equivalent and we have P N (ŵ N (t)) = P(ŵ N (t)).
Thus, the inequality (3.4) holds for P N and for all w 1 , w 2 ∈ S −1 m−1 (I h ) 2ν . However, Q N is different from Q in this space and we need further computing. Actually, Q N can be described by N p(s, q(s)) (a − s) n−α ds where q is the first ν element of w. More precisely Therefore, Theorem 6.2. Let Then, the numerical solution of system (1.1) obtained by (6.1) exists and is unique.
is a contractor by (6.7). Thus, by using contraction mapping theorem, Π m−1,N T N admits a unique fixed-point.
holds, then Eq (6.8) holds too. For m = 1 we have Λ 1 = 1, and Eq (3) matches with Eq (3.7). However, our estimate may not be optimal for higher degrees of approximations. Also, we guess using the convergence properties of Π m−1,N when N → ∞, we may obtain a convergence result without dependency on Λ 1 .

Detailed implementation: Solving nonlinear system
An important question is how to solve the nonlinear Eq (6.1). There are many methods available in literature that we can employ. The Newton iteration method is one of them [14]. The advantage is that it is fast. The disadvantage is that it needs computation of Jacoby as well as a good initial value for an iteration. It increases competition cost for higher dimension (computation of ν 2 function for each iteration) and complexity of the implementation. However, strict restriction of λ Λ and analytic discussion of the previous section is constructive and suggests the use of an iterative method. Since the operator Π m−1,N T N is a contractor operator, beginning with an arbitrary initial value grantees the convergence of the related iterative method. In our algorithm, we initialize the iteration bŷ Then, we estimate the solution by the recursive formulâ The computation of iteration (6.10) continues with the smallest κ such that where T ol is a desired user tolerance. We report such κ as a number of iteration to achieve the given tolerance.

Convergence analysis
Proof. Let w be the solution of the coupled systems (2.2) and (2.4). Define the residue operator by and let e =ŵ N − w. Subtracting Eq (3.3) from (6.1), we obtain Taking maximum norm, we obtain Applying Eq (6.8), we obtain This is a good news since we can obtain the asymptotic behavior of R(w) ∞ for the known w. Thanks to Theorem 8 of [20], one can easily prove that which completes the proof.

Examples
All numerical experiments is implemented in MATLAB software. The tolerance number for solving related nonlinear equations is T ol = 10 −14 . This tolerance is close to the floating-point relative accuracy 2.2204e − 16 for data type of double-precision (64 bit) in standard machine. The numerical examples show theoretically the obtained order of convergence is sharp and can not be improved more. An estimated norm of the given method can be determined by (6.2) and we have Taking into account (6.8), we obtain λ Λ 1 = max{0.78663, 0.9} = 0.9.
Consequently, by using Theorem 6.2 the nonlinear Eq (6.1) has a unique solution and the iteration process (6.10) converges to that solution. In Tables 3 and 4, we provided the maximum error estimated order of proposed method (O N ) and the number of iteration for required tolerance. As we see O N = 0.5, 1, 1.5, 2 for r = 1, 2, 3, 4 respectively, which is in agreement with theoretical result (see Eq (7.2)).   Remark 4. Theorem 1 provides a sharp result when we consider all the components of the state of the system 1 (w = [q, z]). However, it can be improved for the desired state of the solution q. According to Corollary 1, q has better regularity for n ≥ 1. Thus, we guess that the order of convergence for q component with uniform mesh r = 1, can be greater than or equal to n − 1. Thus, for that component, a greater graded mesh exponent is not necessary. The error E q (r, N) = q N − q ∞ and the estimated order of convergence of Example 2 is reported in Table 5. The order 2 is achieved with r = 1.5.
Consequently, by using Theorem 6.2 the nonlinear Eq (6.1) related to m = 1 has a unique solution and the iteration process (6.10) converges to the exact solution. By Theorem 6.3 we expect Order = 0.5r, r ≤ 2, 1, r > 2. (7.4) The sharpness of this result can be seen in Table 6. Table 5. Numerical results of the proposed collocation method with c = [0, 1], m = 2 and r = 1, 2 for y component. Remark 5. The low order method have the advantage of supporting larger class of terminal value problem with large value of a. The reason is that with m = 1, we have Λ 1 = 1 which is three times smaller than Λ 1 for the case m = 2.

Conclusions
TVPs for higher-order (greater than one) fractional differential equations are rarely studied in the literature. In this paper, we tried to have a comprehensive study with a simple analysis to cover all general interests. Many questions in this topic deserve to study with more scrutiny. The regularity of nonlinear cases as well as optimal bound for obtaining well-posed problems are among them. We think the terminal value problem is important in applied since it monitors the past evolution of a dynamical system. Also, it can be regarded as a control problem for having the desired future by finding suitable initial value. We think this topic will catch more interest similar to the initial value problem.