# A Meshfree Time-Splitting Approach for the Time-Fractional Burgers’ Equation

Authors: Erdal Korkmaz; Kenan Yildirim
Journal: Journal of Mathematics
 (2023)

---

## Abstract

We consider to represent an algorithm for time-fractional Burgers’ equation utilizing the multiquadric-radial basis functions with the time-splitting technique. This algorithm is performed on the three examples. The numerical results indicated that the algorithm yields accurate approximate solutions of the time-fractional Burgers’ equation. The physical behavior of the obtained solutions is given with graphics, and the results showed that the obtained solutions are in good match with the solutions reported in the literature. The algorithm is accurate, flexible, and easy to implement.


---

## Body

## 1. Introduction

We address the numerical solution of the one-dimensional fractional Burgers’ equation as follows:(1)∂αu∂tα+∂∂xu22=ν∂2u∂x2+gx,t,x,t∈0,1×0,T,with the initial and boundary conditions(2)ux,0=u0x,0≤x≤1,u0,t=u1t,u1,t=u2t,0≤t≤T,respectively, where 0<α≤1 is Caputo fractional derivative; t and x are time and space parameters, respectively; ν>0 is the viscosity constant; and u0x, u1t, and u2t are known functions, and u0x is sufficiently smooth. Here, ∂αu/∂tα are fractional derivatives in Caputo sense, defined by(3)∂αux,t∂tα≔1Γ1−α∫0tutx,τt−ταdτ.The steady-state solutions of the classical Burgers’ equation are primarily computed in [1]. Burger [2] proposed to utilize classical Burgers’ equation for the turbulence model. Since then, Burgers’ equation has always been the focus of attention for researchers. Particularly, the development of the fractional derivative in recent years has required the fractional modeling of classical mathematical models. Burgers’ (1) is important one of such models. The reason why fractional derivative has been a hot topic lately is that it finds application in science and engineering in fields such as physics, biology, and control theory [3, 4]. When it was difficult for researchers to find analytical solutions for fractional differential equation models, error analyzes have been made by searching numerical solutions in general. The authors of [5] presented an approximate scheme for fractional differential equations which is based on the generalization of the Adams–Bashforth–Moulton method and then provided an exhaustive error analysis of this scheme. The authors of [6] tailored a fully discrete difference scheme for a diffusion-wave system and it is based on the transforming the governing equation into a low-order system of equations by defining two new variables. In [7], the author proves that the smooth solutions exist for the Caputo-type fractional differential equations. The author in [8] presented a scheme for the fractional differential equation. This scheme is tailored by combining the predictor-corrector scheme with the short memory principle. The authors of [9] analyzed the stability of the explicit multistep methods and derived novel schemes which enable a larger stability region. The local truncation error and the stability of the high-order compact finite difference scheme for the fractional diffusion equation are investigated in [10]. The fractional linear multistep and the finite element techniques are performed for time and space discretization, respectively, for the time-fractional subdiffusion equation in [11]. The convergence and stability analysis of the finite difference schemes for fractional differential equations are investigated in [12]. A numerical scheme based on the combination of the Alikhanov scheme and the Galerkin finite element methods for solving nonlinear time-fractional parabolic equations is proposed in [13]. A linear finite difference scheme is tailored for the numerical solutions of the generalized time-fractional Burgers equation in [14]. In [15, 16], the authors used lower- and higher-order classical time-splitting methods to solve Burgers’ equation. The authors of [17] presented three different time-splitting algorithms with the convergence and stability analysis for nonlinear time-fractional differential equations. In [12–16, 18–27], the authors proposed numerical schemes such as finite difference, finite elements, spectral methods, and radial basis functions to analyze the numerical solutions of the fractional Burgers’-type equations.In this work, we propose an algorithm depending on the utilization of the multiquadric-radial basis functions and a splitting method for solving time-fractional Burgers’ equation numerically. By using advantage of a meshfree scheme in this algorithm, the goal is to get an approach that yields accurate results and exhibits the physical dynamics of the time-fractional Burgers’ equation.The remainder of this article is organized as follows. In Section2, the multiquadric-radial basis functions (MQ-RBF) used in this paper were introduced. Section 3 contains the time-splitting technique. In Section 4, the numerical results are presented to show the performance of the proposed algorithm. Finally, we give some concluding remarks in Section 5.


## 2. The Multiquadric-Radial Basis Functions (MQ-RBFs)

In order to obtain a fractional-order differential equation system, we discretize equation (1) in space by performing the radial basis functions method [28, 29]. To do this, one can represent the numerical solution of equation (1), ux,t, as follows:(4)ux,t≈∑j=0kβjtθjx,x∈ℝl,where θjx is continuously differentiable radial basis functions, k represents the number of distributed nodes, x=x0,x1,…,xl−1 is l dimensional vector, and βj are needed to be defined by utilizing the collocation methods. On the other hand, we use the following MQ-RBF function:(5)θjx=ξj2+c2,where c is the shape constant and ξj=x−xj denoted by the Euclidean norm.We define an approximate solutionux,t by ukx,t and then write(6)ukx,t=∑j=0kβjtθjx=ΘTxβ,where θjx=x−xj2+c2, xj=j/kj=0,1,…,k are uniformly distributed nodes on the range [0, 1], and(7)Θx=θ0x,θ1x,…,θkxT,β=β0t,β1t,…,βktT.Equation (6) can be rewritten in the following form:(8)Λβ=w,where w=u0t,u1t,…,ukt with ukxi,t=uit for i=1,…,k, and(9)Λ=ΘTx0ΘTx1⋮ΘTxk=θ0x0θ1x0…θkx0θ0x1θ1x1…θkx1⋮⋮⋱⋮θ0xkθ1xk…θkxk.By considering equations (4), (6), and (8), it follows(10)ux,t≈ukx,t=∑j=0kβjθjx=ΘTxβ=ΘTxΛ−1w=Φxw.If one replaces equation (10) and its derivatives with the corresponding terms in equation (1) at the collocation nodes xi=i/ki=0,1,…,k in [0, 1], then the following is obtained:(11)D0αCwt+wt∗Φxwt=νΦxx2Awt+gt,where D0αCwt, α∈0,1 is the Caputo derivative, w=u0,u1,…,uk,(12)gt=gx0,t,gx1,t,…,gxk,t,Φx=∂∂xΦjxik+1×k+1,Φxx=∂2∂x2Φjxik+1×k+1.Notice that, the operator “∗” is defined for point-wise multiplication. Furthermore, equations (11) read in the form as follows:(13)D0αCwt=Awt+ft,wt,where A=νΦxx,ft,wt=gt−wt∗Φxwt.


## 3. Time-Splitting Method

Splitting methods have received much attention due to their advantages such as speed, accuracy, and stability. These methods are explicit and easy to implement. Furthermore, splitting methods are designed as geometric numerical schemes. These schemes preserve some qualitative properties of the systems [30]. Let us consider following fractional differential equation:(14)D0αCwt=βwt+ρwt,w0=w0,0<t≤T,β,ρ∈ℂ,where D0αCwt,α∈0,1 is the Caputo derivative. The classical time-splitting method for one time step 0≤t≤h can be considered for equation (14) as follows:(15)D0αCw¯¯t=βw¯¯t,w¯¯0=w0,D0αCwt=ρw˜t,w˜0=w¯¯h.That is, equation (14) is split into two equations which are solvable and easy to deal with separately. These two equations are solved successively. Let us denote w˜ as approximate solution of fractional equation (14) and then the error of the classical splitting methods is given as follows:(16)wh−w˜h=w0∑k=0∞β+ρkhkαΓ1+kα−∑k=0∞βhαkΓ1+kα∑k=0∞ρhαkΓ1+kα=2βρΓ1+2α−βρΓ1+α2h2α+Oh3α.The order of convergence of the classical splitting method isOh2α. Thus, it can be said that the convergence order of the time-splitting scheme obtained for one step decreases as α approaches zero. In order to solve this problem, the author of [17] developed a division method for the fractional nonlinear systems of ordinary differential equations. The order of the convergence of this new division method is Oh1+α and we perform this method to solve equation (13), with α∈0,1, A is k+1×k+1 real matrices, f:0,T×Rk+1⟶Rk+1. Here, the continuous function f fulfills the Lipschitz inequality, that is, on suitable domain G,(17)ft,w1−ft,w2≤Kw1−w2,∀w1,w2∈G.The Caputo derivative is defined by(18)D0αCgt=Ia1−αg′t,Iaαgt=1Γα∫atgτt−τ1−αdτ,t>a.For convenience, the following notations are used [17]:(19)DaαsCgt=1Γ1−α∫atg′τs−ταdτ,Iaαsgt=1Γα∫atgτs−τ1−αdτ,Iaαsft,wt=1Γα∫atfτ,wτs−τ1−αdτ.To obtain splitting schemes, the time-fractional differential system equation (13) is written in the integral form as follows:(20)wt=w0+AI0αwt+I0αft,wt,0≤t≤T.If one definestn=nh, 0≤n≤N, as uniformly distributed points on the range 0,T with h=T/N, then equation (20) is written over tn−1,tnn≥1 as follows:(21)wt=wt0+Itn−1αAwt+It0αtAwtn−1+It0αft,wt.In order to tailor a splitting scheme, the author of [17] defined(22)w¯¯t=w˜tn−1+It0αtAw˜tn−1+Itn−1αAw¯¯t,tn−1<t≤tn,where the solution on 0,tn−1 is given and denoted with w˜t for 0≤t≤tn−1.The initial value ofw¯¯ at tn−1 is defined as follows:(23)w˜ctn−1=w˜tn−1+It0αtAw˜tn−1.Then, equation (22) reads(24)w¯¯t=w˜ctn−1+Itn−1αAw¯¯t,tn−1<t≤tn.With equations (21) and (22), one has(25)w˜t=w¯¯ctn+It0αft,w˜t,tn−1<t≤tn,where(26)w¯¯ctn=w¯¯t−w˜tn−1+wt0−Itn−1αAw¯¯−w˜t.The solutionw˜t for t∈tn,tn+1 can be computed with repetitions of above processes.Notice that from equations (23)–(26), one can reach equation (21). Thus, there is no splitting error.With equations (23)–(26), one writes(27)w¯¯tn=w˜ctn−1+Itn−1αAw¯¯tn,(28)w˜ctn−1=w˜tn−1+It0αtnAw˜tn−1,w˜0=w0,(29)w˜tn=w¯¯ctn+It0αftn,w˜tn,(30)w¯¯ctn=w¯¯tn−w˜tn−1+wt0−Itn−1αAw¯¯−w˜tn.The fully discrete time-splitting scheme can be tailored by discretization of equations (27)–(30). If one performs the weighted right-rectangle rule to equation (30), then the following is obtained:(31)w¯¯ctn≈w¯¯tn−w˜tn−1+wt0−1Γα∫tn−1tnAw¯¯tn−w˜tntn−τ1−αdτ,=w¯¯tn−w˜tn−1+wt0−AhαΓ1+αw¯¯tn−w˜tn.Using equation (31) with equation (29) yields(32)Chαw˜n=Chαw¯¯n+w0−w˜n−1+It0αftn,w˜tn,where w¯¯tn = w¯¯n, w˜tn = w˜n, and Chα=I−Ahα/Γ1+α. If one defines all known information in equation (32) by w¯¯nc, then the following is obtained:(33)Chαw˜n=w¯¯nc+Itn−1αftn,w˜tn,(34)w¯¯nc=Chαw¯¯n+w0−w˜n−1+It0αtnftn−1,w˜tn−1.On the other hand, the integrals in equations (28) and (34) are computed by performing the weighted trapezoidal rule and the weighted midpoint rule with w˜tj−1/2≈1/2w˜tj−1−w˜tj, respectively. However, the integrals in equations (27) and (33) are computed by utilizing the weighted right-rectangle rule and the weighted left-rectangle rule, respectively. Then, the following time-splitting scheme is obtained [17]:(35)w¯¯n=w˜n−1c+AhαΓ1+αw¯¯n,n=1,2,3,…,N,(36)w˜n−1c=w˜n−1+Ahα2Γ1+α∑j=1n−1zn,jαw˜j−1+w˜j,w˜0=w0,(37)Chαw˜n=w¯¯nc+hαΓ1+αftn−1,w˜n−1,n=1,2,3,…,N,(38)w¯¯nc=Chαw¯¯n+w0−w˜n−1+hαΓ1+α∑j=1n−1zn,jαftj−1/2,w˜j−1+w˜j2,where(39)zn,jα=αhα∫tj−1tj1tn−τ1−αdτ=n−j+1α−n−jα.The following approach is performed to get equation (36):(40)It0αtnAw˜tn−1=1Γα∫t0tn−1Aw˜τtn−τ1−αdτ=1Γα∑j=1n−1∫tj−1tjAw˜τtn−τ1−αdτ≈AΓα∑j=1n−1∫tj−1tjw˜j+w˜j+12tn−τ1−αdτ=Ahα2Γ1+α∑j=1n−1zn,jαw˜j−1+w˜j.This time-splitting scheme is called as TS-I and the following convergence theorem is valid [17].Theorem 1 (see [17]).
Letwt and w˜n, n∈0,N, be the solution of equation (13) and the schemes equations (35)–(38), respectively. Let h=T/N and Ahα≠Γ1+αI.Ifw∈C20,T and ft,wt fulfills the Lipschitz inequality equation (17), then following equation holds:(41)wtn−w˜n≤K1+w′0tn−αh1+α,1≤n≤N,with a constants K>0 that is independent of h. If f∈C2G, then(42)wtn−w˜n≤K1+tnα−1h1+α,1≤n≤N.


## 4. Numerical Results

To compare the performance of the mentioned algorithms for each test examples, the following absolute error andL2 and L∞ error norms are used:(43)unum−uL2=1k∑j=0kuj,N−uxj,t2,unum−uL∞=max0≤j≤kuj,N−uxj,t,Abs.Error=uj,N−uxj,t,0≤j≤k.We consider to assign the shape constant using trial and error. We have used errors between exact and numerical solutions to determine the best value of shape parameters. As shown in Figure1, the selection of its best value c=0.72 with corresponding smallest error. We compute L∞ error norm between exact and numerical solutions at the interval [0.01, 10] and compare them to specify the shape parameter.Figure 1
The best value of shape parameterc for Example 1 at T=0.1 with h=0.0005, k=10, and ν=0.1.Example 1.
We first compute the approximate solutions of problem (1) satisfying the conditions as follows:(44)ux,t=0=0,t≥0,ux=0,t=t2,ux=1,t=−t2,0≤x≤1,where(45)gx,t=2t2−αΓ3−α−π4t4sinπx+νπ2t2cosπx.
The exact solution is [19](46)ux,t=t2cosπx.
A comparison of the errors between presented method and [19, 20] for Example 1 is given in Table 1 at T=0.1 with h=0.0005. The results are obtained with k=10 and ν=1, ν=0.5, and ν=0.1 with corresponding shape parameter values c=0.7887, c=5.71, and c=0.72, respectively, for the present method. It is clear that as the viscosity value ν decreases, the errors become smaller as well as better than [19, 20]. In Table 2, we compare the numerical results with those of [31] for Example 1 at T=1 with h=0.02. For the presented scheme, the computation was handled using c=0.05, c=0.0517, c=0.0536 for α=0.25, α=0.75, α=0.9, respectively, and k=10. It indicates that the results obtained with the presented algorithm are in good match with the exact solutions and absolute errors are better than those given in [31]. We have tabulated the numerical solutions for Example 1 at T=1 with h=0.01, ν=0.1, and α=0.9 in Table 3. The shape parameters are c=3.933, c=5.138, and c=5.289 for k=30, k=50,and k=70, respectively. It shows that the approximate solutions are in good match with exact solutions and the absolute errors are small enough. Furthermore, Figures 2 and 3 show the behavior of the approximate solutions at T=0.1 and the absolute errors, respectively, with h=0.0005, c=0.72,v=0.1, k=10, and α=0.5. Figure 2 shows that the obtained solutions are symmetric.Table 1
Error norms of Example1 at T=0.1 for h=0.0005 and α=0.5.
ν=1ν=0.5ν=0.1L2L∞L2L∞L2L∞Reference [19] m′=106.980e − 069.547e − 066.492e − 068.854e − 064.288e − 065.714e − 06Reference [20] k=806.528e − 069.164e − 065.835e − 068.250e − 063.105e − 064.847e − 06Presentk=103.601e − 056.175e − 051.386e − 062.410e − 061.428e − 061.940e − 06Table 2
Results of Example1 at T=1 for h=0.02.
αx[31]PresentExactAbs. errorAbs. errork=50k=10[31]Present0.250.100.98121760.95791270.95105653.016115e − 026.856230e − 030.300.64640590.56801750.58778535.862068e − 021.976778e − 020.50−3.075e − 59−4.5661e − 166.1232e − 173.304634e − 595.178393e − 160.70−0.664067−0.5680175−0.58778535.862068e − 021.976778e − 020.90−9.812136−0.9579127−0.95105653.016115e − 026.856230e − 030.750.100.93198930.96433780.95105651.906739e − 021.328133e − 020.300.54188910.57467420.58778534.589585e − 021.311105e − 020.506.120e − 59−6.1287e − 156.1232e − 175.890384e − 596.189963e − 150.70−0.541856−0.5746742−0.58778534.589585e − 021.311105e − 020.90−0.931992−0.9643378−0.95105651.906739e − 021.328133e − 020.90.100.91464010.96286470.95105653.641245e − 021.180823e − 020.300.50614820.57574880.58778538.163693e − 021.203643e − 020.50−1.950e − 591.9096e − 156.1232e − 172.179615e − 591.848389e − 150.70−0.506129−0.5757488−0.58778538.163693e − 021.203643e − 020.90−0.914626−0.9628647−0.95105653.641245e − 021.180823e − 02Table 3
Results of Example1 at T=1, h=0.01, ν=0.1, and α=0.9.
kxPresentExactAbs. error300.100.95103830.95105651.823706e − 050.300.58776610.58778531.917656e − 050.70−0.5877274−0.58778535.787549e − 050.90−0.9510221−0.95105653.439243e − 05500.100.95101470.95105654.184683e − 050.300.58771170.58778537.358647e − 050.70−0.5877305−0.58778535.470573e − 050.90−0.9509895−0.95105656.705913e − 05700.100.95104080.95105651.568811e − 050.300.58773690.58778534.836005e − 050.70−0.5877776−0.58778537.668952e − 060.90−0.9510778−0.95105652.125912e − 05Figure 2
Numerical solutions of Example1 at T=0.1 with h=0.0005, c=0.72,v=0.1, k=10, and α=0.5.Figure 3
Absolute errors of Example1 at T=0.1 with h=0.0005, c=0.72, v=0.1, k=10, and α=0.5.Example 2.
In the next example, we deal with the time-fractional Burgers’ equation (1) subject to following condition:(47)ux,t=0=1−1+5coshx/2−sinhx/2,t≥0.
The exact solution is given as [32](48)ux,t=1−1+5coshx/2+tα/4Γ1+α−sinhx/2+tα/4Γ1+α,with force term(49)gx,t=0.
The boundary conditions are computed by using the exact solution.
In Table4, numerical solutions and absolute errors at T=0.02 are obtained by the present method for Example 2 with c=3.88, k=50, α=0.8, and h=0.0001. Table 4 shows that the numerical solutions of the presented methods are highly accurate than the solutions tabulated in [32] with the same parameter values.Table 4
Comparison of approximate solutions and absolute errors of Example2 with [32] at T=0.02 for h=0.0001, v=1, and α=0.8.
x[32]PresentExactAbs. error [32]Abs. error0.020.2563210.25690830.25690585.84566e − 042.49983e − 060.040.2595660.26016210.26015955.93809e − 042.62925e − 060.060.2628600.26346590.26346306.03243e − 042.96258e − 060.080.2662040.26682070.26681736.12874e − 043.42370e − 060.100.2696010.27022750.27022356.22707e − 043.96545e − 060.120.2730500.27368710.27368256.32749e − 044.53857e − 06Example 3.
As a final example, we investigate the numerical solutions of (1) satisfying the conditions as follows:(50)ux,t=0=sinπx,ux=0,t=ux=1,t=0,t≥0.
The series solutions are obtained in [33] by performing the Hopf-Cole technique [34, 35] as follows:(51)ux,t=2νπ∑m=1∞bme−m2π2νtαmsinmπxb0+∑m=1∞bme−m2π2νtαcosmπx,where(52)b0=∫01e−2πν−11−cosπxdx,bm=2∫01e−2πν−11−cosπxcosmπxdx,m=1,2,3….
In Table5, numerical solutions and absolute errors at T=0.1 are computed by the present scheme for Example 3 with k=10 and α=0.9. The shape parameters are chosen as c=0.76, c=0.05, and c=0.06 for h=0.01, h=0.001, and h=0.0001, respectively. Table 5 shows a good match between the obtained and the analytical results.Table 5
Approximate solutions and absolute errors of Example3 at T=0.1, for values of ν=1 and α=0.9.
hxPresentExactAbs. error0.010.100.10106470.07479992.62648e − 020.300.23878900.19842504.03640e − 020.700.23775800.20733673.04213e − 020.900.04040550.08031223.99067e − 020.0010.100.08029660.07479995.49670e − 030.300.19826250.19842501.62500e − 040.700.20409260.20733673.24410e − 030.900.08423720.08031223.92490e − 030.00010.100.07138700.07479993.41300e − 030.300.20048050.19842502.05560e − 030.700.20700950.20733673.27100e − 040.900.07515440.08031225.15780e − 03


## 5. Conclusions

A numerical algorithm based on the combination of the MQ-RBF and splitting technique is presented to investigate the solutions of the fractional Burgers’ equation, numerically. Instead of classical splitting methods, a division method for the fractional nonlinear systems of ordinary differential equations is performed. This method hasOh1+α order of convergence which is different from Oh2α. That is, when α tends to zero, the error does not increase significantly. Furthermore, this algorithm also takes advantage of a meshfree scheme. It is concluded that the present approach yields accurate results and exhibits physical dynamics of the time-fractional Burgers’ equation. This proposed scheme can be used for the other nonlinear fractional PDEs.



---

*Source: 5493451-2023-01-21.xml*