# Direct Computation of Operational Matrices for Polynomial Bases

Authors: Osvaldo Guimarães; José Roberto C. Piqueira; Marcio Lobo Netto
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

Several numerical methods for boundary value problems use integral and differential operational matrices, expressed in polynomial bases in a Hilbert space of functions. This work presents a sequence of matrix operations allowing a direct computation of operational matrices for polynomial bases, orthogonal or not, starting with any previously known reference
matrix. Furthermore, it shows how to obtain the reference matrix for a chosen polynomial base. The results presented here can be applied not only for integration and differentiation, but also for any linear operation.


---

## Body

## 1. Introduction

One of the main characteristics of the use of polynomial bases is to reduce the solving process of differential or integral equations to systems of algebraic equations, expressing the solutionf(x) by truncated series approximations, up to order n [1–4], such that(1.1)f(x)=fN(x)=∑i=1NciPi(x),x∈[a,b]⊂ℜ.The choice of the polynomial basis is normally one of the orthogonal bases belonging to the Hilbert space of functions, in order to ensure that the expansion of the series to a higher order does not affect the coefficients previously calculated, being applicable to classical methods, as the Runge-Kutta, for instance [5]. However, it is also possible to use nonorthogonal bases, as done in [6, 7], where Pi(x)=xi,  i=1,2,…,n.ConsideringC the line vector that contains the coefficients ci and B the column vector that contains the base polynomials Pi, expression (1.1) can be written as considering: fN(x)=[c0c1,…,cn][P0P1,…,Pn]T=CB [8].The central idea when working with operational matrices is to write the integral or differential of the elements of the basis as a linear combination of the same base elements, transforming the integral and differential operations offN(x) into matrix operations in a Hilbert space [8].So, definingMI as the operational integration matrix and MD as the operational differential matrix, it is possible to obtain the line vector V containing the coefficients of the series that represent the integrated function or the differentiated function by VIT=MICT and VDT=MDCT.Consequently,(1.2)∫axfN(x)dx=VIB,d[fN(x)]dx=fN(1)(x)=VDB.Recently, Doha and Bhrawy [1] presented a method to obtain the operational matrices of integration considering the Jacobi polynomials.Here, a simpler and more direct way to get the operational matrix, by using Theorem2.1 from the next section is presented. Additionally, a way to extend it to any polynomial basis, by using Theorem 3.1, presented in Section 3, is also developed.In spite of the fact that those theorems are applied to integration and differentiation operations, the result is valid to any linear operation, as shown ahead.


## 2. Obtaining the Operational Matrix

Theorem 2.1.
ConsideringZ=Zn+1,n+1 a square matrix describing the resulting coefficients of a linear operation α in the generic basis BG=[b0b1,…,bn]T as a function of the same basis and fN(x) a series the coefficients vector of which is C=[c0c1,…,cn], and V=[v0v1,…,vn] the line vector of the resulting coefficients of the linear operation α applied to the series, then:
(2.1)VT=ZTCT.Proof.
ConsideringfN(x)=∑i=0nciBGi(x), since α is a linear operation: α[fN(x)]=∑i=0nciαBGi(x).
On the other hand, considering the finite matrix spaceMij, i∈{0,1,…,n} and j∈{0,1,…,n}, α is a scalar, depending on x.
Considering thatαBG=ZBG and VBG=C(ZBG), observing the matrix and vector dimensions, this last result is a scalar and, consequently BGTVT=(ZBG)TCT, implying that:(2.2)BGT(VT-ZTCT)=0.
AsBGT is a generic basis, (2.2) implies (2.1) directly.So, in order to build matrices representing actions of linear operations, as derivative and integration, the main task is to determine matrixZ and transpose it.


### 2.1. Example: Integration Matrix of the Legendre Polynomials

Considering that the polynomial basis is used to describe a function to be composed of Legendre polynomials, in the interval[a,b]=[-1,1], one can observe that, for Legendre polynomials [9](2.3)∫-1xPi(x)dx=Pi+1(x)-Pi-1(x)2i+1,i=1,2,…,n,∫-1xP0(x)dx=P1(x)+P0.Defining(2.4)P=[P0P1⋯0·Pn+1],P*=[P0P1⋯Pn+1],
one can write ∫-1xPdx=ZIP*, with the integral acting over to the elements of the vector. The null coefficient on the (n+1)-order term assures equal dimensions for the vector to be integrated and the vector that results from the process.Following (2.3), ZI is a square matrix Zn+2, n+2 and:(i)
Z1,1=Z1,2=1;(ii)
Zk,k+1=-1/(2k-1),  Zk,k-1=1/(2k-1),  k=2,3,…,n+1;(iii)
Zn+2,n+1=-1/(2(n+2)-1)=-1/(2n+3);(iv)
there is noZn+2,n+3 term, because Pn+1 is not integrated.Since the operational matrixMI is the transpose of ZI, when writing a code to implement a computational algorithm, it is required to exchange the indices in the expressions above, obtaining MI. directly.


## 2.1. Example: Integration Matrix of the Legendre Polynomials

Considering that the polynomial basis is used to describe a function to be composed of Legendre polynomials, in the interval[a,b]=[-1,1], one can observe that, for Legendre polynomials [9](2.3)∫-1xPi(x)dx=Pi+1(x)-Pi-1(x)2i+1,i=1,2,…,n,∫-1xP0(x)dx=P1(x)+P0.Defining(2.4)P=[P0P1⋯0·Pn+1],P*=[P0P1⋯Pn+1],
one can write ∫-1xPdx=ZIP*, with the integral acting over to the elements of the vector. The null coefficient on the (n+1)-order term assures equal dimensions for the vector to be integrated and the vector that results from the process.Following (2.3), ZI is a square matrix Zn+2, n+2 and:(i)
Z1,1=Z1,2=1;(ii)
Zk,k+1=-1/(2k-1),  Zk,k-1=1/(2k-1),  k=2,3,…,n+1;(iii)
Zn+2,n+1=-1/(2(n+2)-1)=-1/(2n+3);(iv)
there is noZn+2,n+3 term, because Pn+1 is not integrated.Since the operational matrixMI is the transpose of ZI, when writing a code to implement a computational algorithm, it is required to exchange the indices in the expressions above, obtaining MI. directly.


## 3. Obtaining the Operational Matrix for Any Polynomial Basis:  The Sandwich Matrix

From an operational matrix expressed in a generic reference basisX, it is possible to obtain the corresponding operational matrix in another basis G, also generic, by using a sequence of simple matrix operations, as Theorem 3.1 states.Theorem 3.1 (“Sandwich Matrix”(Ω)).
Given a generic polynomial basisBG in the interval [a,b], the matrix ZG of the operational matrix theorem is obtained by GΩG-1, where the “sandwich matrix” is Ω=X-1ZXX, with Gn+1,n=1 and Xn+1,n+1 being the matrices that describe the generic polynomials (G) and the reference ones (X) in the canonic basis, respectively.Proof.
Considering thatα is a linear operation: αBG=α[GBC]=Gα[BC], where BC is the canonic base BC=[1xx2,…,xn]T, since the canonic base can be written as a function of the reference base as BC=X-1BX, it can be concluded that αBG=α[GBC]=G(X-1α[BX]).
Asα[BX]=ZXBX, it implies that α[BX]=G(X-1)ZXBX. Now, it is necessary to transform BX into the generic basis BG, in order to express the result as a function of the generic input base.
From the defined matricesBX=XBC and BC=G-1BG, the expression given above is written as: αBG=G(X-1ZXX)G-1BG=GΩG-1BG and G must be nonsingular.
Finally, sinceαBG=ZGBG, the generic matrix ZG is obtained by: ZG=GΩG-1, where Ω=(X-1)ZXX.


### 3.1. Comments

(i) No orthogonality condition has been imposed during the proof, thus, the result is valid for all polynomial bases, orthogonal or not.(ii) The operational matrix of a generic polynomial basis is given by:MG=(G-1)TΩTGT. Indeed, for any linear operation: ZG=G(X-1ZXX)G-1 and, since MG=ZGT, it can be concluded that: MG=(G-1)TΩTGT.(iii) SinceX is arbitrary, the canonic basis can be used, and thus, the matrix describing the elements of X as a function of the canonic base is the identity. So, Ω=I-1ZCI=ZC.(iv) Taking the previous comment into account, the matrixZG of the generic basis can be obtained by: ZG=GZCG-1, that is, considering the resemblance definition [9], the matrix of the generic basis ZG is similar to or resembling matrix ZC. To summarizing, the generating matrices of the operational matrices are similarity classes related to each linear operation in the interval [a,b].(v) Considering the uniqueness of the result of the linear operation of differentiation on continuous functions, the differentiation matrixΩ is invariant for all polynomial bases.(vi) The elements of the last column of the integration matrix are arbitrary, since they multiply a null element. Therefore, if instead of the Legendre basis another one is used to build the “sandwich matrix"Ω, these elements may be different, but the result of the integration remains the same. So, the uniqueness of the “sandwich matrix" for the integration is ensured but for the last line, the impact of which will be on the last column of the integration matrix.


### 3.2. Example: “Sandwich Matrix" for Integration

According to the last comment about Theorem3.1, any polynomial basis can be used to build the “sandwich matrix" Ω, because they are all similar. For the sake of simplicity, the canonic base will be chosen.Considering the interval[-1,1], the matrix Z is built, in order to represent the integrals of the basis as a function of the basis itself. By defining(3.1)X=[x0x1⋯0·xn+1],X*=[x0x1⋯xn+1],
one can write ∫-1xXdx=ZCIX*, with the integral acting over to the elements of the vector, that is, ∫-1xξidξ=(1/(i+1))xi+1-(-1)i+1/(i+1),  i=0,1,…,n.Consequently,∫-1xξidξ=(1/(i+1))Pi+1-(-P0)i+1/(i+1)=(1/(i+1))Pi+1-(-1)iP0/(i+1)=(1/(i+1))Pi+1-(-1)i/(i+1),  i=0,1,…,n.The nonnull elementsZCI=Zn+2,n+2 are, therefore: Zk,k+1=1/k and Zk,1=(-1)k-1/k, with k=1,2,…,n+1. As apparent in comments presented in the former subsection, for each linear operation, in the canonic basis, ZC=Ω.Considering, for instance,n=4, the “sandwich matrix” for integration (ΩI) is:(3.2)ΩI=[110000-1/201200013001300-1400014015100015000000].For the interval[0,b], the elements Zk,1, with k=1,2,…,n+1 will be equal to zero, that is,(3.3)ΩI=[0100000012000000130000001400100015000000].Considering the process using the shifted Legendre base up to order 4, known as shifted-Legendre,Ω is: ΩLI=L-1ZLIL, where ZLIT=MLI is the integration operational matrix in the Legendre basis. For polynomials up to order 5, the matrix describing the Legendre polynomials as a function of the canonic basis is:(3.4)L=[100000010000-120320000-3205200380-1540358001580-3540638].CalculatingΩLI=L-1ZLIL(3.5)ΩLI=[110000-1201200013001300-1400014015100015-1136930-56605220].As mentioned, this matrix differs from the one previously obtained with the canonic basis just by the last line elements. It does not modify the integration operation, since these elements are being multiplied by a null coefficient.The wavelets, orthogonal Jacobi polynomials shifted to the interval (u∈ℜ:0≤u≤1), were used in the Galerkin processes [10], with the appropriate domain transformation. Matrix Ω can also be applied in these cases, either transforming the equation or obtaining the matrix that describes the base polynomials in the chosen interval related to the canonic base.


### 3.3. Example: “Sandwich Matrix" for Differentiation

Starting with the canonic base, matrixZD can be built, describing the derivative elements of the basis as a function of the basis itself. Therefore, rewriting the vector X* defined in the former subsection: dX*/dt=ZX*, with the derivative acting over the components of vector X*,  d[x0]/dx=0,…,d[xn]/dx=nxn+1, for any possible considered domain.Thus, the nonnull terms of the matrixZD=Zn+1,n+1 are Zi,i-1=i-1, i=2,…,n and, according to the last comment from Section 3.1, ZCD=ΩD.Consideringn=4, the differentiation matrix is given by(3.6)ΩD=[0000010000020000030000040].In order to obtain this matrix from the Legendre basis, the process is analogous to the previous one. ApplyingΩD=L-1ZLDL, where ZLDT is the described differentiation operational matrix. As highlighted in the third from Theorem 3.1, presented in Section 3.1, the obtained matrix Ω is identical.


### 3.4. Example: Chebyshev Operational Matrices for Integration and Differentiation

Some works present the Galerkin method supported by ChebyshevI expansions [11–13], when solving differential equations. In order to help in this task, integration and differentiation Chebyshev I matrices will be obtained from the operational matrices on the Legendre basis, even though it would be easier to conduct this process using the canonic base.MatrixT describing the Chebyshev I polynomials, as a function of the canonic basis, up to order 4 is:(3.7)T=[1000001000-102000-304010-808].By using Theorem3.1, the Legendre integration matrix (MLI) is known, and Chebyshev integration matrix can be obtained, by using the second comment from Section 3.1.Consider, for instance,f(x)=-2-10x+6x2+16x3, x∈[-1,1] to be integrated in the interval [-1,x]. Written as a Chebyshev I series: f(x)=[1234][T0T1T2T3]T and MTI=(T-1)TΩITTT, observing that Ω is basis independent.Calculating the several matrices is the case of the example(3.8)T=[1000001000-102000-304010-808],ΩIT=[1-1213-14151000-12101200000130290001/40],MTI=[0-14-1318-1510-120-1210140-14000160-290001/80].
Performing V=MTI[1234]T, with the coefficients of the series representing the integral of this function, that is, F(x)=∫-1xf(x)dx=VBT, where V=[0-1/2-1/21/21/2] and, consequently, in the canonic form: F(x)=1-2x-5x2+2x3+4x4.To obtain the Chebyshev differential matrix, the procedure is analogous, giving:(3.9)ΩDT=[0100000200000300000410000],MTD=[0103000408000600000800000].As expected,MTDVT=[12340]T, with the nonnull terms of the invariant ΩD(n+1,n+1)T given by ΩD(i,i+1)T=i,  i=1,…,n.


## 3.1. Comments

(i) No orthogonality condition has been imposed during the proof, thus, the result is valid for all polynomial bases, orthogonal or not.(ii) The operational matrix of a generic polynomial basis is given by:MG=(G-1)TΩTGT. Indeed, for any linear operation: ZG=G(X-1ZXX)G-1 and, since MG=ZGT, it can be concluded that: MG=(G-1)TΩTGT.(iii) SinceX is arbitrary, the canonic basis can be used, and thus, the matrix describing the elements of X as a function of the canonic base is the identity. So, Ω=I-1ZCI=ZC.(iv) Taking the previous comment into account, the matrixZG of the generic basis can be obtained by: ZG=GZCG-1, that is, considering the resemblance definition [9], the matrix of the generic basis ZG is similar to or resembling matrix ZC. To summarizing, the generating matrices of the operational matrices are similarity classes related to each linear operation in the interval [a,b].(v) Considering the uniqueness of the result of the linear operation of differentiation on continuous functions, the differentiation matrixΩ is invariant for all polynomial bases.(vi) The elements of the last column of the integration matrix are arbitrary, since they multiply a null element. Therefore, if instead of the Legendre basis another one is used to build the “sandwich matrix"Ω, these elements may be different, but the result of the integration remains the same. So, the uniqueness of the “sandwich matrix" for the integration is ensured but for the last line, the impact of which will be on the last column of the integration matrix.


## 3.2. Example: “Sandwich Matrix" for Integration

According to the last comment about Theorem3.1, any polynomial basis can be used to build the “sandwich matrix" Ω, because they are all similar. For the sake of simplicity, the canonic base will be chosen.Considering the interval[-1,1], the matrix Z is built, in order to represent the integrals of the basis as a function of the basis itself. By defining(3.1)X=[x0x1⋯0·xn+1],X*=[x0x1⋯xn+1],
one can write ∫-1xXdx=ZCIX*, with the integral acting over to the elements of the vector, that is, ∫-1xξidξ=(1/(i+1))xi+1-(-1)i+1/(i+1),  i=0,1,…,n.Consequently,∫-1xξidξ=(1/(i+1))Pi+1-(-P0)i+1/(i+1)=(1/(i+1))Pi+1-(-1)iP0/(i+1)=(1/(i+1))Pi+1-(-1)i/(i+1),  i=0,1,…,n.The nonnull elementsZCI=Zn+2,n+2 are, therefore: Zk,k+1=1/k and Zk,1=(-1)k-1/k, with k=1,2,…,n+1. As apparent in comments presented in the former subsection, for each linear operation, in the canonic basis, ZC=Ω.Considering, for instance,n=4, the “sandwich matrix” for integration (ΩI) is:(3.2)ΩI=[110000-1/201200013001300-1400014015100015000000].For the interval[0,b], the elements Zk,1, with k=1,2,…,n+1 will be equal to zero, that is,(3.3)ΩI=[0100000012000000130000001400100015000000].Considering the process using the shifted Legendre base up to order 4, known as shifted-Legendre,Ω is: ΩLI=L-1ZLIL, where ZLIT=MLI is the integration operational matrix in the Legendre basis. For polynomials up to order 5, the matrix describing the Legendre polynomials as a function of the canonic basis is:(3.4)L=[100000010000-120320000-3205200380-1540358001580-3540638].CalculatingΩLI=L-1ZLIL(3.5)ΩLI=[110000-1201200013001300-1400014015100015-1136930-56605220].As mentioned, this matrix differs from the one previously obtained with the canonic basis just by the last line elements. It does not modify the integration operation, since these elements are being multiplied by a null coefficient.The wavelets, orthogonal Jacobi polynomials shifted to the interval (u∈ℜ:0≤u≤1), were used in the Galerkin processes [10], with the appropriate domain transformation. Matrix Ω can also be applied in these cases, either transforming the equation or obtaining the matrix that describes the base polynomials in the chosen interval related to the canonic base.


## 3.3. Example: “Sandwich Matrix" for Differentiation

Starting with the canonic base, matrixZD can be built, describing the derivative elements of the basis as a function of the basis itself. Therefore, rewriting the vector X* defined in the former subsection: dX*/dt=ZX*, with the derivative acting over the components of vector X*,  d[x0]/dx=0,…,d[xn]/dx=nxn+1, for any possible considered domain.Thus, the nonnull terms of the matrixZD=Zn+1,n+1 are Zi,i-1=i-1, i=2,…,n and, according to the last comment from Section 3.1, ZCD=ΩD.Consideringn=4, the differentiation matrix is given by(3.6)ΩD=[0000010000020000030000040].In order to obtain this matrix from the Legendre basis, the process is analogous to the previous one. ApplyingΩD=L-1ZLDL, where ZLDT is the described differentiation operational matrix. As highlighted in the third from Theorem 3.1, presented in Section 3.1, the obtained matrix Ω is identical.


## 3.4. Example: Chebyshev Operational Matrices for Integration and Differentiation

Some works present the Galerkin method supported by ChebyshevI expansions [11–13], when solving differential equations. In order to help in this task, integration and differentiation Chebyshev I matrices will be obtained from the operational matrices on the Legendre basis, even though it would be easier to conduct this process using the canonic base.MatrixT describing the Chebyshev I polynomials, as a function of the canonic basis, up to order 4 is:(3.7)T=[1000001000-102000-304010-808].By using Theorem3.1, the Legendre integration matrix (MLI) is known, and Chebyshev integration matrix can be obtained, by using the second comment from Section 3.1.Consider, for instance,f(x)=-2-10x+6x2+16x3, x∈[-1,1] to be integrated in the interval [-1,x]. Written as a Chebyshev I series: f(x)=[1234][T0T1T2T3]T and MTI=(T-1)TΩITTT, observing that Ω is basis independent.Calculating the several matrices is the case of the example(3.8)T=[1000001000-102000-304010-808],ΩIT=[1-1213-14151000-12101200000130290001/40],MTI=[0-14-1318-1510-120-1210140-14000160-290001/80].
Performing V=MTI[1234]T, with the coefficients of the series representing the integral of this function, that is, F(x)=∫-1xf(x)dx=VBT, where V=[0-1/2-1/21/21/2] and, consequently, in the canonic form: F(x)=1-2x-5x2+2x3+4x4.To obtain the Chebyshev differential matrix, the procedure is analogous, giving:(3.9)ΩDT=[0100000200000300000410000],MTD=[0103000408000600000800000].As expected,MTDVT=[12340]T, with the nonnull terms of the invariant ΩD(n+1,n+1)T given by ΩD(i,i+1)T=i,  i=1,…,n.


## 4. Solving a Boundary Value Problem

In this section, an application of the method presented to build operational matrices is shown, considering the boundary value problem related to convection-diffusion equation [14], given by:(4.1)-u′′(v)+bu′(v)=0,u:[0,1]→[0,1],
with u(0)=0 and u(1)=1.Firstly, in order to have a Jacobi interval domain, we change variables,v=(x+1)/2, with x∈[-1,1] obtaining the transformed equation:(4.2)-4u′′(x)+b2u′(x)=0;u:[-1,1]→[0,1],
with u(-1)=0 and u(1)=1.Ifu*(x)=∑i=0nciPi is the series that approximates u(x) and MD the operational differentiation matrix, one can write: (-4MD2CT+2bMDCT)T·P(xk)=0, defining:(4.3)P(xk)=[P0(xk)P1xk⋯Pn(xk)],
with C=[c0c1⋯cn] and k=1,2,…,n-1.This matrix equation is applied ton-2 domain points generating a linear algebraic equations system with n-2 equations and n unknown variables. The two missing equations are obtained from the boundary conditions: x0=-1 and xn=1. To avoid the Runge phenomenon [15], xk are chosen as nodes of the polynomial basis.Figure1 shows the exact solution and the obtained by using Legendre approximation and considering b=1 and n=5. The two solutions are too closed that, in Figure 2, the error is shown for comparison. Figure 3 shows the exact solution and the obtained by using Chebyshev approximation and considering b=20 and n=25 The two solutions are too closed that, in Figure 4, the error is shown for comparison.Figure 1
Solution of convection-diffusion equation with Legendre approximation (b=1,n=5).Figure 2
Error in the solution of convection-diffusion equation with Legendre approximation (b=1,n=5).Figure 3
Solution of convection-diffusion equation with Chebyshev approximation (b=20; n=25).Figure 4
Error in the solution of convection-diffusion equation with Chebyshev approximation (b=20,n=25).


## 5. Conclusion

All operational matrices applied to polynomial bases in linear operations may be obtained directly from a central matrix(Ω) placed between the matrix product involving the matrix describing the chosen base from the canonical base and its inverse.Considering the available computational facilities, this method may turn the calculation of these matrices easier and quicker, on different bases, and various applications, as the Galerkin process, for instance. Furthermore, “sandwich matrix" allows for directly obtaining the recurrence relations for the derivative and integral of an element of any polynomial basis as a function of other basis elements.



---

*Source: 139198-2010-09-22.xml*