An Operational Matrix Method Based on Poly-Bernoulli Polynomials for Solving Fractional Delay Differential Equations

: In this work, we derive the operational matrix using poly-Bernoulli polynomials. These polynomials generalize the Bernoulli polynomials using a generating function involving a polylogarithm function. We first show some new properties for these poly-Bernoulli polynomials; then we derive new operational matrix based on poly-Bernoulli polynomials for the Atangana–Baleanu derivative. A delay operational matrix based on poly-Bernoulli polynomials is derived. The error bound of this new method is shown. We applied this poly-Bernoulli operational matrix for solving fractional delay differential equations with variable coefficients. The numerical examples show that this method is easy to use and yet able to give accurate results.


Introduction
The classical fractional derivatives, such as Caputo, Riemann-Liouville and Grünwal, are based on singular kernels. On the other hand, recently, some new definitions of fractional derivative were derived based on nonsingular kernels. These derivatives include the Caputo-Fabrizio [1] and Atangana-Baleanu fractional derivatives [2]. These new definitions have been applied in many areas, including the groundwater flow within a confined aquifer [3], the magnetohydrodynamic electroosmotic flow of Maxwell fluids [4], the modeling of a financial system [5,6], the modeling of various types of diseases or epidemics [7][8][9] and many other fields [10,11].
On the other hand, solving a fractional differential equation which is defined as an Atangana-Baleanu fractional derivative always is not an easy task. Furthermore, an analytical solution may not always possible, especially when involving variable coefficients. Hence, we need some reliable numerical methods. To overcome this problem, in [12], the authors used fractional Adams-Bashforth methods to solve the fractional differential equation in a Atangana-Baleanu sense. The Crank-Nicholson difference method and reproducing kernel function had been applied for solving third order fractional differential equations in the sense of Atangana-Baleanu Caputo derivative [13]. An implicit, linear and unconditionally stable finite difference method was derived for solving fractional differential equations in the Atangana-Baleanu sense [14]. Other methods for solving special kinds of problems defined in the Atangana-Baleanu fractional derivative include the quasi wavelet approach of non-linear reaction diffusion and integro reaction-diffusion equation [15], the modified homotopy analysis transform method for solving gas dynamics equations of arbitrary order [16] and reproducing kernel functions [17].
In this research direction, operational matrix related methods are getting considerable attention for solving fractional calculus problems; among them is the Legendre operational matrix for solving fractional differential equations in the Caputo-Fabrizio sense [18]. Apart from that, for solving fractional differential equations defined in the Atangana-Baleanu fractional derivative, fifth-kind Chebyshev polynomials have been applied to derive operational matrices for multi-variable orders differential equations with non-singular kernels [19]. Besides that, the shifted Legendre cardinal functions operational matrix method was developed for solving the nonlinear time fractional Schrödinger equation with variable-order defined in the Atangana-Baleanu derivative [20]. Differently than the previous published results, here, for the first time, poly-Bernoulli polynomials which can be considered as a family of Bernoulli polynomials have been employed to derive the operational matrix. The poly-Bernoulli polynomials are derived using a generating function involving a polylogarithm function. Bernoulli polynomials have been used in solving various type fractional calculus problems, such as in [21,22]. These kinds of semi-orthogonal polynomials, including Genocchi polynomials and Euler polynomials, are widely used in solving differential equations up to arbitrary orders [23][24][25]. More specifically, we derive the new operational matrix based on poly-Bernoulli polynomials for solving the fractional differential equation in an Atangana-Baleanu sense. We also derive the delay operational matrix based on poly-Bernoulli polynomials to tackle the fractional delay differential equation in Atangana-Baleanu sense. In short, our proposed approach is able to overcome the difficulty of finding the solutions for variable coefficients of the fractional delay differential equation in an Atangana-Baleanu derivative via transforming the problem into solving a system of algebraic equations, which greatly reduces the complexity of the problem. This paper is organized as follows. We briefly explain some preliminarie concepts, including the Atangana-Baleanu fractional derivative and our main tool used in this work-poly-Bernoulli polynomials-in Section 2. Section 3 discusses the derivation of new operational matrix based on poly-Bernoulli polynomials for the Atangana-Baleanu fractional derivative. In Section 4, we explain the new scheme and give some examples for solving fractional delay differential equations defined in the Atangana-Baleanu fractional derivative using our propose method via the new operational matrix based on poly-Bernoulli polynomials. Conclusions and some recommendations are highlighted in Section 5.

The Atangana-Baleanu Derivative
In this subsection, we present some basic concepts related to the Atangana-Baleanu derivative.

Lemma 1.
The ABC derivative of order 0 < α < 1 for f (x) = x β can be defined as follows: Proof. By using Definition in Equations (1) and (2), knowing that , we obtain This expression is equivalent with those in [27] where the authors presented it in Mittag-Leffler form.

The Properties of Poly-Bernoulli Polynomials
According to [28][29][30], the poly-Bernoulli polynomials, B (k) n (x) can be obtained via using the generating function in Equation (6).
x n n k is the k th polylogarithm function. Recently, this polylogarithm function has been widely used in deriving new polynomials or number sequences such as in [31]. For x = 0, we obtain the poly-Bernoulli numbers, b By using generating function in Equation (6), the first few poly-Bernoulli polynomials of index k, n (x) are as follows: Some of the properties of Bernoulli polynomials can be extended to poly-Bernoulli polynomials, such as For function f (x) in an arbitrary element of L 2 [0, 1], one can approximate the function in terms of poly-Bernoulli polynomials as in Equation (9), where N (x)] and C or c n can be calculated using the following equation.
where denotes inner product and

The Poly-Bernoulli Delay Operational Matrix
The poly-Bernoulli delay vector B (k) (x − a) can be expressed as where R is the (N + 1) × (N + 1) delay operational matrix given by For any delay function f (x − a), one can approximate using poly-Bernoulli polynomials as: where C = [c 0 , c 1 , c 2 , · · · , c N ].

An Operational Matrix Based on Poly-Bernoulli Polynomials for an ABC-Fractional Derivative
In this section, we derive the new operational matrix based on poly-Bernoulli polynomials for the ABC-fractional derivative. This new operational matrix includes all the k positive integers for the family of poly-Bernoulli polynomials; when k = 1, one will obtain the Bernoulli operational matrix defined in an ABC-sense.
where P α is an (N + 1) × (N + 1) operational matrix of the fractional derivative of order α in an ABC-sense and is defined as follows: where θ n,r,j is given by: and b (k) n is the poly-Bernoulli number, c j can be obtained from inner product via Equation (10) and is the ceiling function. (8), we can write the poly-Bernoulli polynomials in analytical form and the fractional derivative is given as in Equation (17).

Proof. From Equation
The ABC D α x x r can be calculate using Lemma 1.
Let f (x) = x iα+r ; by using truncated poly-Bernoulli polynomials, we have f (

Error Bound
In this subsection, we present the error bound for the approximate ABC-derivative using an operational matrix based on poly-Bernoulli polynomials. The error bound basically follows the operational matrix using Genocchi polynomials as in [32].

Theorem 2. The error for approximating ABC D α
x f (x) by using poly-Bernoulli operational matrix, P α , as in Equation (14), is bounded as follows: Proof. Let us approximate the function f (x) by using infinite series poly-Bernoulli polynomials as follows: The ABC-derivative for the function f (x) can be approximated by using the operational matrix as shown in Equation (14).
Practically, we only use N + 1 terms of poly-Bernoulli polynomials to approximate the function; hence, Equation (24) will become This leads to approximate error as follows: Equation (26) can be written in matrix form as follows: If we consider the absolute error, we have

An Application to Solving Variable Coefficients of Fractional Delay Differential Equations in the ABC-Derivative
In this section, we apply the new poly-Bernoulli operational matrix and delay operational matrix for solving variable coefficients of a fractional differential equation in the ABC-derivative. We briefly explain the numerical scheme that we used. We achieved it by using the collocation scheme. All the numerical computations were carried out using Maple.

Collocation Scheme
We intended to solve the following variable coefficients of fractional delay differential equation in the ABC-derivative: with initial and boundary condition y(0) = y 0 and y(1) = y 1 . To solve Equation (29), let the solution of Equation (29), i.e., y(x), be approximated by the first N + 1 terms of poly-Bernoulli polynomials as follows: where the poly-Bernoulli coefficient vector C and the poly-Bernoulli vector B (k) (x) are given by The procedure for the collocation scheme is as follows: Step 1: Approximate each terms in Equation (29) by using poly-Bernoulli polynomials. For the derivative in an ABC-sense, i.e., ABC D α x y(x), we use Equations (14) and (15). The delay term y(x − a) is approximate via Equations (11) and (12).
Step 2: To find the solution y N (x); we collocate Equation (31) at the collocation points x j = j N+1 , j = 1, 2, · · · , N − 1 to obtain for j = 1, 2, · · · , N − 1. Together with initial and boundary condition, we obtain system of algebraic equations with N + 1 equations; by solving this system of equations, we obtain the value for poly-Bernoulli coefficient vector C.
Step 3: The solution for Equation (29) is obtained by substituting the poly-Bernoulli coefficient

Numerical Examples
In this subsection, some examples are presented to illustrate the applicability and accuracy of this poly-Bernoulli polynomial operational matrix and delay operational matrix for solving variable coefficients of the fractional delay differential equation defined in an ABC-derivative. All the computations were done by using Maple. For the first example, we solved a benchmark problem in the literature [33,34]. Example 1. Consider a simple fractional differential equation defined in the ABC-derivative as in [33,34].
The exact solution is given by y(x) = . As in [34], we consider M(α) This problem is solved by using a collocation scheme with N = 4 and by using poly-Bernoulli polynomials for k = 1 and k = 2. The absolute errors for the Example 1 are shown in Table 1. From the table, although the numerical scheme is simple with small N and easy to use, the solution is accurate. Here, we present the results when using different k (i.e., k = 1, 2) of poly-Bernoulli polynomials.
The exact solution is given by y( M(α)Γ(α)(α 3 +6α 2 +11α+6) . When α = 1 2 , the boundary condition is y(1) = 1 This problem is solved by using the collocation scheme with N = 4. By using poly-Bernoulli polynomials for k = 5, i.e., B (5) n (x), we obtain the approximate solution, as in Figure 1. The approximate solution is fixed on the exact solution. Example 3. Consider variable coefficients of the fractional delay differential equation defined in the ABC-derivative as follows.
The absolute errors for the Example 3 are shown in Table 2. This problem is solved by using the collocation scheme with N = 4. We present the results when using different k (i.e., k = 2, 5) of poly-Bernoulli polynomials. The solution is accurate with small absolute errors.

Conclusions
In this work, we achieved the following results: • A new operational matrix based on poly-Bernoulli polynomials for ABC-derivative.

•
A new delay operational matrix based on poly-Bernoulli polynomials.

•
A collocation scheme via operational matrix and delay operational matrix based on poly-Bernoulli polynomials for fractional differential equations in an ABC-sense.
The poly-Bernoulli operational matrix is the general case for any integer k, while when k = 1, it reduces to a Bernoulli operational matrix. The numerical examples show that the proposed method is a powerful tool for obtaining the numerical solutions for variable coefficients of fractional delay differential equations in an ABC-sense. Hence, we suggest that the following points may be able to apply:

•
The proposed scheme can be modified to solve various other problems that use an ABC-derivative.

•
The operational matrix can be extend to other poly-type polynomials.
In this work, our delay operational matrix was limited to constant delay. Hence, for future work, it would also be interesting to study a similar problem but with proportional delay or variable delay.