A new stable collocation method for solving a class of nonlinear fractional delay differential equations

In this paper, a stable collocation method for solving the nonlinear fractional delay differential equations is proposed by constructing a new set of multiscale orthonormal bases of W2,01\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$W^{1}_{2,0}$\end{document}. Error estimations of approximate solutions are given and the highest convergence order can reach four in the sense of the norm of W2,01\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$W_{2,0}^{1}$\end{document}. To overcome the nonlinear condition, we make use of Newton’s method to transform the nonlinear equation into a sequence of linear equations. For the linear equations, a rigorous theory is given for obtaining their ε-approximate solutions by solving a system of equations or searching the minimum value. Stability analysis is also obtained. Some examples are discussed to illustrate the efficiency of the proposed method.


Introduction
Nowadays, fractional differential equations have been a hot topic in the field of differential equations for their widespread applications in many science fields [1,2]. Among them, fractional delay differential equations begin to arouse attentions of many researchers. These equations have also many applications in various areas such as control theory, biology, and economy [3,4]. Since some models have a great deal to do with past condition, the insertion of a time delay makes these models more realistic. Therefore, the development of theory and numerical algorithms about fractional delay differential equations is of importance. For example, Hu and Zhao [5] give the condition of asymptotical stability of nonlinear fractional system with distributed time delay by utilizing the function monotonous properties and the stability theorem of the fractional linear system. Pimenov et al. [6] use a BDF difference scheme based on approximations of Clenshaw and Curtis type to get the numerical solution of the following equation: For (1), Moghaddam and Mostaghim [4] use the fractional finite difference method to obtain its numerical solutions. Saeed et al. [7] draw on the steps method and Chebyshev wavelet method to solve the nonlinear fractional delay differential equation: and get their approximate solutions. However, to our best knowledge, there are few articles about the study of fractional delay differential equations especially with regard to nonlinear fractional delay differential equations.
Newton's iterative method is a very powerful tool to solve the nonlinear problems. Many researchers have been studying and generalizing Newton's iterative method and they use it to solve nonlinear problems. Deuflhard [8] in his monograph constructs adaptive Newton algorithms for some specific nonlinear problems. Krasnosel'skii et al. [9] in their monograph study the Newton-Kantorovich's method, give a modified Newton-Kantorovich's method and solve the problem of the choice of initial approximations. Xu et al. [10] use quasi-Newton's method to linearize a nonlinear operator equation, and so on. Theoretical analysis shows that the convergence order of Newton's iterative formula is order 2.
Motivated greatly by the abovementioned excellent works, in this paper, we deal with the nonlinear fractional delay differential (2) in the condition p = 1, that is the following equation: where f has continuous second order partial derivative, g(x) ∈ C[0, 1] and y 0 (x) ∈ C 2 [−τ, 0] which are known functions. τ > 0 is a constant delay. The fractional derivative is in the sense of Caputo and y(x) ∈ C 1 [−τ, 1] is the unknown function.
In this article, we develop a stable and effective collocation method to solve (3). The collocation method is one of the most efficient methods for obtaining accurate numerical solutions of differential equations including variable coefficients and nonlinear differential equations [11][12][13][14]. The stability of collocation methods has always been an important topic. At present, the definition of its stability is that when there are many collocation nodes, the resulting equations are not ill conditioned and the results are still valid. For a stable collocation method, we can improve the accuracy by increasing the number of approximation terms. Accordingly, it is particularly important to establish a high-precision and stable collocation method.
The choice of bases of a space is important for a collocation method. Approximate solutions with different accuracy can be derived by using different bases for the same equation. For obtaining higher accuracy solution of the nonlinear fractional differential equation, we construct a new set of multiscale orthonormal bases of W 1 2,0 , give error estimations of approximate solutions, and prove in Section 3 that the highest convergence order can reach four in the sense of W 1 2,0 . Noting that the problem of initial value selection of the Newton's iterative method has been well solved in [9], so in this paper we transform (3) into a list of linear equations by using Newton's iterative method. Then a new stable collocation method is proposed to solve these equations. Compared with [15], the final numerical experiments show that our method is better in dealing with this kind of equations.
The remainder of this paper is organized as follows: In Section 2, some relevant definitions and properties of the fractional calculus and the space W 1 2 and W 1 2,0 are introduced. In Section 3, we construct a set of multiscale orthonormal bases of W 1 2,0 and give an error estimation of the approximate solution. In Section 4, we construct the ε−approximate solution method and apply this method to solve the linear fractional delay differential equation. In Section 5, we analyze the stability of the ε−approximate solution method. In Section 6, we make use of Newton's method to transform the nonlinear equation into a sequence of linear equations. In Section 7, we give the algorithm implementation of Newton's iterative formula for solving the nonlinear fractional delay differential equation. In Section 8, three numerical examples are given to clarify the effectiveness of the algorithm. In the last section, the conclusions are prensented.

Preliminaries and notations
In this section, some preliminary results about the fractional calculus operators and reproducing kernel spaces are recalled [16][17][18].

Definition 2.1
The Riemann-Liouville (R-L) fractional integral operator J α 0 is given by

Definition 2.2
The Caputo fractional differential operator D α C is given by Theorem 2.1 Let u(x) be the exact solution of (3), then D α C u(x) ∈ C[0, 1]. Proof By the assumption, one has Since integrability of a function is equivalent to its absolute integrability, one has

Definition 2.3
Assume that n is a positive integer, then we denote The inner product (·, ·) 1 and the norm 1 of W 1 2 are given by It is easy to see that W 1 2 ⊂ C[0, 1]. Similar to [17], we can prove that W 1 2 is not only a Hilbert space, but also a reproducing kernel space with the reproducing kernel
Proof For any u(x) ∈ W 1 2 , we have Thus The following definition will be used in Section 6.

Definition 2.5 The inner product space
And the inner product and the norm of W α 2 are given by Remark 2.1 Similar to [19], one can prove that W α 2 is not only a Hilbert space, but also a reproducing kernel space.

Construction of the multiscale orthonormal bases of W 1 2,0 and error estimations
In this section, we construct the multiscale orthonormal bases of W 1 2,0 by the famous Legendre multiwavelets, and give error estimations of approximate solutions.
Define the cubic Legendre scaling functions in the interval [0, 1] [20], and then cubic Legendre multiwavelets are given as Proof For the proof one can refer to [21].

Remark 3.1 Noting that
and we can see that we call the Vanishing Moment of μ(x) is r + 1.

Property 3.1 The Vanishing Moment of
Proof By specific calculation, we can reach that So the conclusion holds.
Remark 3.2 According to Property 3.1, we obtain that the Vanishing Moment of

Property 3.2 For any
Let y(x) and y J (x) be an exact solution and an approximate solution of (3), respectively. Denote Proof Without loss of generality, we first consider |c 0 ik |.
Thus, we have The first integral of the above inequality Similarly, one can obtain that In the same way, one can derive that |c l ik | ≤ 2 −(i−1)(j −1/2) AM, l = 1, 2, 3. So the conclusion holds.
Proof Using Property 3.2 and Lemma 3.2, we can derive

Linear fractional delay differential equations
In this section, we construct the ε-approximate solution method for solving (3) with f being linear, i.e., and analyze the convergence order of the ε-approximate solutions.
are known functions and y(x) is an unknown function. τ > 0 is a constant delay. Similar to Theorem 2.1, Then (5) is equivalently transformed into . So we use the denotation D α C and J α 0 below.
then (6) is equivalent to Therefore, the solution y(x) of (5) can be obtained by We need the following Lemmas for our aims.
Proof It is easy to see that So we can derive Thus,

Lemma 4.2 Let α > 0. Then
Proof We can find that hold for any ω(x) ∈ W 1 2,0 . On the other hand, we obtain by using the integration by parts. Hence, we have 2,0 and J α 0 ω 1 ≤M α ω 1 .

Lemma 4.3 Let a(x)
holds.
Proof Using Lemmas 2.1 and 4.2, we have So the conclusion holds. (7) is a bounded linear operator from W 1
Proof Suppose ω(x) is the exact solution of (8). By Theorem 3.1, there exists a positive integer N such that for any n > N, there exists ω n (x) = n i=0 c i ψ i (x) such that So we can derive That is, ω * n (x) is an ε−approximate solution of (8).

Stability analysis
In this section, we consider the stability of our proposed method. Assume λ is a eigenvalue of the matrix G which is defined in the proof of Theorem 4.4, that is, there exists an X = [x 0 , x 1 , · · · , x n ] T ∈ R n+1 , X = 0, such that GX = λX. Hence, x j L(ψ j (x))) 1 .
Thus, we have i=0 is a set of orthonormal bases of W 1 2,0 , we have Hence, we can obtain that Thus, the largest eigenvalue λ max of G satisfies λ max L 2 . On the other hand, we prove that the least eigenvalue of G satisfies λ min ≥ 1 L −1 2 . Otherwise, there exists {u n }, u n 1 = 1 such that L(u n ) 1 < 1 L −1 . Denote ω n = L(u n ). Hence, ω n 1 < 1 L −1 and u n = L −1 (ω n ). So we can derive It is a contradiction. So we obtain λ min ≥ 1 L −1 2 . Therefore, we have That is, Cond(G) 2 ≤ ( L L −1 ) 2 , which implies that the spectral condition number [23] of the matrix G is bounded, and the algorithm is stable.

Nonlinear fractional delay differential equations
In this section, we solve (3) when f is nonlinear. In order to obtain high-accuracy approximate solutions, we employ F-derivative and Newton's iterative formula.
Similar to the linear case, we denote

Then, (3) can be written as
Hence, there is no harm in supposing that the equation to be solved is Let u be the exact solution of the above equation and assume that D α C u(x) = ω(x). At the beginning of Section 4, we have proved that u(x) = J α 0 ω(x) and ω ∈ W 1 2,0 . So u(x) ∈ W α 2 . Define an operator F : (14) is equivalent to the following equation (15), then

Theorem 6.1 Suppose that F is defined by
where F (u) refers to the Fréchet derivative [10].
According to the definition of F-derivative and property of Caputo derivatives, we have for y, z, p, q and u(x), u (x), u(x −τ ), u (x −τ ) for y 0 , z 0 , p 0 , q 0 in the above equation, we get Therefore, ∃M > 0, Squaring on both sides of the inequality, and the inequality h 2 1 ≤ 2 The Newton's iterative formula for solving F (u) = 0 defined by (16) is with initial selection u 0 = 0, here F (u) is defined by (17). Equation (18) can be transformed into which are linear fractional delay differential (5), so we can solve them by the ε−approximate solution method constructed in Section 4.

Remark 6.2
One can refer to [9] about the method of selecting a initial value and the convergence of Newton's iterative formula.

Algorithm implementation of Newton's iterative formula
In this section, we will concretely show algorithm implementation of Newton's iterative formula. Assume where {c k,i }, {d k,i } are unknown. It is not difficult to see that . (23) Using (17) and (20), we have F (u k )(u k+1 − u k )  Substituting the above equations (19), (23), and (24) into the (25), we solve the following equations to obtain {d k,i }: To sum up, the algorithm is as follows: Step 1 Homogenize the equation as (14).
Step 3 Input α. Input the number of approximate items n. Set c 0,i = 0.
Using the software Mathematica, we can get the approximate solution u k+1 (x). We don't need too much iterations to reach desired approximate solutions because of high efficiency of Newton's iterative method.

Remark 7.1
The algorithm can be extended to the case of m−1 < α ≤ m (m ∈ N + ) for the fractional delay differential (3).

Numerical examples
In this section, the algorithm presented above is applied to solve linear and nonlinear fractional delay differential equations. Three examples are considered to illustrate the where g(x) is chosen such that the exact solution is x 2.3 or x 4.1 . When the exact solution is x 2.3 , we solve the equation with n = 65, 257, 513, respectively, and obtain the corresponding approximate solutions. The relative errors between the exact solution and the approximate solutions are displayed in Fig. 1. When the exact solution of this example is chosen to be x 4.1 which is of higher smoothness, the relative errors are displayed in Fig. 2. One can see that the relative errors are decreased with the increase of n. The results show that the error is decreased fast when the smoothness of solution is increased. In Table 1, we provide some numerical results illustrating the fact that the convergence order can be improved when smoothness of the solution is improved. Even if the solution of the equation is less smooth, the computing error is still acceptable for the engineering. Therefore, the algorithm is stable, reliable, and adaptive.   In Table 2, we compare the infinity-norm of relative errors at the discrete points and estimates of convergence order for both solution cases. The results show that our method is still valid, stable, and adaptive when solving nonlinear problems.

Example 8.3
Consider the linear fractional differential equation with delay [15] The exact solution of this example is x 2 . The relative errors between the exact solution and the approximate solutions are displayed in Fig. 5 with n = 65, 257, 513, 1025, respectively. One can see that the relative errors are decreased with the increase of n . We compare the infinity-norm of absolute errors at the discrete points with that from Ref. [15] in Table 3. The numerical results show that our method has higher accuracy in this example and the computing errors may be acceptable for engineering.

Conclusion
In this paper, we construct a new stable collocation method for solving a kind of nonlinear fractional delay differential equations. More suitable multiscale orthonormal bases of W 1 2,0 are constructed, error estimations of approximate solutions are given, and the highest convergence order can reach four in the sense of the norm of W 1 2 . Newton's iterative formula is used to linearize the nonlinear equation, and for the obtained linear equations we develop an ε-approximate solution method based on multiscale orthonormal bases to solve them. A concrete algorithm implementation is given. Numerical examples show that compared with [15], the presented method is more accurate in dealing with this kind of equations. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.