A Novel Hybrid Computational Technique to Study Conformable Burgers’ Equation

Abstract

A fully discrete computational technique involving the implicit finite difference technique and cubic Hermite splines is proposed to solve the non-linear conformable damped Burgers’ equation with variable coefficients numerically. The proposed scheme is capable of solving the equation having singularity at $t=0$. The space direction is discretized using cubic Hermite splines, whereas the time direction is discretized using an implicit finite difference scheme. The convergence, stability and error estimates of the proposed scheme are discussed in detail to prove the efficiency of the technique. The convergence of the proposed scheme is found to be of order $h^2$ in space and order ${(\Delta t)}^{\alpha }$ in the time direction. The efficiency of the proposed scheme is verified by calculating error norms in the Eucledian and supremum sense. The proposed technique is applied on conformable damped Burgers’ equation with different initial and boundary conditions and the results are presented as tables and graphs. Comparison with results already in the literature also validates the application of the proposed technique.

1. Introduction

Many real-world phenomena are described by non-linear partial differential equations. As a result, these equations have received a great deal of attention from scholars in the field of science. Some of these non-linear equations involve the interaction between convection and diffusion or the interplay between reaction and diffusion [1,2,3]. These non-linear equations become difficult to solve in cases where an exact solution is either non-existent or too complicated to be analyzed thoroughly. This may be due to the complex nature of the solution involving exponential, hyperbolic or cylindrical functions and complex integrals.

Burgers’ equation, the Burgers–Huxley equation, the Burgers–Fisher equation, the Kuramoto–Sivashinsky equation, the Fisher–Kolmogorov equation, the Zakharov–Kuznetsov–Burgers (ZKB) equation, the Hasegawa–Mima equation, the Korteweg–de Vries (KdV) equation, the KdV–Burgers equation, the Kolmogrov–Petrovskii–Piskunov (KPP) equation, the Telegraph equation, the Benjamin–Bona–Mohany–Burgers (BBMB) equation, the regularized long-wave (RLW) equation, the modified regularized long-wave (MRLW) equation, etc., are some of the equations that are widely used by scientists to explain different physical phenomena. Burgers’ equation has been solved analytically by the homotopy analysis method (HAM) [4] and homotopy perturbation method [5], and numerically by finite element [6], finite difference [7,8] and spline collocation [9,10,11] approaches. The Burgers–Huxley and Burgers–Fisher equations have also been solved analytically by the homotopy perturbation method [12,13] and numerically using cubic B-splines and quintic Hermite splines [14,15], respectively. The Newell–Whitehead–Seagel and FitzHugh–Nagumo equations have been solved analytically by the homotopy perturbation method, homotopy analysis method and integral transforms [16,17,18] and numerically by the differential quadrature method and cubic Hermite splines [19,20,21]. The ZKB equation has been solved analytically by the homotopy perturbation method as well as numerically by using the finite difference technique [22].

Among all these equations, the modified Burgers’ equation has wide applications in fluid dynamics, aeronautics, astronomy, eddies, ocean engineering, gas dynamics, traffic flow and non-linear acoustics, etc. [23,24]. However, most of the researchers have solved the Burgers’ equation for constant coefficients, ignoring the dissipation effect. The study of Burgers’ equation with variable coefficients is an arduous task as, in this situation, both the dispersion and dissipation factors go side by side, and sometimes the damping effect also interferes. Hence, to study the simultaneous effects of dispersion, dissipation and damping, the conformable damped Burgers’ equation has been studied. The time derivative has been taken to be of conformable type.

Many real-world problems are characterized by fractional-order differential equations as compared to integer-order differential equations. The exact solutions of these models are often difficult to find and as such, approximate methods are often more convenient for studying the behavior of these equations. Burgers’ equation has been studied by a number of researchers such as [25,26,27,28]. Analytic solutions of fractional Burgers’ equation have been pursued by different researchers as well. A new approximate method called the approximate analytical method (AAM) has been applied by [29] to solve the time-fractional damped Burgers’ equation. An iterative method has been applied by [30] to solve the fractional model of damped Burgers’ equation associated with the Caputo–Fabrizio fractional derivative. Approximate solutions of time-fractional Burgers’ equation and time-fractional Sharma–Tasso–Olver equations have been obtained using a new iterative method [31]. Also, by using the variational iteration method (VIM), Abd AL-Hussein et al. [32] obtained the approximate solution of fractional Burgers’ equation alongside its statistical concepts.

The present work is organized into nine sections. Section 1 contains the introduction of progressive non-linear systems and a brief literature survey on approximation techniques. The preliminaries of conformable derivative and its association with damped Burgers’ equation is included in Section 2. Section 3 elaborates on the conformable implicit finite difference method for time discretization. Section 4 describes the cubic Hermite collocation method, and its application is demonstrated in Section 5. Stability and error analysis of the proposed scheme are discussed in Section 6 and Section 7, respectively. Section 8 discusses the performance of the obtained results. Conclusions are drawn in Section 9.

2. Formulation of Conformable Damped Burgers’ Equation

Conformable derivatives can be considered as the special case of Hausdorff fractal derivative and the beta fractional derivative. Regardless of the argument of whether the conformable derivative is a fractional derivative or not, its fractional exponent allows it to regulate the order of the derivative. This local fractional derivative or conformable derivative can be viewed as a natural extension of classical derivative, satisfying most of the properties of classical calculus [33,34,35]. In recent times, fractional-order derivatives have been utilized in a wide range of applications in mathematical and applied sciences such as electricity, signal and image processing, biology, chemistry, economics and other areas [36,37,38,39,40,41]. Conformable derivatives have received prominent attention within the research community due to their suitability in describing real-life phenomena. The enhancement of integer-order derivatives to fractional-order derivatives allows one to control the rate of change to have a better view of evolving dynamics. This property enables one to recognize how past states have affected the present and future states of the concerned system. Therefore, the analytic and numerical analysis of conformable models enables one to perceive the macroscopic and microscopic mechanisms of the system.

Differential equations involving a differential operator in the conformable sense often result in singular equations. To manage this singularity, various advancements in existing methods have been proposed by researchers such as the conformable finite element method [42], the conformable Euler finite difference method [43], the conformable integral transforms, non-polynomial spline methods [44,45,46] and many more. With the motivation of evolving research for conformable differential equations, the present study proposes a new hybrid technique that combines an implicit finite difference formulation for the time direction and cubic Hermite splines for spatial interpolation.

Definition 1

([47]). For a function $u:[0,\infty )\to \mathbb{R}$, the αth conformable derivative is defined as

$$T^{\alpha }\left(u\right)\left(t\right)=\underset{\omega \to 0}{lim}{\displaystyle \frac{u(t+\omega t^{1-\alpha })-u\left(t\right)}{\omega }},$$

(1)

for all $\omega >0,t\ge 0,\alpha \in (0,1]$. If u is α-differentiable in some interval $(0,a)$, $a\ge 0$ and ${lim}_{t\to 0^{+}}{u}^{\left(\alpha \right)}\left(t\right)$ exist, then we define

$$u^{\left(\alpha \right)}\left(0\right)=\underset{t\to 0^{+}}{lim}{u}^{\left(\alpha \right)}\left(t\right).$$

Theorem 1

([47]). Let the function $u:[0,\infty )\to \mathbb{R}$ be α-differentiable at $t_0\ge 0$, $\alpha \in (0,1]$; then, u is continuous at $t_0$.

Theorem 2

([47]). Let $\alpha \in (0,1]$ and u be α-differentiable at a point $t\ge 0$; then, $T^{\alpha }\left(u\right)\left(t\right)=t^{1-\alpha }{\displaystyle \frac{du}{dt}}$.

From the above property of the conformable derivative and the definition of the conformable derivative, it is clear that it is an extension of the classical derivative with local properties and does not have a memory effect. It can be viewed as a special case of the classical derivative in the following form:

$${\displaystyle \frac{1}{2-\alpha }}{\displaystyle \frac{d}{dt}}u\left(t^{2-\alpha }\right)=t^{1-\alpha }{\displaystyle \frac{du}{dt}}.$$

(2)

Most of the properties of conformable fractional derivative are derived from the classical derivative. It helps to study the system having the effect of time being multiplied to its derivative.

The conformable damped Burgers’ equation [31,35] can be presented as follows:

$$T_t^{\alpha }u+\beta \left(t\right)u{u}_{\zeta }-\epsilon \left(t\right)u_{\zeta \zeta }=\lambda \left(t\right)u,0<\alpha \le 1,(\zeta ,t)\in (a,b)\times (0,T],$$

(3)

subject to the Dirichlet boundary conditions at $\zeta =a$ and $\zeta =b$, and the initial condition is as follows:

$$u(\zeta ,0)=\psi \left(\zeta \right),$$

where $u(\zeta ,t)$ is a function of two variables $\zeta$ and t indicating displacement; $\zeta$ is a space variable; t is a time variable; $\epsilon \left(t\right)$ is variable viscosity; $\beta \left(t\right)$ is a variable coefficient of non-linear advection term; and $\lambda \left(t\right)$ is a coefficient of damping term, defined on $[0,T]$ and is continuous on it. $T_t^{\alpha }u$ is the conformable derivative of the function $u(\zeta ,t)$, and $\alpha$ is the order of conformable derivative.

After implementing Theorem 2 on Equation (3) and rearranging the terms, the following equation is obtained:

$$t^{1-\alpha }{\displaystyle \frac{\partial u}{\partial t}}=\epsilon \left(t\right)u_{\zeta \zeta }-\beta \left(t\right)u{u}_{\zeta }+\lambda \left(t\right)u,0<\alpha \le 1,(\zeta ,t)\in (a,b)\times (0,T],$$

(4)

Now, Equation (4) is the required conformable damped Burgers’ equation.

3. Time Discretization

The proposed conformable damped Burgers’ equation will be discretized in the time direction first and then in the space direction. To implement the proposed conformable implicit finite difference scheme in the time direction of the proposed equation, we define the uniform partitioning in the time direction as follows:

$${\Pi }_T:0=t_0\le t_1\le \cdots \le t_M=T\mathrm{with}\Delta t=t_j-t_{j-1};j=1,\dots ,M.$$

The proposed finite difference technique is applied on Equation (4) and the following equation is obtained:

$$\begin{array}{c}\hfill {\displaystyle \frac{u(\zeta ,t_j)-u(\zeta ,t_{j-1})}{{(\Delta t)}^{\alpha }}}={\epsilon }_j{u}_{\zeta \zeta }(\zeta ,t_j)-{\beta }_{j}u(\zeta ,t_j)u_{\zeta }(\zeta ,t_j)+{\lambda }_{j}u(\zeta ,t_j);j=1,2,\dots ,M.\end{array}$$

(5)

For the sake of convenience, we write $u(\zeta ,t_j)=u^j$, and Equation (5) can be written as

$${\displaystyle \frac{u^j-u^{j-1}}{{(\Delta t)}^{\alpha }}}={\epsilon }_j{u}_{\zeta \zeta }^j-{\beta }_j{\left(u{u}_{\zeta }\right)}^j+{\lambda }_j{u}^j;j=1,2,\dots ,M.$$

(6)

The non-linear advection term ${\left(u{u}_{\zeta }\right)}^j$ is quasi-linearized using a formula given by [48], which is explicit in nature. This form is linearly convergent but computationally convenient to implement.

$${\left(u{u}_{\zeta }\right)}^j={\displaystyle \frac{1}{3}}\left(u^j{u}_{\zeta }^{j-1}+2u^{j-1}{u}_{\zeta }^j\right).$$

(7)

After applying the discretization on the non-linear term, the following equation is obtained:

$${\displaystyle \frac{u^j-u^{j-1}}{{(\Delta t)}^{\alpha }}}={\epsilon }_j{u}_{\zeta \zeta }^j-{\displaystyle \frac{{\beta }_j}{3}}\left(u^j{u}_{\zeta }^{j-1}+2u^{j-1}{u}_{\zeta }^j\right)+{\lambda }_j{u}^j;j=1,2,\dots ,M.$$

(8)

After rearranging the terms, the following equation is obtained:

$$\left(1+{\displaystyle \frac{{\beta }_j{(\Delta t)}^{\alpha }}{3}}{u}_{\zeta }^{j-1}-{\lambda }_j{(\Delta t)}^{\alpha }\right)u^j-{\epsilon }_j{(\Delta t)}^{\alpha }{u}_{\zeta \zeta }^j+{\displaystyle \frac{2{\beta }_j{(\Delta t)}^{\alpha }}{3}}{u}^{j-1}{u}_{\zeta }^j=u^{j-1};j=1,2,\dots ,M,$$

(9)

where $u(\zeta ,t_0)=\psi \left(\zeta \right)$.

Lemma 1.

Let $u(\zeta ,t_j)$ be the exact solution at $t=t_j$ and $u^j$ be the approximate finite difference solution at the jth iteration; then, the order of convergence of the proposed technique is ${(\Delta t)}^{\alpha }$ for $0\le \alpha <1$.

$$e^j\le C{(\Delta t)}^{\alpha };0\le \alpha <1$$

(10)

where $e^j=u(\zeta ,t_j)-u^j$; $j=1,2,\dots ,M$.

After the discretization of the proposed equation in the time direction, orthogonal collocation is applied with cubic Hermite splines as the base function. This reduces the non-linear equation into a system of algebraic equations, which are then solved using MATLAB 2021.

In operator form, Equation (9) can be written as follows:

$$\mathfrak{T}\left(u^j\right)={\mathfrak{G}}_{\Delta t},$$

(11)

where $\mathfrak{T}\left(u^j\right)=\left(1+{\displaystyle \frac{{\beta }_j{(\Delta t)}^{\alpha }}{3}}{u}_{\zeta }^{j-1}-{\lambda }_j{(\Delta t)}^{\alpha }\right)u^j-{\epsilon }_j{(\Delta t)}^{\alpha }{u}_{\zeta \zeta }^j+{\displaystyle \frac{2{\beta }_j{(\Delta t)}^{\alpha }}{3}}{u}^{j-1}{u}_{\zeta }^j$ and ${\mathfrak{G}}_{\Delta t}=u^{j-1}$.

4. Cubic Hermite Collocation Method

Weighted residual methods are tools for solving the boundary value problems in which the residual is set orthogonal to the weight function. Depending upon the weight functions, the Galerkin method, collocation method, least square method, etc., are defined. Collocation techniques are one of the weighted residual methods. In this technique, the residual is set equal to zero at the collocation points. The choice of trial function makes the collocation technique different from other residual methods. Various investigators have applied different base functions such as Legendre polynomials [37,49], Chebyshev polynomials [50,51,52], Bessel polynomials [53], Lagrangian interpolating polynomials [54,55], Hermite splines [21,56,57], B-splines [58,59,60,61], cubic Hermite B-splines [62,63], exponential splines [64,65], radial basis functions [66], etc., to study the behavior of non-linear partial differential equations.

In the present study, Hermite splines of the order of three are considered to approximate the trial function in the spatial direction. Hermite interpolating polynomials passes through a set of points and their tangents at some point. These cubic splines achieve high accuracy with fewer collocation points than orthogonal collocation with Lagrangian basis for stiff systems of equations [56]. These cubic splines interpolate the boundary as well and, therefore, it handles the boundary conditions easily and preserves certain properties of the original equation, such as symmetry and conservation laws.

Cubic Hermite splines interpolate the function as well as its first-order derivative and thus convert the non-linear equation into a system of ordinary differential equations. Also, these splines interpolate the derivative of the function; therefore, there is no requirement of fictitious nodes as in B-splines. Hermite splines inherit the continuity property at the node points [62] and produce a diagonal system that allows for easy storage in digital computers, yields high computer speed and reduces the cost of storage. The details of cubic Hermite splines are given in [21].

Consider an interval $\mathcal{I}=(a,b)$, and we define the uniform partitioning ${\Pi }_{\zeta }$ of $(a,b)$ such that

$${\Pi }_{\zeta }:a={\zeta }_0\le {\zeta }_1\le {\zeta }_2\le \dots \le {\zeta }_n=b.$$

Let $P_3$ be the space of all cubic polynomials defined on $[{\zeta }_i,{\zeta }_{i+1}]$. Let $\nu$ be the continuously differentiable function and cubic Hermite approximation of u defined on $\overline{\mathcal{I}}$.

Define the space of all continuously differentiable functions on $[{\zeta }_i,{\zeta }_{i+1}]$ by ${\mathbb{M}}_{\rho }$ such that

$${\mathbb{M}}_{\rho }=\{\nu |\nu \in P_{3}on[{\zeta }_i,{\zeta }_{i+1}],i=0,1,2,\dots ,n-1\}.$$

(12)

If boundary conditions are of homogeneous type, then the dimension of ${\mathbb{M}}_{\rho }=2n$, and if the boundary conditions are of non-homogeneous type, then the dimension of ${\mathbb{M}}_{\rho }=2n+2$ [67].

The cubic Hermite polynomials for $\rho =1$ have been defined by ${\mathfrak{P}}_{\rho }\left(\zeta \right)$ and ${\mathfrak{Q}}_{\rho }\left(\zeta \right)$ and are presented below.

$${\mathfrak{P}}_{\rho }\left(\zeta \right)=\left\{\begin{array}{cc}3{\left({\displaystyle \frac{\zeta -{\zeta }_{\rho -1}}{{\zeta }_{\rho }-{\zeta }_{\rho -1}}}\right)}^2-2{\left({\displaystyle \frac{\zeta -{\zeta }_{\rho -1}}{{\zeta }_{\rho }-{\zeta }_{\rho -1}}}\right)}^3;\hfill & {\zeta }_{\rho -1}\le \zeta \le {\zeta }_{\rho }\hfill \\ 3{\left({\displaystyle \frac{{\zeta }_{\rho +1}-\zeta }{{\zeta }_{\rho +1}-{\zeta }_{\rho }}}\right)}^2-2{\left({\displaystyle \frac{{\zeta }_{\rho +1}-\zeta }{{\zeta }_{\rho +1}-{\zeta }_{\rho }}}\right)}^3;\hfill & {\zeta }_{\rho }\le \zeta \le {\zeta }_{\rho +1}\hfill \\ 0\hfill & \mathrm{elsewhere}\hfill \end{array}\right.$$

(13)

$${\mathfrak{Q}}_{\rho }\left(\zeta \right)=\left\{\begin{array}{cc}-{\displaystyle \frac{{(\zeta -{\zeta }_{\rho -1})}^2}{{\zeta }_{\rho }-{\zeta }_{\rho -1}}}+{\displaystyle \frac{{(\zeta -{\zeta }_{\rho -1})}^3}{{({\zeta }_{\rho }-{\zeta }_{\rho -1})}^2}};\hfill & {\zeta }_{\rho -1}\le \zeta \le {\zeta }_{\rho }\hfill \\ {\displaystyle \frac{{({\zeta }_{\rho +1}-\zeta )}^2}{{\zeta }_{\rho +1}-{\zeta }_{\rho }}}-{\displaystyle \frac{{({\zeta }_{\rho +1}-\zeta )}^3}{{({\zeta }_{\rho +1}-{\zeta }_{\rho })}^2}};\hfill & {\zeta }_{\rho }\le \zeta \le {\zeta }_{\rho +1}\hfill \\ 0\hfill & \mathrm{elsewhere}\hfill \end{array}\right.$$

(14)

The piecewise cubic Hermite splines possess the following properties:

$${\mathfrak{P}}_{\rho }\left({\zeta }_i\right)={\delta }_{\rho i};{\mathfrak{P}}_{\rho }^{\prime }\left({\zeta }_i\right)=0;{\mathfrak{Q}}_{\rho }\left({\zeta }_i\right)=0;{\mathfrak{Q}}_{\rho }^{\prime }\left({\zeta }_i\right)={\delta }_{\rho i}.$$

In orthogonal collocation, the global variable $\zeta$ varies in the pth element, where $p=1,2,\dots ,n$. A new variable $\xi$ is introduced in the pth element as $\xi =(\zeta -{\zeta }_p)/h_p$ in such a way that as $\zeta$ varies from ${\zeta }_p$ to ${\zeta }_{p+1}$, $\xi$ varies from 0 to 1, where $h_p=({\zeta }_{p+1}-{\zeta }_p)$. Orthogonal collocation is then applied on the local variable $\xi$. It converts the polynomials defined in Equations (13) and (14) into the form defined in

Table 1. The cubic Hermite splines and its first and second derivatives.

i	${\mathit{H}}_{\mathit{i}}$	${\mathit{H}}_{\mathit{i}}^{\prime }$	${\mathit{H}}_{\mathit{i}}^{″}$
1	$1-3{\xi }^2-2{\xi }^3$	$6\xi (\xi -1)$	12$\xi -6$
2	$h{\xi }^2(\xi -1)$	$h\xi (3\xi -2)$	$2h(3\xi -1)$
3	$\xi {(\xi -1)}^2$	$2\xi (\xi -1)+{(\xi -1)}^2$	$6\xi -4$
4	$h{\xi }^2(3-2\xi )$	$6h\xi (1-\xi )$	$6h(1-2\xi )$

.

The grid points ${\zeta }_k$ are often called the ’knots’ of the piecewise polynomial since they are the points where the polynomials are “tied together”. The Hermite polynomials do not require the extra condition to make the first derivative continuous. This helps to reduce the number of equations by $(n-1)$, where n is the number of elements.

The next step in the application of cubic Hermite splines on the proposed conformable damped Burgers’ equation is the choice of collocation points. In this study, the zeros of shifted Legendre polynomials have been taken as collocation points. The detailed study of collocation points is given in [21] and the details of Legendre polynomials are given in [68,69].

5. Application of Cubic Hermite Splines Collocation

To fully discretize the system of equations within each sub-interval $[{\zeta }_i,{\zeta }_{i+1}]$, a one-to-one correspondence between $[{\zeta }_{i-1},{\zeta }_i]$ and [0, 1] is introduced, and then the collocation technique is applied with cubic Hermite splines as the basis. A linear combination of four cubic Hermite splines with time-dependent coefficients is applied in the pth element as follows:

$$u(\xi ,t_j)=\sum _{i=1}^4{H}_i\left(\xi \right){\chi }_i^{\left(j\right)p};j=0,1,2,\dots ,M;p=1,2,\dots ,n,$$

(15)

where ${\chi }_i^j$’s are the unknown coefficients which are to be determined. A system of equations is obtained in the pth element after the substitution of Equation (15) in Equation (9) in the following form.

$$\begin{array}{c}\hfill \left(1+{\displaystyle \frac{{\beta }_j{(\Delta t)}^{\alpha }}{3h}}\sum _{i=1}^4{H}_i^{\prime }\left(\xi \right){\chi }_i^{(j-1)p}-{\lambda }_j{(\Delta t)}^{\alpha }\right)\sum _{i=1}^4{H}_i\left(\xi \right){\chi }_i^{\left(j\right)p}-{\displaystyle \frac{{\epsilon }_j{(\Delta t)}^{\alpha }t}{h^2}}\sum _{i=1}^4{H}_i^{″}\left(\xi \right){\chi }_i^{\left(j\right)p}\\ \hfill +{\displaystyle \frac{2{\beta }_j{(\Delta t)}^{\alpha }}{3h}}\sum _{i=1}^4{H}_i\left(\xi \right){\chi }_i^{(j-1)p}\sum _{i=1}^4{H}_i^{\prime }\left(\xi \right){\chi }_i^{\left(j\right)p}=\sum _{i=1}^4{H}_i\left(\xi \right){\chi }_i^{(j-1)p};j=1,2,\dots ,M;p=1,2,\dots ,n.\end{array}$$

(16)

At the kth collocation point, Equation (16) reduces to the following form in the pth element:

$$\begin{array}{c}\hfill \left(1+{\displaystyle \frac{{\beta }_j{(\Delta t)}^{\alpha }}{3h}}\sum _{i=1}^4{H}_{ki}^{\prime }{\chi }_i^{(j-1)p}-{\lambda }_j{(\Delta t)}^{\alpha }\right)\sum _{i=1}^4{H}_{ki}{\chi }_i^{\left(j\right)p}-{\displaystyle \frac{{\epsilon }_j{(\Delta t)}^{\alpha }}{h^2}}\sum _{i=1}^4{H}_{ki}^{″}{\chi }_i^{\left(j\right)p}+{\displaystyle \frac{2{\beta }_j{(\Delta t)}^{\alpha }}{3h}}\sum _{i=1}^4{H}_{ki}{\chi }_i^{(j-1)p}\sum _{i=1}^4{H}_{ki}^{\prime }{\chi }_i^{\left(j\right)p}\\ \hfill =\sum _{i=1}^4{H}_{ki}{\chi }_i^{(j-1)p};j=1,2,\dots ,M;p=1,2,\dots ,n,\end{array}$$

(17)

where $H_{ki}$, $H_{ki}^{\prime }$ and $H_{ki}^{″}$ are Hermite spline interpolants, and first- and second-order derivatives of Hermite spline interpolants at the kth collocation point, respectively. The system of equations is converted into the matrix form as given below:

$$\mathbb{A}{\chi }^j=\mathbb{H}{\chi }^{j-1}+F^{j-1},$$

(18)

where

$$\mathbb{A}=\left\{\begin{array}{c}\left(1+{\displaystyle \frac{{\beta }_j{(\Delta t)}^{\alpha }}{3h}}{\sum }_{i=2}^4{H}_{ki}^{\prime }{\chi }_i^{(j-1)p}-{\lambda }_j{(\Delta t)}^{\alpha }\right)H_{ki}-{\displaystyle \frac{{\epsilon }_j{(\Delta t)}^{\alpha }}{h^2}}{H}_{ki}^{″}+{\displaystyle \frac{2{\beta }_j{(\Delta t)}^{\alpha }}{3h}}{\sum }_{i=2}^4{H}_{ki}{\chi }_i^{(j-1)p}{H}_{ki}^{\prime };\hfill \\ \mathrm{for}{1}^{st}\mathrm{element}i=2,3,4\hfill \\ \left(1+{\displaystyle \frac{{\beta }_j{(\Delta t)}^{\alpha }}{3h}}{\sum }_{i=1}^4{H}_{ki}^{\prime }{\chi }_i^{(j-1)p}-{\lambda }_j{(\Delta t)}^{\alpha }\right)H_{ki}-{\displaystyle \frac{{\epsilon }_j{(\Delta t)}^{\alpha }}{h^2}}{H}_{ki}^{″}+{\displaystyle \frac{2{\beta }_j{(\Delta t)}^{\alpha }}{3h}}{\sum }_{i=1}^4{H}_{ki}{\chi }_i^{(j-1)p}{H}_{ki}^{\prime };\hfill \\ \mathrm{for}2,3,\dots ,(n-1)\mathrm{th}\mathrm{element}\hfill \\ \left(1+{\displaystyle \frac{{\beta }_j{(\Delta t)}^{\alpha }}{3h}}{\sum }_{i=1}^3{H}_{ki}^{\prime }{\chi }_i^{(j-1)p}-{\lambda }_j{(\Delta t)}^{\alpha }\right)H_{ki}-{\displaystyle \frac{{\epsilon }_j{(\Delta t)}^{\alpha }}{h^2}}{H}_{ki}^{″}+{\displaystyle \frac{2{\beta }_j{(\Delta t)}^{\alpha }}{3h}}{\sum }_{i=1}^3{H}_{ki}{\chi }_i^{(j-1)p}{H}_{ki}^{\prime };\hfill \\ \mathrm{for}{n}^{th}\mathrm{element}i=1,2,3\hfill \end{array}\right.$$

and

$$\mathbb{H}=\left\{\begin{array}{c}{H}_{ki};\mathrm{for}{1}^{st}\mathrm{element}i=1,2,3\hfill \\ H_{ki};\mathrm{for}2,3,\dots ,{(n-1)}^{th}\mathrm{element}\hfill \\ H_{ki};\mathrm{for}{n}^{th}\mathrm{element}i=1,2,3\hfill \end{array}\right.$$

The matrices $\mathbb{A}$ and $\mathbb{H}$ are almost quad-diagonal dominant and are of order $2n\times 2n$. ${\chi }^j={[{\chi }_1^j,{\chi }_2^j,⋮,\dots ,{\chi }_{2n}^j]}^T$; the matrix $F^j$ is obtained from the non-homogeneous boundary conditions $u(a,t)=g_1\left(t\right)$ and $u(b,t)=g_2\left(t\right)$ if it exists; otherwise, this will be the zero vector. The resulting system of matrix equations is solved iteratively. In compact form, Equation (18) can be written as follows:

$$\mathbb{A}\chi =\mathcal{H},$$

(19)

where $\mathcal{H}=\mathbb{H}{\chi }^{j-1}+F^{j-1}$.

6. Stability Analysis

The development of every numerical technique depends upon its stability. To discuss the stability of the proposed finite difference scheme, the local bound of $u_{\zeta }$ has been considered for all $(\zeta ,t)$. Let $3\eta$ be the local bound of $u_{\zeta }(\zeta ,t)$ for all $(\zeta ,t)$ in the domain $[a,b]\times [0,T]$. Let ${\gamma }_j=1+(\eta {\beta }_j-{\lambda }_j){(\Delta t)}^{\alpha }$, ${\overline{\epsilon }}_j={\epsilon }_j{(\Delta t)}^{\alpha }$ and ${\overline{\gamma }}_j=1-2\eta {\beta }_j{(\Delta t)}^{\alpha }$. After substituting the local error bound for $u_{\zeta }$, Equation (17) takes the following form:

$$\begin{array}{c}\hfill {\gamma }_j{u}^j-{\overline{\epsilon }}_j{u}_{\zeta \zeta }^j={\overline{\gamma }}_j{u}^{j-1}.j=1,2,\dots ,M;p=1,2,\dots ,n,\end{array}$$

(20)

Let $u(\zeta ,t)$ be the solution of Equation (3) defined on the domain $\Omega$; then, the inner product of two functions u and v over the domain $\Omega$ is defined as follows:

$${\langle u,v\rangle }_{\Omega }={\int }_{\Omega }uv,$$

$${\parallel u\parallel }_{L_2(\Omega )}=\sqrt{{\int }_{\Omega }{u}^{2}d\Omega },$$

$${\langle u,u\rangle }_{\Omega }={\parallel u\parallel }_{L_2(\Omega )}^2,$$

$${\parallel u\parallel }_{L_{\infty (\Omega )}}={max}_{(\zeta ,t)\in \Omega }|u(\zeta ,t)|.$$

Theorem 3.

Cauchy–Schwarz Inequality: Let u and v be any two functions defined on Ω; then,

$$\mid \langle u,v\rangle {\mid }_{\Omega }\le {\parallel u\parallel }_{L_2(\Omega )}{\parallel v\parallel }_{L_2(\Omega )}$$

Lemma 2

([48]). Discrete Sobolev Inequality: For any discrete function $u_p$, there exists constants $K_1$ and $K_2$ such that

$$\parallel u^j{\parallel }_{\infty }\le K_1\parallel u^j\parallel +K_2\parallel u_{\zeta }^j\parallel .$$

(21)

Lemma 3.

Let us suppose that $u_0\in$ $H_2[0,\Omega ]$. Then, there exists a constant C such that

$$\parallel u^j{\parallel }_2\le C.$$

(22)

Proof. 

Consider Equation (20)

$$\begin{array}{c}\hfill {\gamma }_j{u}^j-{\overline{\epsilon }}_j{u}_{\zeta \zeta }^j={\overline{\gamma }}_j{u}^{j-1}.j=1,2,\dots ,M;p=1,2,\dots ,n,\end{array}$$

(23)

Taking the inner product with $u^j$ on both sides,

$${\gamma }_j\langle u^j,u^j\rangle -{\tilde{\epsilon }}_j\langle u_{\zeta \zeta }^j,u^j\rangle \le {\overline{\gamma }}_j\langle u^{j-1},u^j\rangle ,j=1,2,\dots ,M.$$

(24)

Using the result $\langle u_{\zeta \zeta },u\rangle \le 0$, the following inequality is obtained:

$${\gamma }_j\langle u^j,u^j\rangle \le {\overline{\gamma }}_j\langle u^{j-1},u^j\rangle ,j=1,2,\dots ,M.$$

(25)

Using the Cauchy–Schwarz inequality implies the following:

$$\mid {\gamma }_j\mid \parallel u^j{\parallel }_2^2\le \mid {\overline{\gamma }}_j\mid \parallel u^{j-1}{\parallel }_2{\parallel u^j\parallel }_2,j=1,2,\dots ,M.$$

(26)

Rearranging the terms yields the following:

$$\parallel u^j{\parallel }_2\le {\displaystyle \frac{\mid {\overline{\gamma }}_j\mid }{\mid {\gamma }_j\mid }}{\parallel u^{j-1}\parallel }_2,j=1,2,\dots ,M.$$

(27)

For $j=1$, Equation (27) is as follows:

$$\parallel u^1{\parallel }_2\le {\displaystyle \frac{\mid {\overline{\gamma }}_1\mid }{\mid {\gamma }_1\mid }}{\parallel u^0\parallel }_2.$$

(28)

Using the Principle of induction and Equation (28), the upper bound on $u^j$ can be defined as follows:

$$\parallel u^j{\parallel }_2\le C{\parallel u^0\parallel }_2,j=1,2,\dots ,M,$$

(29)

where C is the generic constant. Hence, the upper bound on $u^j$ is independent of the jth iteration, which shows that the proposed implicit scheme is unconditionally stable. □

Theorem 4.

Let us suppose that $u\in$ $H_2[(0,T),\Omega ]$ is the exact solution of Equation (3) and $u^j$ is the approximate finite difference solution at the jth iteration. Let $e^j=u(\zeta ,t_j)-u^j$ be the error estimate at the jth iteration. Then,

$$\parallel e^j{\parallel }_2\le C{\parallel e^0\parallel }_2.$$

(30)

Proof. 

$u(\zeta ,t_j)$ is the exact solution of Equation (3) and $u^j$ is the approximate finite difference solution of Equation (3) at $t=t_j$, therefore, Equation (20) can be written as follows:

$${\gamma }_{j}u(\zeta ,t_j)\le {\tilde{\epsilon }}_j{u}_{\zeta \zeta }(\zeta ,t_j)+{\overline{\gamma }}_{j}u(\zeta ,t_{j-1}),j=1,2,\dots ,M.$$

(31)

$${\gamma }_j{u}^j-{\tilde{\epsilon }}_j{u}_{\zeta \zeta }^j\le {\overline{\gamma }}_j{u}^{j-1},j=1,2,\dots ,M.$$

(32)

Substituting Equation (32) from Equation (31), the following equation is obtained:

$${\gamma }_j{e}^j-{\tilde{\epsilon }}_j{e}_{\zeta \zeta }^j\le {\overline{\gamma }}_j{e}^{j-1},j=1,2,\dots ,M.$$

(33)

Now, taking the inner product with $e^j$ on both sides and proceeding with Lemma 3 and the principle of induction, the following inequality is obtained:

$$\parallel e^j{\parallel }_2\le C{\parallel e^0\parallel }_2,j=0,1,2,\dots ,M-1,$$

(34)

where C is the generic constant. □

Definition 2

([48]). A numerical scheme is said to be globally stable if there exists a constant C

$$\parallel e^j\parallel \le C\parallel e^0\parallel ,$$

(35)

where $e^j=u(\zeta ,t_j)-u^j$, for $j=1,2,\dots ,M$, such that $u(\zeta ,t_j)$ be the exact solution at the $t=t_j$ and $u^j$ be the approximate finite difference solution at the jth iteration.

From Theorem 4 and the definition of global stability, it is clear that the proposed scheme is unconditionally globally stable.

7. Error Analysis

It is fundamental to perform an error analysis in the application of any numerical technique. This forms the basis for finding an approximate solution to the given equations. Numerical values obtained after the discretization are subjected to some form of error analysis. The cubic Hermite splines are bounded and can be applied on two-point boundary value problems [21,70].

Lemma 4

([71]). Let $u(\zeta ,t)$ be the solution of Equation (3) in $\overline{\Omega }$; then, the bound on the derivative of u w.r.t. ζ is given by the following:

$$\left|{\displaystyle \frac{{\partial }^{i}u}{\partial {\zeta }^i}}\right|\le C\forall (\zeta ,t)\in \overline{\Omega },i=0,1,2.$$

(36)

Lemma 5

([48]). Let $u(\zeta ,\tau )$ be the solution of Equation (3) in $\overline{\Omega }$ and $\mid u\mid \le C$; then, the bound on the derivative of u w.r.t. t is given by the following:

$$\left|{\displaystyle \frac{{\partial }^{j}u}{\partial t^j}}\right|\le C\forall (\zeta ,t)\in \overline{\Omega },0\le j<1.$$

(37)

Lemma 6

([71]). Let $u(\zeta ,\tau )$ be the solution of Equation (3) in $\overline{\Omega }$; then, the bound on the derivative of u w.r.t. ζ and t is given by the following:

$$\left|{\displaystyle \frac{{\partial }^{i+j}u}{\partial {\zeta }^i\partial t^j}}\right|\le C\forall (\zeta ,t)\in \overline{\Omega },i=0,1,2and0\le j<1.$$

(38)

Theorem 5

([72,73]). Let H $\in H_{\Delta \zeta }^4\left(\zeta \right)$ be the piecewise cubic Hermite splines over the sub-interval $[{\zeta }_i,{\zeta }_{i+1}]$. Let H approximate $\mathit{U}$ $\in C^4$[a,b], and then

$$\parallel {\mathit{H}}^{\left(s\right)}-{\mathit{U}}^{\left(s\right)}{\parallel }_{\infty }\le C{\mu }_s{h}^{4-s},s=0,1,2,3.$$

where ${\mathit{H}}^{\left(s\right)}$ is the sth-order derivative of $\mathit{H}$. The values of ${\mu }_s$ are given in [72,73].

Theorem 6.

Let $u\left(\zeta \right)$ and $\overline{u}\left(\zeta \right)$ be the exact solution and approximate solution, respectively, of Equation (3) in the space ${\mathbb{H}}^3$ of cubic Hermite splines such that $\overline{u}\left(\zeta \right)\in C^4[a,b]$. Then, the uniform error estimate is as follows:

$$\parallel u\left(\zeta \right)-\overline{u}{\left(\zeta \right)\parallel }_{\infty }\le C{h}^2\left({\mu }_0{h}^2\left(\eta |\overline{{\beta }^j}|+|\overline{{\lambda }^j}|\right)+{\mu }_{1}h|\overline{{\beta }^j}|+{\mu }_2|\overline{{\epsilon }^j}|\right).$$

(39)

Proof. 

Let $\mathbf{H}\left(\zeta \right)$ be the unique cubic Hermite spline interpolate of $u\left(\zeta \right)$. Let

$${\Pi }_{\zeta }=a={\zeta }_0\le {\zeta }_1\le ,\dots ,\le {\zeta }_n=b,$$

be the partition of $[a,b]$ with uniform step size h. Let $\overline{{\lambda }^j}=(1-{\lambda }^j{(\Delta t)}^{\alpha })$, $\overline{{\epsilon }^j}={\epsilon }^j{(\Delta t)}^{\alpha }$ and $\overline{{\beta }^j}={\beta }^j{(\Delta t)}^{\alpha }.$

Consider

$$\begin{array}{c}\hfill \parallel \mathfrak{T}\left(u\left({\zeta }_k\right)\right)-\mathfrak{T}\left(H\left({\zeta }_k\right)\right)){\parallel }_{\infty }=\parallel \overline{{\lambda }^j}(u\left({\zeta }_k\right)-H\left({\zeta }_k\right))-\overline{{\epsilon }^j}(u^{″}\left({\zeta }_k\right)-H^{″}\left({\zeta }_k\right))+\overline{{\beta }^j}(u\left({\zeta }_k\right)u^{\prime }\left({\zeta }_k\right)-H\left({\zeta }_k\right)H^{\prime }\left({\zeta }_k\right)){\parallel }_{\infty }\\ \hfill \le |\overline{{\lambda }^j}|\parallel u\left({\zeta }_k\right)-H\left({\zeta }_k\right){\parallel }_{\infty }+|\overline{{\epsilon }^j}|\parallel u^{″}\left({\zeta }_k\right)-H^{″}\left({\zeta }_k\right){\parallel }_{\infty }+|\overline{{\beta }^j}|\parallel u\left({\zeta }_k\right)u^{\prime }\left({\zeta }_k\right)-u^{\prime }\left({\zeta }_k\right)H\left({\zeta }_k\right)+u^{\prime }\left({\zeta }_k\right)H\left({\zeta }_k\right)-H\left({\zeta }_k\right)H^{\prime }\left({\zeta }_k\right){\parallel }_{\infty }\\ \hfill \le |\overline{{\lambda }^j}|C{\mu }_0{h}^4+|\overline{{\epsilon }^j}|C{\mu }_2{h}^2+|\overline{{\beta }^j}|\left(C{\mu }_0{h}^4\parallel u^{\prime }\left({\zeta }_k\right){\parallel }_{\infty }+C{\mu }_1{h}^3\parallel H\left({\zeta }_k\right){\parallel }_{\infty }\right)\\ \hfill \le C{h}^2\left({\mu }_0{h}^2\left(\eta |\overline{{\beta }^j}|+|\overline{{\lambda }^j}|\right)+{\mu }_{1}h|\overline{{\beta }^j}|+{\mu }_2|\overline{{\epsilon }^j}|\right),\end{array}$$

(40)

where $\eta =max|u_{\zeta }|$ for all $(\zeta ,t)$ lying in the domain $(a,b)\times (0,T)$, and the Hermite splines are bounded by unity. Using Theorem 5, the following inequality is obtained:

$$\parallel \mathfrak{T}\left(u\left({\zeta }_k\right)\right)-\mathfrak{T}\left(H\left({\zeta }_k\right)\right)){\parallel }_{\infty }\le C{h}^2\left({\mu }_0{h}^2\left(\eta |\overline{{\beta }^j}|+|\overline{{\lambda }^j}|\right)+{\mu }_{1}h|\overline{{\beta }^j}|+{\mu }_2|\overline{{\epsilon }^j}|\right).$$

(41)

where C is the generic constant and the values of ${\mu }_0$ and ${\mu }_2$ are given in [72,73].

Now, consider

$$\begin{array}{c}\hfill \parallel \mathfrak{T}\left(\overline{u}\left({\zeta }_k\right)\right)-\mathfrak{T}\left(H\left({\zeta }_k\right)\right){\parallel }_{\infty }\hspace{0.17em}=\hspace{0.17em}\parallel {\mathfrak{G}}_{\Delta t}\left({\zeta }_k\right)-\mathfrak{T}\left(H\left({\zeta }_k\right)\right){\parallel }_{\infty }\\ \hfill \le \hspace{0.17em}\parallel \mathfrak{T}\left(u\left({\zeta }_k\right)\right)-\mathfrak{T}\left(H\left({\zeta }_k\right)\right){\parallel }_{\infty }\\ \hfill \le C{h}^2\left({\mu }_0{h}^2\left(\eta |\overline{{\beta }^j}|+|\overline{{\lambda }^j}|\right)+{\mu }_{1}h|\overline{{\beta }^j}|+{\mu }_2|\overline{{\epsilon }^j}|\right).\end{array}$$

(42)

Now, we use $\mathfrak{T}\left(\overline{u}\left({\zeta }_k\right)\right)={\mathfrak{G}}_{\Delta t}\left({\zeta }_k\right)$ and $\mathbb{A}\chi =\mathcal{H}$.

This implies $\mathbb{A}(\chi -\overline{\chi })=\mathcal{H}-\overline{\mathcal{H}}$, and using Equation (42),

$$\begin{array}{c}\hfill \parallel \mathcal{H}-\overline{\mathcal{H}}{\parallel }_{\infty }\hspace{0.17em}\le max\parallel {\mathfrak{G}}_{\Delta t}\left({\zeta }_k\right)-{\overline{\mathfrak{G}}}_{\Delta t}\left({\zeta }_k\right){\parallel }_{\infty },\\ \hfill \le max\parallel \mathfrak{T}\left(\overline{u}\left({\zeta }_k\right)\right)-\mathfrak{T}\left(H\left({\zeta }_k\right)\right){\parallel }_{\infty },\\ \hfill \le C{h}^2\left({\mu }_0{h}^2\left(\eta |\overline{{\beta }^j}|+|\overline{{\lambda }^j}|\right)+{\mu }_{1}h|\overline{{\beta }^j}|+{\mu }_2|\overline{{\epsilon }^j}|\right).\end{array}$$

(43)

Also, $\mathbb{A}$ is bounded being the matrix of collocation coefficients. Therefore, $\parallel {\mathbb{A}}^{-1}\parallel \hspace{0.17em}\le C$

$$\Rightarrow \parallel \chi -\overline{\chi }{\parallel }_{\infty }\le C{h}^2\left({\mu }_0{h}^2\left(\eta |\overline{{\beta }^j}|+|\overline{{\lambda }^j}|\right)+{\mu }_{1}h|\overline{{\beta }^j}|+{\mu }_2|\overline{{\epsilon }^j}|\right).$$

(44)

Now, $\overline{u}\left(\zeta \right)-H\left(\zeta \right)={\sum }_{i=1}^4{H}_i\left(\zeta \right)({\chi }_i-{\overline{\chi }}_i)$

$$\begin{array}{c}\hfill \Rightarrow \parallel \overline{u}\left(\zeta \right)-H\left(\zeta \right){\parallel }_{\infty }=\parallel \sum _{i=1}^4{H}_i\left(\zeta \right)({\chi }_i-{\overline{\chi }}_i){\parallel }_{\infty },\\ \hfill \le \sum _{i=1}^4|H_i\left(\zeta \right)|\parallel ({\chi }_i-{\overline{\chi }}_i){\parallel }_{\infty },\\ \hfill \le C{h}^2\left({\mu }_0{h}^2\left(\eta |\overline{{\beta }^j}|+|\overline{{\lambda }^j}|\right)+{\mu }_{1}h|\overline{{\beta }^j}|+{\mu }_2|\overline{{\epsilon }^j}|\right).\end{array}$$

(45)

According to [72],

$$\parallel u\left(\zeta \right)-H\left(\zeta \right){\parallel }_{\infty }=C{\mu }_0{h}^4.$$

(46)

Now, using the triangle inequality,

$$\begin{array}{cc}\hfill \parallel u\left(\zeta \right)-\overline{u}\left(\zeta \right){\parallel }_{\infty }& \le \hspace{0.17em}\parallel u\left(\zeta \right)-H\left(\zeta \right){\parallel }_{\infty }+\parallel \overline{u}\left(\zeta \right)-H\left(\zeta \right){\parallel }_{\infty }\hfill \\ \hfill & \le C{\mu }_0{h}^4+C{h}^2\left({\mu }_0{h}^2\left(\eta |\overline{{\beta }^j}|+|\overline{{\lambda }^j}|\right)+{\mu }_{1}h|\overline{{\beta }^j}|+{\mu }_2|\overline{{\epsilon }^j}|\right)\hfill \\ \hfill & \le C{h}^2\left({\mu }_0{h}^2\left(1+\eta |\overline{{\beta }^j}|+|\overline{{\lambda }^j}|\right)+{\mu }_{1}h|\overline{{\beta }^j}|+{\mu }_2|\overline{{\epsilon }^j}|\right).\hfill \end{array}$$

(47)

□

Theorem 7.

Let $u(\zeta ,t)$ and $\overline{u}(\zeta ,t)\in C^4[a,b]$ be the exact solution and approximate solution, respectively, of Equation (3). Then, the following uniform error estimate is obtained:

$$\parallel u-\overline{u}{\parallel }_{\infty }\hspace{0.17em}\le C\left(h^2\left({\mu }_0{h}^2\left(\eta |\overline{{\beta }^j}|+|\overline{{\lambda }^j}|\right)+{\mu }_{1}h|\overline{{\beta }^j}|+{\mu }_2|\overline{{\epsilon }^j}|\right)+{(\Delta t)}^{\alpha }\right).$$

(48)

Proof. 

Equation (48) is the direct consequence of Theorem 6 and Lemma 1. □

8. Numerical Results and Discussion

The demonstration of any technique is incomplete without examination of error norms.

Define the $L_2$-norm and $L_{\infty }$-norm as follows:

$${\parallel u\parallel }_{\infty }={\max }_{\zeta ϵ\Omega }\mid u(\zeta ,t_j)\mid ,$$

$${\parallel u\parallel }_2=\sqrt{h\sum _k^{n+1}u{({\zeta }_k,t_j)}^2},$$

$${\parallel e\parallel }_{\infty }=\max {\parallel u-u_{exact}\parallel }_{\infty },$$

$${\parallel e\parallel }_2={\parallel u-u_{exact}\parallel }_2.$$

Example 1.

Consider the conformable damped Burgers’ equation. The proposed equation is defined below [35,74]:

$$T^{\left(\alpha \right)}u\left(t\right)=\epsilon \left(t\right){\displaystyle \frac{{\partial }^{2}u}{\partial {\zeta }^2}}-\beta \left(t\right)u{\displaystyle \frac{\partial u}{\partial \zeta }}+\lambda \left(t\right)u.$$

(49)

where $\epsilon \left(t\right)$, $\beta \left(t\right)$ and $\lambda \left(t\right)$ are viscosity, advection and damping coefficients, respectively.

Case 1. With $\beta =1$ and $\lambda =0$, boundary conditions are taken to be Dirichlet-type, i.e., $u=0$ at $\zeta =0$ and $\zeta =1$, and the initial condition is taken as $u(\zeta ,0)=sin\left(\pi \zeta \right)$.

In the present case, the exact solution of the problem is unknown for $\alpha <1$. For $\alpha =1$, the exact solution for the considered initial condition is given as [74]:

$$u(\zeta ,t)=2\pi \epsilon {\displaystyle \frac{{\sum }_{n=1}^{\infty }{c}_{n}exp(-n^2{\pi }^2\epsilon t)nsin\left(n\pi \zeta \right)}{c_0+{\sum }_{n=1}^{\infty }{c}_{n}exp(-n^2{\pi }^2\epsilon t)cos\left(n\pi \zeta \right)}}$$

(50)

with the Fourier coefficients defined as

$$c_0={\int }_0^{1}exp\{-{\left(2\pi \epsilon \right)}^{-1}[1-cos\left(\pi \zeta \right)]\}d\zeta ,$$

$$c_n=2{\int }_0^{1}exp\{-{\left(2\pi \epsilon \right)}^{-1}[1-cos\left(\pi \zeta \right)]\}cos\left(n\pi \zeta \right)d\zeta ,n=1,2,3,\dots .$$

The internal dynamics of the given problem is obtained for suitable discretization of the spatial domain with $\Delta t=0.001$ and is presented in Figure 1, Figure 2, Figure 3 and Figure 4. Figure 1 and Figure 2 exhibit 3D visuals of the system for small viscosity coefficients $\epsilon =0.005,0.0005$ with variation in $\alpha =0.25,0.5,0.75,1$. Figure 3 exhibits the slower rate of declining amplitude of the propagation front with the progress of time, as $\alpha$ increases. Figure 4 shows that the propagation front turns steeper as viscosity coefficient decreases from $\epsilon =0.05$ to $\epsilon =0.001$ and leads to shock formation.

Table 2. Comparison of the proposed scheme for Example 1 with HWQA [75] and CHCM [76] for $\alpha =1$ with $\Delta t=0.001$ and $h=0.0125$.

		$\mathit{\epsilon }=0.004$				$\mathit{\epsilon }=0.005$
$\mathit{\zeta }$	$\mathit{t}$	Proposed Scheme	HWQA	CHCM	Exact	Proposed Scheme	HWQA	CHCM	Exact
0.25	1	0.1889046	0.18886	0.18887	0.1889040	0.1887889	0.18874	0.18875	0.1887881
	5	0.0469723	0.04696	0.04695	0.0469723	0.0469635	0.04695	0.04696	0.0469635
	10	0.0242194	0.02421	0.02419	0.0242194	0.0242169	0.02421	0.0242	0.0242168
	15	0.0163154	0.01631	0.01630	0.0163154	0.0163076	0.01631	0.01631	0.0163076
0.5	1	0.3759728	0.37591	0.37594	0.3759762	0.3757202	0.37565	0.37569	0.3757228
	5	0.0939379	0.09393	0.09393	0.0939378	0.0939202	0.09391	0.09393	0.0939201
	10	0.0484372	0.04843	0.04843	0.0484372	0.0484214	0.04842	0.04841	0.0484214
	15	0.0325946	0.03259	0.03259	0.0325946	0.0324388	0.03244	0.03244	0.0324388
0.75	1	0.5588145	0.55875	0.5588	0.5588341	0.5583661	0.55831	0.55836	0.5583839
	5	0.1408869	0.14088	0.14089	0.1408869	0.1408317	0.14083	0.14083	0.1408316
	10	0.0722025	0.07221	0.0721	0.0722025	0.0711341	0.07114	0.07112	0.0711338
	15	0.0467756	0.04679	0.04678	0.0467753	0.0441337	0.04415	0.04414	0.0441329

confirms the competitive nature of the proposed scheme in comparison to the Haar wavelet quasilinearization approach (HWQA) [75] and the cubic Hermite collocation method (CHCM) [76]. The ${\Vert u\Vert }_2$ and ${\Vert u\Vert }_{\infty }$ norms are tabulated in

Table 3. ${\Vert u\Vert }_2$-norm and ${\Vert u\Vert }_{\infty }$-norm of solution of Example 1 with $\Delta t=0.001$.

	$\mathit{\epsilon }=0.0005$ $\mathbf{(}\mathit{h}\mathbf{=}\mathbf{0.002}\mathbf{)}$
	$\mathit{\alpha }\mathbf{=}\mathbf{0.25}$		$\mathit{\alpha }\mathbf{=}\mathbf{0.5}$		$\mathit{\alpha }\mathbf{=}\mathbf{0.75}$
$\mathit{t}$	${\mathbf{\Vert }\mathit{u}\mathbf{\Vert }}_{\mathbf{2}}$	${\mathbf{\Vert }\mathit{u}\mathbf{\Vert }}_{\mathbf{\infty }}$	${\mathbf{\Vert }\mathit{u}\mathbf{\Vert }}_{\mathbf{2}}$	${\mathbf{\Vert }\mathit{u}\mathbf{\Vert }}_{\mathbf{\infty }}$	${\mathbf{\Vert }\mathit{u}\mathbf{\Vert }}_{\mathbf{2}}$	${\mathbf{\Vert }\mathit{u}\mathbf{\Vert }}_{\mathbf{\infty }}$
0.02	4.6075 × ${10}^{-1}$	7.7772 × ${10}^{-1}$	7.0280 × ${10}^{-1}$	9.9316 × ${10}^{-1}$	7.0666 × ${10}^{-1}$	9.9927 × ${10}^{-1}$
0.04	3.7910 × ${10}^{-1}$	6.4728 × ${10}^{-1}$	6.9818 × ${10}^{-1}$	9.9156 × ${10}^{-1}$	7.0637 × ${10}^{-1}$	9.9885 × ${10}^{-1}$
0.06	3.3787 × ${10}^{-1}$	5.7833 × ${10}^{-1}$	6.6847 × ${10}^{-1}$	9.9051 × ${10}^{-1}$	7.0612 × ${10}^{-1}$	9.9849 × ${10}^{-1}$
0.08	3.1140 × ${10}^{-1}$	5.3386 × ${10}^{-1}$	6.3145 × ${10}^{-1}$	9.8806 × ${10}^{-1}$	7.0586 × ${10}^{-1}$	9.9818 × ${10}^{-1}$
0.1	2.9237 × ${10}^{-1}$	5.0164 × ${10}^{-1}$	5.9735 × ${10}^{-1}$	9.6527 × ${10}^{-1}$	7.0558 × ${10}^{-1}$	9.9790 × ${10}^{-1}$
	$\mathit{\epsilon }\mathbf{=}\mathbf{0.005}$ $\mathbf{(}\mathit{h}\mathbf{=}\mathbf{0.005}\mathbf{)}$
0.02	4.5174 × ${10}^{-1}$	7.4829 × ${10}^{-1}$	6.9401 × ${10}^{-1}$	9.8293 × ${10}^{-1}$	7.0451 × ${10}^{-1}$	9.9626 × ${10}^{-1}$
0.04	3.7053 × ${10}^{-1}$	6.1833 × ${10}^{-1}$	6.8152 × ${10}^{-1}$	9.7634 × ${10}^{-1}$	7.0263 × ${10}^{-1}$	9.9374 × ${10}^{-1}$
0.06	3.2946 × ${10}^{-1}$	5.5012 × ${10}^{-1}$	6.5436 × ${10}^{-1}$	9.7127 × ${10}^{-1}$	7.0084 × ${10}^{-1}$	9.9154 × ${10}^{-1}$
0.08	3.0309 × ${10}^{-1}$	5.0598 × ${10}^{-1}$	6.1923 × ${10}^{-1}$	9.6244 × ${10}^{-1}$	6.9902 × ${10}^{-1}$	9.8950 × ${10}^{-1}$
0.1	2.8412 × ${10}^{-1}$	4.7384 × ${10}^{-1}$	5.8606 × ${10}^{-1}$	9.3623 × ${10}^{-1}$	6.9705 × ${10}^{-1}$	9.8761 × ${10}^{-1}$

for $\epsilon =0.005,0.0005$ with $\alpha =0.25,0.5$, and $0.75$, respectively.

Table 4. Comparison of obtained solution for Example 1 with $\alpha =1,\epsilon =0.004$ and $h=0.0125$ at $t=1$.

$\mathit{\zeta }$	$\mathbf{\Delta }\mathit{t}=0.1$	$\mathbf{\Delta }\mathit{t}=0.01$	$\mathbf{\Delta }\mathit{t}=0.001$	Exact
0.2	0.1512765	0.1512136	0.1512080	0.1512074
0.4	0.3014956	0.3014908	0.3014981	0.3014990
0.6	0.4493716	0.4497232	0.4497877	0.4497952
0.8	0.5930046	0.5943589	0.5945706	0.5947340
1	0.0000000	0.0000000	0.0000000	0.0000000

signifies the increasing precision of obtained approximate solution with consideration of small time steps. Further,

Table 5. Numerical simulation of Example 1 for $\Delta t=0.001$ at $t=0.1$.

	$\mathit{\epsilon }=0.0003$ $\mathbf{(}\mathit{h}\mathbf{=}\mathbf{0.002}\mathbf{)}$			$\mathit{\epsilon }=0.0005$ $\mathbf{(}\mathit{h}\mathbf{=}\mathbf{0.002}\mathbf{)}$
$\mathit{\zeta }$	$\mathit{\alpha }\mathbf{=}\mathbf{0.25}$	$\mathit{\alpha }\mathbf{=}\mathbf{0.5}$	$\mathit{\alpha }\mathbf{=}\mathbf{0.75}$	$\mathit{\alpha }\mathbf{=}\mathbf{0.25}$	$\mathit{\alpha }\mathbf{=}\mathbf{0.5}$	$\mathit{\alpha }\mathbf{=}\mathbf{0.75}$
0.1	0.051050	0.110252	0.180629	0.051047	0.110236	0.180601
0.2	0.102076	0.22001	0.357714	0.10207	0.219979	0.357656
0.3	0.153054	0.328748	0.527307	0.153045	0.328700	0.527216
0.4	0.203958	0.435861	0.684527	0.903245	0.435792	0.684394
0.5	0.254759	0.540607	0.822682	0.254792	0.540514	0.822495
0.6	0.305428	0.642014	0.931559	0.305407	0.64189	0.931292
0.7	0.35593	0.738706	0.993458	0.355905	0.738539	0.993063
0.8	0.406226	0.828559	0.971871	0.406195	0.828328	0.971227
0.9	0.456268	0.907874	0.766651	0.456232	0.907535	0.765397
	$\mathit{\epsilon }\mathbf{=}\mathbf{0.003}$ $\mathbf{(}\mathit{h}\mathbf{=}\mathbf{0.005}\mathbf{)}$			$\mathit{\epsilon }\mathbf{=}\mathbf{0.005}$ $\mathbf{(}\mathit{h}\mathbf{=}\mathbf{0.005}\mathbf{)}$
0.1	0.0510075	0.1100448	0.1802532	0.1019247	0.1098897	0.1799736
0.2	0.1019903	0.2195870	0.3569364	0.1019247	0.2192694	0.3563581
0.3	0.1529229	0.3280877	0.5260710	0.1528229	0.3275913	0.5251514
0.4	0.2037786	0.4349270	0.6827295	0.2036421	0.4342244	0.6813914
0.5	0.2545286	0.5393417	0.8201426	0.2543528	0.5383881	0.8182503
0.6	0.3051411	0.6403229	0.9279448	0.3049219	0.6390448	0.9252489
0.7	0.3555800	0.7364290	0.9880943	0.3553120	0.7346995	0.9840868
0.8	0.4058036	0.8253964	0.9631312	0.4054797	0.8229734	0.9565905
0.9	0.4557622	0.9031732	0.7497448	0.4553306	0.8995028	0.7373023

presents numerical results on varying diffusion or viscosity coefficient along with changes in $\alpha$. Increasing $\alpha \to 1$ describes the evolution of the propagation front, and for lower viscosity, it attains its high values or amplitude earlier.

Case 2. When damping effect is introduced in the present shock-exhibiting system, the energy dissipation hinders the formation of the shock and significantly impacts the extremum of the propagation front. Figure 5 depicts the interplay between fractional order, damping and diffusion effects for $\lambda =3$ with variation of $\alpha$ and $\epsilon$.

Example 2.

Consider the conformable damped Burgers’ equation. The proposed equation can be written as follows [35]:

$$T^{\left(\alpha \right)}u\left(t\right)=\epsilon \left(t\right){\displaystyle \frac{{\partial }^{2}u}{\partial {\zeta }^2}}-\beta \left(t\right)u{\displaystyle \frac{\partial u}{\partial \zeta }}+\lambda \left(t\right)u.$$

(51)

where $\epsilon \left(t\right)$, $\beta \left(t\right)$ and $\lambda \left(t\right)\ne 0$ are viscosity, advection and damping coefficients, respectively. Since the exact solution of this equation does not exist as per the knowledge of the authors, the initial and boundary conditions can be extracted from the solution of the particular case of this problem for $\epsilon =1$, $\beta =1$ and $\alpha =1$ given in [27] and is defined by $u(\zeta ,t)={\displaystyle \frac{\lambda \zeta }{2e^{\lambda t}-1}}$.

Numerical results obtained from the cubic Hermite splines collocation method and the exact solution for Equation (51) till time $T=4$ reflect the effect of varying the incorporated coefficients through numerical and graphical simulations. The solution profiles given in Figure 6 and Figure 7 are obtained for $n=200$, $\Delta t=0.001$ for different fractional orders. In the context of the damped Burgers’ equation, the damping coefficient governs the dissipation of energy from the system. Figure 6 shows the increment in the upper bound of the solution with increases in damping coefficient. This suggests a reduction in oscillations and perturbations in the system, which leads the system to achieve a higher upper bound initially before attenuation. In Figure 7, the upper bound of the solution leads to attenuation with the passage of time signifying the pronounced dissipation of energy because of damping, i.e., $\lambda =1$.

The results obtained from the proposed scheme, with $n=200$ and $\Delta t=0.001$, are compared to the results obtained from the homotopy analysis method (HAM) [4], conformable homotopy analysis method (C-HAM) [35] and conformable fractional reduced differential transform method (C-FRDTM) [35] in

Table 6. Comparison of the proposed scheme for Example 2 with HAM [4], C-FRDTM [35] and C-HAM [35] for $\alpha =\epsilon =1$ with $n=200$ and $\Delta t=0.001$.

		$\mathit{\lambda }=1$		$\mathit{\lambda }=1/5$
		Absolute Error [4]	Absolute Error Proposed Scheme	Absolute Error [35]	Absolute Error Proposed Scheme
0.5	0.2	6.2000 × ${10}^{-5}$	5.9074 × ${10}^{-5}$	2.1848 × ${10}^{-5}$	6.3393 × ${10}^{-6}$
	0.4	1.2500 × ${10}^{-4}$	1.1417 × ${10}^{-4}$	4.3696 × ${10}^{-5}$	1.2620 × ${10}^{-5}$
	0.6	1.8700 × ${10}^{-4}$	1.6103 × ${10}^{-4}$	6.5544 × ${10}^{-5}$	1.8769 × ${10}^{-5}$
	0.8	2.4900 × ${10}^{-4}$	1.9490 × ${10}^{-4}$	8.7392 × ${10}^{-5}$	2.4685 × ${10}^{-5}$
	1.0	3.1300 × ${10}^{-4}$	0.000000	1.0924 × ${10}^{-5}$	0.000000
1	0.2	7.9000 × ${10}^{-5}$	1.7931 × ${10}^{-5}$	3.1043 × ${10}^{-4}$	4.5495 × ${10}^{-6}$
	0.4	1.6100 × ${10}^{-4}$	3.4205 × ${10}^{-5}$	6.2086 × ${10}^{-4}$	9.0263 × ${10}^{-6}$
	0.6	2.4200 × ${10}^{-4}$	4.7158 × ${10}^{-5}$	9.3129 × ${10}^{-4}$	1.3357 × ${10}^{-5}$
	0.8	3.2300 × ${10}^{-4}$	5.5109 × ${10}^{-5}$	1.2417 × ${10}^{-3}$	1.7467 × ${10}^{-5}$
	1.0	4.0300 × ${10}^{-4}$	0.000000	1.5522 × ${10}^{-3}$	0.000000

and

Table 7. Comparison of the proposed scheme for Example 2 with C-FRDTM [35] and C-HAM [35] for $\epsilon =1$ and $\lambda =1/5$ with $n=200$ and $\Delta t=0.001$.

		$\mathit{\alpha }=0.8$		$\mathit{\alpha }=0.9$
$\mathit{t}$	$\mathit{\zeta }$	[35]	Proposed Scheme	[35]	Proposed Scheme
0.5	0.2	0.030474	0.032724	0.031882	0.032875
	0.4	0.060948	0.065537	0.063764	0.065800
	0.6	0.091422	0.098521	0.095647	0.098818
	0.8	0.121896	0.131748	0.127529	0.131966
	1.0	0.152370	0.165243	0.159411	0.165242
1	0.2	0.024792	0.027690	0.026246	0.027709
	0.4	0.049583	0.055393	0.052492	0.055424
	0.6	0.074375	0.083118	0.078738	0.083151
	0.8	0.099167	0.110870	0.104984	0.110890
	1.0	0.123958	0.138619	0.131230	0.138619

. The results obtained from the present scheme are found to be better than the results obtained from the other schemes given in the literature. Table 7 shows the comparison of numerical and semi-analytic results for different fractional orders. In

Table 8. ${\Vert e\Vert }_2$ and ${\Vert e\Vert }_{\infty }$ error norms for Example 2 when $\alpha =\epsilon =1$ and $\lambda =1$ with $\Delta t=0.001$.

	$\mathit{n}=100$		$\mathit{n}=200$		$\mathit{n}=400$
$\mathit{t}$	${\mathbf{\Vert }\mathit{e}\mathbf{\Vert }}_{\mathbf{2}}$	${\mathbf{\Vert }\mathit{e}\mathbf{\Vert }}_{\mathbf{\infty }}$	${\mathbf{\Vert }\mathit{e}\mathbf{\Vert }}_{\mathbf{2}}$	${\mathbf{\Vert }\mathit{e}\mathbf{\Vert }}_{\mathbf{\infty }}$	${\mathbf{\Vert }\mathit{e}\mathbf{\Vert }}_{\mathbf{2}}$	${\mathbf{\Vert }\mathit{e}\mathbf{\Vert }}_{\mathbf{\infty }}$
0.05	6.8261 × ${10}^{-4}$	1.7816 × ${10}^{-3}$	3.4998 × ${10}^{-4}$	9.0431 × ${10}^{-4}$	1.8126 × ${10}^{-4}$	4.5554 × ${10}^{-4}$
0.25	5.2863 × ${10}^{-4}$	8.9567 × ${10}^{-4}$	2.7589 × ${10}^{-4}$	4.5094 × ${10}^{-4}$	1.4859 × ${10}^{-4}$	2.2624 × ${10}^{-4}$
0.5	2.7049 × ${10}^{-4}$	4.1801 × ${10}^{-4}$	1.4273 × ${10}^{-4}$	2.1029 × ${10}^{-4}$	7.8402 × ${10}^{-5}$	1.0762 × ${10}^{-4}$
1	7.5265 × ${10}^{-5}$	1.1209 × ${10}^{-4}$	4.0807 × ${10}^{-5}$	5.6794 × ${10}^{-5}$	2.3471 × ${10}^{-5}$	3.1060 × ${10}^{-5}$
4	3.2095 × ${10}^{-7}$	4.2808 × ${10}^{-7}$	2.6716 × ${10}^{-7}$	3.6273 × ${10}^{-7}$	2.4075 × ${10}^{-7}$	3.3091 × ${10}^{-7}$

, ${\Vert e\Vert }_2$-norm and ${\Vert e\Vert }_{\infty }$-norm are presented for $\epsilon =1,\lambda =1$ with successive refinements of the space domain. It signifies the importance of precise refinement.

Table 9. ${\Vert u\Vert }_2$-norm and ${\Vert u\Vert }_{\infty }$-norm of solution of Example 2 for $\epsilon =1$ and $\lambda =1$ with $\Delta t=0.001$ and $n=200$.

	$\mathit{\alpha }=0.25$		$\mathit{\alpha }=0.5$		$\mathit{\alpha }=0.75$
$\mathit{t}$	${\mathbf{\Vert }\mathit{u}\mathbf{\Vert }}_{\mathbf{2}}$	${\mathbf{\Vert }\mathit{u}\mathbf{\Vert }}_{\mathbf{\infty }}$	${\mathbf{\Vert }\mathit{u}\mathbf{\Vert }}_{\mathbf{2}}$	${\mathbf{\Vert }\mathit{u}\mathbf{\Vert }}_{\mathbf{\infty }}$	${\mathbf{\Vert }\mathit{u}\mathbf{\Vert }}_{\mathbf{2}}$	${\mathbf{\Vert }\mathit{u}\mathbf{\Vert }}_{\mathbf{\infty }}$
0.05	4.7338 × ${10}^{-1}$	9.0699 × ${10}^{-1}$	4.7820 × ${10}^{-1}$	9.0699 × ${10}^{-1}$	4.9665 × ${10}^{-1}$	9.0699 × ${10}^{-1}$
0.25	3.4455 × ${10}^{-1}$	6.3773 × ${10}^{-1}$	3.4885 × ${10}^{-1}$	6.3773 × ${10}^{-1}$	3.5592 × ${10}^{-1}$	6.3773 × ${10}^{-1}$
0.5	2.4184 × ${10}^{-1}$	4.3527 × ${10}^{-1}$	2.4425 × ${10}^{-1}$	4.3527 × ${10}^{-1}$	2.4751 × ${10}^{-1}$	4.3527 × ${10}^{-1}$
1	1.2991 × ${10}^{-1}$	2.2540 × ${10}^{-1}$	1.3013 × ${10}^{-1}$	2.2540 × ${10}^{-1}$	1.3038 × ${10}^{-1}$	2.2540 × ${10}^{-1}$
4	6.0887 × ${10}^{-3}$	9.2425 × ${10}^{-3}$	5.7361 × ${10}^{-3}$	9.2425 × ${10}^{-3}$	5.5071 × ${10}^{-3}$	9.2425 × ${10}^{-3}$

shows computed ${\Vert u\Vert }_2$ and ${\Vert u\Vert }_{\infty }$ norms of numerical solution for $\epsilon =1,\lambda =1$ at different values of $\alpha$. In

Table 10. Numerical simulation of Example 2 for $\epsilon =1$ and $\alpha =1$ with $\Delta t=0.001$ and $n=200$.

	$\mathit{\lambda }=0.001$			$\mathit{\lambda }=0.1$			$\mathit{\lambda }=1$
$\mathit{\zeta }$	$\mathit{t}\mathbf{=}\mathbf{0.25}$	$\mathit{t}\mathbf{=}\mathbf{0.5}$	$\mathit{t}\mathbf{=}\mathbf{1}$	$\mathit{t}\mathbf{=}\mathbf{0.25}$	$\mathit{t}\mathbf{=}\mathbf{0.5}$	$\mathit{t}\mathbf{=}\mathbf{1}$	$\mathit{t}\mathbf{=}\mathbf{0.25}$	$\mathit{t}\mathbf{=}\mathbf{0.5}$	$\mathit{t}\mathbf{=}\mathbf{1}$
0.1	9.9950 × ${10}^{-5}$	9.9900 × ${10}^{-5}$	9.9800 × ${10}^{-5}$	9.5189 × ${10}^{-3}$	9.0709 × ${10}^{-3}$	8.2629 × ${10}^{-3}$	6.3822 × ${10}^{-2}$	4.3556 × ${10}^{-2}$	2.2549 × ${10}^{-2}$
0.2	1.9990 × ${10}^{-4}$	1.9980 × ${10}^{-4}$	1.9960 × ${10}^{-4}$	1.9038 × ${10}^{-2}$	1.8142 × ${10}^{-2}$	1.6526 × ${10}^{-2}$	1.2764 × ${10}^{-1}$	8.7112 × ${10}^{-2}$	4.5098 × ${10}^{-2}$
0.3	2.9985 × ${10}^{-4}$	2.9970 × ${10}^{-4}$	2.9940 × ${10}^{-4}$	2.8557 × ${10}^{-2}$	2.7213 × ${10}^{-2}$	2.4789 × ${10}^{-2}$	1.9147 × ${10}^{-1}$	1.3067 × ${10}^{-1}$	6.7646 × ${10}^{-2}$
0.4	3.9980 × ${10}^{-4}$	3.9960 × ${10}^{-4}$	3.9920 × ${10}^{-4}$	3.8076 × ${10}^{-2}$	3.6283 × ${10}^{-2}$	3.3052 × ${10}^{-2}$	2.5529 × ${10}^{-1}$	1.7422 × ${10}^{-1}$	9.0194 × ${10}^{-2}$
0.5	4.9975 × ${10}^{-4}$	4.9950 × ${10}^{-4}$	4.9900 × ${10}^{-4}$	4.7595 × ${10}^{-2}$	4.5354 × ${10}^{-2}$	4.1315 × ${10}^{-2}$	3.1911 × ${10}^{-1}$	2.1777 × ${10}^{-1}$	1.1274 × ${10}^{-1}$
0.6	5.9970 × ${10}^{-4}$	5.9940 × ${10}^{-4}$	5.9880 × ${10}^{-4}$	5.7114 × ${10}^{-2}$	5.4425 × ${10}^{-2}$	4.9577 × ${10}^{-2}$	3.8294 × ${10}^{-1}$	2.6132 × ${10}^{-1}$	1.3529 × ${10}^{-1}$
0.7	6.9965 × ${10}^{-4}$	6.9930 × ${10}^{-4}$	6.9860 × ${10}^{-4}$	6.6633 × ${10}^{-2}$	6.3496 × ${10}^{-2}$	5.7840 × ${10}^{-2}$	4.4676 × ${10}^{-1}$	3.0487 × ${10}^{-1}$	1.5783 × ${10}^{-1}$
0.8	7.9960 × ${10}^{-4}$	7.9920 × ${10}^{-4}$	7.9840 × ${10}^{-4}$	7.6153 × ${10}^{-2}$	7.2567 × ${10}^{-2}$	6.6103 × ${10}^{-2}$	5.1057 × ${10}^{-1}$	3.4841 × ${10}^{-1}$	1.8037 × ${10}^{-1}$
0.9	8.9955 × ${10}^{-4}$	8.9910 × ${10}^{-4}$	8.9820 × ${10}^{-4}$	8.5672 × ${10}^{-2}$	8.1638 × ${10}^{-2}$	7.4366 × ${10}^{-2}$	5.7438 × ${10}^{-1}$	3.9195 × ${10}^{-1}$	2.0292 × ${10}^{-1}$

, the numerical simulation of effect of varying $\lambda$ again demonstrates the attentuation of the solution with the passage of time.

9. Conclusions

In this article, a cubic Hermite splines collocation method in space with an implicit finite difference scheme in time is proposed for solving a modified form of time-fractional Burgers’ equation with linear damping known as conformable damped Burgers’ equation with time-dependent variable coefficients. The use of the conformable time-fractional operator provided a convenient way to study the fractional dynamics of the system. For non-linear advection term, a linear approximation is estimated by quasi-linearization. The behavioral analysis of viscosity or diffusion coefficient and damping coefficient highlights the importance of understanding the interplay between damping, diffusion and fractional order in non-linear partial differential equations and their implications for physical systems. In the present case, damping and viscosity play crucial roles in controlling the amplitude and attenuation of the solution profile. A small viscosity can potentially lead to the shock formation. The error analysis shows that one needs a sufficiently small partition of the respective domains to efficiently capture the dynamics of the solution profile. Moreover, graphical analysis suggests the significance of fine partitions to draw abrupt changes in the propagation front that may be caused as $\epsilon \to 0$.