Application of the B-spline Galerkin approach for approximating the time-fractional Burger's equation

1.   Introduction

Fractional calculus (FC) theory was proposed by N. H. Abel and J. Liouville, and a description of their work is presented in ^[1]. By using FC, integer derivatives, and integrals can be generalized to real or variable derivatives and integrals. FC is studied since fractional differential equations (FDEs) are better suited to modeling natural physics processes and dynamic systems than integer differential equations. Furthermore, FDEs that incorporate memory effects are better suited to describing natural processes that have memory and hereditary properties. In other words, because fractional derivatives have memory effects, FDEs are more accurate in describing physical phenomena with memory or hereditary characteristics. There was a trend to consider FC to be an esoteric theory with no application until the last few years. Now, more and more researchers are investigating how it can be applied to economics, control system and finance. As a result, many fractional order differential operators were developed, such as Hadamard, Riemann-Liouville, Caputo, Riesz, Grünwald-Letnikov, and variable order differential operators. The researchers have devoted considerable effort to solving FDEs numerically so that they can be applied to a variety of problems ^{[2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]}. Several numerical approaches have been proposed in the literature, including eigenvector expansion, the fractional differential transform technique ^[21], the homotopy analysis technique ^[22], the homotopy perturbation transform technique ^[23], the generalized block pulse operational matrix technique ^[24] and the predictor-corrector technique ^[25]. In addition, the use of Legendre wavelets to integrate and differentiate fractional order matrices has been suggested as a numerical method ^[26,27].

In this paper, we study the numerical solution of the time-fractional Burger's equation (TFBE) ^[28] as follows:

which is subject to the following boundary conditions (BCs):

and the following initial condition (IC):

in which $0 < \gamma \le 1$ is a parameter representing the order of the fractional time, $v$ denotes a viscosity parameter and $g\left(x\right), {l}_{1}\left(t\right)\;\mathrm{a}\mathrm{n}\mathrm{d}\;{l}_{2}\left(t\right)$ are given functions of their argument. The TFBE is a kind of sub-diffusion convection, which is widely adopted to describe many physical problems such as unidirectional propagation of weakly nonlinear acoustic waves, shock waves in flow systems, viscous media, compressible turbulence, electromagnetic waves and weak shock propagation ^[29,30,31]. In recent years, there has been some technique development in the study of Burger's equation: an implicit difference scheme and algorithm implementation ^[32], pointwise error analysis of the third-order backward differentiation formula (BDF3) ^[33], pointwise error estimates of a compact difference scheme ^[34], efficient (BDF3) finite-difference scheme ^[35], semi-analytical methods ^[36], composite spectral methods ^[37], least-squares methods ^[38], geometric analysis methods ^[39], error and stability estimate techniques ^[40].

Definition 1. Suppose that $m$ is the smallest integer exceeding $\gamma$; the Caputo time fractional derivative operator of order $\gamma > 0$ can be defined as follows ^[41]:

where $u\left(x, t\right)$ is the unknown function that is $(m-1)$ times continuously differentiable and $\mathit{\Gamma} \left(.\right)$ denotes the usual gamma function. The finite-element method has been an important method for solving both ordinary and partial differential equation, therefore, in recent research, it has been applied to solve the TFBE. In what follows, we describe the solution process by using the finite-element scheme for solving the TFBE.

2.   Cubic B-Spline Galerkin method (CBSGM) with a quadratic weight function:

To discretize the TFBE (1), first let us define the cubic B-spline base function. We partition the interval $[a, b]$, which represents the solution domain of (1) into $M$ uniformly spaced points ${x}_{m}$ such that $a = {x}_{0} < {x}_{1} < \dots < {x}_{M-1} < {x}_{M} = b$ and $h = ({x}_{m+1}-{x}_{m})$. Then, the cubic B-spline ${C}_{m}\left(x\right), (m = -1\left(1\right)(M+1)$, at the knots ${x}_{m}$ which form basis on the solution interval $[a, b]$, is defined as follows ^[42]:

where the set of cubic B-splines ${({C}_{-1}\left(x\right), C}_{0}\left(x\right), \dots , {C}_{M}\left(x\right), {C}_{M+1}\left(x\right))$ is a basis for the functions defined over interval $[a, b]$. Thus, the numerical solution ${U}_{M}(x, t)$ to the analytic solution $U(x, t)$ can be illustrated as

where ${\sigma }_{m}\left(t\right)$ are unknown time-dependent parameters to be determined from the initial, boundary and weighted residual conditions. Since each cubic B-spline covers four consecutive elements, each element ${[x}_{m}, {x}_{m+1}]$ is also covered by four cubic B-splines. So, the nodal values ${U}_{m}$ and its first and second derivatives ${U\mathrm{\text{'}}}_{m}$, ${U"}_{m}$ can be respectively computed in terms of the element parameter ${\sigma }_{m}\left(t\right)$, at the knot ${x}_{m}$ as follows:

and by means of the local coordinate transformation ^[43] as follows:

A cubic B-spline shape function in terms of $\eta$ over the element ${[x}_{m}, {x}_{m+1}]$ is formulated as:

and the variation of ${U}_{M}(\eta , t)$ over the typical element ${[x}_{m}, {x}_{m+1}]$ is represented as

in which B-splines ${C}_{m-1}\left(\eta \right), {C}_{m}\left(\eta \right), {C}_{m+1}\left(\eta \right)$, ${C}_{m+2}\left(\eta \right)$ and ${\sigma }_{m-1}\left(t\right)$, ${\sigma }_{m}\left(t\right), {\sigma }_{m+1}$ and ${\sigma }_{m+2}\left(t\right)$ are element shape functions and element parameters, respectively.

Based on the Galerkin's method with weight function $W\left(x\right) > 0,$ we get the following weak formula of (1):

using transformation (8) and by apply partial integration we obtain:

where $\mathrm{\lambda } = \frac{1}{h}{\rm{Û}}$, $\mathit{\Phi} = \frac{v}{{h}^{2}}$ and ${\rm{Û}} = U(\eta , t)$ which is considered to be a constant on an element to simplify the integral ^[43]; replace the weight function $W$ by quadratic B-spline ${B}_{\boldsymbol{m}}\left(\boldsymbol{x}\right), \mathrm{m} = -1\left(1\right)M$, at the knots ${x}_{m}$, which forms a basis on the solution interval $[a, b]$, introduced as follows ^[44]:

where ${({B}_{-1}\left(x\right), B}_{0}\left(x\right), \dots , {B}_{M}\left(x\right))$ is the set of splines for the basis of functions introduced on $[a, b]$. The numerical solution ${U}_{M}(x, t)$ to the analytic solution $U(x, t)$ is expanded by

where ${\vartheta }_{m}$ are unknown time-dependent parameters, and by using local coordinate transformation (8), the quadratic B-spline shape functions for the typical element [${x}_{m}, {x}_{m+1}$] are given as

The variation of the function $U(\eta , t)$ is approximated by

where ${\vartheta }_{m-1}\left(t\right)$, ${\vartheta }_{m}\left(t\right)$ and ${\vartheta }_{m+1}\left(t\right)$ act as element parameters and B-splines ${B}_{m-1}\left(\eta \right), {B}_{m}\left(\eta \right)$ and ${B}_{m+1}\left(\eta \right)$ as element shape functions based on the above; (12) will be in the following form:

in which "Dot" represents the $\sigma \mathrm{t}\mathrm{h}$ fractional derivative with respect to time. We can write (17) in matrix notation as follows:

in which ${\sigma }^{e} = {({\sigma }_{m-1}, {\sigma }_{m}, {\sigma }_{m+1}, {\sigma }_{m+2})}^{T}$ are the element parameters. The element matrices ${X}_{ij}^{e}, {Y}_{ij}^{e}, {Z}_{ij}^{e}, {Q}_{ij}^{e}$ and ${E}_{i}^{e}$ are rectangular $3\times 4$ matrices introduced through the following integrals:

where $i$ and $j$ take only the values $(m-1, m, m+1)$ and $(m-1, m, m+1, m+2)$ respectively, and a lumped value for λ is defined by $\mathrm{\lambda } = \frac{1}{2h}({\sigma }_{m-1}+{5\sigma }_{m}+5{\sigma }_{m+1}+{\sigma }_{m+2})$.

By assembling all contributions from all elements, we get the following matrix equation:

where $\sigma = {({\sigma }_{-1}, {\sigma }_{0}, {\sigma }_{1}, \dots , {\sigma }_{M}, {\sigma }_{\mathrm{M}+1})}^{T}$ denotes a global element parameter. The matrices $\mathrm{X}, \;\mathrm{Z}$ and $Y$ represent rectangular, septa-diagonal and every sub-diagonal matrices, which include the following forms:

in which,

Following ^[45], we can approximate the temporal Caputo derivative with the help of the $L1$ formula:

where ${b}_{k}^{\gamma } = ({k+1)}^{1-\gamma }-{k}^{1-\gamma }\;\mathrm{a}\mathrm{n}\mathrm{d}\;\Delta t = \frac{{t}_{f}-0}{N}$, and ${t}_{f} = n\left(\Delta t\right), n = \mathrm{0, 1}, \dots N, \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}N$ represents a positive integer. Now, we recall the following lemma.

Lemma 1: Suppose that $0 < \gamma < 1and{b}_{k}^{\gamma } = ({k+1)}^{1-\gamma }-{k}^{1-\gamma }, k = \mathrm{0, 1}, \dots ;then, 1 = {b}_{0}^{\gamma } > {b}_{1}^{\gamma } > \dots > {b}_{k}^{\gamma }\to 0, ask\to \mathrm{\infty }$ ^[46].

Then, we can we write the parameter ${\sigma }_{m}^{ \cdot }$ as follows:

while the parameter $\sigma$ by the Crank-Nicolson scheme, is as follows:

Substitution both parameters above into $\left(18\right)$, we obtain the $(M+2)\times (M+3)$ matrix system:

where $\sigma = {({\sigma }_{m-2}+{\sigma }_{m-1}+{\sigma }_{m}+{\sigma }_{m+1}+{\sigma }_{m+1}+{\sigma }_{m+2}+{\sigma }_{m+3})}^{T}$; to make the matrix equation be square, we need to find an additional constraint of BC (2) and their second derivatives and we obtain discard ${\sigma }_{-1}$ from system (20) as follows:

the variables ${\sigma }_{-1}^{n}$ and ${\sigma }_{M+1}^{n}$ can be ignored from system (20) and then the system can be converted to an $(M+1)\times (M+1)$ matrix system. The initial vector of parameter ${\sigma }^{0} = ({\sigma }_{0}^{0}, {\sigma }_{1}^{0}, \dots , {\sigma }_{M}^{0})$ should be obtained to iterate system (20); the approximation of (6) has been reformulated on the interval $[a, b]$ when time $t = 0$ as follows:

where $U\left(x, 0\right)$ fulfills the following equation at node ${x}_{m}:$

Therefore, we can obtain the following system:

and we solve this identity matrix by applying the Jain algorithm ^[47].

3.   Stability analysis

This section adopts the von Neumann stability analysis to investigate the stability of approximation obtained by scheme (20). First, we introduce the recurrence relationship between successive time levels relating unknown element parameters ${\sigma }_{m}^{n+1}\left(t\right)$, as follows:

where

and $\alpha = {\left(\Delta t\right)}^{-\gamma }\mathit{\Gamma} \left(2-\gamma \right)$.

The growth factor of the typical Fourier mode is defined as

where, $i = \sqrt{-1, }\beta$ is a mode number and $h$ is the element size. Substitution of (22) into (21) yields

let ${\xi }^{n+1} = Ϋ{\xi }^{n}$ and assume that $Ϋ\equiv Ϋ\left(\theta \right)$ is independent of time, therefore, we can write $Ϋ$ as follows:

where

Obviously note that $\left|\mathrm{Ϋ}\right|\le 1$. Therefore, according to the Fourier condition, the scheme (20) is unconditionally stable.

4.   Numerical results

This section introduces two numerical examples, which highlight numerical results for the TFBE with different IC and BCs given by the CBSGM with quadratic weight function. In this section, we use the ${L}_{2}$ and ${L}_{\infty }$ to calculate the accuracy of the CBSGM with a quadratic weight function, which has been employed in this study; we will also show how the analytical results and the numerical results are close to each other. To do this, first we will find the exact solutions to the problem (1) by applying the following problems; then, we compare the results with the numerical solution obtained from the given method. To this aim, the ${\boldsymbol{L}}_{{\infty }}$ and ${\boldsymbol{L}}_{2}$ error norms are respectively defined as ^[48]

where $U$ and ${U}_{M}$ represent the exact solution and numerical solution, respectively.

Example 1: Let us consider the TFBE (1) with the BCs

and IC

such that the forcing term $f(x, t)$ is achieved as ^[45]

where the analytic solution is obtained as

Numerical results are reported in Tables 1–3 and Figure 1. Table 1 lists the numerical solutions and the ${L}_{2}$ and ${L}_{\infty }$ error norms with $\gamma = 0.5, \Delta t = 0.0025, {t}_{f} = 0.05\;\mathrm{a}\mathrm{n}\mathrm{d}\;v = 1$ for various numbers of partitions $M.$ As seen in Table 1, we notice that when the number of partitions M are increased, the ${L}_{\infty }$ and ${L}_{2}$ error norms will decrease considerably. Table 2 displays the numerical solutions with $\gamma = 0.5, M = 40, t = 1, {t}_{f} = 0.05\;\mathrm{a}\mathrm{n}\mathrm{d}\;v = 1$ for various values of $\Delta t.$ In view of Table 2, we can see that when $\Delta t$ decreases, the ${L}_{\infty }$ and ${L}_{2}$ error norms decrease, as was expected. Table 3 shows the numerical solutions with $\Delta t = 0.00025, M = 40, t = 1, {t}_{f} = 0.05, v = 1$ for various values of $\gamma$. As observed in Table 3, the ${L}_{\infty }$ and ${L}_{2}$ error norms decrease when γ increases. A comparison between the results of our proposed strategy and two other methods is demonstrated in detail, the researchers of which relied on their work on a weight function corresponding to the spline function in terms of degree; see ^[44,45]. Figure 1 represents the surfaces of the exact and numerical solutions of the TFBE in Example (1).

Example 2: Finally, we consider the TFBE (1) with the BCs

and IC

where the source term $f\left(x, t\right)$ can be obtained as ^[44]

The exact solution is

Numerical results are represented in Tables 4 and 5 and Figure 2. Tables 4 and 5 report the numerical solutions for various numbers of partitions $M$ and values of $\Delta t$. As seen in Tables 4 and 5, when the number of partitions M increased, the error norms ${L}_{\infty }$ and ${L}_{\bf{2}}$ will decrease considerably, while, in Table 5, we can see that when ∆t decrease, the error norms ${L}_{\infty }$ and ${L}_{\bf{2}}$ decrease. Figure 2 demonstrates the surfaces of the exact and numerical solutions of the TFBE in Example (2).

5.   Conclusions

This paper presented a numerical approach based on the CBSGM with a quadratic weight function for the TFBE including the time Caputo derivative. Numerical results have shown that the proposed method is an appropriate and efficient scheme for solving such problems.

Conflict of interest

The authors declare no conflict of interest.