An efficient higher-order trigonometric cubic B-spline collocation method for timefractional Burgers equations

Murat Önal; Berat Karaagac; Alaattin Esen

doi:10.2478/ijmce-2026-0013

Introduction

The history of fractional analysis, which is a generalized form of the derivative and integral topics discussed in classical analysis to non-fractional order derivatives and integrals, dates back to the 1695 notes between Leibniz and L’hospital discussing the (1/2)^th order derivative of a function. Despite the early interest of mathematicians such as Euler, Laplace and Fourier, fractional analysis remained a purely theoretical field for more than two centuries. However, during the last few centuries many researchers have shown that fractional order analysis is more adequate than classical analysis in describing many real materials. The most fundamental and important property of fractional order operators is their ability to describe the memory hereditary effects of many materials and physical processes [1]. For example, in health science the resistance of the brain to drug on tumor cells and stages of a tumor can be described by fractional calculus [2]. Similarly, in the study of deforestation and its effects on the biodiversity of wildlife species, fractional order differential equations are used due to their ability to deal with complex systems [3]. Over the years, a large number of mathematicians have presented various definitions that fit the idea of a non-integer integral or derivative, using their own notations and approaches. Grünwald and Letnikov introduced a limit-based formulation in the 1860s, now known as the Grünwald–Letnikov fractional derivative [4]. Shortly thereafter, Riemann formalized what became known as the Riemann-Liouville derivative [1]. Since then, many applications of the fractional derivatives and integrals of the Riemann–Liouville type have been demonstrated in various domains of science and technology. In the early twentieth century, Caputo provided an alternative definition of a fractional derivative related to the Riemann–Liouville integral, which gained popularity due to its compatibility with classical initial conditions [5]. The Riemann–Liouville derivative has its foundation in theory, though the Caputo derivative is frequently used in physical modeling due to its suitability with conventional initial and boundary conditions. In addition to the definitions presented, the concept of fractional analysis has evolved over the centuries due to the industrialization of the complexities associated with a heterogeneous formation, the ability of fractional operators to capture the behavior of multi-dimensional environments due to their diffusion operations, and the rapid development of mathematical techniques using computer software, allowing researchers to present their perspectives while analyzing complex formations, providing new definitions and theorems, hundreds of books, and scientific publications. Also, significant progress has been made in the numerical analysis of nonlinear fractional differential equations, including studies on fractional Riccati equations, complex spatiotemporal dynamics in fractional systems, fractional ecological models, and fractional Burger-Huxley equations, highlighting the growing interest in accurate and efficient numerical schemes for nonlinear fractional problems [6–9]. In this study, we are going to consider the fractional Burgers equation with the boundary and initial conditions given as follows, respectively; 1 $D_{t}^{γ} u + u u_{x} - v u_{x x} = f (x, t),$ D_t^\gamma u + u{u_x} - v{u_{xx}} = f(x,t), and, 2 $\begin{array}{l} u (a, t) = h_{1} (t), u (b, t) = h_{2} (t), \\ u_{x} (a, t) = h_{3} (t), u_{x} (b, t) = h_{4} (t), \\ u (x, 0) = g (x), \end{array}$ \matrix{ {u(a,t) = {h_1}(t),\,\,u(b,t) = {h_2}(t),} \hfill \cr {{u_x}(a,t) = {h_3}(t),\,\,{u_x}(b,t) = {h_4}(t),} \hfill \cr {u(x,0) = g(x),} \hfill \cr } where x ∈ [a,b], t ∈ [0,T] and v is the viscosity parameter, D^γ symbolize the fractional order derivative. The fractional Burgers equation, which belongs to the class of fractional order partial differential equations, has a wide range of applications in various fields such as heat transfer processes, dispersed water phenomena, sound and shallow water wave propagation analysis. Due to its wide range of applications in modeling nonlinear and anomalous diffusion phenomena, the fractional Burgers equation has garnered significant attention from researchers. Various numerical methods have been developed to solve the equation, as evidenced by numerous studies in literature. For instance, Esen and Tasbozan [10] have obtained the numerical solutions of the fractional order Burgers equation using the cubic B-spline Collocation method. Singh and Gupta [11], in their study which is based on the trigonometric tension B-spline collocation method, have addressed the fractional-order Burgers equation. Duangpan et al. [12] have obtained results for the fractional Burgers equation by a finite integration approach based on Chebyshev polynomials. In [13], Tamsir et al. have applied Hybrid B-spline collocation approach to obtain numerical results of the equation. And, Jeon and Bu [14] have obtained solutions of the Burgers equation using Adams–Moulton and linearized technique. A direct numerical method for solving time-fractional Burger’s equation using a modified hybrid B-spline basis function is presented by Hashmi et al. in [15]. In [16], Sadri et al. developed a novel pseudo-operational collocation method using bivariate Jacobi polynomials and investigated parameter effects, error bounds, and existence-uniqueness results for the time-fractional Burgers equation. In [17], Yaseen and Abbas introduced a cubic trigonometric B-spline method combined with a linearization technique to efficiently solve the time-fractional Burgers equation and established unconditional stability. In [18], Al-saedi and Rashindinia have applied the Galerkin Finite element method based on cubic B-spline basis function to the fractional-order Burgers equation. Maji and Natesan [19] have investigated the numerical results of the equation using the discontinuous Galerkin method. Raul and Rohil [20] have obtained numerical solutions of the problem with the optimal quantic B-spline method. Additionally, Taghipour and Aminikhah have used Pell collocation method for numerical results in the paper referred as [21]. Danaf and Hadhoud [22] have used parametric spline functions and obtained numerical solutions.

In this paper, we are going to consider the fractional Burgers equation using efficient higher-order collocation approach based on trigonometric cubic B-splines. Although various spline-based collocation methods have been proposed for the time-fractional Burgers equation, most existing approaches are limited to the secondor third-order spatial accuracy, or rely on polynomial B-spline bases. In contrast, the present study employs a fourth-order accurate collocation framework constructed with cubic trigonometric B-splines, which provides improved resolution of nonlinear convective effects. The proposed formulation preserves the local support and computational efficiency of spline methods while achieving enhanced accuracy without increasing the stencil width.

First of all, some definitions are going to present in Section 2. Then, the discretization of the time-fractional Burgers equation using L1 algorithm and Crank-Nicolson approach with the processes obtaining numerical algorithm for the equation is presented in Section 3. In Section 4, the stability of the numerical algorithm which is obtained in Section 3 will be discussed. Lastly, the numerical solutions of the three examples related to timefractional Burgers equation are investigated and numerical solutions of the problems are presented in Section 5 with tables, error norms and graphs. In Section 6, the conclusion of the study is given.

Definitions and preliminaries

In this section, a brief overview of essential definitions and results from fractional calculus is presented.

Definition 1. Riemann-Liouville fractional derivative

Let f ∈ L₁(a,t) and γ > 0. The left-sided RiemannLiouville fractional derivative of order γ is defined by $^{R L} D_{a}^{γ} f (t) = \frac{1}{Γ (n - γ)} \frac{d^{n}}{d t^{n}} \int_{a}^{t} \frac{f (τ)}{{(t - τ)}^{γ - n + 1}} d τ,$ ^{RL}D_a^\gamma f(t) = {1 \over {\Gamma (n - \gamma )}}{{{d^n}} \over {d{t^n}}}\int_a^t {{{f(\tau )} \over {{{(t - \tau )}^{\gamma - n + 1}}}}} d\tau , where $n = ⌈ γ ⌉$ n = \left\lceil \gamma \right\rceil and Γ(·) denotes the Gamma function [1].

Definition 2. Caputo fractional derivative

Let f ∈ Cⁿ(a,t) and γ > 0. The left-sided Caputo fractional derivative of order γ is defined by $^{C} D_{a}^{γ} f (t) = \frac{1}{Γ (n - γ)} \int_{a}^{t} \frac{f^{(n)} (τ)}{{(t - τ)}^{γ - n + 1}} d τ,$ ^CD_a^\gamma f(t) = {1 \over {\Gamma (n - \gamma )}}\int_a^t {{{{f^{(n)}}(\tau )} \over {{{(t - \tau )}^{\gamma - n + 1}}}}} d\tau , where $n = ⌈ γ ⌉, f^{(n)}$ n = \left\lceil \gamma \right\rceil ,{f^{(n)}} denotes the n-th derivative of f, and Γ(·) denotes the Gamma function [1].

2.1

Trigonometric B-spline functions

A fundamental part of the applied mathematics, approximation problems are the representation of a known or unknown function by means of a well-defined set of special functions that are easy to compute and have certain analytic properties. Algebraic and trigonometric polynomials, exponential functions, multinomial and trigonometric splines are examples of special functions [23]. Trigonometric spline functions were first introduced in 1964 by Schoenberg [24]. Due to their geometric characteristics and shape-preserving properties, their continuous nature, positive definiteness, and partition of unity properties, trigonometric spline functions have attracted considerable attention in image processing, computer graphics, CAM design, trajectory generation, and many other mathematical and physical fields [25]. The following recursive formula defines trigonometric B-spline functions of degree k. Let {x_i} denote a uniform partition of the spatial domain with step size h = x_(i+1)–x_i.

T_{i}^{k} (x) = \frac{\sin (\frac{x - x_{i}}{2})}{\sin (\frac{x_{i + k} - x_{i}}{2})} T_{i}^{k - 1} (x) + \frac{\sin (\frac{x_{i + k + 1} - x}{2})}{\sin (\frac{x_{i + k + 1} - x_{i + 1}}{2})} T_{i + 1}^{k - 1} (x) .

T_i^k(x) = {{\sin \left( {{{x - {x_i}} \over 2}} \right)} \over {\sin \left( {{{{x_{i + k}} - {x_i}} \over 2}} \right)}}T_i^{k - 1}(x) + {{\sin \left( {{{{x_{i + k + 1}} - x} \over 2}} \right)} \over {\sin \left( {{{{x_{i + k + 1}} - {x_{i + 1}}} \over 2}} \right)}}T_{i + 1}^{k - 1}(x).

Starting with k = 3 in the Eq.(3) trigonometric cubic B-spline functions can be obtained as follows 4 $T_{i}^{3} (x) = \frac{p (x_{i})}{\sin (\frac{3 h}{2})} T_{i}^{2} (x) + \frac{p (x_{i + 4})}{\sin (\frac{3 h}{2})} T_{i + 1}^{2} (x),$ T_i^3(x) = {{p\left( {{x_i}} \right)} \over {\sin \left( {{{3h} \over 2}} \right)}}T_i^2(x) + {{p\left( {{x_{i + 4}}} \right)} \over {\sin \left( {{{3h} \over 2}} \right)}}T_{i + 1}^2(x), where $p (x_{i}) = \sin (\frac{x - x_{i}}{2})$ p\left( {{x_i}} \right) = \sin \left( {{{x - {x_i}} \over 2}} \right). By performing the necessary operations, we have $T_{i}^{3} (x) = \frac{1}{θ} {\begin{array}{l} p^{3} (x_{i}), & x_{i} \leq x < x_{i + 1} \\ - p^{2} (x_{i}) p (x_{i + 2}) - p (x_{i}) p (x_{i + 3}) p (x_{i + 1}) - p (x_{i + 4}) p^{2} (x_{i + 1}), & x_{i + 1} \leq x < x_{i + 2} \\ p (x_{i}) p^{2} (x_{i + 3}) + p (x_{i + 4}) p (x_{i + 3}) p (x_{i + 1}) + p^{2} (x_{i + 4}) p (x_{i + 2}), & x_{i + 2} \leq x < x_{i + 3} \\ - p^{3} (x_{i + 4}) & x_{i + 3} \leq x < x_{i + 4}, \\ 0 & other cases, \end{array}$ T_i^3(x) = {1 \over \theta }\left\{ {\matrix{ {{p^3}\left( {{x_i}} \right),} \hfill & {{x_i} \le x < {x_{i + 1}}} \hfill \cr { - {p^2}\left( {{x_i}} \right)p\left( {{x_{i + 2}}} \right) - p\left( {{x_i}} \right)p\left( {{x_{i + 3}}} \right)p\left( {{x_{i + 1}}} \right) - p\left( {{x_{i + 4}}} \right){p^2}\left( {{x_{i + 1}}} \right),} \hfill & {{x_{i + 1}} \le x < {x_{i + 2}}} \hfill \cr {p\left( {{x_i}} \right){p^2}\left( {{x_{i + 3}}} \right) + p\left( {{x_{i + 4}}} \right)p\left( {{x_{i + 3}}} \right)p\left( {{x_{i + 1}}} \right) + {p^2}\left( {{x_{i + 4}}} \right)p\left( {{x_{i + 2}}} \right),} \hfill & {{x_{i + 2}} \le x < {x_{i + 3}}} \hfill \cr { - {p^3}\left( {{x_{i + 4}}} \right)} \hfill & {{x_{i + 3}} \le x < {x_{i + 4}},} \hfill \cr 0 \hfill & {{\rm{ other cases, }}} \hfill \cr } } \right. where $θ = \sin (\frac{h}{2}) \sin (h) \sin (\frac{3 h}{2})$ \theta = \sin \left( {{h \over 2}} \right)\sin (h)\sin \left( {{{3h} \over 2}} \right). Generally, an approximate solution of a given problem using cubic trigonometric cubic B-spline functions can be defined as 5 $u (x, t) \approx S (x, t) = \sum_{i = - 1}^{N + 1} T_{i}^{3} (x) δ_{i} (t) .$ u(x,t) \approx S(x,t) = \sum\limits_{i = - 1}^{N + 1} {T_i^3} (x){\delta _i}(t).

Here, the u(x,t) symbolize the exact solution of the problem and the coefficients δ_i are time-dependent variables, $T_{i}^{3} (x)$ T_i^3(x) functions represent trigonometric cubic B-spline functions. Additionally, the approximate solution and its first and second order derivatives at nodes x_j are 6 $\begin{array}{l} S (x, t) = a_{1} δ_{j - 1} + a_{2} δ_{j} + a_{3} δ_{j + 1} \\ S^{'} (x, t) = b_{1} δ_{j - 1} + b_{2} δ_{j + 1} \\ S^{″} (x, t) = c_{1} δ_{j - 1} + c_{2} δ_{j} + c_{3} δ_{j + 1} . \end{array}$ \matrix{ {S(x,t) = {a_1}{\delta _{j - 1}} + {a_2}{\delta _j} + {a_3}{\delta _{j + 1}}} \hfill \cr {S'(x,t) = {b_1}{\delta _{j - 1}} + {b_2}{\delta _{j + 1}}} \hfill \cr {S''(x,t) = {c_1}{\delta _{j - 1}} + {c_2}{\delta _j} + {c_3}{\delta _{j + 1}}.} \hfill \cr }

Here; 7 $\begin{array}{l} a_{1} = \sin^{2} (\frac{h}{2}) \csc (h) \csc (\frac{3 h}{2}), a_{2} = \frac{2}{1 + 2 \cos (h)}, a_{1} = a_{3} \\ b_{1} = - \frac{3 \csc (\frac{3 h}{2})}{4}, b_{2} = \frac{3 \csc (\frac{3 h}{2})}{4} \\ c_{1} = \frac{3 (3 \cos^{2} (\frac{h}{2}) - 1)}{4 (\sin (h) \sin (\frac{3 h}{2}))}, c_{2} = - \frac{3 \cot^{2} (\frac{h}{2})}{(2 + 4 \cos (h))}, c_{3} = c_{1} . \end{array}$ \matrix{ {{a_1} = {{\sin }^2}\left( {{h \over 2}} \right)\csc (h)\csc \left( {{{3h} \over 2}} \right),{a_2} = {2 \over {1 + 2\cos (h)}},\,\,\,\,{a_1} = {a_3}} \hfill \cr {{b_1} = - {{3\csc \left( {{{3h} \over 2}} \right)} \over 4},\;\;\;{\mkern 1mu} {b_2} = {{3\csc \left( {{{3h} \over 2}} \right)} \over 4}} \hfill \cr {{c_1} = {{3\left( {3{{\cos }^2}\left( {{h \over 2}} \right) - 1} \right)} \over {4\left( {\sin (h)\sin \left( {{{3h} \over 2}} \right)} \right)}},{c_2} = - {{3{{\cot }^2}\left( {{h \over 2}} \right)} \over {(2 + 4\cos (h))}},{c_3} = {c_1}.} \hfill \cr }

2.2

Higher-order approach

As defined in the previous subsection, the approximate solution is given. Now, in this subsection, a higherorder nodal approximation of the solution and its derivatives is presented. Suppose that S(x,t) satisfies the following interpolation and boundary conditions for the collocation nodes x_j of the uniform spatial grid, with j = 0,1,···,N.8 $\begin{array}{l} S (x_{j}, t) = u (x_{j}, t), & 0 \leq j \leq N \\ S^{″} (x_{j}, t) = u^{″} (x_{j}, t) - \frac{1}{12} h^{2} u^{(4)} (x_{j}, t), & j = 0, N . \end{array}$ \[\begin{array}{*{35}{l}} S\left( {{x}_{j}},t \right)=u\left( {{x}_{j}},t \right), & 0\le j\le N \\ {S}''\left( {{x}_{j}},t \right)={u}''\left( {{x}_{j}},t \right)-\frac{1}{12}{{h}^{2}}{{u}^{(4)}}\left( {{x}_{j}},t \right), & j=0,N. \\ \end{array}\]

If u(x,t) is sufficiently smooth function in [a,b] and S(x,t) satisfies the above conditions, then the following relations hold [26,27] 9 $\begin{array}{l} S^{'} (x_{j}, t) = u^{'} (x_{j}, t) + O (h^{4}), & 0 \leq j \leq N \\ S^{″} (x_{j}, t) = u^{″} (x_{j}, t) - \frac{1}{12} h^{2} u^{(4)} (x_{j}, t), + O (h^{4}) & 0 \leq j \leq N . \end{array}$ \[\begin{array}{*{35}{l}} {S}'\left( {{x}_{j}},t \right)={u}'\left( {{x}_{j}},t \right)+O\left( {{h}^{4}} \right), & 0\le j\le N \\ {S}''\left( {{x}_{j}},t \right)={u}''\left( {{x}_{j}},t \right)-\frac{1}{12}{{h}^{2}}{{u}^{(4)}}\left( {{x}_{j}},t \right),+O\left( {{h}^{4}} \right) & 0\le j\le N. \\ \end{array}\]

Using finite differences and Taylor expansion, the approximate solution and its derivatives at the nodes can be defined as following way 10 $\begin{array}{l} S (x_{j}) = u (x_{j}), & 0 \leq j \leq N \\ S^{″} (x_{j}) = u^{″} (x_{j}), & j = 0, N \end{array}$ \matrix{ {S\left( {{x_j}} \right) = u\left( {{x_j}} \right),} \hfill & {0 \le j \le N} \hfill \cr {S''\left( {{x_j}} \right) = u''\left( {{x_j}} \right),} \hfill & {j = 0,N} \hfill \cr } and 11 $u^{(4)} (x_{j}) = \frac{S^{″} (x_{j - 1}) - 2 S^{″} (x_{j}) + S^{″} (x_{j + 1})}{h^{2}} + O (h^{2}) .$ {u^{(4)}}\left( {{x_j}} \right) = {{S''\left( {{x_{j - 1}}} \right) - 2S''\left( {{x_j}} \right) + S''\left( {{x_{j + 1}}} \right)} \over {{h^2}}} + O\left( {{h^2}} \right).

Here, one can obtain 12 $\begin{array}{l} j = 1, u^{(4)} (x_{1}) = \frac{S^{″} (x_{0}) - 2 S^{″} (x_{1}) + S^{″} (x_{2})}{h^{2}} + O (h^{2}), \\ j = 2, u^{(4)} (x_{2}) = \frac{S^{″} (x_{1}) - 2 S^{″} (x_{2}) + S^{″} (x_{3})}{h^{2}} + O (h^{2}), \end{array}$ \matrix{ {j = 1,\;\;\;{\mkern 1mu} {u^{(4)}}\left( {{x_1}} \right) = {{S''\left( {{x_0}} \right) - 2S''\left( {{x_1}} \right) + S''\left( {{x_2}} \right)} \over {{h^2}}} + O\left( {{h^2}} \right),} \hfill \cr {j = 2,\;\;\;{\mkern 1mu} {u^{(4)}}\left( {{x_2}} \right) = {{S''\left( {{x_1}} \right) - 2S''\left( {{x_2}} \right) + S''\left( {{x_3}} \right)} \over {{h^2}}} + O\left( {{h^2}} \right),} \hfill \cr } and 13 $u^{(4)} (x_{0}) = 2 u^{(4)} (x_{1}) - u^{(4)} (x_{2}) .$ {u^{(4)}}\left( {{x_0}} \right) = 2{u^{(4)}}\left( {{x_1}} \right) - {u^{(4)}}\left( {{x_2}} \right).

When we use the equalities in (12) in the equation given in (13), it can easily be obtained 14 $\begin{array}{l} u^{(4)} (x_{0}) & = \frac{2 S^{″} (x_{0}) - 4 S^{″} (x_{1}) + 2 S^{″} (x_{2}) - S^{″} (x_{1}) + 2 S^{″} (x_{2}) - S^{″} (x_{3})}{h^{2}} \\ = \frac{2 S^{″} (x_{0}) - 5 S^{″} (x_{1}) + 4 S^{″} (x_{2}) - S^{″} (x_{3})}{h^{2}} + O (h^{2}) . \end{array}$ \matrix{ {{u^{(4)}}\left( {{x_0}} \right)} \hfill & { = {{2S''\left( {{x_0}} \right) - 4S''\left( {{x_1}} \right) + 2S''\left( {{x_2}} \right) - S''\left( {{x_1}} \right) + 2S''\left( {{x_2}} \right) - S''\left( {{x_3}} \right)} \over {{h^2}}}} \hfill \cr {} \hfill & { = {{2S''\left( {{x_0}} \right) - 5S''\left( {{x_1}} \right) + 4S''\left( {{x_2}} \right) - S''\left( {{x_3}} \right)} \over {{h^2}}} + O\left( {{h^2}} \right).} \hfill \cr }

Similarly, using (11), it yields 15 $\begin{array}{l} j = N - 1, & u^{(4)} (x_{N - 1}) = \frac{S^{″} (x_{N - 2}) - 2 S^{″} (x_{N - 1}) + S^{″} (x_{N})}{h^{2}}, \\ j = N - 2, & u^{(4)} (x_{N - 2}) = \frac{S^{″} (x_{N - 3}) - 2 S^{″} (x_{N - 2}) + S^{″} (x_{N - 1})}{h^{2}}, \end{array}$ \matrix{ {j = N - 1,} \hfill & {{u^{(4)}}\left( {{x_{N - 1}}} \right) = {{S''\left( {{x_{N - 2}}} \right) - 2S''\left( {{x_{N - 1}}} \right) + S''\left( {{x_N}} \right)} \over {{h^2}}},} \hfill \cr {j = N - 2,} \hfill & {{u^{(4)}}\left( {{x_{N - 2}}} \right) = {{S''\left( {{x_{N - 3}}} \right) - 2S''\left( {{x_{N - 2}}} \right) + S''\left( {{x_{N - 1}}} \right)} \over {{h^2}}},} \hfill \cr } and using this equalities for 16 $u^{(4)} (x_{N}) = 2 u^{(4)} (x_{N - 1}) - u^{(4)} (x_{N - 2}) .$ {u^{(4)}}\left( {{x_N}} \right) = 2{u^{(4)}}\left( {{x_{N - 1}}} \right) - {u^{(4)}}\left( {{x_{N - 2}}} \right).

We get 17 $u^{(4)} (x_{N}) = \frac{2 S^{″} (x_{N}) - 5 S^{″} (x_{N - 1}) + 4 S^{″} (x_{N - 2}) - S^{″} (x_{N - 3})}{h^{2}} + O (h^{2}) .$ {u^{(4)}}\left( {{x_N}} \right) = {{2S''\left( {{x_N}} \right) - 5S''\left( {{x_{N - 1}}} \right) + 4S''\left( {{x_{N - 2}}} \right) - S''\left( {{x_{N - 3}}} \right)} \over {{h^2}}} + O\left( {{h^2}} \right).

Lastly, we can obtain as following approximations for fourth order derivative; 18 $\begin{array}{l} u^{(4)} (x_{j}, t) = \frac{2 S^{″} (x_{0}, t) - 5 S^{″} (x_{1}, t) + 4 S^{″} (x_{2}, t) - S^{″} (x_{3}, t)}{h^{2}} + O (h^{2}), & j = 0 \\ u^{(4)} (x_{j}, t) = \frac{S^{″} (x_{j - 1}, t) - 2 S^{″} (x_{j}, t) + S^{″} (x_{j + 1}, t)}{h^{2}} + O (h^{2}), & 1 \leq j \leq N - 1 \\ u^{(4)} (x_{j}, t) = \frac{2 S^{″} (x_{N}, t) - 5 S^{″} (x_{N - 1}, t) + 4 S^{″} (x_{N - 2}, t) - S^{″} (x_{N - 3}, t)}{h^{2}} + O (h^{2}), & j = N . \end{array}$ \matrix{ {{u^{(4)}}\left( {{x_j},t} \right) = {{2S''\left( {{x_0},t} \right) - 5S''\left( {{x_1},t} \right) + 4S''\left( {{x_2},t} \right) - S''\left( {{x_3},t} \right)} \over {{h^2}}} + O\left( {{h^2}} \right),} \hfill & {j = 0} \hfill \cr {{u^{(4)}}\left( {{x_j},t} \right) = {{S''\left( {{x_{j - 1}},t} \right) - 2S''\left( {{x_j},t} \right) + S''\left( {{x_{j + 1}},t} \right)} \over {{h^2}}} + O\left( {{h^2}} \right),} \hfill & {1 \le j \le N - 1} \hfill \cr {{u^{(4)}}\left( {{x_j},t} \right) = {{2S''\left( {{x_N},t} \right) - 5S''\left( {{x_{N - 1}},t} \right) + 4S''\left( {{x_{N - 2}},t} \right) - S''\left( {{x_{N - 3}},t} \right)} \over {{h^2}}} + O\left( {{h^2}} \right),} \hfill & {j = N.} \hfill \cr }

Using equations given by (9), we get following approximations for seconds order derivative; 19 $\begin{array}{l} u^{″} (x_{0}, t) = \frac{14 S^{″} (x_{0}, t) - 5 S^{″} (x_{1}, t) + 4 S^{″} (x_{2}, t) - S^{″} (x_{3}, t)}{12} + O (h^{4}), & j = 0 \\ u^{″} (x_{j}, t) = \frac{S^{″} (x_{j - 1}, t) + 10 S^{″} (x_{j}, t) + S^{″} (x_{j + 1}, t)}{12} + O (h^{4}), & 1 \leq j \leq N - 1 \\ u^{″} (x_{N}, t) = \frac{14 S^{″} (x_{N}, t) - 5 S^{″} (x_{N - 1}, t) + 4 S^{″} (x_{N - 2}, t) - 5 S^{″} (x_{N - 3}, t)}{12} + O (h^{4}), & j = N . \end{array}$ \matrix{ {u''\left( {{x_0},t} \right) = {{14S''\left( {{x_0},t} \right) - 5S''\left( {{x_1},t} \right) + 4S''\left( {{x_2},t} \right) - S''\left( {{x_3},t} \right)} \over {12}} + O\left( {{h^4}} \right),} \hfill & {j = 0} \hfill \cr {u''\left( {{x_j},t} \right) = {{S''\left( {{x_{j - 1}},t} \right) + 10S''\left( {{x_j},t} \right) + S''\left( {{x_{j + 1}},t} \right)} \over {12}} + O\left( {{h^4}} \right),} \hfill & {1 \le j \le N - 1} \hfill \cr {u''\left( {{x_N},t} \right) = {{14S''\left( {{x_N},t} \right) - 5S''\left( {{x_{N - 1}},t} \right) + 4S''\left( {{x_{N - 2}},t} \right) - 5S''\left( {{x_{N - 3}},t} \right)} \over {12}} + O\left( {{h^4}} \right),} \hfill & {j = N.} \hfill \cr }

Application of the method to time-fractional Burgers equation

This section presents the application of a higher-order trigonometric cubic B-spline collocation method to solve time-fractional Burgers equation with initial and boundary conditions. Let us reconsider the equation with related boundary and initial condition given in (1)-(2). The numerical implementation of the method to aforementioned equation can be decomposed into 5 key steps. The first one is discretization of the equation employing Crank-Nicolson approach for spatial derivatives and L1 algorithm for the time fractional order derivative given in [28]. When we apply the first step, we have the following equation 20 $\begin{array}{l} \frac{{(Δ t)}^{- γ}}{Γ (2 - γ)} \sum_{k = 0}^{n - 1} b_{k} [u^{n - k} - u^{n - k - 1}] + \frac{{(u u_{x})}^{n + 1} + {(u u_{x})}^{n}}{2} - v \frac{(u_{x x}^{n + 1}) + (u_{x x}^{n})}{2} = \frac{f (x_{j}, t_{n + 1}) + f (x_{j}, t_{n})}{2}, & 1 \leq j \leq N - 1. \end{array}$ \matrix{ {{{{{(\Delta t)}^{ - \gamma }}} \over {\Gamma (2 - \gamma )}}\sum\limits_{k = 0}^{n - 1} {{b_k}} \left[ {{u^{n - k}} - {u^{n - k - 1}}} \right] + {{{{\left( {u{u_x}} \right)}^{n + 1}} + {{\left( {u{u_x}} \right)}^n}} \over 2} - v{{\left( {u_{xx}^{n + 1}} \right) + \left( {u_{xx}^n} \right)} \over 2} = {{f\left( {{x_j},{t_{n + 1}}} \right) + f\left( {{x_j},{t_n}} \right)} \over 2},} \hfill & {1 \le j \le N - 1.} \hfill \cr }

The source term f(x,t) is discretized in time using the Crank-Nicolson scheme and evaluated at the collocation points x_j. As the second step, the nonlinear term (uu_x)⁽ⁿ⁺¹⁾ is linearized using the Rubin-Graves type technique according to the following formula [29,30] 21 ${(u^{p} u_{x})}^{n + 1} = {(u^{p})}^{n} u_{x}^{n + 1} + p {(u^{p - 1})}^{n} u_{x}^{n} u^{n + 1} - p {(u^{p})}^{n} u_{x}^{n} .$ {\left( {{u^p}{u_x}} \right)^{n + 1}} = {\left( {{u^p}} \right)^n}u_x^{n + 1} + p{\left( {{u^{p - 1}}} \right)^n}u_x^n{u^{n + 1}} - p{\left( {{u^p}} \right)^n}u_x^n.

For p = 1, the Rubin–Graves linearization yields 22 ${(u u_{x})}^{n + 1} \approx u^{n} u_{x}^{n + 1} + u_{x}^{n} u^{n + 1} - u^{n} u_{x}^{n} .$ {\left( {u{u_x}} \right)^{n + 1}} \approx {u^n}u_x^{n + 1} + u_x^n{u^{n + 1}} - {u^n}u_x^n.

Substituting (22) into (20) leads to 23 $\begin{array}{l} u^{n + 1} (1 + \frac{S ψ_{1}}{2}) + u_{x}^{n + 1} (\frac{S ψ_{2}}{2}) - u_{x x}^{n + 1} (\frac{S v}{2}) = u^{n} + u_{x x}^{n} (\frac{S v}{2}) \\ + \frac{S}{2} (f (x, t_{n + 1}) + f (x, t_{n})) - \sum_{k = 1}^{n} b_{k} [u^{n + 1 - k} - u^{n - k}] . \end{array}$ \matrix{ {{u^{n + 1}}\left( {1 + {{S{\psi _1}} \over 2}} \right) + u_x^{n + 1}\left( {{{S{\psi _2}} \over 2}} \right) - u_{xx}^{n + 1}\left( {{{Sv} \over 2}} \right) = {u^n} + u_{xx}^n\left( {{{Sv} \over 2}} \right)} \hfill \cr {\;\;\;{\mkern 1mu} \,\,\,\,\,\,\, + {S \over 2}\left( {f\left( {x,{t_{n + 1}}} \right) + f\left( {x,{t_n}} \right)} \right) - \sum\nolimits_{k = 1}^n {{b_k}} \left[ {{u^{n + 1 - k}} - {u^{n - k}}} \right].} \hfill \cr }

Here, $S = Δ t^{γ} Γ (2 - γ), ψ_{1} = u_{x}^{n}$ S = \Delta {t^\gamma }\Gamma (2 - \gamma ),\,\,{\psi _1} = u_x^n and ψ² = uⁿ. Accordingly, the coefficients multiplying u⁽ⁿ⁺¹⁾ and $u_{x}^{n + 1}$ u_x^{n + 1} are $u_{x}^{n}$ u_x^n and uⁿ, respectively. Now, the third key step is using approximate solutions in (5) into the discretized equation given in (23). In the light of previous section, and using approximate solutions, the third step of the numerical implementation is resulted as following algebraic equation; 24 $j = 0 : {\begin{array}{l} δ_{- 1}^{n + 1} [a_{1} (1 + \frac{S}{2} ψ_{1}) + \frac{S b_{1}}{2} ψ_{2} - \frac{S v}{24} (14 c_{1})] + δ_{0}^{n + 1} [a_{2} (1 + \frac{S}{2} ψ_{1}) - \frac{S v}{24} (14 c_{2} - 5 c_{1})] \\ + δ_{1}^{n + 1} [a_{3} (1 + \frac{S}{2} ψ_{1}) + \frac{S b_{2}}{2} ψ_{2} - \frac{S v}{24} (14 c_{3} - 5 c_{2} + 4 c_{1})] + δ_{2}^{n + 1} [- \frac{S v}{24} (- 5 c_{3} + 4 c_{2} - c_{1})] \\ + δ_{3}^{n + 1} [- \frac{S v}{24} (4 c_{3} - c_{2})] + δ_{4}^{n + 1} [- \frac{S v}{24} (- c_{3})] \\ = δ_{- 1}^{n} [a_{1} + \frac{S v}{24} (14 c_{1})] + δ_{0}^{n} [a_{2} + \frac{S v}{24} (14 c_{2} - 5 c_{1})] + δ_{1}^{n} [a_{3} + \frac{S v}{24} (14 c_{3} - 5 c_{2} + 4 c_{1})] \\ + δ_{2}^{n} [\frac{S v}{24} (- 5 c_{3} + 4 c_{2} - c_{1})] + δ_{3}^{n} [\frac{S v}{24} (4 c_{3} - c_{2})] + δ_{4}^{n} [\frac{S v}{24} (- c_{3})] + \frac{S}{2} (f (x_{0}, t_{n + 1}) + f (x_{0}, t_{n})) \\ - \sum_{k = 1}^{n} b_{k} [(a_{1} δ_{- 1}^{n + 1 - k} + a_{2} δ_{0}^{n + 1 - k} + a_{3} δ_{1}^{n + 1 - k}) - (a_{1} δ_{- 1}^{n - k} + a_{2} δ_{0}^{n - k} + a_{3} δ_{1}^{n - k})], \end{array}$ j = 0:\left\{ {\matrix{ {\delta _{ - 1}^{n + 1}\left[ {{a_1}\left( {1 + {S \over 2}{\psi _1}} \right) + {{S{b_1}} \over 2}{\psi _2} - {{Sv} \over {24}}\left( {14{c_1}} \right)} \right] + \delta _0^{n + 1}\left[ {{a_2}\left( {1 + {S \over 2}{\psi _1}} \right) - {{Sv} \over {24}}\left( {14{c_2} - 5{c_1}} \right)} \right]} \hfill \cr { + \delta _1^{n + 1}\left[ {{a_3}\left( {1 + {S \over 2}{\psi _1}} \right) + {{S{b_2}} \over 2}{\psi _2} - {{Sv} \over {24}}\left( {14{c_3} - 5{c_2} + 4{c_1}} \right)} \right] + \delta _2^{n + 1}\left[ { - {{Sv} \over {24}}\left( { - 5{c_3} + 4{c_2} - {c_1}} \right)} \right]} \hfill \cr { + \delta _3^{n + 1}\left[ { - {{Sv} \over {24}}\left( {4{c_3} - {c_2}} \right)} \right] + \delta _4^{n + 1}\left[ { - {{Sv} \over {24}}\left( { - {c_3}} \right)} \right]} \hfill \cr {\;\;\;{\mkern 1mu} \,\,\,\,\,\,\, = \delta _{ - 1}^n\left[ {{a_1} + {{Sv} \over {24}}\left( {14{c_1}} \right)} \right] + \delta _0^n\left[ {{a_2} + {{Sv} \over {24}}\left( {14{c_2} - 5{c_1}} \right)} \right] + \delta _1^n\left[ {{a_3} + {{Sv} \over {24}}\left( {14{c_3} - 5{c_2} + 4{c_1}} \right)} \right]} \hfill \cr {\;\;\;\,\,\,\,\,\,{\mkern 1mu} + \delta _2^n\left[ {{{Sv} \over {24}}\left( { - 5{c_3} + 4{c_2} - {c_1}} \right)} \right] + \delta _3^n\left[ {{{Sv} \over {24}}\left( {4{c_3} - {c_2}} \right)} \right] + \delta _4^n\left[ {{{Sv} \over {24}}\left( { - {c_3}} \right)} \right] + {S \over 2}\left( {f\left( {{x_0},{t_{n + 1}}} \right) + f\left( {{x_0},{t_n}} \right)} \right)} \hfill \cr {\;\;\;\,\,\,\,\,\,{\mkern 1mu} - \sum\limits_{k = 1}^n {{b_k}} \left[ {\left( {{a_1}\delta _{ - 1}^{n + 1 - k} + {a_2}\delta _0^{n + 1 - k} + {a_3}\delta _1^{n + 1 - k}} \right) - \left( {{a_1}\delta _{ - 1}^{n - k} + {a_2}\delta _0^{n - k} + {a_3}\delta _1^{n - k}} \right)} \right],} \hfill \cr } } \right. and 25 $j = 1, 2, .., N - 1 : {\begin{array}{l} δ_{j - 2}^{n + 1} [- \frac{S v}{24} c_{1}] + δ_{j - 1}^{n + 1} [a_{1} (1 + \frac{S}{2} ψ_{1}) + \frac{S b_{1}}{2} ψ_{2} - \frac{S v}{24} (c_{2} + 10 c_{1})] \\ + δ_{j}^{n + 1} [a_{2} (1 + \frac{S}{2} ψ_{1}) - \frac{S v}{24} (c_{3} + 10 c_{2} + c_{1})] \\ + δ_{j + 1}^{n + 1} [a_{3} (1 + \frac{S}{2} ψ_{1}) + \frac{S b_{2}}{2} ψ_{2} - \frac{S v}{24} (10 c_{3} + c_{2})] + δ_{j + 2}^{n + 1} [- \frac{S v}{24} c_{3}] \\ = δ_{j - 2}^{n} [\frac{S v}{24} c_{1}] + δ_{j - 1}^{n} [a_{1} + \frac{S v}{24} (c_{2} + 10 c_{1})] + δ_{j}^{n} [a_{2} + \frac{S v}{24} (c_{3} + 10 c_{2} + c_{1})] \\ δ_{j + 1}^{n} [a_{3} + \frac{S v}{24} (10 c_{3} + c_{2})] + δ_{j + 2}^{n} [\frac{S v}{24} (c_{3})] + \frac{S}{2} (f (x_{j}, t_{n + 1}) + f (x_{j}, t_{n})) \\ - \sum_{k = 1}^{n} b_{k} [(a_{1} δ_{m - 1}^{n + 1 - k} + a_{2} δ_{m}^{n + 1 - k} + a_{3} δ_{m + 1}^{n + 1 - k}) - (a_{1} δ_{m - 1}^{n - k} + a_{2} δ_{m}^{n - k} + a_{3} δ_{m + 1}^{n - k})], \end{array}$ j = 1,2,..,N - 1:\left\{ {\matrix{ {\delta _{j - 2}^{n + 1}\left[ { - {{Sv} \over {24}}{c_1}} \right] + \delta _{j - 1}^{n + 1}\left[ {{a_1}\left( {1 + {S \over 2}{\psi _1}} \right) + {{S{b_1}} \over 2}{\psi _2} - {{Sv} \over {24}}\left( {{c_2} + 10{c_1}} \right)} \right]} \hfill \cr { + \delta _j^{n + 1}\left[ {{a_2}\left( {1 + {S \over 2}{\psi _1}} \right) - {{Sv} \over {24}}\left( {{c_3} + 10{c_2} + {c_1}} \right)} \right]} \hfill \cr { + \delta _{j + 1}^{n + 1}\left[ {{a_3}\left( {1 + {S \over 2}{\psi _1}} \right) + {{S{b_2}} \over 2}{\psi _2} - {{Sv} \over {24}}\left( {10{c_3} + {c_2}} \right)} \right] + \delta _{j + 2}^{n + 1}\left[ { - {{Sv} \over {24}}{c_3}} \right]} \hfill \cr {\;\;\;{\mkern 1mu} \,\,\,\,\,\,\, = \delta _{j - 2}^n\left[ {{{Sv} \over {24}}{c_1}} \right] + \delta _{j - 1}^n\left[ {{a_1} + {{Sv} \over {24}}\left( {{c_2} + 10{c_1}} \right)} \right] + \delta _j^n\left[ {{a_2} + {{Sv} \over {24}}\left( {{c_3} + 10{c_2} + {c_1}} \right)} \right]} \hfill \cr {\;\;\;{\mkern 1mu} \,\,\,\,\,\,\delta _{j + 1}^n\left[ {{a_3} + {{Sv} \over {24}}\left( {10{c_3} + {c_2}} \right)} \right] + \delta _{j + 2}^n\left[ {{{Sv} \over {24}}\left( {{c_3}} \right)} \right] + {S \over 2}\left( {f\left( {{x_j},{t_{n + 1}}} \right) + f\left( {{x_j},{t_n}} \right)} \right)} \hfill \cr {\;\;\;{\mkern 1mu} \,\,\,\,\, - \sum\limits_{k = 1}^n {{b_k}} \left[ {\left( {{a_1}\delta _{m - 1}^{n + 1 - k} + {a_2}\delta _m^{n + 1 - k} + {a_3}\delta _{m + 1}^{n + 1 - k}} \right) - \left( {{a_1}\delta _{m - 1}^{n - k} + {a_2}\delta _m^{n - k} + {a_3}\delta _{m + 1}^{n - k}} \right)} \right],} \hfill \cr } ,} \right. and similarly, 26 $j = N : {\begin{array}{l} δ_{N - 4}^{n + 1} [\frac{S v}{24} c_{1}] + δ_{N - 3}^{n + 1} [- \frac{S v}{24} (- c_{2} + 4 c_{1})] + δ_{N - 2}^{n + 1} [- \frac{S v}{24} (- c_{3} + 4 c_{2} - 5 c_{1})] \\ + δ_{N - 1}^{n + 1} [a_{1} (1 + \frac{S}{2} ψ_{1}) + \frac{S b_{1}}{2} ψ_{2} - \frac{S v}{24} (4 c_{3} - 5 c_{2} + 14 c_{1})] + δ_{N}^{n + 1} [a_{2} (1 + \frac{S}{2} ψ_{1}) - \frac{S v}{24} (- 5 c_{3} + 14 c_{2})] \\ + δ_{N + 1}^{n + 1} [a_{3} (1 + \frac{S}{2} ψ_{1}) + \frac{S b_{2}}{2} ψ_{2} - \frac{S v}{24} (14 c_{3})] \\ = δ_{N - 4}^{n} [- \frac{S v}{24} c_{1}] + δ_{N - 3}^{n} [\frac{S v}{24} (- c_{2} + 4 c_{1})] δ_{N - 2}^{n} [\frac{S v}{24} (- c_{3} + 4 c_{2} - 5 c_{1})] \\ + δ_{N - 1}^{n} [a_{1} + \frac{S v}{24} (4 c_{3} - 5 c_{2} + 14 c_{1})] + δ_{N}^{n} [a_{2} + \frac{S v}{24} (- 5 c_{3} + 14 c_{2})] \\ + δ_{N + 1}^{n} [a_{3} + \frac{S v}{24} (14 c_{3})] + \frac{S}{2} (f (x_{N}, t_{n + 1}) + f (x_{N}, t_{n})) \\ - \sum_{k = 1}^{n} b_{k} [(a_{1} δ_{N - 1}^{n + 1 - k} + a_{2} δ_{N}^{n + 1 - k} + a_{3} δ_{N + 1}^{n + 1 - k}) - (a_{1} δ_{N - 1}^{n - k} + a_{2} δ_{N}^{n - k} + a_{3} δ_{N + 1}^{n - k})] . \end{array}$ j = N:\left\{ {\matrix{ {\delta _{N - 4}^{n + 1}\left[ {{{Sv} \over {24}}{c_1}} \right] + \delta _{N - 3}^{n + 1}\left[ { - {{Sv} \over {24}}\left( { - {c_2} + 4{c_1}} \right)} \right] + \delta _{N - 2}^{n + 1}\left[ { - {{Sv} \over {24}}\left( { - {c_3} + 4{c_2} - 5{c_1}} \right)} \right]} \hfill \cr { + \delta _{N - 1}^{n + 1}\left[ {{a_1}\left( {1 + {S \over 2}{\psi _1}} \right) + {{S{b_1}} \over 2}{\psi _2} - {{Sv} \over {24}}\left( {4{c_3} - 5{c_2} + 14{c_1}} \right)} \right] + \delta _N^{n + 1}\left[ {{a_2}\left( {1 + {S \over 2}{\psi _1}} \right) - {{Sv} \over {24}}\left( { - 5{c_3} + 14{c_2}} \right)} \right]} \hfill \cr { + \delta _{N + 1}^{n + 1}\left[ {{a_3}\left( {1 + {S \over 2}{\psi _1}} \right) + {{S{b_2}} \over 2}{\psi _2} - {{Sv} \over {24}}\left( {14{c_3}} \right)} \right]} \hfill \cr {\;\;\;\,\,\,\,\,\,\,{\mkern 1mu} = \delta _{N - 4}^n\left[ { - {{Sv} \over {24}}{c_1}} \right] + \delta _{N - 3}^n\left[ {{{Sv} \over {24}}\left( { - {c_2} + 4{c_1}} \right)} \right]\delta _{N - 2}^n\left[ {{{Sv} \over {24}}\left( { - {c_3} + 4{c_2} - 5{c_1}} \right)} \right]} \hfill \cr {\;\;\;{\mkern 1mu} \,\,\,\,\,\, + \delta _{N - 1}^n\left[ {{a_1} + {{Sv} \over {24}}\left( {4{c_3} - 5{c_2} + 14{c_1}} \right)} \right] + \delta _N^n\left[ {{a_2} + {{Sv} \over {24}}\left( { - 5{c_3} + 14{c_2}} \right)} \right]} \hfill \cr {\;\;\;{\mkern 1mu} \,\,\,\,\,\, + \delta _{N + 1}^n\left[ {{a_3} + {{Sv} \over {24}}\left( {14{c_3}} \right)} \right] + {S \over 2}\left( {f\left( {{x_N},{t_{n + 1}}} \right) + f\left( {{x_N},{t_n}} \right)} \right)} \hfill \cr {\;\;\;\,\,\,\,\,\,{\mkern 1mu} - \sum\limits_{k = 1}^n {{b_k}} \left[ {\left( {{a_1}\delta _{N - 1}^{n + 1 - k} + {a_2}\delta _N^{n + 1 - k} + {a_3}\delta _{N + 1}^{n + 1 - k}} \right) - \left( {{a_1}\delta _{N - 1}^{n - k} + {a_2}\delta _N^{n - k} + {a_3}\delta _{N + 1}^{n - k}} \right)} \right].} \hfill \cr } } \right.

The algebraic equations given in (24) – (26) produces an algebraic equation system consisting of N + 1 equations and N + 3 unknowns for the values of j = 0, 1, ⋯, N. In order to obtain a solvable system 2 equations must be added, or 2 unknowns must be eliminated from the algebraic equation system. The fourth step of implementation is involves the reducing to system into N + 1 equations and N + 1 unknowns by utilizing the boundary conditions. Eliminating the parameters δ₋₁ and δ_{N + 1} yields an algebraic equation system is to ready to be solved by iteratively. In the sake of seeing clearly, the discretization boundary conditions and removing process can be expressed by following equations; 27 $\begin{array}{l} u (a, t) = a_{1} δ_{- 1} + a_{2} δ_{0} + a_{3} δ_{1} = h_{1} (t) \\ δ_{- 1} = \frac{h_{1} (t)}{a_{1}} - \frac{a_{2}}{a_{1}} δ_{0} - \frac{a_{3}}{a_{1}} δ_{1} \\ u (b, t) = a_{1} δ_{N - 1}^{n + 1 - k} + a_{2} δ_{N}^{n + 1 - k} + a_{3} δ_{N + 1}^{n + 1 - k} = h_{2} (t) \\ δ_{N + 1} = \frac{h_{2} (t)}{a_{3}} - \frac{a_{2}}{a_{3}} δ_{N} - \frac{a_{1}}{a_{3}} δ_{N - 1} . \end{array}$ \matrix{ {u(a,t) = {a_1}{\delta _{ - 1}} + {a_2}{\delta _0} + {a_3}{\delta _1} = {h_1}(t)} \hfill \cr {\,\,\,\,\,\,\,{\delta _{ - 1}} = {{{h_1}(t)} \over {{a_1}}} - {{{a_2}} \over {{a_1}}}{\delta _0} - {{{a_3}} \over {{a_1}}}{\delta _1}} \hfill \cr {u(b,t) = {a_1}\delta _{N - 1}^{n + 1 - k} + {a_2}\delta _N^{n + 1 - k} + {a_3}\delta _{N + 1}^{n + 1 - k} = {h_2}(t)} \hfill \cr {\,\,\,\,\,\,{\delta _{N + 1}} = {{{h_2}(t)} \over {{a_3}}} - {{{a_2}} \over {{a_3}}}{\delta _N} - {{{a_1}} \over {{a_3}}}{\delta _{N - 1}}.} \hfill \cr }

The equation system obtained above can be written in matrix form as 28 $A δ^{n + 1} = B δ^{n} .$ A{\delta ^{n + 1}} = B{\delta ^n}.

As the last step, the initial vector $[δ_{0}^{0}, δ_{1}^{0} δ_{2}^{0}, \dots, δ_{N}^{0}]$ \left[ {\delta _0^0,\delta _1^0\delta _2^0, \cdots ,\delta _N^0} \right] will be constructed using the initial condition u(x,0) = g(x) of the equation and approximate solution. For t = 0 in S(x,t), it follows that 29 $\begin{array}{l} j = 0 & a_{1} δ_{- 1} + a_{2} δ_{0} + a_{3} δ_{1} = g (x_{0}) \\ j = 1 & a_{1} δ_{0} + a_{2} δ_{1} + a_{3} δ_{2} = g (x_{1}) \\ ⋮ & ⋮ \\ j = N & a_{1} δ_{N - 1} + a_{2} δ_{N} + a_{3} δ_{N + 1} = g (x_{N}) \end{array}$ \matrix{ {j = 0} \hfill & {{a_1}{\delta _{ - 1}} + {a_2}{\delta _0} + {a_3}{\delta _1} = g\left( {{x_0}} \right)} \hfill \cr {j = 1} \hfill & {{a_1}{\delta _0} + {a_2}{\delta _1} + {a_3}{\delta _2} = g\left( {{x_1}} \right)} \hfill \cr \vdots \hfill & \vdots \hfill \cr {j = N} \hfill & {{a_1}{\delta _{N - 1}} + {a_2}{\delta _N} + {a_3}{\delta _{N + 1}} = g\left( {{x_N}} \right)} \hfill \cr } the system in (29) is in the form N + 1 equations and N + 3 unknowns. The system given in (29) can be solved after the elimination of the parameters δ₋₁ and δ_{N + 1} like as previous step using the Neumann boundary conditions u_x(a,t) = h₃(t) and u_x(a,t) = h₄(t). The implementation of the method will be completed via obtaining initial vector as follows 30 $\begin{array}{l} \begin{array}{l} j = 0; & a_{2} δ_{0} + (a_{3} - \frac{b_{2}}{b_{1}} a_{1}) δ_{1} = g (x_{0}) - \frac{h_{3} (0)}{b_{1}} a_{1} \end{array} \\ 1 \leq j \leq N - 1; a_{1} δ_{j - 1} + a_{2} δ_{j} + a_{3} δ_{j + 1} = g (x_{j}) \\ \begin{array}{l} j = N; & (a_{1} - \frac{b_{1}}{b_{2}} a_{3}) δ_{N - 1} + a_{2} δ_{N} = g (x_{N}) - \frac{h_{4} (0)}{b_{2}} a_{3} \end{array} \end{array}$ \eqalign{ & \matrix{ {j = 0;} \hfill & {{a_2}{\delta _0} + \left( {{a_3} - {{{b_2}} \over {{b_1}}}{a_1}} \right){\delta _1} = g\left( {{x_0}} \right) - {{{h_3}(0)} \over {{b_1}}}{a_1}} \hfill \cr } \cr & 1 \le j \le N - 1;{a_1}{\delta _{j - 1}} + {a_2}{\delta _j} + {a_3}{\delta _{j + 1}} = g\left( {{x_j}} \right) \cr & \matrix{ {j = N;} \hfill & {\left( {{a_1} - {{{b_1}} \over {{b_2}}}{a_3}} \right){\delta _{N - 1}} + {a_2}{\delta _N} = g\left( {{x_N}} \right) - {{{h_4}(0)} \over {{b_2}}}{a_3}} \hfill \cr } \cr}

Consequently, to obtain numerical solutions, the system given in (24) – (26) is going to be solved by iteratively after implementation of the boundary conditions and iteration will be started by the initial vector which is obtained after solving the system given in (30). More clearly, the parameters $δ_{j}^{n + 1}$ \delta _j^{n + 1} are going to be obtained by the parameters $δ_{j}^{n}$ \delta _j^n at desired time level n. And, using the known values of $δ_{j}^{n + 1}$ \delta _j^{n + 1} in the (5) will produce approximate solutions for time-fractional Burgers equation.

Stability analysis

In this section, we are going to perform the stability analysis of the scheme obtained the application of suggested method for time-fractional Burgers equation using von-Neumann criterion and with the help of mathematical induction. According to the von-Neumann criterion, if the amplification factor |λ|≤1, then the error difference between numerical and exact solutions remains bounded as the number of temporal iteration increases. For the stability analysis, the homogenous part of the equation is considered as used in [31]. Consider the numerical scheme which is obtained using suggested method of related with linearized form of the Eq. (1) by 31 $\begin{array}{l} δ_{j - 2}^{n + 1} (\frac{- S v c_{1}}{24}) + δ_{j - 1}^{n + 1} (a_{1} + \frac{S}{2} (τ b_{1} - \frac{v}{12} (10 c_{1} + c_{2}))) + δ_{j}^{n + 1} (a_{2} - \frac{S v}{12} (c_{1} + 5 c_{2})) \\ + δ_{j + 1}^{n + 1} (a_{1} + \frac{S}{2} (τ b_{2} - \frac{v}{12} (10 c_{1} + c_{2}))) + δ_{j + 2}^{n + 1} (\frac{- S v c_{1}}{24}) \\ = δ_{j - 2}^{n} (\frac{S v c_{1}}{24}) + δ_{j - 1}^{n} (a_{1} - \frac{S}{2} (τ b_{1} - \frac{v}{12} (10 c_{1} + c_{2}))) + δ_{j}^{n} (a_{2} + \frac{S v}{12} (c_{1} + 5 c_{2})) \\ + δ_{j + 1}^{n} (a_{1} - \frac{S}{2} (τ b_{2} - \frac{v}{12} (10 c_{1} + c_{2}))) + δ_{j + 2}^{n} (\frac{S v c_{1}}{24}) \\ - \sum_{k = 0}^{n} b_{k}^{α} ((a_{1} δ_{j - 1}^{n + 1 - k} + a_{1} δ_{j}^{n + 1 - k} + a_{1} δ_{j + 1}^{n + 1 - k}) - (a_{1} δ_{j - 1}^{n - k} + a_{1} δ_{j}^{n - k} + a_{1} δ_{j + 1}^{n - k})) \end{array}$ \matrix{ {\delta _{j - 2}^{n + 1}\left( {{{ - Sv{c_1}} \over {24}}} \right) + \delta _{j - 1}^{n + 1}\left( {{a_1} + {S \over 2}\left( {\tau {b_1} - {v \over {12}}\left( {10{c_1} + {c_2}} \right)} \right)} \right) + \delta _j^{n + 1}\left( {{a_2} - {{Sv} \over {12}}\left( {{c_1} + 5{c_2}} \right)} \right)} \hfill \cr { + \delta _{j + 1}^{n + 1}\left( {{a_1} + {S \over 2}\left( {\tau {b_2} - {v \over {12}}\left( {10{c_1} + {c_2}} \right)} \right)} \right) + \delta _{j + 2}^{n + 1}\left( {{{ - Sv{c_1}} \over {24}}} \right)} \hfill \cr {\,\,\,\,\,\,\,\,\,\,\, = \delta _{j - 2}^n\left( {{{Sv{c_1}} \over {24}}} \right) + \delta _{j - 1}^n\left( {{a_1} - {S \over 2}\left( {\tau {b_1} - {v \over {12}}\left( {10{c_1} + {c_2}} \right)} \right)} \right) + \delta _j^n\left( {{a_2} + {{Sv} \over {12}}\left( {{c_1} + 5{c_2}} \right)} \right)} \hfill \cr {\,\,\,\,\,\,\,\,\,\,\,\,\, + \delta _{j + 1}^n\left( {{a_1} - {S \over 2}\left( {\tau {b_2} - {v \over {12}}\left( {10{c_1} + {c_2}} \right)} \right)} \right) + \delta _{j + 2}^n\left( {{{Sv{c_1}} \over {24}}} \right)} \hfill \cr {\,\,\,\,\,\,\,\,\,\,\,\,\, - \sum\limits_{k = 0}^n {b_k^\alpha } \left( {\left( {{a_1}\delta _{j - 1}^{n + 1 - k} + {a_1}\delta _j^{n + 1 - k} + {a_1}\delta _{j + 1}^{n + 1 - k}} \right) - \left( {{a_1}\delta _{j - 1}^{n - k} + {a_1}\delta _j^{n - k} + {a_1}\delta _{j + 1}^{n - k}} \right)} \right)} \hfill \cr } where $τ = u_{j}^{n}$ \tau = u_j^n. Assume that amplification factor of the Fourier mode is defined as 32 $δ_{j}^{n} = λ^{n} e^{i j h β}$ \delta _j^n = {\lambda ^n}{e^{ijh\beta }} where $i = \sqrt{- 1}$ i = \sqrt { - 1} and β is the number of wavelengths. Using (32) in Eq. (31) and dividing e^ijhβ yields following equation 33 $\begin{array}{l} λ^{n + 1} {(\frac{- S v c_{1}}{24}) (e^{- 2 i h β} + e^{2 i h β}) + (a_{1} - \frac{S v}{24} (10 c_{1} + c_{2})) (e^{- i h β} + e^{i h β}) \\ + a_{2} - \frac{S v}{12} (c_{1} + 5 c_{2}) + \frac{S}{2} τ (b_{1} + b_{2})} \\ = λ^{n} {(\frac{S v c_{1}}{24}) (e^{- 2 i h β} + e^{2 i h β}) + (a_{1} + \frac{S v}{24} (10 c_{1} + c_{2})) (e^{- i h β} + e^{i h β}) \\ + a_{2} + \frac{S v}{12} (c_{1} + 5 c_{2}) - \frac{S}{2} τ (b_{1} + b_{2})} \\ - λ^{n} \sum_{k = 0}^{n} b_{k}^{α} (λ^{1 - k} (a_{2} + a_{1} (e^{- i h β} + e^{i h β})) - λ^{- k} (a_{2} + a_{1} (e^{- i h β} + e^{i h β}))) . \end{array}$ \matrix{ {{\lambda ^{n + 1}}\left\{ {\left( {{{ - Sv{c_1}} \over {24}}} \right)\left( {{e^{ - 2ih\beta }} + {e^{2ih\beta }}} \right) + \left( {{a_1} - {{Sv} \over {24}}\left( {10{c_1} + {c_2}} \right)} \right)\left( {{e^{ - ih\beta }} + {e^{ih\beta }}} \right)} \right.} \hfill \cr {\left. { + {a_2} - {{Sv} \over {12}}\left( {{c_1} + 5{c_2}} \right) + {S \over 2}\tau \left( {{b_1} + {b_2}} \right)} \right\}} \hfill \cr {\,\,\,\,\,\,\,\,\, = {\lambda ^n}\left\{ {\left( {{{Sv{c_1}} \over {24}}} \right)\left( {{e^{ - 2ih\beta }} + {e^{2ih\beta }}} \right) + \left( {{a_1} + {{Sv} \over {24}}\left( {10{c_1} + {c_2}} \right)} \right)\left( {{e^{ - ih\beta }} + {e^{ih\beta }}} \right)} \right.} \hfill \cr {\left. {\,\,\,\,\,\,\,\,\,\,\, + {a_2} + {{Sv} \over {12}}\left( {{c_1} + 5{c_2}} \right) - {S \over 2}\tau \left( {{b_1} + {b_2}} \right)} \right\}} \hfill \cr {\,\,\,\,\,\,\,\,\,\,\, - {\lambda ^n}\sum\limits_{k = 0}^n {b_k^\alpha } \left( {{\lambda ^{1 - k}}\left( {{a_2} + {a_1}\left( {{e^{ - ih\beta }} + {e^{ih\beta }}} \right)} \right) - {\lambda ^{ - k}}\left( {{a_2} + {a_1}\left( {{e^{ - ih\beta }} + {e^{ih\beta }}} \right)} \right)} \right).} \hfill \cr }

We know that $e^{- i h β} + e^{i h β} = 2 \cos (h β), e^{i h β} - e^{- i h β} = 2 i \sin (h β), b_{2} = - b_{1}$ {e^{ - ih\beta }} + {e^{ih\beta }} = 2\cos (h\beta ),\,\,{e^{ih\beta }} - {e^{ - ih\beta }} = 2i\sin (h\beta ),\,\,{b_2} = - {b_1}, grouping the resembling terms, we get 34 $\begin{array}{l} λ^{n + 1} {2 (\frac{- S v c_{1}}{24}) \cos (2 h β) + 2 (a_{1} - \frac{S v}{24} (10 c_{1} + c_{2})) \cos (h β) \\ + a_{2} - \frac{S v}{12} (c_{1} + 5 c_{2}) - i S τ b_{1} \sin (h β)} \\ = λ^{n} {2 (\frac{S v c_{1}}{24}) \cos (2 h β) + 2 (a_{1} + \frac{S v}{24} (10 c_{1} + c_{2})) \cos (h β) \\ + a_{2} + \frac{S v}{12} (c_{1} + 5 c_{2}) + i S τ b_{1} \sin (h β)} \\ - λ^{n} \sum_{k = 0}^{n} b_{k}^{α} (λ^{1 - k} (2 a_{1} \cos (h β) + a_{2}) - λ^{- k} (2 a_{1} \cos (h β) + a_{2})) . \end{array}$ \matrix{ {{\lambda ^{n + 1}}\left\{ {2\left( {{{ - Sv{c_1}} \over {24}}} \right)\cos (2h\beta ) + 2\left( {{a_1} - {{Sv} \over {24}}\left( {10{c_1} + {c_2}} \right)} \right)\cos (h\beta )} \right.} \hfill \cr {\left. { + {a_2} - {{Sv} \over {12}}\left( {{c_1} + 5{c_2}} \right) - iS\tau {b_1}\sin (h\beta )} \right\}} \hfill \cr {\,\,\,\,\,\,\,\,\,\, = {\lambda ^n}\left\{ {2\left( {{{Sv{c_1}} \over {24}}} \right)\cos (2h\beta ) + 2\left( {{a_1} + {{Sv} \over {24}}\left( {10{c_1} + {c_2}} \right)} \right)\cos (h\beta )} \right.} \hfill \cr {\left. {\,\,\,\,\,\,\,\,\,\,\, + {a_2} + {{Sv} \over {12}}\left( {{c_1} + 5{c_2}} \right) + iS\tau {b_1}\sin (h\beta )} \right\}} \hfill \cr {\,\,\,\,\,\,\,\,\,\,\,\, - {\lambda ^n}\sum\limits_{k = 0}^n {b_k^\alpha } \left( {{\lambda ^{1 - k}}\left( {2{a_1}\cos (h\beta ) + {a_2}} \right) - {\lambda ^{ - k}}\left( {2{a_1}\cos (h\beta ) + {a_2}} \right)} \right).} \hfill \cr }

After some arrangements, we have 35 $\begin{array}{l} λ^{n + 1} {2 a_{1} \cos (h β) + a_{2} - \frac{S v}{24} (2 c_{1} \cos (2 h β) + 2 (10 c_{1} + c_{2}) \cos (h β) + (2 c_{1} + 10 c_{2})) - i S τ b_{1} \sin (h β)} \\ = λ^{n} {2 a_{1} \cos (h β) + a_{2} + \frac{S v}{24} (2 c_{1} \cos (2 h β) + 2 (10 c_{1} + c_{2}) \cos (h β) + (2 c_{1} + 10 c_{2})) + i S τ b_{1} \sin (h β)} \\ - λ^{n} \sum_{k = 0}^{n} b_{k}^{α} (λ^{1 - k} (2 a_{1} \cos (h β) + a_{2}) - λ^{- k} (2 a_{1} \cos (h β) + a_{2})) . \end{array}$ \matrix{ {{\lambda ^{n + 1}}\left\{ {2{a_1}\cos (h\beta ) + {a_2} - {{Sv} \over {24}}\left( {2{c_1}\cos (2h\beta ) + 2\left( {10{c_1} + {c_2}} \right)\cos (h\beta ) + \left( {2{c_1} + 10{c_2}} \right)} \right) - iS\tau {b_1}\sin (h\beta )} \right\}} \hfill \cr {\,\,\,\,\,\,\,\,\,\,\,\, = {\lambda ^n}\left\{ {2{a_1}\cos (h\beta ) + {a_2} + {{Sv} \over {24}}\left( {2{c_1}\cos (2h\beta ) + 2\left( {10{c_1} + {c_2}} \right)\cos (h\beta ) + \left( {2{c_1} + 10{c_2}} \right)} \right) + iS\tau {b_1}\sin (h\beta )} \right\}} \hfill \cr {\,\,\,\,\,\,\,\,\,\,\, - {\lambda ^n}\sum\limits_{k = 0}^n {b_k^\alpha } \left( {{\lambda ^{1 - k}}\left( {2{a_1}\cos (h\beta ) + {a_2}} \right) - {\lambda ^{ - k}}\left( {2{a_1}\cos (h\beta ) + {a_2}} \right)} \right).} \hfill \cr }

Without the loss of generality, we can assume that β = 0 as in [32,33], it yields 36 $\begin{array}{l} λ^{n + 1} {2 a_{1} + a_{2} - \frac{S v}{24} (2 c_{1} + 2 (10 c_{1} + c_{2}) + (2 c_{1} + 10 c_{2}))} \\ = λ^{n} {2 a_{1} + a_{2} + \frac{S v}{24} (2 c_{1} + 2 (10 c_{1} + c_{2}) + (2 c_{1} + 10 c_{2}))} \end{array}$ \matrix{ {{\lambda ^{n + 1}}\left\{ {2{a_1} + {a_2} - {{Sv} \over {24}}\left( {2{c_1} + 2\left( {10{c_1} + {c_2}} \right) + \left( {2{c_1} + 10{c_2}} \right)} \right)} \right\}} \hfill \cr {\,\,\,\,\,\,\,\,\,\,\,\, = {\lambda ^n}\left\{ {2{a_1} + {a_2} + {{Sv} \over {24}}\left( {2{c_1} + 2\left( {10{c_1} + {c_2}} \right) + \left( {2{c_1} + 10{c_2}} \right)} \right)} \right\}} \hfill \cr } and 37 $\frac{λ^{n + 1}}{λ^{n}} = \frac{2 a_{1} + a_{2} + \frac{S v}{12} (2 c_{1} + c_{2})}{2 a_{1} + a_{2} - \frac{S v}{12} (2 c_{1} + c_{2})} - \frac{(2 a_{1} + a_{2})}{2 a_{1} + a_{2} - \frac{S v}{12} (2 c_{1} + c_{2})} \sum_{k = 0}^{n} b_{k}^{α} (λ^{1 - k} - λ^{- k}) .$ {{{\lambda ^{n + 1}}} \over {{\lambda ^n}}} = {{2{a_1} + {a_2} + {{Sv} \over {12}}\left( {2{c_1} + {c_2}} \right)} \over {2{a_1} + {a_2} - {{Sv} \over {12}}\left( {2{c_1} + {c_2}} \right)}} - {{\left( {2{a_1} + {a_2}} \right)} \over {2{a_1} + {a_2} - {{Sv} \over {12}}\left( {2{c_1} + {c_2}} \right)}}\sum\limits_{k = 0}^n {b_k^\alpha } \left( {{\lambda ^{1 - k}} - {\lambda ^{ - k}}} \right).

Now, we apply mathematical induction. For n = 0, equation (37) becomes $\frac{λ^{1}}{λ^{0}} = \frac{2 a_{1} + a_{2} + \frac{S v}{12} (2 c_{1} + c_{2})}{2 a_{1} + a_{2} - \frac{S v}{12} (2 c_{1} + c_{2})} .$ {{{\lambda ^1}} \over {{\lambda ^0}}} = {{2{a_1} + {a_2} + {{Sv} \over {12}}\left( {2{c_1} + {c_2}} \right)} \over {2{a_1} + {a_2} - {{Sv} \over {12}}\left( {2{c_1} + {c_2}} \right)}}.

For $| \frac{λ^{n + 1}}{λ^{n}} | \leq 1, | 2 a_{1} + a_{2} + \frac{S v}{12} (2 c_{1} + c_{2}) | \leq | 2 a_{1} + a_{2} - \frac{S v}{12} (2 c_{1} + c_{2}) | and 2 c_{1} + c_{2} \leq 0$ \left| {{{{\lambda ^{n + 1}}} \over {{\lambda ^n}}}} \right| \le 1,\left| {2{a_1} + {a_2} + {{Sv} \over {12}}\left( {2{c_1} + {c_2}} \right)} \right| \le \left| {2{a_1} + {a_2} - {{Sv} \over {12}}\left( {2{c_1} + {c_2}} \right)} \right|{\rm{ and }}2{c_1} + {c_2} \le 0. Using Taylor expansion for c₁ and c₂, we get $\begin{array}{l} c_{1} \approx \frac{1}{h^{2}} - \frac{1}{10} + \frac{113 h^{2}}{16800} - \frac{197 h^{4}}{739200} + O (h^{6}) \\ c_{2} \approx - \frac{2}{h^{2}} + \frac{1}{6} - \frac{29 h^{2}}{1080} + \frac{19 h^{4}}{10080} + O (h^{6}) \end{array}$ \matrix{ {{c_1} \approx {1 \over {{h^2}}} - {1 \over {10}} + {{113{h^2}} \over {16800}} - {{197{h^4}} \over {739200}} + O\left( {{h^6}} \right)} \hfill \cr {{c_2} \approx - {2 \over {{h^2}}} + {1 \over 6} - {{29{h^2}} \over {1080}} + {{19{h^4}} \over {10080}} + O\left( {{h^6}} \right)} \hfill \cr } and it is clear that for sufficiently small $h (i . e ., h \to 0), 2 c_{1} + c_{2} \leq 0$ h(i.e{\rm{., }}h \to 0),\,\,2{c_1} + {c_2} \le 0. Thus, for $n = 0, | λ^{1} | \leq | λ^{0} |$ n = 0,\,\,\left| {{\lambda ^1}} \right| \le \left| {{\lambda ^0}} \right|. Now, assume that $| λ^{1} | \leq | λ^{0} |$ \left|\lambda^{\ell}\right| \leq\left|\lambda^0\right| for $ℓ = 1, 2, \dots, n$ \ell = 1,2, \cdots ,n. We now show that it also holds for $ℓ = n + 1$ \ell = n + 1. Using the assumption $| λ^{1} | \leq | λ^{0} |$ |{\lambda ^1}| \le |{\lambda ^0}| and $b_{k}^{α} = {(k + 1)}^{2 - α} - k^{2 - α} > 0$ b_k^\alpha = {(k + 1)^{2 - \alpha }} - {k^{2 - \alpha }} > 0 for 0 < α < 1, we can conclude from Eq. (37) 38 $| \frac{λ^{n + 1}}{λ^{n}} | \leq | \frac{2 a_{1} + a_{2} + \frac{S v}{12} (2 c_{1} + c_{2})}{2 a_{1} + a_{2} - \frac{S v}{12} (2 c_{1} + c_{2})} | + 2 | λ^{0} | | \frac{(2 a_{1} + a_{2})}{2 a_{1} + a_{2} - \frac{S v}{12} (2 c_{1} + c_{2})} | \sum_{k = 0}^{n} b_{k}^{α} .$ \left| {{{{\lambda ^{n + 1}}} \over {{\lambda ^n}}}} \right| \le \left| {{{2{a_1} + {a_2} + {{Sv} \over {12}}\left( {2{c_1} + {c_2}} \right)} \over {2{a_1} + {a_2} - {{Sv} \over {12}}\left( {2{c_1} + {c_2}} \right)}}} \right| + 2\left| {{\lambda ^0}} \right|\left| {{{\left( {2{a_1} + {a_2}} \right)} \over {2{a_1} + {a_2} - {{Sv} \over {12}}\left( {2{c_1} + {c_2}} \right)}}} \right|\sum\limits_{k = 0}^n {b_k^\alpha } .

Thus, since $2 c_{1} + c_{2} \leq 0$ 2{c_1} + {c_2} \le 0, we get $| λ^{n + 1} | \leq L | λ^{0} |$ \left| {{\lambda ^{n + 1}}} \right| \le L\left| {{\lambda ^0}} \right| where $L = | \frac{2 a_{1} + a_{2} + \frac{S v}{12} (2 c_{1} + c_{2})}{2 a_{1} + a_{2} - \frac{S v}{12} (2 c_{1} + c_{2})} | + 2 | \frac{(2 a_{1} + a_{2})}{2 a_{1} + a_{2} - \frac{S v}{12} (2 c_{1} + c_{2})} | \sum_{k = 0}^{n} b_{k}^{α}$ L = \left| {{{2{a_1} + {a_2} + {{Sv} \over {12}}\left( {2{c_1} + {c_2}} \right)} \over {2{a_1} + {a_2} - {{Sv} \over {12}}\left( {2{c_1} + {c_2}} \right)}}} \right| + 2\left| {{{\left( {2{a_1} + {a_2}} \right)} \over {2{a_1} + {a_2} - {{Sv} \over {12}}\left( {2{c_1} + {c_2}} \right)}}} \right|\sum\limits_{k = 0}^n {b_k^\alpha } . Hence, the proposed higher-order trigonometric cubic B-spline collocation scheme is unconditionally stable in the sense of the von Neumann criterion.

Numerical results

In this section, three numerical examples are presented to validate the accuracy and performance of the proposed higher-order trigonometric cubic B-spline collocation method for the time-fractional Burgers equation. To quantitatively assess the accuracy of the numerical solutions and to compare the obtained results with those reported in the literature, the following error norms are employed: 39 $\begin{array}{l} L_{2} = ‖ u (x, t) - S (x, t) ‖_{2} = \sqrt{h \sum_{j = 1}^{N} {| u (x_{j}, t_{n}) - S (x_{j}, t_{n}) |}^{2}} \\ L_{\infty} = ‖ u (x, t) - S (x, t) ‖_{\infty} = \max_{j} | u (x_{j}, t_{n}) - S (x_{j}, t_{n}) | . \end{array}$ \matrix{ {{L_2} = \,u(x,t) - S(x,t){_2} = \sqrt {h\sum\limits_{j = 1}^N {{{\left| {u\left( {{x_j},{t_n}} \right) - S\left( {{x_j},{t_n}} \right)} \right|}^2}} } } \hfill \cr {{L_\infty } = \,u(x,t) - S(x,t){_\infty } = \mathop {\max }\limits_j \left| {u\left( {{x_j},{t_n}} \right) - S\left( {{x_j},{t_n}} \right)} \right|.} \hfill \cr }

Additionally, the convergence rate of the method is calculated with the formula 40 $R o C = \frac{l n (E r r o r (N_{2}) / E r r o r (N_{1}))}{l n (N_{1} / N_{2})},$ RoC = {{ln\left( {Error\left( {{N_2}} \right)/Error\left( {{N_1}} \right)} \right)} \over {ln\left( {{N_1}/{N_2}} \right)}}, where Error(N₁) and Error(N₂) represent the errors at number of partition N₁ and N₂, respectively. The numerical results are organized into three parts: accuracy tests, comparative studies, and convergence analysis for each problems.

5.1

Example 1

As the first example; consider the time-fractional Burgers equation given by (1) having the related boundary conditions given in below. As the first test problem, we consider the time-fractional Burgers equation given by (1) subject to the following initial and boundary conditions: 41 $\begin{array}{l} u (0, t) = u (1, t) = 0, & t \in [0, T], \\ u_{x} (0, t) = u_{x} (1, t) = 2 π t^{2}, & t \in [0, T] \\ u (x, 0) = 0, & x \in [0, 1] \end{array}$ \matrix{ {u(0,t) = u(1,t) = 0,} \hfill & {t \in [0,T],} \hfill \cr {{u_x}(0,t) = {u_x}(1,t) = 2\pi {t^2},} \hfill & {t \in [0,T]} \hfill \cr {u(x,0) = 0,} \hfill & {x \in [0,1]} \hfill \cr } and the corresponding source term is 42 $f (x, t) = \frac{2 t^{2 - γ} \sin (2 π x)}{Γ (3 - γ)} + 2 π t^{4} \sin (2 π x) \cos (2 π x) + 4 v π^{2} t^{2} \sin (2 π x)$ f(x,t) = {{2{t^{2 - \gamma }}\sin (2\pi x)} \over {\Gamma (3 - \gamma )}} + 2\pi {t^4}\sin (2\pi x)\cos (2\pi x) + 4v{\pi ^2}{t^2}\sin (2\pi x) and exact solution of the problem is $u (x, t) = t^{2} \sin (2 π x)$ u(x,t) = {t^2}\sin (2\pi x) [20]. With the aim of analyzing the first problem, firstly, the time step and viscosity parameter are taken as Δt = 0.0005, v = 1 for varying values of space partition number of N and fractional order derivative. Specifically, the values of fractional derivatives are chosen as γ = 0.25, 0.5, 0.75, 0.9, respectively and presented in Table 1. Additionally, for a fixed value of N = 120, the time step size is varied from Δt = 0.002 to 0.0005 , and the corresponding numerical results are presented in Table 2. Upon examining Table 1, it is evident that increasing the number of spatial partitions has a significant effect on reducing both of the error norms. Notably, for N = 10, the infinity norm of the error decreases from 1.55586455×10⁻³ to 1.71989088×10⁻⁵, demonstrating a clear reduction trend. Similarly, Table 2 shows that a decrease in time step size generally leads to a reduction in the error norms. Subsequently, comparison tables (Tables 3 and 4) are presented to observe the contribution and effectiveness of the proposed method in solving the problem. In Table 3, for γ = 0.5 Δt = 0.00025, the results obtained using the proposed method with different numbers of spatial partitions (N = 40, 80 and 100) are compared with those from the cubic B-spline finite element method [10] and the trigonometric tension B-spline collocation method [11]. Table 4 presents a comparison of results for γ = 0.5, N = 120 and t_f = 1 using different time step sizes Δt = 0.002, 0.001 and 0.0005, again with references [10] and [11]. Examination of Table 3 reveals that the error norms obtained using the proposed higher-order trigonometric cubic B-spline collocation method are lower than those obtained with the cubic Bspline finite element method and are consistent with the results from [11]. Similarly, Table 4 indicates that for N = 120, the error norms obtained using a higher-order trigonometric cubic B-spline collocation method are lower than those obtained using both the cubic B-spline finite element method and the trigonometric tension B-spline collocation method. Moreover, the reduction in error norms becomes more significant as the time step size decreases.

Table 1

The error norms for Δt = 0.0005 and different values of N and γ of Problem 1.

N	Norm	γ = 0.25	γ = 0.5	γ = 0.75	γ = 0.9
10	L₂L_∞	1.17608481 × 10⁻³1.56691086 × 10⁻³	1.16808997 × 10⁻³1.55586455 × 10⁻³	1.16122058 × 10⁻³1.54633448 × 10⁻³	1.15977847 × 10⁻³1.54428531 × 10⁻³
20	L₂L_∞	3.24966594 × 10⁻⁴4.58600285 × 10⁻⁴	3.22510891 × 10⁻⁴4.55139112 × 10⁻⁴	3.21299945 × 10⁻⁴4.53437127 × 10⁻⁴	3.23144788 × 10⁻⁴4.56047172 × 10⁻⁴
40	L₂L_∞	7.78030654 × 10⁻⁵1.09833692 × 10⁻⁴	7.69389503 × 10⁻⁵1.08615283 × 10⁻⁴	7.73445392 × 10⁻⁵1.09189833 × 10⁻⁴	8.01209788 × 10⁻⁵1.13110804 × 10⁻⁴
80	L₂L_∞	1.25916352 × 10⁻⁵1.78288364 × 10⁻⁵	1.21469751 × 10⁻⁵1.71989088 × 10⁻⁵	1.29784520 × 10⁻⁵1.83758881 × 10⁻⁵	1.60002218 × 10⁻⁵2.26544495 × 10−5

Table 2

The error norms for N = 120 and different values of Δt and γ at T = 1 for Problem 1.

∆t	Norm	γ = 0.25	γ = 0.5	γ = 0.75	γ = 0.9
0.002	L₂L_∞	2.79799013 × 10⁻⁵3.96292257 × 10⁻⁵	2.89453357 × 10⁻⁵4.09955552 × 10⁻⁵	2.39205481 × 10⁻⁵3.38827254 × 10⁻⁵	1.09172130 × 10⁻⁵1.54759469 × 10⁻⁵
0.001	L₂L_∞	9.09662936 × 10⁻⁶1.28824323 × 10⁻⁵	9.70022354 × 10⁻⁶1.37368370 × 10⁻⁵	7.53455829 × 10⁻⁶1.06714152 × 10⁻⁵	1.19234428 × 10⁻⁶1.69358924 × 10⁻⁶
0.0005	L₂L_∞	3.57839140 × 10⁻⁷5.05989927 × 10⁻⁷	1.18088684 × 10⁻⁸2.09713397 × 10⁻⁸	9.03171600 × 10⁻⁷1.27789156 × 10⁻⁶	3.97091243 × 10⁻⁶5.62041476 × 10⁻⁶

Table 3

Comparison of error norms for γ = 0.5, Δt = 0.00025, t_f = 1 and different values of N for Problem 1.

N	Norm	Proposed method	Ref. [10]	Ref. [11]
40	L₂L_∞	8.18 × 10⁻⁵1.15 × 10⁻⁴	1.22 × 10⁻³1.73 × 10⁻³	1.60 × 10⁻⁵2.63 × 10⁻⁵
80	L₂L_∞	1.70 × 10⁻⁵2.40 × 10⁻⁵	1.78 × 10⁻⁴2.53 × 10⁻⁴	7.72 × 10⁻⁶1.34 × 10⁻⁵
100	L₂L_∞	9.14 × 10⁻⁶1.29 × 10⁻⁵	5.23 × 10⁻⁵7.65 × 10⁻⁵	7.24 × 10⁻⁶1.19 × 10⁻⁵

Table 4

Comparison of error norms for γ = 0.5, N = 120, t_f = 1 and different values of Δt for Problem 1.

∆t	Norm	Proposed method	Ref. [10]	Ref. [11]
0.002	L₂L_∞	2.89 × 10⁻⁵4.09 × 10⁻⁵	1.22 × 10⁻³1.72 × 10⁻³	5.55 × 10⁻⁵8.01 × 10⁻⁵
0.001	L₂L_∞	9.70 × 10⁻⁶1.37 × 10⁻⁵	5.32 × 10⁻⁴7.53 × 10⁻⁴	2.77 × 10⁻⁵3.92 × 10⁻⁵
0.0005	L₂L_∞	1.18 × 10⁻⁸2.09 × 10⁻⁸	1.89 × 10⁻⁴2.68 × 10⁻⁴	1.39 × 10⁻⁵2.05 × 10⁻⁵

To investigate the convergence rate of the method for Problem 1, the convergence ratios are presented in Table 5 and 6 for γ = 0.5, v = 1 and t_f = 1, using halved values of both (b–a)/N and Δt, Since both the spatial and temporal discretizations are refined simultaneously, the reported convergence rates represent the overall observed convergence of the fully discrete scheme. In the first table, the number of spatial partitions N is doubled while the Δt values are reduced by a factor of two. Thus, the convergence rates are 4.78, 4.51, 3.31 and 4.27, respectively. In the second table, the number of N = 10,20,40 is doubled again and similarly the Δt values are decreased by 1/4, 1/32 and 1/256. It is evident from both of the tables that the convergence rate is approximately 4.

Table 5

Maximum errors and convergence rates for γ = 0.5, v = 1, t_f = 1 for different values of N and Δt.

N	∆t	L_∞	RoC
6	1/5	4.51974778 × 10⁻³	−
12	1/20	1.64465660 × 10⁻⁴	4.78
24	1/80	7.17815902 × 10⁻⁶	4.51
48	1/320	7.21429488 × 10⁻⁷	3.31
96	1/1280	3.72153127 × 10⁻⁸	4.27

Table 6

Maximum errors and convergence rates for γ = 0.5, v = 1, t_f = 1 for different values of N and Δt.

N	∆t	L_∞	RoC
10	1/4	5.26044936 × 10⁻³	−
20	1/32	3.50923812 × 10⁻⁴	3.90
40	1/256	1.65112329 × 10⁻⁵	4.40

For example 5.1, we present five figures to illustrate the behavior of the numerical solutions. In Figure 1 (a) involves the graphs of numerical solutions and exact one for γ = 0.5, N = 50 and Δt = 0.001 at time t_f = 1. It can be seen from the Figure 1 (a), the numerical solutions exhibit the same behavior as the exact solution and two solutions are almost overlapping. In Figure 1 (b), the space partition number and step size are kept as N = 100 and Δt = 0.001 at t_f = 1, and the fractional order chosen as γ = 0.2, 0.5, 0.75, 0.97, respectively. It is clear that a very small change is observed between the solutions with increasing values of γ. For this reason, a zoomed subplot has been inserted in the figure so that the change can be seen. Figure 2 (a) and (b) depict the evolution of the numerical solutions of the time-fractional Burgers equation for N = 100, γ = 0.5, ∆t = 0.001 at time levels t = 0(0.1)1. It is seen that with the progression of time, the form of the wave also adapts to its structure. And, lastly in Figure 3, we present the graph of absolute error between the exact and numerical solutions across the spatial domain for N = 100, γ = 0.5, ∆t = 0.001 at t_f = 1. It is clear that from the Figure 3, the absolute error remains on the order of 10⁻⁶, which show high numerical accuracy. One can see a minimum error at x = 0.5 and maximum errors are occurred at the wave peaks.

5.2

Example 2

As the second example, we consider the time-fractional Burgers equation subject to the following initial and boundary conditions $\begin{matrix} u (0, t) = t^{2}, u (1, t) = \exp (1) t^{2}, t \geq 0 \\ u (x, 0) = 0, 0 \leq x \leq 1. \end{matrix}$ \matrix{ {u(0,t) = {t^2},u(1,t) = \exp (1){t^2},t \ge 0} \cr {u(x,0) = 0,0 \le x \le 1.} \cr }

The corresponding source term; $f (x, t) = \frac{2 t^{2 - γ}}{Γ (3 - γ)} \exp (x) + t^{4} \exp (2 x) - v t^{2} \exp (x)$ f(x,t) = {{2{t^{2 - \gamma }}} \over {\Gamma (3 - \gamma )}}\exp (x) + {t^4}\exp (2x) - v{t^2}\exp (x) and the exact solution of the problem 2 is $u (x, t) = t^{2} \exp (x) .$ u(x,t) = {t^2}\exp (x).

With the similar idea of first example, first of all, to show the effect of fractional order and partition number of space on the error norms, Δt is fixed to 0.0005, viscosity parameter v, and final time t_f are chosen equal to 1. Then, to show the effect of fractional order and space step size on the error norms L₂ and L_∞, the partition number of space domain is fixed to N = 120 for v = 1, t_f = 1. The obtained results of the error norms L₂ and L_∞ are presented in Tables 7 and 8, respectively. Regarding to Table 7, the error norms decrease with increasing values of N for each value of fractional order γ. Additionally, it is also observed that for N = 80 and N = 100, while the error values have an increasing behaviour at first, then it can be observed a decreasing behavior on the error norms as the derivative order γ increases. In other words, the derivative order has less negative effect on the error at values close to 0 and 1, whereas the errors increase around γ = 0.5. In Table 8, as it is expected, for each fixed fractional order γ, the error norms decrease as the number of spatial partitions increases. And, if we consider a fixed value of Δt, there can be seen a peak on error norms, however, for different values of fractional order γ, the decrease is again noticeable.

Table 7

The error norms L₂ and L_∞ for varying values of γ and N for ∆t = 0.0005 of Problem 2.

N	Norm	γ = 0.25	γ = 0.5	γ = 0.75	γ = 0.9
20	L₂L_∞	1.61504086 × 10⁻⁴3.08801388 × 10⁻⁴	1.44824109 × 10⁻⁴1.98520497 × 10⁻⁴	1.38338571 × 10⁻⁴1.89620989 × 10⁻⁴	1.15834339 × 10⁻⁴1.58724458 × 10⁻⁴
40	L₂L_∞	8.55740423 × 10⁻⁵1.17333257 × 10⁻⁴	8.82849416 × 10⁻⁵1.21056287 × 10⁻⁴	8.27624785 × 10⁻⁵1.13520309 × 10⁻⁴	6.08117897 × 10⁻⁵8.34162090 × 10⁻⁵
80	L₂L_∞	7.11925979 × 10⁻⁵9.76772696 × 10⁻⁵	7.41488030 × 10⁻⁵1.01737613 × 10⁻⁴	6.88664654 × 10⁻⁵9.45207070 × 10⁻⁵	4.70544770 × 10⁻⁵6.45980625 × 10⁻⁵
100	L₂L_∞	6.94667643 × 10⁻⁵9.53080411 × 10⁻⁵	7.24524065 × 10⁻⁵9.94124999 × 10⁻⁵	6.71988959 × 10⁻⁵9.22404235 × 10⁻⁵	4.54035631 × 10⁻⁵6.23420557 × 10⁻⁵

Table 8

The error norms L₂ and L_∞ for varying values of γ and Δt for N = 120 of Problem 2.

∆t	Norm	γ = 0.25	γ = 0.5	γ = 0.75	γ = 0.9
0.002	L₂L_∞	2.66079805 × 10⁻⁴3.65024346 × 10⁻⁴	2.75969290 × 10⁻⁴3.78638010 × 10⁻⁴	2.46238628 × 10⁻⁴3.37996968 × 10⁻⁴	1.54248522 × 10⁻⁴2.11803042 × 10⁻⁴
0.001	L₂L_∞	1.34644789 × 10⁻⁴1.84721757 × 10⁻⁴	1.40213575 × 10⁻⁴1.92384206 × 10⁻⁴	1.27686564 × 10⁻⁴1.75274357 × 10⁻⁴	8.28012104 × 10⁻⁵1.13702401 × 10⁻⁴
0.0005	L₂L_∞	6.85292690 × 10⁻⁵9.40182449 × 10⁻⁵	7.15309025 × 10⁻⁵9.81468420 × 10⁻⁵	6.62930520 × 10⁻⁵9.09997581 × 10⁻⁵	4.45067678 × 10⁻⁵6.11154536 × 10⁻⁵

Table 9 presents a comparison between the proposed method and the previously published method reported in Ref. [10]. for increasing values of the discretization number N and γ = 0.5, Δt = 0.00025, t_f = 1. From the Table 9, it can be seen that as N increases from 40 to 80, both method exhibit a decreasing behaviour in error norms which is expected. Even so, for both partition number N = 40 and 80, the proposed method yields lower errors than method used in [10] for the norms L₂ and L_∞. The second comparison table is Table 10 which provides a pointwise comparisons of numerical solutions at selected grids, obtained using the collocation method and [12], alongside the exact solution using the error L_∞ for γ = 0.5, ∆t = 0.00025, t_f = 1, ν = 1 and N = 40. It can be seen from the Table 10, across all data points, proposed collocation method based on trigonometric cubic B-splines yields lower errors than method used in [12]. We can conclude that proposed method offers more accurate approximations across the interval.

Table 9

Comparison of errors for γ = 0.5, ∆t = 0.00025, t_f = 1 and different values of N for Problem 2.

N	Method	L₂	L_∞
40	Proposed methodRef. [10]	5.36 × 10⁻⁵6.77 × 10⁻⁵	7.36 × 10⁻⁵2.09 × 10⁻⁴
80	Proposed methodRef. [10]	3.95 × 10⁻⁵4.57 × 10⁻⁵	5.42 × 10⁻⁵6.92 × 10⁻⁵

Table 10

Comparison of numerical solutions and errors for γ = 0.5, ∆t = 0.00025, t_f = 1, ν = 1 and N = 40 for Problem 2.

x	Proposed method	L_∞	Ref. [12]	L_∞	Exact solution
0.2	1.221366	3.66 × 10⁻⁵	1.221462	5.95 × 10⁻⁵	1.221402
0.4	1.491762	6.24 × 10⁻⁵	1.491934	1.09 × 10⁻⁴	1.491824
0.6	1.822045	7.36 × 10⁻⁵	1.822258	1.39 × 10⁻⁴	1.822118
0.8	2.225480	6.05 × 10⁻⁵	2.225666	1.025 × 10⁻⁴	2.225540

Tables 11 and 12 present the maximum error norms L_∞ and the rates of convergence (RoC) for various combinations of spatial discretization number N and time step size Δt for the values of the parameters γ = 0.9, ν = 1, and t_f = 1. In Table 11, while N increases from 30 to 120 and Δt is diminished from 1/4 to 1/1024, the error norms L_∞ decreases from 2.20295625×10⁻² to 1.11304212×10⁻⁴, and the RoC values are approximately 4.19 and 3.46. In Table 12, one can see a different set of partition number as N = 100, 200 and 400. For selected values of N and Δt, the error norms L_∞ decreases from 2.20820268×10⁻² to 1.08765049×10⁻⁴ and RoC values remain approximately 4.19 and 3.46. From the tables, it can be concluded that the results obtained demonstrate that the proposed method achieves an observed convergence rate of around 4, and the method is both accurate as well as economical, displaying high-order convergence properties for various spatial and temporal discretizations.

Table 11

The error norms L_∞ and convergence rates for varying values of N and Δt for γ = 0.9, ν = 1, t_f = 1 for Problem 2.

N	∆t	L_∞	RoC
30	1/4	2.20295625 × 10⁻²	−
60	1/64	1.21511503 × 10⁻³	4.18
120	1/1024	1.11304212 × 10⁻⁴	3.44

Table 12

The error norms L_∞ and convergence rates for varying values of N and Δt for γ = 0.9, ν = 1, t_f = 1 for Problem 2.

N	∆t	L_∞	RoC
100	1/4	2.20820268 × 10⁻²	−
200	1/64	1.20495216 × 10⁻³	4.19
400	1/1024	1.08765049 × 10⁻⁴	3.46

For example 5.2, to share the properties of the behaviour of numerical solutions, we present the Figures 4–6. In Figure 4 (a), the values are chosen as γ = 0.75, N = 50 and ∆t = 0.0025 at time t_f = 1 and numerical solutions and exact solutions are plotted on the same axis to show how well they are in agreement with each other. In Figure 4 (b), the partition number of x axis are chosen as N = 100, time step size Δt = 0.001 and the fractional orders are chosen as γ = 0.2,0.5,0.75,0.97, and numerical solutions of example 5.2 are depicted with a zoomed graphs to show the changing values of solutions regarding to fractional order. The behaviour of solutions at different time levels with 3D graphs for γ = 0.5, N = 100 and ∆t = 0.001 are presented in Figure 5 (a) and (b). Figure 5 (a) is a clear representation of the temporal evolution of the numerical solution profile at different time levels. Figure 5 (b), a 3D surface representation of the numerical solution, provides a helpful view of both spatial and temporal dynamics. And, Figure 6 shows the absolute error distribution between the exact and numerical solutions. It is clear from Figure 6 which shows a maximum absolute error of approximately 1.95×10⁻⁴ and a peak around x ≈ 0.6. Since, the error vanishes at the boundaries, it shows that boundary conditions are respected well.

5.3

Example 3

Finally, as the last example, we are going to consider time-fractional Burgers equation given in (1) via following initial and boundary conditions $\begin{array}{l} u (x, 0) = 0, 0 \leq x \leq 1 \\ u (0, t) = t^{2}, u (1, t) = - t^{2}, & t \geq 0. \end{array}$ \matrix{ {u(x,0) = 0,\;\;\;{\mkern 1mu} 0 \le x \le 1} \hfill & {} \hfill \cr {u(0,t) = {t^2},u(1,t) = - {t^2},{\mkern 1mu} } \hfill & {t \ge 0.} \hfill \cr }

For the problem 5.3, the corresponding source term and exact solution are given as $\begin{array}{l} f (x, t) = \frac{2 t^{2 - γ} \cos (π x)}{Γ (3 - γ)} - π t^{4} \cos (π x) \sin (π x) + v π^{2} t^{2} \cos (π x) \\ u (x, t) = t^{2} \cos (π x) . \end{array}$ \matrix{ {f(x,t) = {{2{t^{2 - \gamma }}\cos (\pi x)} \over {\Gamma (3 - \gamma )}} - \pi {t^4}\cos (\pi x)\sin (\pi x) + v{\pi ^2}{t^2}\cos (\pi x)} \hfill \cr {u(x,t) = {t^2}\cos (\pi x).} \hfill \cr }

In order to follow the same idea in previous problems, Table 13 and Table 14 are presented to follow the effect of changing fractional order derivative and space partition number and time step size on numerical solutions for fixed values of ν = 1, t_f = 1, ∆t = 0.00025 and N = 120, respectively. Table 13 demonstrates a distinct decrease in errors corresponding to an increase in the number of spatial partitions. Furthermore, when analyzing any row with a constant number of divisions, a comparable decline in error is noted as the value of γ rises. Table 14 illustrates an overall trend of diminishing error norms as the time step lowers and the value of γ approaches 1.

Table 13

The error norms L₂ and L_∞ for varying values of γ and N for Δt = 0.00025 of Problem 3.

N	Norm	γ = 0.1	γ = 0.2	γ = 0.4	γ = 0.6
10	L₂L_∞	3.01229397 × 10⁻⁴4.14249903 × 10⁻⁴	3.00321184 × 10⁻⁴4.13024518 × 10⁻⁴	2.98502812 × 10⁻⁴4.10575706 × 10⁻⁴	2.96792891 × 10⁻⁴4.08282315 × 10⁻⁴
20	L₂L_∞	7.28223128 × 10⁻⁵1.00695720 × 10⁻⁴	7.24982978 × 10⁻⁵1.00239858 × 10⁻⁴	7.19129216 × 10⁻⁵9.94134006 × 10⁻⁵	7.14959653 × 10⁻⁵9.88181649 × 10⁻⁵
40	L₂L_∞	1.52632425 × 10⁻⁵2.12934859 × 10⁻⁵	1.50860659 × 10⁻⁵2.10460160 × 10⁻⁵	1.48106479 × 10⁻⁵2.06611232 × 10⁻⁵	1.47187945 × 10⁻⁵2.05322390 × 10⁻⁵
80	L₂L_∞	8.18623973 × 10⁻⁷1.14201244 × 10⁻⁶	6.78279439 × 10⁻⁷9.46133375 × 10⁻⁷	4.80627584 × 10⁻⁷6.70126490 × 10⁻⁷	4.70380184 × 10⁻⁷6.55416123 × 10⁻⁷

Table 14

The error norms L₂ and L_∞ for varying values of γ and Δt for N = 120 of Problem 3.

∆t	Norm	γ = 0.25	γ = 0.5	γ = 0.75	γ = 0.9
0.002	L₂L_∞	3.14766794 × 10⁻⁵4.39151250 × 10⁻⁵	3.22741925 × 10⁻⁵4.50291619 × 10⁻⁵	2.77834353 × 10⁻⁵3.87661038 × 10⁻⁵	1.62793905 × 10⁻⁵2.27160596 × 10⁻⁵
0.001	L₂L_∞	1.46307138 × 10⁻⁵2.04121181 × 10⁻⁵	1.51160676 × 10⁻⁵2.10899969 × 10⁻⁵	1.31564294 × 10⁻⁵1.83571994 × 10⁻⁵	7.53330880 × 10⁻⁶1.05120674 × 10⁻⁵
0.0005	L₂L_∞	6.23879040 × 10⁻⁶8.70409358 × 10⁻⁶	6.51771718 × 10⁻⁶9.09364603 × 10⁻⁶	5.66868928 × 10⁻⁶7.90976865 × 10⁻⁶	2.93671313 × 10⁻⁶4.09826535 × 10⁻⁶

Tables 15 and 16 serve as comparative tables. In Table 15, γ = 0.5, ∆t = 0.00025, and t_f = 1 are kept constant, and the spatial partition numbers are taken as 40,80, and 100. The resultant error norms are compared to those in references [10] and [11]. For N = 40, the approach utilized in this study results in error norms comparable to [10] but less than those in [11]. Also, Table 15 clearly shows that as the number of spatial divisions rises, the suggested technique outperformed the other studies. Table 16 uses fixed values of γ = 0.5, N = 120, and t_f = 1 for comparing decreasing Δt values. The method produces reduced error norms for each specified Δt compared to those described in [10] and [11].

Table 15

Comparison of errors for γ = 0.5, ∆t = 0.00025, t_f = 1 and different values of N (Problem 3).

N Method	Method	L₂	L_∞
40	Proposed methodRef. [10]Ref. [11]	1.47 × 10⁻⁵1.22 × 10⁻³1.60 × 10⁻⁵	2.05 × 10⁻⁵1.73 × 10⁻³2.63 × 10⁻⁵
80	Proposed methodRef. [10]Ref. [11]	4.41 × 10⁻⁷1.78 × 10⁻⁴7.72 × 10⁻⁶	6.15 × 10⁻⁷2.53 × 10⁻⁴1.34 × 10⁻⁵
100	Proposed methodRef. [10]Ref. [11]	1.27 × 10⁻⁶5.23 × 10⁻⁵7.24 × 10⁻⁶	1.77 × 10⁻⁶7.65 × 10⁻⁵1.19 × 10⁻⁵

Table 16

Comparison of errors for γ = 0.5, N = 120, t_f = 1 and different values of Δt (Problem 3).

∆t	Method	L₂	L_∞
0.002	Proposed methodRef. [10]Ref. [11]	3.22 × 10⁻⁵1.22 × 10⁻³5.55 × 10⁻⁵	4.50 × 10⁻⁵1.72 × 10⁻³8.01 × 10⁻⁵
0.001	Proposed methodRef. [10]Ref. [11]	1.51 × 10⁻⁵5.32 × 10⁻⁴2.77 × 10⁻⁵	2.10 × 10⁻⁵7.53 × 10⁻⁴3.92 × 10⁻⁵
0.0005	Proposed methodRef. [10]Ref. [11]	6.51 × 10⁻⁶1.89 × 10⁻⁴1.39 × 10⁻⁵	9.09 × 10⁻⁶2.68 × 10⁻⁴2.05 × 10⁻⁵

Tables 17–19 demonstrate how changing spatial and temporal step sizes influences the convergence order for γ = 0.5, v = 1, and, t_f = 1. The results demonstrate that the maximum error decreases as the number of spatial divisions increases. Furthermore, achieving a convergence rate close to 4 with finer spatial and temporal discretizations indicates that the use of higher-order basis functions in the construction of the algorithm allows for a better alignment with the solution behavior of the problem.

Table 17

Maximum errors and convergence rates for γ = 0.5, ν = 1, t_f = 1 and different values of Δt and N (Problem 3).

N	∆t	L_∞	RoC
8	1/5	6.43241077 × 10⁻³	−
16	1/40	4.49763147 × 10⁻⁴	3.83
32	1/320	3.33582433 × 10⁻⁵	3.75

Table 18

Maximum errors and convergence rates for γ = 0.5, ν = 1, t_f = 1 and different values of Δt and N (Problem 3).

N	∆t	L_∞	RoC
5	1/4	7.90193068 × 10⁻³	−
10	1/32	3.54754064 × 10⁻⁴	4.47
20	1/256	1.22189859 × 10⁻⁵	4.85

Table 19

Maximum errors and convergence rates for γ = 0.5, ν = 1, t_f = 1 and different values of Δt and N (Problem 3).

N	∆t	L_∞	RoC
12	1/5	1.64465660 × 10⁻⁴	−
24	1/40	7.17815902 × 10⁻⁶	3.64
48	1/320	7.21429488 × 10⁻⁷	3.25
96	1/2560	3.72153127 × 10⁻⁸	3.55

For the problem 5.3, in order to reiterate the visual aims of problem 1 and problem 2, Figure 7 (a) is given to show the correspondence between the numerical solution and the exact solution for problem 3 for the values γ = 0.5, ν = 1, N = 100 and ∆t = 0.001. Then, in Figure 7 (b), the change of the fractional derivative on the curve is shown with a zoomed image by choosing N = 100, ∆t = 0.001, t_f = 1 for the change of the fractional order γ = 0.2, 0.5, 0.75, and, 0.97. The temporal and temporal-spatial evolutions of the nature of the solution are shown in Figure 8 (a) and (b) for γ = 0.5, N = 100 and ∆t = 0.001 in 3D plots. At the end, absolute error distribution between the exact and numerical solutions are presented in Figure 9. The absolute error shows a similar behaviour to the absolute error in the first problem, and the vanishing of the error at the boundaries again indicates that the boundary conditions are well provided. From a computational perspective, the proposed method retains the local support property of cubic trigonometric B-splines, leading to sparse and banded system matrices. Therefore, the computational cost per time step can be of the same order as that of existing spline collocation methods. Moreover, the higher-order spatial accuracy enables higher accuracy on relatively coarse grids, resulting in improved efficiency in terms of accuracy versus computational effort.

Conclusions

In the current article, we have designed and implemented a numerical approach for solving the time-fractional Burgers equation using a higher-order trigonometric cubic B-spline collocation method. The numerical experiments showed that the suggested technique has a convergence rate around four, suggesting a high degree of accuracy. The method’s efficacy was further tested by doing error analysis and comparing the results to those of existing methodologies. The smoothness and local support of the trigonometric cubic B-spline basis functions contributed to the scheme’s stability and accuracy, making it well-suited for dealing with fractional-order differential equations, which are nonlinear and nonlocal. These properties emphasize the potential of the method for larger applications in complicated physical models driven by fractional dynamics. In future work, we hope to investigate the applicability of our technique to various types of fractional partial differential equations. Furthermore, we have intended to study how other types of basis functions such as hyperbolic, polynomial, or wavelet-based splines affect the method’s accuracy, computing efficiency, and resilience. These investigations will help to further the development of efficient numerical tools for fractional calculus and its applications.

An efficient higher-order trigonometric cubic B-spline collocation method for timefractional Burgers equations

Full Article

Paradigm

My account