Direct Computation of Operational Matrices for Polynomial Bases

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.

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 solution 𝑓(𝑥) by truncated series approximations, up to order 𝑛 [1–4], such that𝑓(𝑥)=𝑓_𝑁(𝑥)=^𝑁_𝑖=1𝑐_𝑖𝑃_𝑖(𝑥),𝑥∈[𝑎,𝑏]⊂ℜ. (1.1)

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 𝑃_𝑖(𝑥)=𝑥^𝑖,  𝑖=1,2,…,𝑛.

Considering 𝐶 the line vector that contains the coefficients 𝑐_𝑖 and 𝐵 the column vector that contains the base polynomials 𝑃_𝑖, expression (1.1) can be written as considering: 𝑓_𝑁(𝑥)=[𝑐₀𝑐₁,…,𝑐_𝑛][𝑃₀𝑃₁,…,𝑃_𝑛]^𝑇=𝐶𝐵 [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 of 𝑓_𝑁(𝑥) into matrix operations in a Hilbert space [8].

So, defining 𝑀_𝐼 as the operational integration matrix and 𝑀_𝐷 as the operational differential matrix, it is possible to obtain the line vector 𝑉 containing the coefficients of the series that represent the integrated function or the differentiated function by 𝑉^𝑇_𝐼=𝑀_𝐼𝐶^𝑇 and 𝑉^𝑇_𝐷=𝑀_𝐷𝐶^𝑇.

Consequently, ^𝑥_𝑎𝑓_𝑁(𝑥)𝑑𝑥=𝑉_𝐼𝐵,𝑑[𝑓_𝑁(𝑥)]𝑑𝑥=𝑓_𝑁⁽¹⁾(𝑥)=𝑉_𝐷𝐵. (1.2)

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 Theorem 2.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. Considering 𝑍=𝑍_{𝑛+1,𝑛+1} a square matrix describing the resulting coefficients of a linear operation 𝛼 in the generic basis 𝐵_𝐺=[𝑏₀𝑏₁,…,𝑏_𝑛]^𝑇 as a function of the same basis and 𝑓_𝑁(𝑥) a series the coefficients vector of which is 𝐶=[𝑐₀𝑐₁,…,𝑐_𝑛], and 𝑉=[𝑣₀𝑣₁,…,𝑣_𝑛] the line vector of the resulting coefficients of the linear operation 𝛼 applied to the series, then: 𝑉^𝑇=𝑍^𝑇𝐶^𝑇. (2.1)

Proof. Considering 𝑓_𝑁(𝑥)=∑^𝑛_𝑖=0𝑐_𝑖𝐵_𝐺𝑖(𝑥), since 𝛼 is a linear operation: 𝛼[𝑓_𝑁(𝑥)]=∑^𝑛_𝑖=0𝑐_𝑖𝛼𝐵_𝐺𝑖(𝑥).
On the other hand, considering the finite matrix space 𝑀_𝑖𝑗, 𝑖∈{0,1,…,𝑛} and 𝑗∈{0,1,…,𝑛}, 𝛼 is a scalar, depending on 𝑥.
Considering that 𝛼𝐵_𝐺=𝑍𝐵_𝐺 and 𝑉𝐵_𝐺=𝐶(𝑍𝐵_𝐺), observing the matrix and vector dimensions, this last result is a scalar and, consequently 𝐵^𝑇_𝐺𝑉^𝑇=(𝑍𝐵_𝐺)^𝑇𝐶^𝑇, implying that:𝐵^𝑇_𝐺(𝑉^𝑇−𝑍^𝑇𝐶^𝑇)=0. (2.2)
As 𝐵^𝑇_𝐺 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 matrix 𝑍 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 [𝑎,𝑏]=[−1,1], one can observe that, for Legendre polynomials [9]^𝑥₋₁𝑃_𝑖𝑃(𝑥)𝑑𝑥=^𝑖+1(𝑥)−𝑃_𝑖−1(𝑥)2𝑖+1,𝑖=1,2,…,𝑛,_𝑥−1𝑃₀(𝑥)𝑑𝑥=𝑃₁(𝑥)+𝑃₀. (2.3)

Defining𝑃𝑃=⁰𝑃₁⋯0⋅𝑃_𝑛+1,𝑃^∗=𝑃₀𝑃₁⋯𝑃_𝑛+1, (2.4) one can write ∫^𝑥₋₁𝑃𝑑𝑥=𝑍_𝐼𝑃^∗, with the integral acting over to the elements of the vector. The null coefficient on the (𝑛+1)-order term assures equal dimensions for the vector to be integrated and the vector that results from the process.

Following (2.3), 𝑍_𝐼 is a square matrix 𝑍𝑛+2, 𝑛+2 and:(i)𝑍_1,1=𝑍_1,2=1;(ii)𝑍_𝑘,𝑘+1=−1/(2𝑘−1),  𝑍_{𝑘,𝑘−1}=1/(2𝑘−1),  𝑘=2,3,…,𝑛+1;(iii)𝑍_{𝑛+2,𝑛+1}=−1/(2(𝑛+2)−1)=−1/(2𝑛+3); (iv)there is no 𝑍_{𝑛+2,𝑛+3} term, because 𝑃_𝑛+1 is not integrated.

Since the operational matrix 𝑀_𝐼 is the transpose of 𝑍_𝐼, when writing a code to implement a computational algorithm, it is required to exchange the indices in the expressions above, obtaining 𝑀_𝐼. directly.

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

From an operational matrix expressed in a generic reference basis 𝑋, it is possible to obtain the corresponding operational matrix in another basis 𝐺, also generic, by using a sequence of simple matrix operations, as Theorem 3.1 states.

Theorem 3.1 (“Sandwich Matrix” (Ω)). Given a generic polynomial basis 𝐵_𝐺 in the interval [𝑎,𝑏], the matrix 𝑍_𝐺 of the operational matrix theorem is obtained by 𝐺Ω𝐺⁻¹, where the “sandwich matrix” is Ω=𝑋⁻¹𝑍_𝑋𝑋, with 𝐺_{𝑛+1,𝑛=1} and 𝑋_{𝑛+1,𝑛+1} being the matrices that describe the generic polynomials (𝐺) and the reference ones (𝑋) in the canonic basis, respectively.

Proof. Considering that 𝛼 is a linear operation: 𝛼𝐵_𝐺=𝛼[𝐺𝐵_𝐶]=𝐺𝛼[𝐵_𝐶], where 𝐵_𝐶 is the canonic base 𝐵_𝐶=[1𝑥𝑥²,…,𝑥^𝑛]^𝑇, since the canonic base can be written as a function of the reference base as 𝐵_𝐶=𝑋⁻¹𝐵_𝑋, it can be concluded that 𝛼𝐵_𝐺=𝛼[𝐺𝐵_𝐶]=𝐺(𝑋⁻¹𝛼[𝐵_𝑋]).
As 𝛼[𝐵_𝑋]=𝑍_𝑋𝐵_𝑋, it implies that 𝛼[𝐵_𝑋]=𝐺(𝑋⁻¹)𝑍_𝑋𝐵_𝑋. Now, it is necessary to transform 𝐵_𝑋 into the generic basis 𝐵_𝐺, in order to express the result as a function of the generic input base.
From the defined matrices 𝐵_𝑋=𝑋𝐵_𝐶 and 𝐵_𝐶=𝐺⁻¹𝐵_𝐺, the expression given above is written as: 𝛼𝐵_𝐺=𝐺(𝑋⁻¹𝑍_𝑋𝑋)𝐺₋₁𝐵_𝐺=𝐺Ω𝐺⁻¹𝐵_𝐺 and 𝐺 must be nonsingular.
Finally, since 𝛼𝐵_𝐺=𝑍_𝐺𝐵_𝐺, the generic matrix 𝑍_𝐺 is obtained by: 𝑍_𝐺=𝐺Ω𝐺⁻¹, where Ω=(𝑋⁻¹)𝑍_𝑋𝑋.

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: 𝑀_𝐺=(𝐺⁻¹)^𝑇Ω^𝑇𝐺^𝑇. Indeed, for any linear operation: 𝑍_𝐺=𝐺(𝑋⁻¹𝑍_𝑋𝑋)𝐺⁻¹ and, since 𝑀_𝐺=𝑍^𝑇_𝐺, it can be concluded that: 𝑀_𝐺=(𝐺⁻¹)^𝑇Ω^𝑇𝐺_𝑇.

(iii) Since 𝑋 is arbitrary, the canonic basis can be used, and thus, the matrix describing the elements of 𝑋 as a function of the canonic base is the identity. So, Ω=𝐼⁻¹𝑍_𝐶𝐼=𝑍_𝐶.

(iv) Taking the previous comment into account, the matrix 𝑍_𝐺 of the generic basis can be obtained by: 𝑍_𝐺=𝐺𝑍_𝐶𝐺⁻¹, that is, considering the resemblance definition [9], the matrix of the generic basis 𝑍_𝐺 is similar to or resembling matrix 𝑍_𝐶. To summarizing, the generating matrices of the operational matrices are similarity classes related to each linear operation in the interval [𝑎,𝑏].

(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 Theorem 3.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 𝑍 is built, in order to represent the integrals of the basis as a function of the basis itself. By defining𝑥𝑋=⁰𝑥¹⋯0⋅𝑥^𝑛+1,𝑋^∗=𝑥⁰𝑥¹⋯𝑥^𝑛+1, (3.1) one can write ∫^𝑥₋₁𝑋𝑑𝑥=𝑍_𝐶𝐼𝑋^∗, with the integral acting over to the elements of the vector, that is, ∫^𝑥₋₁𝜉^𝑖𝑑𝜉=(1/(𝑖+1))𝑥^𝑖+1−(−1)^𝑖+1/(𝑖+1),  𝑖=0,1,…,𝑛.

Consequently, ∫^𝑥₋₁𝜉^𝑖𝑑𝜉=(1/(𝑖+1))𝑃_𝑖+1−(−𝑃₀)^𝑖+1/(𝑖+1)=(1/(𝑖+1))𝑃_𝑖+1−(−1)^𝑖𝑃₀/(𝑖+1)=(1/(𝑖+1))𝑃_𝑖+1−(−1)^𝑖/(𝑖+1),  𝑖=0,1,…,𝑛.

The nonnull elements 𝑍_𝐶𝐼=𝑍_{𝑛+2,𝑛+2} are, therefore: 𝑍_𝑘,𝑘+1=1/𝑘 and 𝑍_𝑘,1=(−1)^𝑘−1/𝑘, with 𝑘=1,2,…,𝑛+1. As apparent in comments presented in the former subsection, for each linear operation, in the canonic basis, 𝑍_𝐶=Ω.

Considering, for instance, 𝑛=4, the “sandwich matrix” for integration (Ω_𝐼) is:Ω_𝐼=1110000−1/202100031003−100410004015110005000000. (3.2)

For the interval [0,𝑏], the elements 𝑍_𝑘,1, with 𝑘=1,2,…,𝑛+1 will be equal to zero, that is,Ω_𝐼=1010000002100000031000000401010005000000. (3.3)

Considering the process using the shifted Legendre base up to order 4, known as shifted-Legendre, Ω is: Ω_𝐿𝐼=𝐿⁻¹𝑍_𝐿𝐼𝐿, where 𝑍^𝑇_𝐿𝐼=𝑀_𝐿𝐼 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:−1𝐿=100000010000203230000−205230080−1540358001580−3540638. (3.4)

Calculating Ω_𝐿𝐼=𝐿⁻¹𝑍_𝐿𝐼𝐿Ω_𝐿𝐼=−11100002012100031003−100410004015110005−11356930−0566022. (3.5)

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 (𝑢∈ℜ∶0≤𝑢≤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, matrix 𝑍_𝐷 can be built, describing the derivative elements of the basis as a function of the basis itself. Therefore, rewriting the vector 𝑋^∗ defined in the former subsection: 𝑑𝑋^∗/𝑑𝑡=𝑍𝑋^∗, with the derivative acting over the components of vector 𝑋^∗,  𝑑[𝑥⁰]/𝑑𝑥=0,…,𝑑[𝑥^𝑛]/𝑑𝑥=𝑛𝑥^𝑛+1, for any possible considered domain.

Thus, the nonnull terms of the matrix 𝑍_𝐷=𝑍_{𝑛+1,𝑛+1} are 𝑍_{𝑖,𝑖−1}=𝑖−1, 𝑖=2,…,𝑛 and, according to the last comment from Section 3.1, 𝑍_𝐶𝐷=Ω_𝐷.

Considering 𝑛=4, the differentiation matrix is given byΩ_𝐷=0000010000020000030000040. (3.6)

In order to obtain this matrix from the Legendre basis, the process is analogous to the previous one. Applying Ω_𝐷=𝐿⁻¹𝑍_𝐿𝐷𝐿, where 𝑍^𝑇_𝐿𝐷 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 Chebyshev 𝐼 expansions [11–13], when solving differential equations. In order to help in this task, integration and differentiation Chebyshev 𝐼 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.

Matrix 𝑇 describing the Chebyshev 𝐼 polynomials, as a function of the canonic basis, up to order 4 is:𝑇=1000001000−102000−304010−808. (3.7)

By using Theorem 3.1, the Legendre integration matrix (𝑀_𝐿𝐼) is known, and Chebyshev integration matrix can be obtained, by using the second comment from Section 3.1.

Consider, for instance, 𝑓(𝑥)=−2−10𝑥+6𝑥²+16𝑥³, 𝑥∈[−1,1] to be integrated in the interval [−1,𝑥]. Written as a Chebyshev 𝐼 series: 𝑓(𝑥)=[1234][𝑇₀𝑇₁𝑇₂𝑇₃]^𝑇 and 𝑀_𝑇𝐼=(𝑇⁻¹)^𝑇Ω^𝑇_𝐼𝑇^𝑇, observing that Ω is basis independent.

Calculating the several matrices is the case of the example,Ω𝑇=1000001000−102000−304010−808_𝑇𝐼=11−213−141511000−012121000003029,𝑀0001/40_𝑇𝐼=10−4−1318−15110−210−0121410−40100620−9.0001/80 (3.8) Performing 𝑉=𝑀_𝑇𝐼[1234]^𝑇, with the coefficients of the series representing the integral of this function, that is, 𝐹(𝑥)=∫−1^𝑥𝑓(𝑥)𝑑𝑥=𝑉𝐵_𝑇, where 𝑉=[0−1/2−1/21/21/2] and, consequently, in the canonic form: 𝐹(𝑥)=1−2𝑥−5𝑥²+2𝑥³+4𝑥⁴.

To obtain the Chebyshev differential matrix, the procedure is analogous, giving:Ω^𝑇_𝐷=,𝑀0100000200000300000410000_𝑇𝐷=.0103000408000600000800000 (3.9)

As expected, 𝑀_𝑇𝐷𝑉^𝑇=[12340]^𝑇, with the nonnull terms of the invariant Ω^𝑇_{𝐷(𝑛+1,𝑛+1)} given by Ω^𝑇_{𝐷(𝑖,𝑖+1)}=𝑖,  𝑖=1,…,𝑛.

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:−𝑢^′′(𝑣)+𝑏𝑢^′(𝑣)=0,𝑢∶[0,1]→[0,1], (4.1) with 𝑢(0)=0 and 𝑢(1)=1.

Firstly, in order to have a Jacobi interval domain, we change variables, 𝑣=(𝑥+1)/2, with 𝑥∈[−1,1] obtaining the transformed equation:−4𝑢^′′(𝑥)+𝑏2𝑢^′(𝑥)=0;𝑢∶[−1,1]→[0,1], (4.2) with 𝑢(−1)=0 and 𝑢(1)=1.

If 𝑢^∗(𝑥)=∑^𝑛_𝑖=0𝑐_𝑖𝑃_𝑖 is the series that approximates 𝑢(𝑥) and 𝑀_𝐷 the operational differentiation matrix, one can write: (−4𝑀²_𝐷𝐶^𝑇+2𝑏𝑀_𝐷𝐶^𝑇)^𝑇⋅𝑃(𝑥_𝑘)=0, defining:𝑃(𝑥_𝑘𝑃)=⁰(𝑥_𝑘)𝑃₁𝑥_𝑘⋯𝑃_𝑛(𝑥_𝑘), (4.3) with 𝐶=[𝑐₀𝑐₁⋯𝑐_𝑛] and 𝑘=1,2,…,𝑛−1.

This matrix equation is applied to 𝑛−2 domain points generating a linear algebraic equations system with 𝑛−2 equations and 𝑛 unknown variables. The two missing equations are obtained from the boundary conditions: 𝑥₀=−1 and 𝑥_𝑛=1. To avoid the Runge phenomenon [15], 𝑥_𝑘 are chosen as nodes of the polynomial basis.

Figure 1 shows the exact solution and the obtained by using Legendre approximation and considering 𝑏=1 and 𝑛=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 𝑏=20 and 𝑛=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,\mathrm{\hspace{0.17em}\hspace{0.17em}}n=5$).

Figure 2 

Error in the solution of convection-diffusion equation with Legendre approximation ($b=1,\mathrm{\hspace{0.17em}\hspace{0.17em}}n=5$).

Figure 3 

Solution of convection-diffusion equation with Chebyshev approximation ($b=20\text{;}$$n=25$).

Figure 4 

Error in the solution of convection-diffusion equation with Chebyshev approximation ($b=20,\mathrm{\hspace{0.17em}\hspace{0.17em}}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.