Numerical approximation for solving time-fractional Benjamin-Bona-Mahony-Burger model via cubic B-spline functions

1.   Introduction

It is well-established that fractional derivatives offer greater flexibility and accuracy compared to integer-order derivatives, particularly when modeling nonclassical engineering and scientific phenomena. Examples include fractional-order recurrent neural networks and various dynamical topics, as discussed in ^[1], a new extension of the fractality concept introduced in financial mathematics by Laskin et al. ^[2], frequency dependance of cell rheological behavior studied in ^[3], and the HIV/AIDS model and the present scenario of COVID-19 discussed in ^[4,5]. Fractional derivatives enable accurate modeling of systems requiring accurate modeling of damping. Numerous examples demonstrate the superiority of fractional calculus over integer-order calculus in such applications. Examples illustrating the application of fractional derivatives include: modeling non-linear earthquake fluctuations and developing fluid dynamic traffic models that mitigate problems caused by high traffic flow ^[6,7]. Numerous mathematical equations are increasingly being formulated as fractional partial differential equations, finding applications in diverse fields such as mathematical biology, transport, finance, visco-elasticity, particle chemistry, and population dynamics ^[8,9,10]. Wave breaking models have two most important equations, which are Kortewegde Varies (KdV) and the Banjamin-Bona-Mahony-Burger (BBM-Burger) equation ^[11]. In mathematics, KdV is a model of wave and shallow water surface. These equations are not applicable to certain long-wave physical systems, which led to the introduction of the BBM-Burger equations. The BBM-Burger equations serve as a refinement of the KdV equations. Specifically, the integer-order nonlinear BBM-Burger equation, which accounts for both dissipative and dispersive effects, is presented in ^[12].

The BBM-Burger equations have numerous applications across various fields, including cracked rock, acoustic-gravity waves in fluids, acoustic waves in a harmonic crystal, and thermodynamics ^[13]. Several methods have been proposed to solve fractional partial-integro differential equations and fractional differential equations, including finite element method with cubic B-spline (CBS) functions for numerical solutions of Burger and Fisher equations ^[14], CBS collocation method for numerical solutions of nonlinear inhomogeneous time-fractional (TF) Burger-Huxley equations ^[15], CBS functions for solving TF diffusion equations involving Caputo-Fabrizio fractional derivatives ^[16], and Atangana-Baleanu fractional derivative (ABFD) with CBS functions for solving TF Burger's equations ^[17]. Many fractional equations have been solved by using the quartic B-spline collocation method ^[18], an extended CBS technique ^[19], shifted Legendre polynomials involved in orthogonal basis function method ^[20,21], homotopy analysis method ^[22], and the fifth-kind Chebyshev polynomial collocation method ^[23]. For discussing dynamic physical structures, a nonhomogeneous fractional BBM-Burger model with a nonlocal viscous term is proposed ^[24],

with initial condition (IC) and boundary conditions (BCs)

where $\tau$, ${\phi}$, $\rho$ and $\chi$ are positive parameters. Due to the complexity of the nonhomogeneous BBM-Burger equation, finding an analytical solution is often challenging. Consequently, recent studies have focused on developing numerical methods to solve this equation, exploring alternative approaches to approximate its solution. Linearized difference schemes for the BBM-Burger equation incorporating a fractional nonlocal viscous term have been investigated in ^[25]. The homotopy analysis method was employed to investigate the BBM-Burger equation by Fakhari ^[26]. Song et al. ^[27] obtained an approximate solution for the fractional BBM-Burgers equation using the homotopy analysis method. Kumar et al. ^[28] solved non-integer order BBM-Burger equations using the novel homotopy analysis transform method. Shakeel et al. ^[29] obtained the exact solution of fractional BBM-Burger equation by using the $(G'/G)$ expansion method. The generalized Atangana-Baleanu BBM-Burgers equation involving dissipative term has been solved by using the modified sub-equation method and new $G'/(b G'+G+a)$ expansion schemes ^[30]. The residual power series technique for the series solution of the TF BBM-Burger equation has been used by Zhang et al. ^[31]. Salih et al. ^[32] applied cubic trigonometric B-spline functions to solve the BBM-Burger equation. A hybrid numerical technique has been used to solve the BBM-Burger equation in ^[33]. Karakoc et al. ^[34] obtained the exact traveling wave solution and computational solutions of the BBM-Burger equation by using the modified Kudryashov method and septic B-spline finite element method, respectively. Omrani ^[35] employed the Crank-Nicolson finite difference method to obtain a numerical solution for the BBM-Burger equation. Karakoc et al. ^[36] applied the CBS finite element method, in conjunction with a lumped Galerkin scheme, to spatially approximate the solitary-wave solution of the nonlinear BBM-Burger equation. Majeed et al. ^[37] employed the Caputo fractional derivative and CBS functions to solve the nonhomogeneous TF BBM-Burger equation. Lin et al. suggested an iterative RBF-based approach in order to solve the nonlinear TF BBM-Burger equation in any domain ^[38]. Mohsin Kamran et al. solved the BBM-Burger equation numerically using Caputo derivative and B-spline basis functions in ^[39]. Hamad Salih ^[40] solved the one-dimensional nonlinear BBM-Burger problem using the novel quartic trigonometric B-spline method based on finite difference. Lalit Mohan and Amit prakash provided a very effective method for analyzing the diffusion wave equation and the fractional BBM-Burger equation in ^[41]. Vineesh Kumar used the ansatz method and Adomian decomposition method on the BBM-Burger equation in ^[42]. Atallah El-shenawy et al. ^[43] uses the collocation approach, which is based on the cubic trigonometric B-Spline methodology, to numerically investigate Troesch's problem. Atallah El-shenawy et al. ^[44] developed a numerical method based on B-spline to solve the time-dependent Emden-Fow-ler-type equations.

The given study is inspired by recent developments in the investigation of computational solution of the nonhomogeneous time-fractional BBM-Burger equation. The objective of this study is to apply the $\theta$-weighted scheme with CBS functions to the nonhomogeneous BBM-Burger equation for obtaining the numerical solution. In the past, many researchers have used B-spline techniques to solve non-homogeneous BBM-Burger equations. However, no one has ever used the CBS method together with the ABFD on the TF term involved in the BBM-Burger equation. The usage of non-singular kernel operator in B-spline methods is novel. Convergence and stability of the proposed problem are analyzed. By presenting few numerical examples, the efficiency and applicability of the proposed scheme is also analyzed. By contrasting analytical and numerical solutions, one can find that the present approach has provided the more effective results. The present scheme is novel for the approximate solution of the nonhomogeneous time-fractional BBM-Burger equation involving ABFD, and, as far as authors are aware, it has never been employed for this purpose before.

This paper has been presented in this way: Basic definitions utilized in this work are given in Section 2. This section contains an important definition of ABFD, which is the core of our work. Section 3 contains the description of the proposed method. Section 4 gives the answer of "how to get the values of the initial vector". Stability and convergence of this problem are given in Sections 5 and 6, respectively. Section 7 contains numerical examples and their discussions. Concluding remarks are presented in Section 8.

2.   Basic definitions

The fundamental definitions including the ABFD with several properties ^[45] are a part of this section.

Definition 1. The ABFD of order $\psi$ and $v\in H'(0, 1)$ is defined as ^[46]:

where $E_{\psi}(\nu)$ is the Mittag-Leffler function given as:

and the standardization function $J(\psi)$ can be described as:

For $0 < \psi \leq 1$, it becomes

Definition 2. If $\hat{q}\in L^{2} [\mathfrak{a}, \mathfrak{b}]$, Parseval's identity is given as ^[47]:

where $\tilde{q}(\acute{n}) = \int_{\mathfrak{a}}^{\mathfrak{b}} \hat{q}(r) e^{2 \pi i \acute{n} r}$ is the Fourier transform for every integer $\acute{n}$.

2.1.   Basis functions

The spatial domain $[\mathfrak{a}, \mathfrak{b}]$ is partitioned as $\mathfrak{a} = \mathfrak{u}_{0} < \mathfrak{u}_{1} < \cdots < \mathfrak{u}_{\hat{N}} = \mathfrak{b}$, where $\mathfrak{u}_{j} = \mathfrak{u}_{0}+j h, \; j = 0(1)\hat{N}$. Domain is divided into $\hat{N}$ equal subintervals of length $h = \frac{\mathfrak{b}-\mathfrak{a}}{\hat{N}}$. Now, assume that $\Psi(\mathfrak{u}, s)$ is the CBS approximation for $v(\mathfrak{u}, s)$ s.t.

where $F_{j}(\mathfrak{u})$ are CBS functions and $y_{j}^{\hat{m}}$ are control points that will be computed at each time interval. The CBS functions are defined as ^[48]:

Domain for time $[0, S]$ can break using knots $0 = s_{0} < s_{1} < \cdots < s_{\hat{M}} = S$ in $\hat{M}$ uniform subintervals $[s_{\hat{m}}, s_{\hat{m}+1}]: s_{\hat{m}} = \hat{m}\Delta s, \; \hat{m} = 0, 1, \cdots, \hat{M}$, where $\Delta s = S/\hat{M}$ implements numerical scheme. The values of $\Psi$, $\Psi_{\mathfrak{u}}$, and $\Psi_{\mathfrak{u} \mathfrak{u}}$ at nodal points can be expressed in terms of the parameter $y_{j}$ by combining Eq (2.2) with expression $F_{j}(\mathfrak{u})$. These values are summarized in the following Table 1:

3.   Description of the method

The ABFD used in (1.1) is discretized at $s = s_{\hat{m}+1}$ as ^[49]:

where $l_{\zeta} = (\zeta+1)E_{\zeta+1}-\zeta E_{\zeta}$ and $E_{\zeta} = E_{\psi, 2}[-\frac{\psi}{1-\psi}(\zeta \Delta s)^{\psi}]$. It is easy to see that

● $l_{\zeta} > 0$ and $l_{0} = E_{1}$, $\zeta = 1:1: \hat{m}$,

● $l_{0} > l_{1} > l_{2} > \cdots > l_{\zeta}$, $l_{\zeta}\rightarrow 0$ as $\zeta \rightarrow \infty$,

● $\sum\limits^{\hat{m}}_{\zeta = 0}(l_{\zeta}-l_{\zeta+1})+l_{\hat{m}+1} = (E_{1}-l_{1})+\sum\limits^{\hat{m}-1}_{\zeta = 1}(l_{\zeta}-l_{\zeta+1})+ l_{\hat{m}} = E_{1}$.

Moreover, truncation error $\lambda^{\hat{m}+1}_{\Delta s}$ is given by ^[49]:

$c_{1}$ is constant, and

where $\hat{\vartheta}$ is a constant. The nonlinear term in Eq (1.1) can be linearized by the generalized formula used in ^[50].

where $\varrho$ is a positive integer. Applying the $\theta$-weighted scheme for spatial discretization, forward difference for temporal discretization, and Eq (3.1) for the ABFD to Eq (1.1) yields:

where $\gamma = \frac{J(\psi)}{1-\psi}$ and $v^{\hat{m}} = v(\mathfrak{u}, s_{\hat{m}})$. For $\theta = 0$, this is an explicit scheme. When $\theta = 0.5$, it becomes the Crank-Nicolson scheme, and for $\theta = 1$, it is fully implicit. This problem is solved here for $\theta = 1$, which yields better results. Substituting Eq (3.3) into (3.4) and simplifying yields:

Combining Eqs (2.2) and (3.5) with the data in Table 1 results in the following system of equations:

where,

The summation term from the right side of Eq (3.6) is dropped for $\zeta = 0$, and this equation is considered for the iteration process for $\zeta \geq 1$. Here, a system of $\hat{N}+1$ linear equations can be obtained for $\hat{N}+3$ unknown parameters $(y_{-1}, y_{0}, y_{1}, y_{2}, y_{3}, \cdots, y_{\hat{N}+1})^{T}$. To get $\hat{N}+3$ equations, boundary conditions will be used. Boundary conditions are

Equation (3.6) in matrix form can be written as:

where

Insert the vector $Y^{\hat{m}+1}$ in (2.2), and the numerical solution at $(\hat{m}+1)^{th}$ time stage can be obtained for $\hat{m} = 0, 1, \cdots, \hat{M}$.

4.   Initial vector

By using the IC, the initial vector $Y^{0} = (y^{0}_{-1}, y^{0}_{0}, y^{0}_{1}, ..., y^{0}_{\hat{N}-1}, y^{0}_{\hat{N}}, y^{0}_{\hat{N}+1})^{T}$ of the given problem can be calculated. Derivatives are approximated at boundary points and the IC is used to know the value of $Y^{0}$ as:

This system gives a matrix of dimension $(\hat{N}+3) \times (\hat{N}+3)$,

where,

5.   The stability analysis

The von Neumann stability analysis is based on the decomposition of numerical errors of numerical approximations into Fourier series ^[51]. This is employed in this section to examine the stability of the proposed numerical scheme. For this, suppose $y^{\hat{m}}_{j}$ symbolizes the growth factor in Fourier mode as:

where $\mu = \mathfrak{p} h$, $\acute{A}$ is the harmonic amplitude, and $\mathfrak{p}$ is the mode number. Now using Eq (5.1) in Eq (3.6), the result is

After further simplification, Eq (5.2) becomes

Using the values of $\kappa_{0}, \; \kappa_{1}, \; \kappa_{2}, \; \varpi_{1}, \; \varpi_{2}$, and $\varpi_{3}$ in Eq (5.3), we have

where

Denominator value is greater than the numerator value. Thus, $| \xi|^2 \leq 1$.

Since the modulus of the eigenvalues must be less than one, the suggested approach for the time-fractional BBM-Burger equation is, therefore, unconditionally stable from (5.4). This indicates that the grid size $h$ and step size $\Delta s$ in the time level are not limited; rather, we should favor the values of $h$ and $\Delta s$ that yield the highest scheme accuracy.

6.   The convergence analysis

The convergence of the proposed technique is examined using the approach presented in ^[52]. We first introduce Theorem 1 and Lemma 6.1, which are based on the work of Hall ^[53] and Boor ^[54].

Theorem 1. Assume that $g(\mathfrak{u}, s)$ and $v(\mathfrak{u}, s)$ belong to $C^{2}(\mathfrak{a}, \mathfrak{b})$ and $C^{4}(\mathfrak{a}, \mathfrak{b})$ respectively. The equidistance partition of $[\mathfrak{a}, \mathfrak{b}]$ is $\Upsilon = [\mathfrak{a} = \mathfrak{u}_{0}, \mathfrak{u}_{1}, \mathfrak{u}_{2}, \cdots, \mathfrak{u}_{\hat{N}} = \mathfrak{b}]$ with stepsize $h$. If $\hat{\Psi}$ is the unique spline interpolation of the given problem at knots $\mathfrak{u}_{0}, \mathfrak{u}_{1}, \mathfrak{u}_{2}, \cdots, \mathfrak{u}_{\hat{N}}\in \Upsilon$, then there exists a constant $\Im_{j}$ independent of $h$ in which $\mathfrak{u}_{j} = \mathfrak{a}+j h$, $j = 0, 1, 2, \cdots, \hat{N}$, then for every $s\geq 0$, it can be obtained that

Lemma 6.1. The cubic B-spline set ${F_{-1}, F_{0}, F_{1}, \cdots, F_{\hat{N}+1}}$ in Eq (2.3) satisfies the inequality

Theorem 2. For the BBM-Burger Eq (1.1) and BCs (1.2), there is a computational approximation $\Psi(\mathfrak{u}, s)$ to the analytical solution $v(\mathfrak{u}, s)$. Moreover, if $g\in C^{2}[\mathfrak{a}, \mathfrak{b}]$, then

where $h$ is relatively small and $\tilde{\Im} > 0$ is free of $h$.

Proof. Assume $\Psi(\mathfrak{u}, s)$ is approximated as $\hat{\Psi}(\mathfrak{u}, s) = \sum\limits^{\hat{N}+1}_{j = -1} w_{j}^{\hat{m}}(t)F_{j}(\mathfrak{u})$. From triangular inequality:

By using Theorem 1, Eq (6.2) becomes:

The proposed scheme has collocation conditions as: $L v(\mathfrak{u}_{j}, s) = L \Psi(\mathfrak{u}_{j}, s) = g(\mathfrak{u}_{j}, s), \; \; j = 0(1) \hat{N}$. Suppose that $L \hat{\Psi}(\mathfrak{u_{j}}, s) = \hat{g}(\mathfrak{u}_{j}, s)$, $j = 0(1)\hat{N}.$ Thus, for any temporal level, the difference $\hat{\Psi}(\mathfrak{u}_{j}, s)-\Psi(\mathfrak{u}_{j}, s)$ with $\theta = 1$ in linear form can be given as:

where

BCs can be written as:

where

and

From inequality (6.1), it is clear that

Define $\beta^{\hat{m}} = max \{ { \mid \beta^{\hat{m}}_{j} \mid; 0\leq j \leq \hat{N} } \}, \; e^{\hat{m}}_{j} = \mid \Omega^{\hat{m}}_{j}\mid$ and $e^{\hat{m}} = max \{\mid e^{\hat{m}}_{j} \mid; 0\leq j \leq \hat{N}\}$. When $\hat{m} = 0$, Eq (6.4) becomes

From IC, $e^{0} = 0$. For sufficiently small mesh spacing $h$ and norms of $\beta^{1}_{i}, \; \Omega^{1}_{i}$, Eq (6.5) yields:

From BCs $e^{1}_{-1}$ and $e^{1}_{\hat{N}+1}$, it can be written as:

From above conditions, it is carried out that

where $\Im_{1}$ is not depending on $h$. To prove this theorem, mathematical induction is applied to $\hat{m}$. It is considered that the term $e^{z}_{j}\leq \Im_{z}h^{2}$ is true, when $z = 1, 2, 3, \cdots, \hat{m}$ and $\Im = max \{\Im_{z}: z = 0, 1, 2, \cdots, \hat{m} \}$, then Eq (6.4) becomes:

Apply the norm again on $\Omega^{\hat{m}+1}_{j}$ and $\beta^{\hat{m}+1}_{j}$, and Eq (6.6) gives

Similarly, from BCs $e^{\hat{m +1}}_{-1}$ and $e^{\hat{m}+1}_{\hat{N}+1}$, it can be taken as:

For all $\hat{m}$, we get

Particularly,

Therefore, the inequality (6.7) and Lemma 6.1 gives

From inequalities (6.3) and (6.8), we obtain:

where $\check{\Im} = \Im_{0}h^{2}+\frac{5}{3}\Im$. □

Theorem 3. The BBM-Burger equation is convergent with the initial and BCs.

Proof. Consider the BBM-Burger equation has analytical solution $v(\mathfrak{u}, s)$ and numerical solution $\Psi(\mathfrak{u}, s)$. Therefore, the previous theorem and inequality (3.2) validate that

where $\check{\Im}$ and $\hat{\vartheta}$ are arbitrary constants. Consequently, the present scheme is second order convergent in both the spatial and temporal directions. □

7.   Examples and discussion

This section includes numerical results and their consistency to show how accurate our results are. Error norms $L_{2}(\hat{N})$ and $L_{\infty}(\hat{N})$ are defined as:

The convergence order of the present scheme can be calculated as ^[55]:

Normalization function in all examples is considered as $R(\psi) = 1$ and $\theta = 1$. Mathematica $9$ is used for numerical calculations on an Intel(R)Core(TM) i5-3437U CPU@2.60GHz, 2712Mhz with 16GB RAM, SSD and 64-bit operating system (Windows 11 pro). Processing to compute numerical results takes less than a minute for all computed results given in tables.

Example 7.1. Consider the nonhomogeneous time-fractional BBM-Burger Eq (1.1) with specific parameter values, namely, ${\phi} = \rho = \tau = \chi = 1.$

with IC and BCs

and the source term is $g(\mathfrak{u}, s) = 2 \bigg(\frac{R(\psi)}{1-\psi} \bigg)s^{2} e^{\mathfrak{u}} E_{\psi, 3}\bigg[\frac{-\psi}{1-\psi} s^{\psi} \bigg]+s^{4} e^{2\mathfrak{u}}$.

The exact solution is $v(\mathfrak{u}, s) = s^{2} e^{\mathfrak{u}}$. Tables 2 and 3 compare exact and computational solutions, along with absolute errors, for different choices of $\psi$ and $s = 0.25, 0.75$, $\psi = 0.5$ with $\hat{N} = 80$ and $\hat{N} = 60$, $\Delta s = 0.001$, respectively, demonstrating good agreement between the exact and approximate solutions. Error norms for various $\psi$ values are displayed in Table 4 at different time stages. The comparison of error norms and the convergence order for different choices of $\Delta s$ is represented in Table 5. Table 6 illustrates convergence order and error norm for different values of $h$. A very good agreement between exact and spline solutions at various time levels for fixed $\hat{N} = 40$ with $\psi = 0.3$ and $\Delta s = 0.01$ is shown in Figure 1. This figure also demonstrates good results for $\hat{N} = 80$, $\psi = 0.7$, and $\Delta s = 0.002$. The correctness of the suggested approach is demonstrated in Figure 2 using 3D graphs comparing analytical and numerical answers at $\hat{N} = 60$, $\psi = 0.5$, $\Delta s = 0.001$, and $s = 1$. At $s = 1$, 2D and 3D error profiles are demonstrated in Figure 3. All tables and figures guarantee the exactness of the scheme.

Example 7.2. Consider the nonhomogeneous time-fractional BBM-Burger Eq (1.1) with specific parameter values, namely, ${\phi} = \rho = \tau = \chi = 1.$

with IC and BCs

and the source term is $g(\mathfrak{u}, s) = 6 \bigg(\frac{R(\psi)}{1-\psi} \bigg)s^{3} \sin(\mathfrak{u}) E_{\psi, 4}\bigg[\frac{-\psi}{1-\psi} s^{\psi} \bigg]+3 s^{2} \sin (\mathfrak{u})+s^{3} cos (\mathfrak{u})+3 s^{2} \sin(\mathfrak{u})+(s^{3} \sin(\mathfrak{u}))(s^{3} \cos(\mathfrak{u}))+ s^{3} \sin(\mathfrak{u})$.

The exact solution is $\psi(\mathfrak{u}, s) = s^{3} sin(\mathfrak{u})$. The absolute errors for various values of $\mathfrak{u}$ setting $\Delta s \; = \; 0.001$, $s = 1$, $\psi = 0.5$, and $\hat{N} = 30$ are reported in Table 7 of Example 7.2. Table 8 presents the exact solutions, approximate results, and absolute errors for the proposed problem when $\Delta s \; = \; 0.001$, $s = 0.25$, $\psi = 0.2$, $\hat{N} = 50$. Comparison of error norms is tabulated in Table 9. Tables 10 and 11 exhibit error norms and their convergence orders for various choices of parameters in temporal and spatial directions. Figure 4 shows the significant agreement between the numerical findings and exact solutions of the suggested scheme at different time levels. A 3D representation of exact and computational solutions at distinct values of parameters is given in Figure 5. The 2D and 3D error profile is exhibited in Figure 6, which demonstrates procedure of best accuracy.

8.   Conclusions

An accurate numerical solution of the nonhomogeneous time-fractional BBM-Burger equation based on CBS functions has been determined in this study. For this purpose, CBS functions have been used to build up a collocation technique for the nonhomogeneous time-fractional BBM-Burger equation. The TF derivative has been approximated by the typical finite difference scheme and Atangana-Baleanu fractional operator, whereas the spatial derivative was discretized by using the $\theta$-weighted scheme with CBS functions. Two test problems have been solved, and their graphical and numerical comparison uncovers that the proposed method is computationally very effective. The scheme has possessed the second order temporal and spatial convergence, as well as unconditionally stable. This method can be used to get efficient approximations for a large number of fractional differential equations. It has provided the results in the form of an improved solution for the problems for which there is no precise solution. This method offers effective solutions to numerous problems. In the future, we may consider the solution of higher dimensional and higher order BBM-Burger equations by using spline functions of higher degrees.

Author contributions

Muserat Shaheen: Methodology, writing-original draft; Muhammad Abbas: Supervision, methodology, writing-original draft; Miguel Vivas-Cortez: Software, formal analysis, writing-review & editing; M. R. Alharthi: Visualization, writing-review & editing; Y. S. Hamed: Software, formal analysis, writing-review & editing. All authors have read and approved the final version of the manuscript for publication.

Use of Generative-AI tools declaration

The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgements

The authors would like to acknowledge the Deanship of Graduate Studies and Scientific Research, Taif University for funding this work.

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.