R E S E A R C H Time accurate solution to Benjamin–Bona–Mahony–Burgers equation via Taylor–Boubaker series scheme

The object of this paper is to develop an accurate combined spectral collocation approach to numerically solve the generalized nonlinear Benjamin–Bona–Mahony–Burgers equation. The ﬁrst stage is devoted to discretization in time, which is carried out with the aid of the well-known Taylor series expansions. Then the spectral collocation procedure based on the Boubaker polynomials is applied for the resulting discretized spatial operator in each time step. A detailed error analysis of the presented technique is carried out with regard to the space variable. The advantages of the hybrid technique are shown via performing several simulations through four test examples. Comparisons between our numerical results and the outcomes of some existing schemes indicate that the proposed technique is not only simple and easy-to-implement, but also suﬃciently accurate using a moderate number of bases and a large time step. 65M70; 65M12; 65M06


Introduction
Diverse important physical phenomena in nature are mathematically described by means of (nonlinear) partial differential equations (PDEs) with appropriate initial and boundary conditions. Except for some particular cases, the closed-form solutions to the most PDEs either do not exist or are not intractable in practice. Thus developing an accurate approximate or numerical solution for nonlinear PDEs becomes practically important in the filed. In the past decades the subject has attracted many authors, and difference equations have been appeared as a promising research field on both applied and theoretical levels (see [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]).
The generalized Benjamin-Bona-Mahony-Burgers (BBMB) model (1) was first considered in [1]. Due to the importance and vast applications, many researchers have been considered the BBMB-type equations over the past decades. In this respect, diverse computational and approximation techniques were proposed in the literature. The Galerkin approach was applied to the Benjamin-Bona-Mahony (BBM) equations in [18][19][20]. The finite difference schemes with different combined methods in space and time are developed in [6,7]. The B-spline collocation approaches are studied in [21][22][23]. The meshlessbased methods were investigated in the previously published papers [24,25]. The spectral collocation procedures based upon (orthogonal) basis functions were discussed in [26][27][28][29]. Additionally, a discussion of some other methods for solving PDEs can be found in [30][31][32][33][34][35][36][37][38][39][40][41]. In most of the proposed numerical solution algorithms for (13), the implemented time marching procedure has the first-order accuracy. In addition, to achieve a reasonable accuracy, the time step cannot be taken rather large. Alternatively, in this research work, we construct a second-order time advancement technique for the numerical treatment of the BBMB equation. Our efficient and accurate approach based upon a hybrid of Taylor method for the time discretization and exponential spectral Boubaker collocation scheme for the spatial operator. The presented hybrid technique is straightforward in implementation compared to other existing numerical models; see also [42][43][44].
Recently, the authors in [45] obtained new exact solitary solutions for a version of (3 + 1)dimensional Wazwaz-Benjamin-Bona-Mahony equation formulated in the sense of conformable derivative with the help of two novel techniques: the generalized Kudryashov method and exp(-φ(ℵ)) method. The general fractional formulation of the Wazwaz-Benjamin-Bona-Mahony equation can be expressed as follows: where D ζ is the fractional operator of order 0 < ζ ≤ 1, and the function : [0, ∞) → R is ζ -differentiable at a point t. Matar et al. [46], studied the following fractional differential equation: ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ d dt (ϕ p ( C D α,ρ x(t))) = q(t, x(t), C D γ ,ρ x(t)), 0 ≤ t ≤ 1, x(0) + μx(1) = θ 1 (x(0), x(1)), x (1) = θ 2 (x(0), x(1)), where C D α,ρ and C D γ ,ρ are the generalized Caputo fractional derivatives of orders 0 < α < 1 and 0 < γ < 1, respectively, ρ > 1, and s q : [0, 1] × R 2 → R and θ i : R 2 → R, i = 1, 2, are given nonlinear continuous functions. They investigated possible solutions to the following fractional boundary value problem for an implicit nonlinear an implicit nonlinear fractional differential equation: where C D ς 1 ,ρ and C D ς 2 ,ρ ℘(t) are generalized Caputo fractional derivatives of orders 1 < σ 1 < 2 and 0 < σ 2 < 1, respectively, ρ > 1, ϕ p (p > 1) is a p-Laplacian operator, and w : [0, 1] × R 2 → R and v 1 , v 2 : R 2 → R are given nonlinear continuous functions. Also, the authors in [49] found a theoretical method to investigate the existence of solutions for the strongly singular fractional model of thermostat control given as x (j) = 0, j = 1, 2, . . . , n -1, where (ω -k -1)p(1) > aη k , with j = k coupled, q : [0, 1] → R is singular or strongly singular at some points of [0, 1], p : , and C D ω is the Caputo derivative of order ω. For more studied related application, see [50,51]. Rezapour et al. investigated the multisingular integro-differential pointwise equation , D ς is the Caputo fractional q-derivative of order ς , and w : [0, 1] × R 4 → R is a function such that w(t, ., ., ., .) is singular at some points 0 ≤ t ≤ 1 [52]. In 2021, Izadi [42] presented an effective approximation algorithm to solve the nonlinear Hunter-Saxton equation subjected to the boundary condition lim x→∞ w(x, t) = 0. Abdelwahed and Chorf [53] consider the nonlinear heat equation: Find a solution ϕ of where is a bounded simply connected domain of R d (d = 1, 2, 3), ∂ is its connected Lipschitz continuous boundary, and T is a positive constant. They considered families of large eddy simulation models, which are variants of the classical Smagorinsky model, and similarly to the model of Cottet, Jiroveanu, and Michaux, proposed a selective model based on the local behavior of the angle of the vorticity direction [54]. The chief goal of this research paper is developing an effective time-accurate computational procedure for finding the solutions of a PDE model problem arising in broad branches of science and engineering. We consider a hybrid approximation technique to numerically treat the following BBMB equation: where x ∈ [x L , x R ] and t ∈ [0, T f ]. Also, two coefficient parameters μ, γ , σ are positive constants, η ∈ R, and h(x, t) is a familiar real-valued function. Along with this equation, an initial condition is given as Moreover, with the initial-value problem (1)-(2), the following boundary conditions are supplemented: for t ∈ [0, T f ], where w L (t) and w R (t) are two prescribed functions. If γ = 0, then equation (1) is called the BBM equation, which first studied as a model for the propagation of long waves in nonlinear dispersive systems [2]. Some facts on physical significance of this model were given in [3,5]. For μ = 0 and η = -1, the model (1) reduces to (generalized) Burgers equation [4]. Our plan in the rest of the paper is as follows. In the next Sect. 2, we illustrate the time advancement procedure for BBMB equation (1), which is relied on the Taylor series expansion. Afterward, in Sect. 3, we give an overview of the Boubaker functions together with some important their properties. In this section, we also discuss the convergence for this class of polynomials. The hybrid strategy of Taylor and Boubaker functions is illustrated in Sect. 4 in detail. The results of numerical simulations presented through tables and figures are given in Sect. 5. We end the paper with the conclusion Sect. 6.

Time-marching procedure
To gain an accurate time solution of the nonlinear BBMB equation (1), we first consider the Taylor expansion series to discretize it in time. For this purpose, we uniformly partition the time interval [0, T f ] into J uniform subdivisions. The corresponding grid points and time-step t are given by where t = t jt j-1 , and j = 1, 2, . . . , J. In what follows, by w j we denote the approximate solution at time level t j , that is, We then evaluate the original BBMB equation (1) at time level t j to arrive at Let us for convenience define Application of the Taylor formula to u j yields Differentiating (4) with regard to t reveals Now we replace all first-order derivatives ∂w j ∂t by the forward difference relation After multiplying both sides by t and some manipulations, the resultant equation becomes where h j := h(x, t j ). Our next task is to place (6) into the right-hand side of (5) followed by equating to (4). After multiplying by 2 t, we get the time-discretized equation for the nonlinear model (1) with second-order accuracy in time. Again for simplicity of notations, we introduce Similarly, we define Therefore the resulting second-order equation with regard to space variable can be written as for j = 0, 1, . . . , J -1. In each time level, we must solve the linear equation (7). So we need the given initial condition w 0 = w 0 (x) and its first-and second-order derivatives as they appear in the coefficients of (7). In other words, for j = 0, we have At each time level t j+1 , we exploit the boundary conditions obtained from (3) at two endpoints x = x L , x R as follows: for j = 0, 1, . . . , J -1.

Boubaker functions: an overview
The family of nonorthogonal Boubaker polynomials B (x) naturally appeared in the study of heat transfer equation when its solution can be expanded in terms of Bessel functions of the first kind [55]. These polynomials are very closely related to the Chebyshev polynomials U (x) of the second kind denoted by the relation These polynomials are obtained via the following recursive formula: where the first three Boubaker functions are B 0 (x) = 1, B 1 (x) = x, and For > 0, these functions are the unique solutions of the second-order differential equation where y(x) = B (x). Furthermore, the explicit series solution of Boubaker polynomial B (x) of degree is given by The next theorem establishes the distribution of the zeros of these polynomials on [-2, 2]; see [56,Theorem 2.4].

Theorem 3.1
The function B (x) of degree has has -2 real zeros and two nonreal, purely imaginary complex zeros. All zeros of B (x) are located in [-2, 2], and the two purely imaginary zeros lie outside the unit circle.
Remark 3.1 Based on the observation established in Theorem 3.1, the shifted version of these polynomials were previously considered on an arbitrary domain [x L , x R ] in [57]. This can be accomplished trough the change of variable . Thus we get the shifted version of these polynomial by using B * (η) = B (x). In the numerical examples, we may appropriately use the shifted Boubaker polynomials.

Convergence results
To proceed, we define = [x L , x R ] and the related weighted L 2 space as [57] with ω(x) = 1 x R -x L denotes the induced norm resulting from the following inner product of the space L 2,ω ( ): In practice, we consider a finite-dimensional subspace of L 2,ω ( ) of the form Clearly, dim(L M ) = M + 1, and L M is a closed and thus complete subspace of L 2,ω ( ). Thus any function u ∈ L 2,ω ( ) has a unique best approximation u * ∈ L M in the sense that We further invoke the following result from the approximation theory [58].
Now assume that an arbitrary function u ∈ L 2,ω ( ) can be written in terms of shifted Boubaker functions as By restricting our attention to the finite subsets L M of L 2,ω ( ) we may write a truncated series for u as The next theorem provides an error bound for Let u M be the best approximation to u in the space L M in the sense of (12). Then we have Proof Since u M is the best approximation belonging to L M , we obtain The foregoing inequality is valid in particular for v = P ∈ L M . We therefore conclude that According to Theorem 3.2 with (M + 1) nodes, we get The proof is finished by tending M to infinity.

The hybrid procedure
The discretization of the BBMB equation (1) with regard to time is already carried out by relation (7). Also, the given boundary conditions (3) are converted to the boundary conditions (8) for the second-order differential equation (7). So the main objective is to treat the initial-boundary value problem (7)- (8) numerically with respect to the space variable x. Suppose that the solutions χ j+1 (x) of (7) at each time level j can be written as combinations of the (shifted) Boubaker functions defined by (11). Let assume that X j,M (x) are the computed Boubaker approximations to χ j (x) at the time steps t j . To start, we need χ 0 (x), which can be derived from the given initial condition w 0 (x). In the next time step t j+1 , we look for the approximate solution X j+1,M (x) for j = 0, 1, . . . , J -1 as follows: Here our ultimate goal is to find the unknowns λ j+1 for = 0, 1, . . . , M at time level j ≥ 0.
For convenience, we set With the aid of these vectors, the (M +1)-term finite expansion series (13) can be expressed compactly in the matrix form as After defining the vector of monomials we further decompose the vector of Boubaker bases as where the upper-triangular matrix has the following elements: where the only nonzero elements on the first and second columns are n 0,0 , n 1,1 = 1. For example, for M = 6, we get We note that determinant of N N N M is equal to unity. Ultimately, we requires a set of collocation points on [x L , x R ] to acquire an approximate solution of the discretized model problem (7) in the form (9). In this respect, we use the zeros of the shifted Chebyshev polynomials on [x L , x R ] given by where s τ are roots of Chebyshev functions of degree (M + 1) on (-1, 1). Further, we combine relationships (14) and (15) to express the approximate solution X j+1,M (x) in (13) in a compact form as Using points of collocation (16) and putting them into (17), we obtain The next objective would be finding a connection between V V V M (x) and its derivatives. We can show that the derivatives of V V V M (x) can be stated in terms of differentiation matrix D D D as for k = 1, 2. Twice differentiation of relation (17) with respect to x and using the former relation (19) reveal that Now it is sufficient to insert the points of collocations (16) into the last formulas to obtain the following matrix expressions for the first and second derivatives in (7): j+1,M (x 1 ) . . .
If we substitute the approximate solution X j+1,M (x) and its two derivatives (7), then the resulting equation for j = 0, 1, . . . , J -1 is Inserting the collocation points inserted into the foregoing equation followed by writing it in a matrix representation yields where we have used the coefficient matrices A A A j (B B B j or C C C j ) of size (M + 1) × (M + 1) and the vectors F F F j of size (M + 1) × 1 defined as Based on the matrix forms derived in relations (18)-(22), we get the fundamental matrix equation where To obtain the unknowns λ  . . .
Ultimately, once the above linear system is solved, we get the unknown Boubaker coefficients in (13) or (17).

The given periodic initial and boundary conditions are
The exact analytical solution takes the form We first consider t = 0.1 and T f = 1. We use (13) with M = 8 to obtain the following approximate solutions at the first and last time steps t = t and t = T f for 0 ≤ x ≤ 1: The profile of the approximate solution using these parameters are presented in Fig. 1, left. In Fig. 1, right, we further present the achieved absolute errors E j,M (x) for x ∈ [0, 1], t = 0.1, and T f = 1, evaluated at various time levels t j = j t, j = 1, 2, . . . , 10.
In Tables 1 and 2, some comparisons are done to validate our computational results. For this purpose, the error norms in the L ∞ evaluated at the final time t = T f are calculated. In Table 1 Table 1. We can clearly see that the performance of the proposed hybrid technique is better than the FODM with less computational efforts and the number of bases.
The second-order accuracy of the present technique with respect to time variable is investigated in Table 2. In this respect, we use a fixed M = 10 and various The related outcomes of the existing FODM [7] are also presented in Table 2 for comparison. Obviously, our method with considerably larger time steps t shows its order of accuracy as O( t 2 ) in comparison with FODM. Finally, for this example, we go beyond the unit time interval and consider T f = 3π . In this case the approximate solutions at time t = T f is given by X 94,10 (x) = 0.000214278x 10 -0.00364302x 9 + 0.0137945x 8    The initial and boundary conditions are taken for x ∈ [x L , x R ] = [0, π] and t ∈ [0, T f ] as The exact solution for this model problem is Let us first consider t = 0.1, T f = 1, 2, 4, and T f = 10. The approximated solutions using M = 8 obtained at these t = T f are given as follows:  this purpose, we exploit the quartic B-spline collocation method (QSCM) [21], the nonpolynomial spline method (NPSM) [22], and the improvised cubic B-spline collocation method (ICSCM) [23].
The results using fixed M = 15 and t = 0.01, and diverse values of T f = 1, 2, 4, 10 are listed in Table 3. Clearly, our approach with less computational efforts is more accurate than the QSCM/ICSCM and is comparable with the NPSM.
We next fix M = 12 and vary t as 1/2, 1/4, 1/8, 1/16. The estimated order of accuracy in time (ROCt ∞ ) is investigated at time T f = 10 for this Test Problem as well. The results of achieved E ∞ error norms, together with related ROCt ∞ , are shown in Table 4, where a further comparison with some other available computational methods is done.
We have used the combined method based on a finite difference procedure and a new class of polynomials (FDCP) [28], the FDM [6], and the method based on Lucas polynomials [27]. Obviously, the presented results confirm that the order of accuracy of two is obtainable for the used time-marching algorithm and with relatively large time steps in comparison with other approaches. Page 18 of 29

Test Problem 5.3
The third test problem is devoted to the following nonhomogeneous model with nonunity coefficients [27,28], μ = 0.1, γ = 1, σ = 0.01, η = 1: h(x, t) = 9 5 + 189 100 The computational domain is taken as and the initial and boundary conditions are w(x R , t) = 0.  The exact solution of this example is We first set M = 8 and T f = 1 and use t = 0.1 in the computations. The approximate solution at t = T f on [x L , x R ] is obtained as  t FDCP [28] (N = 10) FDM [6] L u c a s [ 27] τ  are shown in Fig. 6, left. The related absolute errors are further shown in Fig. 6, right. In the next experiments, we indicate the second-order convergence of our approach with respect to time discretization. In Table 5, using M = 10 and various we present the achieved E ∞ and calculate the corresponding ROCt ∞ . Similar results reported by other the existing schemes but with larger number of resources are also presented in Table 5. These approaches are the FDCP [28] with N = 18, the Chebyshev-Legendre method (CLM) [26] with N = 36, and the method based on Lucas polynomials [27] with N = 18. Clearly, the TBM using a few bases, together with relatively large time steps, confirms the order of convergence two. The spectral accuracy with regard to the achieved E 2 /E ∞ error norms for the space variable is also investigated in Fig. 7. In these results, we have used t = 10 -3 and diverse M = 2 i , i = 1, 2, 3, 4.

Test Problem 5.4
In the last test problem, we pay attention to [24,25] with μ = 1, γ = 1, σ = 1, and η = ±1. For η = ±1, we have the right-hand side functions h(x, t), respectively,   As we can see, both approximate solutions are very close to each other. We thus only plot the the graphs of X -10,7 (x) and the related absolute errors for η = -1 in Fig. 9. Next, we examine the behavior of E ∞ error norms when M = 8 is fixed and t varies as 1/2, 1/4, 1/8, 1/16. In addition, the final time is taken as T f = 1, and [x L , x R ] = [0, 1]. These results for η = -1 are reported in Table 6. For comparison, the outcomes of the previously existing computational procedures are also displayed in Table 5. To this end, we use the Liegroup approach based on radial basis functions (LG-RBFs) [25], the meshless method [24], and the Legendre spectral element method (LSEM) [29]. We can observe from the results shown in Table 6 that our numerical results provide second-order accuracy in time while employing a smaller number of bases and time steps compared to other well-established numerical procedures in the literature.
Similar results for η = +1 are further reported in Table 7. However, here we have used T f = 0.1, and the values of t are 1/10, 1/20, 1/40 as used in the results presented in the ICSCM [23] and meshless methods [24] Table 6. The aforementioned results are shown in Table 8, which indicate the superiority of our procedure compared to the LG-RBFs, LSEM, and meshless approaches. Besides, the second-order accuracy in time is visible from the given results if Table 8. For completeness,  For the simulations, we employ both η = ±1 and compute the results at T f = 1 and the spatial domain [-1, 1]. The results are shown in Table 9.

Conclusions
We have developed an accurate time numerical solution algorithm (using a large time step) for the generalized BBMB-type equations (1) arising in diverse disciplines of engineering science. The proposed technique is constructed based on a combination of the Boubaker collocation procedure for the spatial variable and the Taylor series formula for the temporal discretization. The main characteristic of the presented work is that we need to solve an algebraic system of equations at each time step rather than solving a global system obtained in the spectral collocation methods developed in the past. The convergence analysis of the hybrid technique is discussed. The numerical results shown in tables and  figures justify the second-order accuracy in time and the high-order accuracy in the space of the presented technique in comparison with some existing well-established computational schemes.