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

Abstract

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 coefficients 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.

1. 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 monograph [4]. Podlubny [5] wrote another book that focuses on the solution methods of the fractional differential equations and some applications. The reader can also find, in [6,7], other applications of the fractional calculus and fractional differential equations in model physical and engineering processes. Therein, the authors of [6] investigated the numerical solutions of non-linear multi-term of fractional differential equations, while the authors of [7] extended the applications of Legendre polynomials to solve different fractional differential equations. It is worth noting that the uses of the shifted Legendre polynomials of solving fractional differential equations proved to be much more effective. Hence, several research papers have been published recently to study the applications of shifted Legendre polynomials (see [8,9] and the cited references therein). Srivastava et al. [10] studied the global exponential stability of the fractional-order cellular neural network. Different applications of fractional calculus in engineering can be found in [11], where Tenreiro et al. explored fractional calculus uses in the field of controller tuning, legged robots, redundant robots, heat diffusion, and digital circuit synthesis.

The fractional-order analogs of the standard Rodrigues’ formulas were recently used to either define new classes of special functions or to study and classify the special functions of fractional calculus as generalized calculus operators of some basic functions [12,13,14]. Such an approach allows for establishing several new integral and differential representations that can be useful in applications and for examining related numerical algorithms. In this respect, a valuable use concerning fractional percolation on random graphs was introduced by Shang [15] in dealing with the vulnerability of networks. In [12,16], the authors gave several types of fractional generalizations of Rodrigues-type formulas for a class of special functions that deal with Laguerre polynomials of two fractional (arbitrary) orders. In [17], Rajković and Kiryakova used the fractional-order Rodrigues’ formula to investigate Legendre type special functions.

In the past few years, matrix functions and polynomials defined in terms of power matrix series [18,19,20,21] attracted the attention to new and exciting problems, specifically that involve the generalized hypergeometric matrix functions and future directions. It is worth mentioning the recent work of Shang [22] who uncovered some of the information of topological properties hidden in the smaller eigenvalues of A in matrix functions of the form $f\left(A\right)={\sum }_{k=0}^{\infty }{c}_k{A}^k.$ The graph-theoretic index introduced in [22] with the information encoded in the eigenvalues of the Hamiltonian near-zero plays a fundamental role in molecular magnetic/stability properties when A delineates the tight-binding Hamiltonian in the Hücked molecular orbited theory. This study provides some key tools to connect the function theory of approximation with graph theory.

Analogously to the classical case, the Rodrigues matrix formulas [23,24,25,26] are useful to define sequences of orthogonal matrix polynomials and the establishment of their mathematical properties. Furthermore, the Rodrigues matrix formulas allow defining new classes of special matrix functions and their fractional-order differentiation.

Recently, many authors (see, e.g., [18,19,20,21]) proposed the extension of the classical polynomials to the matrix framework. The operational matrix of fractional derivatives has been studied for various types of orthogonal polynomials, including Chebyshev polynomials [27]. Operational matrices for solving fractional-order differential equations can be found in [28]. In [29], the authors introduced a new operational matrix for shifted Legendre polynomials and fractional differential equations with variable coefficients. The authors of [30,31] also contributed to the study of the Legendre matrix polynomials and their properties. It is worth mentioning that, in [32], the authors used a shifted Legendre spectral method to deal with fractional order multi-point boundary value problems.

Now, owing to the significance of the earlier work regarding to Legendre type functions and other important special functions, one should record that many researchers became interested in analyzing the scalar cases of the classical special functions in the matrix setting. Much of these recent studies that signify the fractional setting and the matrix approach furnish a road to investigate the shifted Legendre matrix polynomials further using fractional-order Rodrigues matrix formula [17,29,33].

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 ${\mathbb{C}}^{N\times 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.

2. Preliminaries

Throughout this paper, ${\mathbb{C}}^{N\times N}$ is a vector space of all $N\times N$ square matrices with complex number entries z where the real and imaginary parts are $Re\left(z\right)$ and $Im\left(z\right)$, respectively. For any matrix A in ${\mathbb{C}}^{N\times N}$, the spectrum $\sigma \left(A\right)$ is the set of all eigenvalues of A, for which we set

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \alpha \left(A\right)=max\{Re(z):z\in \sigma (A\left)\right\}\mathrm{and}\beta \left(A\right)=min\{Re(z):z\in \sigma (A\left)\right\},\end{array}\end{array}$$

(1)

where $\alpha \left(A\right)$ refer to the spectral abscissa of A and for which $\alpha (-A)=-\beta \left(A\right)$. By definition, a matrix $A\in {\mathbb{C}}^{N\times N}$ is said to be positive stable if and only if $\beta \left(A\right)>0$.

The gamma and the beta matrix function are defined, respectively [26], as follows:

$$\Gamma \left(A\right)={\int }_0^{\infty }{e}^{-t}{t}^{A-I}dt;t^{A-I}=exp((A-I)lnt)$$

(2)

and

$$\mathbf{B}(A,B)={\int }_0^1{t}^{A-I}{(1-t)}^{B-I}dt.$$

(3)

Following the authors of [24,26], we always consider in the sequel that, for any matrix $A\in {\mathbb{C}}^{N\times N}$, the following Pochhammer symbol, or shifted factorial, denoted by ${\left(A\right)}_n$ is

$$\begin{array}{cc}\hfill {\left(A\right)}_0=I,{\left(A\right)}_n& =A(A+I)\dots (A+(n-1)I)=\Gamma (A+nI){\Gamma }^{-1}\left(A\right),n\ge 1,\hfill \end{array}$$

where ${\Gamma }^{-1}\left(A\right)$ is a well defined matrix, and invertible (see [23], p. 273) as well as

$$A+nI\mathrm{is}\mathrm{invertible}\mathrm{for}\mathrm{all}\mathrm{integers}n\ge 0.$$

(4)

Note that [2], if $A=-jI$ and j is a positive integer, then ${\left(A\right)}_n=\mathbf{0}$ whenever $n>j$. Further, by definition, if

$$A-nI\mathrm{is}\mathrm{invertible}\mathrm{for}\mathrm{all}\mathrm{integers}n\ge 0,$$

(5)

then

$$\begin{array}{cc}\hfill {(-A)}_n& ={(-1)}^{n}A(A-I)\cdots (A-(n-1)I)\Gamma (A-(n-1)I){\Gamma }^{-1}(A-(n-1)I)\hfill \\ \hfill & ={(-1)}^n\Gamma (A+I){\Gamma }^{-1}(A-(n-1)I).\hfill \end{array}$$

In accordance with [24], if A and B are commuting positive stable matrices in ${\mathbb{C}}^{N\times N}$, then the Beta matrix function defined in Equation (3) will satisfy $\mathbf{B}(A,B)=\mathbf{B}(B,A)$; besides, if $A+nI,$$B+nI$, and $A+B+nI$ are invertible for all integers $n\ge 0$, then

$$\mathbf{B}(A,B)=\Gamma \left(A\right)\Gamma \left(B\right){[\Gamma (A+B)]}^{-1}.$$

(6)

Definition 1

([24]). For the positive integers p and q, the generalized hypergeometric matrix function is defined by the matrix power series

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {}_p{F}_q(A_i;B_j;z)=\sum _{n\ge 0}\prod _{i=1}^p{\left(A_i\right)}_n{\left[\prod _{j=1}^q\left[{\left(B_j\right)}_n\right]\right]}^{-1}\frac{z^n}{n!},(p=q+1,|z|<1),\end{array}\end{array}$$

(7)

where $A_i,1\le i\le p$ and $B_j,1\le j\le q$ in ${\mathbb{C}}^{N\times N}$ are commutative matrices and

$$\begin{array}{c}\hfill B_j+nIareinvertibleforallintegern\ge 0.\end{array}$$

(8)

In the case of $p=2$ and $q=1$, Equation (7) reduces to the Gauss hypergeometric matrix function [24]

$$\begin{array}{c}\hfill {}_2{F}_1(A_1,A_2;B_1;z)=\sum _{n=0}^{\infty }{\left(A_1\right)}_n{\left(A_2\right)}_n{\left[{\left(B_1\right)}_n\right]}^{-1}\frac{z^n}{n!},\left|z\right|<1.\end{array}$$

(9)

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 $\left|z\right|<1$, the following is given by ([24], Theorem 5, p. 216):

$$\begin{array}{cc}\hfill {}_2{F}_1(A_1,A_2;B_1;z)=& \left({\int }_0^1{(1-zt)}^{-A_1}{t}^{A_2-I}{(1-t)}^{B_1-A_2-I}dt\right){\Gamma }^{-1}\left(A_2\right){\Gamma }^{-1}(B_1-A_2)\Gamma \left(B_1\right).\hfill \end{array}$$

(10)

The fractional order integral and derivatives of Riemann–Liouville operator of order $\mu \in \mathbb{C}$, $Re\left(\mu \right)>0$, are given (see [14,37,38]), respectively, as

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \left({\mathbf{I}}_{}^{\mu }f\right)\left(z\right)=\frac{1}{\Gamma \left(\mu \right)}{\int }_0^z{(z-t)}^{\mu -1}f\left(t\right)dt\end{array}\end{array}$$

(11)

and

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\mathbf{D}}_{}^{\mu }f\left(z\right)={\mathbf{I}}_{}^{n-\mu }{\mathbf{D}}^{n}f\left(z\right),\mu \in (n-1,n),n\in {\mathbb{Z}}_{+},Re\left(z\right)>0,\mathbf{D}=\frac{d}{dz}.\end{array}\end{array}$$

(12)

Definition 2

([19]). Suppose that $A\in {\mathbb{C}}^{N\times N}$ is a positive stable matrix with $\mu \in \mathbb{C}$ and $Re\left(\mu \right)>0$; the Riemann–Liouville fractional integrals of order μ is defined by

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\mathbf{I}}^{\mu }\left(z^A\right)=\frac{1}{\Gamma \left(\mu \right)}{\int }_0^z{(z-t)}^{(\mu -1)I}{t}^{A}dt,A\in {\mathbb{C}}^{N\times N}.\end{array}\end{array}$$

(13)

Lemma 1

([19]). Suppose that $A\in {\mathbb{C}}^{N\times N}$ is a positive stable matrix with $\mu \in \mathbb{C}$ and $Re\left(\mu \right)>0$; the Riemann–Liouville fractional integrals of order μ is given by

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\mathbf{I}}^{\mu }\left(z^{A-I}\right)=\Gamma \left(A\right){\Gamma }^{-1}(A+\mu I)z^{A+(\mu -1)I}.\end{array}\end{array}$$

(14)

This Lemma can be extended to:

Lemma 2.

Let A be a positive stable matrix in ${\mathbb{C}}^{N\times N}$, and let $\mu ,\lambda \in \mathbb{C}$ such that $Re(\mu -\lambda )>-1$. Then, the Riemann–Liouville fractional integrals of order μ can be written by

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\mathbf{I}}^{\mu }\left(z^{A-\lambda I}\right)=\Gamma (A-(\lambda -1)I){\Gamma }^{-1}(A+(\mu -\lambda +1)I)z^{A+(\mu -\lambda )I}.\end{array}\end{array}$$

(15)

Definition 3

([12,39,40]). The Caputo type fractional derivative of a function $f\left(z\right)$ of order $\nu \in {\mathbb{R}}^{+}$ is defined by:

$$\begin{array}{c}\hfill \left({\mathcal{D}}^{\nu }f\right)\left(z\right)=\left({\mathbf{I}}^{{\epsilon }_{\nu }}{D}^{\lceil \nu \rceil }f\right)\left(z\right),\end{array}$$

(16)

where $\lceil \nu \rceil$ is the smallest integer greater than or equal to ν and ${\epsilon }_{\nu }\ge 0$, refer to the difference ${\epsilon }_{\nu }=\lceil \nu \rceil -\nu .$

If $\left({\mathcal{D}}^{\nu }f\right)\left(z\right)$ and $\left({\mathcal{D}}^{\nu }g\right)\left(z\right)$ exist for the differential functions $f\left(z\right)$ and $g\left(z\right)$, 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\in {\mathbb{C}}^{N\times N}$, and satisfy the spectral conditions

$$Re\left(z\right)>-1\forall z\in \sigma \left(A\right),Re\left(z\right)>-1\forall z\in \sigma \left(B\right)\mathrm{and}AB=BA,$$

(17)

then for the fractional order $\nu ,$ the shifted Legendre matrix functions given by the Rodrigues’ formula as

$${\mathbf{P}}_{\nu }(A,B,z)=\frac{1}{\Gamma (\nu +1)}{\mathcal{D}}^{\nu }\left\{z^A{(1-z)}^B\right\}.$$

(18)

where ${\mathcal{D}}^{\nu }{z}^A=\frac{\Gamma (A+I)}{\Gamma (A-(\nu -1\left)I\right)}{z}^{A-\nu I}$ (in the Caputo sense). Noting that, for $N=1,$ taking $A=B=\nu =n\in {\mathbb{N}}_0,$ these matrix functions reduce [17] to the shifted Legendre polynomials $\left\{L_n(z;0,1)\right\}$ 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.

3. The Recurrence Relations

Theorem 1.

The shifted Legendre matrix function $\left\{{\mathit{P}}_{\nu }(A,B,z)\right\}$ given by Equation (18) with the assumption in Equation (17) satisfies the following recurrence relations

$$\begin{array}{cc}\hfill {\mathit{P}}_{\nu }(A,B,z)& ={\mathit{P}}_{\nu }(A,B-I,z)-{\mathit{P}}_{\nu }(A+I,B-I,z),\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill {\mathit{P}}_{\nu }(A,B,z)& ={\mathit{P}}_{\nu }(A-I,B,z)-{\mathit{P}}_{\nu }(A-I,B+I,z),\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill \nu {\mathit{P}}_{\nu }(A,B,z)& =A{\mathit{P}}_{\nu -1}(A-I,B,z)-B{\mathit{P}}_{\nu -1}(A,B-I,z).\hfill \end{array}$$

(21)

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

$$\begin{array}{cc}\hfill {\mathbf{P}}_{\nu }(A,B,z)& =\frac{1}{\Gamma (\nu +1)}{\mathcal{D}}^{\nu }\left(z^A{(1-z)}^{B-I}-z^{A+I}{(1-z)}^{B-I}\right).\hfill \end{array}$$

Similarly, using

$${\mathbf{P}}_{\nu }(A,B,z)=\frac{1}{\Gamma (\nu +1)}{\mathcal{D}}^{\nu }\left(z^{A-I}{(1-z)}^B-z^{A-I}{(1-z)}^{B+I}\right),$$

the assertion in Equation (20) follows directly using Equation (18). Finally, using the composition property of derivative and Equation (18) again, it follows that:

$$\begin{array}{c}\hfill {\mathbf{P}}_{\nu }(A,B,z)=\frac{1}{\Gamma (\nu +1)}{\mathcal{D}}^{\nu -1}\left(A{z}^{A-I}{(1-z)}^B-B{z}^A{(1-z)}^{B-I}\right),\end{array}$$

which complete the proof of the theorem. □

Corollary 1.

Let A and B be matrices in ${\mathbb{C}}^{N\times N}$ satisfying Equation (18) with the assumption in Equation (17), then

$$\begin{array}{cc}\hfill {\mathit{P}}_{\nu }(A+nI,B,z)& =\sum _{k=0}^n{(-1)}^k\left(\genfrac{}{}{0pt}{}{n}{k}\right){\mathit{P}}_{\nu }(A,B+kI,z),n\ge 0,\hfill \end{array}$$

(22)

$${\mathit{P}}_{\nu }(A,B+mI,z)=\sum _{k=0}^m{(-1)}^k\left(\genfrac{}{}{0pt}{}{m}{k}\right){\mathit{P}}_{\nu }(A+kI,B,z),m\ge 0,$$

(23)

$$\begin{array}{cc}\hfill {\mathit{P}}_{\nu }(A+nI,B+mI,z)& =\sum _{k=0}^n\sum _{j=0}^{m+k}{(-1)}^{k+j}\left(\genfrac{}{}{0pt}{}{n}{k}\right)\left(\genfrac{}{}{0pt}{}{m+k}{j}\right){\mathit{P}}_{\nu }(A+jI,B,z),\hfill \\ \hfill & =\sum _{k=0}^m\sum _{j=0}^{n+k}{(-1)}^{k+j}\left(\genfrac{}{}{0pt}{}{m}{k}\right)\left(\genfrac{}{}{0pt}{}{n+k}{j}\right){\mathit{P}}_{\nu }(A,B+jI,z),n,m\ge 0\hfill \end{array}$$

(24)

Proof. 

By induction on n, shifting B to $B+I$, Equation (19) can be written as

$${\mathbf{P}}_{\nu }(A+I,B,z)={\mathbf{P}}_{\nu }(A,B,z)-{\mathbf{P}}_{\nu }(A,B+I,z)$$

which proves Equation (22) when $n=1$. Assume that Equation (22) is true for $n=\ell$, that is to say

$${\mathbf{P}}_{\nu }(A+\ell I,B,z)=\sum _{k=0}^{\ell }{(-1)}^k\left(\genfrac{}{}{0pt}{}{\ell }{k}\right){\mathbf{P}}_{\nu }(A,B+kI,z)$$

Using the binomial identity $\left(\genfrac{}{}{0pt}{}{\ell +1}{k}\right)=\left(\genfrac{}{}{0pt}{}{\ell }{k}\right)+\left(\genfrac{}{}{0pt}{}{\ell }{k-1}\right)$, one can write for $n=\ell +1$ that

$$\begin{array}{cc}\hfill \sum _{k=0}^{\ell +1}{(-1)}^k\left(\genfrac{}{}{0pt}{}{\ell +1}{k}\right){\mathbf{P}}_{\nu }(A,B+kI,z)& =\sum _{k=0}^{\ell +1}{(-1)}^k\left(\genfrac{}{}{0pt}{}{\ell }{k}\right){\mathbf{P}}_{\nu }(A,B+kI,z)+\sum _{k=0}^{\ell +1}{(-1)}^k\left(\genfrac{}{}{0pt}{}{\ell }{k-1}\right){\mathbf{P}}_{\nu }(A,B+kI,z)\hfill \\ \hfill & =\sum _{k=0}^{\ell }{(-1)}^k\left(\genfrac{}{}{0pt}{}{\ell }{k}\right){\mathbf{P}}_{\nu }(A,B+kI,z)+\sum _{k=1}^{\ell +1}{(-1)}^k\left(\genfrac{}{}{0pt}{}{\ell }{k-1}\right){\mathbf{P}}_{\nu }(A,B+kI,z)\hfill \\ \hfill & =\sum _{k=0}^{\ell }{(-1)}^k\left(\genfrac{}{}{0pt}{}{\ell }{k}\right){\mathbf{P}}_{\nu }(A,B+kI,z)-\sum _{k=0}^{\ell }{(-1)}^k\left(\genfrac{}{}{0pt}{}{\ell }{k}\right){\mathbf{P}}_{\nu }(A,B+I+kI,z)\hfill \\ \hfill & ={\mathbf{P}}_{\nu }(A+\ell I,B,z)-{\mathbf{P}}_{\nu }(A+\ell I,B+I,z)\hfill \\ \hfill & ={\mathbf{P}}_{\nu }(A+(\ell +1)I,B,z),\hfill \end{array}$$

as required. Here, we use the well-known formula for gamma function $1/\Gamma (-m)=0,m=0,1,2,\cdots .$ Equation (23) follows similarly. For Equation (24), we note that

$$\begin{array}{cc}\hfill {\mathbf{P}}_{\nu }(A+nI,B+mI,z)& =\sum _{k=0}^n{(-1)}^k\left(\genfrac{}{}{0pt}{}{n}{k}\right)P_{\nu }(A,B+mI+kI,z)\hfill \\ \hfill & =\sum _{k=0}^n{(-1)}^k\left(\genfrac{}{}{0pt}{}{n}{k}\right)\sum _{j=0}^{m+k}{(-1)}^j\left(\genfrac{}{}{0pt}{}{m+k}{j}\right)P_{\nu }(A+jI,B,z),\hfill \end{array}$$

which complete the proof of the Corollary. □

4. Differential Relations

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

Theorem 2.

For $\nu ,{\nu }_1,{\nu }_2\in (-1,n-1)$, $n,m\in \mathbb{N}$, we have the following differential formulas:

$$\begin{array}{cc}\hfill D{\mathit{P}}_{\nu }(A,B,z)& =(\nu +1){\mathit{P}}_{\nu +1}(A,B,z),\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill D^n{\mathit{P}}_{\nu }(A,B,z)& ={(\nu +1)}_n{\mathit{P}}_{\nu +n}(A,B,z),\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill D^{{\nu }_1+n}{\mathit{P}}_{{\nu }_2}(A,B,z)& ={({\nu }_2+1)}_n{D}^{{\nu }_1}{\mathit{P}}_{{\nu }_2+n}(A,B,z),\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill {\mathcal{D}}^{{\nu }_1}{\mathit{P}}_{{\nu }_2}(A,B,z)& =\frac{\Gamma ({\nu }_1+1)}{\Gamma ({\nu }_2+1)}{D}^{{\nu }_2}{\mathit{P}}_{{\nu }_1}(A,B,z),\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill D^{{\nu }_1+n}{\mathit{P}}_{{\nu }_2+m}(A,B,z)& =\frac{{({\nu }_2+m+1)}_n}{{({\nu }_2+n+1)}_m}{D}^{{\nu }_1+m}{\mathit{P}}_{{\nu }_2+n}(A,B,z),\hfill \end{array}$$

(29)

Proof. 

Equation (25) follows from Equation (18) as

$${\mathbf{P}}_{\nu +1}(A,B,z)=\frac{1}{\Gamma (\nu +2)}{\mathcal{D}}^{\nu +1}\left\{z^A{(1-z)}^B\right\}=\frac{\Gamma (\nu +1)}{\Gamma (\nu +2)}D\left\{\frac{{\mathcal{D}}^{\nu }\left\{z^A{(1-z)}^B\right\}}{\Gamma (\nu +1)}\right\},$$

while Equation (26) follows by applying Equation (25) recursively n times. Equation (27) is easily deduced from Equation (26) by writing

$$D^{{\nu }_1+n}{\mathbf{P}}_{{\nu }_2}(A,B,z)=D^{{\nu }_1}\left(D^n{\mathbf{P}}_{{\nu }_2}(A,B,z)\right).$$

Equation (28) follows using Equation (18) (with $\nu$ is replaced with ${\nu }_2$) followed by reversing the order of ${\nu }_1$ and ${\nu }_2$, then applying Equation (18) once again (with $\nu$ is now replaced by ${\nu }_1$), thus

$$\begin{array}{cc}\hfill {\mathcal{D}}^{{\nu }_1}{\mathbf{P}}_{{\nu }_2}(A,B,z)& ={\mathcal{D}}^{{\nu }_1}\left\{\frac{{\mathcal{D}}^{{\nu }_2}\left\{z^A{(1-z)}^B\right\}}{\Gamma ({\nu }_2+1)}\right\}\hfill \\ \hfill & =\frac{{\mathcal{D}}^{{\nu }_1+{\nu }_2}\left\{z^A{(1-z)}^B\right\}}{\Gamma ({\nu }_2+1)}=\frac{\Gamma ({\nu }_1+1)}{\Gamma ({\nu }_2+1)}{\mathcal{D}}^{{\nu }_2}\left\{\frac{{\mathcal{D}}^{{\nu }_1}\left\{z^A{(1-z)}^B\right\}}{\Gamma ({\nu }_1+1)}\right\}.\hfill \end{array}$$

Finally, Equation (29) follow similarly. □

5. Hypergeometric Representation

In this section, we present two interesting results as follows.

Theorem 3.

The matrix function ${\mathbf{P}}_{\nu }(A,B,z)$ has the power series representation by means of the Gauss hypergeometric matrix function in the form:

$$\begin{array}{cc}\hfill {\mathit{P}}_{\nu }(A,B,z)=& {(-1)}^n\frac{\Gamma (A+I)\Gamma (B+I)}{\Gamma (\nu +1)}{z}^{A+\epsilon I}\sum _{k=0}^n\left(\genfrac{}{}{0pt}{}{n}{k}\right){\Gamma }^{-1}(A-(k-\epsilon -1)I){\Gamma }^{-1}(B-(n-k-1)I)\hfill \\ \hfill & \times {(-z)}^{-kI}{}_2{F}_1(-B+(n-k)I,A-(k-1)I;A+(\epsilon -k+1)I;z),\hfill \end{array}$$

(30)

where $n=\lceil \nu \rceil$ and $\epsilon =n-\nu .$

Proof. 

Since

$$D^n{z}^A=D^{n-1}\left(D{z}^A\right)=\Gamma (A+I){[\Gamma (A-(k-1)I)]}^{-1}{D}^{n-k}{z}^{A-kI},n\ge k,$$

we have, for $n=k$,

$$D^k{z}^A=\Gamma (A+I){[\Gamma (A-(k-1)I)]}^{-1}{z}^{A-kI},$$

and, similarly,

$$D^{n-k}{(1-z)}^B={(-1)}^{n-k}\Gamma (B+I){[\Gamma (B-(n-k-1)I)]}^{-1}{(1-z)}^{A-(n-k)I}.$$

By means of Leibniz formula, we have

$$\begin{array}{cc}\hfill {\mathbf{P}}_{\nu }(A,B,z)& =\frac{1}{\Gamma (\nu +1)}{\mathbf{I}}^{\epsilon }{D}^n\left\{z^A{(1-z)}^B\right\}\hfill \\ \hfill & =\frac{1}{\Gamma (\nu +1)}{\mathbf{I}}^{\epsilon }\sum _{k=0}^n\left(\genfrac{}{}{0pt}{}{n}{k}\right)\left(D^k{z}^A\right)\left(D^{n-k}{(1-z)}^B\right)\hfill \\ \hfill & =\frac{1}{\Gamma (\nu +1)}\sum _{k=0}^n{(-1)}^{n-k}\left(\genfrac{}{}{0pt}{}{n}{k}\right)\Gamma (A+I){[\Gamma (A-(k-1)I)]}^{-1}\hfill \\ \hfill & \times \Gamma (B+I){[\Gamma (B-(n-k-1)I)]}^{-1}{\mathbf{I}}^{\epsilon }\left[z^{A-kI}{(1-z)}^{B-(n-k)I}\right].\hfill \end{array}$$

However, using the Binomial theorem [20],

$${(1-z)}^{B-(n-k)I}=\sum _{s=0}^{\infty }\frac{{(-B+(n-k)I)}_s}{s!}{z}^s,$$

and Lemma 1, we have

$$\begin{array}{cc}\hfill {\mathbf{I}}^{\epsilon }\left(z^{A-kI}{(1-z)}^{B-(n+k)I}\right)& =\sum _{s=0}^{\infty }\frac{{(-B+(n-k)I)}_s}{s!}{\mathbf{I}}^{\epsilon }\left(z^{A-(k-s)I}\right)\hfill \\ \hfill & =\sum _{s=0}^{\infty }\frac{{(-B+(n-k)I)}_s}{s!}\hfill \\ \hfill & \times \Gamma (A-(k-s-1)I){\Gamma }^{-1}(A-(k-s-\epsilon -1)I)z^{A-(k-s-\epsilon )I}\hfill \\ \hfill & =\sum _{s=0}^{\infty }\frac{{(-B+(n-k)I)}_s}{s!}{(A-(k-1)I)}_s\Gamma (A-(k-1)I)\hfill \\ \hfill & \times {[{(A-(k-\epsilon -1)I)}_s\Gamma (A-(k-\epsilon -1)I)]}^{-1}{z}^{A-(k-s-\epsilon )I}\hfill \\ \hfill & =\Gamma (A-(k-1)I){\Gamma }^{-1}(A-(k-\epsilon -1)I)z^{A-(k-\epsilon )I}\hfill \\ \hfill & \times \sum _{s=0}^{\infty }\frac{{(-B+(n-k)I)}_s}{s!}{(A-(k-1)I)}_s{\left[{(A-(k-\epsilon -1)I)}_s\right]}^{-1}{z}^{sI}\hfill \\ \hfill & =\Gamma (A-(k-1)I){\Gamma }^{-1}(A-(k-\epsilon -1)I)z^{A-(k-\epsilon )I}\hfill \\ \hfill & \times {}_2{F}_1(-B+(n-k)I,A-(k-1)I;A-(k-\epsilon -1)I;z).\hfill \end{array}$$

Thus,

$$\begin{array}{cc}\hfill {\mathbf{P}}_{\nu }(A,B,z)& =\frac{\Gamma (A+I)\Gamma (B+I)}{\Gamma (\nu +1)}\sum _{k=0}^n{(-1)}^{n-k}\left(\genfrac{}{}{0pt}{}{n}{k}\right){\Gamma }^{-1}(B-(n-k-1)I){\Gamma }^{-1}(A-(k-\epsilon -1)I)z^{A-(k-\epsilon )I}\hfill \\ \hfill & \times {}_2{F}_1(-B+(n-k)I,A-(k-1)I;A-(k-\epsilon -1)I;z),\hfill \end{array}$$

as required. □

Remark 1.

Using the identities

$$\begin{array}{c}\hfill {\Gamma }^{-1}(A-(k-\epsilon -1)I)={(-1)}^k{\Gamma }^{-1}(A+(\epsilon +1)I){(-A-\epsilon )}_k\end{array}$$

(31)

and

$$\begin{array}{c}\hfill {\Gamma }^{-1}(B-(n-k-1)I)={(-1)}^n{\Gamma }^{-1}(B+I){(-B)}_n{\left[{(B-nI+I)}_k\right]}^{-1},\end{array}$$

(32)

Equation (30) can be written as

$$\begin{array}{cc}\hfill {\mathbf{P}}_{\nu }(A,B,z)& =\frac{\Gamma (A+I)}{\Gamma (\nu +1)}{(-B)}_n{\Gamma }^{-1}(A+(\epsilon +1)I)\hfill \\ \hfill & \times \sum _{k=0}^n\frac{{(-1)}^k{(-n)}_k}{k!}{\left[{(B-(n-1)I)}_k\right]}^{-1}{(-A-\epsilon )}_k{z}^{A-(k-\epsilon )I}\hfill \\ \hfill & \times {}_2{F}_1(-B+(n-k)I,A-(k-1)I;A-(k-\epsilon -1)I;z).\hfill \end{array}$$

(33)

Remark 2.

When $\nu =n$, a nonnegative integer, $\epsilon =0$ and, since ${}_2{F}_1(A,B;B;z)={}_1{F}_0(A;z)$, Equation (30) gives

$$\begin{array}{cc}\hfill {\mathit{P}}_n(A,B,z)=& \frac{\Gamma (A+I)\Gamma (B+I)}{\Gamma (n+1)}\sum _{k=0}^n{(-1)}^{n-k}\left(\genfrac{}{}{0pt}{}{n}{k}\right){\Gamma }^{-1}(A-(k-1)I){\Gamma }^{-1}(B-(n-k-1)I)\hfill \\ \hfill & \times z^{A-kI}{}_1{F}_0(-B+(n-k)I;z),\hfill \end{array}$$

(34)

However, by means of the identity

$$\begin{array}{c}\hfill \Gamma (B-(n-k-1)I)={(-1)}^n{(B-(n-1)I)}_k{\left[{(-B)}_n\right]}^{-1}\Gamma (B+I),\end{array}$$

(35)

it follow that

$$\begin{array}{c}\hfill {\mathit{P}}_n(A,B,z)=\frac{{(-B)}_n}{n!}\sum _{k=0}^n\frac{{(-1)}^k{(-n)}_k}{k!}{\left[{(B-(n-1)I)}_k\right]}^{-1}{(-A)}_k{z}^{A-kI}{}_1{F}_0(-B+(n-k)I;-;z).\end{array}$$

(36)

The following lemma is useful in the forthcoming result.

Lemma 3.

For$A\in {\mathbb{C}}^{N\times N}$,

$$\begin{array}{c}\hfill {}_1{F}_0\left(-A;z\right)={(1-z)}^A.\end{array}$$

Proof. 

Using Equation (7) with the matrix parameter $A\in {\mathbb{C}}^{N\times N},$ we get the following reduced hypergeometric function

$$\begin{array}{c}\hfill {}_1{F}_0\left(-A;z\right)=\sum _{n=0}^{\infty }{\left(-A\right)}_n\frac{z^n}{n!}.\end{array}$$

However, by the Pochhammer symbol definition ${\left(-A\right)}_n={(-1)}^n\Gamma (A+I){\Gamma }^{-1}(A-(n-1)),$ the lemma follows. □

Theorem 4.

For nonnegative integer n, and $A,B\in {\mathbb{C}}^{N\times N}$ satisfying the assumption in Equation (17), if $\nu =n$, Equation (18) yields the representation

$${\mathit{P}}_n(A,B,z)=\frac{1}{n!}{(-B)}_n{z}^A{(1-z)}^{B-nI}{}_2{F}_1\left(-nI,-A;B-(n-1)I;1-\frac{1}{z}\right).$$

(37)

Proof. 

For $\nu =n$, $\epsilon =0$ and due to lemma 3, Equation (36) can be written as

$$\begin{array}{cc}\hfill {\mathbf{P}}_n(A,B,z)=& {(-1)}^n\frac{\Gamma (A+I)\Gamma (B+I)}{n!}{z}^A{(1-z)}^{B-nI}\hfill \\ \hfill & \times \sum _{k=0}^n\left(\genfrac{}{}{0pt}{}{n}{k}\right){\Gamma }^{-1}(A-(k-1)I){\Gamma }^{-1}(B-(n-k-1)I){\left(1-\frac{1}{z}\right)}^{kI}.\hfill \end{array}$$

Using the following identities

$$\Gamma (B-(n-1)I+kI)=\Gamma (B-(n-1)I){(B-(n-1)I)}_k$$

and

$$\Gamma (A-kI+I)={(A+I)}_{-k}\Gamma (A+I)={(-1)}^k{\left[{(-A)}_k\right]}^{-1}\Gamma (A+I),$$

we have

$$\begin{array}{cc}\hfill {\mathbf{P}}_n& (A,B,z)={(-1)}^n\frac{\Gamma (A+I)\Gamma (B+I)}{n!}{z}^A{(1-z)}^{B-nI}\hfill \\ \hfill & \times \sum _{k=0}^n\left(\genfrac{}{}{0pt}{}{n}{k}\right){(-1)}^k{\Gamma }^{-1}(A+I){(-A)}_k{\Gamma }^{-1}(B-(n-1)I){\left[{(B-(n-1)I)}_k\right]}^{-1}{\left(1-\frac{1}{z}\right)}^{kI}\hfill \\ \hfill & ={(-1)}^n\frac{\Gamma (B+I){\Gamma }^{-1}(B+I-nI)}{n!}{z}^A{(1-z)}^{B-nI}\sum _{k=0}^n\left(\genfrac{}{}{0pt}{}{n}{k}\right){(-1)}^k{(-A)}_k{\left[{(B-(n-1)I)}_k\right]}^{-1}{\left(1-\frac{1}{z}\right)}^{kI}\hfill \\ \hfill & ={(-1)}^n\frac{\Gamma (B+I){\Gamma }^{-1}(B+I-nI)}{n!}{z}^A{(1-z)}^{B-nI}{}_2{F}_1(-nI,-A;B-(n-1)I;1-\frac{1}{z})\hfill \\ \hfill & ={(-1)}^n\frac{{\left[{(B+I)}_{-n}\right]}^{-1}}{n!}{z}^A{(1-z)}^{B-nI}{}_2{F}_1(-nI,-A;B-(n-1)I;1-\frac{1}{z})\hfill \\ \hfill & =\frac{{(-B)}_n}{n!}{z}^A{(1-z)}^{B-nI}{}_2{F}_1\left(-nI,-A;B-(n-1)I;1-\frac{1}{z}\right),\hfill \end{array}$$

and the conclusion follows. □

6. Integral Representation

Theorem 5.

The matrix function ${\mathit{P}}_{\nu }(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$,

$$\begin{array}{cc}\hfill {\mathit{P}}_{\nu }(A,B,z)& =\frac{{(-B)}_n{\Gamma }^{-1}\left(\epsilon \right)}{\Gamma (\nu +1)}{z}^{A+\epsilon I}\hfill \\ \hfill & \times {\int }_0^1{(1-tz)}^{B-nI}{t}^A{(1-t)}^{(\epsilon -1)I}{}_2{F}_1\left(-A,-nI;B-(n-1)I;1-\frac{1}{tz}\right)dt,\hfill \end{array}$$

(38)

where $n=\lceil \nu \rceil$ and $\epsilon =n-\nu >0.$

Proof. 

Using the integral representation in Equation (10), we have

$$\begin{array}{cc}\hfill {}_2{F}_1(-B+(n-k)I,A-(k-1)I;& A+(\epsilon -k+1)I;z)\hfill \\ \hfill & =\left({\int }_0^1{(1-tz)}^{B-(n-k)I}{t}^{A-kI}{(1-t)}^{(\epsilon -1)I}dt\right)\hfill \\ \hfill & \times {\Gamma }^{-1}(A-(k-1)I){\Gamma }^{-1}\left(\epsilon \right)\Gamma (A-(k-\epsilon -1)I).\hfill \end{array}$$

Now, applying Equation (30), it follows that

$$\begin{array}{cc}\hfill {\mathbf{P}}_{\nu }& (A,B,z)=\frac{\Gamma (A+I)\Gamma (B+I)}{\Gamma (\nu +1)}\sum _{k=0}^n{(-1)}^{n-k}\left(\genfrac{}{}{0pt}{}{n}{k}\right){\Gamma }^{-1}(B-(n-k-1)I)\hfill \\ \hfill & \times z^{A+(\epsilon -k)I}\left({\int }_0^1{(1-tz)}^{B-(n-k)I}{t}^{A-kI}{(1-t)}^{(\epsilon -1)I}dt\right){\Gamma }^{-1}(A-(k-1)I){\Gamma }^{-1}\left(\epsilon \right),\hfill \end{array}$$

from which we can write

$$\begin{array}{cc}\hfill {\mathbf{P}}_{\nu }(A,B,z)& ={(-1)}^n\frac{\Gamma (A+I)\Gamma (B+I){\Gamma }^{-1}\left(\epsilon \right)}{\Gamma (\nu +1)}{z}^{A+\epsilon I}{\int }_0^1{(1-tz)}^{B-nI}{t}^A{(1-t)}^{(\epsilon -1)I}\hfill \\ \hfill & \times \sum _{k=0}^n\left(\genfrac{}{}{0pt}{}{n}{k}\right){\Gamma }^{-1}(B-(n-k-1)I){\Gamma }^{-1}(A-(k-1)I){(-zt)}^{-kI}{(1-tz)}^{kI}dt.\hfill \end{array}$$

However, it is straightforward to show, using the series representation of ${}_2{F}_1$, that

$$\begin{array}{cc}\hfill \sum _{k=0}^n\left(\genfrac{}{}{0pt}{}{n}{k}\right){\Gamma }^{-1}(B-(n-k-1)I)& {\Gamma }^{-1}(A-(k-1)I){(-zt)}^{-kI}{(1-tz)}^{kI}\hfill \\ \hfill & ={\Gamma }^{-1}(A+I){\Gamma }^{-1}(B-(n-1)I){}_2{F}_1\left(-A,-n;B-n+1;1-\frac{1}{tz}\right).\hfill \end{array}$$

Thus,

$$\begin{array}{cc}\hfill {\mathbf{P}}_{\nu }(A,B,z)& ={(-1)}^n\frac{\Gamma (A+I)\Gamma (B+I){\Gamma }^{-1}\left(\epsilon \right)}{\Gamma (\nu +1)}{z}^{A+\epsilon I}{\int }_0^1{(1-tz)}^{B-nI}{t}^A{(1-t)}^{(\epsilon -1)I}\hfill \\ \hfill & \times {\Gamma }^{-1}(A+I){\Gamma }^{-1}(B-(n-1)I){}_2{F}_1\left(-A,-n;B-n+1;1-\frac{1}{tz}\right)dt\hfill \\ \hfill & =\frac{{(-B)}_n{\Gamma }^{-1}\left(\epsilon \right)}{\Gamma (\nu +1)}{z}^{A+\epsilon I}{\int }_0^1{(1-tz)}^{B-nI}{t}^A{(1-t)}^{(\epsilon -1)I}{}_2{F}_1\left(-A,-nI;B-(n-1)I;1-\frac{1}{tz}\right)dt\hfill \end{array}$$

which complete the proof of the theorem. A similar technique in proving these results has been applied in [41]. □

Corollary 2.

For $n\in \mathbb{N}$, the matrix function ${\mathit{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$,

$$\begin{array}{cc}\hfill {\mathit{P}}_{n-1}(A,B,z)& =\frac{{(-B)}_n}{\Gamma \left(n\right)}{z}^{A+I}{\int }_0^1{t}^A{(1-tz)}^{B-nI}{}_2{F}_1\left(-A,-nI;B-(n-1)I;1-\frac{1}{tz}\right)dt.\hfill \end{array}$$

(39)

Proof. 

The proof follows by shifting $\epsilon$ to $\epsilon +1$ in Equation (38), and then taking the limit as $\epsilon \to 0$. □

Remark 3.

The change of variable $tz=u$ allows writing Equation (39) as

$$\begin{array}{cc}\hfill {\mathit{P}}_{n-1}(A,B,z)& =\frac{{(-B)}_n}{\Gamma \left(n\right)}{\int }_0^z{u}^A{(1-u)}^{B-nI}{}_2{F}_1\left(-A,-nI;B-(n-1)I;1-\frac{1}{u}\right)du\hfill \end{array}$$

(40)

which confirms Equation (37) using the fundamental theorem of calculus and Equation (25).

Remark 4.

For positive integer n, and $A,B\in {\mathbb{C}}^{N\times N}$ satisfying the assumption in Equation (17), we have, by means of Equations (37) and (40), for $\left|z\right|<1$, the following identity

$$\begin{array}{cc}\hfill {\int }_0^z{u}^A& {(1-u)}^{B-nI}{}_2{F}_1\left(-A,-nI;B-(n-1)I;1-\frac{1}{u}\right)du\hfill \\ \hfill & ={(-B)}_{n-1}{\left[{(-B)}_n\right]}^{-1}{z}^A{(1-z)}^{B-(n-1)I}{}_2{F}_1\left(-A,-(n-1)I;B-(n-2)I;1-\frac{1}{z}\right).\hfill \end{array}$$

(41)

Remark 5.

By the representation in Equation (37),

$${\mathit{P}}_n(A,B,z)=\frac{1}{n!}{(-B)}_n{z}^A{(1-z)}^{B-nI}{}_2{F}_1\left(-nI,-A;B-(n-1)I;1-\frac{1}{z}\right),$$

(42)

we have the following integral for the shifted Legendre type polynomials ${\mathit{P}}_n(A,B,z)$:

$$\begin{array}{cc}\hfill {\int }_0^z{\mathit{P}}_n(A,B,u)du& =-{(B-(n-1)I)}^{-1}(1-z){\mathit{P}}_n(A,B,z).\hfill \end{array}$$

(43)

7. 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

$$\begin{array}{cc}\hfill {\mathbf{P}}_{\nu }(A,B,z)& =\frac{\Gamma (A+I){(-B)}_n{\Gamma }^{-1}(A+(\epsilon +1)I)}{\Gamma (\nu +1)}\sum _{k=0}^n\frac{{(-1)}^k{(-n)}_k{(-A-\epsilon )}_k{\left[{(B-(n-1)I)}_k\right]}^{-1}}{k!}\hfill \\ \hfill & \times z^{A-kI+\epsilon I}{}_2{F}_1(-B+(n-k)I,A-(k-1)I;A-(k-\epsilon -1)I;z)\hfill \end{array}$$

where $\epsilon =n-\nu$ and $n=\lceil \nu \rceil .$ The orthogonality is studied by means of

$${\int }_0^1{\mathbf{P}}_{\nu }(A,B,z){\mathbf{P}}_{{\nu }^{\prime }}(A,B,z)dz=I_{\nu {\nu }^{\prime }}.$$

From this,

$$\begin{array}{cc}\hfill I_{\nu {\nu }^{\prime }}& =\frac{\Gamma (A+I){(-B)}_n{\Gamma }^{-1}(A+(\epsilon +1)I)}{\Gamma (\nu +1)}\frac{\Gamma (A+I){(-B)}_m{\Gamma }^{-1}(A+({\epsilon }^{\prime }+1)I)}{\Gamma ({\nu }^{\prime }+1)}\hfill \\ \hfill & \times \sum _{k=0}^n\sum _{\ell =0}^m\frac{{(-1)}^{k+\ell }{(-n)}_k{(-A-\epsilon )}_k{\left[{(B-(n-1)I)}_k\right]}^{-1}}{k!}\frac{{(-m)}_{\ell }{(-A-{\epsilon }^{\prime })}_{\ell }{\left[{(B-(m-1)I)}_{\ell }\right]}^{-1}}{\ell !}\hfill \\ \hfill & \times {\int }_0^1{z}^{2A-(k+\ell )I+(\epsilon +{\epsilon }^{\prime })I}{}_2{F}_1\left(\begin{array}{cc}-B+(n-k)I,& A-(k-1)I\\ A-(k-\epsilon -1)I\end{array};z\right){}_2{F}_1\left(\begin{array}{cc}-B+(m-\ell )I,& A-(\ell -1)I\\ A-(\ell -{\epsilon }^{\prime }-1)I\end{array};z\right)dz\hfill \end{array}$$

This equation can be written as

$$\begin{array}{cc}\hfill I_{\nu {\nu }^{\prime }}& =\frac{\Gamma (A+I){(-B)}_n{\Gamma }^{-1}(A+(\epsilon +1)I)}{\Gamma (\nu +1)}\frac{\Gamma (A+I){(-B)}_m{\Gamma }^{-1}(A+({\epsilon }^{\prime }+1)I)}{\Gamma ({\nu }^{\prime }+1)}\hfill \\ \hfill & \times \sum _{k=0}^n\sum _{\ell =0}^m\frac{{(-1)}^{k+\ell }{(-n)}_k{(-A-\epsilon )}_k{\left[{(B-(n-1)I)}_k\right]}^{-1}}{k!}\frac{{(-m)}_{\ell }{(-A-{\epsilon }^{\prime })}_{\ell }{\left[{(B-(m-1)I)}_{\ell }\right]}^{-1}}{\ell !}\hfill \\ \hfill & \times {\int }_0^1{z}^{2A-(k+\ell )I+(\epsilon +{\epsilon }^{\prime })I}{(1-z)}^{2B+(k+\ell )I+(\epsilon +{\epsilon }^{\prime })I-(m+n)I}\hfill \\ \hfill & \times {}_2{F}_1\left(\begin{array}{cc}A+B+(\epsilon +1-n)I,& \epsilon \\ A-(k-\epsilon -1)I\end{array};z\right){}_2{F}_1\left(\begin{array}{cc}A+B+({\epsilon }^{\prime }+1-m)I,& {\epsilon }^{\prime }\\ A-(\ell -{\epsilon }^{\prime }-1)I\end{array};z\right)dz\hfill \end{array}$$

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

$$\begin{array}{cc}\hfill I_{\nu {\nu }^{\prime }}& =\frac{\Gamma (A+I){(-B)}_n{\Gamma }^{-1}(A+(\epsilon +1)I)}{\Gamma (\nu +1)}\frac{\Gamma (A+I){(-B)}_m{\Gamma }^{-1}(A+({\epsilon }^{\prime }+1)I)}{\Gamma ({\nu }^{\prime }+1)}\hfill \\ \hfill & \times \sum _{k=0}^n\sum _{\ell =0}^m\frac{{(-1)}^{k+\ell }{(-n)}_k{(-A-\epsilon )}_k{\left[{(B-(n-1)I)}_k\right]}^{-1}}{k!}\frac{{(-m)}_{\ell }{(-A-{\epsilon }^{\prime })}_{\ell }{\left[{(B-(m-1)I)}_{\ell }\right]}^{-1}}{\ell !}\hfill \\ \hfill & \times \sum _{s=0}^{\infty }\frac{{(A+B+({\epsilon }^{\prime }+1-m)I)}_s{\left({\epsilon }^{\prime }\right)}_s{\left[{(A-(\ell -{\epsilon }^{\prime }-1)I)}_s\right]}^{-1}}{s!}\hfill \\ \hfill & \times {\int }_0^1{z}^{2A-(k+\ell )I+(\epsilon +{\epsilon }^{\prime })I+sI}{(1-z)}^{2B+(k+\ell )I+(\epsilon +{\epsilon }^{\prime })I-(m+n)I}{}_2{F}_1\left(\begin{array}{cc}A+B+(\epsilon +1-n)I,& \epsilon \\ A+(\epsilon +1-k)I\end{array};z\right)dz\hfill \end{array}$$

The definite integral can be computed, in terms of the beta function B, to yield

$$\begin{array}{cc}\hfill I_{\nu {\nu }^{\prime }}& =\frac{\Gamma (A+I){(-B)}_n{\Gamma }^{-1}(A+(\epsilon +1)I)}{\Gamma (\nu +1)}\frac{\Gamma (A+I){(-B)}_m{\Gamma }^{-1}(A+({\epsilon }^{\prime }+1)I)}{\Gamma ({\nu }^{\prime }+1)}\hfill \\ \hfill & \times \sum _{k=0}^n\sum _{\ell =0}^m\frac{{(-1)}^{k+\ell }{(-n)}_k{(-A-\epsilon )}_k{\left[{(B-(n-1)I)}_k\right]}^{-1}}{k!}\frac{{(-m)}_{\ell }{(-A-{\epsilon }^{\prime })}_{\ell }{\left[{(B-(m-1)I)}_{\ell }\right]}^{-1}}{\ell !}\hfill \\ \hfill & \times \sum _{s=0}^{\infty }\frac{{(A+B+({\epsilon }^{\prime }+1-m)I)}_s{\left({\epsilon }^{\prime }\right)}_s{\left[{(A-(\ell -{\epsilon }^{\prime }-1)I)}_s\right]}^{-1}}{s!}\hfill \\ \hfill & \times \mathrm{B}(2A-(k+\ell )I+(\epsilon +{\epsilon }^{\prime })I+sI+I,2B+(k+\ell )I+(\epsilon +{\epsilon }^{\prime })I-(m+n)I+I)\hfill \\ \hfill & \times {}_3{F}_2\left(\begin{array}{c}A+B+(\epsilon +1-n)I,\epsilon ,2A-(k+\ell )I+(\epsilon +{\epsilon }^{\prime })I+sI\\ A+(\epsilon +1-k)I,2(A+B)-(n+m)I+2(\epsilon +{\epsilon }^{\prime })I+(s+2)I\end{array};1\right).\hfill \end{array}$$

In the case $\nu =n,{\nu }^{\prime }=m$, $\epsilon ={\epsilon }^{\prime }=0$, it follow that

$$\begin{array}{cc}\hfill I_{nm}& =\frac{{(-B)}_n}{\Gamma (\nu +1)}\frac{{(-B)}_m}{\Gamma (m+1)}{\Gamma }^{-1}(2(A+B)-(m+n)I+2I)\hfill \\ \hfill & \times \sum _{k=0}^n\sum _{\ell =0}^m\frac{{(-1)}^{k+\ell }{(-n)}_k{(-A)}_k{\left[{(B-(n-1)I)}_k\right]}^{-1}}{k!}\frac{{(-m)}_{\ell }{(-A)}_{\ell }{\left[{(B-(m-1)I)}_{\ell }\right]}^{-1}}{\ell !}\hfill \\ \hfill & \times \Gamma (2A-(k+\ell )I+I)\Gamma (2B+(k+\ell )I-(m+n)I+I)\hfill \end{array}$$

Using the Pochhammer identities

$$\begin{array}{cc}\hfill \Gamma (2A-(k+\ell )I+I)& ={(-1)}^{k+\ell }\Gamma (2A+1){\left[{(k-2A)}_{\ell }\right]}^{-1}{\left[{(-2A)}_k\right]}^{-1},\hfill \\ \hfill \Gamma (2B-(n+m)I+(k+\ell +1\left)I\right)& ={(2B-(n+m)I+(k+I)I)}_{\ell }{(2B-(n+m)+I)}_k\Gamma (2B-(n+m)+I)\hfill \end{array}$$

we have

$$\begin{array}{cc}\hfill I_{nm}& =\frac{{(-B)}_n}{\Gamma (\nu +1)}\frac{{(-B)}_m}{\Gamma (m+1)}{\Gamma }^{-1}(2(A+B)-(m+n)I+2I)\Gamma (2A+I)\Gamma (2B-(n+m)+I)\hfill \\ \hfill & \times \sum _{k=0}^n\sum _{\ell =0}^m\frac{{(-n)}_k{(-A)}_k{\left[{(B-(n-1)I)}_k\right]}^{-1}}{k!}\frac{{(-m)}_{\ell }{(-A)}_{\ell }{\left[{(B-(m-1)I)}_{\ell }\right]}^{-1}}{\ell !}\hfill \\ \hfill & \times {\left[{(k-2A)}_{\ell }\right]}^{-1}{\left[{(-2A)}_k\right]}^{-1}{(2B-(n+m)I+(k+I)I)}_{\ell }{(2B-(n+m)+I)}_k\hfill \end{array}$$

Finally, we obtain

$$\begin{array}{cc}\hfill I_{nm}(A,B)& =\frac{\Gamma (2A+I){(-B)}_m{(-B)}_n\Gamma (2B-(m+n)I+I)}{m!n!}{(\Gamma (2(A+B)-(m+n)I+2I))}^{-1}\hfill \\ \hfill & \times \sum _{k=0}^n\frac{{(-A)}_k{(-n)}_k{(2B-(m+n)I+I)}_k}{k!}{\left[{(-2A)}_k\right]}^{-1}{\left[{(B-nI+I)}_k\right]}^{-1}\hfill \\ \hfill & \times {}_3{F}_2(-A,-m,2B+kI-(m+n)I+I;kI-2A,B-mI+I;1).\hfill \end{array}$$

Conjecture. 

The shifted Legendre matrix functions $\left\{P_n(A,B,z)\right\}$ are orthogonal on $(0,1)$ in the following sense that:

$$\begin{array}{c}\hfill I_{nm}(A,B)=0.\end{array}$$

when $n+m$ is odd.

Particular Case

It is not known even in the scalar case if the set of polynomials $\left\{{\mathbf{P}}_{\nu }(A,B,z)\right\}$ are orthogonal in general but a trial for a particular case will be given through the following subsection. We prove the following result.

Theorem 6.

The shifted Legendre matrix functions $\left\{{\mathit{P}}_{\nu }(A,B,z)\right\}$ are orthogonal on $(0,1)$ in the following way:

$$\begin{array}{cc}\hfill \langle {\mathit{P}}_n(nI+\Phi ,nI+\Phi ;z),& {\mathit{P}}_m(mI+\Theta ,mI+\Theta ;z)\rangle \hfill \\ \hfill & ={\int }_0^1\left[{\mathit{P}}_n(nI+\Phi ,nI+\Phi ;z)\right]{\left[{\mathit{P}}_m(mI+\Theta ,mI+\Theta ;z)\right]}^{*}dz=0\hfill \end{array}$$

when $m+n$ is odd where $\Phi ,\Theta$ are $N\times N$ matrices.

Proof. 

We consider the case where $\nu =n$ for which $\epsilon =0$, then by definition

$$\begin{array}{cc}\hfill {\mathbf{P}}_n(A,B,z)=& \frac{\Gamma (A+I)\Gamma (B+I)}{n!}\hfill \\ \hfill & \times \sum _{k=0}^n{(-1)}^{n+k}\left(\genfrac{}{}{0pt}{}{n}{k}\right){\Gamma }^{-1}(A-(k-1)I){\Gamma }^{-1}(B-(n-k-1)I)z^{A-kI}{(1-z)}^{B-nI+kI}.\hfill \end{array}$$

In particular, for $A=B$, we have

$$\begin{array}{cc}\hfill {\mathbf{P}}_n(A,A;z)=& \frac{{[\Gamma (A+I)]}^2}{\Gamma (n+1)}\hfill \\ \hfill & \times \sum _{k=0}^n{(-1)}^{n+k}\left(\genfrac{}{}{0pt}{}{n}{k}\right){\Gamma }^{-1}(A-(k-1)I){\Gamma }^{-1}(A-(n-k-1)I)z^{A-kI}{(1-z)}^{A-(n-k)I}.\hfill \end{array}$$

For $A=nI+\Phi$, which is invertible, where $\Phi$ is a matrix of type $N\times N$, that is $\Phi \in {\mathbb{C}}^{N\times N}$, we get

$$\begin{array}{cc}\hfill {\mathbf{P}}_n(nI+\Phi ,nI+\Phi ;z)=& \frac{{\Gamma }^2((n+1)I+\Phi )}{\Gamma (n+1)}\sum _{k=0}^n\left(\genfrac{}{}{0pt}{}{n}{k}\right)\frac{{(-1)}^{n+k}{z}^{\Phi +(n-k)I}{(1-z)}^{(\Phi +k)I}}{\Gamma (\Phi -(-n+k-1\left)I\right)\Gamma (\Phi +(k+1\left)I\right)}\hfill \end{array}$$

(44)

Similarly, for $\mathbf{P}=mI+\Theta$ where $\Theta$ is an $N\times N$ matrix, we have

$$\begin{array}{cc}\hfill {\mathbf{P}}_m(mI+\Theta ,mI+\Theta ;z)=& \frac{{\Gamma }^2((m+1)I+\Theta )}{\Gamma \left(\right(m+1\left)I\right)}\sum _{j=0}^m\left(\genfrac{}{}{0pt}{}{m}{j}\right)\frac{{(-1)}^{m+j}{z}^{\Theta +(m-j)I}{(1-z)}^{(\Theta +k)I}}{\Gamma (\Theta -(-m+j-1\left)I\right)\Gamma (\Theta +(j+1\left)I\right)}\hfill \end{array}$$

(45)

Noting that ${\mathbf{P}}_n(A,A;z)$ is Hermitian matrix, then following [42], define the matricial inner product on ${\mathbb{C}}^{N\times N}$ for the two matrix functions $P_n(nI+\Phi ,nI+\Phi ;z)$ and ${\mathbf{P}}_m(mI+\Theta ,mI+\Theta ;z)$ defined by Equations (44) and (45) as follows.

$$\begin{array}{cc}\hfill \langle {\mathbf{P}}_n(nI+\Phi ,nI+\Phi ;z),& {\mathbf{P}}_m(mI+\Theta ,mI+\Theta ;z)\rangle \hfill \\ \hfill & ={\int }_0^1\left[{\mathbf{P}}_n(nI+\Phi ,nI+\Phi ;z)\right]{\left[{\mathbf{P}}_m(mI+\Theta ,mI+\Theta ;z)\right]}^{*}dz\hfill \end{array}$$

where ${\left[{\mathbf{P}}_m(mI+\Theta ,mI+\Theta ;z)\right]}^{*}=\overline{{\left[{\mathbf{P}}_m(mI+\Theta ,mI+\Theta ;z)\right]}^T}$ is the conjugate transpose matrix of ${\mathbf{P}}_n(nI+\Phi ,nI+\Phi ;z)$. Then,

$$\mathcal{J}={\int }_0^1{P}_n(nI+\Phi ,nI+\Phi ;z){\mathbf{P}}_m(mI+\Theta ,mI+\Theta ;z)dz,$$

where

$$\begin{array}{c}{\mathbf{P}}_n(nI+\Phi ,nI+\Phi ;z){\mathbf{P}}_m(mI+\Theta ,mI+\Theta ;z)={(-1)}^{m+n}\frac{{\Gamma }^2(n+\Phi +I){\Gamma }^2(m+\Theta +I)}{n!m!}\hfill \\ \hfill \times \sum _{k=0}^n\sum _{j=0}^m{(-1)}^{k+j}\left(\genfrac{}{}{0pt}{}{m}{j}\right)\left(\genfrac{}{}{0pt}{}{n}{k}\right)\frac{z^{\Phi +\Theta +(n+m-j-k)I}{(1-z)}^{(\Phi +\Theta +j+k)I}}{\Gamma (\Phi -(-n+k-1\left)I\right)\Gamma (\Phi +(k+1\left)I\right)\Gamma (\Theta -(-m+j-1\left)I\right)\Gamma (\Theta +(j+1\left)I\right)}\end{array}$$

Now, using the properties of the Gamma and Beta matrix functions, we conclude that the integral $\mathcal{J}$ gives

$$\begin{array}{cc}\hfill \mathcal{J}& ={(-1)}^{m+n}\frac{{\Gamma }^2(n+\Phi +I){\Gamma }^2(m+\Theta +I)}{n!m!\Gamma (2\Phi +2\Theta +(m+n+2\left)I\right)}\hfill \\ \hfill & \times \sum _{k=0}^n\sum _{j=0}^m{(-1)}^{k+j}\left(\genfrac{}{}{0pt}{}{m}{j}\right)\left(\genfrac{}{}{0pt}{}{n}{k}\right)\frac{\Gamma (\Phi +\Theta +(n+m-j-k+1\left)I\right)\Gamma (\Phi +\Theta +(j+k+1\left)I\right)}{\Gamma (\Phi -(-n+k-1\left)I\right)\Gamma (\Phi +(k+1\left)I\right)\Gamma (\Theta -(-m+j-1\left)I\right)\Gamma (\Theta +(j+1\left)I\right)},\hfill \end{array}$$

that allows us to write

$$\begin{array}{c}\hfill \mathcal{J}={(-1)}^{m+n}\frac{{\Gamma }^2(nI+\Phi +I){\Gamma }^2(mI+\Theta +I)}{n!m!\Gamma (2\Phi +2\Theta +(m+n+2\left)I\right)}\sum _{k=0}^n\sum _{j=0}^m{C}_{j,k}\end{array}$$

(46)

where

$$\begin{array}{c}\hfill C_{j,k}={(-1)}^{k+j}\left(\genfrac{}{}{0pt}{}{m}{j}\right)\left(\genfrac{}{}{0pt}{}{n}{k}\right)\frac{\Gamma (\Phi +\Theta +(n+m-j-k+1\left)I\right)\Gamma (\Phi +\Theta +(j+k+1\left)I\right)}{\Gamma (\Phi -(-n+k-1\left)I\right)\Gamma (\Phi +(k+1\left)I\right)\Gamma (\Theta -(-m+j-1\left)I\right)\Gamma (\Theta +(j+1\left)I\right)}\end{array}$$

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

$$\begin{array}{cc}\hfill C_{m-j,n-k}& ={(-1)}^{n-k+m-j}\left(\genfrac{}{}{0pt}{}{m}{m-j}\right)\left(\genfrac{}{}{0pt}{}{n}{n-k}\right)·\hfill \\ \hfill & \times \frac{\Gamma (\Phi +\Theta +(j+k+1\left)I\right)\Gamma (\Phi +\Theta +(n+m-j-k+1\left)I\right)}{\Gamma (\Phi +(k+1\left)I\right)\Gamma (\Phi -(-n+k-1\left)I\right)\Gamma (\Theta -(j+1\left)I\right)\Gamma (\Theta -(-m+j-1\left)I\right)}\hfill \\ \hfill & ={(-1)}^{m+n}{(-1)}^{k+j}\left(\genfrac{}{}{0pt}{}{m}{j}\right)\left(\genfrac{}{}{0pt}{}{n}{k}\right)·\hfill \\ \hfill & \times \frac{\Gamma (\Phi +\Theta +(n+m-j-k+1\left)I\right)\Gamma (\Phi +\Theta +(j+k+1\left)I\right)}{\Gamma (\Phi -(-n+k-1\left)I\right)\Gamma (\Phi +(k+1\left)I\right)\Gamma (\Theta -(-m+j-1\left)I\right)\Gamma (\Theta +(j+1\left)I\right)}\hfill \\ \hfill & ={(-1)}^{m+n}{C}_{j,k},j=0,1,\cdots ,m;k=0,1,\dots ,n\hfill \end{array}$$

Note that, for $m+n$ odd, we get, for instance, if $m+n=1$, ${\sum }_{k=0}{\sum }_{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 ${\sum }_{k=0}^n{\sum }_{j=0}^m{C}_{j,k}$ has its opposite signed companion, therefore ${\sum }_{k=0}^n{\sum }_{j=0}^m{C}_{j,k}=0$ and the orthogonality property is verified. □

Remark 6.

Note that, for $m+n$ is odd, which means $m\ne n$ and for $m=n=0$, the relation in Equation (46) gives

$$\mathcal{J}=\beta (2\Theta +I,2\Theta +I).$$

(47)

8. 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).

8.1. Technical Scheme

The explicit analytical form of the shifted Legendre matrix polynomials ${\mathbf{P}}_n(A,B;z)$ of degree n may be written as:

$$\begin{array}{c}\hfill {\mathbf{P}}_n(A,B;z)=\sum _{k,s=0}^n{Q}_{k,s}{z}^{A-(k-s)I}\end{array}$$

(48)

where

$$\begin{array}{c}\hfill Q_{k,s}=\left(\genfrac{}{}{0pt}{}{n}{k}\right){(-1)}^{n+k+s}\frac{\Gamma (A+I)\Gamma (B+I)\Gamma (A+(s-k\left)I\right)}{s!\Gamma (n+1)\Gamma (B-(n+k-s-1\left)I\right)\Gamma (A-(k-1\left)I\right)\Gamma (A-kI)\Gamma (A-(k-s\left)I\right)}\end{array}$$

such that $\parallel A\parallel <\parallel (k-s)I\parallel$. Theorem 6 leads to the orthogonality condition in the form:

$$\mathcal{J}={\int }_0^1{\mathbf{P}}_i(A,B;z){\mathbf{P}}_j(A,B;z)dz=\left\{\begin{array}{cc}0,\hfill & \mathrm{for}i+j\mathrm{is}\mathrm{odd},(i\ne j)\hfill \\ \beta (2\Theta +I,2\Theta +I),\hfill & \mathrm{for}i=j=0.\hfill \end{array}\right.$$

Thus, any square integrable function $y\left(z\right)$ defined on $(0,1)$ may be expressed in terms of shifted Legendre matrix polynomials as

$$y\left(z\right)=\sum _{j=0}^{\infty }{c}_j{\mathbf{P}}_j(A,B;z)$$

from which the coefficients $c_j$ are given by

$$\begin{array}{c}\hfill c_j=\beta (2\Theta +I,2\Theta +I){\int }_0^{1}y\left(z\right){\mathbf{P}}_j(A,B;z)dz\end{array}$$

If we approximate $y\left(z\right)$ by truncating the first $(m+1)$ terms, we have

$$y_m\left(z\right)\simeq \sum _{j=0}^m{c}_j{\mathbf{P}}_j(A,B,z),$$

which alternatively may be written in the matrix form as

$$y\left(z\right)=C^T{\Delta }_m(A,B;z),$$

where the shifted Legendre coefficient $(m+1)\times 1$ vector C and the shifted Legendre vector ${\Delta }_m(A,B,z)$ are given by

$$\begin{array}{c}\hfill C^T=[c_0,c_1,\cdots ,c_m]\end{array}$$

(49)

and

$$\begin{array}{c}\hfill {\Delta }_m(A,B;z)={\left[{\mathbf{P}}_0(A,B;z),{\mathbf{P}}_1(A,B;z),\cdots ,{\mathbf{P}}_m(A,B;z)\right]}^T,\end{array}$$

(50)

respectively.

8.2. 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 ${\Delta }_m(A,B;z)$ can be expressed by

$$\begin{array}{c}\hfill \frac{d}{dz}{\Delta }_m(A,B;z)={\mathcal{D}}^{\left(1\right)}{\Delta }_m(A,B;z)\end{array}$$

(51)

where ${\mathcal{D}}^{\left(1\right)}$ is the $(m+1)\times (m+1)$ operational matrix of derivatives given by

$${\mathcal{D}}^{\left(1\right)}=\left(d_{ij}\right)=\left\{\begin{array}{cc}2\beta (2\Theta +I,2\Theta +I),\hfill & \mathrm{for}j=i-k,k=1,3,\cdots ,m\mathrm{if}m\mathrm{odd}\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

8.3. Legendre Operational Matrix in Fractional Setting

It is clear using Equation (51) that

$$\begin{array}{c}\hfill \frac{d^n{\Delta }_m(A,B;z)}{d{z}^n}={\left({\mathcal{D}}^{\left(1\right)}\right)}^n{\Delta }_n(A,B;z)\end{array}$$

(52)

where $n\in \mathbb{N}$ and the subscript, in ${\mathcal{D}}^{\left(1\right)}$, denotes matrix powers. Thus,

$$\begin{array}{c}\hfill {\mathcal{D}}^{\left(n\right)}={\left({\mathcal{D}}^{\left(1\right)}\right)}^n,n=1,2,\cdots \end{array}$$

(53)

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:

$$\begin{array}{c}\hfill {\mathcal{D}}^{\alpha }y\left(x\right)=\sum _{s=1}^k{A}_s{D}^{{\beta }_s}y\left(x\right)+A_{k+1}y\left(x\right)+A_{k+2}y\left(x\right)\end{array}$$

(54)

with initial conditions

$$\begin{array}{c}\hfill y^{\left(i\right)}\left(0\right)=d_i,i=0,1,\cdots ,n\end{array}$$

(55)

where $A_s$ for $s=1,2,\cdots ,k+2$ are real coefficients and $n<\alpha \le n+1,$$0<{\beta }_1<{\beta }_2<\cdots <{\beta }_k<\alpha$ and ${\mathcal{D}}^{\alpha }$ refers to the Caputo fractional derivative of order $\alpha$. To solve the problem in Equations (54) and (55), it is required to approximate $y\left(x\right)$ and $g\left(x\right)$ by the shifted Legendre matrix polynomials as follows.

$$\begin{array}{c}\hfill y\left(z\right)\simeq \sum _{i=0}^m{c}_i{\mathbf{P}}_i(A,B;z)=C^T{\Delta }_m(A,B;z)\end{array}$$

(56)

$$\begin{array}{c}\hfill g\left(z\right)\simeq \sum _{i=0}^m{g}_i{\mathbf{P}}_i(A,B;z)=G^T{\Delta }_m(A,B;z)\end{array}$$

(57)

where the vector $G={\left[g_0,\cdots ,g_m\right]}^T$ is known but the vector $C={\left[c_0,\cdots ,c_m\right]}^T$ is unknown. Following similar procedure as in the scalar case [7] leads to the following systems.

$$\begin{array}{cc}\hfill D^{\alpha }y\left(x\right)& \simeq C^T{D}^{\alpha }{\Delta }_m(A,B;z)\simeq C^T{D}^{\left(\alpha \right)}{\Delta }_m(A,B;z)\end{array}$$

(58)

$$\begin{array}{cc}\hfill D^{{\beta }_j}y\left(x\right)& \simeq C^T{D}^{{\beta }_j}{\Delta }_m(A,B;z)\simeq C^T{D}^{\left({\beta }_j\right)}{\Delta }_m(A,B;z),j=1,\cdots ,k\end{array}$$

(59)

Using Equations (56)–(59), the residual $R_m\left(z\right)$ for Equation (54) can be written as

$$\begin{array}{c}\hfill R_m(A,B;z)\simeq \left[C^T{\mathcal{D}}^{\left(\alpha \right)}-C^T\sum _{j=1}^k{A}_j{\mathcal{D}}^{\left({\beta }_j\right)}-A_{k+1}{C}^T-a_{k+2}{G}^T\right]{\Delta }_m(A,B;z)\end{array}$$

(60)

Owing to Tau method [46], it follows that $m-n$ linear equations can be generated by using

$$\begin{array}{c}\hfill \langle R_m(A,B;z),{\mathbf{P}}_j(A,B;z)\rangle ={\int }_0^1{R}_m(A,B;z){\mathbf{P}}_j(A,B;z)dz=0,j=0,1,\cdots ,m-n-1\end{array}$$

(61)

Using Equations (53) and (56) in Equation (55), we obtain

$$\begin{array}{cc}\hfill y\left(0\right)& =C^T{\Delta }_m\left(0\right)=d_0,\hfill \\ \hfill y{}^{\left(1\right)}\left(0\right)& =C^T{\mathcal{D}}^{\left(1\right)}{\Delta }_m\left(0\right)=d_1\hfill \\ \hfill & ⋮\hfill \\ \hfill y^{\left(n\right)}\left(0\right)& =C^T{\mathcal{D}}^{\left(n\right)}{\Delta }_m\left(0\right)=d_n.\hfill \end{array}$$

(62)

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\left(z\right)$ 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=\nu =n\in {\mathbb{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].

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill D^{2}y\left(z\right)& +D^{3/2}y\left(z\right)+y\left(z\right)=1+z\hfill \\ \hfill & y\left(0\right)=1,y^{\prime }\left(0\right)=1\hfill \end{array}\end{array}$$

(63)

This problem has the exact solution in the form $y\left(z\right)=1+z$. By applying the technical scheme described above with $m=2,$$A=B=\nu =n\in {\mathbb{N}}_0$ and $N=1$, one can approximate solution as

$$y\left(z\right)=c_0{\mathit{P}}_0\left(z\right)+c_1{\mathit{P}}_1\left(z\right)+c_2{\mathit{P}}_2\left(z\right)=C^T{\Delta }_m\left(z\right)$$

Acting by fractional derivative, we get

$${\mathcal{D}}^{\left(1\right)}=\left[\begin{array}{ccc}0& 0& 0\\ 2& 0& 0\\ 0& 6& 0\end{array}\right],{\mathcal{D}}^{\left(2\right)}=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 12& 0& 0\end{array}\right],{\mathcal{D}}^{(3/2)}=\frac{16}{\sqrt{\pi }}\left[\begin{array}{ccc}0& 0& 0\\ 2& 0& 0\\ 1& \frac{3}{5}& \frac{-1}{7}\end{array}\right],G=\left[\begin{array}{c}\frac{3}{2}\\ \frac{1}{2}\\ 0\end{array}\right]$$

From this, one can deduce

$$c_0+\left(12+\frac{16}{\sqrt{\pi }}\right)c_2-\frac{3}{2}=0$$

Applying Equation (62), we get

$$\begin{array}{cc}\hfill c_0-c_1+c_2-1& =0\hfill \\ \hfill 2c_1-6c_2-1& =0\hfill \end{array}$$

Therefore, $c_0=\frac{3}{2}$, $c_1=\frac{1}{2},$$c_2=0$. Thus, we can write

$$y\left(z\right)=\left[\begin{array}{ccc}\frac{3}{2}& \frac{1}{2}& 0\end{array}\right]\left[\begin{array}{c}1\\ 2z-1\\ 6z^2-6z+1\end{array}\right]=z+1$$

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

$$\begin{array}{cc}\hfill \frac{{\partial }^{{\alpha }_1}}{\partial t^{{\alpha }_1}}u(x,t)& =\frac{-{\partial }^3}{\partial x^3}u(x,t)-6u(x,t)\frac{\partial }{\partial x}u(x,t)+3v(x,t)\frac{\partial }{\partial x}v(x,t)\hfill \\ \hfill \frac{{\partial }^{{\alpha }_2}}{\partial t^{{\alpha }_2}}u(x,t)& =\frac{-{\partial }^3}{\partial x^3}u(x,t)-3u(x,t)\frac{\partial }{\partial x}v(x,t)\hfill \end{array}$$

with the initial condition $u(x,0)=\frac{4c^2{e}^{cx}}{{(1+e^{cx})}^2},$$v(x,0)=\frac{4c^2{e}^{cx}}{{(1+e^{cx})}^2}$. The exact solution of this problem, for ${\alpha }_1={\alpha }_2=1$, is

$$u(x,t)=v(x,t)=\frac{4c^2{e}^{c(x-c^{2}t)}}{{(1+e^{c(x-c^{2}t)})}^2}.$$

9. Conclusions

Classically, the Rodrigues differential formulas are often used as definitions of the orthogonal polynomials and their generalizations in special functions. Many of the essential mathematical properties of the orthogonal polynomials can be derived using these definitions.

Notably, the fractional-order analogs of the Rodrigues’ formulas were recently used to either define new classes of special functions or to study and classify many of the special functions in fractional calculus (see, e.g., [12,16]). The study developed in Section 3, Section 4, Section 5, Section 6 and Section 7 provides several mathematical properties of the shifted Legendre-type matrix polynomials under the assumption in Equation (5) in terms of the fractional-order Rodrigues’ formula.

It should be noted that the integral representation for the shifted Legendre matrix function has been used in studying the delayed, scaled consensus problems in the work of Shang [41].

By adopting a matrix setting with the notion of the Hermitian matrices, we introduce the notion of a matricial inner product on ${\mathbb{C}}^{N\times N}$. This approach enabled us to establish an orthogonality property in a matrix case for some particular cases. Although a general orthogonality condition in the matrix setting was not entirely resolved in the present study, we hope the present work will encourage further investigations in that direction.

To the best of our knowledge, the introduced fractional setting and matrix framework for shifted Legendre-Type matrix polynomials is newly presented. Based on the earlier work related to the shifted Legendre polynomials, we may suggest some possible applications to be studied using the approach introduced in the present work.

• Referring to the work given in [32], where the authors used the shifted Legendre collocation method for non-linear multi-order fractional differential equations, a similar problem may be treated through the matrix framework. More precisely, we may suggest the analog of FDE stated as follows

  $$D^{\nu }f\left(z\right)+A{D}^{{\nu }_1}f\left(z\right)+B{D}^{{\nu }_2}f\left(z\right)+Cf\left(z\right)=g\left(z\right),z\in \mathbb{C},A,B,C\in {\mathbb{C}}^{N\times N}$$

  subject to some boundary conditions.
• Another application is possible based on [34], where the authors applied shifted Legendre polynomials to the state equations of a linear system. It is possible to investigate similar problems utilizing fractional shifted Legendre polynomials in the sense of matrix framework, as discussed in the present work. We may impose the following question: “Can the state and forcing functions be expressed in terms of a series of fractional shifted Legendre polynomials with expansion coefficients?”
• It is also an interesting problem to investigate the solution of a linear dynamic system with initial or boundary value conditions by the proposed fractional shifted Legendre-type matrix polynomials as a generalization of the work given by Wang et al. [35].
• As the operational matrix of shifted Legendre polynomials in the scalar case has been employed to approach the identification problem (see [36]), a similar problem may be investigated in the sense of our proposed approach.

The results given in Section 7 are obtained for a particular case. A general orthogonality relation in the matrix setting involving fractional sense is still an open problem. We hope that the present study will encourage further investigation to solve this problem. Obtaining such tools, would shed light for the numerical analysis field with applications.