Modified HPMs Inspired by Homotopy Continuation Methods

Abstract

Nonlinear differential equations have applications in the modelling area for a broad variety of phenomena and physical processes; having applications for all areas in science and engineering. At the present time, the homotopy perturbation method (HPM) is amply used to solve in an approximate or exact manner such nonlinear differential equations. This method has found wide acceptance for its versatility and ease of use. The origin of the HPM is found in the coupling of homotopy methods with perturbation methods. Homotopy methods are a well established research area with applications, in particular, an applied branch of such methods are the homotopy continuation methods, which are employed on the numerical solution of nonlinear algebraic equation systems. Therefore, this paper presents two modified versions of standard HPM method inspired in homotopy continuation methods. Both modified HPMs deal with nonlinearities distribution of the nonlinear differential equation. Besides, we will use a calcium-induced calcium released mechanism model as study case to test the proposed techniques. Finally, results will be discussed and possible research lines will be proposed using this work as a starting point.

1. Introduction

Many important physical phenomena on the engineering and science fields are frequently modelled by nonlinear differential equations. Such equations are often difficult or impossible to solve in a closed way. Nevertheless, analytic methods to obtain approximate solutions have gained importance in recent years. Among the most employed methods we have the homotopy perturbation method (HPM) [1–9]; it was proposed by the Chinese mathematician Ji-Huan He, and it was introduced as a powerful tool to approach various kinds of nonlinear problems. As is well known, nonlinear phenomena appear in a wide variety of scientific fields, such as applied mathematics, physics, and engineering. Scientists in those disciplines are constantly faced with the task of finding solutions for linear and nonlinear ordinary differential equations, partial differential equations, and systems of nonlinear ordinary differential equations. In fact, there are several methods used to find approximate solutions to nonlinear problems like variational approaches [10–14], Tanh Method [15], exp-function [16], Adomian’s decomposition method [17, 18], parameter expansion [19], and so on. Nevertheless, the HPM method is powerful, relatively simple to use, and has been tested successfully in a variety of applications [13, 20–33]. The HPM method can be considered as combination of the classical perturbation technique and the homotopy (whose origin is in the topology), but has eliminated the limitations of the traditional perturbation methods. For example, the method does not need a small parameter or linearisation, in fact, only requires few iterations to obtain accurate solutions.

The basic idea of HPM method is to introduce a homotopy parameter 𝑝, which takes values from 0 to 1. When parameter 𝑝=0, the equation usually reduces to a simple, or trivial, equation to solve. Then 𝑝 gradually is increase to 1, producing a sequence of deformations, where every solution is close to the last one. Eventually at 𝑝=1, the system takes the original form of the equation and the final stage of deformation provides the desired solution. However, only a few iterations are needed to achieve a good accuracy.

Contributions to the HPM have been done in order to solve, in the most efficient way, certain types of nonlinear problems [21, 34–36]. Nevertheless, we will see it is possible to modify HPM by using analogies to the homotopy continuation methods (HCM) [37–52]. The HCM methods are a mathematical technique that allows to locate one or multiple roots from nonlinear algebraic equation systems (NAESs). For the HCM methods, the homotopy parameter plays a role similar to the homotopy parameter in the HPM methods; it serves to create a continuous deformation from, say 𝜆=0, where the homotopy equation has a trivial solution or simple to solve, up to 𝜆=1, where the homotopy equation recovers the form of the original nonlinear problem. In every deformation step for 𝜆 the previous solution is employed to find the next. Nevertheless, for some cases the homotopy simulation HCM is unable to locate a solution; in the HCM the nonlinearities in the equation or system of equations to solve play a determinant role on the probability to locate one or more solutions. Therefore, the HCM have been developed with nonlinear mitigation techniques by embedding the homotopy parameter in the equation to solve, generating a mitigation on the nonlinearities capable to increase the probability to find the solution even for initial points located far from the solution. Hence, due to the similarity of the HPM and HCM, we propose that it is possible to distribute the nonlinearities of the differential equation, in the HPM to create iterations easier to solve.

In this work, we have applied the new proposed modification of HPM to find an approximate solution of the model of calcium-induced calcium release (CICR) mechanism. A detailed physical interpretation of the CICR mechanism is given by [53, 54].

This paper is organized as follows. In Section 2, we provide a brief review of the basic idea for the HPM and HCM. In Section 3, we propose a variant of the HPM based on the nonlinear distributions of the differential equation between subsequent iterations of the method. In Section 4, we propose a multiparameter version of the HPM. Section 5 provides a brief explanation of the CICR mechanism, then we solve it by the HPM and the proposed modifications. In Section 6, we discuss the results, summarize our findings, and suggest possible directions for future investigations. Finally, a brief conclusion is given in Section 7.

2. Basic Idea of HPMs and HCMs

HPM methods are a mathematical tool which bases its origin coupling the perturbation and homotopy methods. Therefore, in this section we will explore the basic concepts for the HPM and HCM in order to find analogies or similarities that allows to propose modifications to the HPM inspired on the HCMs.

2.1. Basic Idea of HPMs

In the HPM is considered that a nonlinear differential equation can be expressed as𝐴(𝑢)−𝑓(𝑟)=0,where𝑟∈Ω. (2.1) With the boundary condition𝐵𝑢,𝜕𝑢𝜕𝜂,where𝑟∈Γ, (2.2) where 𝐴 is a general differential operator, 𝑓(𝑟) is a known analytic function, 𝐵 is a boundary operator, and Γ is the boundary of the domain Ω. The 𝐴 operator, generally, can be divided into two operators, 𝐿 and 𝑁, where 𝐿 is the linear operator and 𝑁 is the nonlinear operator. Hence, (2.1) can be rewritten as𝐿(𝑢)+𝑁(𝑢)−𝑓(𝑟)=0. (2.3) Now, the homotopy function is𝐻𝑣,𝑝=1−𝑝𝐿(𝑣)−𝐿𝑢₀+𝑝𝐿(𝑣)+𝑁(𝑣)−𝑓(𝑟)=0where𝑝∈[0,1], (2.4) where 𝑢₀ is the initial approximation of (2.3) which satisfies the boundary conditions and 𝑝 is known as the perturbation homotopy parameter. Analysing (2.4) can be concluded that𝐻(𝑣,0)=𝐿(𝑣)−𝐿𝑢₀,𝐻(𝑣,1)=𝐿(𝑣)+𝑁(𝑣)−𝑓(𝑟). (2.5) We assume that the solution of (2.4) can be written as a power series of 𝑝 [1]𝑣=𝑝⁰𝑣₀+𝑝¹𝑣₁+𝑝²𝑣₂+⋯. (2.6) Adjusting 𝑝=1 results that the approximate solution for (2.1) is𝑢=lim_𝑝→1𝑣=𝑣₀+𝑣₁+𝑣₂+⋯. (2.7)

The series (2.7) is convergent on most cases, nevertheless, the convergence depends on the nonlinear operator 𝐴(𝑣) [1, 2].

2.2. Basic Idea of HCMs

In the area of the HCMs, the first step to formulate a homotopy is to establish the nonlinear algebraic equation to be solved; which is being defined as𝑓(𝑥)=0,where𝑓∶∈ℜ^𝑛⟶ℜ^𝑛, (2.8) where 𝑥 represents the variables of problem and 𝑛 is the number of variables. This equation plays the same role as the nonlinear differential equation (2.3) for the HPM.

The homotopy function can be represented as𝐻𝑓(𝑥),𝑝=0, (2.9) where 𝑓(𝑥) is the equilibrium equation and 𝑝 is the homotopy parameter just like the HPM. For this case, (2.9) represents any homotopy formulation that fulfils the following basic conditions:(i)for 𝑝=0 the solution of 𝐻⁻¹(0) is known or easy to obtain using numerical methods. This point is known as Homotopy’s initial point (𝑥₀),(ii)for 𝑝=1,𝐻(𝑓(𝑥),1)=𝑓(𝑥). Means that at 𝑝=1 all solutions for 𝑓(𝑥) are located,(iii)the path for 𝐻⁻¹(0) is a continuous function of 𝑝 within the range of 0≤𝑝≤1.

The homotopy path is the solution set of 𝐻⁻¹(0), representing a continuous curve that can be traced by some numerical continuation technique or path following method [48, 55]. An example of HCM is Newton’s homotopy. The formulation is𝐻𝑓(𝑥),𝑝=𝑓(𝑥)−1−𝑝𝑓𝑥₀=0, (2.10) where 𝑓(𝑥) is the equilibrium equation, 𝑥₀ is the initial point of the homotopy, and 𝑝 is the homotopy parameter.

If 𝑝=0,𝐻𝑓(𝑥),0=𝑓(𝑥)−𝑓𝑥₀=0, (2.11) where 𝑥₀ represents the initial point. It is simple to obtain or is chosen arbitrarily.

If 𝑝=1𝐻𝑓(𝑥),1=𝑓(𝑥)=0, (2.12) in this case the solution of 𝐻 is, precisely, the desired solution of the NA equation 𝑓(𝑥). That is, a continuous path is created between 𝑝=0 and 𝑝=1 which can be described as in Figure 1. In case that there is more than one solution: 𝑥^∗₁,𝑥^∗₂,𝑥^∗₃, and 𝑥^∗₄, for 𝑓(𝑥), there is a chance that those are connected to the same path.

Figure 1 

Homotopy path.

3. HPM with Nonlinearities Distributions

For the HCM methods the selection of the initial point is vital, since it may affect the probabilities to find a solution, even producing a failure on the search for the solution of the NAES. It is known that nonlinearities of the NAES usually directly affect the probabilities to find a solution, so a common strategy is to embed the homotopy parameter inside the NAES, multiplying 𝑝 by the nonlinear terms responsible for the nonlinearities that affect the NAES. This is equivalent to soften the nonlinearities inherent to the NAES and increasing the probabilities of success for the HCM. In other words, the nonlinearities of the NAES are distributed along the homotopy path from 𝑝=0 to 𝑝=1, creating a smooth path for the homotopy. An example of this type of homotopy can be formulated as𝐻𝑓(𝑥),𝑝=1−𝑝𝑥−𝑥₀+𝑝𝑓𝑥,𝑝=0, (3.1) where 𝑓(𝑥) is the NA equation, 𝑥₀ is the initial point of the homotopy, and 𝑝 is the homotopy parameter. It can be seen that when 𝑝=0, then𝐻𝑓(𝑥),0=𝑥−𝑥₀=0, (3.2) when 𝑝=1, then𝐻𝑓(𝑥),1=𝑓(𝑥,1)=𝑓(𝑥)=0, (3.3) where for 𝑝 between zero and one, the homotopy parameter reduce the nonlinear effect from nullifying it at 𝑝=0 until giving them their original state at 𝑝=1. This kind of strategy may provide better convergence to the solutions.

Therefore, using the previously exposed for the HCMs, it is proposed to distribute nonlinearities of 𝑁 and 𝑓 among sequential iterations of the HPM.𝐻𝑣,𝑝=1−𝑝𝐿(𝑣)−𝐿𝑢₀+𝑝𝐿(𝑣)+𝑁𝑣,𝑝−𝑓𝑟,𝑝=0,where𝑝∈[0,1]. (3.4)

The homotopy function (3.4) is essentially the same as (2.4), except for the nonlinear operator 𝑁 and function 𝑓 which contain embedded the homotopy parameter. We propose, first, to subdivide 𝑁 and 𝑓 into a series of terms and multiply 𝑝^𝑖 by the most nonlinear terms, where 𝑖 is an integer number greater or equal to zero. The 𝑖 power is selected according to how much displacement is desired in the iterations for the corresponding nonlinear term of 𝑁 or 𝑓.

As mentioned before, embedding the homotopy parameter 𝑝 within the differential equation is a strategy to redistribute the nonlinearities between the successive iterations of the HPM, we expect that this strategy helps to increase the probabilities to find the solution sought. The rest of the method is exactly the same as the standard procedure for the HPM. Therefore, we establish that𝑣=𝑝⁰𝑣₀+𝑝¹𝑣₁+𝑝²𝑣₂+⋯, (3.5) adjusting 𝑝=1 turns out that the approximate solution of (2.1) is𝑢=lim_𝑝→1𝑣=𝑣₀+𝑣₁+𝑣₂+⋯. (3.6)

An advantage of this procedure is that given the distribution of nonlinearities from the differential equation over the successive iterations of (3.5), less complex analytic approximations may be obtained than the generated by the original standard of the HPM. This method will be designated as nonlinear distribution HPM (NDHPM).

The convergence of the method will be tested on the next subsection taking as reference the exposed in [56, 57].

3.1. Convergence

We construct a homotopy 𝑣(𝑟,𝑝)∶Ω×[0,1]→ℜ^𝑛 which satisfies𝐻𝑣,𝑝=𝐿(𝑣)−𝐿𝑢₀+𝑝𝐿𝑢₀+𝑝𝑁𝑣,𝑝−𝑓𝑟,𝑝=0, (3.7) where 𝑢₀ is the initial approximation of (2.1).

Let one writes (3.7) in the following form:𝐿(𝑣)=𝐿𝑢₀+𝑝𝑓𝑟,𝑝−𝑁𝑣,𝑝−𝐿𝑢₀=0. (3.8) Applying the inverse operator, 𝐿⁻¹, to both sides of (3.8) we obtain𝑣=𝑢₀+𝑝𝐿⁻¹𝑓𝑟,𝑝−𝐿⁻¹𝑁𝑣,𝑝−𝑢₀. (3.9) Suppose that𝑣=^∞_𝑖=0𝑝^𝑖𝑣_𝑖, (3.10) substituting (3.10) in the right-hand side of (3.9) in the following form:𝑣=𝑢₀+𝑝𝐿⁻¹𝑓𝑟,𝑝−𝐿⁻¹𝑁^∞_𝑖=0𝑝^𝑖𝑣_𝑖−𝑢₀. (3.11) If 𝑝→1, the exact solution may be obtained by using𝑢=lim_𝑝→1𝑣=𝐿⁻¹𝑓(𝑟)−𝐿⁻¹𝑁^∞_𝑖=0𝑣_𝑖=𝐿⁻¹𝑓(𝑟)−^∞_𝑖=0𝐿⁻¹𝑁𝑣_𝑖. (3.12)

Now, to study the convergence of the method we use the Banach’s fixed point theorem as done in [56]. Besides, in [56] a study of convergence for the standard HPM is performed, giving as result exactly the same equation (3.12). Therefore, it can be concluded that the NDHPM, essentially, has the same convergence as the standard method; except that it allows the redistribution of the nonlinearities between the iterations of the method.

4. Multiparameter HPM

In the HCM method there are multiparameter versions [52]. Those multiparameter homotopy methods have certain advantages like the distribution of nonlinearities, bifurcations evasion, sharp folds, singularities, location of complex solutions, among others [55, 58]. We expect the creation of a multiparameter HPM may offer practical advantages, like the strategic distribution of nonlinearities between successive iterations of the HPM. The homotopy perturbation function is reformulated as𝐻𝑣,𝑝=1−𝑝₁𝐿(𝑣)−𝐿𝑢₀+𝑝₁𝐿(𝑣)+𝑁𝑣,𝑝₂−𝑓𝑟,𝑝₂=0, (4.1) where 𝑝₁∈[0,1] and 𝑝₂∈[0,1], both of them being homotopy parameters.

The homotopy formulation (4.1) is similar to (2.4), except that now there are two homotopy parameters 𝑝₁ and 𝑝₂; the first homotopy parameter 𝑝₁ formulates the basic homotopy (2.4) and the second homotopy parameter 𝑝₂ helps to distribute the nonlinearities of the differential equation between successive iterations of the HPM. Hence, the main benefit is that the nonlinearity degree can be redistributed for the resultant approximation, with respect to the obtained results with the original HPM (see (2.4)).

Just like in the homotopy (3.4), we propose to subdivide 𝑁 and 𝑓 in a set of terms and multiply 𝑝^𝑖₂ by the most nonlinear terms, where 𝑖 is an integer number greater or equal to zero. The 𝑖 power is selected according to how much displacement is desired in the iterations the corresponding term 𝑁 or 𝑓.

The rest of the method is the basically the standard HPM method, except that the power series are done in a multivariable way. Therefore, it is established that𝑣=𝑣₀+𝑣₁𝑝₁+𝑣₂𝑝₂+𝑣₃𝑝₁𝑝₂+𝑣₄𝑝²₁+𝑣₅𝑝²₂+𝑣₆𝑝₁𝑝²₂+𝑣₇𝑝²₁𝑝₂+⋯, (4.2) adjusting 𝑝=1 results that the approximate solution for (2.1) is𝑢=lim_{𝑝1→1𝑝2→1}𝑣=𝑣₀+𝑣₁+𝑣₂+⋯. (4.3)

This method will be designated as multiparameter HPM (MHPM).

4.1. Convergence

We construct a homotopy 𝑣(𝑟,𝑝₁,𝑝₂)∶Ω×[0,1]→ℜ^𝑛 which satisfies𝐻𝑣,𝑝=𝐿(𝑣)−𝐿𝑢₀+𝑝₁𝐿𝑢₀+𝑝₁𝑁𝑣,𝑝₂−𝑓𝑟,𝑝₂=0, (4.4) where 𝑢₀ is the initial approximation of (2.1).

Let one writes (4.4) in the following form:𝐿(𝑣)=𝐿𝑢₀+𝑝₁𝑓𝑟,𝑝₂−𝑁𝑣,𝑝₂−𝐿𝑢₀=0. (4.5) Applying the inverse operator 𝐿⁻¹ to both sides of (4.5), we obtain𝑣=𝑢₀+𝑝₁𝐿⁻¹𝑓𝑟,𝑝₂−𝐿⁻¹𝑁𝑣,𝑝₂−𝑢₀. (4.6) Substituting (4.2) in the right-hand side of (4.6) in the following form:𝑣=𝑢₀+𝑝₁𝐿⁻¹𝑓𝑟,𝑝₁,𝑝₂−𝐿⁻¹𝑁𝑣₀+𝑣₁𝑝₁+𝑣₂𝑝₂+𝑣₃𝑝₁𝑝₂+𝑣₄𝑝²₁+𝑣₅𝑝²₂+𝑣₆𝑝₁𝑝²₂+⋯−𝑢₀. (4.7) If 𝑝₁→1 and 𝑝₂→1, the exact solution may be obtained by using𝑢=lim_{𝑝1→1𝑝2→1}𝑣=𝐿⁻¹𝑓(𝑟)−𝐿⁻¹𝑁^∞_𝑖=0𝑣_𝑖=𝐿⁻¹𝑓(𝑟)−^∞_𝑖=0𝐿⁻¹𝑁𝑣_𝑖. (4.8)

Now, to study the convergence of the method we use the Banach’s fixed point theorem as done in [56]. Besides, in [56] a study of the convergence for the standard HPM is performed, resulting in exactly the same equation (4.8). Therefore, it can be concluded that the NDHPM has, essentially, the same convergence as the standard method except that allows the redistribution of the nonlinearities between the iterations of the method.

5. Solution of the CICR Model

There are a number of phenomena in biological sciences, where the precursor of a particular process is the appearance of a travelling wave of chemical concentration, mechanical deformation, electrical signals, and so on [53, 54]. There are, for instance, both chemical and mechanical waves that propagate on the surface of many vertebrate eggs. In the case of the egg of Medaka fish, a Calcium (Ca⁺⁺) wave sweeps over the surface; it emanates from the point of sperm entry. Another example, related to interacting populations, is the progressing wave of an epidemic, on which, for instance, the rabies epizootic spreading a country. Another example is the movement of microorganisms moving into a food source chemotactically directed. The existence of wave phenomena in biomedical sciences requires a detailed study of travelling waves, and the search for analytic solutions of the equations that govern them. A kind of biological wave is the calcium-induced calcium release mechanism (CICR). This mechanism is relevant, for instance, to understand how the membrane enclosing certain fertilized amphibian eggs works. Due to the importance of the CICR mechanism, this work proposes to implement an approximate solution with good accuracy that describes the behaviour of such process.

During the development of living systems there is almost continual interchange of information at both inter and intra cellular level. Embryogenesis is an example on how such communication is necessary for a sequential development. Propagating waveforms of varied biochemical concentrations are the transmission medium of such information.

A biochemical switch is a mechanism whereby sufficiently large perturbation from one steady state can move a system to another steady state. An important example, which arises experimentally, is known as the calcium stimulated, calcium release mechanism. In this process, if calcium (Ca⁺⁺) is perturbed above a given threshold concentration, it causes the further release of sequestered calcium, that is, the system moves to another steady state. This happens, for instance, in certain calcium sites on the membrane enclosing fertilised amphibian eggs.

The model represents the kinetics model equations of the calcium resequestration and the autocatalytic release of calcium [53, 54]:𝑦^′(𝑥)𝑘=𝐴+¹𝑦²1+𝑦²−𝑘₂𝑦, (5.1) where 𝑥 represents time, 𝑦 represents calcium concentration, 𝑘₁ and 𝑘₂ are constants, and 𝐴 represents a small leakage into the membrane.

In this work the mechanism for the CICR is solved just as a mathematical case to explore the variants of the HPM proposed here.

5.1. Using Basic HPM

In order to analyse (5.1), it should be reordered𝑦′(𝑥)+𝑦²𝑦′(𝑥)+𝑘₂𝑦³−𝐴+𝑘₁𝑦²+𝑘₂𝑦−𝐴=0. (5.2)

To exemplify, we assign the following values 𝐴=1,𝑘1=3,𝑘2=1, and the initial condition 𝑦(0)=0.9. Therefore, (5.2) is updated becoming𝑦^′(𝑥)+𝑦²𝑦^′(𝑥)+𝑦³−4𝑦²+𝑦−1=0,where𝑦(0)=0.9. (5.3)

Equation (5.3) represents the transition from a state given by the initial condition 𝑦(0)=0.9 until a steady state for the calcium concentration. Now, the linear part is separated𝐿𝑦=𝑦′(𝑥)+𝑦−1, (5.4) and its nonlinear part𝑁𝑦=𝑦²(𝑥)𝑦^′(𝑥)+𝑦³−4𝑦². (5.5)

Adjustment constants are added to (5.4), the result is𝐿_𝑐𝑦=𝑎𝑦^′(𝑥)+𝑏𝑦+𝑐−1, (5.6) here 𝑎,𝑏, and 𝑐 are adjustment constants.

Now, we establish the homotopy equation1−𝑝𝐿_𝑐𝑦+𝑝𝐿𝑦+𝑁𝑦=0, (5.7) where 𝑝 is the homotopy parameter.

We suppose that solution for (5.7) has the form𝑦(𝑥)=𝑦₀(𝑥)+𝑝𝑦₁(𝑥)+𝑝²𝑦₂(𝑥)+⋯. (5.8)

Adjusting 𝑝=1, we can obtain an approximate solution𝑦(𝑥)=𝑦₀(𝑥)+𝑦₁(𝑥)+𝑦₂(𝑥)+⋯. (5.9)

Substituting (5.8) into (5.7) and equaling terms having potentials in the same order as 𝑝, it can be solved for 𝑦₀(𝑥),𝑦₁(𝑥),𝑦₂(𝑥), and so on (in order to fulfil initial conditions from 𝑦(0)=0.9, it follows that 𝑦₀(0)=0.9,𝑦₁(0)=0,𝑦₂(0)=0, and so on). Therefore, the following differential equations are established:𝑎𝑦^′₀(𝑥)+𝑏𝑦₀(𝑥)+𝑐−1=0,𝑎𝑦^′₁(𝑥)+𝑏𝑦₁(𝑥)+𝑦²₀(𝑥)−𝑎+1𝑦^′₀(𝑥)+𝑦³₀(𝑥)−4𝑦²₀(𝑥)+(1−𝑏)𝑦₀(𝑥)⋮−𝑐=0, (5.10) where terms 𝑦₂ and subsequent are ignored because the first two terms (𝑦₀(𝑥) and 𝑦₁(𝑥)) contribute to obtain a good accurate solution. In a similar process presented in [20, 59], we use NonlinearFit (given a total of 𝑘-samples from the exact model, the NonlinearFit command finds values of the approximate model parameters such that the sum of the squared 𝑘-residuals is minimized). command from Maple (Release 15) and command “convert” (with option “rational”), to set the adjustment parameters (this process will repeat for the rest of the examples in this paper), giving as a result:𝑦(𝑥)=−73−97exp86𝑥+33115−36exp33𝑥+1953+−14293+5537𝑥−53exp33𝑥38 (5.11)

5.2. Using the NDHPM

The NDHPMs establish that1−𝑝𝐿_𝑐𝑦+𝑝𝐿𝑦+𝑁_𝑖𝑦,𝑝−𝑓(𝑟)=0, (5.12) where 𝑝 is the homotopy parameter, 𝐿_𝑐 is established in (5.6), 𝑓(𝑟)=0, and the nonlinear part 𝑁_𝑖(𝑦,𝑝) (𝑖=1,2,3,4) is one of the following representative cases:𝑁₁𝑦,𝑝=𝑦²𝑦^′(𝑥)+𝑦³−4𝑦²𝑁𝑝,₂𝑦,𝑝=𝑦²𝑦^′(𝑥)+𝑦³𝑝−4𝑦²𝑁𝑝,₃𝑦,𝑝=𝑦²𝑦^′(𝑥)𝑝+𝑦³−4𝑦²𝑁𝑝,₄𝑦,𝑝=𝑦²𝑦^′(𝑥)𝑝+𝑦³𝑝−4𝑦². (5.13)

Consider that when lim_𝑝→1𝑁_𝑖=𝑁.

Now, the exact basic steps for the HPM method are followed to obtain the following differential equation system for cases 𝑁₁,𝑁₂,𝑁₃, and 𝑁₄:𝑝⁰∶𝑎𝑦^′₀(𝑥)+𝑏𝑦₀(𝑥)𝑝+𝑐−1=0,₁∶𝑎𝑦^′₁(𝑥)+𝑏𝑦₁(𝑥)+(−𝑎+1)𝑦^′₀(𝑥)+(1−𝑏)𝑦₀(𝑥)⋮𝑝−𝑐=0,₀∶𝑎𝑦^′₀(𝑥)+𝑏𝑦₀(𝑥)𝑝+𝑐−1=0,₁∶𝑎𝑦^′₁(𝑥)+𝑏𝑦₁(𝑥)+𝑦²₀(𝑥)−𝑎+1𝑦^′₀(𝑥)+(1−𝑏)𝑦₀(𝑥)⋮𝑝−𝑐=0,₀∶𝑎𝑦^′₀(𝑥)+𝑏𝑦₀(𝑥)𝑝+𝑐−1=0,₁∶𝑎𝑦^′₁(𝑥)+𝑏𝑦₁(𝑥)+(−𝑎+1)𝑦^′₀(𝑥)+𝑦³₀(𝑥)+(1−𝑏)𝑦₀(𝑥)⋮𝑝−𝑐=0,₀∶𝑎𝑦^′₀(𝑥)+𝑏𝑦₀(𝑥)𝑝+𝑐−1=0,₁∶𝑎𝑦^′₁(𝑥)+𝑏𝑦₁(𝑥)+(−𝑎+1)𝑦^′₀(𝑥)−4𝑦²₀(𝑥)+(1−𝑏)𝑦₀(𝑥)⋮−𝑐=0, (5.14) with initial conditions 𝑦₀(0)=0.9,𝑦₁(0)=0,𝑦₂(0)=0,…. This is done to fulfil initial conditions of (5.3).

In this case, the terms 𝑦₂ and subsequent are ignored because the first two terms (𝑦₀(𝑥) and 𝑦₁(𝑥)) provide good accuracy to the solution. Then, we obtain the following accurate approximations for all cases 𝑁₁,𝑁₂,𝑁₃ and 𝑁₄:𝑦(𝑥)=53+−14269+16𝑥−29exp15𝑥,𝑦(𝑥)=129−1264exp43𝑥−1211−118exp43𝑥+18195+−45−1637𝑥−19exp37𝑥𝑦(𝑥)=531,−54exp43𝑥−1185−71exp99𝑥+38195+−79−4494𝑥−29exp43𝑥𝑦(𝑥)33,=−13−41exp39𝑥+16195+−31−1229𝑥−13exp39𝑥,32 (5.15)(5.16)(5.17)(5.18) respectively.

5.3. Using the Multiparameter HPM

The homotopy equation is established according to the basic guidelines for the multiparameter HPM:1−𝑝₁𝐿_𝑐𝑦+𝑝₁𝐿𝑦+𝑁_𝑖𝑦,𝑝₂=0, (5.19) where 𝑝₁ is the first homotopy parameter, 𝐿_𝑐  is defined in (5.6), and we propose the nonlinearoperator 𝑁_𝑖(𝑦,𝑝₂)(𝑖=1,2,3,4) as the following representative case:𝑁₁𝑦,𝑝₂=𝑦²𝑦^′(𝑥)+𝑦³−4𝑦²𝑝₂,𝑁₂𝑦,𝑝₂=𝑦²𝑦^′(𝑥)𝑝₂+𝑦³𝑝²₂−4𝑦²𝑝₂,𝑁₃𝑦,𝑝₂=𝑦²𝑦^′(𝑥)𝑝₂+𝑦³𝑝₂−4𝑦²𝑝²₂,𝑁₄𝑦,𝑝₂=𝑦²𝑦^′(𝑥)𝑝²₂+𝑦³𝑝²₂−4𝑦²𝑝₂. (5.20)

From (4.2) we obtain𝑦(𝑥)=𝑦₀(𝑥)+𝑦₁(𝑥)𝑝₁+𝑦₂(𝑥)𝑝₂+𝑦₃(𝑥)𝑝₁𝑝₂+𝑦₄(𝑥)𝑝²₁+𝑦₅(𝑥)𝑝²₂+⋯. (5.21)

Now, we substitute (5.21) into (5.19) and follow the same basic steps to obtain the following differential equations for the first six terms of the multivariable power series (with initial conditions 𝑦₀(0)=0.9,𝑦₁(0)=0,𝑦₂(0)=0,𝑦₃(0)=0,𝑦₄(0)=0, and 𝑦₅(0)=0) for cases 𝑁₁,𝑁₂,𝑁₃, and 𝑁₄, we obtain𝑝⁰₁𝑝⁰₂∶𝑎𝑦^′₀(𝑥)+𝑏𝑦₀(𝑥)𝑝+𝑐−1=0,₁₁𝑝⁰₂∶𝑎𝑦^′₁(𝑥)+𝑏𝑦₁(𝑥)+(1−𝑎)𝑦^′₀(𝑥)+(1−𝑏)𝑦₀(𝑥)𝑝−𝑐=0,₀₁𝑝¹₂∶𝑎𝑦^′₂(𝑥)+𝑏𝑦₂(𝑥)𝑝=0,₁₁𝑝¹₂∶𝑎𝑦^′₃(𝑥)+𝑏𝑦₃(𝑥)+(1−𝑎)𝑦^′₂(𝑥)+(1−𝑏)𝑦₂(𝑥)+𝑦²₀(𝑥)𝑦^′₀(𝑥)−4𝑦²₀(𝑥)+𝑦³₀(𝑥)𝑝=0,₂₁𝑝⁰₂∶𝑎𝑦^′₄(𝑥)+𝑏𝑦₄(𝑥)+(1−𝑎)𝑦^′₁(𝑥)+(1−𝑏)𝑦₁(𝑥)𝑝=0,₀₁𝑝²₂∶𝑎𝑦^′₅(𝑥)+𝑏𝑦₅(𝑥)𝑝=0,₀₁𝑝⁰₂∶𝑎𝑦^′₀(𝑥)+𝑏𝑦₀(𝑥)𝑝+𝑐−1=0,₁₁𝑝⁰₂∶𝑎𝑦^′₁(𝑥)+𝑏𝑦₁(𝑥)+(1−𝑎)𝑦^′₀(𝑥)+(1−𝑏)𝑦₀(𝑥)𝑝−𝑐=0,₀₁𝑝¹₂∶𝑎𝑦^′₂(𝑥)+𝑏𝑦₂(𝑥)𝑝=0,₁₁𝑝¹₂∶𝑎𝑦^′₃(𝑥)+𝑏𝑦₃(𝑥)+(1−𝑎)𝑦^′₂(𝑥)+(1−𝑏)𝑦₂(𝑥)+𝑦²₀(𝑥)𝑦^′₀(𝑥)−4𝑦²₀(𝑥)𝑝=0,₂₁𝑝⁰₂∶𝑎𝑦^′₄(𝑥)+𝑏𝑦₄(𝑥)+(1−𝑎)𝑦^′₁(𝑥)+(1−𝑏)𝑦₁(𝑥)𝑝=0,₀₁𝑝²₂∶𝑎𝑦^′₅(𝑥)+𝑏𝑦₅(𝑥)𝑝=0,₀₁𝑝⁰₂∶𝑎𝑦^′₀(𝑥)+𝑏𝑦₀(𝑥)𝑝+𝑐−1=0,₁₁𝑝⁰₂∶𝑎𝑦^′₁(𝑥)+𝑏𝑦₁(𝑥)+(1−𝑎)𝑦^′₀(𝑥)+(1−𝑏)𝑦₀(𝑥)𝑝−𝑐=0,₀₁𝑝¹₂∶𝑎𝑦^′₂(𝑥)+𝑏𝑦₂(𝑥)𝑝=0,₁₁𝑝¹₂∶𝑎𝑦^′₃(𝑥)+𝑏𝑦₃(𝑥)+(1−𝑎)𝑦^′₂(𝑥)+(1−𝑏)𝑦₂(𝑥)+𝑦²₀(𝑥)𝑦^′₀(𝑥)+𝑦³₀(𝑥)𝑝=0,₂₁𝑝⁰₂∶𝑎𝑦^′₄(𝑥)+𝑏𝑦₄(𝑥)+(1−𝑎)𝑦^′₁(𝑥)+(1−𝑏)𝑦₁(𝑥)𝑝=0,₀₁𝑝²₂∶𝑎𝑦^′₅(𝑥)+𝑏𝑦₅(𝑥)𝑝=0,₀₁𝑝⁰₂∶𝑎𝑦^′₀(𝑥)+𝑏𝑦₀(𝑥)𝑝+𝑐−1=0,₁₁𝑝⁰₂∶𝑎𝑦^′₁(𝑥)+𝑏𝑦₁(𝑥)+(1−𝑎)𝑦^′₀(𝑥)+(1−𝑏)𝑦₀(𝑥)𝑝−𝑐=0,₀₁𝑝¹₂∶𝑎𝑦^′₂(𝑥)+𝑏𝑦₂(𝑥)𝑝=0,₁₁𝑝¹₂∶𝑎𝑦^′₃(𝑥)+𝑏𝑦₃(𝑥)+(1−𝑎)𝑦^′₂(𝑥)+(1−𝑏)𝑦₂(𝑥)−4𝑦²₀(𝑥)𝑝=0,₂₁𝑝⁰₂∶𝑎𝑦^′₄(𝑥)+𝑏𝑦₄(𝑥)+(1−𝑎)𝑦^′₁(𝑥)+(1−𝑏)𝑦₁(𝑥)𝑝=0,₀₁𝑝²₂∶𝑎𝑦^′₅(𝑥)+𝑏𝑦₅(𝑥)=0. (5.22)(5.23)(5.24)(5.25) respectively.

In all cases, solution for 𝑦(𝑥)=𝑦₀(𝑥) and 𝑦(𝑥)=𝑦₀(𝑥)+𝑦₁(𝑥) are𝑦(𝑥)=236−44−15exp13𝑥,𝑦(𝑥)=17580−153269exp−17𝑥−1435−18exp17𝑥14𝑥, (5.26)(5.27) respectively.

For the 𝑁₁,𝑁₂,𝑁₃, and 𝑁₄ cases, the results for 𝑦(𝑥)=𝑦₀(𝑥)+𝑦₁(𝑥)+𝑦₂(𝑥)+𝑦₃(𝑥) are𝑦(𝑥)=68579−1810486964−34529exp24𝑥−191199−488exp24𝑥𝑥+319−448exp1225𝑥−323104−277exp48𝑥,𝑦(𝑥)=19259257−68476316829−130380exp14𝑥−115557−2160exp14𝑥𝑥−911−19exp28𝑥+311−176exp42𝑥,𝑦(𝑥)=11758584−199745157979−46956exp41𝑥−36181−140exp41𝑥𝑥−536−127exp41𝑥+1240−79exp41𝑥,𝑦(𝑥)=183863−10201322−525exp24𝑥−1922697−9238exp24𝑥19𝑥−17−46exp48𝑥,19 (5.28)(5.29)(5.30)(5.31) respectively.

6. Result and Discussion

The graph showing the exact solution for (5.3) (cross) against the approximate solution for (5.11) (solid line), obtained by the standard HPM, is shown in Figure 2. In this figure can be seen a good approximation to the real behaviour of the exact solution. Also, the rest of the located solutions using the proposed methods in this work accomplished similar results to the presented in Figure 2. Therefore, in Table 1 is shown the numerical comparison between the exact result for (5.3) using numerical methods against the results obtained by the HPM standard method, NDHPM, and MHPM. Also, the average absolute relative error (A.A.R.E) is calculated for the five samples taken between the exact numerical solution and the obtained approximations. In all cases the results showed good accuracy (See Table 1). Nevertheless, the approximate solutions (5.15), (5.26), and (5.27), turned to be a less complex mathematical expression compared to the solutions from the basic HPM.

Table 1 

Comparison between the obtained solution against standard HPM, NHPM, and MHPM.

Figure 2 

Exact solution for (5.3) (cross) and its approximate solution (5.11) (solid line).

We suggested a two modified versions NDHPM y MHPM of standard HPM inspired in HCMs. Also, the convergence of the NDHPM and MHPM were analysed; it was determined that those methods are capable to converge according to the Banach’s fixed point theorem. Then, we applied those modified HPMs for solving the CICR model obtaining acceptable results (see Table 1). Additionally, we showed that it is possible to find different solutions to the ones obtained by the standard HPM, with the added value that for the CICR model was possible to find solutions with less exponential terms ((5.15), (5.26), and (5.27)). Approximations containing few exponential terms may considerably reduce computing time, for the case of extensive simulations. Future research should be done to develop the criteria to establish how to place the parameter within the nonlinear part of the differential equation.

For the case of the HPMs just one solution is searched, although is known that there are nonlinear differential equations with multiple critical points that represent different solutions that depends on the initial conditions. Therefore, one possible research line for the HPMs would be to locate multiple solutions at the same time in analogy with Figure 1.

Homotopy methods are a well-established research area with applications, in particular, an applied branch of such methods are the homotopy continuation methods [38–45, 45–52, 60–73]. Therefore, it is possible to research on the properties, characteristics, or variations in such methods in order to be reused focusing on the modification of the HPM.

7. Conclusions

This work presented a series of proposals to create two modified HPMs: NDHPM and MHPM. Both methods were tested by obtaining, successfully, the approximate solution of the CICR model. The main advantage of both methods is the distribution of the nonlinear terms between iterations of the HPM, could generate solutions with good accuracy without high mathematical complexity, compared to the standard HPM. As future work for this paper, other nonlinear equations should be solved based on the proposals made in this paper.