Rational Biparameter Homotopy Perturbation Method and Laplace-Padé Coupled Version

The fact that most of the physical phenomena are modelled by nonlinear differential equations underlines the importance of having reliable methods for solving them. This work presents the rational biparameter homotopy perturbation method RBHPM as a novel tool with the potential to find approximate solutions for nonlinear differential equations. The method generates the solutions in the form of a quotient of two power series of different homotopy parameters. Besides, in order to improve accuracy, we propose the Laplace-Padé rational biparameter homotopy perturbationmethod LPRBHPM , when the solution is expressed as the quotient of two truncated power series. The usage of the method is illustrated with two case studies. On one side, a Ricatti nonlinear differential equation is solved and a comparison with the homotopy perturbation method HPM is presented. On the other side, a nonforced Van der Pol Oscillator is analysed and we compare results obtained with RBHPM, LPRBHPM, and HPM in order to conclude that the LPRBHPM and RBHPMmethods generate the most accurate approximated solutions.


Introduction
Solving nonlinear differential equations is an important issue in sciences because many physical phenomena are modelled using such classes of equations. One of the most powerful methods to approximately solve nonlinear differential equations is the homotopy perturbation method HPM 1-45 . The HPM method is based on the use of a power series, which transforms the original nonlinear differential equation into a series of linear differential equations. In this work, we propose a generalization of the aforementioned concept by using

2.5
For the HPM method 11-14 , we assume, without loss of generality, that the solution for 2.4 can be expressed as a power series of p v p 0 v 0 p 1 v 1 p 2 v 2 · · ·. 2.6 In the limit, when p → 1, the approximate solution for 2.1 is give as where v 0 , v 1 , v 2 , . . . are unknown functions to be determined by the HPM method. The series in 2.7 is convergent in most cases 1, 2, 11, 14, 36 . For the RBHPM method, the homotopy in 2.4 can be rewritten as Q ap 1 − a q, 0 < a < 1, 2.9 where p and q are the homotopy parameters, and a is a weight factor that affects them in complementary proportions. Now, we assume that the solution for 2.8 can be written as the quotient of the power series of both homotopy parameters: v p 0 v 0 p 1 v 1 p 2 v 2 · · · q 0 w 0 q 1 w 1 q 2 w 2 · · · , 2.10 where v 0 , v 1 , v 2 , . . . and w 1 , w 2 , . . . are unknown functions to be determined by the RBHPM method. In addition, w 0 is an arbitrary trial function, which is chosen in order to improve the RBHPM convergence in the same way as the trial function u 0 does for the HPM method 50 . On one side, the order of the approximation for the HPM method is determined by the highest power of p considered in the formulation. On the other side, the order for the RBHPM method is given as i, k , where i and k are the highest power of p and q employed in the numerator and denominator of 2.10 , respectively.
The limit of 2.10 , when p → 1 and q → 1, provides an approximate solution for 2.3 given as Journal of Applied Mathematics The limit above exists in the case that both limits exist.

Convergence of RBHPM Method
In order to analyse the convergence of the RBHPM method, 2.8 is rewritten as After applying the inverse operator, L −1 , on both sides of 3.1 , we obtain By assuming that see 2. 10 and after substituting 3.3 in the right-hand side of 3.2 in the following form: The exact solution of 2.3 is obtained in the limit when p → 1 and q → 1 in 3.4 , which results in Then, according to Banach Fixed Point Theorem, N has a unique fixed point u, that is, N u u. Assume that the sequence generated by the RBHPM method can be written as 3.8 and suppose that W 0 v 0 /β ∈ B r u , where B r u {w * ∈ X| w * − u < r}, then, under these conditions: Proof. i By inductive approach, for n 1 we have Assuming that W n−1 − u ≤ γ n−1 w 0 − u , as induction hypothesis, then Using i , we have ii Because of W n − u ≤ γ n w 0 − u and lim n → ∞ γ n 0, lim n → ∞ W n − u 0, that is, lim n → ∞ W n u. 3.12

Padé Approximants
A rational approximation to f x on a, b is the quotient of two polynomials P N x and Q M x of degrees N and M, respectively. We use the notation R N,M x to denote this quotient. The R N,M x Padé approximations 51, 52 to a function f x are given by The method of Padé requires f x and its derivative to be continuous at x 0. The polynomials used in 4.1 are expressed as

4.2
The polynomials in 4.2 are constructed so that f x and R N,M x agree at x 0 as well as their derivatives up to N M agree at x 0. A special case occurs for Q 0 x 1, wherein the approximation in 4.1 becomes the Maclaurin expansion for f x . For a fixed value of N M the error is smallest when P N x and Q M x have the same degree or when P N x has degree one higher than Q M x .
Notice that the constant coefficient of Q M is q 0 − 1. This is permissible, because it can be noted that 0 and R N,M x are not changed when both P N x and Q M x are divided by the same constant. Hence the rational function R N,M x has N M 1 unknown coefficients. Assume that f x is analytic and has the Maclaurin expansion f x a 0 a 1 x a 2 x 2 · · · a k x k · · · . 4.3 And The lower index j N M 1 in the summation on the right side of 4.4 is chosen because the first N M derivatives of f X and R N,M x should agree at x 0.
When the left side of 4.4 is multiplied out and the coefficients of the powers of x i are set equal to zero for k 0, 1, 2, . . . , N M, the result is a system of N M 1 linear equations: Notice that the sum of the subscripts of the terms of each product is the same in each equation, and it increases consecutively from 0 to N M. The M equations in 4.6 involve only the unknowns q 1 , q 2 , . . . , q M and must be firstly solved. Then the equations in 4.5 are used successively to find p 1 , p 2 , . . . , p N 51 .

Laplace-Padé Transform and RBHPM Method Coupling
The coupling of Laplace transform and Padé approximant 46 is used in order to recover part of the lost information due to the truncated power series 51, 53-60 . The process can be recast as follows.
1 First, Laplace transformation is applied to power series.
2 Next, s is substituted by 1/x in the resulting equation.
3 After that, we convert the transformed series into a meromorphic function by forming its Padé approximant of order N/M . N and M are arbitrarily chosen, but they should be of smaller value than the order of the power series. In this step, the Padé approximant extends the domain of the truncated series solution to obtain better accuracy and convergence. 4 Then, x is substituted by 1/s. 5 Finally, by using the inverse Laplace s transformation, we obtain the modified approximate solution.
In order to improve the approximate rational solutions generated by the RBHPM method, we propose to apply the Laplace-Padé method, separately, to the denominator and the numerator in 2.11 , only when they are expressed as power series. We will denominate to this process as the Laplace-Padé rational biparameter homotopy perturbation method LPRBHPM .

Solution Calculated by RBHPM
We establish the following homotopy equation: where we have chosen a 0.25. Besides, the linear operator L is given as and the trial function is We suppose that the solution of 6.3 is of order 2, 2 , which is expressed as follows where the trial w-function is chosen as w 0 1. After substituting 6.6 into 6.3 , regrouping, and equating the terms having the following powers: p 0 q 0 , p 1 , p 2 , q 1 , and q 2 , it can be solved for v 0 , v 1 , v 2 , w 1 , and w 2 . In order to fulfil the initial conditions of 6.1 , it follows that v 0 0 0, v 1 0 0, v 2 0 0, w 1 0 0, and w 2 0 0. The results are recast in the following system of differential equations: By substituting 6.8 into 6.6 , and calculating the limits when p → 1 and q → 1, we obtain 6.9

Solution by HPM
We establish the following homotopy equation: where the linear operator L and the trial function u 0 are 6.4 and 6.5 , respectively. Substituting 2.6 into 6.10 , regrouping, and equating the terms with identical powers of p, it can be solved for v 0 , v 1 , v 2 , and so on in order to fulfil initial conditions from v 0 y 0 0, it follows that v 0 0 0, v 1 0 0, v 2 0 0 and so on . The result is recast in the following system of differential equations: . . ..

6.12
By using 6.12 and 2.6 and calculating the limit when p → 1, we obtain the second and four order approximations 14 respectively.

Case Study 2: Van Der Pol Oscillator
Consider the Van der Pol Oscillator problem 8, 48 without external forcing u u u u 2 u 0, u 0 0, u 0 1, 7.1

Solution by the RBHPM Method
We establish the following homotopy equation: where the linear operator L is and the trial function is We assume that the solution for 7.2 is order 2, 2 , which is expressed as follows: where the trial w-function is chosen as w 0 1. Substituting 7.5 into 7.2 , regrouping, and equating the terms having the following powers: p 0 q 0 , p 1 , p 2 , q 1 , and q 2 , it can be solved for v 0 , v 1 , v 2 , w 1 , and w 2 . In order to fulfil the initial conditions of 7.1 , it follows that v 0 0 0, v 0 0 1 and the rest are v i 0 0, v i 0 0, w i 0 0 and w i 0 0, i 1, 2, 3, . . . . Then, considering that w 0 1, we establish the following system: By solving 7.6 , we obtain v 0 t, v 1 − 1 12 at 2 t 2 2t 6 , v 2 1 2520 a 2 t 3 30t 4 77t 3 315t 2 210t 420 , where k 1 and k 2 are integration constants. Substituting 7.7 into 7.5 , and calculating the limits when p → 1 and q → 1, we obtain u t lim We select the parameters as a 27/100, k 1 22/53, and k 2 −17/39 by using the numerical procedure reported in 35-37 . In order to guarantee the validity of the approximate solution 7.8 for large t, the quotient of series solutions is transformed by the Laplace-Padé after-treatment. The procedure is applied separately to numerator and denominator of the expression in 7.8 . First, Laplace transformation is applied to numerator of 7.8 and then 1/t is written in place of s in the equation. Then, the Padé approximant 3/3 is applied and 1/s is written in place of t. Finally, by using the inverse Laplace s transformation, we obtain the modified approximate solution for numerator

Solution by HPM
We establish the following homotopy equation: where the linear operator L and the trial function u 0 are 7.3 and 7.4 , respectively. Substituting 2.6 into 7.12 , regrouping, and equating the terms with identical powers of p. In order to fulfil initial conditions of 7.1 , it follows that v 0 0 0, v 0 0 1,

Journal of Applied Mathematics
13 and the rest are v i 0 0 and v i 0 0 i 1, 2, 3, . . . . Then, we establish the following system:

7.13
By solving 7.13 , we obtain 7.14 By substituting the solutions from 7.14 into 2.6 and calculating the limit when p → 1, we can obtain the second order approximation u t lim In the same way, we obtain the four order approximation RBHPM (6.9) HPM (6.14) Exact (6.2) HPM (6.13) Figure 1: Exact solution 6.2 diagonal cross for 6.1 and its approximate solutions 6.9 solid line , 6.13 solid circles , and 6.14 solid diamonds .

Numerical Simulation and Discussion
On one hand, Figure 1 and Table 1 show a comparison between the exact solution 6.2 for the Riccati nonlinear differential equation 6.1 and the analytic approximations 6.9 , 6.13 , and 6.14 . The RBHPM method of order 2, 2 gives the smallest average absolute error A.A.E. of all solutions, followed by the solutions obtained with HPM of order 4 and 2. It is remarkable to observe that the RBHPM method generates an accurate rational approximation that successfully replicates the asymptotic behaviour of 6.2 .
On the other hand, Figure 2 and Table 2  In this case study, both LPRBHPM and RBHPM generate more accurate solutions than HPM. The coupling of Laplace and Padé with RBHPM was a key factor for improving the accuracy in the Van der Pol problem. Both numerator and denominator of the rational expression 7.8 were considered as truncated power series. Laplace-Padé after-treatment 46 was successfully applied, and it yields better accuracy and it allows for larger ranges of the domain. Additionally, the trial functions w 0 and u 0 play an important role in the behaviour of the RBHPM and LPRBHPM methods; therefore, future research must be done in order to understand what kinds of trial functions produce better results; in this context, the operators L and N will be a key aspect to consider. Furthermore, it is important to point out that the RBHPM and LPRBHPM methods do not resort to linearization, a perturbation parameter, or assumptions of weak nonlinearity, and it clearly results that the generated solution has a general character and it is more realistic compared to the method of simplifying the physical problems.
Finally, the homotopy in 2.4 can be replaced by homotopy schemes such as those reported in the literature 63-93 . This future line of research can lead us to improve the performance of RBHPM or LPRBHPM methods.

Conclusions
This paper presented the RBHPM and LPRBHPM methods as a novel tool with high potential to solve nonlinear differential equations. Furthermore, a comparison between the results of applying the proposed methods and HPM was given. For the Riccati nonlinear asymptotic problem, a comparison between the RBHPM and HPM methods was presented, showing how the RBHPM method generates highly accurate approximate solutions, similar to the results obtained when using the HPM method. Additionally, a Van der Pol Oscillator problem was tackled with the proposed methods and compared to solutions obtained by HPM. It resulted that the RBHPM and LPRBHPM methods produced the most accurate approximated solutions. Because RBHPM, LPRBHPM, and HPM are closely related, it is possible that differential equations solved by HPM can be solved by using RBHPM or LPRBHPM. Besides, further research should be done to apply the proposed methods to the calculation of approximate solutions of nonlinear partial differential equations, nonlinear fractional equations, and boundary value problems, among others. An important remark is that the RBHPM and LPRBHPM methods do not resort to any kind of linearization procedure or perturbation parameter. Thereupon, these methods promise to become important mathematical tools, useful for scientist and engineers working in the area of mathematical modelling and computer simulation.