Pseudospectral method for fourth-order fractional Sturm-Liouville problems

1.   Introduction

The Sturm-Liouville problems with integer derivatives have been a significant field of research for centuries, due to their importance in science, engineering, and mathematics ^[3]. This equation has many applications in various fields of science and engineering, including sea-breeze flow ^[23], unidirectional pressure-driven flow of a second-grade fluid in a plane channel with impermeable solid walls ^[10], and flow of the antarctic circumpolar current in rotating spherical coordinates ^[24]. It is crucial to determine both the eigenvalues and the corresponding eigenfunctions for this equation, as they play a very significant role in theory and applications. It is often difficult, if not impossible, to accurately determine the eigenvalue. Numerical methods can be extremely useful to achieve this objective, including boundary value methods ^[18], Legendre-Galerkin-Chebyshev collocation method ^[15], Numerov's method ^[5], the Chebyshev collocation method ^[36], etc.

When dealing with a fourth-order Sturm-Liouville equation, implementing four boundary conditions makes the problem more complex. Dealing with two boundary conditions at each end of the computational domain can be a challenging problem in certain applications. For instance, an eigenvalue problem of a fourth-order differential equation is derived by analyzing the free lateral vibration of a homogeneous beam using the Euler-Bernoulli beam theory ^[35]. A fourth-order Sturm-Liouville equation is generally considered as

appropriate boundary conditions on $w$, $w'$, $pw''$ or$/$and $(pw'')'-qy'$. Here, $p(t)$, $q(t)$, $r(t)$ and $w(t)$ are piecewise continuous functions, and $p(t), r(t) \geq 0$ ^[8].

Several numerical solutions have been found for calculating eigenvalues of fourth-order Sturm-Liouville problems, including the extend sampling method ^[12], the Adomian decomposition method ^[8], the Boubaker polynomial expansion scheme ^[37], variational iteration methods ^[31], the Haar wavelet method ^[34].

Fractional calculus, which generalizes the derivative of a function to non-integer order, has remained a mystery to mathematicians for 300 years. The branch's origin dates back to a 1695 letter from Leibniz to L'Hospital. It was reported in the nineteenth century that some theoretical work was related to fractional calculus. In recent decades, there has been an increasing interest in fractional calculus due to its applications in various fields of physics and engineering. For comprehensive reviews, refer to Podlubny ^[29], Oldham and Spanier ^[28], and Miller and Ross ^[27]. Fractional differential equations arise in various physical phenomena, including mathematical biology, fluid mechanics, electrochemistry, and viscoelasticity ^[6,13,25,26]. Several analytical and numerical methods have been implemented and developed in recent years to solve the fractional differential equation, such as the multiwavelet Galerkin method ^[7,33], Adomian decomposition ^[16], the Kuratowski MNC technique ^[9], the B-spline collocation method ^[22], least-squares finite element ^[17], etc.

The objective of this paper is to present a simple and accurate method, based on the pseudospectral method, to find the eigenvalues of

with the boundary conditions

where $s_{k, j}$ for $k, j = 0, 1, \ldots, 3$ are given real constants, $q_i(t)$, $i = 1, \ldots, 3$, and $r(t)$ are in $L^1[0, 1]$, and $\mu_i$ for $i = 1, \ldots, 4$ are real numbers such that $\mu_i\in (i-1, i]$. The constant $S$ indicates the number of conditions in point 0. $\mathcal{D}$ and ${}^c\mathcal{D}_{0}^{\mu}$ indicate the derivative operator and the Caputo fractional derivative, respectively. In this problem, $\lambda$ is called an eigenvalue.

It is important to note that fractional Sturm-Liouville problems (1.1)-(2.13) arise frequently when dealing with separable linear fractional partial differential equations ^[4]. In ^[19], the authors utilized the series solution to solve this equation. This is the first report on fourth-order fractional Sturm-Liouville problems with varying coefficients. The presented work can be considered the second paper in this field, which solves the problem with a simpler method and more accurately.

The rest of the paper is structured in the following manner: Chebyshev Cardinal polynomials and their properties are reviewed and introduced in Section 2. In Section 3, the pseudospectral method is applied to solve fourth-order fractional Sturm-Liouville problems using Chebyshev cardinal polynomials. Section 4 is dedicated to illustrate the applicability and accuracy of the method. To sum up our work, we have included a conclusion in Section 5.

2.   Chebyshev cardinal polynomials

Given $N\geq 0$. Let $\mathcal{R}: = \{r_j:T_{N+1}(r_j) = 0, j\in \mathcal{N}\}$ be the set of the roots of the TChebyshev polynomial $T_{N+1}$ in which $\mathcal{N}: = \{1, 2, \ldots, N+1\}$. Recall that the TChebyshev polynomials are defined on $[-1, 1]$ by

and their roots are specified by

Shifted TChebyshev polynomials for generic intervals $[a, b]$ are related to the TChebyshev polynomials by

and the roots of $T^*_{N+1}$ in its turn are obtained by $t_j = \frac{(r_j+1)(b-a)}{2}+a$, $j\in \mathcal{N}$.

The Chebyshev cardinal function (CCF) is one of the orthogonal polynomials' most notable cardinal functions ^[1,11,32]. Considering $T^*_{N+1, t}(t_j)$ as the derivative of function $T^*_{N+1}(t)$ with respect to the variable $t$, Chebyshev cardinal functions can be denoted by

The most striking feature of these polynomials is their cardinality, i.e.,

in which $\delta_{ji}$ indicates the Kronecker delta. This property is mostly important as it enables us to approximate any function $w\in \mathcal{H}^\alpha([a, b])$ (the Sobolev space $\mathcal{H}^\alpha([a, b])$ will be briefly introduced) easily and without integration in finding the coefficients, viz,

In what follows, since we will need the definition of Sobolev spaces and their norm, we will provide a brief definition of it. For $\alpha \in \mathbb{N}$, we denote by $\mathcal{H}^\alpha([a, b])$ the sobolev space of functions $w(t)$ which have continuous derivatives up to order $\alpha$ such that $\mathcal{D}^{\beta}w\in L^2([a, b])$:

with the norm

and the semi-norm

Lemma 2.1. (cf ^[14]) Given $N\geq 0$, the error of approximation (2.5), obtained using the shifted Chebyshev nodes $\{t_j\}_{j\in\mathcal{N}}$, can be bounded

where the constant $C$ is independent of $N$. Furthermore, it can be verified that

2.1.   Operational matrix of derivative

Let $\Psi(t)$ be a vector function with entries $\{\psi_j\}_{j\in \mathcal{N}}$. We specify the operational matrix of derivatives for CCFs as

To evaluate the elements of $D$, they can be obtained via the following process using the approximation (2.5). It follows from (2.5) that

It is worth noting that there is another presentation of CCFs ^[2]

where $\varrho = 2^{2N+1}/((b-a)^{N+1}T^*_{N+1, t}(t_j))$. When the operator $\mathcal{D}$ acts on both sides of (2.12) and taking into account (2.3), we obtain

It can be shown by (2.11) and (2.13) that

2.2.   Operational matrix of fractional integration

Considering the interval $[0, 1]$, the fractional integral is defined as

where $\Gamma(\mu)$ denotes the Gamma function.

Note that there is a square matrix $I_{\mu}$ such that the acting of the fractional integral operator on $\Psi(x)$ can be represented by it, viz,

It is straightforward to show that the elements of this matrix can be obtained by

After performing some simple calculations, it can be inferred from ^[30] that

in which

and

Motivated by (2.12), the CCFs can be determined by

Using this definition of CCFs, (2.16) leads to

So it can be concluded from (2.16) that

2.3.   Operational matrix of fractional derivative

Definition 2.1. ^[21] Let $\mu\in \mathbb{R}^+$ and $m: = \lceil \mu \rceil \in \mathbb{N}$ ($\lceil.\rceil$ denotes the ceiling function). The Caputo fractional derivative is denoted by

where $\mathcal{D}^m: = \frac{d^m}{dt^m}$.

Lemma 2.2. [21, cf Corollary 2.3(a)] Let $\mu\in \mathbb{R}^+$, $m: = \lceil \mu \rceil \in \mathbb{N}$ and $\mu\not \in \mathbb{N}_0$. Then we have

Taking into account Definition 2.1 and the operational matrices of derivative $D$ and fractional integral $I_{\mu}$, when the Caputo derivative operator acts on $\Psi(x)$, it follows that

So, the operational matrix for the Caputo operator is specified by

3.   Materials and methods

The present chapter will be focused on solving fourth-order fractional Sturm-Liouville equations (SLEs) with the Caputo operator using an efficient and accurate scheme based on the pseudospectral method. As mentioned above, we consider the fourth-order fractional Sturm-Liouville equation (SLE)

equipped with the boundary conditions

where $s_{k, j}$ for $k, j = 0, 1, \ldots, 3$ are given real constants, $q_i(t)$, $i = 1, \ldots, 3$, and $r(t)$ are in $L^1[0, 1]$, and $\mu_i$ for $i = 1, \ldots, 4$ are real numbers such that $\mu_i\in (i-1, i]$. The constant $S$ indicates the number of conditions in point 0. In this problem, $\lambda$ is called an eigenvalue and is not given in advance. Our objective while solving the Sturm-Liouville equation is to determine the eigenvalues.

One of the common schemes that is used to solve a differential equation is based on converting it to an integral equation. To give rise to such a conversion for Eq (3.1)-(3.2), we reformulate it as

Equation (3.3) is written due to the following formula ^[21]:

In this equation, we either have the values of the $w(t)$ function and its derivatives at zero or assume them as unknowns and add them to our unknowns.

Pseudospectral scheme

To obtain the pseudospectral discretization of Eq (3.3), the unknown solution $w$ is approximated using Chebyshev cardinal polynomials as

where the $(N+1)$-dimensional vector $W$ consists of the unknowns $(w_j)_{j = 1}^{N+1}$. Substituting $w_N$ instead of $w$ in (3.3), we have

in which $\bar{w}(t) = \sum_{i = 0}^{3}\frac{w^{(i)}(0)}{i!}(t)^i$. Now, we approximate all terms in (3.4) as follows:

● Using the Chebyshev cardinal polynomials, $\bar{w}$ can be approximated as

where $\bar{W}$ is a $(N+1)$-dimensional vector whose elements may be known or unknown according to the boundary conditions.

● Let us put $g_i(t): = ^c\mathcal{D}_{0}^{\mu_i}(w_N)(t)$, $i = 1:3$. Taking into account the operational matrices of derivative or the operational matrix for the Caputo operator, we obtain

By approximating $q_i(t)\tilde{g}_i(t)$ using CCFs, one can write

where $Q_i$ is a $(N+1)$-dimensional vector whose elements consist of unknowns $w_j$, $j = 1, \ldots, N+1$. Finally, taking the fractional integral from both sides of (3.8) and taking into account the operational matrix of the fractional integral, we obtain

in which $G_i$ is an $(N+1)$-dimensional matrix with constant entries.

● Putting $g_0: = q_0w_N$ and then approximating it using CCFs leads to

in which $Q_0$ is a $(N+1)$-dimensional vector whose elements consist of unknowns $w_j$, $j = 1, \ldots, N+1$, and $G_0$ is an $(N+1)$-dimensional matrix with constant entries. Taking the fractional integral from both sides of (3.10) and motivated by the operational matrix of the fractional integral, it is easy to achieve

● Similar to the previous one, the term $\mathcal{I}_0^{ \mu_4 }(rw_N)(t)$ can be approximated by

where $\tilde{G}$ is an $(N+1)$-dimensional matrix with constant entries.

Let us consider $G: = \sum_{i = 0}^{3} G_i$. By this assumption, the residual function can be determined by

Our goal is to make $r(t)$ as close to zero as possible. By utilizing the shifted Chebyshev nodes $\{t_j\}_{j\in\mathcal{N}}$ and taking into account Eq (2.4), the pseudospectral method results in

Setting $\Upsilon^T: = \left(I+G+\lambda\tilde{G}\right)I_{\mu_4}$, we have

or equivalently, we have

To determine the $\bar{w}$ function, $4-s$ unknowns are undetermined. So, (3.16) consists of $N+1$ equations with $N+5-s$ unknowns. We add $4-s$ equations to this system of equations according to the boundary conditions (3.2), and then we obtain a new system

in which $\hat{W}$ consists of all unknowns, including $w_j$, $j = 1, \ldots, N+1$, and $4-s$ unknowns of function $\bar{w}$. Also, $A(\lambda)$ is a matrix whose elements are functions of $\lambda$. In order for Eq (3.17) to have a non-zero eigenvector, it is necessary that the matrix $A(\lambda)$ is singular when $\lambda$ is an eigenvalue. Equivalently, we have

Note that $det(A(\lambda))$ refers to the characteristic polynomial, and $\lambda$ represents its root. To calculate the roots, the Maple software can be used.

It is worth realizing that the eigenvector $\hat{W}$ associated with $\lambda$ belongs to

Since $\hat{W}$ must be nonzero, the matrix $A(\lambda)$ has a nonzero kernel. On the other hand, since we obtain the eigenvalues approximately, the characteristic polynomial $det(A(\lambda))$ does not become exactly zero in these eigenvalues. Still, it will have a value very close to zero. To find the eigenvector we are looking for, we need to identify the eigenvalue of $A(\lambda)$ that is the smallest and then select the eigenvector corresponding to that eigenvalue. Finally, we obtain the eigenvector corresponding to the eigenvalue $\lambda$, via

4.   Numerical results

Example 4.1. As the first illustrative example, the fourth-order fractional SLE is considered

with conditions

The exact eigenvalues can be calculated by $\lambda_l = (l\pi)^4$ ^[19,20] for $\mu = 4$.

Table 1 is tabulated to demonstrate the approximation of the first three eigenvalues for different values of $\mu$ when the value of $N$ is fixed. As we vary $\mu$ from $3.7$ to 4, we aim to verify that the eigenvalues approach the eigenvalues for $\mu = 4$ when $\mu$ values approach the exact ones for 4. What is important is that eigenvalues approach the exact ones, when $\mu$ tends to 4.

Table 2 illustrates the convergence of the method. To this end, we simulate the method for different values of $N$ when the $\mu$ value is fixed. Motivated by this, we can deduce that the presented scheme is accurate.

For more illustration, the approximate eigenfunctions corresponding to the first 3 eigenvalues are plotted in Figure 1. This figure shows that for each eigenvalue $\lambda_l$, the corresponding eigenfunction $w_l$ has $l-1$ zeros.

Example 4.2. For the second example, the fourth-order fractional Sturm-Liouville equation

with conditions

is considered ^[20,37].

Table 3 presents the numerical results for the first three eigenvalues. In this table, we can find the results of the differential quadrature method ^[37] and the polynomial expansion and integral technique ^[20] with the presented technique results and compare them.

Table 4 illustrates the convergence of the method. To this end, we simulate the method for different values of $N$ when the $\mu$ value is fixed.

Finally, for more illustration, the approximate eigenfunctions corresponding to the first 3 eigenvalues are plotted in Figure 2. This figure shows that for each eigenvalue $\lambda_l$, the corresponding eigenfunction $w_l$ has $l-1$ zeros.

5.   Conclusions

The Sturm-Liouville problem is a significant ordinary differential equation with numerous applications in various fields of science. Thus, solving it and obtaining its eigenvalues and eigenfunctions can be an intriguing task. In this study, we propose a method for solving the fourth-order fractional Sturm-Liouville equation using Chebyshev cardinal polynomials and the pseudospectral method, which is highly efficient. This field has only seen a few numerical methods to solve this type of equation. This paper is the second of its kind. The numerical examples that have been solved confirm the accuracy and efficiency of the method used.

According to the experimental observations, the following can be concluded:

(1) The proposed schemes are effective for solving these types of equations.

(2) The eigenvalues approach the exact ones when $\mu$ tends to 4.

(3) The presented method is convergent.

(4) The abilities of the presented method are simplicity, high accuracy, and reduction of computational cost by avoiding integration in finding coefficients.

In the future, we plan to extend our numerical approaches for solving higher-order fractional Sturm-Liouville equations and fractional Dirac equations.

Author contributions

Haifa Bin Jebreen: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Writing–original draft, Writing–review & editing, Funding acquisition; Beatriz Hernández-Jiménez: Formal analysis, Invesigation, Software, Validation, Investigation, Writing–original draft, Writing–review & editing. All authors have read and agreed to the published version of the manuscript.

Acknowledgments

This project was supported by Researchers Supporting Project number (RSP2024R210), King Saud University, Riyadh, Saudi Arabia.

Conflict of interest

The authors declare that they have no conflicts of interest.