Generalized squared remainder minimization method for solving multi-term fractional differential equations

Hashemia Mir Sajjad hashemi_math396@yahoo.com

University of Bonab, República Islámica del Irán

Inc Mustafa minc@firat.edu.tr

Firat University, Turquía

Hajikhah Somayeh hashemi_math396@yahoo.com

University of Bonab, República Islámica del Irán

Esta obra está bajo una Licencia Creative Commons Atribución 4.0 Internacional.

Fractional integration and differentiation are generalizations of integer-order calculus to noninteger ones. It is demonstrated in literature that fractional calculus can play a justifiable and beneficial role in the modeling of various phenomena e.g. science and engineering [2, 3, 6, 25, 29, 30, 35–37, 42]. The classical differential operators are defined as local operators, whereas the fractional differential operators are nonlocal. This significant property makes studying fractional differential equations an active area of research. There are a several number of definitions of fractional derivatives involving the kernel of the special functions, such as Mittag–Leffler function, Prabhakar function and so on. Some important ones are, Riemann–Liouville and Caputo fractional derivatives in [14, 41], Caputo–Fabrizio fractional derivative in [8, 13], and Atangana and Baleanu suggested another version of fractional-order derivative, which uses the generalized Mittag–Leffler function with strong memory as nonlocal and nonsingular kernel in [5]. Multi-point fractional differential equations appear in different types of visco-elastic damping [40]. The multi-point boundary conditions in the fractional differential equations can be understood as the controllers at the end points dissipate or add energy according to censors located at intermediate points. Model equations proposed so far are almost always linear. Therefore, we concentrate on the following multi-point fractional differential equation:

where j , j=1; : : : ; k+1 are real constant coefficients and 0 < B1 < B2 <_ _ _< Bk < . Here, we consider this equation with most usual fractional derivative in Caputo sense. Equation (1) permits us to describe the model more accurately than the classical integer equation. The nonlocality of Caputo fractional derivative means that the next state of a system depends not only upon its current state but also upon all of its historical states.

Bhrawy et al. [11] applied the spectral algorithm based on generalized Laguerre tau (GLT) method with generalized Laguerre–Gauss (GQ) and generalized Laguerre–Gauss– Radau (GRQ) quadrature methods for Eq. (1). The shifted Chebyshev spectral tau (SCT) method based on the integrals of shifted Chebyshev polynomials is utilized to construct the approximate solutions of such problems [12]. Spectral tau method combined with the shifted Chebyshev polynomials are proposed to solve the Eq. (1) in [15]. Three alternative decomposition approaches are introduced for the approximate solution of Eq. (1) by Ford et al. [16]. More discussion about the approximate solution of Eq. (1) can be found in [27, 34, 38].

The paper is organized as follows: In Section 2, we briefly discuss some necessary definitions and mathematical preliminaries of fractional calculus, which will be needed in the forthcoming sections. The third section deals with a generalization of squared remainder minimization (GSRM) method for the multi-point fractional differential equations. Section 4 is devoted to the study of convergence analysis for the proposed method. Some illustrative examples are investigated in Section 5. Conclusion remarks provide our final section.

In this section, we review some preliminaries and properties of well-known fractional derivatives and integrals for the purpose of acquainting with sufficient fractional calculus theory. Among various definitions of fractional derivatives e.g. Caputo–Fabrizio [4,13,19], Atangana-Baleanu [7,17,18] and conformable fractional derivative [20,21,26], we introduce two most commonly used definitions, namely, the Riemann–Liouville and Caputo derivatives [9, 14, 33]. Various analytical and numerical methods are utilized to consider the fractional differential equations e.g. [1, 10, 22–24, 31, 32, 39].

Definition 1

Let _ 2 R+. The operator J_ x+ 0 E L1[x0; xf ] defined by

is called the Riemann–Liouville fractional integral operator.

Now, we discuss the interchange of the Riemann–Liouville fractional integration and limit operation, which is useful in the convergence analysis of further discussed method.

Lemma 1

Let _ 2 R+. Assume that fung11 is a uniformly convergent sequence of functions on C[x0; xf ]. Then

Definition 2

Let _ 2 R+. The operator RLD_ x+ 0 defined by

is called the Riemann–Liouville fractional differential operator, where the ceiling function d_e denotes the smallest integer greater than or equal to ∞

Definition 3

Let _ 2 R+. The operator CD_ x+ 0 defined by

For problem (1), we specify the linear operator

Suppose ψ0(x), ψ1(x), . . . , ψn(x) are linearly independent functions on [x0, xf ] and Ψn = span ψ0(x), ψ1(x), . . . , ψn(x) . Let un(x) Ψn is the approximate solution of Eq. (1). Then there exist unknown constants c0, c1, . . . , cn such that

Obviously, the approximate solution un(x) needs to satisfy the following conditions:

Substituting (3) into Eq. (2) concludes

Were

For any xE[x0, xf ], if en(x) = un(x)Eu(x), then the nth-order remaining terms can be given by

Remark 2.

If the linearly independent functions are supposed as

then Eq. (5) becomes

The important point to note here is the minimization problem in the GSRM method by the fact (4). That is, the problem is minimizing (Oun)(x) in such a way that

To do this, we introduce the real functions

and

Therefore in order to find the unknown vector C = (c0, c1, . . . , cn), we have to solve the minimization problem

Remark 3.

If we set polynomials (6) as basis functions, then constrains in minimization

problem (7) become

We use the Lagrange-multiplier method [28] to minimize problem (7). By using this method, we solve the following system of algebraic equations2:

Where and Equivalently, we can define the vector

where

For simplicity of notations, if we use

then

Therefore, the linear system of equations (8) becomes

or, in the abstract form,

where

Therefore, our main goal is finding the unknown vector C from Eq. (9).

This section is a discussion about the convergence of the GSRM method for the multipoint

fractional differential equations of the form (1).

Theorem 1

Suppose that u(x) is the exact solution of fractional differential equation (1) defined on [x0, xf ], and un(x) is the corresponding approximate solution of problem given by the GSRM method. If there exists a polynomial pn(x) ∈ Pn[x0, xf ] such that for any x ∈ [x0, xf ], pn(x) → u(x), then

Proof. One can easily conclude that for any n ∈ N, we have

Therefore

Moreover, from the continuity of norms and pn(x) → u(x) we get

From Eqs. (1), (2), obviously, we have ( u)(x) = 0 and therefore ( u)(x) L2 [x0 ,xf ] = 0. Hence, from (10) and (11) we obtain

which completes the proof.

n this section, we present the results of the GSRM method on five test problems. We perform our computations using Maple 18 software with 30 digits.

Example 1. Consider the equation

with initial condition u(0) = 0. Exact solution of this equation is given by u(x) = x4 − x3/2. Figure 1 shows the absolute errors for various α by the GSRM method

Figure 1
Absolute errors of example 1 for (a) α = 0.01, (b) α = 0.5 and (c) α = 0.99

Table 1
Comparison of maximum absolute values for Example 1 and various differential orders

with polynomial bases as (6) and n = 4. In Table 1, the absolute errors obtained by present method are compared with GLT method [11] with N = 50 and SCT method [12] with N = 64. From this table it may be concluded that the GSRM method is more accurate than the mentioned approaches. The CPU time used in this example is 0.842, and condition number of coefficient matrix in (9) w.r.t. α = 0.01 is 2.61

Example 2. Let us consider the following fractional equation

with initial conditions u(0) = 0 when 0 < α ≤ 1 and u(0) = u′(0) = 0 when 1 < α ≤ 2. The reported exact solution of this fractional equation is u(x) = xα+1. In this example, we suppose Ψ3 = span xl/2+1, l = 0, . . . 3 . Maximum absolute error in the best case E 03 is reported w.r.t. α = 1.5 by Bhrawy et al. [11], whereas 5.93E 18 and 1.31E 21 are extracted by the present method for α = 0.5 and α = 1.5, respectively. High accuracy of the GSRM method for both choices of α can be found from Fig. 2. In Fig. 3, the approximate solutions obtained by the proposed method are plotted with respect to different values of fractional-order α. The CPU time used in this example w.r.t. α = 0.5 and α = 1.5 are 0.375s and 0.421s, respectively.

Figure 2
Absolute errors of example 2 for (a) α = 0.5 and (b) α = 1.5.

Figure 3
Approximate solutions of example 2 for various values of α.

Example 3. Consider the Bagley–Torvik equation

with initial conditions u(0) = 0 and u′(0) = 0. Exact solution of this equation is given by u(x) = x2. For this problem, we choose the polynomial linear independent functions with n = 2, and we obtain

By solving the resultant system, we get c0 = c1 = 0 and c2 = 1. Therefore we obtain the exact solution for this example by using the GSRM method. The best maximum absolute errors by GLT(GQ) and GLT(GRQ) reported in [11] w.r.t. N = 90, are 7.06E 06 and 1.75E 05, respectively. The Haar wavelet collocation method obtained the error 4.23E 22 for h = 1/512 in [38]. By using an operational matrix method [34] it can be found the maximum absolute error 1.89E 12. The best reached error value among the utilized methods in [16] is 1.54E 05, which is reported for method 2 w.r.t. h = 1/512. So, the present method is most reliable and powerful approach for this example

Example 4. Let us consider the equation

with initial conditions u(0) = u′(0) = 0. Exact solution of this equation is given by u(x) = x3. Numerical results will not be presented since the exact solution is achieved by choosing n = 3

Regarding Example 4 and in [16], the best result is attained with 512 steps, and the maximum absolute errors are 6.93E 05, 1.18E 04 and 3.10E 06 by using method 1, method 2 and method 3, respectively. Obtained results by GLT(GQ) and GLT(GRQ) methods [11] with N = 64 are 1.43E 05 and 1.80E 05, respectively. Moreover, in [38], the absolute error 1.86E 09 is reported by HWCM method, and 3.39E 13 is reported by SCT method in [15].

Example 5. As the last test problem, we consider

with initial conditions u(0) = u′(0) = u′′(0) = 0 and f (x) = x3/3 + 2x0.8/Γ(1.8) + 2x1.75/Γ(2.75) + 2x2.25/Γ(3.25). Exact solution of this equation is reported as u(x) = x3/3. By choosing n = 4 we obtain the absolute error plot of this equation in Fig. 4.

In [27], the maximum absolute error by the Haar wavelet operational matrix method is 1.12E 02, and reported error in [38] is 2.91E 03, w.r.t. N = 64. Whereas by the GSRM method, we obtain 4.99E 28. The CPU time used in this example is 0.811s, and condition number of coefficient matrix in (9) is 849.54

Figure 4
Absolute error of Example 5.

In the present paper, the squared remainder minimization method is developed to the multi-term fractional differential equations. A minimization problem is manifested and it considered by the Lagrange-multiplier method. Convergence of the GSRM method is theoretically proved. Five test problems are investigated. For some of the given examples, exact solutions by the present method are extracted. Accuracy and reliability of the GSRM method is revealed by the reported figures and tables.

The authors are pleased to acknowledge the helpful comments of

the referees on an earlier version of this paper.