On a New Reliable Algorithm

Abstract

An analytical expression for the solution of the prey-predator problem by an adaptation of the homotopy analysis method (HAM) is presented. The HAM is treated as an algorithm for approximating the solution of the problem in a sequence of time intervals; that is, HAM is converted into a hybrid numeric-analytic method called the multistage HAM (MSHAM). Comparisons between the MSHAM solutions and the fourth-order Runge-Kutta (RK4) numerical solutions reveal that the new technique is a promising tool for the nonlinear systems of ordinary differential equations.

1. Introduction

Most modelling of biological problems is characterized by systems of ordinary differential equations (ODEs). For example, the relationship of increasing and decreasing in the population of two kind of animals (such as rabbits and foxes) can be described by the so-called mathematical model of the prey-predator problem which is given by a system of nonlinear equations:

𝑑𝑥(𝑡)𝑑𝑡=𝑥(𝑡)(𝑎−𝑏𝑦(𝑡)),𝑥(𝑡₀)=𝑐₁,𝑑𝑦(𝑡)𝑑𝑡=−𝑦(𝑡)(𝑐+𝑑𝑥(𝑡)),𝑦(𝑡₀)=𝑐₂, (1.1) where 𝑥(𝑡) and 𝑦(𝑡) are, respectively, the populations of rabbits and foxes at the time 𝑡 and 𝑎,𝑏,𝑐, and 𝑑 are known coefficients (for more details, see [1]). Problems of this nature have been solved using classical numerical techniques such as Runge-Kutta and are now published in most textbooks on differential equations.

Authors in [2, 3] used the Adomian decomposition method (ADM) to handle the systems of prey-predator problem. Yusufoglu and Erbas [4] and Rafei et al. [5] employed the variational iteration method (VIM) to compute an approximation to the solution of the system of nonlinear differential equations governing the problem. Biazar et al. [6] used the power series method (PSM) to handle the systems. All the solutions above are in the form of convergent power series with polynomial base function.

In recent years, a great deal of attention has been devoted to study HAM, which was first envisioned by Liao in his Ph.D. thesis [7] for solving a wide range of problems whose mathematical models yield differential equation or system of differential equations [8–12]. HAM has successfully been applied to many situations. For example, Hayat and Sajid [13], Hayat and Khan [14], and Hayat et al. [15, 16] used HAM to solve many kinds of modelling in fluid problems. Xu and Liao [17] applied HAM to give the dual solutions of boundary layer flow over an upstream moving plate. Abbasbandy [18] used HAM to present solitary wave solutions to the Kuramoto-Sivashinsky equation. Sami Bataineh et al. [19] modified the HAM (MHAM) to show that Taylor series converge to the exact solution by expanding the coefficients variables using Taylor series. Alomari et al. [20] introduced a new reliable algorithm based on an adaptation of the standard HAM to solve Chen system. In recent years, this method has been successfully employed to solve many types of problems in science and engineering [21–28].

In this paper, we are interested to find the approximate analytic solution of the system of coupled nonlinear ODEs (1.1) by treating the HAM as an algorithm for approximating the solution of the problem in a sequence of time intervals. We shall call this technique as the multistage homotopy analysis method (MSHAM). The freedom of choosing the linear operator and the auxiliary parameter ℏ is still present in this modification. Different from the series solution in [2–5], the solution we present here uses a combination of base functions, namely, the polynomial and exponential functions, and it is effective for longer time. Comparison with the classical fourth-order Runge-Kutta (RK4) shall be made.

2. Solution Procedure

Consider the following general system of first-order ordinary differential equations (ODEs):

𝑦^′_𝑖=𝑓_𝑖(𝑡,𝑦₁,𝑦₂,…,𝑦_𝑛),𝑦_𝑖(𝑡₀)=𝑦_𝑖,0,𝑖=1,2,…,𝑛, (2.1) where 𝑓_𝑖 are (linear or nonlinear) real-valued functions, 𝑡₀∈ℜ is the initial condition, and 𝑦_𝑖,0∈ℜ.

2.1. Solution by HAM

In HAM [8], system (2.1) is first written in the form

𝑁_𝑖[𝑦_𝑖(𝑡)]=0,𝑖=1,2,…,𝑛, (2.2) where 𝑁_𝑖 are nonlinear operators, 𝑡 denotes the independent variable, and 𝑦_𝑖(𝑡) are the unknown functions. By means of generalizing the traditional homotopy method, Liao [8] constructs the so-called zeroth-order deformation equation:

(1−𝑞)𝐿[𝜙_𝑖(𝑡;𝑞)−𝑦_𝑖,0]=𝑞ℏ_𝑖𝑁_𝑖[𝜙_𝑖(𝑡;𝑞)], (2.3) where 𝑞∈[0,1] is an embedding parameter, ℏ_𝑖 are the nonzero auxiliary parameters, 𝐿 is an auxiliary linear operator, 𝑦_𝑖,0(𝑡) are the initial guesses of 𝑦_𝑖(𝑡), and 𝜙_𝑖(𝑡;𝑞) are the unknown functions. It is important to note that one has great freedom to choose auxiliary objects such as ℏ_𝑖 and 𝐿 in HAM. We note that, in the framework of HAM, the solution 𝑦_𝑖(𝑡) can be represented by many different base functions such as the polynomial functions, exponential functions and rational functions, and so forth. Obviously, when 𝑞=0 and 𝑞=1, both 𝜙_𝑖(𝑡;0)=𝑦_𝑖,0(𝑡) and 𝜙_𝑖(𝑡;1)=𝑦_𝑖(𝑡) hold. Thus, as 𝑞 increases from 0 to 1, the solution 𝜙_𝑖(𝑡;𝑞) varies from the initial guess 𝑦_𝑖,0(𝑡) to the solution 𝑦_𝑖(𝑡). Expanding 𝜙_𝑖(𝑡;𝑞) in Taylor series with respect to 𝑞, one has 𝜙_𝑖(𝑡;𝑞)=𝑦_𝑖,0(𝑡)+^∞_𝑚=1𝑦_𝑖,𝑚(𝑡)𝑞^𝑚, (2.4) where

𝑦_𝑖,𝑚1(𝑡)=𝜕𝑚!^𝑚𝜙_𝑖(𝑡;𝑞)𝜕𝑞^𝑚|_𝑞=0. (2.5) If the auxiliary linear operators, the initial guesses, the auxiliary parameters ℏ_𝑖, and the auxiliary function are so properly chosen, then the series (2.4) converges at 𝑞=1 and

𝜙_𝑖(𝑡;1)=𝑦_𝑖,0(𝑡)+^∞_𝑚=1𝑦_𝑖,𝑚(𝑡), (2.6) which must be one of the solutions of the original nonlinear equations, as proved by Liao [8].

Define the vector

𝑦_𝑖,𝑛(𝑡)={𝑦_𝑖,0(𝑡),𝑦_𝑖,1(𝑡),…,𝑦_𝑖,𝑛(𝑡)}. (2.7) Differentiating (2.3) 𝑚 times with respect to the embedding parameter 𝑞 and then setting 𝑞=0 and finally dividing them by 𝑚!, we have the so-called mth-order deformation equation:

𝐿[𝑦_𝑖,𝑚(𝑡)−𝜒_𝑚𝑦_{𝑖,𝑚−1}(𝑡)]=ℏ_𝑖𝐻_𝑖(𝑡)𝑅_𝑖,𝑚(𝑦_{𝑖,𝑚−1}), (2.8) where

𝑅_𝑖,𝑚(𝑦_{𝑖,𝑚−1}1)=𝜕(𝑚−1)!^𝑚−1𝑁_𝑖[𝜙_𝑖(𝑡;𝑞)]𝜕𝑞^𝑚−1|_𝑞=0,𝜒_𝑚=0,𝑚≤1,1,𝑚>1. (2.9) It should be emphasized that 𝑦_𝑖,𝑚(𝑡)(𝑚≥1) is governed by the linear (2.8) with the linear boundary conditions that come from the original problem, which can be easily solved by symbolic computation softwares such as Maple and Mathematica. Also,

𝑦_𝑖(𝑡)≃𝜑_𝑖,𝑘(𝑡)=^𝑘_𝑚=0𝑦_𝑖,𝑚(𝑡), (2.10)

2.2. Solution by MSHAM

The approximate solutions (2.10) are generally, as shall be shown in the numerical experiments of this paper, not valid for large 𝑡. A simple way of ensuring validity of the approximations for large 𝑡 is to treat (2.3) and (2.8) as an algorithm for approximating the solutions of (2.1) in a sequence of intervals: the solution from [𝑡₀,𝑡) will be derived by subdividing this interval into [𝑡₀,𝑡₁),[𝑡₁,𝑡₂),…,[𝑡_𝑛−1,𝑡) and applying the HAM solution on each subinterval. The initial approximation in each interval is taken from the solution on the previous interval:

𝑦_𝑖,0(𝑡)=𝑦_𝑖(𝑡^∗)=𝑐^∗₁, (2.11) where 𝑡^∗ is the left-end point of each subinterval.

Now we solve (2.8) for unknowns 𝑦_𝑖,𝑗(𝑡)(𝑖=1,2,…;𝑗=1,2…,𝑚). In order to carry out the iteration in every subinterval of equal length ℎ, [𝑡₀,𝑡₁), [𝑡₁,𝑡₂), [𝑡₂,𝑡₃),…,[𝑡_𝑗−1,𝑡), we need to know the values of the following:

𝑦^∗_𝑖,0=𝑦_𝑖(𝑡^∗),𝑖=1,2,…. (2.12) But, in general, we do not have these information at our disposal except at the initial point 𝑡^∗=𝑡₀. A simple way for obtaining the necessary values could be by means of the previous 𝑛-term approximations 𝜑_1,𝑘,𝜑_2,𝑘,…,𝜑_𝑖,𝑘 of the preceding subinterval given by (2.10), that is,

𝑦^∗_1,0=𝜑_1,𝑘(𝑡^∗),𝑦^∗_2,0=𝜑_2,𝑘(𝑡^∗),…,𝑦^∗_𝑖,0=𝜑_𝑖,𝑘(𝑡^∗). (2.13)

3. Application

In this part, we apply the MSHAM for the prey-predator problem subject to the initial conditions

𝑥(𝑡^∗)=𝑐₁,𝑦(𝑡^∗)=𝑐₂. (3.1) We note that when 𝑡^∗=0, we have the initial condition of (1.1). Let the set of the base functions be

{𝑡^𝑛𝑒^−st∣𝑛≥0,𝑠≥0}. (3.2) So the solutions are

𝑥(𝑡)=𝑑_0,0+^+∞_𝑠=0^+∞_𝑛=0𝑑^𝑛_𝑛,𝑠𝑡^𝑛𝑒^−st,𝑦(𝑡)=𝑏_0,0+^+∞_𝑠=0^+∞_𝑛=0𝑏^𝑛_𝑛,𝑠𝑡^𝑛𝑒^−st, (3.3) where 𝑑^𝑛_𝑛,𝑠 and 𝑏^𝑛_𝑛,𝑠 are the coefficients. It is straightforward to choose

𝑥₀(𝑡)=𝑐₁,𝑦₀(𝑡)=𝑐₂, (3.4) as our initial approximations of 𝑥(𝑡) and 𝑦(𝑡), and the linear operator should be

𝐿[𝜙(𝑡;𝑞)]=𝜕𝜙(𝑡;𝑞)𝜕𝑡+𝜙(𝑡;𝑞) (3.5) with the property

𝐿[𝐴𝑒^−𝑡]=0, (3.6) where 𝐴 is the integration constant, which will be determined by the initial condition.

If 𝑞∈[0,1] and ℏ indicate the embedding and nonzero auxiliary parameters, respectively, then the zeroth-order deformation problems are of the following form:

(1−𝑞)𝐿[𝑥(𝑡;𝑞)−𝑥₀(𝑡)]=𝑞ℏ𝑁_𝑥[𝑥(𝑡;𝑞),𝑦(𝑡;𝑞)],(1−𝑞)𝐿[𝑦(𝑡;𝑞)−𝑦₀(𝑡)]=𝑞ℏ𝑁_𝑦[𝑥(𝑡;𝑞),𝑦(𝑡;𝑞)], (3.7) subject to the initial conditions

𝑥(𝑡^∗;𝑞)=𝑐₁,𝑦(𝑡^∗;𝑞)=𝑐₂, (3.8) in which we define the nonlinear operators 𝑁_𝑥 and 𝑁_𝑦 as

𝑁_𝑥[𝜕𝑥(𝑡;𝑞),𝑦(𝑡;𝑞)]=𝑥(𝑡;𝑞)−𝑁𝜕𝑡𝑥(𝑡;𝑞)(𝑎−𝑏𝑦(𝑡;𝑞)),_𝑦[𝜕𝑥(𝑡;𝑞),𝑦(𝑡;𝑞)]=𝑦(𝑡;𝑞)+𝜕𝑡𝑦(𝑡;𝑞)(𝑐−𝑑𝑥(𝑡;𝑞)). (3.9) For 𝑞=0 and 𝑞=1, the above zeroth-order equations (3.7) have the solutions

𝑥(𝑡;0)=𝑥₀(𝑡),𝑦(𝑡;0)=𝑦₀(𝑡),𝑥(𝑡;1)=𝑥(𝑡),𝑦(𝑡;1)=𝑦(𝑡). (3.10)(3.11) When 𝑞 increases from 0 to 1, then 𝑥(𝑡;𝑞) and 𝑦(𝑡;𝑞) vary from 𝑥₀(𝑡) and 𝑦₀(𝑡) to 𝑥(𝑡) and 𝑦(𝑡). Expanding 𝑥 and 𝑦 in Taylor series with respect to 𝑞, we have

𝑥(𝑡;𝑞)=𝑥₀(𝑡)+^∞_𝑚=1𝑥_𝑚(𝑡)𝑞^𝑚,𝑦(𝑡;𝑞)=𝑦₀(𝑡)+^∞_𝑚=1𝑦_𝑚(𝑡)𝑞^𝑚, (3.12) in which

𝑥_𝑚1(𝑡)=𝜕𝑚!^𝑚𝑥(𝑡;𝑞)𝜕𝑞^𝑚|_𝑞=0,𝑦_𝑚1(𝑡)=𝜕𝑚!^𝑚𝑦(𝑡;𝑞)𝜕𝑞^𝑚|_𝑞=0, (3.13) where ℏ is chosen in such a way that these series are convergent at 𝑞=1. Therefore, we have through (3.11) that

𝑥(𝑡)=𝑥₀(𝑡)+^∞_𝑚=1𝑥_𝑚(𝑡),𝑦(𝑡)=𝑦₀(𝑡)+^∞_𝑚=1𝑦_𝑚(𝑡). (3.14)

Define the vectors

𝑥(𝑡)={𝑥₀(𝑡),𝑥₁(𝑡),…,𝑥_𝑛(𝑡)},𝑦(𝑡)={𝑦₀(𝑡),𝑦₁(𝑡),…,𝑦_𝑛(𝑡)}. (3.15) Differentiating the zeroth-order equations (3.7) 𝑚 times with respect to 𝑞, then setting 𝑞=0, and finally dividing by 𝑚!, we have the mth-order deformation equations:

𝐿[𝑥_𝑚(𝑡)−𝜒_𝑚𝑥_𝑚−1(𝑡)]=ℏ𝑅_𝑥,𝑚(𝑥(𝑡),𝑦(𝑡)),𝐿[𝑦_𝑚(𝑡)−𝜒_𝑚𝑦_𝑚−1(𝑡)]=ℏ𝑅_𝑦,𝑚(𝑥(𝑡),𝑦(𝑡)), (3.16) with the following boundary conditions:

𝑥_𝑚(𝑡^∗)=0,𝑦_𝑚(𝑡^∗)=0, (3.17) for all 𝑚≥1, where

𝑅_𝑥,𝑚(𝑥(𝑡),𝑦(𝑡))=𝑥^′_𝑚−1−𝑎𝑥_𝑚−1+𝑏^𝑚−1_𝑖=0𝑥_𝑖(𝑡)𝑦_{𝑚−1−𝑖}𝑅(𝑡),_𝑦,𝑚(𝑥(𝑡),𝑦(𝑡))=𝑦^′_𝑚−1+𝑐𝑥_𝑚−1−𝑑^𝑚−1_𝑖=0𝑦_𝑖(𝑡)𝑥_{𝑚−1−𝑖}(𝑡). (3.18) This way, it is easy to solve the linear nonhomogeneous Equations (3.16) at general initial conditions by using Maple, one after the other in the order 𝑚=1,2,3,…. Thus we successfully have 𝑥₁(𝑡)=−ℏ𝑐₁𝑎+ℏ𝑐₁𝑏𝑐₂−𝑒^−𝑡(−ℏ𝑐₁𝑎+ℏ𝑐₁𝑏𝑐₂)𝑒^−𝑡∗,𝑦₁(𝑡)=ℏ𝑐₂𝑐−ℏ𝑐₂𝑑𝑐₁−𝑒^−𝑡(ℏ𝑐₂𝑐−ℏ𝑐₂𝑑𝑐₁)𝑒^−𝑡∗,𝑥₂(𝑡)=(−ℏ²𝑐₁𝑒^𝑡∗𝑏𝑐₂𝑡^∗−ℏ𝑐₁𝑒^𝑡𝑎+ℏ𝑐₁𝑒^𝑡𝑏𝑐₂+ℏ²𝑐₁𝑎²𝑒^𝑡−ℏ²𝑐₁𝑎²𝑒^𝑡∗−2ℏ²𝑐₁𝑎𝑏𝑐₂𝑒^𝑡+2ℏ²𝑐₁𝑎𝑒^𝑡∗𝑏𝑐₂+ℏ²𝑐₁𝑏²𝑐²₂𝑒^𝑡−ℏ²𝑐₁𝑏²𝑐²₂𝑒^𝑡∗+𝑐₁𝑏ℏ²𝑐₂𝑐𝑒^𝑡−𝑐₁𝑏ℏ²𝑐₂𝑐𝑒^𝑡∗−𝑐²₁𝑏ℏ²𝑐₂𝑑𝑒^𝑡+𝑐²₁𝑏ℏ²𝑐₂𝑑𝑒^𝑡∗+ℏ𝑐₁𝑒^𝑡∗𝑎−ℏ²𝑐₁𝑎𝑒^𝑡∗𝑡+ℏ²𝑐₁𝑒^𝑡∗𝑎𝑡^∗−ℏ𝑐₁𝑏𝑐₂𝑒^𝑡∗+ℏ²𝑐₁𝑒^𝑡∗𝑎²𝑡^∗−2ℏ²𝑐₁𝑒^𝑡∗𝑎𝑏𝑐₂𝑡^∗+ℏ²𝑐₁𝑒^𝑡∗𝑏²𝑐²₂𝑡^∗+ℏ²𝑐₁𝑒^𝑡∗𝑏𝑐₂𝑐𝑡^∗−ℏ²𝑐²₁𝑒^𝑡∗𝑏𝑐₂𝑑𝑡^∗+ℏ²𝑐₁𝑏𝑐₂𝑒^𝑡∗𝑡−ℏ²𝑐₁𝑎²𝑒^𝑡∗𝑡+2ℏ²𝑐₁𝑎𝑒^𝑡∗𝑏𝑐₂𝑡−ℏ²𝑐₁𝑏²𝑐²₂𝑒^𝑡∗𝑡−ℏ²𝑐₁𝑏𝑐₂𝑐𝑒^𝑡∗𝑡+ℏ²𝑐²₁𝑏𝑐₂𝑑𝑒^𝑡∗𝑡)𝑒^−𝑡,𝑦₂(𝑡)=(−2ℏ²𝑐₂𝑐𝑑𝑐₁𝑒^𝑡+ℏ²𝑐₂𝑑²𝑐²₁𝑒^𝑡+2ℏ²𝑐₂𝑒^𝑡∗𝑐𝑑𝑐₁−ℏ²𝑐₂𝑒^𝑡∗𝑑²𝑐²₁+𝑐₂𝑑ℏ²𝑐₁𝑎𝑒^𝑡−𝑐²₂𝑑ℏ²𝑐₁𝑏𝑒^𝑡−𝑐₂𝑑ℏ²𝑐₁𝑒^𝑡∗𝑎+𝑐²₂𝑑ℏ²𝑐₁𝑒^𝑡∗𝑏+𝑒^𝑡ℏ𝑐₂𝑐−𝑒^𝑡ℏ𝑐₂𝑑𝑐₁+ℏ²𝑐₂𝑒^𝑡∗𝑑²𝑐²₁𝑡^∗+ℏ𝑐₂𝑒^𝑡∗𝑑𝑐₁−ℏ²𝑐₂𝑒^𝑡∗𝑐𝑡^∗+ℏ²𝑐₂𝑒^𝑡∗𝑑𝑐₁𝑡^∗+ℏ²𝑐₂𝑒^𝑡∗𝑐²𝑡^∗−2ℏ²𝑐₂𝑒^𝑡∗𝑐𝑑𝑐₁𝑡^∗+ℏ²𝑐₂𝑒^𝑡∗𝑑𝑐₁𝑎𝑡^∗−ℏ²𝑐²₂𝑒^𝑡∗𝑐₁𝑏𝑑𝑡^∗+ℏ²𝑐₂𝑒^𝑡∗𝑐𝑡−ℏ²𝑐₂𝑒^𝑡∗𝑑𝑐₁𝑡−ℏ²𝑐₂𝑒^𝑡∗𝑐²𝑡+2ℏ²𝑐₂𝑒^𝑡∗𝑐𝑑𝑐₁𝑡−ℏ²𝑐₂𝑒^𝑡∗𝑑²𝑐²₁𝑡−ℏ²𝑐₂𝑑𝑐₁𝑒^𝑡∗𝑎𝑡+ℏ²𝑐²₂𝑐₁𝑏𝑑𝑒^𝑡∗𝑡−ℏ𝑐₂𝑒^𝑡∗𝑐+ℏ²𝑐₂𝑐²𝑒^𝑡−ℏ²𝑐₂𝑒^𝑡∗𝑐²)𝑒^−𝑡,⋮ (3.19) By the same way we can get the first fourth term to be as analytical approximate solution as 𝑥(𝑡)≃∑⁴_𝑖=0𝑥_𝑖(𝑡) and 𝑦(𝑡)≃∑⁴_𝑖=0𝑦_𝑖(𝑡) terms. Now we divide the interval [0,𝑇] to subintervals by time step Δ𝑡=0.01; Then we start from the initial conditions and we get the solution on the interval [0,0.01). Further, we take 𝑐₁=𝑥(0.01) and 𝑐₂=𝑦(0.01) and 𝑡^∗=0.01, so we get the solution on the new interval [0.01,0.02), and so on. Therefore, by choosing this initial approximation on the starting of each interval, the solution on the whole interval should be continuous. It is worth mentioning that if we take 𝑡^∗=0 and we fixed 𝑐₁ and 𝑐₂, then the solution will be the standard HAM solution which is not effective at large value of 𝑡.

4. Analysis of Results

In this section, we compare the fourth term of the MSHAM solution using step size 0.01 with RK4 using step size Δ𝑡=0.001 for the following cases.

Case 1. 𝑎=1, 𝑏=1, 𝑐=0.1, 𝑑=1, and the initial conditions 𝑥(0)=14 and 𝑦(0)=18.

Case 2. 𝑎=0.1, 𝑏=1, 𝑐=1, 𝑑=1, and the initial conditions 𝑥(0)=14 and 𝑦(0)=18.

Case 3. 𝑎=1, 𝑏=1, 𝑐=0.1, 𝑑=1, and the initial conditions 𝑥(0)=16 and 𝑦(0)=10.

Case 4. 𝑎=0.1, 𝑏=1, 𝑐=1, 𝑑=1, and the initial conditions 𝑥(0)=16 and 𝑦(0)=10.

Case 5. 𝑎=1, 𝑏=1, 𝑐=1, 𝑑=1, and the initial conditions 𝑥(0)=3 and 𝑦(0)=2.

Figure 1 presents the ℏ-curves for 5th-order of approximations of the prey-predator problem at different cases. It is clear that ℏ=−1 is in the convergent region which is parallel to the 𝑥-axis. Figure 2 shows the MSHAM solutions in comparison with the RK4 and HAM solutions for the time span 𝑡∈[0,30] using 16 Digits in Maple software. This 4th-term HAM solution is not accurate enough when compared with the RK4 solution but the 4th-term MSHAM solution works very well. The results obtained in [2–5] are effective for 𝑡≤1 but the result in MSHAM is effective for 𝑡∈[0,30], that means the MSHAM is more effective than the methods in [2–5]. Moreover, the absolute errors between the new algorithm and RK4 using the same benchmark for Cases 2 and 5 are presented in Tables 1 and 2, respectively. These show high accuracy of the new method since the absolute error is up to 10⁻¹¹ in Case 2 and 10⁻⁷ in Case 5, which is not possible to get if the linear operator is 𝜕𝜙(𝑡;𝑞)/𝜕𝑡.

Table 1 

The absolute error from the MSHAM0.001 and RK40.001 for Case 2 of the prey-predator problem.

Table 2 

The absolute error from the MSHAM0.001 and RK40.001 for Case 5 of the prey-predator problem.

(a)

(b)

(c)

(d)

Figure 1 

The ℏ-curves for 5th-order of approximations. (a) and (b) The curves of 𝑥 for Cases 1 and 5. (c) and (d) The curves of 𝑦 for Cases 1 and 5.

(a)

(b)

(c)

(d)

(e)

Figure 2 

Comparisons between the populations of rabbits 𝑥 and foxes 𝑦 at time 𝑡 for (a) Case 1, (b) Case 2, (c) Case 3, (d) Case 4, and (e) Case 5, using the 4th-term classical HAM, 4th-term MSHAM (Δ𝑡=0.01), and RK4 (Δ𝑡=0.001).

5. Conclusions

In this paper, a prey-predator problem was simulated accurately by MSHAM. MSHAM has the advantages of giving an analytical form of the solution within each time interval which is not possible in purely numerical techniques like RK4, which provides solution only at the two ends of a given time interval, and provided that the interval is chosen small enough for convergence. The method also gives the freedom to choose the auxiliary linear operator and the auxiliary parameter ℏ. The present technique offers an explicit time marching algorithm that works accurately over such a bigger time span. The results presented in this paper suggest that MSHAM is also readily applicable to more complex cases of nonlinear ordinary differential equations.

Acknowledgments

The authors would like to acknowledge the financial supports received from Universiti Kebangsaan Malaysia under the Science Fund UKM-ST-06-FRGS0008-2008, and UKM-GUP-BTT-07-25-174. The authors would also like to thank the reviewers for their valuable comments and suggestions.