Numerical Solution of Third-Order Rosenau–Hyman and Fornberg–Whitham Equations via B-Spline Interpolation Approach

Abstract

This study aims to find the numerical solution of the Rosenau–Hyman and Fornberg–Whitham equations via the quintic B-spline collocation method. Quintic B-spline, along with finite difference and theta-weighted schemes, is used for the discretization and approximation purposes. The effectiveness and robustness of the procedure is assessed by comparing the computed results with the exact and available results in the literature using absolute and relative error norms. The stability of the proposed scheme is studied using von Neumann stability analysis. Graphical representations are drawn to analyze the behavior of the solution.

1. Introduction

Partial differential equation (PDE) arises in the modeling of various physical phenomenon and is considered the backbone of mathematical simulations of real-life applications. Many models have been introduced to describe different problems in engineering, physics, electrochemistry, computational biology, and fluid mechanics. Examples of such models appear in mathematical biology [1], Keller–Segel models in pattern formation in biological tissues [2], health estimation via the electrochemical model [3], the fluid dynamics model for supply chains [4], population genetics [5], hormone-dependent cancer [6], and to estimate the internal temperature of a battery [7].

Mostly non-linear PDEs arise in the modeling of any physical systems. Closed-form or analytical solutions of such types of models are difficult to find. Perturbation techniques have made significant contributions to tackling such non-linear phenomena; however, these techniques are based on some small parameters and fail when there are strongly non-linear problems. Semi-analytical methods such as the Adomian decomposition method and variational iteration method are valid for such problems but are again computationally expensive. Keeping all these difficulties in mind, we are moving towards numerical techniques.

Numerical techniques are simple, fast, accurate, and easy to analyze. Most of the numerical techniques have been used to solve non-linear PDEs, including the radial basis function method (RBF). It is a meshless method that has been used for a wide range of problems. The main disadvantage is that the selection of the shape parameter in RBF is a difficult task. Another method is the Haar wavelet, which involves integration; sometimes, it is difficult to evaluate that integration and may lead to computational costs while dealing with higher-order DEs. The finite difference method (FDM) is one of the oldest and finest techniques. In FDM, derivatives are approximated at the mesh points through finite differences. Although this method is very fruitful, it fails when solving problems with irregular regions.

The spline-based method has gained significant consideration in the approximation of different physical problems. Initially, splines were used for curve fitting. Among the family of splines, B-splines are a simple and effective method due to their high degree of smoothness. A B-spline is a piecewise polynomial function that gives local compact support, meaning that a slight change in any of values will result in a change in that particular part of the curve. The theory of the B-spline technique is an active field of approximation theory and BVPs. B-spline curves are valuable because of the suitable shape representation and analysis [8,9,10,11]. Several problems have been tackled using the proposed method such as the telegraphic equation [12], Fokker–Planck equation [13], Fisher’s equation [14], regularized long wave (RLW) equation [15], generalized Kuramoto–Sivashinsky equation [16], coupled Burger’s equation [17], generalized Black–Scholes equation [18], incompressible Navier–Stokes equations [19], and Rosenau–KdV–RLW equations [20]. The order approximation of quintic B-spline was briefly discussed and proved error bound, given in [21] and the reference therein.

RH and FW equations are special cases of the Gilson–Pickering (GP) equation, introduced by Gilson and Pickering [22], given by

$${\mathrm{U}}_{\mathrm{t}}-\nu {\mathrm{U}}_{\varrho \varrho \mathrm{t}}+2\kappa {\mathrm{U}}_{\varrho }-\mathrm{U}{\mathrm{U}}_{3\varrho }-\alpha \mathrm{U}{\mathrm{U}}_{\varrho }-\beta {\mathrm{U}}_{\varrho }{\mathrm{U}}_{\varrho \varrho }=0,c\le \varrho \le d.$$

(1)

The classical RH are given by [23]:

$${\mathrm{U}}_{\mathrm{t}}+\alpha {\left({\mathrm{U}}^n\right)}_{\varrho }+\beta {\left({\mathrm{U}}^n\right)}_{3\varrho }=0.$$

(2)

From Equation (1), when $\nu =\kappa =0,\alpha =1\mathrm{and}\beta =3$, we obtain a RH equation of the form:

$${\mathrm{U}}_{\mathrm{t}}-\mathrm{U}{\mathrm{U}}_{3\varrho }-\mathrm{U}{\mathrm{U}}_{\varrho }-3{\mathrm{U}}_{\varrho }{\mathrm{U}}_{\varrho \varrho }=0.$$

(3)

In Equation (2), when $n=4,\alpha =0.2\mathrm{and}\beta ={\displaystyle \frac{1}{2}}$, we obtain a RH equation of the form:

$${\mathrm{U}}_{\mathrm{t}}+0.8{\mathrm{U}}^3{\mathrm{U}}_{\varrho }+2{\mathrm{U}}^3{\mathrm{U}}_{3\varrho }+18{\mathrm{U}}^2{\mathrm{U}}_{\varrho }{\mathrm{U}}_{\varrho \varrho }+12\mathrm{U}{\left({\mathrm{U}}_{\varrho }\right)}^3=0.$$

(4)

Using Equation (1) with $\nu =1,\kappa ={\displaystyle \frac{1}{2}},\alpha =-1\mathrm{and}\beta =3$, we have a FW equation of the form:

$${\mathrm{U}}_{\mathrm{t}}-{\mathrm{U}}_{\varrho \varrho \mathrm{t}}+{\mathrm{U}}_{\varrho }-\mathrm{U}{\mathrm{U}}_{3\varrho }+\mathrm{U}{\mathrm{U}}_{\varrho }-3{\mathrm{U}}_{\varrho }{\mathrm{U}}_{\varrho \varrho }=0,$$

(5)

along with initial and boundary conditions as

$$\begin{array}{cc}\hfill & \mathrm{U}(\varrho ,0)=f\left(\varrho \right),\hfill \\ \hfill & \mathrm{U}(c,\mathrm{t})=g\left(\mathrm{t}\right),{\mathrm{U}}_{\varrho }(c,\mathrm{t})=g_1\left(\mathrm{t}\right),\hfill \\ \hfill & \mathrm{U}(d,\mathrm{t})=g_2\left(\mathrm{t}\right),{\mathrm{U}}_{\varrho }(d,\mathrm{t})=g_3\left(\mathrm{t}\right),\hfill \end{array}$$

(6)

where $f,g,g_1,g_2$ and $g_3$ are differentiable functions. GP equation is the KdV type equation, which is a well-known equation in the field of non-linear waves. Philip Rosenau and James M. Hyman used this equation for the study of compactons and was named after that as RH in 1993. This equation has been utilized for the modeling of non-linear dispersion in droplets of liquid patterns [23]. It is used to calculate the growth rate of particles in a liquid suspension. FW equation was first introduced by Whitham in 1967. These equations have been utilized for analyzing the qualitative characteristics of wave breaking. Because of their use in physics and engineering, many authors have tackled these equations through different numerical techniques; some of the techniques to solve RH equations are ADM [24], the hybrid and reduce differential transform method (RDTM) [25], the Hermite wavelet method [26], VIM and HPM [27]. Similarly, different techniques for FW equations are HAM and HPM [28,29] and modified VIM [30].

The rest of the article is organized in the following way. In Section 2, the definition of quintic B-splines is discussed, followed by the methodology. In Section 3, the stability of the technique is analyzed using von Neumann stability analysis. In Section 4, the effectiveness of the proposed method is studied by comparing the numerical and graphical results in the existing methods in the literature and finally, the article is concluded in Section 5.

2. Definition

In this section, we construct a method for the numerical solution of the third-order RH and FW equations using the quintic B-spline collocation method. We have to find a solution in the domain $[a,b]$, where $a={\varrho }_0<{\varrho }_1<{\varrho }_2\cdots <{\varrho }_N=b$. The partition of the space domain is assumed to be distributed uniformly with $h={\varrho }_i-{\varrho }_{i-1}$, where $i=1,2,\cdots ,N$. The quintic B-spline is given by

$$B_i\left(\varrho \right)={\displaystyle \frac{1}{h^5}}\left\{\begin{array}{cc}{\left(\varrho -{\varrho }_{i-3}\right)}^5,\hfill & \left[{\varrho }_{i-3},{\varrho }_{i-2}\right)\hfill \\ {\left(\varrho -{\varrho }_{i-3}\right)}^5-6{\left(\varrho -{\varrho }_{i-2}\right)}^5,\hfill & \left[{\varrho }_{i-2},{\varrho }_{i-1}\right)\hfill \\ {\left(\varrho -{\varrho }_{i-3}\right)}^5-6{\left(\varrho -{\varrho }_{i-2}\right)}^5+15{\left(\varrho -{\varrho }_{i-1}\right)}^5,\hfill & \left[{\varrho }_{i-1},{\varrho }_i\right)\hfill \\ {\left(\varrho -{\varrho }_{i-3}\right)}^5-6{\left(\varrho -{\varrho }_{i-2}\right)}^5+15{\left(\varrho -{\varrho }_{i-1}\right)}^5-20{\left(\varrho -{\varrho }_i\right)}^5,\hfill & \left[{\varrho }_i,{\varrho }_{i+1}\right)\hfill \\ {\left(\varrho -{\varrho }_{i-3}\right)}^5-6{\left(\varrho -{\varrho }_{i-2}\right)}^5+15{\left(\varrho -{\varrho }_{i-1}\right)}^5-20{\left(\varrho -{\varrho }_i\right)}^5+15{\left(\varrho -{\varrho }_{i+1}\right)}^5,\hfill & \left[{\varrho }_{i+1},{\varrho }_{i+2}\right)\hfill \\ {\left(\varrho -{\varrho }_{i-3}\right)}^5-6{\left(\varrho -{\varrho }_{i-2}\right)}^5+15{\left(\varrho -{\varrho }_{i-1}\right)}^5-20{\left(\varrho -{\varrho }_i\right)}^5+\hfill \\ 15{\left(\varrho -{\varrho }_{i+1}\right)}^5-6{\left(\varrho -{\varrho }_{i+2}\right)}^5,\hfill & \left[{\varrho }_{i+2},{\varrho }_{i+3}\right)\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(7)

The function approximation in terms of the quintic B-spline (QBS) is defined as

$$\mathrm{U}(\varrho ,\mathrm{t})=\sum _{i=-2}^{N+2}{\mathrm{\Lambda }}_i\left(\mathrm{t}\right)B_i\left(\varrho \right),$$

(8)

where ${\mathrm{\Lambda }}_i$ are time-dependent unknowns, yet to be calculated. The QBS has zero contribution outside the interval except for $B_{i-2},B_{i-1},B_i,B_{i+1}$, and $B_{i+2}$. The values of function and their derivatives at grid points are given by

$$\left\{\begin{array}{l}{\displaystyle \mathrm{U}({\varrho }_i,\mathrm{t})=\sum _{j=i-2}^{i+2}{\mathrm{\Lambda }}_j\left(\mathrm{t}\right)B_j\left(\varrho \right)=({\mathrm{\Lambda }}_{i-2}\left(\mathrm{t}\right)+26{\mathrm{\Lambda }}_{i-1}\left(\mathrm{t}\right)+66{\mathrm{\Lambda }}_i\left(\mathrm{t}\right)+26{\mathrm{\Lambda }}_{i+1}\left(\mathrm{t}\right)+{\mathrm{\Lambda }}_{i+2}\left(\mathrm{t}\right)),}\hfill \\ {\displaystyle {\mathrm{U}}_{\varrho }({\varrho }_i,\mathrm{t})=\sum _{j=i-2}^{i+2}{\mathrm{\Lambda }}_j\left(\mathrm{t}\right)B_j^{\prime }\left(\varrho \right)=\frac{5}{h}(-{\mathrm{\Lambda }}_{i-2}\left(\mathrm{t}\right)-10{\mathrm{\Lambda }}_{i-1}\left(\mathrm{t}\right)+10{\mathrm{\Lambda }}_{i+1}\left(\mathrm{t}\right)+{\mathrm{\Lambda }}_{i+2}\left(\mathrm{t}\right)),}\hfill \\ {\displaystyle {\mathrm{U}}_{\varrho \varrho }({\varrho }_i,\mathrm{t})=\sum _{j=i-2}^{i+2}{\mathrm{\Lambda }}_j\left(\mathrm{t}\right)B_j^{\prime \prime }\left(\varrho \right)=\frac{20}{h^2}({\mathrm{\Lambda }}_{i-2}\left(\mathrm{t}\right)+2{\mathrm{\Lambda }}_{i-1}\left(\mathrm{t}\right)-6{\mathrm{\Lambda }}_i\left(\mathrm{t}\right)+2{\mathrm{\Lambda }}_{i+1}\left(\mathrm{t}\right)+{\mathrm{\Lambda }}_{i+2}\left(\mathrm{t}\right)),}\hfill \\ {\displaystyle {\mathrm{U}}_{3\varrho }({\varrho }_i,\mathrm{t})=\sum _{j=i-2}^{i+2}{\mathrm{\Lambda }}_j\left(\mathrm{t}\right)B_j^{\prime \prime \prime }\left(\varrho \right)=\frac{60}{h^3}(-{\mathrm{\Lambda }}_{i-2}\left(\mathrm{t}\right)+2{\mathrm{\Lambda }}_{i-1}\left(\mathrm{t}\right)-2{\mathrm{\Lambda }}_{i+1}\left(\mathrm{t}\right)+{\mathrm{\Lambda }}_{i+2}\left(\mathrm{t}\right)).}\hfill \end{array}\right.$$

(9)

2.1. Description of the Method

Rewriting Equation (1), we have

$${\mathrm{U}}_{\mathrm{t}}-\nu {\mathrm{U}}_{\varrho \varrho \mathrm{t}}+2\kappa {\mathrm{U}}_{\varrho }-\mathrm{U}{\mathrm{U}}_{3\varrho }-\alpha \mathrm{U}{\mathrm{U}}_{\varrho }-\beta {\mathrm{U}}_{\varrho }{\mathrm{U}}_{\varrho \varrho }=0.$$

(10)

The discretization of PDEs is essential for applying numerical methods, as it transforms continuous equations into a system of algebraic equations that can be solved computationally. This approach enables the approximation of solutions for complex PDEs that are otherwise analytically intractable. Therefore, for time discretization, the forward difference method is used and by applying the $\theta$-weighted scheme to Equation (10), we obtain

$$\begin{array}{cc}\hfill & \left[{\displaystyle \frac{{\mathrm{U}}^{n+1}-{\mathrm{U}}^n}{\delta \mathrm{t}}}\right]-\nu {\displaystyle \frac{{\mathrm{U}}_{\varrho \varrho }^{n+1}-{\mathrm{U}}_{\varrho \varrho }^n}{\delta \mathrm{t}}}+2\theta \kappa {\mathrm{U}}_{\varrho }^{n+1}-\theta {\left(\mathrm{U}{\mathrm{U}}_{3\varrho }\right)}^{n+1}-\theta \alpha {\left(\mathrm{U}{\mathrm{U}}_{\varrho }\right)}^{n+1}\hfill \\ \hfill & -\beta \theta {\left({\mathrm{U}}_{\varrho }{\mathrm{U}}_{\varrho \varrho }\right)}^{n+1}=-(1-\theta )\left(2\kappa {\mathrm{U}}_{\varrho }^n-{\left(\mathrm{U}{\mathrm{U}}_{3\varrho }\right)}^n-\alpha {\left(\mathrm{U}{\mathrm{U}}_{\varrho }\right)}^n-\beta {\left({\mathrm{U}}_{\varrho }{\mathrm{U}}_{\varrho \varrho }\right)}^n\right),\hfill \end{array}$$

(11)

where $\delta t$ represents the time step size, n represents the time level, and the time grid points are $t_n$. The terms ${\left(\mathrm{U}{\mathrm{U}}_{3\varrho }\right)}^{n+1}$, ${\left(\mathrm{U}{\mathrm{U}}_{\varrho }\right)}^{n+1}\mathrm{and}{\left({\mathrm{U}}_{\varrho }{\mathrm{U}}_{\varrho \varrho }\right)}^{n+1}$ are lineraized using the quasi-linearization technique:

$$\begin{array}{cc}\hfill & {\left(\mathrm{U}{\mathrm{U}}_{3\varrho }\right)}^{n+1}={\mathrm{U}}^{n+1}{\mathrm{U}}_{3\varrho }^n+{\mathrm{U}}^n{\mathrm{U}}_{3\varrho }^{n+1}-{\left(\mathrm{U}{\mathrm{U}}_{3\varrho }\right)}^n,\hfill \\ \hfill & {\left(\mathrm{U}{\mathrm{U}}_{\varrho }\right)}^{n+1}={\mathrm{U}}^{n+1}{\mathrm{U}}_{\varrho }^n+{\mathrm{U}}^n{\mathrm{U}}_{\varrho }^{n+1}-{\left(\mathrm{U}{\mathrm{U}}_{\varrho }\right)}^n,\hfill \\ \hfill & {\left({\mathrm{U}}_{\varrho }{\mathrm{U}}_{\varrho \varrho }\right)}^{n+1}={\mathrm{U}}_{\varrho }^{n+1}{\mathrm{U}}_{\varrho \varrho }^n+{\mathrm{U}}_{\varrho }^n{\mathrm{U}}_{\varrho \varrho }^{n+1}-{\left({\mathrm{U}}_{\varrho }{\mathrm{U}}_{\varrho \varrho }\right)}^n.\hfill \end{array}$$

(12)

Using Equation (12) in (11), we obtain

$$\begin{array}{cc}\hfill & \left(1-\theta \delta \mathrm{t}{\mathrm{U}}_{3\varrho }^n-\alpha \theta \delta \mathrm{t}{\mathrm{U}}_{\varrho }^n\right){\mathrm{U}}^{n+1}+\left(2\kappa -\alpha {\mathrm{U}}^n-\beta {\mathrm{U}}_{\varrho \varrho }^n\right)\theta \delta \mathrm{t}{\mathrm{U}}_{\varrho }^{n+1}+\left(-\nu -\beta \theta \delta \mathrm{t}{\mathrm{U}}_{\varrho }^n\right){\mathrm{U}}_{\varrho \varrho }^{n+1}\hfill \\ \hfill & -\theta \delta \mathrm{t}{\mathrm{U}}^n{\mathrm{U}}_{3\varrho }^{n+1}={\mathrm{U}}^n-\nu {\mathrm{U}}_{\varrho \varrho }^n-\theta \delta \mathrm{t}\left({\left(\mathrm{U}{\mathrm{U}}_{3\varrho }\right)}^n+\alpha {\left(\mathrm{U}{\mathrm{U}}_{\varrho }\right)}^n+\beta {\left({\mathrm{U}}_{\varrho }{\mathrm{U}}_{\varrho \varrho }\right)}^n\right)-\hfill \\ \hfill & (1-\theta )\delta \mathrm{t}\left(2\kappa {\mathrm{U}}_{\varrho }^n-{\left(\mathrm{U}{\mathrm{U}}_{3\varrho }\right)}^n-\alpha {\left(\mathrm{U}{\mathrm{U}}_{\varrho }\right)}^n-\beta {\left({\mathrm{U}}_{\varrho }{\mathrm{U}}_{\varrho \varrho }\right)}^n\right).\hfill \end{array}$$

(13)

Using Equation (9), we have

$$\begin{array}{cc}\hfill & d_i{\mathrm{\Lambda }}_{i-2}^{n+1}+o_i{\mathrm{\Lambda }}_{i-1}^{n+1}+p_i{\mathrm{\Lambda }}_i^{n+1}+q_i{\mathrm{\Lambda }}_{i+1}^{n+1}+r_i{\mathrm{\Lambda }}_{i+2}^{n+1}\hfill \\ \hfill & =D_i{\mathrm{\Lambda }}_{i-2}^n+O_i{\mathrm{\Lambda }}_{i-1}^n+p_i{\mathrm{\Lambda }}_i^n+P_i{\mathrm{\Lambda }}_{i+1}^n+Q_i{\mathrm{\Lambda }}_{i+2}^n,\hfill \end{array}$$

(14)

where

$$\begin{array}{cc}\hfill & d_i=\left(a_i+{\displaystyle \frac{5\theta \delta \mathrm{t}{b}_i}{h}}+{\displaystyle \frac{20c_i}{h^2}}+{\displaystyle \frac{60e_i}{h^3}}\right),o_i=\left(26a_i+{\displaystyle \frac{50\theta \delta \mathrm{t}{b}_i}{h}}+{\displaystyle \frac{40c_i}{h^2}}-{\displaystyle \frac{120e_i}{h^3}}\right),p_i=\left(66a_i-{\displaystyle \frac{120c_i}{h^2}}\right),\hfill \\ \hfill & q_i=\left(26a_i-{\displaystyle \frac{50\theta \delta \mathrm{t}{b}_i}{h}}+{\displaystyle \frac{40c_i}{h^2}}+{\displaystyle \frac{120e_i}{h^3}}\right),r_i=\left(a_i-{\displaystyle \frac{5\theta \delta \mathrm{t}{b}_i}{h}}+{\displaystyle \frac{20c_i}{h^2}}-{\displaystyle \frac{60e_i}{h^3}}\right),\hfill \\ \hfill & D_i=\left(a_i-{\displaystyle \frac{5(1-\theta )\delta \mathrm{t}{b}_i}{h}}+{\displaystyle \frac{20c_i}{h^2}}-{\displaystyle \frac{60f_i}{h^3}}\right),O_i=\left(26a_i-{\displaystyle \frac{5(1-0)\theta \delta \mathrm{t}{b}_i}{h}}+{\displaystyle \frac{40c_i}{h^2}}+{\displaystyle \frac{120f_i}{h^3}}\right),\hfill \\ \hfill & P_i=\left(26a_i+{\displaystyle \frac{50(1-\theta )\delta \mathrm{t}{b}_i}{h}}+{\displaystyle \frac{40c_i}{h^2}}-{\displaystyle \frac{120f_i}{h^3}}\right),Q_i=\left(a_i+{\displaystyle \frac{5(1-\theta )\delta \mathrm{t}{b}_i}{h}}+{\displaystyle \frac{20c_i}{h^2}}+{\displaystyle \frac{60f_i}{h^3}}\right),\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill & a_i=\left(1-\theta \delta \mathrm{t}\left({\displaystyle \frac{60}{h^3}}(-{\mathrm{\Lambda }}_{i-2}^n+2{\mathrm{\Lambda }}_{i-1}^n-2{\mathrm{\Lambda }}_{i+1}^n+{\mathrm{\Lambda }}_{i+2}^n)\right)-\alpha \theta \delta \mathrm{t}\left({\displaystyle \frac{5}{h}}(-{\mathrm{\Lambda }}_{i-2}^n-10{\mathrm{\Lambda }}_{i-1}^n+10{\mathrm{\Lambda }}_{i+1}^n+{\mathrm{\Lambda }}_{i+2}^n)\right)\right),\hfill \\ \hfill & b_i=\left(2\kappa -\alpha ({\mathrm{\Lambda }}_{i-2}^n+26{\mathrm{\Lambda }}_{i-1}^n+66{\mathrm{\Lambda }}_i^n+26{\mathrm{\Lambda }}_{i+1}^n+{\mathrm{\Lambda }}_{i+2}^n)-\beta \left({\displaystyle \frac{20}{h^2}}({\mathrm{\Lambda }}_{i-2}^n+2{\mathrm{\Lambda }}_{i-1}^n-6{\mathrm{\Lambda }}_i^n+2{\mathrm{\Lambda }}_{i+1}^n+{\mathrm{\Lambda }}_{i+2}^n)\right)\right),\hfill \\ \hfill & c_i=\left(-\nu -\beta \theta \delta \mathrm{t}\left({\displaystyle \frac{5}{h}}(-{\mathrm{\Lambda }}_{i-2}^n-10{\mathrm{\Lambda }}_{i-1}^n+10{\mathrm{\Lambda }}_{i+1}^n+{\mathrm{\Lambda }}_{i+2}^n)\right)\right),e_i=\theta \delta \mathrm{t}({\mathrm{\Lambda }}_{i-2}^n+26{\mathrm{\Lambda }}_{i-1}^n+66{\mathrm{\Lambda }}_i^n+26{\mathrm{\Lambda }}_{i+1}^n+{\mathrm{\Lambda }}_{i+2}^n),\hfill \\ \hfill & f_i=(1-\theta )\delta \mathrm{t}({\mathrm{\Lambda }}_{i-2}^n+26{\mathrm{\Lambda }}_{i-1}^n+66{\mathrm{\Lambda }}_i^n+26{\mathrm{\Lambda }}_{i+1}^n+{\mathrm{\Lambda }}_{i+2}^n).\hfill \end{array}$$

The corresponding BCs are evaluated at the grid points $t_n$, which are given by

$$\begin{array}{cc}\hfill & {\mathrm{\Lambda }}_{-2}^n+26{\mathrm{\Lambda }}_{-1}^n+66{\mathrm{\Lambda }}_0^n+26{\mathrm{\Lambda }}_1^n+{\mathrm{\Lambda }}_2^n=g\left({\mathrm{t}}_{\mathrm{n}}\right),\hfill \\ \hfill & {\displaystyle \frac{5}{h}}(-{\mathrm{\Lambda }}_{-2}^n-10{\mathrm{\Lambda }}_{-1}^n+10{\mathrm{\Lambda }}_1^n+{\mathrm{\Lambda }}_2^n)=g_1\left({\mathrm{t}}_{\mathrm{n}}\right),\hfill \\ \hfill & {\mathrm{\Lambda }}_{N-2}^n+26{\mathrm{\Lambda }}_{N-1}^n+66{\mathrm{\Lambda }}_N^n+26{\mathrm{\Lambda }}_{N+1}^n+{\mathrm{\Lambda }}_{N+2}^n=g_2\left({\mathrm{t}}_{\mathrm{n}}\right),\hfill \\ \hfill & {\displaystyle \frac{5}{h}}(-{\mathrm{\Lambda }}_{N-2}^n-10{\mathrm{\Lambda }}_{N-1}^n+10{\mathrm{\Lambda }}_{N+1}^n+{\mathrm{\Lambda }}_{N+2}^n)=g_3\left({\mathrm{t}}_{\mathrm{n}}\right).\hfill \end{array}$$

Using these BCs in (14), we obtain

$$S_1{\mathrm{\Lambda }}^{n+1}=S_2{\mathrm{\Lambda }}^n,$$

(15)

where $S_1$ and $S_2$ are square matrices of dimensions $(N+5)\times (N+5)$, given by

$$S_1=\left[\begin{array}{cccccccccc}1& 26& 66& 26& 1& \cdots & \cdots & \cdots & \cdots & 0\\ {\displaystyle \frac{-5}{h}}& {\displaystyle \frac{-50}{h}}& 0& {\displaystyle \frac{50}{h}}& {\displaystyle \frac{5}{h}}& \ddots & & & & ⋮\\ d_0& o_0& p_0& q_0& r_0& \ddots & \ddots & & & ⋮\\ 0& d_1& o_1& p_1& q_1& r_1& \ddots & \ddots & & ⋮\\ ⋮& \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ ⋮& & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ ⋮& & & \ddots & d_{N-1}& o_{N-1}& p_{N-1}& q_{N-1}& r_{N-1}& 0\\ ⋮& & & & \ddots & d_N& o_N& p_N& q_N& r_N\\ 0& \cdots & \cdots & \cdots & \cdots & {\displaystyle \frac{-5}{h}}& {\displaystyle \frac{-50}{h}}& 0& {\displaystyle \frac{50}{h}}& {\displaystyle \frac{5}{h}}\\ 0& \cdots & \cdots & \cdots & \cdots & 1& 26& 66& 26& 1\end{array}\right],$$

(16)

$${\mathrm{\Lambda }}^{n+1}={({\mathrm{\Lambda }}_{-2}^{n+1},{\mathrm{\Lambda }}_{-1}^{n+1},\cdots ,{\mathrm{\Lambda }}_{N+2}^{n+1})}^{\prime }.$$

Similarly, one can construct matrix for $S_2$. For $n=0$, the unknown $\mathrm{\Lambda }$ can be determined as:

$$K_1{\mathrm{\Lambda }}^0=K_2,$$

(17)

where

$$K_1=\left[\begin{array}{cccccccccc}{\displaystyle \frac{-5}{h}}& {\displaystyle \frac{-50}{h}}& 0& {\displaystyle \frac{50}{h}}& {\displaystyle \frac{5}{h}}& \cdots & \cdots & \cdots & \cdots & 0\\ {\displaystyle \frac{20}{h^2}}& {\displaystyle \frac{40}{h^2}}& {\displaystyle \frac{-120}{h^2}}& {\displaystyle \frac{40}{h^2}}& {\displaystyle \frac{20}{h^2}}& \ddots & & & & ⋮\\ 1& 26& 66& 26& 1& \ddots & \ddots & & & ⋮\\ 0& 1& 26& 66& 26& 1& \ddots & \ddots & & ⋮\\ ⋮& \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ ⋮& & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ ⋮& & & \ddots & 1& 26& 66& 26& 1& 0\\ ⋮& & & & \ddots & 1& 26& 66& 26& 1\\ 0& \cdots & \cdots & \cdots & \cdots & {\displaystyle \frac{20}{h^2}}& {\displaystyle \frac{40}{h^2}}& {\displaystyle \frac{-120}{h^2}}& {\displaystyle \frac{40}{h^2}}& {\displaystyle \frac{20}{h^2}}\\ 0& \cdots & \cdots & \cdots & \cdots & {\displaystyle \frac{-5}{h}}& {\displaystyle \frac{-50}{h}}& 0& {\displaystyle \frac{50}{h}}& {\displaystyle \frac{5}{h}}\end{array}\right],$$

(18)

$${\mathrm{\Lambda }}^0=\left[\begin{array}{c}{\mathrm{\Lambda }}_{-2}^0\\ {\mathrm{\Lambda }}_{-1}^0\\ {\mathrm{\Lambda }}_0^0\\ {\mathrm{\Lambda }}_1^0\\ ⋮\\ {\mathrm{\Lambda }}_{N-1}^0\\ {\mathrm{\Lambda }}_N^0\\ {\mathrm{\Lambda }}_{N+1}^0\\ {\mathrm{\Lambda }}_{N+2}^0\end{array}\right],K_2=\left[\begin{array}{c}{f}^{\prime }\left({\varrho }_0\right)\\ f^{\prime \prime }\left({\varrho }_0\right)\\ f\left({\varrho }_0\right)\\ f\left({\varrho }_1\right)\\ f\left({\varrho }_2\right)\\ ⋮\\ f\left({\varrho }_{N-1}\right)\\ f\left({\varrho }_N\right)\\ f^{\prime \prime }\left({\varrho }_N\right)\\ f^{\prime }\left({\varrho }_N\right)\end{array}\right].$$

By solving Equation (15), one can obtain the required unknowns. Once these unknowns are determined, one can put them into Equation (8) to obtain the required solution.

2.2. Algorithm of the Method

To find the solution of the given equation using QBS, the following steps were followed:

• Input $N;{\varrho }_0,\cdots ,{\varrho }_N\mathrm{and}n;t_0,\cdots ,t_n$
• Output To calculate the unknown ${\mathrm{\Lambda }}^{\prime }s$.
• Step 1 for $i=0,1,\cdots ,N\mathrm{and}j=0,1,\cdots ,n,h={\varrho }_{i+1}-{\varrho }_i,\delta \mathrm{t}=t_{j+1}-t_j$.
• Step 2 Input values of $B\left({\varrho }_i\right),B^{\prime }\left({\varrho }_i\right),B^{\prime \prime }\left({\varrho }_i\right),B^{\prime \prime \prime }\left({\varrho }_i\right).$
• Note ${\displaystyle \mathrm{U}(\varrho ,\mathrm{t})=\sum _{i=-2}^{N+2}{\mathrm{\Lambda }}_i\left(\mathrm{t}\right)B_i\left(\varrho \right).}$
• Set $i=0,1,\cdots ,N$, $U({\varrho }_i,0)=f\left({\varrho }_i\right)$.
• Calculate ${\mathrm{\Lambda }}^0=K_1^{-1}{K}_2$.

Use (15) and repeat the process to calculate the ${\mathrm{\Lambda }}^{n}s$ for each time level and the solution.

3. Stability Analysis

This section explores the von Neumann method to assess the stability of the approach. Rewriting Equation (10), we have:

$$\begin{array}{cc}\hfill & {\displaystyle \frac{{\mathrm{U}}^{n+1}-{\mathrm{U}}^n}{\delta \mathrm{t}}}-\nu {\displaystyle \frac{{\mathrm{U}}_{\varrho \varrho }^{n+1}-{\mathrm{U}}_{\varrho \varrho }^n}{\delta \mathrm{t}}}+\kappa ({\mathrm{U}}_{\varrho }^{n+1}+{\mathrm{U}}_{\varrho }^n)-{\displaystyle \frac{1}{2}}({\left(\mathrm{U}{\mathrm{U}}_{3\varrho }\right)}^{n+1}+{\left(\mathrm{U}{\mathrm{U}}_{3\varrho }\right)}^n)\hfill \\ \hfill & -{\displaystyle \frac{\alpha }{2}}({\left(\mathrm{U}{\mathrm{U}}_{\varrho }\right)}^{n+1}+{\left(\mathrm{U}{\mathrm{U}}_{\varrho }\right)}^n)-{\displaystyle \frac{\beta }{2}}({\left({\mathrm{U}}_{\varrho }{\mathrm{U}}_{\varrho \varrho }\right)}^{n+1}+{\left({\mathrm{U}}_{\varrho }{\mathrm{U}}_{\varrho \varrho }\right)}^n)=0.\hfill \end{array}$$

(19)

The non-linear terms are transformed into linear ones by taking $\mathrm{U}$ as a local constant $\lambda$, and by further simplification, they give:

$$\begin{array}{cc}\hfill & {\mathrm{U}}^{n+1}-\nu {\mathrm{U}}_{\varrho \varrho }^{n+1}+\kappa \delta \mathrm{t}{\mathrm{U}}_{\varrho }^{n+1}-{\displaystyle \frac{\lambda \delta \mathrm{t}}{2}}{\mathrm{U}}_{3\varrho }^{n+1}-{\displaystyle \frac{\alpha \lambda \delta \mathrm{t}}{2}}{\mathrm{U}}_{\varrho }^{n+1}=\hfill \\ \hfill & {\mathrm{U}}^n-\nu {\mathrm{U}}_{\varrho \varrho }^n-\kappa \delta \mathrm{t}{\mathrm{U}}_{\varrho }^n+{\displaystyle \frac{\lambda \delta \mathrm{t}}{2}}{\mathrm{U}}_{3\varrho }^m+{\displaystyle \frac{\alpha \lambda \delta \mathrm{t}}{2}}{\mathrm{U}}_{\varrho }^n.\hfill \end{array}$$

(20)

Evaluating (20) at a grid point by substituting the values of $\mathrm{U}$, ${\mathrm{U}}_{\varrho }$, ${\mathrm{U}}_{\varrho \varrho }$ and ${\mathrm{U}}_{3\varrho }$ from (14), we obtain:

$$\begin{array}{cc}\hfill & c_1{\mathrm{\Lambda }}_{j-2}^{n+1}+c_2{\mathrm{\Lambda }}_{j-1}^{n+1}+c_3{\mathrm{\Lambda }}_j^{n+1}+c_4{\mathrm{\Lambda }}_{j+1}^{n+1}+c_5{\mathrm{\Lambda }}_{j+2}^{n+1}\hfill \\ \hfill & =c_5{\mathrm{\Lambda }}_{j-2}^n+c_4{\mathrm{\Lambda }}_{j-1}^n+c_3{\mathrm{\Lambda }}_j^n+c_2{\mathrm{\Lambda }}_{j+1}^n+c_1{\mathrm{\Lambda }}_{j+2}^n,\hfill \end{array}$$

(21)

where

$$\begin{array}{cc}\hfill & c_1=\left(1-{\displaystyle \frac{20\nu }{h^2}}-{\displaystyle \frac{5}{h}}(\kappa \delta \mathrm{t}-{\displaystyle \frac{\alpha \lambda \delta \mathrm{t}}{2}})+{\displaystyle \frac{30\lambda \delta \mathrm{t}}{h^3}}\right),c_2=\left(26-{\displaystyle \frac{40\nu }{h^2}}-{\displaystyle \frac{50}{h}}(\kappa \delta \mathrm{t}-{\displaystyle \frac{\alpha \lambda \delta \mathrm{t}}{2}})-{\displaystyle \frac{60\lambda \delta \mathrm{t}}{h^3}}\right),\hfill \\ \hfill & c_3=\left(66+{\displaystyle \frac{120\nu }{h^2}}\right),c_4=\left(26-{\displaystyle \frac{40\nu }{h^2}}+{\displaystyle \frac{50}{h}}(\kappa \delta \mathrm{t}-{\displaystyle \frac{\alpha \lambda \delta \mathrm{t}}{2}})+{\displaystyle \frac{60\lambda \delta \mathrm{t}}{h^3}}\right),\hfill \\ \hfill & c_5=\left(1-{\displaystyle \frac{20\nu }{h^2}}+{\displaystyle \frac{5}{h}}(\kappa \delta \mathrm{t}-{\displaystyle \frac{\alpha \lambda \delta \mathrm{t}}{2}})-{\displaystyle \frac{30\lambda \delta \mathrm{t}}{h^3}}\right).\hfill \end{array}$$

(22)

Substituting

$${\mathrm{\Lambda }}_j^n={\eta }^n{e}^{ijkh},$$

(23)

where $\eta$ is the amplification factor, i is the imaginary unit, and k and h are the mode number and mesh width, respectively. By inserting Equation (23) in (21), we obtain

$$|x+iy|{\eta }^{n+1}=|x-iy|{\eta }^n,$$

(24)

where

$$\begin{array}{c}\hfill \left\{\begin{array}{c}x=(2-{\displaystyle \frac{40\nu }{h^2}})cos\left(2kh\right)+(52-{\displaystyle \frac{80\nu }{h^2}})cos\left(kh\right)+(66+{\displaystyle \frac{120\nu }{h^2}}),\hfill \\ y=({\displaystyle \frac{10}{h}}(\kappa \delta \mathrm{t}-{\displaystyle \frac{\alpha \lambda \delta \mathrm{t}}{2}})-{\displaystyle \frac{60\lambda \delta \mathrm{t}}{h^3}})isin\left(2kh\right)+({\displaystyle \frac{100}{h}}(\kappa \delta \mathrm{t}-{\displaystyle \frac{\alpha \lambda \delta \mathrm{t}}{2}})+{\displaystyle \frac{120\lambda \delta \mathrm{t}}{h^3}})isin\left(kh\right).\hfill \end{array}\right.\end{array}$$

(25)

It is clear from the values of x and y that the amplification factor $\left|\eta \right|=1$, and the scheme is unconditionally stable for the linearized equation.

4. Numerical Examples

To validate the precision of the method, we employ absolute and relative error norms.

• Absolute error:

  $$A.E=|{\mathrm{U}}_e({\varrho }_i,{\mathrm{t}}_{\mathrm{i}})-{\mathrm{U}}_a({\varrho }_i,{\mathrm{t}}_{\mathrm{i}})|,$$
• Relative error norm:

  $$R.E={\displaystyle \frac{|{\mathrm{U}}_e({\varrho }_i,{\mathrm{t}}_{\mathrm{i}})-{\mathrm{U}}_a({\varrho }_i,{\mathrm{t}}_{\mathrm{i}})|}{|{\mathrm{U}}_e({\varrho }_i,{\mathrm{t}}_{\mathrm{i}})|}}.$$
• Rate of convergence:

  $$\mathrm{Rate}={\displaystyle \frac{L_{\infty }(\delta t_i)}{L_{\infty }(\delta t_{i+1})}}.$$

The software that is used for the computational purposes is MATLAB (R2020a), also used to compute the convergence rate (CR) and spectral radius (SR).

Example 1.

Consider the RH equation:

$$\begin{array}{cc}\hfill & {\mathrm{U}}_{\mathrm{t}}-\mathrm{U}{\mathrm{U}}_{3\varrho }-\mathrm{U}{\mathrm{U}}_{\varrho }-3{\mathrm{U}}_{\varrho }{\mathrm{U}}_{\varrho \varrho }=0,\hfill \\ \hfill & \mathrm{t}>0,\varrho \in [0,1],\hfill \end{array}$$

(26)

along with the initial and boundary conditions:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\mathrm{U}(\varrho ,0)={\displaystyle \frac{-8c}{3}}{cos}^2\left({\displaystyle \frac{\varrho }{4}}\right),\hfill \\ \mathrm{U}(0,\mathrm{t})={\displaystyle \frac{-8c}{3}}{cos}^2\left({\displaystyle \frac{-\mathrm{t}}{4}}\right),\mathrm{U}(1,\mathrm{t})={\displaystyle \frac{-8c}{3}}{cos}^2\left({\displaystyle \frac{1-\mathrm{t}}{4}}\right),\hfill \\ {\mathrm{U}}_{\varrho }(0,\mathrm{t})={\displaystyle \frac{16c}{3}}sin2\left({\displaystyle \frac{-\mathrm{t}}{4}}\right),{\mathrm{U}}_{\varrho }(1,\mathrm{t})={\displaystyle \frac{16c}{3}}sin2\left({\displaystyle \frac{1-\mathrm{t}}{4}}\right).\hfill \end{array}\right.\end{array}$$

(27)

The exact solution is $\mathrm{U}(\varrho ,\mathrm{t})={\displaystyle \frac{-8c}{3}}{cos}^2\left({\displaystyle \frac{\varrho -\mathrm{ct}}{4}}\right)$. The results are given in

Table 1. Comparison of absolute error of the given method with the hybrid and RDTM methods of Example 1 at $\delta \mathrm{t}=0.00001$, $h=0.1$, $t=0.0001$, $c=1$ and $\theta =1/2$.

$\mathit{\varrho }$	Quintic	[25]	[25]
0	4.441 × 10${}^{-16}$	9.000 × 10${}^{-9}$	2.504 × 10${}^{-4}$
0.1	3.908 × 10${}^{-14}$	1.111 × 10${}^{-5}$	2.501 × 10${}^{-4}$
0.2	5.196 × 10${}^{-14}$	2.17 × 10${}^{-5}$	2.495 × 10${}^{-4}$
0.3	9.104 × 10${}^{-14}$	3.314 × 10${}^{-5}$	2.486 × 10${}^{-4}$
0.4	1.101 × 10${}^{-13}$	4.398 × 10${}^{-5}$	2.475 × 10${}^{-4}$
0.5	1.492 × 10${}^{-13}$	5.465 × 10${}^{-5}$	2.460 × 10${}^{-4}$
0.6	1.652 × 10${}^{-13}$	6.509 × 10${}^{-5}$	2.443 × 10${}^{-4}$
0.7	1.847 × 10${}^{-13}$	7.528 × 10${}^{-5}$	2.423 × 10${}^{-4}$
0.8	2.802 × 10${}^{-13}$	8.518 × 10${}^{-5}$	2.401 × 10${}^{-4}$
0.9	2.385 × 10${}^{-13}$	9.474 × 10${}^{-5}$	2.375 × 10${}^{-4}$
1	4.441 × 10${}^{-16}$	1.039 × 10${}^{-4}$	2.347 × 10${}^{-4}$

and

Table 2. Comparison of error norms at $\delta t=0.01$ and $h=0.1$.

t	c	${\mathit{L}}_2$	[31]	${\mathit{L}}_{\mathit{\infty }}$	[31]
0.1	0.5	4.86 × 10${}^{-9}$	2.17 × 10${}^{-4}$	3.22 × 10${}^{-9}$	6.67 × 10${}^{-5}$
0.5	0.5	1.04 × 10${}^{-8}$	1.08 × 10${}^{-3}$	6.87 × 10${}^{-9}$	3.33 × 10${}^3$
1	0.5	1.11 × 10${}^{-8}$	2.17 × 10${}^{-3}$	7.22 × 10${}^{-9}$	6.67 × 10${}^{-4}$
3	0.01	1.09 × 10${}^{-12}$	6.50 × 10${}^{-5}$	7.13 × 10${}^{-13}$	2.00 × 10${}^{-5}$
10	0.01	2.51 × 10${}^{-12}$	2.17 × 10${}^{-4}$	1.67 × 10${}^{-12}$	6.67 × 10${}^{-5}$
50	0.01	3.70 × 10${}^{-12}$	1.08 × 10${}^{-3}$	2.41 × 10${}^{-12}$	3.32 × 10${}^{-4}$

for different values of ${\varrho }_i$, t and $\theta =1/2$. In Table 1, the absolute error norm is used to check the performance of the technique. The proposed method is more efficient than hybrid and RDTM. In Table 2, we compute the simulations for various values of t and c, and the results of the proposed method are better than [25]. In Figure 1, the dashed line and asterisk show the approximate and exact solutions, respectively. Two-dimensional heat maps for the solution profiles are given at $t=1$. From Figure 1, it is concluded that the heat maps plots of the proposed method agree with the exact method.

Example 2.

Consider the RH equation given as

$$\begin{array}{cc}\hfill & {\mathrm{U}}_{\mathrm{t}}+0.8{\mathrm{U}}^3{\mathrm{U}}_{\varrho }+2{\mathrm{U}}^3{\mathrm{U}}_{3\varrho }+18{\mathrm{U}}^2{\mathrm{U}}_{\varrho }{\mathrm{U}}_{\varrho \varrho }+12\mathrm{U}{\left({\mathrm{U}}_{\varrho }\right)}^3=0,\hfill \\ \hfill & \mathrm{t}>0,\varrho \in [0,1].\hfill \end{array}$$

(28)

The initial and boundary conditions are

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\mathrm{U}(\varrho ,0)=2{\left(c{cosh}^2\left({\displaystyle \frac{3i\left(\varrho \right)\sqrt{0.4}}{8}}\right)\right)}^{{\displaystyle \frac{1}{3}}},\hfill \\ \mathrm{U}(0,\mathrm{t})=2{\left(c{cosh}^2\left({\displaystyle \frac{3i(-c\mathrm{t})\sqrt{0.4}}{8}}\right)\right)}^{{\displaystyle \frac{1}{3}}},\mathrm{U}(1,\mathrm{t})=2{\left(c{cosh}^2\left({\displaystyle \frac{3i(1-c\mathrm{t})\sqrt{0.4}}{8}}\right)\right)}^{{\displaystyle \frac{1}{3}}},\hfill \\ {\mathrm{U}}_{\varrho }(0,\mathrm{t})=i{\displaystyle \frac{\sqrt{0.4}}{2}}2{\left(c{cosh}^2\left({\displaystyle \frac{3i(-c\mathrm{t})\sqrt{0.4}}{8}}\right)\right)}^{{\displaystyle \frac{1}{3}}}2{\left(ctanh\left({\displaystyle \frac{3i(-c\mathrm{t})\sqrt{0.4}}{8}}\right)\right)}^{{\displaystyle \frac{1}{3}}},\hfill \\ {\mathrm{U}}_{\varrho }(1,\mathrm{t})=i{\displaystyle \frac{\sqrt{0.4}}{2}}2{\left(c{cosh}^2\left({\displaystyle \frac{3i(1-c\mathrm{t})\sqrt{0.4}}{8}}\right)\right)}^{{\displaystyle \frac{1}{3}}}2{\left(ctanh\left({\displaystyle \frac{3i(1-c\mathrm{t})\sqrt{0.4}}{8}}\right)\right)}^{{\displaystyle \frac{1}{3}}}.\hfill \end{array}\right.\end{array}$$

(29)

The exact solution is $\mathrm{U}(\varrho ,\mathrm{t})=2{\left(c{cosh}^2\left({\displaystyle \frac{3i(\varrho -c\mathrm{t})\sqrt{0.4}}{8}}\right)\right)}^{{\displaystyle \frac{1}{3}}}$, parameter $c=0.1$. The absolute error and R.E error norms are calculated to measure the accuracy of the method.

Table 3. Comparison of the given method with the Genocchi wavelet method using the absolute and relative error at $\delta \mathrm{t}=0.01$, $h=0.1$, and $\theta =1/2$.

		Absolute Error		Relative Error
$\mathbf{\varrho }$	$\mathrm{t}$	Quintic	[32]	Quintic	[32]
0	0	1.11 × 10${}^{-16}$	2.31 × 10${}^{-11}$	1.20 × 10${}^{-16}$	2.48 × 10${}^{-11}$
0	0.2	5.55 × 10${}^{-16}$	2.57 × 10${}^{-10}$	5.98 × 10${}^{-16}$	2.76 × 10${}^{-10}$
0	0.4	1.11 × 10${}^{-16}$	4.90 × 10${}^{-10}$	1.20 × 10${}^{-16}$	5.28 × 10${}^{-10}$
0	0.6	7.77 × 10${}^{-16}$	6.52 × 10${}^{-10}$	8.37 × 10${}^{-16}$	7.03 × 10${}^{-10}$
0	0.8	3.33 × 10${}^{-16}$	7.62 × 10${}^{-10}$	3.59 × 10${}^{-16}$	8.21 × 10${}^{-10}$
0	1	8.88 × 10${}^{-16}$	9.99 × 10${}^{-16}$	9.57 × 10${}^{-16}$	1.08 × 10${}^{-15}$
0.5	0	1.11 × 10${}^{-16}$	7.10 × 10${}^{-11}$	1.20 × 10${}^{-16}$	7.69 × 10${}^{-11}$
0.5	0.2	1.52 × 10${}^{-10}$	5.30 × 10${}^{-9}$	1.64 × 10${}^{-10}$	5.73 × 10${}^{-9}$
0.5	0.4	1.81 × 10${}^{-10}$	9.63 × 10${}^{-9}$	1.95 × 10${}^{-10}$	1.04 × 10${}^{-8}$
0.5	0.6	6.65 × 10${}^{-11}$	1.35 × 10${}^{-8}$	7.19 × 10${}^{-11}$	1.45 × 10${}^{-8}$
0.5	0.8	2.34 × 10${}^{-10}$	1.68 × 10${}^{-8}$	2.53 × 10${}^{-10}$	1.82 × 10${}^{-8}$
0.5	1	2.49 × 10${}^{-11}$	4.39 × 10${}^{-8}$	2.69 × 10${}^{-11}$	4.74 × 10${}^{-8}$
1	0	2.22 × 10${}^{-16}$	9.99 × 10${}^{-16}$	2.44 × 10${}^{-16}$	1.10 × 10${}^{-15}$
1	0.2	2.22 × 10${}^{-16}$	0.00 × 10${}^{00}$	2.44 × 10${}^{-16}$	0.00 × 10${}^{00}$
1	0.4	5.55 × 10${}^{-16}$	9.99 × 10${}^{-16}$	6.08 × 10${}^{-16}$	1.10 × 10${}^{-15}$
1	0.6	0.00 × 10${}^{00}$	1.11 × 10${}^{-16}$	0.00 × 10${}^{00}$	1.22 × 10${}^{-16}$
1	0.8	1.11 × 10${}^{-16}$	1.11 × 10${}^{-15}$	1.22 × 10${}^{-16}$	1.22 × 10${}^{-15}$
1	1	1.11 × 10${}^{-16}$	0.00 × 10${}^{00}$	1.21 × 10${}^{-16}$	0.00 × 10${}^{00}$

shows that the QBS results are better than the Genocchi wavelets method. The temporal and spatial step sizes are $\delta \mathrm{t}=0.01$ and $h=0.1$, respectively. Figure 2 shows the heat map profiles of the exact solutions, QBS and absolute error, and the exact vs. approximate plot at $t=1$. The curve of the proposed method is validated with the exact one. By increasing the values of $\varrho$, the error decreases.

Table 4. Convergence and stability analysis of the proposed method for temporal step size with computational time $t=1$ and $h=0.1$.

$\mathit{\delta }\mathbf{t}$	${\mathit{L}}_{\mathit{\infty }}$	C.R	CPU(s)	S.R
0.050	5.69 × 10${}^{-9}$		0.10462	0.0381
0.025	2.35 × 10${}^{-9}$	2.42	0.111444	0.0436
0.013	1.18 × 10${}^{-9}$	1.98	0.129875	0.0469
0.010	5.47 × 10${}^{-10}$	2.17	0.132137	0.0476

shows the values of the spectral radius, which is less than one for different temporal step sizes, results which lead to stability. Furthermore, the computational time and convergence rate are given, which show that as the temporal step size decreases, the results converge.

Example 3.

Consider a FW equation

$$\begin{array}{cc}\hfill & {\mathrm{U}}_{\mathrm{t}}-{\mathrm{U}}_{\varrho \varrho \mathrm{t}}+{\mathrm{U}}_{\varrho }-\mathrm{U}{\mathrm{U}}_{3\varrho }+\mathrm{U}{\mathrm{U}}_{\varrho }-3{\mathrm{U}}_{\varrho }{\mathrm{U}}_{\varrho \varrho }=0.\hfill \\ \hfill & \mathrm{t}>0,\varrho \in [-5,5].\hfill \end{array}$$

(30)

Subjected to ICs and BCs,

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\mathrm{U}(\varrho ,0)={\displaystyle \frac{4}{3}}exp\left({\displaystyle \frac{1}{2}}\varrho \right),\hfill \\ \mathrm{U}(a,\mathrm{t})={\displaystyle \frac{4}{3}}exp({\displaystyle \frac{1}{2}}a-{\displaystyle \frac{2}{3}}\mathrm{t}),\mathrm{U}(b,\mathrm{t})={\displaystyle \frac{4}{3}}exp({\displaystyle \frac{1}{2}}b-{\displaystyle \frac{2}{3}}\mathrm{t}),\hfill \\ {\mathrm{U}}_{\varrho }(a,\mathrm{t})={\displaystyle \frac{2}{3}}exp({\displaystyle \frac{1}{2}}a-{\displaystyle \frac{2}{3}}\mathrm{t}),{\mathrm{U}}_{\varrho }(b,\mathrm{t})={\displaystyle \frac{2}{3}}exp({\displaystyle \frac{1}{2}}b-{\displaystyle \frac{2}{3}}\mathrm{t}).\hfill \end{array}\right.\end{array}$$

(31)

The exact solution is $\mathrm{U}(\varrho ,\mathrm{t})={\displaystyle \frac{4}{3}}exp({\displaystyle \frac{1}{2}}\varrho -{\displaystyle \frac{2}{3}}\mathrm{t})$. For $\theta ={\displaystyle \frac{1}{2}}$, the absolute error norm is given in

Table 5. Comparison of the given method with different methods using the absolute error at $\delta \mathrm{t}=0.01$, $h=0.1$ and $\theta =1/2$.

	$\mathit{\varrho }=-5$			$\mathit{\varrho }=-2.5$		$\mathit{\varrho }=2.5$		$\mathit{\varrho }=5$
$\mathbf{t}$	Quintic	[30]	[34]	Quintic	[34]	Quintic	[34]	Quintic	[30]	[34]
0.02	4.16 × 10${}^{-17}$	1.43 × 10${}^{-13}$	3.13 × 10${}^{-7}$	2.64 × 10${}^{-8}$	1.10 × 10${}^{-6}$	2.25 × 10${}^{-7}$	1.43 × 10${}^{-5}$	3.55 × 10${}^{-15}$	2.12 × 10${}^{-11}$	5.85 × 10${}^{-5}$
0.04	9.71 × 10${}^{-17}$	4.65 × 10${}^{-13}$	5.47 × 10${}^{-7}$	5.21 × 10${}^{-8}$	1.93 × 10${}^{-6}$	4.49 × 10${}^{-7}$	2.69 × 10${}^{-5}$	3.55 × 10${}^{-15}$	6.90 × 10${}^{-11}$	1.26 × 10${}^{-4}$
0.06	8.33 × 10${}^{-17}$	6.71 × 10${}^{-13}$	7.09 × 10${}^{-7}$	7.69 × 10${}^{-8}$	2.51 × 10${}^{-6}$	6.47 × 10${}^{-7}$	3.78 × 10${}^{-5}$	0.00 × 10${}^{00}$	9.96 × 10${}^{-11}$	2.00 × 10${}^{-4}$
0.08	0	1.97 × 10${}^{-13}$	8.08 × 10${}^{-7}$	1.01 × 10${}^{-7}$	2.88 × 10${}^{-6}$	8.24 × 10${}^{-7}$	4.71 × 10${}^{-5}$	1.78 × 10${}^{-15}$	2.93 × 10${}^{-11}$	2.78 × 10${}^{-4}$
0.1	2.78 × 10${}^{-17}$	7.25 × 10${}^{-13}$	8.52 × 10${}^{-7}$	1.24 × 10${}^{-7}$	3.07 × 10${}^{-6}$	9.84 × 10${}^{-7}$	5.49 × 10${}^{-5}$	1.78 × 10${}^{-15}$	1.08 × 10${}^{-10}$	3.57 × 10${}^{-4}$

and

Table 6. Comparison of the given method with different methods using the absolute error at $\delta \mathrm{t}=0.1$, $h=0.1$ and $\theta =1/2$.

	$\mathit{\varrho }=-4$			$\mathit{\varrho }=0$		$\mathit{\varrho }=4$
t	Quintic	[33]	[33]	Quintic	[33]	Quintic	[33]	[33]
0.2	1.67 × 10${}^{-16}$	6.75 × 10${}^{-8}$	1.72 × 10${}^{-5}$	4.52 × 10${}^{-5}$	1.11 × 10${}^{-4}$	7.11 × 10${}^{-15}$	8.29 × 10${}^{-8}$	3.04 × 10${}^{-5}$
0.4	2.22 × 10${}^{-16}$	2.31 × 10${}^{-7}$	3.35 × 10${}^{-6}$	7.46 × 10${}^{-5}$	5.72 × 10${}^{-5}$	5.33 × 10${}^{-15}$	1.04 × 10${}^{-9}$	5.42 × 10${}^{-6}$
0.6	2.36 × 10${}^{-16}$	6.72 × 10${}^{-9}$	1.09 × 10${}^{-5}$	9.41 × 10${}^{-5}$	9.64 × 10${}^{-5}$	4.44 × 10${}^{-15}$	7.31 × 10${}^{-9}$	1.37 × 10${}^{-5}$
0.8	5.55 × 10${}^{-17}$	3.92 × 10${}^{-7}$	3.09 × 10${}^{-5}$	1.07 × 10${}^{-4}$	4.69 × 10${}^{-5}$	4.44 × 10${}^{-15}$	3.34 × 10${}^{-7}$	1.37 × 10${}^{-5}$
1	1.25 × 10${}^{-16}$	1.39 × 10${}^{-7}$	1.02 × 10${}^{-5}$	1.14 × 10${}^{-4}$	9.43 × 10${}^{-5}$	3.55 × 10${}^{-15}$	5.78 × 10${}^{-8}$	1.23 × 10${}^{-5}$
5	2.60 × 10${}^{-18}$	1.39 × 10${}^{-7}$	3.06 × 10${}^{-5}$	4.88 × 10${}^{-5}$	1.88 × 10${}^{-4}$	2.78 × 10${}^{-16}$	7.54 × 10${}^{-9}$	1.75 × 10${}^{-3}$

, respectively. In Table 5, the values of the proposed scheme in comparison with VIM are given for temporal step size $\delta \mathrm{t}=0.01$, which demonstrates that the proposed methodology exhibits improved accuracy. In Table 6, our results are compared with the RBF method. The results from QBS are superior to those indicated in [33] at final time $\mathrm{t}=1$. Moreover, at the larger time, the proposed method attains good accuracy. In this case, the temporal step size is $\delta \mathrm{t}=0.1$ and the collocation points are $N=20$. The heat map profiles of the absolute error, exact and approximate solutions along with the plot of the exact vs. approximate solution plot at $t=1$ are give in Figure 3, which shows that the QBS results correspond to the exact solution.

Table 7. Convergence and stability analysis of the proposed method for different temporal step sizes with computational time.

$\mathbf{t}$	${\mathit{L}}_{\mathit{\infty }}$	[33]	CPU(s)	[33]	$\mathit{\delta }\mathbf{t}$	${\mathit{L}}_{\mathit{\infty }}$	C.R	S.R
0.2	1.25 × 10${}^{-4}$	2.39 × 10${}^{-4}$	0.193664	0.525753	0.100	2.30 × 10${}^{-4}$	-	0.4392
0.4	1.78 × 10${}^{-4}$	3.66 × 10${}^{-4}$	0.206210	0.535698	0.03	1.56 × 10${}^{-5}$	3.88	0.6216
0.8	2.21 × 10${}^{-4}$	6.00 × 10${}^{-4}$	0.231908	0.542491	0.017	7.66 × 10${}^{-6}$	2.04	0.7015
1	2.30 × 10${}^{-4}$	7.06 × 10${}^{-4}$	0.221241	0.559415	0.013	4.87 × 10${}^{-6}$	1.57	0.7585
5	2.55 × 10${}^{-4}$	1.45 × 10${}^{-3}$	0.567964	0.614328	0.010	3.58 × 10${}^{-6}$	1.36	0.8009

shows the convergence rate and spectral radius values along with the computational time. Furthermore, the $L_{\infty }$ error norm and computational time show that the proposed method results are better than [33].

Example 4.

Consider a FW equation

$$\begin{array}{cc}\hfill & {\mathrm{U}}_{\mathrm{t}}-{\mathrm{U}}_{\varrho \varrho \mathrm{t}}+{\mathrm{U}}_{\varrho }-\mathrm{U}{\mathrm{U}}_{3\varrho }+\mathrm{U}{\mathrm{U}}_{\varrho }-3{\mathrm{U}}_{\varrho }{\mathrm{U}}_{\varrho \varrho }=0,\hfill \\ \hfill & \mathrm{t}>0,\varrho \in [-4,4],\hfill \end{array}$$

(32)

subjected to initial and boundary conditions

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\mathrm{U}(\varrho ,0)=exp\left({\displaystyle \frac{1}{2}}\varrho \right),\hfill \\ \mathrm{U}(-4,\mathrm{t})=exp(-2-{\displaystyle \frac{2}{3}}\mathrm{t}),\mathrm{U}(4,\mathrm{t})=exp(2-{\displaystyle \frac{2}{3}}\mathrm{t}),\hfill \\ {\mathrm{U}}_{\varrho }(-4,\mathrm{t})={\displaystyle \frac{1}{2}}exp(-2-{\displaystyle \frac{2}{3}}\mathrm{t}),{\mathrm{U}}_{\varrho }(4,\mathrm{t})={\displaystyle \frac{1}{2}}exp(2-{\displaystyle \frac{2}{3}}\mathrm{t}).\hfill \end{array}\right.\end{array}$$

(33)

The exact solution is $\mathrm{U}(\varrho ,\mathrm{t})=exp({\displaystyle \frac{1}{2}}\varrho -{\displaystyle \frac{2}{3}}\mathrm{t})$. The absolute error norm of the proposed method in comparison with HAM at different space and time values is given in

Table 8. Comparison of the given method with HAM using the absolute error at $\delta \mathrm{t}=0.01$, $h=0.1$ and $\theta =1/2$.

	$\mathbf{t}=0.2$		$\mathbf{t}=0.4$		$\mathbf{t}=0.6$		$\mathbf{t}=0.8$		$\mathbf{t}=1$
$\mathbf{\varrho }$	Quintic	[28]	Quintic	[28]	Quintic	[28]	Quintic	[28]	Quintic	[28]
−4	1.11 × 10${}^{-16}$	2.22 × 10${}^{-5}$	1.11 × 10${}^{-16}$	9.47 × 10${}^{-6}$	0.00 × 10${}^{00}$	4.84 × 10${}^{-5}$	4.16 × 10${}^{-17}$	6.72 × 10${}^{-5}$	1.11 × 10${}^{-16}$	5.36 × 10${}^{-5}$
−2	1.66 × 10${}^{-7}$	6.04 × 10${}^{-5}$	2.85 × 10${}^{-7}$	2.58 × 10${}^{-5}$	3.68 × 10${}^{-7}$	1.32 × 10${}^{-4}$	4.26 × 10${}^{-7}$	1.83 × 10${}^{-4}$	4.62 × 10${}^{-7}$	1.46 × 10${}^{-4}$
0	4.20 × 10${}^{-7}$	1.64 × 10${}^{-5}$	7.07 × 10${}^{-7}$	7.00 × 10${}^{-5}$	9.08 × 10${}^{-7}$	3.58 × 10${}^{-4}$	1.05 × 10${}^{-6}$	4.96 × 10${}^{-4}$	1.14 × 10${}^{-6}$	3.96 × 10${}^{-4}$
2	7.57 × 10${}^{-7}$	4.46 × 10${}^{-4}$	1.27 × 10${}^{-6}$	1.90 × 10${}^{-4}$	1.64 × 10${}^{-6}$	9.72 × 10${}^{-4}$	1.85 × 10${}^{-6}$	1.35 × 10${}^{-3}$	1.99 × 10${}^{-6}$	1.08 × 10${}^{-3}$
4	5.33 × 10${}^{-15}$	1.21 × 10${}^{-4}$	4.44 × 10${}^{-15}$	5.17 × 10${}^{-4}$	3.55 × 10${}^{-15}$	2.64 × 10${}^{-3}$	2.66 × 10${}^{-15}$	3.67 × 10${}^{-3}$	3.55 × 10${}^{-15}$	2.93 × 10${}^{-3}$

. The results show that the QBS is better than HAM for different ${\varrho }_i$’s and ${\mathrm{t}}_i$’s. Heat map plots of the absolute error, exact and approximate are given in Figure 4. From Figure 4, it is clear that the plots of the proposed solution overlap with the exact solution at time $\mathrm{t}=1$.

5. Conclusions

In the present study, we applied the QBS technique to the computational analysis of RH and FW equations. Initially, the temporal component was discretized using finite differences, and functions along with their derivatives were approximated via the QBS technique. The results were computed using various error norms, presented in tables for comparison. Our findings demonstrated better performance compared to previously published results. Through von Neumann stability analysis, we confirmed the stability of the proposed methodology. Heat map plots, alongside plots of approximate solution, showed excellent agreement with the exact solutions. This method is observed to be simple, fast, and accurate for handling certain non-linear PDEs.