Finite Difference / Collocation Method for a Generalized Time-Fractional KdV Equation

Abstract: In this paper, we studied the numerical solution of a time-fractional Korteweg–de Vries (KdV) equation with new generalized fractional derivative proposed recently. The fractional derivative employed in this paper was defined in Caputo sense and contained a scale function and a weight function. A finite difference/collocation scheme based on Jacobi–Gauss–Lobatto (JGL) nodes was applied to solve this equation and the corresponding stability was analyzed theoretically, while the convergence was verified numerically. Furthermore, we investigated the behavior of solution of the generalized KdV equation depending on its parameter δ, scale function z(t) in fractional derivative. We found that the full discrete scheme was effective to obtain a numerical solution of the new KdV equation with different conditions. The wave number δ in front of the third order space derivative term played a significant role in splitting a soliton wave into multiple small pieces.


Introduction
The study of nonlinear phenomena has always been an active subject in applied science and physics.As a classic nonlinear partial differential equation, the Korteweg-de Vires (KdV) equation is very important in analyzing waves on shallow water surfaces.In history, the KdV equation was first introduced by Boussinesq around 1877, and rediscovered by Diederik Korteweg and Gustav de Vries in 1895 [1].The KdV equation was firstly derived as an evolution equation that governs small-amplitude, long-surface gravity waves propagating in a shallow channel of water.However, the KdV equation can only be solved analytically in selected initial and boundary conditions due to its nonlinearity [2][3][4][5].With the development of computational techniques, the KdV equation gained much attention of mathematicians and physicists in the past few years.In [6], by choosing appropriate base functions, a Legendre-Petrov-Galerkin method is applied for the KdV equation.An optimal rate of convergence in L 2 -norm is derived without the nonlinear part.By using different test and trial function spaces and basis functions, [7] proposes a new dual-Petrov-Galerkin method for high odd-order differential equations.In [8], pseudo-spectral method is applied to solve the KdV equation, in which the spectral method is implemented parallelled with split-step.This is one of the original results of the pseudo-spectral method applied to the KdV equation.Later in [9], the KdV equation is solved by the homotopy perturbation method.By constructing special forms of initial conditions, numerical solutions expressed by the trigonometric functions, hyperbolic functions are obtained.Many authors applied the collocation method to solve KdV equation due to its high accuracy and good performance in the nonlinear case.In [10], the KdV equation is solved by using the radial basis functions.It is a meshless method based on the collocation with radial basis functions.Five standard radial basis functions are used in the method of the collocation.In [11], A KdV-Burgers equation was discussed by using the modified Bernstein polynomials.The B-polynomials were used to expand the desired solution requiring discretization with only the time variable.The Galerkin method was used to determine the expansion coefficients to construct initial trial functions.A fourth-order Runge-Kutta method was used to solve the system of equations for the time variable.In [12], the sinc-collocation method was applied to solve the nonlinear third-order KdV equation in the space direction with periodic forcing at the boundary numerically and a θ-weight finite difference method was used to approximate the first-order time derivative.However, all the work above is on the KdV equation with integer-order derivatives.In the recent forty years, many partial differential equations are reconsidered by changing their classical derivative into fractional derivatives.
Fractional derivative is not a new concept.In fact, it was firstly proposed in the discussion of Leibniz and l'Hospital in 1695.Since 1974, when the first international conference of fractional calculus and its application was held, many researchers from scientific and engineering fields brought fractional derivatives back to their research interests.A large number of physical models were renewed by fractional calculus.Under such a background, fractional KdV equation, as well as many other fractional PDEs, are also studied by different methods, either analytical or numerical techniques [13][14][15][16].In [17], a fractional KdV equation with a Caputo derivative is solved by the Adomian decomposition method.The approximation solution is presented via a convergent series with easily computable components.In [18], based on Maple, the Adomian decomposition method is extended to derive explicit and numerical solutions of the fractional KdV-Burgers equation.In [19], the homotopy analysis method that was developed for integer-order differential equations is directly extended to derive explicit and numerical solutions of the nonlinear fractional KdV equation.In [20], the homotopy perturbation method is implemented to construct solitary solutions for variants of the KdV equations with fractional time derivatives.In this scheme, the solution takes the form of a convergent series, and its terms can be easily calculated via iteration.In a recent work [21], the explicit and approximate solutions of the nonlinear fractional KdV-Burgers equation with time/space-fractional derivatives were presented and discussed.The numerical method is a technique based on the generalized Taylor series formula which constructs an analytical solution in the form of a convergent series.Noting that most of above work is analytical/semianalytical, we are interested in solving the fractional KdV equation numerically by a finite difference/collocation method, and trying to discover some particular features of the soliton solution of the fractional KdV equation.
Fractional derivatives have various definitions, out of which the frequently used ones are defined in Riemann-Liouville and Caputo senses, which can be regarded as a convolution between a fractional power kernel and a function.These kinds of fractional power kernel put more weight in the present than in the past, which is particularly suitable to describe physical processes in the material with short memory, but is incapable of properly describing some behavior observed in materials with huge heterogeneities and structures with different scales.In [22], a new kind of fractional KdV model was established by using the developed Caputo-Fabrizio fractional derivative, and the conditions of existence and uniqueness of the continuous solution were provided.In [23], a time-fractional Shamel-KdV equation with Riesz fractional derivatives was derived to describe nonlinear behavior of dust-ion acoustic waves.In 2012, to unify different existing fractional derivatives, a new fractional calculus was proposed in [24].Taking a scale function and a weight function into consideration, fractional integrals and derivatives were essentially generalized.Partial differential equations with such generalized fractional derivatives are not studied much [25][26][27][28].Therefore, in this paper, we choose the KdV equation to further study this topic.By using a generalized derivative with respect to scale function and weight function, we defined a new time-fractional KdV equation.As there is strong nonlinearity in the structure of the KdV equation, the collocation method is applied to the space direction while the finite difference method is implemented in the time direction.Because of the complexity of generalized fractional derivative, some classic mathematical tools are not available anymore in studying the generalized KdV equation, such as variational principle, energy functional, Hölder estimation and Sobolev-type inequalities.We only have very limited methods of theoretical analysis to detect such KdV equations.However, in the current circumstance, the numerical method shall be a useful and direct way to show the details of the considered KdV equation with generalized fractional derivatives.After applying the collocation method in the space direction, the generalized fractional KdV equation is reduced to a nonlinear generalized fractional differential equation with respect to t. Reference [29] studied the existence and uniqueness of a class of generalized fractional boundary value problems by using topological degree theory.In our paper, we focus on numerical methods to analyze such a KdV equation.Thus, we choose initial and boundary conditions that are sufficiently good to guarantee the existence and uniqueness of the new time-fractional KdV equation.We apply a semi-discretization scheme to solve this new fractional KdV.Based on numerical simulation, many new characteristics of the new fractional KdV equation are presented.The remainder of this article is organized as follows: Section 2 is devoted to the preliminaries of generalized fractional calculus.In Section 3, we discuss the numerical scheme and its stability.Two numerical examples are given in Section 4. The first one is calculated to confirm the effectiveness and accuracy of the proposed numerical method.The second one with soliton-type initial condition is used to show important differences between the conventional KdV equation and our generalized fractional KdV equation.Some nonlinear scale functions are considered, which makes the KdV equation more general.

Generalized Fractional Calculus
In this part, we introduce the definitions of new generalized fractional integrals and derivatives.To guarantee that the new generalized fractional integral and derivative are well-defined, we need to assume that the weight function is positive and the scale function is monotonically increasing on the domain.

Definition 1. ([24]
) The left-generalized fractional integral of order α > 0 of a function u(t) with respect to a scale function z(t) and a weight function w(t) is defined as provided the integral exists.The superscript '*' indicates the integral here is a new generalized fractional integral.

Definition 2. ([24]
) The left first-order generalized derivative of a function u(t) with respect to a scale function z(t) and a weight function w(t) is defined as provided the right side of the equation is finite.

Definition 3. ([24]
) The left-generalized derivative of order m of a function u(t) with respect to a scale function z(t) and a weight function w(t) is defined as provided the right side of the equation is finite, where m is a positive integer.

Definition 4. ([24]
) The left Caputo-type generalized fractional derivative of order α > 0 of a function u(t) with respect to a scale function z(t) and a weight function w(t) is defined as provided the right side of the equation is finite, where m − 1 ≤ α < m, m is a positive integer.
Similarly, for one function containing two variables, that is, u(x, t), its (left) Caputo-type generalized fractional partial derivative with respective to t is defined as [25] ( Remark 1. Left Riemann-Liouville-type generalized fractional derivatives can be defined if we switch the order of integral and differential.When the scale function is selected as z(t) = t and the weight function is a nonzero constant, the operators defined in (1), ( 4) and left Riemann-Liouville-type generalized fractional derivatives are reduced to left Riemann-Liouville fractional integrals, left Caputo fractional derivatives and left Riemann-Liouville fractional derivatives, respectively.

Jacobi Polynomials
In recent years, more concentration has been focused on Jacobi polynomials J β,γ n for solving differential equations because the Jacobi parameters β and γ make such polynomials easy to implement.The well-known Jacobi polynomials J β,γ n with indices β, γ > −1 of degree n are defined on interval [−1, 1] and can be expressed as which are solutions of a class of Sturm-Liouville problems [30].
The integer derivative of J β,γ n on the interval [−1, 1] can be represented in the form Let x j , wj N j=0 be a set of Jacobi-Gauss-Lobatto nodes and weights.The corresponding discrete inner product, denoted by •, • N, w, can be defined as The exactness of the quadratures implies where P 2N−1 is denoted by the space of polynomials whose degree does not exceed 2N − 1.

Numerical Scheme and Stability Analysis
The integer-order KdV equation, is firstly derived to describe one-dimensional shallow-water waves with small but finite amplitudes.Generally, ξ is a positive parameter and δ > 0 indicates the wave number.Replacing the time derivative by the Caputo-type generalized fractional derivative (see ( 5)), we get which is called a generalized fractional KdV equation.Note that the second term of (10) makes it a nonlinear equation.In order to solve this equation, we impose the following initial and boundary conditions In the following parts, as a premise, we make an assumption that the problem ( 10) and ( 11) is well-posed.

Numerical Discretization in the Time Direction
To establish a numerical approximation scheme in the time direction, we divide the time domain into M equal parts, where M is a positive integer.Thus ∆t = T M is the time-grid size.We label the successive time nodes as , where t 0 = 0 is the initial time.Assuming w(t) > 0, 0 < α < 1 and z (t) is strictly monotonic increasing on (0, T], then s = z −1 (η) by setting η = z(s).The Caputo-type generalized fractional derivative of u(t) can be discretized as where and R k is the truncation error given by Proof.According to the mean-value theorem of integrals, there exist w (t k ) = w k and z (t k ) = z k for convenience.By using the Newton interpolation method, we interpolate U (η) and according to the theory of Newton interpolation, there exists Then, (14) can be rewritten as Denote [t i * −1 , t i * ] as the interval satisfying max 1≤i≤k (z i − z i−1 ) = z i * − z i * −1 , then the truncation error (14) has the form when z(t) satisfies the Lipschitz condition on [t i * −1 , t i * ] and L is the Lipschitz constant.
According to the finite difference approximation (12) and replacing the nonlinear part u (x, t k ) with the former step u (x, t k−1 ), we have the semi-discretization form of the generalized fractional KdV Equation (10) as where

Numerical Discretization in the Space Direction
As the homogeneous Neumann condition is included in the boundary conditions, we choose a type of polynomial as p (x) = (b − x) p (x) to approximate to the u (x, t) in the space direction.In our problem, we set that p (x) is a combination of Jacobi polynomials of degree ≤ N, and p (x) satisfies the Equation (10) exactly at x j , j = 1, • • • , N − 1 and boundary condition at x 0 and x N .x j , j = 0, • • • , N are Jacobi-Gauss-Lobatto (JGL) nodes on [a, b] by introducing the change of variable Based on the analysis above, we take as the basis functions in the space direction.Then, the numerical solution of problem (10) has the form: satisfying for any t ∈ [0, T].Substituting (24) to the semi-discretization ( 22), we have (27)  where f k (x) = f (x, t k ) and R s is the error in the space direction caused by replacing u (x, t) with u N (x, t).The classical first-and third-order partial derivatives in ( 27) are derived as [ψ n (x)] = Γ (n + 4 + γ + β) Neglecting the error parts, we obtain the full discretization scheme of Equation (10).For where (31)

Stability Analysis
Firstly, we consider the linear case of the generalized fractional KdV equation as and similarly, we discretize the linear case as Now we state our main theorem as follows.
Theorem 1. Assume that u (x, t) is a continuous function satisfying boundary conditions (11) for any t, and f (x, t) is continuous function.Letting u k N (x) = ∑ N n=0 φ n (t k ) ψ n (x) is the numerical solution of discrete scheme (33), the scheme (33) is stable unconditionally and for any k ≥ 0, it holds that where The discrete form (30) at JGL nodes x j can be rewritten as Multiplying both sides of (35) by v k N−2 x j ŵj and summing for j from 0 to N leads to where ŵj is the corresponding weight.Since the degrees of u k N (x), ∂ x u k N (x) and ∂ xxx u k N (x) are not exceeding N + 1, according to (8), they can be exactly expressed by the Jacobi-Gauss-Lobatto quadrature: which yields by applying Cauchy-Schwartz inequality and Lemma 1.
In the following, we will continue our proof by induction.Obviously, (34) holds when k = 0 according to (42): Assuming that (34) is valid for any i = 0, 1, • • • , k − 1, it is easy to verify from (42) that This completes the proof.
Remark 3.For the nonlinear case (10), at every collocation node x j , we have the 'frozen coefficient' ξu x j .Since Theorem 1 is sustained for any constant ξ, the discretization scheme (30) is stable by regarding ξu k−1 N x j as a constant coefficient.

Numerical Results and Analysis
In Section 3.3, we proved the stability of the proposed scheme in the linear case.Following numerical examples will show that our algorithm is also available in the nonlinear case.In this part, two numerical examples are provided to demonstrate the effectiveness of the finite difference/collocation scheme.The initial and boundary conditions are given to establish an exact solution in the first example to analyze error and convergence numerically.In the second example, the initial condition is taken as a hyperbolic sine function and the boundary condition is homogenous.For convenience, we take the weight function w (t) = 1 in the following part.
Applying the proposed finite difference/collocation scheme, we obtain the numerical solutions over the domain under different situations.Firstly, letting z (t) = t, the generalized fractional derivative in (10) is equivalent to the Caputo fractional derivative.Figure 2 shows the influence of δ and α under conditions of T = 2, ξ = 6.By fixing α = 0.8, diffusion processes with δ = 0.2 and δ = 1 are illustrated.It is shown that when δ is small, the soliton (initial condition) will be decomposed to neighboring solitons, whose amplitudes are affected by the fractional derivative α.This decomposing phenomenon cannot be observed when δ becomes larger.In the case of δ = 1, the soliton wave does not move much with 0 < t < 2 as expected, and with the decreasing of α, the amplitudes are getting smaller.Next, numerical results with different scale functions z (t) = t 0.5 , t , t 2 and t 5 are shown in Figure 3 when T = 1.Denoting z (t) = t s , 0 ≤ t ≤ 1, scale function z (t) behaves like a "contracting function" for the time variable over domain [0, 1] when s > 1, while behaving like a "stretching function" as s < 1 (see [25]).From Figure 3, when s = 0.5, decomposing and diffusion phenomena appeared earlier in contrast with the case when s = 1, just like being contracted.On the contrary, the decomposing process works quite slowly when s > 1.Moreover, the speed of decomposing is slower as s becomes larger.However, comparing with Figure 2a, Figure 3e,f shows that when t > 1, with the increasing of s, the speed of decomposing is getting faster.

Conclusions
We numerically study a new class of KdV equation, in which the partial derivative on the time variable is defined by a new generalized fractional derivative.The finite difference method is applied to solve the generalized fractional KdV equation under different selections of parameter, order of fractional derivatives and initial conditions.The stability of the numerical scheme is verified and the order of convergence is evaluated numerically in cases of different scale functions.Numerical examples are provided to demonstrate the effectiveness and accuracy of the numerical scheme.Furthermore, we investigate the dynamics of soliton solution in the generalized fractional KdV equation via simulation.We find that the number and amplitude of waves are significantly influenced by the parameter in the equation and the order of fractional derivatives.

Figure 1 Figure 1 .
Figure1shows the numerical solutions and the distribution of absolute error when α = 0.8, N = 20 and M = 50.Note that numerical solutions are highly aligned with exact solutions.