Redefined Quintic B-Spline Collocation Method to Solve the Time-Fractional Whitham-Broer-Kaup Equations

Abstract

This article proposes a collocation approach based on a redefined quintic B-spline basis for solving the time-fractional Whitham-Broer-Kaup equations. The presented method involves discretizing the time-fractional derivatives using an -approximation scheme and then approximating the spatial derivatives using the redefined quintic B-spline basis. The von Neumann technique has been used to demonstrate that the proposed method is unconditionally stable. The error estimates are discussed and show that the proposed method is third-order convergent. The results demonstrate the potential of the proposed method as a reliable tool for solving fractional differential equations.

1. Introduction

Fractional calculus has grown significantly in relevance in recent years. The fractional derivatives and integrals have been used in numerous applications in the fields of science and engineering, including but not limited to fluid mechanics, chemical physics, electricity, control theory, epidemic diseases, biomedicine, signal processing, and issues with heat conduction and diffusion [1–5]. There are several definitions of fractional-order derivatives, each with a variety of uses [6–10].

The Whitham-Broer-Kaup (WBK) equations are a set of coupled nonlinear partial differential equations that describe the propagation of shallow water waves in a channel. Its fractional counterpart describes shallow water in a porous medium, which can absorb wave energy and prevent tsunamis. Several analytical and numerical methods have been developed to solve the WBK equations [11–14]. Wang and Chen [15] applied an analytic iterative technique called the residual power series method to solve time-fractional WBK equations. Wang et al. [14] proposed the generalized exponential rational function method to elucidate the basic solution properties of the WBK equation. Wang et al. [16] used the generalized projective Riccati equation method to solve the classical WBK equations. Yasmin [17] used the Yang decomposition method for fractional-order nonlinear WBK equations. Sadat and Kassem [18] used Lie point symmetries for the fractional Riemann-Liouville system to reduce fractional WBK equations to nonlinear fractional ordinary differential equations using the prolongation theorem. Wang and Li [19] provided a streamlined homogeneous balance technique to investigate the shallow water small-amplitude WBK model equations. Nonlaopon et al. [20] used the Laplace homotopy perturbation transform technique to solve the fractional-order WBK equations. Cao et al. [13] used the conformal fractional derivative to transform the nonlinear space-time fraction WBK equation into an ordinary differential equation and then used the complete polynomial discriminant system to find the exact solutions. Ali et al. [21] applied the Laplace Adomian decomposition technique to obtain an approximate solution of the nonlinear coupled system of WBK equations of time-space fractional order. Shah et al. [22] used the -homotopy analysis transform method and the natural decomposition method to solve time-fractional WBK equations. Our paper focuses on the homogeneous and nonhomogeneous fractional WBK equations that take the following form [23]:

subject to the initial conditions and boundary conditions

where , , , , , , , and are supposed to be smooth functions with continuous first-order derivatives and represents the horizontal velocity, while denotes the height that deviates from the equilibrium position. The constants and are real numbers, which are expressed as different diffusion powers, and is the Caputo derivative operator, where . When , the resulting equations are the usual WBK equations. Importantly, setting and yields the fractional-order modified Boussinesq (MB) equation, while setting and produces the fractional-order approximate long wave (ALW) equation.

The collocation method is widely used to obtain solutions for partial differential equations [24–27]. Depending on the situation, it can be useful to find the solution of fractional partial differential equations (FPDEs) at various locations within the given problem domain. In such cases, spline solutions can provide information on spline interpolation between mesh points. The nonpolynomial, cubic, quadratic, trigonometric, and quintic B-spline methods are used to solve many fractional-order partial differential equations [28–32].

In the usual collocation method, the basis functions are required to vanish on the boundary where the Dirichlet-type boundary conditions are specified. However, in the set of quintic B-splines , the basis functions do not vanish at one of the boundary points. Therefore, it is necessary to redefine the basis functions into a new set of basis functions that vanish on the boundary where the Dirichlet-type boundary conditions are specified. The primary goal of this work is to propose an efficient computational approach based on a redefined quintic B-spline (RQBS) algorithm for obtaining the numerical solution of time-fractional WBK equations. RQBS functions are essentially a generalization of typical quintic B-spline functions that include a free parameter that gives the ability to adjust the solution curve. We used the -approximation formula to discretize the Caputo time-fractional derivative, whereas RQBS functions are used to discretize the spatial derivatives. This approach is developed for numerical solutions of fractional-order WBK equations. Moreover, this scheme is equally effective for homogeneous and nonhomogeneous FPDEs. The redefined quintic B-spline collocation discretization for the problem considered leads to a system with the pentadiagonal matrix.

This paper’s brief outline is as follows. In Section 2, we provide some basic definitions and lemmas. In Section 3, we explain the quintic B-spline collocation scheme and its redefinition. Then, we describe the method and apply it to the coupled time-fractional WBK equation. Section 4 discusses the von Neumann technique for ensuring the stability of the method. Error analysis is discussed in Section 5. In Section 6, numerical examples are presented to demonstrate the applicability and accuracy of the proposed method. Finally, we finished this paper with the conclusion.

2. Basic Concepts

In this section, we will introduce some fundamental definitions of the fractional derivative of order where . Numerous definitions of the fractional derivative can be found in the literature, but the Riemann-Liouville and Caputo fractional derivatives are the most widely utilized ones.

Definition 1. The Riemann-Liouville fractional derivative of order of a function , is defined by

The fractional derivative used in this study is in the Caputo meaning, which is defined as follows.

Definition 2 (see [6]). The Caputo fractional derivative of order of a function , is defined by

3. Derivation Method

3.1. Temporal Discretization

Let denote the integration time ; the time-fractional derivative is approximated by the -approximation [30, 33–35], which is valid for . Explicitly, the time Caputo derivative of order is replaced by the -approximation at which is given in the following lemma.

Lemma 3 (see [36]). Suppose and , it holds that where .

Lemma 4 (see [36]). Let and ; then, , as .

Following Lemma 3 and some algebraic simplifications, we can approximate the time Caputo derivative at as follows:

3.2. Quintic B-Spline Method

This section introduces the quintic B-spline collocation method. We begin by dividing the domain into subinterval where with uniform step for . The quintic B-spline functions for are described by the following relationships [32, 37–39]:

The values of and its first three derivatives are given in Table 1.

Table 1 

Quintic B-splines and their corresponding derivatives.

Let and be the quintic B-spline approximations of the exact solutions and , respectively, of the system considered in eqs. (1) and (2). Since the set of quintic B-splines forms a basis over the domain , the approximate solutions and can be written as where and are the time-dependent unknown quantities to be computed and are the quintic B-spline basis function as shown in Table 1.

Let and be the approximate solutions of and , respectively; then, , , and their first three-order derivatives are determined at the th time level and the nodal points in terms of and as

3.3. Redefined Quintic B-Spline Method

To obtain an approximate solution to systems (1) and (2), we have redefined the quintic B-spline basis functions into a new set of basis functions that vanish on the boundary points since , , , , , , , , , and are nonzero at one of the boundary points. The basis functions are redefined as follows. Allowing the approximate solutions and given by Eqs. (11) and (12) to satisfy the boundary conditions (Eq. (4)) and eliminating , , , , , , , and from the resultant equations, we obtain the approximate solutions for and as

where the weight functions and are given by

and the basis functions as

Substituting the values of from Table 1 into Eqs. (16)–(18), we get

and as

3.4. Description of Numerical Method

In order to apply the suggested method using the redefined set of quintic B-splines basis functions to systems (1)–(4), we write the system by approximate solutions and at time level and nodal points as follows: where and .

The Caputo fractional derivatives are discretized using the -approximation as described in Eqs. (8) and (9), and the nonlinear terms , , and are linearized using the linearization form given by Rubin and Graves [40], ; we get where . After simplifying Eqs. (22) and (23), we get

Using Eqs. (14) and (15) in Eqs. (24) and (25), we obtain

This leads to that

where and are the resulting terms from the weight functions

Using Eq. (20) and Table 1 for the coefficients , , , and , we may rewrite systems (27) and (28) as follows:

where and the following coefficients

System (30) to (34) can be written in blocks of five-diagonal matrices as follows: where and , , , and are matrices where

The initial vectors and may be calculated from Eq. (13) and the initial conditions (Eq. (4)) as follows:

Equation (41) creates two systems, each consisting of equations and unknown variables. So, using Eqs. (42) and (43) at and to eliminate each , , , , , , , and from Eq. (41), we get

Systems (44) and (45) can be written in matrix form as follows:

4. Stability Analysis

In this section, the stability of the suggested approach has been investigated through the implementation of the von Neumann technique. In order to execute this technique, the nonlinear terms , , and in Eqs. (1) and (2) have been linearized by regarding and as local constants and , respectively.

By implementing the -approximation method as shown in Eqs. (8) and (9) and subsequently linearizing the nonlinear terms using the linearization form given by Rubin and Graves [40], the substitution of the approximate solutions for , , and their respective derivatives at the knots in the modified equation results in a difference equation with the variables and given as where

Now, we consider the solutions in terms of Fourier series and at a given point , where and are the harmonic amplitude, and are the mode number and element size, respectively, and . Substituting these solutions in Eq. (47) and simplifying the terms, we get

After performing algebraic simplifications on the given terms, we obtain where , , and .

When we take the absolute value on both sides of Eq. (51), we obtain where using Lemma 4. The necessary and sufficient condition for is that . Since and are positive, it follows that and . Hence, and

Using Eq. (53), we get for all .

Similarly, by substitution of the solution in terms of Fourier series into Eq. (48) and simplification, we get where , , and , and then, for all . As a result, the schemes are unconditionally stable.

5. Error Analysis

Theorem 5 (see [41, 42]). Assume that the exact solutions , , and are an equidistant partition, each of length , over the interval such that . Let and be the unique spline approximations to the given problem at the spatial grid points ; then, for all , there exist , independent of , such that

Theorem 6. Let and be the numerical approximations obtained by the redefined quintic B-spline method to the analytical exact solutions and , respectively, for Eqs. (1)–(4). If , then for sufficiently small and , we have

Proof. Let be the calculated spline for the approximate solutions and the exact solution , respectively. Let , be the collocating conditions. Then, .
Using the difference system (32), the th time step of can be written as where , and
It is evident from Eq. (55) that We define and
Now, we can write Eqs. (57) to (61) in the matrix form as where , , , and The submatrices , , , and are defined by Eqs. (39) and (40).
By defining we can rewrite , , and as follows: where , , and are matrices as follows: The pentadiagonal matrices and are invertible and hold the following condition: According to [43], matrix is invertible. Moreover, From Eq. (63) and norm inequalities, we have From the classifications of the matrices , , and defined in Eq. (66) and the truncation error of time-fractional discretization shown in Eqs. (8) and (9) and the fact that , we have Based on the previous analysis, we deduce that as . Moreover, the convergence rate of the proposed method is of third order.

6. Numerical Results

This section provides some illustrated cases to demonstrate the applicability and efficiency of the proposed technique. All the computations associated with the experiments discussed above were performed in Wolfram Mathematica 13.2 on a PC with Windows 11 64-bit OS + processor Intel Core i7~2.4 GHz.

In order to calculate the maximum absolute error and , we use the following formula:

Accordingly [44], the convergence order (CO) of the proposed approach is given by where and are errors that correspond to grids with mesh size and , respectively, and and .

Example 1. Considering the system of the Whitham-Broer-Kaup equations (Eqs. (1)–(4)) with and , we obtain the suggested equation in [13] taking into account that is only fractional for the temporal variable and for the spatial variable follows

According to [13], the exact solutions for this system are where , , and are arbitrary constants. The initial and boundary conditions (Eqs. (3) and (4)) can be obtained from the exact solution.

In this example, we calculate the maximum absolute error between the approximate solution obtained by the proposed method and two exact solutions and computed by Aminikhah et al. [45] by setting , , and different values of as shown in Tables 2 and 3. Table 2 shows the maximum absolute error between the first exact solutions and the solutions obtained by the proposed method where , , , , and . Moreover, the exact and approximate solutions at the same previous values in addition to are represented for both and of Example 1 in Figure 1. The maximum absolute error between the second exact solutions and the approximate solutions for two sets of parameters and are tabulated in Tables 3 and 4, respectively. Table 5 compares the accuracy of the proposed method and the other popular existing methods for Example 1 at .

Example 2. Consider systems (1)–(4) with , , and

Table 2 

Maximum absolute error for Example 1 at , , , , , , and .

Table 3 

Maximum absolute error for Example 1 at , , , , , and .

(a)

(b)

(c)

(d)

Figure 1 

Comparison between (a, c) the exact solutions and (b, d) the approximate solutions of and for Example 1 at , , , , , , , , and .

Table 4 

Maximum absolute error for Example 1 at , , , , , and .

Table 5 

Maximum absolute error for Example 1 at , , , , and .

The exact solution to this problem is as follows:

The initial and boundary conditions (Eqs. (3) and (4)) can be obtained from the exact solution. In this example, we applied the proposed method for solving the nonhomogeneous systems (1)–(4) for two different sets and . The calculated solutions using the suggested approach are compared to the exact solutions at , , , , and , and the maximum absolute errors and convergence rates are tabulated in Table 6. Figure 2 shows the exact and approximate solutions of and of Example 2 at , , , , , , and .

Table 6 

Maximum absolute error for Example 2 at , , , , and .

(a)

(b)

(c)

(d)

Figure 2 

Comparison between (a, c) the exact solutions and (b, d) the approximate solutions of and for Example 2 at , , , , , , and .

7. Conclusion

In this study, the time-fractional WBK equations were successfully solved using the redefined quintic B-spline collocation method. To achieve this, the -approximation technique in time and the redefined quintic B-spline collocation scheme in space have been combined. We conducted a von Neumann stability analysis, which confirmed that the method used for solving the time-fractional WBK equations is unconditionally stable. The order of convergence is shown to be . To evaluate the accuracy of our approach, we compared our solutions to the exact solutions obtained by Aminikhah et al. [45]. The comparison revealed that our method is highly effective in solving the given equations, as it produced results that closely matched the exact solutions.

Data Availability

No underlying data was collected or produced in this study.

Conflicts of Interest

The authors declare no potential conflict of interest.

Authors’ Contributions

All authors contributed equally. All authors read and approved the final manuscript.