GiNaCDE: the high-performance F-expansion and First Integral Methods with C++ library for solving Nonlinear Differential Equations

We present the algorithms for three popular methods: F-expansion, modified F-expansion, and first integral methods to automatically get closed-form traveling-wave solutions of nonlinear partial differential equations (NLPDEs). We generalize and improve the methods. The proposed algorithms are manageable, straightforward, and powerful tools providing a high-performance evaluation of the exact solutions of nonlinear ordinary differential equations (NLODEs) and NLPDEs. For differential equations with parameters, the new algorithms determine the conditions on the parameters to obtain exact solutions. The algorithms show solutions to a wide variety of NLODEs and NLPDEs, both integrable and non-integrable. It can solve NLODEs and NLPDEs that contain complex functions. The algorithms are implemented in a C++ library named GiNaCDE. The efficiency and effectiveness of the algorithms are demonstrated by some examples with the help of GiNaCDE. The output results tally with the previously known results, and in some cases, new exact traveling-wave solutions are explicitly obtained. Use of the library, implementation issues, scope, limitations, and future extensions of the software are addressed.


I. INTRODUCTION
The Nonlinear Ordinary Differential Equations (NLODEs) and Nonlinear Partial Differential Equations (NLPDEs) play an important role in the theoretical sciences to explain many Nonlinear phenomena in various fields of science, such as biology, chemistry, engineer, solid-state physics, plasma physics, optical fibers, and so on. The exact (closedform) traveling-wave solutions (TWS) of such NLPDEs give much extra information, which helps us to study the result more deeply. The knowledge of closed-form solutions of NLODEs and NLPDEs helps to test the degree of accuracy of numerical solvers and also facilitates stability analysis. In the past few decades, many powerful methods have been presented for seeking exact solutions of NLPDEs. Some of them are tanh-expansion method [1][2][3], Jacobi elliptic function expansion method [4,5], multiple exp-function method [6], tanh-sech method [7], extended tanh method [8], hyperbolic function method [9], sine-cosine method [10], inverse scattering method [11], Hirota's direct method [12], F-expansion method [13], first integral method [14,15], etc. We have observed that most of the methods (such as first seven methods mentioned above) are function-expansion methods where the exact solutions of nonlinear differential equations restrict to polynomial solutions in terms of specific functions. For example the exact solutions are obtained in terms of the tanh function for tanh-expansion method [1][2][3] and Jacobi's sn (JacobiSN) or cn (JacobiCN) functions for Jacobi-expansion method [4,5]. These methods are useful for finding explicit traveling solitary wave solutions to nonlinear evolution equations.
On the other hand, F-expansion [13] and first integral methods [14,15] are different kinds of methods, which can overcome the limitations of the function-expansion method. F-expansion method was first introduced by Zhou et al. [13]. Later this method has been further improved [16][17][18][19][20][21][22][23][24], and successfully applied to many nonlinear equations for finding the exact solutions. In the F-expansion method, solutions of NLPDEs are expressed in finite power series of a new function F (ξ) which depends on traveling-wave coordinate ξ. F (ξ) satisfies a first-order NLODE called auxiliary equation (A.E.). By repeated substitutions and using the chain rules arise from A.E., the order of NLPDEs is reduced to first-order NLODEs. As a result, solutions can be easily deduced from known solutions of first-order NLODEs. In the modified F-expansion method [24] (in short mF-expansion) which is very similar to the F-expansion method, a different kind of A.E. is taken. In our proposed algorithms for F-expansion and modified F-expansion methods, one important advantage is that one can set A.E. in the general forms (13), (29) while in [19,24] the A.E. have been taken in the specific forms (C-6), (A-1) respectively. Unlike the function-expansion method, a huge variant of solutions can be obtained using different auxiliary equations. Even the solutions obtained from the function-expansion method can also be recovered from the F-expansion method. On the other hand, the first integral method (FIM) is based on the ring theory of commutative algebra. In FIM, a second-order NLPDE is automatically converted to the first-order NLODE when possible, which is known as the first integral of the original second-order NLPDE. Such automatic conversion is not possible in F-expansion and modified F-expansion methods. Then the solutions of second-order NLPDE are determined from known solutions of the first-order NLODE. This method has been successfully applied to a number of NLPDEs by researchers [25][26][27][28][29][30][31][32][33]. Keeping in mind the popularity of the F-expansion, mF-expansion, and first integral methods in the scientific community in finding the exact solutions of NLODEs and NLPDEs, we have been motivated to introduce the algorithms of these methods that can be easily implemented in a computer package for automatic derivations of closed-form traveling-wave solutions by employing these methods. In our proposed algorithms, we have generalized these methods and introduced some new features in the methods so that one can apply the new algorithms to different kinds of NLODEs and NLPDEs arising in various fields of science.
Searching for traveling-wave solutions of NLPDEs by hand is very tedious, cumbersome, and sometimes takes a long time. However, many modern Computer Algebra Systems (CAS), such as Maple, Mathematica, Reduce, Maxima, etc., help us to perform complicated and tedious algebraic calculations on a computer to find exact solutions to such problems accurately in a short time. In order to solve the NLPDEs, many computer packages are available. In 1996, Parkes and Duffy [34] have implemented tanh-expansion in their Mathematica package ATFM. Later complete implementation of tanh-expansion has been done by Li and Liu (2002) [35] designing the Maple package RATH. Baldwin et al. [36,37] have developed Mathematica package PDESpecialSolutions.m which admits polynomial solutions in tanh, sech, combinations thereof, JacobiSN, JacobiCN. RAEEM [38] is one of the most popular packages written in the Maple programming language, which is a comprehensive and complete implementation of some powerful methods such as the tanh-method, the extended tanh-method, the Jacobi elliptic function method, and the elliptic equation method. All these computer packages have implemented the function-expansion methods. One serious drawback of the function-expansion method is that the solutions which contain functions other than some specific type of functions, such as tanh, sech, JacobiSN, JacobiCN, etc. are not obtained. Additionally, the non-polynomial forms of these specific functions are not obtained also. We note that most packages have been developed within commercially available software frameworks, such as Maple and Mathematica. Besides this, to the best of our knowledge, the computer packages with the implementation of F-expansion and first integral methods are not available so far.
Keeping in mind all the above points of view, we have been motivated to develop a computer package or C++ library called GiNaCDE 1 [39], which implements our proposed algorithms of the high-performance automated F-expansion, modified F-expansion, and first integral methods for the first time. The symbolic manipulations of GiNaCDE depend only on GiNaC [40]. GiNaCDE extends the capability of the GiNaC library towards the solutions of differential equations. There are several advantages to use GiNaC over other CAS. GiNaC is a pure C++ library. It can accept C++ programming language, a general-purpose object-oriented programming (OOP) language, and it is fast like commercially available computer algebra systems. The source code is freely available as this library is written under the terms and conditions of the GNU general public license (GPL). One important property of GiNaC that differentiates it from other computer algebra programs we may have used: GiNaC assigns a unique (hidden) serial number for each newly created symbol object, and GiNaC uses this unique serial number instead of its name for algebraic manipulations. The serial number for the same name of a symbol may be changed in each running session of the GiNaC program. As a result, the symbols in the same algebraic expressions may be ordered differently during each running session of the GiNaC program. This happens because to order the symbols of an algebraic expression GiNaC internally uses a comparison predicate, called ex is less which uses an internal symbol id counter. There is also have another problem with using GiNaC as CAS. One should write a C++ program every time to solve each NLPDE, and compilation is also required for each one. However, many GUI tools like GTK, FLTK, Qt are available to make GUI applications with GiNaCDE, which facilitates users to solve NLPDEs automatically without writing programming and compilation each time. In this context, we have developed a GUI application called GiNaCDE GUI within the GTK framework.
We have noted that in the Maple package RAEEM [38], similar to the F-expansion method, an A.E. has also been considered, whose algebraic form is very similar to the A.E. as taken in the F-expansion method. However, in RAEEM, A.E. is handled with a maximum of four parameters, but in the case of GiNaCDE, one can choose the number of parameters higher than four that can be controlled from the user end. In order to get the exact traveling wave solutions including polynomial, exponential, triangular, hyperbolic, rational, Jacobi elliptic, Weierstrass elliptic type, RAEEM automatically chooses some predefined constant values of the four parameters of A.E., and these values cannot be controlled from the user end. On the other hand, in GiNaCDE software, one can assign any numerical or algebraic values to the parameters of A.E. and choose the upper limit on the number of parameters in A.E. in own choices. All this flexibility in handling the parameters of A.E. from the user end makes the GiNaCDE software more powerful and efficient to investigate more new types of closed-form solutions of NLPDEs.
The paper is organized as follows: In Secs. II, III, IV, we have simply described the algorithms for automated F-expansion, modified F-expansion, and first integral methods respectively and described each step of the algorithms with a working example of NLPDE. Key features of the implementation of the proposed algorithms in GiNaCDE software are discussed in Sec. V. In Sec. VI, to illustrate the efficiency and effectiveness of the GiNaCDE library, we have applied the library to a wider class of NLPDEs and found the exact solutions. Our results are also compared with an existing Maple package RAEEM. Finally, in Sec. VII, we discuss our results, and some conclusions have been made on our works. The exact analytical solutions of some auxiliary equations are given in the Appendix: A, B, C, D.

II. ALGORITHM OF AUTOMATED F-EXPANSION METHOD
This section explains the algorithms for the F-expansion method [19] to obtain closed-form traveling-wave solutions of NLPDEs automatically. However, in the proposed algorithms, we have to guide the solution process initially by providing some initial data. In that sense, the algorithms are not fully automated, but this makes the algorithms more powerful, and we can apply these algorithms to a huge variant of NLPDEs with different types of initial data. We divide the algorithm into five main steps (labeled F1-F5).
Let us consider an NLPDE with independent variables t, x 1 , x 2 , . . . x m = X and dependent variable u in the following general form where α i (i = 1, . . . , l) are the parameters, u = u(t, x 1 , x 2 , . . . x m ) and F is a polynomial about u and its derivatives.
Equation (1) (1), and sometimes the code may failure to give solutions due to the complexity in the problems. Here, it should be noted that our algorithms are also applicable to Eq. (1) without the parameters α i . Equation (1) is called input-NLPDE throughout this paper.
Step F1(Transform the NLPDE into a NLODE): At first we take a transformation in Eq. (1) where Here ξ is called traveling-wave coordinate and U (ξ) is the traveling-wave part of the solutions. The second part of Eq. (2) is e Iθ which is called the phase part of the solutions and θ is called phase coordinate. Usually, this part is present when differential equation (1) has an imaginary part. In our program library we must retain this phase part when Eq. (1) is complex otherwise we only retain traveling-wave part. Taking proper forms or values of the constant coefficients (k i (i = 0 . . . m)) of traveling-wave coordinate ξ and phase angle constants p i (i = 0 . . . m), we can transform the NLPDE (1) into the nonlinear ordinary differential equation (NLODE) with single independent variable ξ and dependent variable U (ξ). Repeatedly applying the chain rule on Eq. (2), Eq. (1) is transformed into NLODE. The form of NLODE in general form is given by where U (n) (n = 1, 2 . . .) indicates n times differentiation with respect to ξ. All the three algorithms (F-expansion, mF-expansion, FIM) work only when Eq. (5) is a polynomial in variable U (ξ) and its derivatives, and does not explicitly depend on the independent variables ξ. Then our algorithm check the integrability of Eq. (5). If Eq. (5) is integrable, our method try to integrate Eq. (5) and if integration is successful then one can assign a value to each integrating constant (ic i , i = 1, 2, . . . , η, η is the number of integration) in own choices. In the case of complex NLPDE, if real and imaginary both parts are present in the transformed NLPDE, one part is taken to solve. To select one part from complex NLODE, following strategies are followed i. Assume all the parameters (α i , k i , p i ) are real.
iii. Check in which part (Re or Im) dependent variable U (ξ) is not present, but minimum one number of parameter is present there. Get the constraint on that parameters contained in this part. Take other part (may be Re or Im) as NLODE whose solutions are to be solved.
Suppose in the real part Re, the dependent variable U (ξ) is not present and there is only present the parameters α i , k i , p i . Then Eq. (6) can be expressed by and we have to solve the NLODE Im α i , k i , p i , U, U (1) , U (2) . . . = 0 subject to the constraint Re (α i , k i , p i ) = 0.
iv. If the above step is failure, compare the expressions of Re and Im to check whether they are the same NLODE for some constraint on the parameters of Re or Im. If they are the same NLODE, take anyone part as NLODE whose solutions are to be solved.
Suppose for the constraint we get from Eq. (6), then the solutions can be determined either from the NLODE Re or Im subject to the constraint relation (8).
If any above criteria are not satisfied, the complex input-NLPDE equation cannot be transformed into single NLODE in our algorithm and the method does not work successfully.
Step F2(Determine the highest power N of finite power series): Now according to the F-expansion method, the solution of the NLODE (5) can be expressed as the finite power series, which is where a i (i = 0 . . . N ) and b i (i = 1 . . . N ) are constants to be determined later. In Eq. (10) first term is positive part and second term is negative part in the solution. The value of N (positive integer or positive noninteger) in Eq. (10) can be determined by considering a homogeneous balance between the highest order nonlinear term with the highest order derivative of U (ξ) in Eq. (5). However, to automate this process, we have employed a method described in [35]. At first, Eq. (5) is expanded in the sum of product (SOP) form. Our aim is to determine the highest possible value of N , and so it is sufficient to replace U with U N . Assuming degree of U (ξ) is D[U (ξ)] = N , we replace U by U N and collect the degrees of each term appearing in Eq. (5) by a variable, say E. To determine the degree of an expression, we use the relations Generally in (max( E),N ) plane we get a turning point and the value of N is taken at this point. Whole procedure is automated in the following order: i. Substitute U by U N and simplify.
ii. Expand and express in SOP form.
iii. Collect the degree of each term appearing in SOP form and store them in a list E.
iv. Replace N in E by a sequence of numbers whose first number is 0. In GiNaCDE, we take three number sequence with common differences-1/2, 1/3, 1/4 up to last number 11.
v. For each number in sequence, calculate max(E).
vi. Calculate differences between the value of max(E) for successive numbers in the sequence.
vii. To get highest power N , take the number in sequence for which differences are not same with previous one. More clearly if i is the current number in the the sequence and (max( It is clear that for a larger value of N the complexity in mathematical operations is increased. To avoid such complexity in derivations, one should set the maximum allowed value of N . However, sometimes in some cases, the auto-evaluation of N may fail. Step F3(Derive the system of Nonlinear Algebraic Equations for the coefficients of F (ξ)): F (ξ) satisfy the firstorder nonlinear ODE (also called auxiliary equation (A.E.)) where F (F (ξ)) is some known functions of F (ξ). The prime over F (ξ) represents differentiation with respect to ξ. In case of F-expansion method where δ is a positive integer and A i (i = 0, 1, . . . δ) are coefficients of A.E. Here one can use any functional form of Eq. (13) by choosing any positive integer value of δ and real values of A i . As a result the solutions of first-order NLODE (13) can be expressed in terms of a large variety of functions such as polynomial, exponential, trigonometric, hyperbolic, rational, Jacobi elliptic etc.. However, taking some well-known functional forms of F in Eq. (13), we have shown some solutions of F (ξ) in the Appendix: C, D. The higher derivatives of F (ξ) using Eq. (12) can be expressed by Now substituting Eq. (10) into Eq. (5) and using Eq. (14) with (13) we get an expression, and the numerator of resulted expression contains F (ξ) j F k (j = 0, 1, 2, . . . ; k = 0, 1) terms. Setting each coefficient of F (ξ) j F k to zero an  (5) is integrable, are present.
The nonlinear system of equations is a set of simultaneous equations in which the unknowns (the constant parameters) appear as variables of a polynomial of degree one or higher than one. Suppose the system of algebraic equations is solved for all parameters that are present in the system of algebraic equations. In that case, it takes a larger time to get solutions for the system of equations, even sometimes solutions are not obtained for a complicated system. To reduce the calculating time and the complexity in derivations, we categorize all the parameters into two different types. They are external parameters and internal parameters. External parameters are α i (i = 1 . . . l), A i (i = 0, 1, . . . δ) which are present in input-NLPDE equation and auxiliary equation respectively. When input-NLPDE equations are integrable, the generated integrating constant(s) ic i (i = 1 . . . η) is also external parameter. On the other hand all the remaining parameters, such as . . m) are called internal parameters as they generate internally. Nonlinear algebraic systems are always solved for internal parameters. But one can have control over external parameters to choose the parameters for which nonlinear algebraic system is to be solved. For this purposes, the programming variables ASolve and paraInDiffSolve (detailed descriptions are given in the software user manual) are used to choose the unknowns from the external parameters in our choice, and it will reduce the calculating time and can handle more complicated algebraic expressions. At the same time, the exact solutions are determined subject to the conditions on chosen external parameters.
Step F4(Solve the system of Nonlinear Algebraic Equations): Analysing and solving the nonlinear algebraic system is a vital and challenging step among all steps of the method. In fact, the number of exact solutions of NLPDE derived by F-expansion, modified F-expansion, or first integral methods is entirely depending on how many solutions are obtained from a nonlinear algebraic system. The executing time of the software mainly depends on this step. Many methods are available to solve nonlinear algebraic systems such that Gröbner basis methods [41], the Ritt-Wu characteristic sets method implemented by Wang [42,43], and the Reduced Involutive Form (Rif) code by Wittkopf and Reid [44]. D. Baldwin et al. [45] have employed a simple algorithm to design a powerful nonlinear solver in Mathematica. We have followed their algorithm to create a nonlinear solver in GiNaC symbolic system. The nonlinear solver implemented in [45] solves the entire system in an automated way using the built-in Mathematica function Reduce. Their solver can solve polynomial and non-polynomial systems both. The nonlinear solver implemented by us can solve only the polynomial system required in this application, and its own C++ function solves the polynomial equations.
The steps used in this algorithm are very much like the steps used to solve a nonlinear algebraic system by hand. In this method, the simplest equation is solved for sorted unknown parameters. Then the solutions are substituted in the remaining equations. Such solving and substitution procedures are repeated until the system is completely solved. We operate the whole procedure in the following order: i. Check whether each equation is polynomial in unknowns.
ii. Factor and simplify each equation.
iii. Measure complexity of each equation by the number of add containers, unknown parameters, and the degree of unknowns. Then, sort the system based on their complexity. If more than one equations have the same complexity, in GiNaCDE they are sorted according to the GiNaC in-built comparison predicate ex is less.
iv. Sort the unknown parameters contained in the simplest equation by their degree.
v. Solve the simplest equation for the lowest degree unknown. If the number of unknown for the lowest degree is greater than one, then GiNaCDE uses the comparison predicate ex is less to choose the unknown. If solutions are absent, solve the unknown for the next higher degree.
vi. Substitute the solutions into the remaining equations and simplify.
vii. Repeat the steps i-vi until all the equations are reduced to zero.
viii. Substitute all the unknowns which are present in the computed solutions with the help of other solutions.
ix. Test the solutions by substituting them into each equation.
x. Finally, collect all solutions branches.
Our solver is powerful and can easily handle nonlinear equations (of course, polynomials in unknowns) with multiparameters. Sometimes, there are risks of missing some solutions due to numerous parameters in the system or if the system is high degree. In this solver, the unknowns from all parameters appearing in the system are chosen in order of complexity. Then the solutions for these unknowns are expressed in terms of other parameters that are to be regarded as arbitrary parameters. Sometimes it is observed that the solutions become simpler where these arbitrary parameters are taken as unknowns.
Step F5(Build solutions with calculating steps): Substitute the solutions obtained in step F4 into Eq. (13) and obtain the solutions of F using the Appendix: C, D. Then, to obtain traveling-wave solutions of Eq. (5), substitute F and the solutions obtained in step F4 into Eq. (10). Finally the explicit solutions in original variables are obtained using Eqs. (2), (3).

A. Explanation of each step with a working example
To illustrate each algorithm steps for F-expansion method described above, we take the one dimensional cubic nonlinear Schrödinger (NLS) equation [46] Iu t − pu xx + q|u| as an example. In Eq. (15) p, q are non-zero real constants and u(x, t) is a complex-valued function depends on the variables x, t.
Step F1: Equation (15) has an imaginary and a real part. So retaining the phase part in Eq. (2) and making traveling-wave transformation Eq. (4) yields Substituting Eqs. (17) into Eq. (15) we get Note that, Eq. (18) has real and imaginary both parts. From imaginary part we get the condition and from real part we obtain the NLODE Step F2: Substituting U → U N , degree of each term in Eq. (20) is collected by a list variable E: 4.5 1.5 2 6 1. 5   TABLE I: Here a number sequence with common difference 1/2 has been used to determine the highest power N . Bold number denotes the turning point and the corresponding value N = 1 is balancing highest power.
Then using balancing highest power N = 1, finite power series is expressed by Step F3: We choose A.E.
Step F4: Resulting nonlinear algebraic equations are solved for internal parameters only and the solutions are In the above solutions the parameters p, q, k 0 , p 0 , A 0 , A 2 are arbitrary.
Step F5: The solutions for F of Eq. (23) are derived with the help of Eq. (D-13) in Appendix: D. Then combining (22) with (19) and substituting (26a), we obtain the closed-form solutions to (15). For brevity the solutions are not shown here. Here we present only one solution which illustrates our algorithm: where θ = x and C is an arbitrary constant. There are two sign combinations in each equation and total 12 solutions are obtained. Similarly, using (26b), we can derive an another exact solution

III. ALGORITHM OF AUTOMATED MODIFIED F-EXPANSION METHOD
The algorithm for the automated modified F-expansion method [24] is very similar to the algorithm of the automated F-expansion method. This algorithm also has five main steps (labeled MF1-MF5), and it has only one difference to the F-expansion method. The difference is that a different form of A.E. is taken in step MF3 in comparison to step F3. Therefore, one can check new exact solutions of the NLPDE applying both methods (F-expansion method and modified F-expansion method) to the same NLPDE with different forms of A.E. Details of all steps are described below: Step MF1(Transform the NLPDE into an NLODE): Same as step F1.
Step MF2(Determine the highest power N of finite power series): Same as step F2.
Step MF3(Derive the system of Nonlinear Algebraic Equations for the coefficients of (F (ξ)): In the modified F-expansion method, the solution of the NLODE (5) is also expressed by a finite series like Eq. (10). In this method, we have generalized the modified F-expansion method [24] taking the A.E. in more general form where δ is a positive integer and A i (i = 0, 1, . . . δ) are coefficients of A.E. One can choose any functional form of Eq.  29) we get an expression appearing the terms F (ξ) j (j = 0, 1, 2, . . .) in the numerator. The equations must vanish identically. Hence, to generate a nonlinear algebraic system, equate to zero the coefficients of the power terms in F .
Step MF4(Solve the system of Nonlinear Algebraic Equation): Similar strategy as in step F4.
Step MF5(Build solutions with calculating steps): Substitute the solutions of step MF4 into Eq. (29). Obtain the solutions of F using Appendix: A, B. Then, substitute F along with the solutions of step MF4 into Eq. (10). To get the explicit solutions in original variables, Eqs. (2), (3) are used.
As the modified expansion method is very much similar to the F-expansion method (the difference with the Fexpansion method is that a different form of A.E. is taken here), we do not explain each step of this method with a working example.

IV. ALGORITHM OF AUTOMATED FIRST INTEGRAL METHOD
In first integral method [14], one important advantage over F-expansion and modified F-expansion methods is that one does not have to choose A.E. to solve NLPDEs; instead, the input-NLPDE is automatically reduced to a suitable first-order NLODE whose solutions have to be calculated.
The algorithm for the automated first integral method has eight main steps (labeled FIM1-FIM8). Now we give an outline of every step as follows: Step FIM1(Transform the NLPDEs into NLODEs): Same as step F1. The condition for applying first integral method to the Eq. (5) is that Eq. (5) must be a second-order NLODE. Therefore Eq. (5) is expressed in the form Step FIM2(Convert into a system of NLODEs): We assume that U (ξ) = X(ξ) and introducing a new independent variable Y (ξ) = X ξ (ξ), Eq. (30) can be rewritten as a system of NLODEs [14] X is the coefficient of the highest derivative term in Eq. (30).
Step FIM3(Apply Division Theorem): If X(ξ), Y (ξ) are nontrivial solutions of Eq. (31), then applying the Division Theorem [14] there exist an irreducible polynomial in the complex domain C[X, Y ] such that where a i (i = 0 . . . N ) are polynomials of X and a N = 0. Equation (32) is called the first integral to Eqs. (31a) and (31b). Using Division Theorem there exists a polynomial (g(X) + h(X)Y ) such that Using Eqs. (31), (32) in the Eq. (33), we get Dot over a i (X) denotes derivative with respect to X.
Step FIM4(Derive the Algebraic System of equations for coefficients of Y i ): Comparing coefficients of Y i (i =  N + 1, N, . . . , 1, 0) on both sides of (35), and for H(X) = 0 canceling H(X) in denominator from both sides we obtain where a i (X) = 0 for i < 0 and i > N .
In the next step FIM5 we take a N = 1 to derive the polynomial forms of h(X), g(X), a i (i = 0.. . . . N − 1). For a N = 1 from Eq. (36a) we obtain If H(X) is not a constant and at the same time degree of P (X(ξ), Y (ξ)) in variable Y is 2 then it is clear from Eqs. (36a), (37) that h(X) will not be polynomial in X. In this case we avoid such non-polynomial form of h(X) by making the transformation [47] dξ = H(X)dη, in Eq. (31) temporarily. Applying the transformation in Eq. (31) we get Consequently the nonlinear algebraic system becomes Y N : H(X)ȧ N −1 (X) + N K 1 (X)a N (X) + (N + 1)K 0 (X)a N +1 (X) = g(X)a N (X) + h(X)a N −1 (X), (40b) . . . . . .
Now for a N = 1, from Eq. (40a) we get h(X) = N K 2 (X) which is polynomial in X. In the following steps, we explain all the procedures with the help of Eq. (36), because the same procedures are applicable when Eq. (40) is considered for the non-polynomial case of h(X).
Step FIM6(Derive the Nonlinear Algebraic System for the parameters g i , a ij ): Substitute a N = 1 and (41) into each equation of (36). The coefficients of the power in X, Y in each equation of (36) must vanish. Collect the coefficients and generate a nonlinear algebraic system of equations parametrized by g i , a ij , K, P , integrating constants (for integrable NLPDE) and parameters appearing in input-NLPDE.
Step FIM7(Solve the Nonlinear Parameterized Algebraic System): Here, the external parameters are parameters appearing in input-NLPDE and integrating constants. Internal parameters are g i , a ij , K, P . The nonlinear algebraic system is solved following a similar process in step F4. Like step F4, here also, the runtime of this algorithm mainly depends on this step.
Step To illustrate every step of the above algorithm for the first integral method, like the F-expansion method, we have also used the NLS equation (15).
Step FIM6: Substituting a 1 = 1 and Eq. (41) into each equation of Eqs. (45) and collecting coefficients of X i Y j (i = 0, 1, 2, 3, j = 0, 1), we obtain the following set of nonlinear algebraic equations: Step FIM7: Here, the internal parameters are g 0 , g 1 , a 00 , a 01 , a 02 , k 0 , k 1 , p 0 , p 1 and external parameters are p, q. Solving the set of equations (49) for the internal parameters only, the solutions are . (50b) Step FIM8: Substituting (50a) into (43) the first integral form to Eq. (20) is Equation (51) is Riccati equation and to find its solutions Appendix: A is used. Finally, using (19) and substituting (50a), 5 exact solutions of Eq. (15) can be obtained. One can note that when g 1 = −2q p , we have the same first integral form as given in [32].

V. ALGORITHM IMPLEMENTATION
We have implemented the proposed algorithms described in Secs. II, III and IV, into GiNaCDE. GiNaCDE is a C++ library that is built on a pure C++ symbolic library GiNaC [40]. Besides this library version, we have also developed a GUI version of GiNaCDE called GiNaCDE GUI. When we solve differential equations using GiNaCDE GUI, we do not have to write any C++ code, and compilation of any code is not required. This GUI version guides us in each step to obtain the output results. However, a complete guide of GiNaCDE GUI has been provided with the GiNaCDE source code. In both versions of GiNaCDE (library and GUI), the output results are saved in a text file with calculating steps. Output results can be saved in Maple or Mathematica programming language by assigning the C++ macros Maple or Mathematica to the programming variable output.
In order to start a solution process for a given NLPDE or NLODE in GiNaCDE, we require some initial data (the options and parameters specified by the user). Some programming variables such as twcPhase, positivePart, negativePart, NValue, degAcoeff, ASolve, paraInDiffSolve (detailed descriptions of these variables are given in the software user manual) are available, which can be initially set up by the users in their own choices for getting better results for a given NLPDE before starting the solution process in the GiNaCDE software. In all three methods, if the input NLPDE or NLODE is complex, then the software tries to separate the real and imaginary parts following the step F1. In all three methods, if the input NLPDE or NLODE is integrable, the software tries to integrate them after starting the solution process. If the integration is successful, then the software gives an option to us to assign a numerical or symbolic value to the integrating constant(s) ic i (i = 1, 2, ..) in our choices. All these options make the software more powerful and flexible, enhancing its ability to find many new exact solutions to huge variants of NLPDEs. Now we shall discuss some implementation details of each method separately: F-expansion and mF-expansion methods: These methods have been implemented in F expns methd.cpp and F expns methd.h files. The F-expansion and mF-expansion methods are chosen by the C++ macros F expansion, mF expansion respectively. One can use the coordinates K, P in own choices with the help of the programming variable twcPhase. We can take any one or both parts in the solutions (10) with the help of the programming variables positivePart, negativePart. In GiNaCDE, we set the maximum allowed value of N at 10. However, sometimes in some cases, the auto-evaluated value of N exceeds 10. Then this step fails to find N . In this case, we can check the solutions of input-NLPDE by specifying the value of N in our choice lower than 10 with the help of the programming variable NValue. If we do not assign any value to the variable NValue, the value of N is auto-evaluated following the criteria in step F2. In these methods, we initially have to input A.E. to start the solution process in GiNaCDE manually. One can choose the parameters A i (i = 0, 1, . . . δ) in own choices with the help of the programming variable degAcoeff. The variable ASolve confirms whether the nonlinear algebraic system will be solved for the parameters contained in A.E. ( i.e., the parameters A i where i = 0, 1, . . . δ) along with other parameters. The parameters appearing in input-NLPDE (these parameters are belong to external parameters) are supplied in the programming variable paraInDiffSolve to solve the nonlinear algebraic system for those parameters also. This will determine the conditions on that external parameters so that exact solutions are obtained.
FIM: This method has been implemented in fim.cpp and fim.h files. This method is chosen by the C++ macro FIM. For the first integral method, we do not have to input A.E. Here, for initializations, we have only three programming variables twcPhase, NValue, paraInDiffSolve. The variables twcPhase, paraInDiffSolve have been discussed above. Like F-expansion method the value of N is also assigned by the variable NValue and the default value is N = 1. However, in our library, the allowed values of NValue are 1 and 2.

VI. APPLICATION
In this paper, we have presented the algorithms that automate three different methods which can be applied to determine exact solutions of the wide variant of NLPDEs. F-expansion and mF-expansion methods can be applied to higher-order NLPDEs. But, FIM method is applicable to an NLPDE when its transformed NLODE (30) is secondorder only. Actually, there are no rules to know in advance the appropriate method among three to solve the given NLPDE. To illustrate the effectiveness of the algorithms, in the following subsections, using the GiNaCDE library, we have applied all these three methods to solve some NLPDEs, including higher-order and complex NLPDEs.
To start a solution process for a given NLPDE, we require some initializations (the options and parameters specified by the user). There is no rule to know in advance a specific type of initialization to get better results. However, we have searched some initializations for solving the following five examples so that the manually chosen initializations give better results. In the following examples, such initializations for each NLPDE are also described. The chosen assigned values of the integrating constant(s) of the integrable NLPDEs among the following examples are also described. In the case of F-expansion and mF-expansion methods, there are no rules to know in advance what type of A.E. can give exact solutions for a given NLPDE. In this context, we have shown the solutions of some well-known A.E. in the Appendix: A, B, C, D. We have checked the solutions for the NLPDEs listed below using some specific types of A.E. These specific forms of A.E. and the corresponding number of solutions with calculating times are shown in Table II.
Here one can note that more new types of solutions of the NLPDEs can be obtained using different forms of A.E.
We have well described the procedure to get the exact solutions of an NLPDE using GiNaCDE in the user manual. We have solved the following examples using the C++ code as given in the examples folder provided with the GiNaCDE source code. However, in this paper, we do not quote the algebraic expressions of the exact solutions obtained in the following examples since their final expressions are too lengthy. We have observed that some of the solutions calculated by us were previously reported in the literature. Interestingly in some cases, new solutions are also obtained, and some new solutions are shown here.
which belongs to a different class of integrable NLPDEs, often termed C-integrable equations, i.e. linearizable through a change of dependent variables, into the linear, free Schrödinger equation [48]. In seeking exact solutions of Eq. (52), many powerful methods have been applied, such as the first-integral method [49], G ′ /G-expansion method [50]. In the all methods (F-expansion, mF-expansion and FIM), after substituting the traveling-wave solution (16) into Eq. (52), it is transformed into NLODE. Then this NLODE is separated into real and imaginary parts. The imaginary part is integrable and after one time integration GiNaCDE itself assigns the integrating constant to 0. Then from imaginary part, following the step F1 in the algorithm of automated F-expansion method, GiNaCDE derives the condition In F-expansion and mF-expansion methods, value of N is 1 2 which is auto-evaluated following the step F2 in our proposed algorithm. With the condition (53), solutions of the real part are evaluated in all three methods. In Fexpansion method, we do not get any solution. In mF-expansion method, using the solutions of Bernouli equation given in (B-2) with condition (53), the solutions to Eq. (52) are obtained.
B. Seventh-order Sawada-Kotara equations J. Feng [51] considered the following Seventh-order Sawada-Kotera equation (sSK) to find out some closed-form trigonometric, hyperbolic, rational solutions applying G ′ /G expansion method. Several methods such as Adomian decomposition method (ADM) [52], He's variational iteration method [53], Reconstruction of Variational Iteration Method (RVIM) [54], has been applied for computing approximated solutions to Eq. (54). The initializations we have used to solve the Eq. (54) by employing each method available in the GiNaCDE library are given below: In the all methods, after transforming Eq. (54) into NLODE using the traveling-wave coordinate ξ given by GiNaCDE integrates this NLODE one time, and we have assigned the integrating constant ic 1 to 0 in our choice. In FIM, no solution is obtained. In the case of the F-expansion method, the library returns some new Jacobian solutions using the known solutions of A.E. given in (C-7). Total 40 new solutions are obtained. For simplicity only 3 new solutions are presented here: where In the above solutions two sign combinations in each solution are to be taken.

C. Fifth-order Generalized KdV equation
Consider the general fifth-order KdV equation (gKdV) [55][56][57][58][59] where p, q, r are arbitrary non-zero and real parameters. The model equation (62) has a wide range of applications in many important physical phenomenon including quantum mechanics, nonlinear optics, plasma physics, one dimensional nonlinear lattice. The values of the parameters p, q, r will drastically change the characteristics of the gKdV equation (62). Many well-known equations can be constructed from the gKdV equation by changing these parameters, such as: the Lax equation q = 2p, r = 3p 2 10 , the Sawada-Kotera (SK) equation q = p, r = p 2 5 , the generalized Kaup-Kupershmidt equation (GKK) q = 5p 2 , r = p 2 5 , the generalized Ito equation (GI) q = 2p, r = 2p 2

9
. A variety of powerful and direct methods have been applied to obtain exact solutions to Eq. (62). Among of them are: the extended tanh method [55], the generalized tanh-coth method [56], an extended Jacobian elliptic function expansion approach [57], Exp-Function method [58,59]. Zhi-bin Li et al. in their works [35] have found out some new conditions among p, q, r for which some new solutions were obtained using their computer package RATH. We ask GiNaCDE to derive solutions of Eq. (62) by applying all three methods available in the library.
We have used the following initializations: In the all methods, GiNaCDE transforms Eq. (62) into traveling-wave coordinate ξ using the Eq. (55) and integrates it one time. We assign the numerical value 0 to integrating constant ic 1 . Interestingly, in the case of F-expansion method, it finds out some new exact solutions with a new condition which is completely different from a list of conditions among p, q, r derived in [35]. For the condition (63), with the help of solutions (C-9) we obtain 10 new Jacobian elliptic solutions. For brevity we quote only following 3 new solutions where

D. Perturbed NLS Equation with Kerr Law Nonlinearity
We consider the perturbed NLS equation with Kerr law nonlinearity [29,[63][64][65][66] where u(t, x) represents the complex function and the parameters G 1 , G 2 and G 3 are the higher order dispersion coefficient, the coefficient of Raman scattering, the coefficient of nonlinear dispersion term respectively, while A represents fiber loss. The model equation (68) has important application in various fields, such as semiconductor materials, optical fiber communications, plasma physics, fluid and solid mechanics. Several methods for finding the exact solutions to (68) have been applied, such as first integral method [29], the improved tan φ(ξ) 2 -expansion method [64], the modified trigonometric function series method [65], the modified mapping method and the extended mapping method [66]. G. Akram et al. [63] recently have successfully applied the extended G ′ /G 2 -expansion method and the first integral method on Eq. (68) to find some new exact solutions, which include hyperbolic function solutions, trigonometric function solutions, rational function solutions and soliton solutions. We run GiNaCDE software applying all the three available methods on Eq. (68). Here we have used the following initializations: In all methods, the software substitutes the traveling-wave solution (16) into Eq. (68), and separates the real part and the imaginary part following the step F1. The imaginary part is integrated one time and the software assigns the constant of integration to zero. Then the algebraic expressions of real part and imaginary part are compared and the software detects that they are same equations subject to the following conditions Same conditions were obtained in [66]. Therefore GiNaCDE evaluates the exact analytical solutions of imaginary part only.
In the case of F-expansion method, using the solutions (D-8) GiNaCDE produces two new traveling-wave solutions given by where

E. Kudryashov-Sinelshchikov Equation
Now, we study the following Kudryashov-Sinelshchikov equation proposed in [67,68]: where g, n, k, d and e are real parameters. Equation (73) models the pressure waves in a liquid and gas bubbles mixture when the viscosity of liquid and the heat transfer are both considered. In [47], authors have found the exact traveling-wave solutions of Eq. (73) by using the first integral method. We have also employed the first integral method to find exact traveling-wave solutions of Eq. (73) with the help of GiNaCDE and compare our result with the result in [47].
Finally, we have also solved the above-listed five examples with the help of RAEEM Maple package [38]. We have seen that RAEEM with Maple v. 8.0 is unable to solve the NLPDEs (52), (62) and (68). RAEEM can only solve the NLPDEs (54), (73), and gives 6, 6 solutions with CPU times 61 s, 250 s respectively. Beside the above examples, we have successfully solved more than 20 NLPDEs (provided in the result and test folders of GiNaCDE source code) using the GiNaCDE library and it gives results not more than 100 S.

VII. CONCLUSIONS
We have presented the algorithms for the high-performance automated F-expansion and First Integral Methods. We have also implemented these algorithms in a C++ library named GiNaCDE. This library is used to find the closed-form traveling-wave solution of the NLPDEs of the form (1). The solution methods described in the Secs. II, III and IV, are very tedious if we apply these three methods on an NLPDE by hand. The program library automates the methods and delivers possible solutions when they exist. However, in this context, in order to start the solution process, we have to provide some initial data (the options and parameters specified by the user) to the library. After running the solution process, there is also scope to assign the values to integrating constants generated from integrable NLPDEs. From these points of view, we can tell that our algorithms are not fully automated. But, these features make the library more efficient and powerful to tackle a large class of NLPDEs. Due to the implementations of three different methods in one software, one can easily check exact traveling-wave solutions of a huge variant of NLPDEs by applying those methods in one place with less labor.
We have introduced an exciting feature in our proposed algorithms, by which one can solve the complex NLPDEs. The program library can integrate the input-NLPDE or the transformed NLODE when possible. The generated integrating constants can be assigned with the values in our choice after running the program. Following [47], we have added another interesting feature in the algorithm for the first integral method. According to this new feature, our proposed algorithm is capable of tackling the non-polynomial form of h(X) in Eq. (37). We have applied the package to a wide class of nonlinear evolution equations. It successfully recovered all previously known solutions that many powerful methods had found. More importantly, we have found new solutions and a more general form of solutions for some of the equations considered.
Many computer packages such as RATH [35], PDESpecialSolutions.m [36,37], RAEEM [38] are currently available which can solve the NLPDEs of the form (1) using some popular methods other than F-expansion, modified Fexpansion, and first integral methods. These currently available packages cannot solve the NLPDEs containing complex conjugate functions, but GiNaCDE can solve such type of NLPDEs. We have compared the performances of GiNaCDE and RAEEM [38] by solving some NLPDEs listed in Sec. VI. We have observed that RAEEM cannot solve 3 NLPDEs among 5 examples, but interestingly GiNaCDE can solve all the NLPDEs. We have also noted that GiNaCDE is much faster than RAEEM. We have tested the GiNaCDE library on more than 20 NLPDEs and it gives results not more than 100 S.
The program library shows the results with calculating steps, and the results are saved in a text file. We have also provided a GUI version of our library that can be used without any programming knowledge. Thus, the GiNaCDE can be easily used to solve a broad class of NLPDEs for obtaining exact solutions of NLPDEs, and it is very efficient and fast to get the solutions.
Since our software performs the computations automatically from start to finish without human intervention except for few steps in the algorithms, it may not return the solutions in the simplest form. If required, the user can further manipulate and graph the solutions. The most important and vital step in the algorithms is analyzing and solving the nonlinear algebraic system. In the algorithms, most of the computation time is spent in solving the nonlinear algebraic system. So for a complicated algebraic system, the library takes more time to solve an NLPDE. Furthermore, sometimes the nonlinear algebraic system may be quintic or higher degree and may contain a large number of parameters, and then it may be unsolvable in analytic form. Therefore, although unlikely, due to the limitations of the algebraic solver, some exact solutions of NLPDEs may be missed.
The GiNaCDE library for its algebraic manipulations depends only on the GiNaC library [40]. GiNaC algebra system assigns a symbol id to each name of a symbol, and unlike other algebra systems, it uses the symbol id instead of its name for algebraic manipulations. Since the symbols id may change for each running session of the program, the same algebraic expressions may appear in different order of symbols. We have noted that the present version of GiNaCDE cannot handle a system of NLPDEs. However, it is expected that these limitations can be removed in the future version of GiNaCDE. For F-expansion and modified F-expansion methods, the solutions were considered in the form (10). Many researchers have considered the solutions in several other forms. Z. Sheng [23] considers solutions of the form Here the parameters a 0 , a i , a −i and b i are not constant, and they depend on traveling-wave coordinate ξ. Y.M. Zhao [19] seeks the solutions of the form So there is a scope to extend the capabilities within GiNaCDE to find the solutions in the forms (74), (75). We will introduce such extensions in a future version of GiNaCDE.

Appendix A: Solutions of Riccati equation
In case of Riccati equation, Eq. (29) take the form The solutions of the equation (A-1) are [? ] where S = A 2 1 − 4A 0 A 2 and C is auxiliary constant. In the above solutions, from Eq. (A-2a) to (A-2e), the condition A 2 1 − 4A 0 A 2 > 0 must be satisfied. Please note that in the following solutions if not mentioned C has to be assumed as an auxiliary constant.
The solutions of the equation (B-1) are Solutions of Eq. (C-6) in terms of Jacobi elliptic functions are listed in below (with the conditions that all the algebraic expressions within square-root must be greater than 0): We note that Jacobian functions in all the above solutions are in the form C 1 JacobiXY(C 2 ξ, C 3 ), where C 1 , C 2 , C 3 are constants which contains only the constant parameters of differential equation (C-6). It will be better to express them in the convenient form JacobiXY(ξ, m) by putting C 1 = 1, C 2 = 1, C 3 = m. We derive the values of the constant parameters A 0 , A 2 , A 4 in terms of modulus m by solving the system of equations {C 1 = 1, C 2 = 1, C 3 = m} simultaneously. There have some advantages to have solutions in the form JacobiXY(ξ, m). We can transform JacobiSN(ξ, m), JacobiCN(ξ, m), JacobiDN(ξ, m) into hyperbolic functions tanh(ξ), sech(ξ), sech(ξ) respectively simply making the approximation m → 1. We can also simply get the trigonometric transformations JacobiSN(ξ, m) → sin(ξ), JacobiCN(ξ, m) → cos(ξ) and also get the transformation JacobiDN(ξ, m) → 1 in the limit case m → 0. So by choosing a single A.E. (C-6) one can obtain a variety of exact solutions of NLPDEs in terms of Jacobian elliptic functions, hyperbolic functions and trigonometric functions. Table III shows all the solutions of the parameters A 0 , A 2 , A 4 solving the system of equations {C 1 = 1, C 2 = 1, C 3 = m} and the corresponding solutions of Eq. (C-6) are given in last column. Interestingly the values of A 0 , A 2 , A 4 in Table III coincide with the ones given in [13].   [13].
Here the values of A 0 ,A 2 and A 4 are taken through modulus (m) in most cases.
The solutions of an another first-order NLODE can also be expressed in terms of Jacobi elliptic functions. The solutions of Eq. (C-8) are (with the conditions that all the algebraic expressions within square-root must be greater than 0) (C-10) In this section, we discuss the solutions of more different types of the first-order NLODE. Among them, we can choose an auxiliary equation suitable for input-NLPDE. Type-1 first-order NLODE: Let us first consider the first-order NLODE as follows (D-1) The above equation admits following special hyperbolic solutions [69]: 4A2 then it has a bell profile solution and a singular solution .

(D-2b)
Type-2 first-order NLODE: For A 0 = 0 Eq. (D-1) is reduced to 4A2 , Eq. (D-3) also admits a kink profile solution and a singular solution Type-3 first-order NLODE: An another type of first-order NLODE is which has different type solitary wave solutions [70] F .

(D-11a)
Type-4 first-order NLODE: We will now consider one more simplified first-order NLODE We find some exact solutions of Eq. (D-12) containing exponential, hyperbolic functions, which are listed below