On the Fractional Order Rodrigues Formula for the Shifted Legendre-Type Matrix Polynomials

: The generalization of Rodrigues’ formula for orthogonal matrix polynomials has attracted the attention of many researchers. This generalization provides new integral and differential representations in addition to new mathematical results that are useful in theoretical and numerical computations. Using a recently studied operational matrix for shifted Legendre polynomials with the variable coefﬁcients fractional differential equations, the present work introduces the shifted Legendre-type matrix polynomials of arbitrary (fractional) orders utilizing some Rodrigues matrix formulas. Many interesting mathematical properties of these matrix polynomials are investigated and reported in this paper, including recurrence relations, differential properties, hypergeometric function representation, and integral representation. Furthermore, the orthogonality property of these polynomials is examined in some particular cases. The developed results provide a matrix framework that generalizes and enhances the corresponding scalar version and introduces some new properties with proposed applications. Some of these applications are explored in the present work.


Introduction
The classic subject of special functions has emerged from a wide variety of practical problems that interest not only mathematicians but also other researchers in science to study their properties, characteristics, and applications.
Traditionally, the special functions of mathematical physics are defined using power series representations. Many orthogonal polynomials are then generated by imposing different termination conditions on these representations [1]. Another standard method to generate a sequence of orthogonal polynomials is to use the Rodrigues' formula (see, e.g., [1,2]). This method allows the attainment of many mathematical properties for these polynomials. For that reason, the generalization of the Rodrigues' formulas for orthogonal matrix polynomials has occupied the researchers' attention for the past few decades (see, e.g., [3] and the references therein) Special functions of fractional calculus, classified as generalized fractional integrals or derivatives, currently receive growing attention due to its enormous applications in various disciplines of science and engineering. For the historical development of fractional calculus, the reader is referred to the The fractional calculus, as a generalization of integer order integration and differentiation to its non-integer (fractional) orders counterpart, has been proved to be a valuable tool in the modeling of many phenomena in the fields of physics, chemistry, engineering, aerodynamic electrodynamics of complex medium, and polymer rheology [32]. This tool allows for describing many problems more accurately than in its classical setting. It is well-known that the integer-order differential operator is a local operator, but the fractional-order differential operator is non-local. The nonlocality property indicates that the next state of the system depends not only on its current state but also on all of its historical states. This makes studying fractional-order systems an active area of research.
The purpose of the present work is to introduce and examine a new fractional-order Rodrigues' formula for the shifted Legendre type matrix polynomials and functions starting by setting up the analog of Rodrigues' formula in matrix framework. Using this established formula, we investigate several properties of these polynomials including the recurrence relations, differential relations, integral representation, and hypergeometric representation. To our knowledge, the matrix framework for the shifted Legendre-type polynomials is reported here for the first time, especially in the fractional setting. The current study may open the door for further investigations concerning the practical applications [32]. By adopting a matrix setting with the notion of the Hermitian matrices, we also introduce a matricial inner product on C N×N that enables us to investigate the orthogonality property of these polynomials under certain constraints. The given results develop and generalize those of [17,33] by placing it in the matrix framework. In addition, we introduce some new results concerning fractional differentiation and integral representations. This study will provide a tool in developing a suitable technique to solve some multi-order fractional differential equations related to the boundary-value problems. More precisely and owing to the work [32,[34][35][36], one can use the presented approach to give an estimation of complex functions in terms of shifted Legendre-type matrix polynomials and then applying it to some problems mentioned in [32].
The paper is organized as follows. In the next section, some basic notations and some algebraic tools that are essential for the present work are introduced. In Section 3, we obtain the recurrence relation and establish some generalizations. In Section 4, the different differential formulas of the polynomials are achieved. In Section 5, a new representation of these matrix functions in terms of the Gauss hypergeometric matrix function is given. A new integral representation is developed in Section 6, and the orthogonality condition is proved in Section 7. Some applications are explored in Section 8. The conclusions and further suggested practical problems are recorded in Section 9. We should also mention that the present work is mostly analytic and designed to develop new properties with proposed algorithm, which are necessary for the forthcoming applications.

Preliminaries
Throughout this paper, C N×N is a vector space of all N × N square matrices with complex number entries z where the real and imaginary parts are Re(z) and Im(z), respectively. For any matrix A in C N×N , the spectrum σ(A) is the set of all eigenvalues of A, for which we set where α(A) refer to the spectral abscissa of A and for which α(−A) = −β(A). By definition, a matrix A ∈ C N×N is said to be positive stable if and only if β(A) > 0. The gamma and the beta matrix function are defined, respectively [26], as follows: and Following the authors of [24,26], we always consider in the sequel that, for any matrix A ∈ C N×N , the following Pochhammer symbol, or shifted factorial, denoted by (A) n is where Γ −1 (A) is a well defined matrix, and invertible (see [23], p. 273) as well as Note that [2], if A = −jI and j is a positive integer, then (A) n = 0 whenever n > j . Further, by definition, if A − nI is invertible for all integers n ≥ 0, In accordance with [24], if A and B are commuting positive stable matrices in C N×N , then the Beta matrix function defined in Equation (3) will satisfy B(A, B) = B(B, A); besides, if A + nI, B + nI, and A + B + nI are invertible for all integers n ≥ 0, then Definition 1 ([24]). For the positive integers p and q, the generalized hypergeometric matrix function is defined by the matrix power series where A i , 1 ≤ i ≤ p and B j , 1 ≤ j ≤ q in C N×N are commutative matrices and B j + n I are invertible for all integer n ≥ 0.
In the case of p = 2 and q = 1, Equation (7) reduces to the Gauss hypergeometric matrix function [24] Assume that B 1 A 2 = A 2 B 1 and B 1 , A 2 , B 1 − A 2 are a positive, stable, integral representation of the Gauss hypergeometric matrix function. for |z| < 1, the following is given by ([24], Theorem 5, p. 216): The fractional order integral and derivatives of Riemann-Liouville operator of order µ ∈ C, Re(µ) > 0, are given (see [14,37,38]), respectively, as (11) and Definition 2 ([19]). Suppose that A ∈ C N×N is a positive stable matrix with µ ∈ C and Re(µ) > 0; the Riemann-Liouville fractional integrals of order µ is defined by Lemma 1 ([19]). Suppose that A ∈ C N×N is a positive stable matrix with µ ∈ C and Re(µ) > 0; the Riemann-Liouville fractional integrals of order µ is given by This Lemma can be extended to: Lemma 2. Let A be a positive stable matrix in C N×N , and let µ, λ ∈ C such that Re(µ − λ) > −1.
Then, the Riemann-Liouville fractional integrals of order µ can be written by Definition 3 ([12,39,40]). The Caputo type fractional derivative of a function f (z) of order ν ∈ R + is defined by: where ν is the smallest integer greater than or equal to ν and ε ν ≥ 0, refer to the difference ε ν = ν − ν.
If (D ν f )(z) and (D ν g)(z) exist for the differential functions f (z) and g(z), it is straightforward to show that the Caputo type fractional derivative is a linear operator. For a detailed survey of the fractional derivatives and integrals, we refer the interested reader to [40].
Following Rajkovic and Kiryakova's [17] approach, we consider the functions of shifted Legendre type in the matrix setting by means of Rodrigues' formula as follows: If A, B ∈ C N×N , and satisfy the spectral conditions then for the fractional order ν, the shifted Legendre matrix functions given by the Rodrigues' formula as where D ν z A = Γ(A+I) Γ(A−(ν−1)I) z A−νI (in the Caputo sense). Noting that, for N = 1, taking A = B = ν = n ∈ N 0 , these matrix functions reduce [17] to the shifted Legendre polynomials {L n (z; 0, 1)} that are orthogonal on the interval (0, 1).
The generalization of the scalar case of recurrence relations of the fractional shifted Legendre-type function into the matrix framework is given in the following result.

The Recurrence Relations
Theorem 1. The shifted Legendre matrix function {P ν (A, B, z)} given by Equation (18) with the assumption in Equation (17) satisfies the following recurrence relations Proof. Relying on Equation (18) with the assumption in Equation (17) on A and B and using the linearity of the Caputo derivative operator, the assertion in Equation (19) follows from Similarly, using the assertion in Equation (20) follows directly using Equation (18). Finally, using the composition property of derivative and Equation (18) again, it follows that: which complete the proof of the theorem. (18) with the assumption in Equation (17), then

Corollary 1. Let A and B be matrices in C N×N satisfying Equation
Proof. By induction on n, shifting B to B + I, Equation (19) can be written as which proves Equation (22) when n = 1. Assume that Equation (22) is true for n = , that is to say Using the binomial identity ( +1 k ) = ( k ) + ( k−1 ), one can write for n = + 1 that as required. Here, we use the well-known formula for gamma function 1/Γ(−m) = 0, m = 0, 1, 2, · · · . Equation (23) follows similarly. For Equation (24), we note that which complete the proof of the Corollary.

Differential Relations
Differentiation of the shifted Legendre matrix functions given by the Rodrigues' formula has the following interesting properties.

Hypergeometric Representation
In this section, we present two interesting results as follows.
Theorem 3. The matrix function P ν (A, B, z) has the power series representation by means of the Gauss hypergeometric matrix function in the form: where n = ν and ε = n − ν.
Proof. Since and, similarly, By means of Leibniz formula, we have However, using the Binomial theorem [20], and Lemma 1, we have Thus, as required.

Remark 1. Using the identities
and Equation (30) can be written as Remark 2. When ν = n, a nonnegative integer, ε = 0 and, However, by means of the identity it follow that The following lemma is useful in the forthcoming result.

Theorem 4.
For nonnegative integer n, and A, B ∈ C N×N satisfying the assumption in Equation (17), if ν = n, Equation (18) yields the representation Proof. For ν = n, ε = 0 and due to lemma 3, Equation (36) can be written as Using the following identities we have and the conclusion follows.

Integral Representation
Theorem 5. The matrix function P ν (A, B, z) given by Equation (18) with assumption in Equation (17) have the following integral representation, by means of the hypergeometric function 2 F 1 , where n = ν and ε = n − ν > 0.
Proof. Using the integral representation in Equation (10) Now, applying Equation (30), it follows that from which we can write However, it is straightforward to show, using the series representation of 2 F 1 , that Thus, which complete the proof of the theorem. A similar technique in proving these results has been applied in [41].

Corollary 2.
For n ∈ N, the matrix function P n (A, B, z) given by Equation (18) with assumption in Equation (17) have the following integral representation, by means of the hypergeometric function 2 F 1 , Proof. The proof follows by shifting ε to ε + 1 in Equation (38), and then taking the limit as ε → 0.

Remark 4.
For positive integer n, and A, B ∈ C N×N satisfying the assumption in Equation (17), we have, by means of Equations (37) and (40), for |z| < 1, the following identity (37),

Remark 5. By the representation in Equation
we have the following integral for the shifted Legendre type polynomials P n (A, B, z):

Orthogonality Relation
It is worth mentioning that orthogonal polynomials comprise an emerging field of study, with important results in both theory and applications continuing to appear in the literature. Therefore, development and generalization of orthogonality from the scalar case to matrix setting, for particular cases, are given through the following result. From Equation (33), we have where ε = n − ν and n = ν . The orthogonality is studied by means of From this, This equation can be written as Using the series representation of one of the 2 F 1 function and assuming the switch of the infinite sum with the integral is permitted, one obtains The definite integral can be computed, in terms of the beta function B, to yield In the case ν = n, ν = m, ε = ε = 0, it follow that Using the Pochhammer identities Finally, we obtain Conjecture. The shifted Legendre matrix functions {P n (A, B, z)} are orthogonal on (0, 1) in the following sense that: when n + m is odd.
Now, using the properties of the Gamma and Beta matrix functions, we conclude that the integral J gives that allows us to write where C j,k = (−1) k+j m j n k For these coefficients, we observe that, if we replace j by m − j and k by n − k, we get a relation between C j,k and C m−j,n−k as follows = (−1) m+n C j,k , j = 0, 1, . . . , m; k = 0, 1, . . . , n Note that, for m + n odd, we get, for instance, if m + n = 1, ∑ k=0 ∑ j=0 C j,k = C 0,0 + C 0,1 = 0 where C 0,0 = −C 0,1 . Similarly, we easily see that C 0,0 + C 1,0 = 0 where C 0,0 = −C 1,0 . This means that, taking m and n to be integers and letting m + n be odd, each summand C j,k in ∑ n k=0 ∑ m j=0 C j,k has its opposite signed companion, therefore ∑ n k=0 ∑ m j=0 C j,k = 0 and the orthogonality property is verified.

Some Applications
Nowadays, fractional differential equations (FDEs) have been stimulated due to their enormous applications in several problems in science and engineering. Among these applications, we mention viscoelasticity; bioengineering; dynamics of interfaces between nanoparticles and substrates; fluid-dynamic traffic; signal processing, solid mechanics, control theory, statistical mechanics, economics; and nonlinear oscillations of earthquakes (see [43] and the references therein for more details). In general, most FDEs do not have an exact solution, which is available only for a considerable small class of FDEs. Sometimes it is challenging to obtain the exact analytic solutions of FDEs. In some cases, it becomes impossible to arrive at the exact analytic solution. The reason for this difficulty is the computational complexities of fractional calculus involving in these equations. Therefore, finding numerical solutions has become the destination of many mathematicians. They introduced many methods such as Laplace transforms, power series methods, collocation methods, Tau method, Legendre operational matrix method, etc. (see [43], and reference therein).
Approaching the problems mentioned using the fractional shifted Legendre matrix polynomials (FSLMPs) may provide new insights on solving some of these problems. With the aid of the given condition of orthogonality and owing to the Legendre operational matrix method [43], a truncated matrix shifted the Legendre series together with the Legendre operational matrix of the fractional derivative are proposed. This approach allows us to deal with the numerical solution of a class of FDEs (see [7]).
As practical applications, we mention two types of problems: the first is the initial and boundary value problems for the fractional Bagley-Torvik differential equations [44] and the second possible application is the time fractional coupled Korteweg-de Vries (KdV) system (see [45]), which is considered to be one of the most important nonlinear evolution systems of equations for its several applications in physics and engineering (see [45] and references therein).

Technical Scheme
The explicit analytical form of the shifted Legendre matrix polynomials P n (A, B; z) of degree n may be written as: respectively.

Overview on Operational Matrix
Recently, various operational matrices for the polynomials have been developed to study the numerical solution of differential equations. The main advantage of these operational matrices is that they replace differential operators with some matrices. Consequently, they reduce such problems to analyzing a system of algebraic equations, thus greatly simplifying the problem (see [43] and references therein). This trend is out of the scope of the present paper. However, we proceed to explore a technical scheme. Thus, the derivative of the vector ∆ m (A, B; z) can be expressed by

Legendre Operational Matrix in Fractional Setting
It is clear using Equation (51) that where n ∈ N and the subscript, in D (1) , denotes matrix powers. Thus, The generalization of Theorem 1 of [7] leads to the following scheme, which can be used to solve multi-order fractional differential equations. Analog to the scalar case, we consider the generalized linear multi-order fractional differential equation in the form: with initial conditions where A s for s = 1, 2, . . . , k + 2 are real coefficients and n < α ≤ n + 1, 0 < β 1 < β 2 < · · · < β k < α and D α refers to the Caputo fractional derivative of order α. To solve the problem in Equations (54) and (55), it is required to approximate y(x) and g(x) by the shifted Legendre matrix polynomials as follows.
where the vector G = [g 0 , . . . , g m ] T is known but the vector C = [c 0 , . . . , c m ] T is unknown. Following similar procedure as in the scalar case [7] leads to the following systems.
Using Equations (56)-(59), the residual R m (z) for Equation (54) can be written as Owing to Tau method [46], it follows that m − n linear equations can be generated by using It is clear that Equations (61) and (62) generate (m − n) and (n + 1) set of linear equations, respectively. These linear equations can be solved for unknown coefficients of the vector C, and hence y(z) stated in Equation (56) can be established.

Remark 7.
For the sake of simplicity, we restrict ourselves to the case when N = 1 (order of the matrices A and B) and taking A = B = ν = n ∈ N 0 which give the the technique used in [7]; however, our proposed approach can be extended for application. The main characteristic behind the approach using this technique of Legendre operational matrix is that it reduces the given problem to those of solving a system of algebraic equations. Thus, it greatly simplifies the problem and this needs hard calculations which would be difficult to be recorded for the moment. Example 1. We consider the following initial value problem in the case of inhomogeneous Bagley-Torvik equation [47]. D 2 y(z) + D 3/2 y(z) + y(z) = 1 + z y(0) = 1, y (0) = 1 (63) This problem has the exact solution in the form y(z) = 1 + z. By applying the technical scheme described above with m = 2, A = B = ν = n ∈ N 0 and N = 1, one can approximate solution as y(z) = c 0 P 0 (z) + c 1 P 1 (z) + c 2 P 2 (z) = C T ∆ m (z) Acting by fractional derivative, we get which is the exact solution. It can be expected that the above problem in Equation (63) may be generalized with matrix parameters and the shifted Legendre matrix polynomials can be used to solve it. Such tools would have made light for the forthcoming work besides numerical analysis technique.

Example 2.
(The time-fractional coupled Kortewg-de Varies (KdV) system) The KdV equation is considered as one of the most important nonlinear evolution equations for its several applications in physics and engineering (see [45]). The presented shifted Legendre matrix polynomials can be considered as a basic function to solve the KdV equation for the two types homogeneous and inhomogeneous. The homogeneous problem is stated as with the initial condition u(x, 0) = 4c 2 e cx (1+e cx ) 2 , v(x, 0) = 4c 2 e cx (1+e cx ) 2 . The exact solution of this problem, for α 1 = α 2 = 1, is u(x, t) = v(x, t) = 4c 2 e c(x−c 2 t) (1 + e c(x−c 2 t) ) 2 .
their genuine thanks to the anonymous referees for valuable suggestions which improved the final manuscript. Finally, they are indebted to the assistant editor Syna Mu for managing the process of reviewing the paper.

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