Mixed Chebyshev and Legendre polynomials differentiation matrices for solving initial-boundary value problems

1.   Introduction

The ordinary differential equations (ODEs) pose many engineering ^[1], biology ^[2,3], physics ^[4,5,6], and fluid mechanics ^[7,8] problems. Its applications have grown in prominence and relevance over the last several decades, owing primarily to its proven applicability in a wide range of seemingly disparate and vast disciplines. Recently, many studies have been introduced to solve such boundary value problems (BVPs) numerically ^[9,10,11], including the quintic non-polynomial spline solutions, the Galerkin algorithm, the quartic spline method, the cubic spline method, the variational iteration technique, and the sinc-Galerkin method ^{[12,13,14,15,16,17]}.

Different types of polynomials are applied via other techniques to approximate the solutions of the differential equations. The Chebyshev-spectral method is used to solve ODEs in ^[18]. In ^[19], the authors solved ODEs using the Legendre wavelet basis. In addition, Ultraspherical polynomials with the pseudo-Galerkin method are used in ^[20]. In ^{[21,22,23,24,25,26]}, the authors investigated the derivatives of Legendre (LPs), Chebyshev polynomials (CHPs) and monic CHPs, respectively. In addition, methods for higher-order ODEs have been developed in ^[27].

The pseudospectral approach, unlike the finite differences and finite element methods, is unique for particular classes of solving ODEs. The pseudospectral method performs well; it saves several orders of magnitude in computer memory se well as time, ^[28,29]. The pseudospectral method has more specialty than the other methods. When we generate the weight, orthogonality, and nodes, we directly construct the differentiation matrices (D-matrices). So, we have a direct method to find the derivative of any order for functions. In ODEs, the unknown function can be found easily without any complicated steps. Thus, the pseudospectral method can be classified as a significantly important method for solving ODEs. D-matrices were expanded based on the relationship between the coefficients of derivatives and those of the function itself. Many authors introduced the D-matrices from the perspective of the pseudospectral method ^{[30,31,32,33]}.

The spirit of the spectral and pseudospectral methods is the choice of the used polynomials. Many authors are trying to construct, introduce, and modify new polynomials for the original polynomials. The authors in ^[34] defined the monic CHPs. While in ^[35], the authors introduced modified shifted LPs. Similarly, modified shifted CPs are investigated in ^[36]. Other proposals are presented in ^[37] as the derivatives of the polynomials.

Consequently, we continued in this direction and investigated modified basis functions. The main focus is on constructing two sets of polynomials generated from CHPs and LPs. As is known, in the spectral methods, the unknown function is considered to be a sum of unknown constants multiple by selected basis functions. Herein, the introduced basis functions will be alternating between CHPs and LPs. First, the spectral's sum is chosen as CHPs of even degrees, and LPs of odd degrees are applied. The newly generated polynomials will be called Chebyshev-Legendre polynomials (CH-L)Ps. Additionally, recurrence and other relations are investigated for the (CH-L)Ps. Moreover, the associated Gauss-Lobatto quadrature (GLQ) zeros and weights dependent on the Chebyshev GLQ (CH-GLQ) and Legendre GLQ (L-GLQ) are determined. Consequently, the new D-matrices ((CH-L) D-matrices) is built. Similarly, D-matrices are constructed for the LPs even terms and CHPs odd terms (i.e., Legendre-Chebyshev polynomials (L-CH)Ps).

Our paper consists of seven sections. In Section (2), we mention the main definitions and relations that we use. In Section (3), we generate the mixed basis functions and their properties, like the recurrence relation, corresponding weights and some relations, that are used in our method. In Section (4), two D-matrices are formed with two GLQs. In Section (5), the error analysis determines the error order of the presented method. The method of obtaining the solution through the use of our techniques to solve ODEs is shown in Section (6). Finally, we solve some problems in Section (7).

2.   Preliminaries

Through this section, a brief of the needed concepts definitions is presented. The CHPs, $T_n(\xi)$; $n = 0, 1, \cdots$, are the solutions of the following Chebyshev's differential equation ^[28,34]:

where, $\xi \in[-1, 1]$.

Similar results can be obtained via the following recurrence formula:

where $T_{0}(\xi) = 1$, and $T_{1}(\xi) = \xi$.

CHPs form an orthogonal set as follows:

where $\delta_{nm}$ is the Kronecker delta, $c_{0} = 2$ and $c_{n} = 1$ for $n > 0$.

CHPs are bounded according to the following property:

The $N+1$ CH-GLQ points are given by

and the quadrature weights are given by

where $\theta_{0} = \theta_{N} = \dfrac{1}{2}$, $\theta_{j} = 1$, and $0 < j < N$. The discrete inner product of the CHPs is defined as:

The LPs, $P_n(\zeta)$ ; $n = 0, 1, 2, \cdots$, are the solutions of the following Legendre differential equation ^[38]:

where $\zeta\in [-1, 1]$.

LPs can be generated by using the following relation:

with the initial $P_0 (\zeta) = 1$ and $P_1 (\zeta) = \zeta$, which are orthogonal by satisfying the following orthogonal relation:

The boundedness of LPs can be given by:

Unlike CH-GLQ points, the $N+1$ L-GLQ points cannot be obtained exactly. Thus, it is necessary to numerically solve the equation $(1-\zeta^2)P'_N(\zeta) = 0$ to get the points. The general form for the GLQ weights is given by

and the discrete inner product is given by

As in the spectral method, the function, $f(\xi)$, in the pseudospectral method can be expanded as follows ^[38]:

such that $\{\alpha_{k}\}_0^N$ denotes arbitrary constants and $\{q_{k}\}_0^N$ is a set of orthogonal polynomials.

Using the discrete inner product with $\{\xi _{j}, w_{j}\}_{j = o}^{N}$ as the associated GLQ points and weights we get

Substituting Eq (2.14) into Eq (2.13) gives

This approximation is actually represented not by its coefficients but by its values of $f(\xi _{j})$ at $N+1$ GLQ points $\xi _{j}, j = 0, 1, 2...N$.

In the pseudospectral D-matrices, we want to evaluate different derivatives of the approximation given by Eq (2.15):

where

Equation (2.16) can be written in the following matrix form:

The matrices $\{D^{(n)}:n = 1, 2, ...\}$ are called the D-matrices.

3.   Chebyshev and Legendre polynomials mixed basis function

In this section, we shall define a mixed basis function constructed from CHPs and LPs.

Definition 3.1. The set $\{\phi_j(z)\}_{j = 0}^{N}$ of mixed polynomials that are constructed by alternating between CHPs and LPs is called a set of Chebyshev-Legendre polynomials (CH-L)Ps if

Definition 3.2. The set $\{\psi_j(z)\}_{j = 0}^{N}$ of mixed polynomials that are constructed by alternating between LPs and CHPs is called a set of Legendre-Chebyshev polynomials (L-CH)Ps if

Now, the recurrence and some essential relations will be investigated.

3.1.   Recurrence relations

Since two polynomials will be used, we shall modify the recurrence relations of CHPs and LPs.

Lemma 3.1. Let $\{T_j(\xi)\}_{j = o}^{n}$ denote CHPs. Then,

Proof. It is straightforward by using Eq (2.2).

Lemma 3.2. Let $\{P_j(\zeta)\}_{j = o}^{N}$ denote LPs. Then,

Proof. The proof is straightforward by using Eq (2.8).

The following will present the recurrence relations for $\{\phi_n(z)\}_{n = 0}^{N}$ and $\{\psi_n(z)\}_{n = 0}^{N}$.

Corollary 3.1. Let $\{\phi_n(z)\}_{n = 0}^{N}$ denote the (CH-L)Ps defined by Definition (3.1). Then

where $g =$GCD $(n, 2)$, $\phi_0(z) = 1$, $\phi_1(z) = z$, $\phi_2(z) = 2z^{2}-1$ and $\phi_3(z) = \dfrac{1}{2}(5z^{3}-3z)$.

Proof. If $n$ is even, then $\delta_{g1} = 0$. That transforms Eq (3.5) to Eq (3.3). For the odd case for $n$, $\delta_{g1} = 1$. With simple calculations Eq (3.5) takes the form of Eq (3.4).

We can estimate the following lemma by using similar steps as for the previous lemma.

Corollary 3.2. Let $\{\psi_n(z)\}_{n = 0}^{N}$ denote the (L-CH)Ps defined by Definition (3.2). Then

where $g =$GCD $(n, 2)$, $\psi_0(z) = 1$, $\psi_1(z) = z$, $\psi_2(z) = \dfrac{1}{2}(3z^{2}-1)$ and $\psi_3(z) = 4z^{3}-3z$.

3.2.   Weights and orthogonal relations

This section will discuss the orthogonal relationship between CHPs and LPs. This requires investigating important integration-type relations between the two novel constructed polynomials $\{\phi_j(z)\}_{j = 0}^{N}$ and $\{\psi_j(z)\}_{j = 0}^{N}$.

Recall that we let $f(\xi)$ be an integral odd function. Then, the integration $\int_{-a}^{a}f(\xi)d\xi$ is equal to zero. Moreover, if $n$ is even (odd), then the CHPs and LPs are even (odd) functions. Consequently, the following corollary will be fulfilled.

Corollary 3.3. Let $\{T_n(z)\}_{n = 0}^{N}$ and $\{P_m(z)\}_{m = 0}^{N}$ denote CHPs and LPs, respectively. Then

where

or

By using the orthogonal properties of CHPs and LPs and Corollary (3.3), it will be an easy task to prove the following theorems.

Theorem 3.1. Let $\{\phi_n(z)\}_{n = 0}^{N}$ denote a set of (CH-L)Ps. Then, $\{\phi_n(z)\}_{n = 0}^{N}$ satisfies the following:

where $i = 1, 2, \cdots$.

Theorem 3.2. Let $\{\psi_n(z)\}_{n = 0}^{N}$ denote a set of (L-CH)Ps. Then, $\{\psi_n(z)\}_{n = 0}^{N}$ satisfies the following:

where $i = 0, 1, 2, \cdots$.

The following section is devoted to establishing and constructing the pseudospectral D-matrices that use the investigated novel (CH-L)Ps, $\{\phi_n(z)\}_{n = 0}^{N}$ and (L-CH)Ps, $\{\psi_n(z)\}_{n = 0}^{N}$, as trial functions.

4.   Pseudospectral Chebyshev and Legendre differentiation matrices

The introduced mixed polynomials are classified as (CH-L)Ps and (L-CH)Ps. Hence, two main D-matrices can be constructed.

4.1.   Chebyshev-Legendre differentiation matrices

Consider the approximation of $f(z)$ as a summation of the basis function $\phi_n(z)$:

In the following subsections, the used zeros will be CH-GLQ points ($\{\xi_i\}_{i = 0}^{N}$) and L-GLQ points ($\{\zeta_i\}_{i = 0}^{N}$).

4.1.1.   Chebyshev-Gauss-Lobatto points

Lemma 4.1. Consider the expansion given by Eq (4.1) associated with $N+1$ CH-GLQ points. Then, the spectral expansion constants can be given as follows:

where

and

Proof. The expansion given by Eq (4.1) can be written as follows:

Then, by using Theorem (3.1), $a_{2n}$ and $a_{2n-1}$ can be determined to get the following:

Comparing the above equation with Eq (4.3) completes the proof.

Theorem 4.1. Let $f(\xi)$ be real valued function that satisfies Lemma (4.1). Then, its derivative can be obtained by

where $D1_{Ch}f(\boldsymbol \xi) = \left(f'(\xi_0), \, f'(\xi_1), \, \cdots, \, f'(\xi_N)\right)^T$, $f(\boldsymbol \xi) = \left(f(\xi_0), \, f(\xi_1), \, \cdots, \, f(\xi_N)\right)^T$, $f(\boldsymbol \zeta) = \left(f(\zeta_0), \, f(\zeta_1), \, \cdots, \, f(\zeta_N)\right)^T$,

Proof. By using Lemma (4.1) we get

Then, by differentiating the above equation and collocating it by CH-GLQ points we get

where:

The abbreviation "Ch" in Eq (4.4) denotes that the points used are CH-GLQ points.

4.1.2.   Legendre-Gauss-Lobatto points

Similar to the above subsection, but using L-GLQ $\{\zeta_j\}_{j = 0}^{N}$ points instead of CH-GLQ points, we have

where $D1_{Lg}f(\boldsymbol \zeta) = \left(f'(\zeta_0), \, f'(\zeta_1), \, \cdots, \, f'(\zeta_N)\right)^T$,

Alternatively, it may be written as follows:

The abbreviation "Lg" denotes that the points used are L-GLQ points.

4.2.   Legendre-Chebyshev differentiation matrices

In this subsection, the (L-CH)Ps will be used to construct another two matrices. To avoid redundant work, we shall only define the symbols for the matrices. The matrices $D2Ch$ and $D2Lg$ represent the (L-Ch) D-matrices using CH-GLQ points and L-GLQ points, respectively.

where,

Also,

where,

4.3.   Test functions

Herein, the investigated matrices will be used as differentiation tools. Consider the following functions:

$f_1(z) = z^2$, $f_2(z) = e^{z}$ and $f_3(z) = sin{\pi z}$. Tables 1–3 show the maximum absolute error (MAE) for the derivatives of the functions $f_1(z)$, $f_2(z)$ and $f_3(z)$ corresponding to different values of $N$. The results were compared with the obtained results by using the CHPs D-matrix (DCh) and LPs D-matrix (DLg).

Before starting the problem formulation, the method of obtaining the solution and the numerical simulation, the expansion's convergence must be guaranteed.

5.   Error analysis

In this section, the investigation of the error analysis of the presented method will be proved. This will prove the accuracy of our approach.

Lemma 5.1. ^[39] Consider the function $f(z)$ that satisfies the following:

(1) $f(k) = a_k$.

(2) $f(z)$ is a positive, continuous, and decreasing function for $z \geqslant k$.

(3) $\sum a_j$ is convergent and $R_k = \sum_{j = k+1}^{\infty}a_n;$

then,

Lemma 5.2. The (CH-L)Ps, "$\phi_{n}(z)$" and the (L-CH)Ps, "$\psi_{n}(z)$" are bounded such that

where $n \ge 0$.

Proof. It is straightforward by using Eqs (2.4) and (2.10).

Theorem 5.1. If a continuous function $f(z)$, $z\in [-1, 1]$, and $|f^{(q)}(z)| < R$ can be expanded as an infinite sum of (CH-L)Ps or (L-CH)Ps as $f(z) = \sum\limits_{n = 0}^{\infty}a_{n}\phi_{n}(z)$ or $f(z) = \sum\limits_{n = 0}^{\infty}a_{n}\psi_{n}(z)$, then

where $q\in \mathrm{N}$.

Proof. Let $f(z)$ be a function that can be expressed as an infinite sum of (CH-L)Ps, i.e., $f(z) = \sum\limits_{n = 0}^{\infty}a_{n}\phi_{n}(z)$. Then according to Eq (3.10),

As special case of Theorem 2 in ^[36] and Theorem 4 in ^[35] we have

A similar result can be proved by using (L-CH)Ps.

Theorem 5.2. Suppose that $f(z)$ follows the assumptions of Theorem (5.1). Then,

where $f_N(z) = \sum\limits_{n = 0}^{N}a_n\, \, \phi_n(z)$.

Proof. Since $f(z) = \sum\limits_{n = 0}^{\infty}a_n\, \, \phi_n(z)$ and $f_N(z) = \sum\limits_{n = 0}^{N}a_n\, \, \phi_n(z)$,

Using Lemma (5.1) and Theorem (5.1) we have

The next theorem concerns the stability of the error, which estimates the error propagation.

Theorem 5.3. Let $f_N$ and $f_{N+1}$ be two successive approximations of the function $f(z)$ that satisfies Theorem (5.2). Then

Proof. Since $f(z)$ satisfies Theorem (5.2), it follows that

and

Thus

6.   Problem formulation and the method of obtaining the solution

The proposed method for solving ODEs with the presented matrices will be introduced in this section. Consider the ODEs as follows:

for $z\in[-1, 1]$, subject to the following $n$ initial and boundary conditions:

where $s = 0, 1, \cdots, l$, $m = 0, 1, \cdots, r$, $k, n, l, r, \in \mathrm{Z}^{+}$, $c_{k}(z)$ denotes real functions of $z$, and $b_s, e_m \in \mathrm{R}$.

The next subsection will be devoted to solving the BVPs given by Eqs (6.1) and (6.2) using the (CH-L) D-matrices.

6.1.   Chebyshev-Legendre differentiation matrices for solving ordinary differential equations

As mentioned in Eqs (4.10) and (4.17), the first derivative of the unknown function can be written as:

where

The second order derivative can be approximated by ${\bf D_{(CH-L)}}^{(2)} = {\bf D_{(CH-L)}}\cdot{\bf D_{(CH-L)}}$. So, the n$^{th}$ derivative can be written as

By applying Eq (6.5) to the BVPs given by Eqs (6.1) and (6.2), we have

Consider Eqs (4.10), (4.17), (6.4) and (6.5); Eqs (6.6)–(6.8) will be transformed into an algebraic system of the unknown function $f(z)$. This system can be solved by any solver.

6.2.   Legendre-Chebyshev differentiation matrices for solving change: ordinary differential equations

Similar procedures can be executed by using (L-CH)Ps. However, instead, we will use the following:

where $d21_{i, j}$, $d22_{i, j}$, $d23_{i, j}$ and $d22_{i, j}$ are defined in Eqs (4.20, 4.21, 4.24, 4.25). Hence, the derivative of the unknown function can be written in the form:

where:

7.   Numerical examples

This section applies the D-matrices to several ODEs. Then, comparisons with exact solutions and other numerical methods are made. These comparisons show the efficiency of our mixed matrices.

These simulations were conducted by using an Intel® Core$^{TM}$ i7-4500 CPU @ 1.80GHz, 2.40 GHz, and with a SSD hard disk. The software used are MATLAB R2013a and Mathematics 11.

Example 7.1. Consider the following fourth-order equation ^[40]:

with the exact solution

One of the applications of ODEs involves studying the transverse vibration of a uniform beam. This specific scenario can be effectively modeled by using a fourth-order ODE, which describes the relationship between the beam's deflection and the forces acting upon it. By applying the presented matrices to the equation, the results were obtained as reported in Table 4. The table shows that the presented matrices are more accurate and more efficient than the method in ^[40]. The stability of the approximate solution is reported in Figures 1 and 2.

Example 7.2. Consider the obtained fourth-order equation ^[41]:

with the exact solution

According to Table 5, the presented matrices achieved higher accuracy at $N = 8, 12$, and almost the same result at $N = 14$. Those results showed the efficiency of the constructed matrices. A comparison of the exact and approximate solutions is reported in Figure 3 for $N = 14$.

Example 7.3. Consider the following nonlinear Ricatti equation ^[42]:

with exact solution

A point-wise absolute error comparison is shown in Table 6 and Figure 4 for this problem at $N = 12$.

Example 7.4. The following nonlinear system describes the effect of COVID-19 ^[26,43,44]:

and

where

Figure 5 shows that the presented method is in agreement with the methods in ^[26,44]. Table 7 reports the maximum residual error for the eight differential equations comprising Eq (7.1) at different values of $N$.

8.   Conclusions

From the perspective of the pseudospectral spectral method, matrices for differentiation have been constructed. A mix between CHPs and LPs is used in these matrices as basis functions. We show that alternating between both polynomials' even and odd degrees generates two novel polynomials. Then, recurrence and some essential relations for the mixed polynomials are presented and proved. However, the study showed that the generated polynomials are not orthogonal. Thus, quasi-orthogonal relations have been established. Hence, two D-matrices have been created for each set of basis functions. Then, as usual for the D-matrices, the derivatives of three test functions have been calculated, and the results have been compared with the exact differentiation, Chebyshev differentiation matrix and Legendre differentiation matrix. Moreover, several theorems for error analysis have been developed to ensure the correctness of the new expansions. Finally, numerical simulations for linear and non-linear BVPs were conducted to examine and verify the correctness of the investigated matrices. The promising results encourage us, as future work, to use the introduced matrices to solve integro-differential equations, partial differential equations and optimal control problems. Also, the shifted matrices may be introduced to deal with the fractional cases.

Use of AI tools declaration

The authors declare that they have not used artificial intelligence tools in the creation of this article.

Acknowledgments

Princess Nourah bint Abdulrahman University Researchers Supporting Project number (PNURSP2023R157), Princess Nourah bint Abdulrahman University, Riyadh, Saudi Arabia. The authors are also sincerely grateful to the Helwan School of Numerical Analysis in Egypt (HSNAE) members for their valuable support during this research. This study is supported via funding from Prince Sattam bin Abdulaziz University project number (PSAU/2023/R/1444). The authors are thankful to the Deanship of Scientific Research at University of Bisha for supporting this work through the Fast-Track Research Support Program.

Conflict of interest

The authors declare that they have no conflict of interest.