An efficient higher-order trigonometric cubic B-spline collocation method for time-fractional 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;

D_{t}^{γ} u + u u_{x} - v u_{x x} = f (x, t),

\begin{equation} D_{t}^{\gamma }u+uu_{x}-\nu u_{xx}=f(x,t), \label{intro1} \end{equation}

and,

\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}

\begin{equation} \begin{array}{cc} u(a,t)=h_{1}(t), & u\left( b,t\right) =h_{2}\left( t\right) , \\ u_{x}(a,t)=h_{3}(t), & u_{x}\left( b,t\right) =h_{4}\left( t\right) , \\ u\left( x,0\right) =g\left( x\right) , & \end{array} \label{intro2} \end{equation}%

where $x \in [a, b], t \in [0, T]$ x\in \left[ a,b\right] ,t\in \left[ 0,T\right] 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 the 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 second- or 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 presented 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 time-fractional 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 Riemann-Liouville 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^{\gamma}_{a} f(t) = \frac{1}{\Gamma(n - \gamma)} \frac{d^n}{dt^n} \int_{a}^{t} \frac{f(\tau)}{(t - \tau)^{\gamma - n + 1}} \, d\tau, \]

where n = [γ] 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 τ,$ \[ {}^{C}D^{\gamma}_{a} f(t) = \frac{1}{\Gamma(n - \gamma)} \int_{a}^{t} \frac{f^{(n)}(\tau)}{(t - \tau)^{\gamma - n + 1}} \, d\tau, \]

where n = [γ], f⁽ⁿ⁾ 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) .

\begin{equation} T_{i}^{k}\left( x\right) =\frac{\sin (\frac{x-x_{i}}{2})}{\sin \left( \frac{% x_{i+k}-x_{i}}{2}\right) }T_{i}^{k-1}\left( x\right) +\frac{\sin (\frac{% x_{i+k+1}-x}{2})}{\sin \left( \frac{x_{i+k+1}-x_{i+1}}{2}\right) }% T_{i+1}^{k-1}\left( x\right). \label{s1} \end{equation}

Starting with k = 3 in the Eq.(3) trigonometric cubic B-spline functions can be obtained as follows

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}\left( x\right) =\frac{p(x_{i})}{\sin(\frac{3h}{2})}T_{i}^{2}\left( x\right) +\frac{p(x_{i+4})}{\sin(\frac{3h}{2})}T_{i+1}^{2}\left( x\right), \label{s2}

where $p (x_{i}) = \sin (\frac{x - x_{i}}{2})$ p(x_{i})=\sin(\frac{x-x_{i}}{2}). 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 (x_{i}) (p (x_{i}) p (x_{i + 2}) + 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 (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 \leq x_{i + 4}, \\ 0, & other cases, \end{array}

T_{i}^{3}\left( x\right) =\frac{1}{\theta }\left\{ \begin{array}{ll}{p}^{3}\left( x_{i}\right) , & 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}\left( x_{i+4}\right) & \ \ x_{i+3}\leq x<x_{i+4}, \\ 0 & \mathrm{other\ cases,\ }% \end{array}% \right.

where $θ = \sin (\frac{h}{2}) \sin (h) \sin (\frac{3 h}{2})$ \theta=\sin(\frac{h}{2})\sin(h)\sin(\frac{3h}{2}). Generally, an approximate solution of given problem using cubic trigonometric B-spline functions can be defined as

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_{i=-1}^{N+1}T_{i}^{3}(x)Δ_{i}(t)

Here, the $u (x, t)$ u(x,t) symbolize the exact solution of the problem and the coefficients $δ_{i}$ Δ_{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}$ x_{j} are

\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}

\begin{aligned}S(x,t)&=a_{1}Δ_{j-1}+a_{2}Δ_{j}+a_{3}Δ_{j+1} \\ S^{\prime}(x,t)&=b_{1}Δ_{j-1}+b_{2}Δ_{j+1} \\ S^{\prime\prime}(x,t)&=c_{1}Δ_{j-1}+c_{2}Δ_{j}+c_{3}Δ_{j+1}.\end{aligned}

Here;

\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}

a_{1}={\sin}^{2}(\frac{h}{2})\csc(h)\csc(\frac{3h}{2}), a_{2}=\frac{2}{1+2\cos(h)},\ \ a_{1}=a_{3} \\ b_{1}=-\frac{3\csc(\frac{3h}{2})}{4},\ \ b_{2}=\frac{3\csc(\frac{3h}{2})}{4}\\ c_{1}=\frac{3(3{\cos}^{2}(\frac{h}{2})-1)}{4(\sin(h)\sin(\frac{3h}{2}))},c_{2}=-\frac{3{\cot}^{2}(\frac{h}{2})}{(2+4\cos(h))},c_{3}=c_{1}.%}

2.2

Higher-order approach

As defined in the previous subsection, the approximate solution is given. Now, in this subsection, a higher-order nodal approximation of the solution and its derivatives is presented. Suppose that (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

\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{aligned} S(x_{j},t)&=u(x_{j},t), & 0\le j\le N \\ S^{\prime\prime}(x_{j},t)&=u^{\prime\prime}(x_{j},t)-\frac{1}{12}h^{2}u^{(4)}(x_{j},t), & j=0,N \end{aligned}

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]

\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{aligned}S^{\prime}(x_{j},t)&=u^{\prime}(x_{j},t)+O(h^{4}), & 0\le j\le N \\ S^{\prime\prime}(x_{j},t)&=u^{\prime\prime}(x_{j},t)-\frac{1}{12}h^{2}u^{(4)}(x_{j},t)+O(h^{4}), & 0\le j\le N.\end{aligned}

Using finite differences and Taylor expansion, the approximate solution and its derivatives at the nodes can be defined as following way

\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}

S(x_{j})=u(x_{j}),\ \ \ & 0\leq j\leq N \\ S^{\prime \prime }(x_{j})=u^{\prime \prime }(x_{j}),\ \ & \ j=0,N%

and

u^{(4)} (x_{j}) = \frac{S^{″} (x_{j - 1}) - 2 S^{″} (x_{j}) + S^{″} (x_{j + 1})}{h^{2}} + O (h^{2}) .

u^{(4)}(x_{j})=\frac{S^{\prime\prime}(x_{j-1})-2S^{\prime\prime}(x_{j})+S^{\prime\prime}(x_{j+1})}{h^{2}}+O(h^{2}).

Here, one can obtain

\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}

\begin{aligned}u^{(4)}(x_{1})&=\frac{S^{\prime\prime}(x_{0})-2S^{\prime\prime}(x_{1})+S^{\prime\prime}(x_{2})}{h^{2}}+O(h^{2}), & j=1 \\ u^{(4)}(x_{2})&=\frac{S^{\prime\prime}(x_{1})-2S^{\prime\prime}(x_{2})+S^{\prime\prime}(x_{3})}{h^{2}}+O(h^{2}), & j=2\end{aligned}

and

u^{(4)} (x_{0}) = 2 u^{(4)} (x_{1}) - u^{(4)} (x_{2}) .

u^{(4)}(x_{0})=2u^{(4)}(x_{1})-u^{(4)}(x_{2}).

When we use the equalities in (12) in the equation given in (13), it can easily be obtained

\begin{aligned} 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{aligned}

u^{(4)}(x_{0}) = \frac{2S^{\prime\prime}(x_{0})-4S^{\prime\prime}(x_{1})+2S^{\prime\prime}(x_{2})-S^{\prime\prime}(x_{1})+2S^{\prime\prime}(x_{2})-S^{\prime\prime}(x_{3})}{h^{2}} = \frac{2S^{\prime\prime}(x_{0})-5S^{\prime\prime}(x_{1})+4S^{\prime\prime}(x_{2})-S^{\prime\prime}(x_{3})}{h^{2}} + O(h^{2}).

Similarly, using (11), it yields

\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}

\begin{aligned} u^{(4)}(x_{N-1})&=\frac{S^{\prime\prime}(x_{N-2})-2S^{\prime\prime}(x_{N-1})+S^{\prime\prime}(x_{N})}{h^2}, & j=N-1 \\ u^{(4)}(x_{N-2})&=\frac{S^{\prime\prime}(x_{N-3})-2S^{\prime\prime}(x_{N-2})+S^{\prime\prime}(x_{N-1})}{h^{2}}, & j=N-2\end{aligned}

and using this equalities for

u^{(4)} (x_{N}) = 2 u^{(4)} (x_{N - 1}) - u^{(4)} (x_{N - 2})

u^{(4)}(x_{N})=2u^{(4)}(x_{N-1})-u^{(4)}(x_{N-2})

We get

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)}(x_{N})=\frac{2S^{\prime\prime}(x_{N})-5S^{\prime\prime}(x_{N-1})+4S^{\prime\prime}(x_{N-2})-S^{\prime\prime}(x_{N-3})}{h^{2}}+O(h^{2}).

Lastly, we can obtain as following approximations for fourth order derivative;

\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}

\begin{aligned}u^{(4)}(x_{j},t)&=\frac{2S^{\prime\prime}(x_{0},t)-5S^{\prime\prime}(x_{1},t)+4S^{\prime\prime}(x_{2},t)-S^{\prime\prime}(x_{3},t)}{h^{2}}+O(h^{2}) , & j=0 \\ u^{(4)}(x_{j},t)&=\frac{S^{\prime\prime}(x_{j-1},t)-2S^{\prime\prime}(x_{j},t)+S^{\prime\prime}(x_{j+1},t)}{h^{2}}+O(h^{2}) , & 1\le j\le N-1 \\ u^{(4)}(x_{j},t)&=\frac{2S^{\prime\prime}(x_{N},t)-5S^{\prime\prime}(x_{N-1},t)+4S^{\prime\prime}(x_{N-2},t)-S^{\prime\prime}(x_{N-3},t)}{h^{2}}+O(h^{2}) , & j=N.\end{aligned}

Using equations given by (9), we get following approximations for seconds order derivative;

\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) - S^{″} (x_{N - 3}, t)}{12} + O (h^{4}), & j = N . \end{array}

\begin{aligned}u^{\prime\prime}(x_{0},t)&=\frac{14S^{\prime\prime}(x_{0},t)-5S^{\prime\prime}(x_{1},t)+4S^{\prime\prime}(x_{2},t)-S^{\prime\prime}(x_{3},t)}{12}+O(h^{4}), & j=0 \\ u^{\prime\prime}(x_{j},t)&=\frac{S^{\prime\prime}(x_{j-1},t)+10S^{\prime\prime}(x_{j},t)+S^{\prime\prime}(x_{j+1},t)}{12}+O(h^{4}), & 1\le j\le N-1 \\ u^{\prime\prime}(x_{N},t)&=\frac{14S^{\prime\prime}(x_{N},t)-5S^{\prime\prime}(x_{N-1},t)+4S^{\prime\prime}(x_{N-2},t)-S^{\prime\prime}(x_{N-3},t)}{12}+O(h^{4}) , & j = N.\end{aligned}

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 the aforementioned equation can be decomposed into 5 key steps. The first one is discretization of the equation employing the 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

\frac{δ t^{- γ}}{γ (2 - γ)} \sum_{k = 0}^{n} b_{k} [u_{j}^{n + 1 - k} - u_{j}^{n - k}] + \frac{1}{2} {(u u_{x})}_{j}^{n + 1} + \frac{1}{2} {(u u_{x})}_{j}^{n} - \frac{v}{2} {(u_{x x})}_{j}^{n + 1} - \frac{v}{2} {(u_{x x})}_{j}^{n} = \frac{f (x_{j}, t_{n + 1}) + f (x_{j}, t_{n})}{2}

\frac{Δ t^{-\gamma}}{\Gamma(2-\gamma)}\sum_{k=0}^{n}b_{k}[u_j^{n+1-k}-u_j^{n-k}]+\frac{1}{2}(uu_x)_j^{n+1}+\frac{1}{2}(uu_x)_j^{n}-\frac{v}{2}(u_{xx})_j^{n+1}-\frac{v}{2}(u_{xx})_j^{n}=\frac{f(x_j,t_{n+1})+f(x_j,t_n)}{2}

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)^{n + 1} is linearized using the Rubin-Graves type technique according to the following formula [29,30]

{(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} .

(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}.

For p = 1 the Rubin-Graves linearization yields

{(u u_{x})}^{n + 1} \approx u^{n} u_{x}^{n + 1} + u_{x}^{n} u^{n + 1} - u^{n} u_{x}^{n}

(uu_{x})^{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

\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}

u^{n+1}(1+\frac{S\psi_{1}}{2})+u_{x}^{n+1}(\frac{S\psi_{2}}{2})-u_{xx}^{n+1}(\frac{Sv}{2})=u^{n}+u_{xx}^{n}(\frac{Sv}{2}) +\frac{S}{2}(f(x,t_{n+1})+f(x,t_{n}))-\Sigma_{k=1}^{n}b_{k}[u^{n+1-k}-u^{n-k}].

Here, $S = δ t^{- γ} / γ (2 - γ)$ S=Δ t^{-\gamma}/\Gamma(2-\gamma), $ψ_{1} = u_{x}^{n}$ \Psi_{1}=u_{x}^{n} and $ψ_{2} = u^{n}$ \psi_{2}=u^{n}. Accordingly, the coefficients multiplying u^{n + 1} 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 the previous section, and using approximate solutions, the third step of the numerical implementation is resulted as the following algebraic equation;

\begin{array}{l} j = 0 : & {\begin{cases} δ_{- 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{cases} \end{array}

j=0: \begin{cases} \delta_{-1}^{n+1} \left[a_{1}(1+\frac{S}{2}\psi_{1}) + \frac{Sb_{1}}{2}\psi_{2} - \frac{Sv}{24}(14c_{1})\right] + \delta_{0}^{n+1} \left[a_{2}(1+\frac{S}{2}\psi_{1}) - \frac{Sv}{24}(14c_{2}-5c_{1})\right] \\ + \delta_{1}^{n+1} \left[a_{3}(1+\frac{S}{2}\psi_{1}) + \frac{Sb_{2}}{2}\psi_{2} - \frac{Sv}{24}(14c_{3}-5c_{2}+4c_{1})\right] + \delta_{2}^{n+1} \left[-\frac{Sv}{24}(-5c_{3}+4c_{2}-c_{1})\right] \\ + \delta_{3}^{n+1} \left[-\frac{Sv}{24}(4c_{3}-c_{2})\right] + \delta_{4}^{n+1} \left[-\frac{Sv}{24}(-c_{3})\right] \\ \qquad = \delta_{-1}^{n} \left[a_{1} + \frac{Sv}{24}(14c_{1})\right] + \delta_{0}^{n} \left[a_{2} + \frac{Sv}{24}(14c_{2}-5c_{1})\right] + \delta_{1}^{n} \left[a_{3} + \frac{Sv}{24}(14c_{3}-5c_{2}+4c_{1})\right] \\ \qquad + \delta_{2}^{n} \left[\frac{Sv}{24}(-5c_{3}+4c_{2}-c_{1})\right] + \delta_{3}^{n} \left[\frac{Sv}{24}(4c_{3}-c_{2})\right] + \delta_{4}^{n} \left[\frac{Sv}{24}(-c_{3})\right] + \frac{S}{2}(f(x_{0},t_{n+1})+f(x_{0},t_{n})) \\ \qquad - \sum_{k=1}^{n} b_{k} \left[ (a_{1}\delta_{-1}^{n+1-k} + a_{2}\delta_{0}^{n+1-k} + a_{3}\delta_{1}^{n+1-k}) - (a_{1}\delta_{-1}^{n-k} + a_{2}\delta_{0}^{n-k} + a_{3}\delta_{1}^{n-k}) \right], \end{cases}

and similarly,

\begin{array}{l} j = 1, 2, . . ., N - 1 : & {\begin{cases} δ_{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{cases} \end{array}

j=1,2,\dots,N-1: \begin{cases} \delta_{j-2}^{n+1} \left[-\frac{Sv}{24}c_{1}\right] + \delta_{j-1}^{n+1} \left[a_{1}(1+\frac{S}{2}\psi_{1}) + \frac{Sb_{1}}{2}\psi_{2} - \frac{Sv}{24}(c_{2}+10c_{1})\right] \\ + \delta_{j}^{n+1} \left[a_{2}(1+\frac{S}{2}\psi_{1}) - \frac{Sv}{24}(c_{3}+10c_{2}+c_{1})\right] \\ + \delta_{j+1}^{n+1} \left[a_{3}(1+\frac{S}{2}\psi_{1}) + \frac{Sb_{2}}{2}\psi_{2} - \frac{Sv}{24}(10c_{3}+c_{2})\right] + \delta_{j+2}^{n+1} \left[-\frac{Sv}{24}c_{3}\right] \\ \qquad = \delta_{j-2}^{n} \left[\frac{Sv}{24}c_{1}\right] + \delta_{j-1}^{n} \left[a_{1} + \frac{Sv}{24}(c_{2}+10c_{1})\right] + \delta_{j}^{n} \left[a_{2} + \frac{Sv}{24}(c_{3}+10c_{2}+c_{1})\right] \\ \qquad + \delta_{j+1}^{n} \left[a_{3} + \frac{Sv}{24}(10c_{3}+c_{2})\right] + \delta_{j+2}^{n} \left[\frac{Sv}{24}(c_{3})\right] + \frac{S}{2}(f(x_{j},t_{n+1})+f(x_{j},t_{n})) \\ \qquad - \sum_{k=1}^{n} b_{k} \left[ (a_{1}\delta_{m-1}^{n+1-k} + a_{2}\delta_{m}^{n+1-k} + a_{3}\delta_{m+1}^{n+1-k}) - (a_{1}\delta_{m-1}^{n-k} + a_{2}\delta_{m}^{n-k} + a_{3}\delta_{m+1}^{n-k}) \right], \end{cases}

$j = N :$ j=N:

\begin{array}{l} j = N : & {\begin{cases} δ_{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{cases} \end{array}

j=N: \begin{cases} \delta_{N-4}^{n+1} \left[\frac{Sv}{24}c_{1}\right] + \delta_{N-3}^{n+1} \left[-\frac{Sv}{24}(-c_{2}+4c_{1})\right] + \delta_{N-2}^{n+1} \left[-\frac{Sv}{24}(-c_{3}+4c_{2}-5c_{1})\right] \\ + \delta_{N-1}^{n+1} \left[a_{1}(1+\frac{S}{2}\psi_{1}) + \frac{Sb_{1}}{2}\psi_{2} - \frac{Sv}{24}(4c_{3}-5c_{2}+14c_{1})\right] + \delta_{N}^{n+1} \left[a_{2}(1+\frac{S}{2}\psi_{1}) - \frac{Sv}{24}(-5c_{3}+14c_{2})\right] \\ + \delta_{N+1}^{n+1} \left[a_{3}(1+\frac{S}{2}\psi_{1}) + \frac{Sb_{2}}{2}\psi_{2} - \frac{Sv}{24}(14c_{3})\right] \\ \qquad = \delta_{N-4}^{n} \left[-\frac{Sv}{24}c_{1}\right] + \delta_{N-3}^{n} \left[\frac{Sv}{24}(-c_{2}+4c_{1})\right] + \delta_{N-2}^{n} \left[\frac{Sv}{24}(-c_{3}+4c_{2}-5c_{1})\right] \\ \qquad + \delta_{N-1}^{n} \left[a_{1} + \frac{Sv}{24}(4c_{3}-5c_{2}+14c_{1})\right] + \delta_{N}^{n} \left[a_{2} + \frac{Sv}{24}(-5c_{3}+14c_{2})\right] \\ \qquad + \delta_{N+1}^{n} \left[a_{3} + \frac{Sv}{24}(14c_{3})\right] + \frac{S}{2}(f(x_{N},t_{n+1})+f(x_{N},t_{n})) \\ \qquad - \sum_{k=1}^{n} b_{k} \left[ (a_{1}\delta_{N-1}^{n+1-k} + a_{2}\delta_{N}^{n+1-k} + a_{3}\delta_{N+1}^{n+1-k}) - (a_{1}\delta_{N-1}^{n-k} + a_{2}\delta_{N}^{n-k} + a_{3}\delta_{N+1}^{n-k}) \right]. \end{cases}

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 involves the reducing to system into N + 1 equations and N + 1 unknown by utilizing the boundary conditions. Eliminating the parameters δ_{- 1} and δ_{N + 1} yields an algebraic equation system ready to be solved by iteratively. In the sake of seeing clearly, the discretization boundary conditions and removing process can be expressed by the following equations;

\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}

u(a,t)=a_{1}{Δ }_{-1}+a_{2}{Δ }_{0}+a_{3}{Δ }_{1}=h_{1}(t)\ \\ \\ \hspace{0.5cm}{Δ }_{-1}=\frac{h_{1}(t)\ }{a_{1}}-\frac{a_{2}}{a_{1}}{% Δ }_{0}-\frac{a_{3}}{a_{1}}{Δ }_{1} \\ \\ u\left( b,t\right) =a_{1}{Δ }_{N-1}^{n+1-k}+a_{2}{Δ }% _{N}^{n+1-k}+a_{3}{Δ }_{N+1}^{n+1-k}=h_{2}\left( t\right) \\ \\ \hspace{0.5cm}{Δ }_{N+1}=\frac{h_{2}\left( t\right) \ }{a_{3}}-\frac{% a_{2}}{a_{3}}{Δ }_{N}-\frac{a_{1}}{a_{3}}{Δ }_{N-1}.%

The equation system obtained above can be written in matrix form as

A δ^{n + 1} = B δ^{n} .

AΔ^{n+1}=BΔ^{n}.

As the last step, the initial vector $[δ_{0}^{0}, δ_{1}^{0}, δ_{2}^{0}, \dots, δ_{N}^{0}]$ [Δ_{0}^{0},Δ_{1}^{0},Δ_{2}^{0},\cdot\cdot\cdot,Δ_{N}^{0}] 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

\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}) \\ \dots \\ j = N a_{1} δ_{N - 1} + a_{2} δ_{N} + a_{3} δ_{N + 1} = g (x_{N}) \end{array}

\begin{matrix}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{matrix}

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 δ_{- 1} and δ_{N + 1} like the previous step using the Neumann boundary conditions u_x(a,t)=h₃(t) and u_x(b,t)=h₄(t). The implementation of the method will be completed via obtaining the initial vector as follows

\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} \\ 1 \leq j \leq N - 1; a_{1} δ_{j - 1} + a_{2} δ_{j} + a_{3} δ_{j + 1} = g (x_{j}) \\ 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}

\begin{matrix}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}\\ 1\le j\le N-1;& a_{1}Δ_{j-1}+a_{2}Δ_{j}+a_{3}Δ_{j+1}=g(x_{j})\\ 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{matrix}

Consequently, to obtain numerical solutions, the system given in (24) - (26) is going to be solved 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}$ Δ_{j}^{n+1} are going to be obtained by the parameters $δ_{j}^{n}$ Δ_{j}^{n} at desired time level n. And, using the known values of $δ_{j}^{n + 1}$ Δ_{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 via the application of the suggested method for time-fractional Burgers equation using the von-Neumann criterion and with the help of mathematical induction. According to the von-Neumann criterion, if the amplification factor $| λ | \leq 1$ |\lambda|\le1, 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 the suggested method related with the linearized form of the Eq. (1) by

\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}

\begin{array}{l} Δ _{j-2}^{n+1}\left( \frac{-S\nu c_{1}}{24}\right) +Δ _{j-1}^{n+1}\left( a_{1}+\frac{S}{2}\left( \tau b_{1}-\frac{\nu}{12}\left( 10c_{1}+c_{2}\right) \right) \right) +Δ _{j}^{n+1}\left( a_{2}-\frac{S\nu}{12}\left( c_{1}+5c_{2}\right) \right) \\ \\ +Δ _{j+1}^{n+1}\left( a_{1}+\frac{S}{2}\left( \tau b_{2}-\frac{\nu}{12}% \left( 10c_{1}+c_{2}\right) \right) \right) +Δ _{j+2}^{n+1}\left( \frac{% -S\nu c_{1}}{24}\right) \\ \\ \hspace{1cm}=Δ _{j-2}^{n}\left( \frac{S\nu c_{1}}{24}\right) +Δ _{j-1}^{n}\left( a_{1}-\frac{S}{2}\left( \tau b_{1}-\frac{\nu}{12}\left( 10c_{1}+c_{2}\right) \right) \right) +Δ _{j}^{n}\left( a_{2}+\frac{S\nu}{% 12}\left( c_{1}+5c_{2}\right) \right) \\ \\ \hspace{1.15cm}+Δ _{j+1}^{n}\left( a_{1}-\frac{S}{2}\left( \tau b_{2}-% \frac{\nu}{12}\left( 10c_{1}+c_{2}\right) \right) \right) +Δ _{j+2}^{n}\left( \frac{S\nu c_{1}}{24}\right) \\ \\ \hspace{1.15cm}-\sum\limits_{k=0}^{n}b_{k}^{\alpha }\left( \left( a_{1}Δ _{j-1}^{n+1-k}+a_{1}Δ _{j}^{n+1-k}+a_{1}Δ _{j+1}^{n+1-k}\right) -\left( a_{1}Δ _{j-1}^{n-k}+a_{1}Δ _{j}^{n-k}+a_{1}Δ _{j+1}^{n-k}\right) \right)% \end{array} \label{stab1}

where $τ = u_{j}^{n}$ \tau=u_{j}^{n} Assume that the amplification factor of the Fourier mode is defined as

δ_{j}^{n} = λ^{n} e^{i j h β}

Δ _{j}^{n}=\lambda ^{n}e^{ijh\beta} \label{stab2}

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

\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}

\lambda^{n+1}\{(\frac{-Svc_{1}}{24})(e^{-2ih\beta}+e^{2ih\beta})+(a_{1}-\frac{Sv}{24}(10c_{1}+c_{2}))(e^{-ih\beta}+e^{ih\beta})+a_{2}-\frac{Sv}{12}(c_{1}+5c_{2})+\frac{s}{2}\tau(b_{1}+b_{2})\} =\lambda^{n}\{(\frac{Svc_{1}}{24})(e^{-2ih\beta}+e^{2ih\beta})+(a_{1}+\frac{Sv}{24}(10c_{1}+c_{2}))(e^{-ih\beta}+e^{ih\beta})+a_{2}+\frac{Sv}{12}(c_{1}+5c_{2})-\frac{S}{2}\tau(b_{1}+b_{2})\} -\lambda^{n}\sum_{k=0}^{n}b_{k}^{a}(\lambda^{1-k}(a_{2}+a_{1}(e^{-ih\beta}+e^{ih\beta}))-\lambda^{-k}(a_{2}+a_{1}(e^{-ih\beta}+e^{ih\beta}))).

We know that $e^{- i h β} + e^{i h β} = 2 \cos (h β)$ e^{-ih\beta} + e^{ih\beta} = 2\cos(h\beta), $e^{i h β} - e^{- i h β} = 2 i \sin (h β)$ e^{ih\beta} - e^{-ih\beta} = 2i\sin(h\beta), $b_{2} = - b_{1}$ b_{2}=-b_{1} grouping the resembling terms, we get

\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}

\lambda^{n+1}\{2(\frac{-Svc_{1}}{24})cos(2h\beta)+2(a_{1}-\frac{Sv}{24}(10c_{1}+c_{2}))cos(h\beta)+a_{2}-\frac{Sv}{12}(c_{1}+5c_{2})-iS\tau b_{1}sin(h\beta)\} =\lambda^{n}\{2(\frac{Svc_{1}}{24})cos(2h\beta)+2(a_{1}+\frac{Sv}{24}(10c_{1}+c_{2}))cos(h\beta)+a_{2}+\frac{Sv}{12}(c_{1}+5c_{2})+iS\tau b_{1}sin(h\beta)\} -\lambda^{n}\sum_{k=0}^{n}b_{k}^{\alpha}(\lambda^{1-k}(2a_{1}cos(h\beta)+a_{2})-\lambda^{-k}(2a_{1}cos(h\beta)+a_{2})).

After some arrangements, we have

\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}

\lambda^{n+1}\{2a_{1}cos(h\beta)+a_{2}-\frac{Sv}{24}(2c_{1}cos(2h\beta)+2(10c_{1}+c_{2})cos(h\beta)+(2c_{1}+10c_{2}))-iS\tau b_{1}sin(h\beta)\} =\lambda^{n}\{2a_{1}cos(h\beta)+a_{2}+\frac{Sv}{24}(2c_{1}cos(2h\beta)+2(10c_{1}+c_{2})cos(h\beta)+(2c_{1}+10c_{2}))+iS\tau b_{1}sin(h\beta)\} -\lambda^{n}\sum_{k=0}^{n}b_{k}^{\alpha}(\lambda^{1-k}(2a_{1}cos(h\beta)+a_{2})-\lambda^{-k}(2a_{1}cos(h\beta)+a_{2})).

Without the loss of generality, we can assume that $β = 0$ \beta=0 as in [32, 33], it yields

\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}

\lambda^{n+1}\{2a_{1}+a_{2}-\frac{Sv}{24}(2c_{1}+2(10c_{1}+c_{2})+(2c_{1}+10c_{2}))\} =\lambda^{n}\{2a_{1}+a_{2}+\frac{Sv}{24}(2c_{1}+2(10c_{1}+c_{2})+(2c_{1}+10c_{2}))\}

and

\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}) .

\frac{\lambda^{n+1}}{\lambda^{n}}=\frac{2a_{1}+a_{2}+\frac{Sv}{12}(2c_{1}+c_{2})}{2a_{1}+a_{2}-\frac{Sv}{12}(2c_{1}+c_{2})}-\frac{(2a_{1}+a_{2})}{2a_{1}+a_{2}-\frac{Sv}{12}(2c_{1}+c_{2})}\sum_{k=0}^{n}b_{k}^{\alpha}(\lambda^{1-k}-\lambda^{-k}).

Now, we apply mathematical induction. For $n = 0$ 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})} .

\frac{\lambda^{1}}{\lambda^{0}}=\frac{2a_{1}+a_{2}+\frac{Sv}{12}(2c_{1}+c_{2})}{2a_{1}+a_{2}-\frac{Sv}{12}(2c_{1}+c_{2})}.

$| \frac{λ^{n + 1}}{λ^{π}} | \leq 1$ |\frac{\lambda^{n+1}}{\lambda^{\pi}}|\le1, $| 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}) |$ |2a_{1}+a_{2}+\frac{Sv}{12}(2c_{1}+c_{2})|\le|2a_{1}+a_{2}-\frac{Sv}{12}(2c_{1}+c_{2})| and $2 c_{1} + c_{2} \leq 0$ 2c_{1}+c_{2}\le0. 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}

\begin{matrix}c_{1}\approx\frac{1}{h^{2}}-\frac{1}{10}+\frac{113h^{2}}{16800}-\frac{197h^{4}}{739200}+O(h^{6})\\ c_{2}\approx-\frac{2}{h^{2}}+\frac{1}{6}-\frac{29h^{2}}{1080}+\frac{19h^{4}}{10080}+O(h^{6})\end{matrix}

and it is clear that for sufficiently small $h (i . e ., h \to 0)$ h\rightarrow0, $2 c_{1} + c_{2} \leq 0$ 2c_{1}+c_{2}\le0 Thus, for $n = 0$ n=0, $| λ^{1} | \leq | λ^{0} |$ |\lambda^{1}|\le|\lambda^{0}|. Now, assume that $| λ^{l} | \leq | λ^{0} |$ |\lambda^{l}|\le|\lambda^{0}| for $l = 1, 2, \dots, n$ l=1,2,\cdot\cdot\cdot,n We now show that it also holds for $l = n + 1$ l=n+1 Using the assumption $| λ^{l} | \leq | λ^{0} |$ |\lambda^{l}|\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$ 0<\alpha<1, we can conclude from Eq. (37)

| \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}^{α}

|\frac{\lambda^{n+1}}{\lambda^{n}}|\le|\frac{2a_{1}+a_{2}+\frac{Sv}{12}(2c_{1}+c_{2})}{2a_{1}+a_{2}-\frac{Sv}{12}(2c_{1}+c_{2})}|+2|\lambda^{0}||\frac{(2a_{1}+a_{2})}{2a_{1}+a_{2}-\frac{5v}{12}(2c_{1}+c_{2})}|\sum_{k=0}^{n}b_{k}^{\alpha}

Thus, since $2 c_{1} + c_{2} \leq 0$ 2c_{1}+c_{2}\le0 we get

| λ^{n + 1} | \leq L | λ^{0} |

|\lambda^{n+1}|\le L|\lambda^{0}|

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=|\frac{2a_{1}+a_{2}+\frac{5v}{12}(2c_{1}+c_{2})}{2a_{1}+a_{2}-\frac{5v}{12}(2c_{1}+c_{2})}|+2|\frac{(2a_{1}+a_{2})}{2a_{1}+a_{2}-\frac{5v}{12}(2c_{1}+c_{2})}|\sum_{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:

\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}

\begin{matrix}<italic>L<sub>2</sub></italic>=||u(x,t)-S(x,t)||_{2}=\sqrt{h\sum_{j=1}^{N}|u(x_{j},t_{n})-S(x_{j},t_{n})|^{2}}\\ <italic>L<sub>∞</sub></italic>=||u(x,t)-S(x,t)||_{\infty}=max_{j}|u(x_{j},t_{n})-S(x_{j},t_{n})|.\end{matrix}

Additionally, the convergence rate of the method is calculated with the formula

R o C = \frac{\ln (E r r o r (N_{2}) / E r r o r (N_{1}))}{\ln (N_{1} / N_{2})}

RoC=\frac{ln(Error(N_{2})/Error(N_{1}))}{ln(N_{1}/N_{2})}

where $E r r o r (N_{1})$ Error(N_{1}) and $E r r o r (N_{2})$ Error(N_{2}) represent the errors at number of partition $N_{1}$ N_{1} and $N_{2}$ N_{2}, 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 below. As the first test problem, we consider the time-fractional Burgers equation given by (1) subject to the following initial and boundary conditions:

\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}

\begin{aligned} u(0,t) &= u(1,t) = 0, & t &\in [0,T], \\ u_x(0,t) &= u_x(1,t) = 2\pi t^2, & t &\in [0,T], \\ u(x,0) &= 0, & x &\in [0,1]. \end{aligned}

and the corresponding source term is

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) = \frac{2t^{2-\gamma} \sin(2\pi x)}{\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²sin(2π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^-3 to 1.71989088×10^-5, 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 B-spline 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₂	1.17608481×10^-3	1.16808997×10^-3	1.16122058×10^-3	1.15977847×10^-3
	L_∞	1.56691086×10^-3	1.55586455×10^-3	1.54633448×10^-3	1.54428531×10^-3
20	L₂	3.24966594×10^-4	3.22510891×10^-4	3.21299945×10^-4	3.23144788×10^-4
	L_∞	4.58600285×10^-4	4.55139112×10^-4	4.53437127×10^-4	4.56047172×10^-4
40	L₂	7.78030654×10^-5	7.69389503×10^-5	7.73445392×10^-5	8.01209788×10^-5
	L_∞	1.09833692×10^-4	1.08615283×10^-4	1.09189833×10^-4	1.13110804×10^-4
80	L₂	1.25916352×10^-5	1.21469751×10^-5	1.29784520×10^-5	1.60002218×10^-5
	L_∞	1.78288364×10^-5	1.71989088×10^-5	1.83758881×10^-5	2.26544495×10^-5

Table 2

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

Δ t	Norm	γ=0.25	γ=0.5	γ=0.75	γ=0.9
0.002	L₂	2.79799013×10^-5	2.89453357×10^-5	2.39205481×10^-5	1.09172130×10^-5
	L_∞	3.96292257×10^-5	4.09955552×10^-5	3.38827254×10^-5	1.54759469×10^-5
0.001	L₂	9.09662936×10^-6	1.28824323×10^-5	7.53455829×10^-6	1.19234428×10^-6
	L_∞	9.70022354×10^-6	1.37368370×10^-5	1.06714152×10^-5	1.69358924×10^-6
0.0005	L₂	1.18088684×10^-8	3.57839140×10^-7	9.03171600×10^-7	3.97091243×10^-6
	L_∞	2.09713397×10^-8	5.05989927×10^-7	1.27789156×10^-6	5.62041476×10^-6

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₂	8.18×10^-5	1.22×10^-3	1.60×10^-5
	L_∞	1.15×10^-4	1.73×10^-3	2.63×10^-5
80	L₂	1.70×10^-5	1.78×10^-4	7.72×10^-6
	L_∞	2.40×10^-5	2.53×10^-4	1.34×10^-5
100	L₂	9.14×10^-6	5.23×10^-5	7.24×10^-6
	L_∞	1.29×10^-5	7.65×10^-5	1.19×10^-5

Table 4

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

Δ t	Norm	Proposed method	Ref. [10]	Ref. [11]
0.002	L₂	2.89×10^-5	1.22×10^-5	5.55×10^-5
	L_∞	4.09×10^-5	1.72×10^-3	8.01×10^-5
0.001	L₂	9.70×10^-6	5.32×10^-4	2.77×10^-5
	L_∞	1.37×10^-5	7.53×10^-4	3.92×10^-5
0.0005	L₂	1.18×10^-8	1.89×10^-4	1.39×10^-5
	L_∞	2.09×10^-8	2.68×10^-4	2.05×10^-5

To investigate the convergence rate of the method for Problem 1, the convergence ratios are presented in Table 5 and Table 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 Ar 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 At 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 = 1 for different values of N and Ar.

N	Δ t	L_∞	ROC
6	1/5	4.51974778×10^-3	-
12	1/20	1.64465660×10^-4	4.78
24	1/80	7.17815902×10^-6	4.51
48	1/320	7.21429488×10^-7	3.31
96	1/1280	3.72153127×10^-8	4.27

Table 6

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

N	Δ t	L_∞	ROC
10	1/4	5.26044936×10^-3	-
20	1/32	3.50923812×10^-4	3.90
40	1/256	1.65112329×10^-5	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 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{array}{c} 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{array}

\begin{matrix} u(0,t) = t^2, u(1,t) = \exp(1)t^2, t \ge 0 \\ u(x,0) = 0, 0 \le x \le 1. \end{matrix}

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) = \frac{2t^{2-\gamma}}{\Gamma(3-\gamma)}\exp(x) + t^4\exp(2x) - vt^2\exp(x)

and the exact solution of 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 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 a decreasing behavior can be observed on the error norms as the derivative order y 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, a peak on error norms can be seen, however, for different values of fractional order y, 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₂	1.61504086×10^-4	1.44824109×10^-4	1.38338571×10^-4	1.15834339×10^-4
20	L_∞	3.08801388×10^-4	1.98520497×10^-4	1.89620989×10^-4	1.58724458×10^-4
40	L₂	8.55740423×10^-5	8.82849416×10^-5	8.27624785×10^-5	6.08117897×10^-5
40	L_∞	1.17333257×10^-4	1.21056287×10^-4	1.13520309×10^-4	8.34162090×10^-5
80	L₂	7.11925979×10^-5	7.41488030×10^-5	6.88664654×10^-5	4.70544770×10^-5
80	L_∞	9.76772696×10^-5	1.01737613×10^-4	9.45207070×10^-5	6.45980625×10^-5
100	L₂	6.94667643×10^-5	7.24524065×10^-5	6.71988959×10^-5	4.54035631×10^-5
100	L_∞	9.53080411×10^-5	9.94124999×10^-5	9.22404235×10^-5	6.23420557×10^-5

Table 8

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

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

Table 9

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

N	Methods	L₂	L_∞
40	Proposed method	5.36×10^-5	7.36×10^-5
40	Ref. [10]	6.77×10^-5	2.09×10^-4
80	Proposed method	3.95×10^-5	5.42×10^-5
80	Ref. [10]	4.57×10^-5	6.92×10^-5

Table 10

Comparison of numerical solutions and errors for γ=0.5, Δt=0.00025, t = 1, v=1 and N=40 for Problem 2.

X	Proposed method	L_∞	Ref. [12]	L_∞	Exact solution
0.2	1.221366	3.66×10^-5	1.221462	5.95×10^-5	1.221402
0.4	1.491762	6.24×10^-5	1.491934	1.09×10^-4	1.491824
0.6	1.822045	7.36×10^-5	1.822258	1.39×10^-4	1.822118
0.8	2.225480	6.05×10^-5	2.225666	1.025×10^-4	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 At for the values of the parameters γ=0.9, v=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^{-2} to 1.11304212×10^-4 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^{-2} to 1.08765049×10^-4 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 Ar for γ=0.9, v=1, t_f=1 for Problem 2.

N	Δ t	L_∞	ROC
30	1/4	2.20295625×10^{-2}	-
60	1/64	1.21511503×10^-3	4.18
120	1/1024	1.11304212×10^-4	3.44

Table 12

The error norms Le and convergence rates for varying values of N and Ar for γ=0.9, v=1, t_f=1 for Problem 2.

N	At	L_∞	RoC
100	1/4	2.20820268×10^{-2}	-
200	1/64	1.20495216×10^-3	4.19
400	1/1024	1.08765049×10^-4	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 is 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 graph 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 Figure 5 (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 well respected.

5.3

Example 3

As a 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}

\begin{aligned} u(x,0) &= 0, & 0 \le x \le 1 \\ u(0,t) &= t^2, u(1,t) = -t^2, & t \ge 0. \end{aligned}

For the problem 5.3, the corresponding source term and exact solution are given as

f (x, t) = \frac{2 t^{2 - γ} \cos (π x)}{Γ (3 - γ)} - π t^{4} \cos (π x) \sin (π x) + v π^{2} t^{2} \cos (π x)

f(x,t) = \frac{2t^{2-\gamma} \cos(\pi x)}{\Gamma(3-\gamma)} - \pi t^4 \cos(\pi x) \sin(\pi x) + v \pi^2 t^2 \cos(\pi x)

u (x, t) = t^{2} \cos (π x) .

u(x,t) = t^2 \cos(\pi x).

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 v=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 y rises. Table 14 illustrates an overall trend of diminishing error norms as the time step lowers and the value of y approaches 1.

Table 13

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

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

Table 14

The error norms L₂ and Los for varying values of y and Ar for N=120 of Problem 3.

Δ t	Norm	γ=0.25	γ=0.5	γ=0.75	γ=0.9
0.002	L₂	3.14766794×10^-5	3.22741925×10^-5	2.77834353×10^-5	1.62793905×10^-5
0.002	L_∞	4.39151250×10^-5	4.50291619×10^-5	3.87661038×10^-5	2.27160596×10^-5
0.001	L₂	1.46307138×10^-5	1.51160676×10^-5	1.31564294×10^-5	7.5330880×10^-6
0.001	L_∞	2.04121181×10^-5	2.10899969×10^-5	1.83571994×10^-5	1.05120674×10^-5
0.0005	L₂	6.23879040×10^-6	6.51771718×10^-6	5.66868928×10^-6	2.93671313×10^-6
0.0005	L_∞	8.70409358×10^-6	9.09364603×10^-6	7.90976865×10^-6	4.09826535×10^-6

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 y = 0.5, Δ t=0.00025, t_f=1 and different values of N (Problem 3).

N	Method	L₂	Loo
40	Proposed method	1.47×10^-5	2.05×10^-5
	Ref. [10]	1.22×10^-3	1.73×10^-3
	Ref. [11]	1.60×10^-5	2.63×10^-5
80	Proposed method	4.41×10^-7	6.15×10^-7
	Ref. [10]	1.78×10^-4	2.53×10^-4
	Ref. [11]	7.72×10^-6	1.34×10^-5
100	Proposed method	1.27×10^-6	1.77×10^-6
	Ref. [10]	5.23×10^-5	7.65×10^-5
	Ref. [11]	7.24×10^-6	1.19×10^-5

Table 16

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

Δ t	Method	L₂	L_∞
0.002	Proposed method	3.22×10^-5	4.50×10^-5
	Ref. [10]	1.22×10^-3	1.72×10^-3
	Ref. [11]	5.55×10^-5	8.01×10^-5
0.001	Proposed method	1.51×10^-5	2.10×10^-5
	Ref. [10]	5.32×10^-4	7.53×10^-4
	Ref. [11]	2.77×10^-5	3.92×10^-5
0.0005	Proposed method	6.51×10^-6	9.09×10^-6
	Ref. [10]	1.89×10^-4	2.68×10^-4
	Ref. [11]	1.39×10^-5	2.05×10^-5

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, v = 1, t t_f=1 and different values of Ar and N (Problem 3).

N	Δ t	L_∞	ROC
8	1/5	6.43241077×10^-3	-
16	1/40	4.49763147×10^-4	3.83
32	1/320	3.33582433×10^-5	3.75

Table 18

Maximum errors and convergence rates for y = 0.5, v = 1, tf = 1 and different values of Ar and N (Problem 3).

N	Δ t	L_∞	ROC
5	1/4	7.90193068×10^-3	-
10	1/32	3.54754064×10^-4	4.47
20	1/256	1.22189859×10^-5	4.85

Table 19

Maximum errors and convergence rates for y = 0.5, v=1, t_f=1 and different values of Ar and N (Problem 3).

N	Δ t	L_∞	ROC
12	1/5	1.64465660×10^-4	-
24	1/40	7.17815902×10^-6	3.64
48	1/320	7.21429488×10^-7	3.25
96	1/2560	3.72153127×10^-8	3.55

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, 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 time-fractional Burgers equations

Full Article

Paradigm

My account