A high order method for numerical solution of time-fractional KdV equation by radial basis functions

A radial basis function method for solving time-fractional KdV equation is presented. The Caputo derivative is approximated by the high order formulas introduced in Buhman (Proc. Edinb. Math. Soc. 36:319–333, 1993). By choosing the centers of radial basis functions as collocation points, in each time step a nonlinear system of algebraic equations is obtained. A fixed point predictor–corrector method for solving the system is introduced. The efficiency and accuracy of our method are demonstrated through several illustrative examples. By the examples, the experimental convergence order is approximately 4-α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$4-\alpha $$\end{document}, where α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} is the order of time derivative.


Introduction
In this paper, we consider the time-fractional KdV equation with the following initial condition: Gaussian (GA) φ(r ) = exp(−cr 2 ) Hardy multiquadric (MQ) φ(r ) = √ r 2 + c 2 Inverse multiquadric (IMQ) φ(r ) = ( √ r 2 + c 2 ) −1 Inverse quadric (IQ) φ(r ) = (r 2 + c 2 ) −1 and boundary conditions: where ε and ν are the positive parameters, and h 1 (t), h 2 (t) and g(x) are known functions. Also, ∂ α u ∂t α is the Caputo fractional derivative and 0 < α < 1. The αth order Caputo fractional derivative of function F is defined as The Korteweg-de Vries equation was initially introduced to describe the lossless propagation of shallow water waves [17]. Also, it has been applied in plasma physics, radar, rheology [25,32], optical-fiber communications [37], super conductors [10], etc. Several numerical and analytical methods for solving fractional KdV equations have been introduced. In [29], the variational iteration method for the space-and time-fractional KdV equation was applied. These equations were investigated by Adomian method in [39]. The homotopy perturbation method has been used for fractional KdV equations in [38]. A numerical technique based on generalized Taylor series formula has been presented in [13]. For more information, for example, see [1,18,28,33,35].
The radial basis functions methods are a family of meshless methods. There are other meshless methods. Some of them are as follows: element-free Galerkin method, reproducing kernel particle method, point interpolation method, boundary element-free method, moving least-square approximation, etc. The meshless methods are new and interesting numerical techniques. They can solve many problems in physics and engineering that are not suited to conventional numerical methods, with a minimum of meshing or no meshing at all [2]. For more study about these methods, for example, see [3,9,14,16,[19][20][21]40].
The radial basis function methods are highly flexible and are efficient especially for solving problems with arbitrary geometry [30]. Furthermore, the RBF methods are performed without any mesh generation and they are usually more accurate than low order methods, such as finite differences, finite volumes and finite elements. To see the applications of RBFs in the numerical solution of partial differential equations (PDEs) and fractional PDEs, for example, see [7,11,15,24,27,36].
In Table 1, some well-known globally supported RBFs have been listed. Let r be the Euclidean distance between x * ∈ R d and any x ∈ R d , i.e., x − x * 2 . A radial function φ * = φ( x − x * 2 ) depends only on the distance between x ∈ R d and fixed point x * ∈ R d . Hence, the RBF φ * is radially symmetric about x * . Clearly, the functions in Table 1 are globally supported, infinitely differentiable and depend on a free parameter c. The parameter c is called shape parameter. It allows us to broaden the solution space outside of the polynomial space. Many researches are done for finding the optimal values of c which produce most accurate interpolation (e.g., see [5,8]). But this is still an open problem. In practice, the shape parameter c is chosen experimentally. In fact, the problem is solved by the different values of c. Those values of c which give convergent sequences of approximations are suitable.
The global infinitely differentiable RBFs can be used for interpolating smooth data with a spectral accuracy [4,5,23,26]. More information about the accuracy of the approximations made by RBFs can be found in [12,41]. Let x 1 , x 2 , . . . , x M be a given set of distinct points in R d . The idea behind the use of RBFs is interpolation by translations of a single function, i.e., where Since all φ of the interest have global support, this method produces a dense matrix A. The matrix A corresponding to GS, IMQ and IQ for distinct interpolation points is positive definite, and therefore nonsingular [34].
In the above cases, the condition number of A is usually a very large number and A is very ill-conditioned. Therefore, in our computations more precision arithmetics than the standard floating point arithmetic must be used. In this work, the high order difference formulas introduced in [6] are applied for discretizing on time variable. In each time step, the solution of Eqs. (1)-(4) is approximated by a linear combination of RBFs with unknown coefficients. For finding these coefficients, we choose the centers of RBFs as collocation points, and consequently a nonlinear system of algebraic equations is obtained. The nonlinear system cannot be easily solved by the ordinary methods (for example, Newton method). So a fixed point iteration method for solving the resulted nonlinear system is proposed. By the fixed point method, the computations of the nonlinear system are reduced to some linear systems of equations.
The organization of the paper is as follows: In Sect. 2, a fixed point iteration method for solving the systems of nonlinear algebraic equations is introduced. In Sect. 3, the solution of Eq. (1) by RBFs is considered. Section 4 is devoted to numerical experiments.

Fixed point method
By the fixed point method first the system of nonlinear equations, is rearranged as for F, G : R n → R n . Then by a suitable initial approximation X 0 ∈ R n and the following recursive formula, successive approximations for the solution of (9) are given We will use the fixed point method as follows. First, we decompose (9) into an invertible linear (F 1 ) and a nonlinear (F 2 ) operators as Then which results Therefore, with the formulation, the desired sequence is obtained. The method is a predictor-corrector method. There is a sufficient condition for the convergence of the iteration formula (11) [31].
Theorem 2.1 Suppose that G : D ⊂ R n → R has a fixed point X * in the interior of D and that G is continuously differentiable in a neighborhood of X * . Denote by J G the Jacobian matrix of G and assume that ρ(J G (X * )) < 1. Then there exists a neighborhood S of X * such that S ⊂ D and for any X 0 ∈ S the iterates defined by (11) all lie in D and converge to X * [31].

Description of the method
First, we discretize Eq. (1) in the time direction as where u(x, t n ) = u n , t n = nτ , n = 1, 2, . . . , N , the time step τ , and the time length N τ . The value of ∂ α u(x,t n ) ∂t α for n = 1, n = 2 and n ≥ 3 is obtained as follows: in which and By Eqs. (16)- (19), the following finite differences equations are obtained: and Now using the radial basis functions, we consider the solution u(x, t n ) as where λ n j , j = 1, . . . , M are unknown. To construct the approximations for u 1 (x), first we substitute (23) in (20) and in boundary conditions (3) and (4), and then we collocate the resulted equations. For suitable collocation points, we choose the centers, M−1 , as collocation points. Thus, a nonlinear system of M equations in M unknowns is obtained. By solving the system, λ 1 (21), (3) and (4), and using the collocation method with the same collocation points. Inductively, u n (x), n = 3, 4, 5, ... is approximated by Eqs. (23), (22), (3) and (4) and with the same technique.
We solve the resulted nonlinear systems by utilizing the fixed point method presented in Sect. 2, as follows: In nth, n = 1, 2, . . . time step, the unknown vector, n = λ n 1 , λ n 2 , . . . , λ n M T , is given bỹ with a suitable choice of n 0 . In (24),F n = F n 1 ,F n 2 , . . . ,F n M T is a function of n and for n = 1, n = 2 and n ≥ 3 is as follows: For i = 2, . . . , M − 1,

Numerical illustrations
In this section, three test examples have been presented to show the efficiency of the mentioned method. We performed our computations using Maple 16 software. In all results, δ = 50 floating point arithmetics are used in our computations. Solution of the resulted nonlinear systems is obtained via the corresponding fixed point iteration method. Here, we consider stop condition of fixed point iteration method as follows: where all results are obtained with ξ = 10 −5 . Also, we have ρ(J G (X * )) < 1 and the fixed point method introduced in Sect. 2, gives convergent sequences. The errors and the experimental convergence order (C − Order) are calculated as follows: with the initial and boundary conditions, The exact solution of this problem is u(x, t) = t 5+α sech(x). We solved the problem by the new method and using GS-RBFs, IQ-RBFs, MQ-RBFs and IMQ-RBFs. Also, we put ε = 6, ν = 1 (except in Table 6). Table 2 presents the E ∞ errors and the experimental convergence orders (C − Order) obtained by our method with M = 21. In this case, ρ(J G (X * )) < 0.37. Table 2 shows the C − Order is approximately 4 − α. Table 2, demonstrates that as τ becomes smaller the smaller errors are obtained. Also, with τ = 0.1, 0.05, 0.025 and 0.0125, and respectively, 10, 20, 40 and 80 iterations the very small errors are obtained. So, we can conclude the method has good stability for solving the problem. In Table 3, we report the E ∞ , E 2 errors and ρ(J G (X * )) of present method using GS and MQ-RBFs for different values of α. Table 3 shows that the smaller the values of α, the more accurate results are obtained.
In Table 4, we report the E ∞ , E 2 , RMSE errors and condition number of matrix B =B n , n ≥ 3 resulted by 160 iterations of the new method, for some values of shape parameter c at t = 1. It seems that the optimal value of c is in the interval [2.3, 2.55]. But it cannot be found exactly. Table 5 presents the E ∞ , E 2 , RMSE errors and κ ∞ (B) for different values of M. In this case ρ(J G (X * )) < 0.26. Table 5 shows that as the number of RBFs becomes larger, the smaller errors are obtained while the condition number of B becomes larger. So, the dimension of matrix B should be small sufficiently to guarantee the stability of the solution. Table 6 presents the E ∞ and RMSE errors for ε = 1 and various values of ν resulted by GS-RBFs and MQ-RBFs. As Table 6 shows, the larger the values of ν the more accurate results are obtained. We plot the approximate solution resulted by 200 iterations of the mentioned method and M = 41 IQ-RBFs with τ = 0.005 and c = 1.25 in Fig. 1. In Fig. 2 Table 6 The The exact solution is u(x, t) = 1 3000 t 5 e −x 2 . Table 7 demonstrates the errors of numerical approximations resulted by our method and M = 31 MQ-RBFs with τ = 0.2, 0.1, 0.05 and 0.025, at t = 1, 2, 3, 4 and 5. The table shows that the method has good stability. Table 8 depicts the errors of the approximations resulted by the new method with different values of M, further to listing the corresponding condition number of the matrix B. Table 8 shows that as the number of RBFs becomes larger, the more accurate approximations are obtained while κ ∞ (B) becomes larger. In Fig. 3, we have plotted the absolute error function in −3 ≤ x ≤ 3 for 0 ≤ t ≤ 1 and 5 ≤ t ≤ 6.    We consider − tanh 2 (x) + 5sech 2 (x) − tanh 2 (y) + sech 2 (y) , and the initial and boundary conditions are as u(a 1 , y, t) = t 5+α sech(a 1 )sech(y), t ≥ 0, u(a 2 , y, t) = t 5+α sech(a 2 )sech(y), t ≥ 0, The exact solution of this problem is u(x, t) = t 6 sech(x)sech(y). We developed our method for solving this two-dimensional KdV problem. We chose the centers of RBFs (and the collocation points) as ( . . , M 1 and j = 1, . . . , M 2 . Table 9 depicts the errors and the experimental convergence orders resulted by our method with M 1 = M 2 = 9 for four different time steps and [a 1 , In these cases, we got ρ(J G (X * )) < 0.48. Table 10 presents the E ∞ , E 2 , RMSE errors, ρ(J G (X * )) and κ 2 (B) resulted by IQ-RBFs for the case α = 0.35 and c = 1.95. In Fig. 4, we plot the error function |u exact − u approx. | and the approximate solution resulted by 50 iterations of our method with MQ-RBFs and τ = 0.02, at t = 1.

Conclusion
In this article, an RBF collocation method was used for finding the solution of the time-fractional KdV equation in Caputo sense. The RBFs were applied for discretization of the spatial variable and we chose the centers of RBFs as collocation points. The three formulas introduced in [6] were applied for approximating ∂ α u(X,t n ) ∂t α for n = 1, n = 2 and n ≥3. These formulas are of orders 2 − α, 3 − α and 4 − α, for, respectively, n=1, n=2 and n ≥3, where α is the order of derivative. Instead of the formulas in [6], we can use for example the difference formula [22], in which a 0 = τ −α (2−α) , b k = (k + 1) 1−α − k α , k = 0, 1, . . . , n. But the formulas in [6] have higher order of accuracy than the above formula, and consequently they lead to more accurate results. By our method, the computations of time-fractional KdV equation are reduced to some systems of nonlinear algebraic equations. A fixed point iteration method for solving the resulted nonlinear systems was introduced. Our method is computationally attractive and it can be generalized easily for two-and three-dimensional cases. In each time step, the method provides a closed form approximate solution. Three illustrative examples were included to demonstrate the accuracy, the stability and the convergence of the method and as we expected the experimental convergence order is approximately 4 − α.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.