# Modified HPMs Inspired by Homotopy Continuation Methods

Authors: Héctor Vázquez-Leal; Uriel Filobello-Niño; Roberto  							Castañeda-Sheissa; Luis Hernández-Martínez; Arturo Sarmiento-Reyes
Journal: Mathematical Problems in Engineering
 (2012)

---

## 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.


---

## Body

## 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 parameterp, which takes values from 0 to 1. When parameter p=0, the equation usually reduces to a simple, or trivial, equation to solve. Then p gradually is increase to 1, producing a sequence of deformations, where every solution is close to the last one. Eventually at p=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 Section2, 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(2.1)A(u)-f(r)=0,wherer∈Ω.
With the boundary condition(2.2)B(u,∂u∂η),wherer∈Γ,
where A is a general differential operator, f(r) is a known analytic function, B is a boundary operator, and Γ is the boundary of the domain Ω. The A operator, generally, can be divided into two operators, L and N, where L is the linear operator and N is the nonlinear operator. Hence, (2.1) can be rewritten as(2.3)L(u)+N(u)-f(r)=0.
Now, the homotopy function is(2.4)H(v,p)=(1-p)[L(v)-L(u0)]+p(L(v)+N(v)-f(r))=0wherep∈[0,1],
where u0 is the initial approximation of (2.3) which satisfies the boundary conditions and p is known as the perturbation homotopy parameter. Analysing (2.4) can be concluded that(2.5)H(v,0)=L(v)-L(u0),H(v,1)=L(v)+N(v)-f(r).
We assume that the solution of (2.4) can be written as a power series of p [1](2.6)v=p0v0+p1v1+p2v2+⋯.
Adjusting p=1 results that the approximate solution for (2.1) is(2.7)u=lim⁡p→1v=v0+v1+v2+⋯.The series (2.7) is convergent on most cases, nevertheless, the convergence depends on the nonlinear operator A(v) [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(2.8)f(x)=0,wheref:∈Rn⟶Rn,
where x represents the variables of problem and n 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(2.9)H(f(x),p)=0,
where f(x) is the equilibrium equation and p is the homotopy parameter just like the HPM. For this case, (2.9) represents any homotopy formulation that fulfils the following basic conditions:(i)
forp=0 the solution of H-1(0) is known or easy to obtain using numerical methods. This point is known as Homotopy’s initial point (x0),(ii)
forp=1,H(f(x),1)=f(x). Means that at p=1 all solutions for f(x) are located,(iii)
the path forH-1(0) is a continuous function of p within the range of 0≤p≤1.The homotopy path is the solution set ofH-1(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(2.10)H(f(x),p)=f(x)-(1-p)f(x0)=0,
where f(x) is the equilibrium equation, x0 is the initial point of the homotopy, and p is the homotopy parameter.Ifp=0,(2.11)H(f(x),0)=f(x)-f(x0)=0,
where x0 represents the initial point. It is simple to obtain or is chosen arbitrarily.Ifp=1(2.12)H(f(x),1)=f(x)=0,
in this case the solution of H is, precisely, the desired solution of the NA equation f(x). That is, a continuous path is created between p=0 and p=1 which can be described as in Figure 1. In case that there is more than one solution: x1*,x2*,x3*, and x4*, for f(x), there is a chance that those are connected to the same path.Figure 1
Homotopy path.


## 2.1. Basic Idea of HPMs

In the HPM is considered that a nonlinear differential equation can be expressed as(2.1)A(u)-f(r)=0,wherer∈Ω.
With the boundary condition(2.2)B(u,∂u∂η),wherer∈Γ,
where A is a general differential operator, f(r) is a known analytic function, B is a boundary operator, and Γ is the boundary of the domain Ω. The A operator, generally, can be divided into two operators, L and N, where L is the linear operator and N is the nonlinear operator. Hence, (2.1) can be rewritten as(2.3)L(u)+N(u)-f(r)=0.
Now, the homotopy function is(2.4)H(v,p)=(1-p)[L(v)-L(u0)]+p(L(v)+N(v)-f(r))=0wherep∈[0,1],
where u0 is the initial approximation of (2.3) which satisfies the boundary conditions and p is known as the perturbation homotopy parameter. Analysing (2.4) can be concluded that(2.5)H(v,0)=L(v)-L(u0),H(v,1)=L(v)+N(v)-f(r).
We assume that the solution of (2.4) can be written as a power series of p [1](2.6)v=p0v0+p1v1+p2v2+⋯.
Adjusting p=1 results that the approximate solution for (2.1) is(2.7)u=lim⁡p→1v=v0+v1+v2+⋯.The series (2.7) is convergent on most cases, nevertheless, the convergence depends on the nonlinear operator A(v) [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(2.8)f(x)=0,wheref:∈Rn⟶Rn,
where x represents the variables of problem and n 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(2.9)H(f(x),p)=0,
where f(x) is the equilibrium equation and p is the homotopy parameter just like the HPM. For this case, (2.9) represents any homotopy formulation that fulfils the following basic conditions:(i)
forp=0 the solution of H-1(0) is known or easy to obtain using numerical methods. This point is known as Homotopy’s initial point (x0),(ii)
forp=1,H(f(x),1)=f(x). Means that at p=1 all solutions for f(x) are located,(iii)
the path forH-1(0) is a continuous function of p within the range of 0≤p≤1.The homotopy path is the solution set ofH-1(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(2.10)H(f(x),p)=f(x)-(1-p)f(x0)=0,
where f(x) is the equilibrium equation, x0 is the initial point of the homotopy, and p is the homotopy parameter.Ifp=0,(2.11)H(f(x),0)=f(x)-f(x0)=0,
where x0 represents the initial point. It is simple to obtain or is chosen arbitrarily.Ifp=1(2.12)H(f(x),1)=f(x)=0,
in this case the solution of H is, precisely, the desired solution of the NA equation f(x). That is, a continuous path is created between p=0 and p=1 which can be described as in Figure 1. In case that there is more than one solution: x1*,x2*,x3*, and x4*, for f(x), 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, multiplyingp 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 p=0 to p=1, creating a smooth path for the homotopy. An example of this type of homotopy can be formulated as(3.1)H(f(x),p)=(1-p)(x-x0)+pf(x,p)=0,
where f(x) is the NA equation, x0 is the initial point of the homotopy, and p is the homotopy parameter. It can be seen that when p=0, then(3.2)H(f(x),0)=x-x0=0,
when p=1, then(3.3)H(f(x),1)=f(x,1)=f(x)=0,
where for p between zero and one, the homotopy parameter reduce the nonlinear effect from nullifying it at p=0 until giving them their original state at p=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 ofN and f among sequential iterations of the HPM.(3.4)H(v,p)=(1-p)[L(v)-L(u0)]+p(L(v)+N(v,p)-f(r,p))=0,wherep∈[0,1].The homotopy function (3.4) is essentially the same as (2.4), except for the nonlinear operator N and function f which contain embedded the homotopy parameter. We propose, first, to subdivide N and f into a series of terms and multiply pi by the most nonlinear terms, where i is an integer number greater or equal to zero. The i power is selected according to how much displacement is desired in the iterations for the corresponding nonlinear term of N or f.As mentioned before, embedding the homotopy parameterp 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)v=p0v0+p1v1+p2v2+⋯,
adjusting p=1 turns out that the approximate solution of (2.1) is(3.6)u=lim⁡p→1v=v0+v1+v2+⋯.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 homotopyv(r,p):Ω×[0,1]→ℜn which satisfies(3.7)H(v,p)=L(v)-L(u0)+pL(u0)+p[N(v,p)-f(r,p)]=0,
where u0 is the initial approximation of (2.1).Let one writes (3.7) in the following form:(3.8)L(v)=L(u0)+p[f(r,p)-N(v,p)-L(u0)]=0.
Applying the inverse operator, L-1, to both sides of (3.8) we obtain(3.9)v=u0+p[L-1f(r,p)-L-1N(v,p)-u0].
Suppose that(3.10)v=∑i=0∞pivi,
substituting (3.10) in the right-hand side of (3.9) in the following form:(3.11)v=u0+p[L-1f(r,p)-(L-1N)[∑i=0∞pivi]-u0].
If p→1, the exact solution may be obtained by using(3.12)u=lim⁡p→1v=L-1f(r)-(L-1N)[∑i=0∞vi]=L-1f(r)-[∑i=0∞(L-1N)(vi)].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.


## 3.1. Convergence

We construct a homotopyv(r,p):Ω×[0,1]→ℜn which satisfies(3.7)H(v,p)=L(v)-L(u0)+pL(u0)+p[N(v,p)-f(r,p)]=0,
where u0 is the initial approximation of (2.1).Let one writes (3.7) in the following form:(3.8)L(v)=L(u0)+p[f(r,p)-N(v,p)-L(u0)]=0.
Applying the inverse operator, L-1, to both sides of (3.8) we obtain(3.9)v=u0+p[L-1f(r,p)-L-1N(v,p)-u0].
Suppose that(3.10)v=∑i=0∞pivi,
substituting (3.10) in the right-hand side of (3.9) in the following form:(3.11)v=u0+p[L-1f(r,p)-(L-1N)[∑i=0∞pivi]-u0].
If p→1, the exact solution may be obtained by using(3.12)u=lim⁡p→1v=L-1f(r)-(L-1N)[∑i=0∞vi]=L-1f(r)-[∑i=0∞(L-1N)(vi)].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(4.1)H(v,p)=(1-p1)[L(v)-L(u0)]+p1(L(v)+N(v,p2)-f(r,p2))=0,
where p1∈[0,1] and p2∈[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 p1 and p2; the first homotopy parameter p1 formulates the basic homotopy (2.4) and the second homotopy parameter p2 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 N and f in a set of terms and multiply p2i by the most nonlinear terms, where i is an integer number greater or equal to zero. The i power is selected according to how much displacement is desired in the iterations the corresponding term N or f.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)v=v0+v1p1+v2p2+v3p1p2+v4p12+v5p22+v6p1p22+v7p12p2+⋯,
adjusting p=1 results that the approximate solution for (2.1) is(4.3)u=lim⁡p1→1p2→1v=v0+v1+v2+⋯.This method will be designated as multiparameter HPM (MHPM).


### 4.1. Convergence

We construct a homotopyv(r,p1,p2):Ω×[0,1]→ℜn which satisfies(4.4)H(v,p)=L(v)-L(u0)+p1L(u0)+p1[N(v,p2)-f(r,p2)]=0,
where u0 is the initial approximation of (2.1).Let one writes (4.4) in the following form:(4.5)L(v)=L(u0)+p1[f(r,p2)-N(v,p2)-L(u0)]=0.
Applying the inverse operator L-1 to both sides of (4.5), we obtain(4.6)v=u0+p1[L-1f(r,p2)-L-1N(v,p2)-u0].
Substituting (4.2) in the right-hand side of (4.6) in the following form:(4.7)v=u0+p1[L-1f(r,p1,p2)-(L-1N)[v0+v1p1+v2p2+v3p1p2+v4p12+v5p22+v6p1p22+⋯]-u0].
If p1→1 and p2→1, the exact solution may be obtained by using(4.8)u=lim⁡p1→1p2→1v=L-1f(r)-(L-1N)[∑i=0∞vi]=L-1f(r)-[∑i=0∞(L-1N)(vi)].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.


## 4.1. Convergence

We construct a homotopyv(r,p1,p2):Ω×[0,1]→ℜn which satisfies(4.4)H(v,p)=L(v)-L(u0)+p1L(u0)+p1[N(v,p2)-f(r,p2)]=0,
where u0 is the initial approximation of (2.1).Let one writes (4.4) in the following form:(4.5)L(v)=L(u0)+p1[f(r,p2)-N(v,p2)-L(u0)]=0.
Applying the inverse operator L-1 to both sides of (4.5), we obtain(4.6)v=u0+p1[L-1f(r,p2)-L-1N(v,p2)-u0].
Substituting (4.2) in the right-hand side of (4.6) in the following form:(4.7)v=u0+p1[L-1f(r,p1,p2)-(L-1N)[v0+v1p1+v2p2+v3p1p2+v4p12+v5p22+v6p1p22+⋯]-u0].
If p1→1 and p2→1, the exact solution may be obtained by using(4.8)u=lim⁡p1→1p2→1v=L-1f(r)-(L-1N)[∑i=0∞vi]=L-1f(r)-[∑i=0∞(L-1N)(vi)].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]:(5.1)y′(x)=A+k1y2(1+y2)-k2y,
where x represents time, y represents calcium concentration, k1 and k2 are constants, and A 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(5.2)y′(x)+y2y′(x)+k2y3-(A+k1)y2+k2y-A=0.To exemplify, we assign the following valuesA=1,k1=3,k2=1, and the initial condition y(0)=0.9. Therefore, (5.2) is updated becoming(5.3)y′(x)+y2y′(x)+y3-4y2+y-1=0,wherey(0)=0.9.Equation (5.3) represents the transition from a state given by the initial condition y(0)=0.9 until a steady state for the calcium concentration. Now, the linear part is separated(5.4)L(y)=y′(x)+y-1,
and its nonlinear part(5.5)N(y)=y2(x)y′(x)+y3-4y2.Adjustment constants are added to (5.4), the result is(5.6)Lc(y)=ay′(x)+by+c-1,
here a,b, and c are adjustment constants.Now, we establish the homotopy equation(5.7)(1-p)Lc(y)+p(L(y)+N(y))=0,
where p is the homotopy parameter.We suppose that solution for (5.7) has the form(5.8)y(x)=y0(x)+py1(x)+p2y2(x)+⋯.Adjustingp=1, we can obtain an approximate solution(5.9)y(x)=y0(x)+y1(x)+y2(x)+⋯.Substituting (5.8) into (5.7) and equaling terms having potentials in the same order as p, it can be solved for y0(x),y1(x),y2(x), and so on (in order to fulfil initial conditions from y(0)=0.9, it follows that y0(0)=0.9,y1(0)=0,y2(0)=0, and so on). Therefore, the following differential equations are established:(5.10)ay0′(x)+by0(x)+c-1=0,ay1′(x)+by1(x)+(y02(x)-a+1)y0′(x)+y03(x)-4y02(x)+(1-b)y0(x)-c=0,⋮
where terms y2 and subsequent are ignored because the first two terms (y0(x) and y1(x)) contribute to obtain a good accurate solution. In a similar process presented in [20, 59], we use NonlinearFit (given a total of k-samples from the exact model, the NonlinearFit command finds values of the approximate model parameters such that the sum of the squared k-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:(5.11)y(x)=(-7397)exp⁡(-8633x)+(11536)exp⁡(-3319x)+5314+(-29355+3753x)exp⁡(-3338x)


### 5.2. Using the NDHPM

The NDHPMs establish that(5.12)(1-p)Lc(y)+p(L(y)+Ni(y,p)-f(r))=0,
where p is the homotopy parameter, Lc is established in (5.6), f(r)=0, and the nonlinear part Ni(y,p) (i=1,2,3,4) is one of the following representative cases:(5.13)N1(y,p)=(y2y′(x)+y3-4y2)p,N2(y,p)=y2y′(x)+y3p-4y2p,N3(y,p)=y2y′(x)p+y3-4y2p,N4(y,p)=y2y′(x)p+y3p-4y2.Consider that whenlimp→1Ni=N.Now, the exact basic steps for the HPM method are followed to obtain the following differential equation system for casesN1,N2,N3, and N4:(5.14)p0:ay0′(x)+by0(x)+c-1=0,p1:ay1′(x)+by1(x)+(-a+1)y0′(x)+(1-b)y0(x)-c=0,⋮p0:ay0′(x)+by0(x)+c-1=0,p1:ay1′(x)+by1(x)+(y02(x)-a+1)y0′(x)+(1-b)y0(x)-c=0,⋮p0:ay0′(x)+by0(x)+c-1=0,p1:ay1′(x)+by1(x)+(-a+1)y0′(x)+y03(x)+(1-b)y0(x)-c=0,⋮p0:ay0′(x)+by0(x)+c-1=0,p1:ay1′(x)+by1(x)+(-a+1)y0′(x)-4y02(x)+(1-b)y0(x)-c=0,⋮
with initial conditions y0(0)=0.9,y1(0)=0,y2(0)=0,…. This is done to fulfil initial conditions of (5.3).In this case, the termsy2 and subsequent are ignored because the first two terms (y0(x) and y1(x)) provide good accuracy to the solution. Then, we obtain the following accurate approximations for all cases N1,N2,N3 and N4:(5.15)y(x)=5314+(-269+1629x)exp⁡(-1529x),(5.16)y(x)=(11264)exp⁡(-4312x)-(11118)exp⁡(-4318x)+195+(-4516-3719x)exp⁡(-3731x),(5.17)y(x)=(554)exp⁡(-4311x)-(8571)exp⁡(-9938x)+195+(-7944-9429x)exp⁡(-4333x),(5.18)y(x)=-(1341)exp⁡(-3916x)+195+(-3112-2913x)exp⁡(-3932x),
respectively.


### 5.3. Using the Multiparameter HPM

The homotopy equation is established according to the basic guidelines for the multiparameter HPM:(5.19)(1-p1)Lc(y)+p1(L(y)+Ni(y,p2))=0,
where p1 is the first homotopy parameter, Lc  is defined in (5.6), and we propose the nonlinearoperator Ni(y,p2)(i=1,2,3,4) as the following representative case:(5.20)N1(y,p2)=(y2y′(x)+y3-4y2)p2,N2(y,p2)=y2y′(x)p2+y3p22-4y2p2,N3(y,p2)=y2y′(x)p2+y3p2-4y2p22,N4(y,p2)=y2y′(x)p22+y3p22-4y2p2.From (4.2) we obtain(5.21)y(x)=y0(x)+y1(x)p1+y2(x)p2+y3(x)p1p2+y4(x)p12+y5(x)p22+⋯.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 y0(0)=0.9,y1(0)=0,y2(0)=0,y3(0)=0,y4(0)=0, and y5(0)=0) for cases N1,N2,N3, and N4, we obtain(5.22)p10p20:ay0′(x)+by0(x)+c-1=0,p11p20:ay1′(x)+by1(x)+(1-a)y0′(x)+(1-b)y0(x)-c=0,p10p21:ay2′(x)+by2(x)=0,p11p21:ay3′(x)+by3(x)+(1-a)y2′(x)+(1-b)y2(x)+y02(x)y0′(x)-4y02(x)+y03(x)=0,p12p20:ay4′(x)+by4(x)+(1-a)y1′(x)+(1-b)y1(x)=0,p10p22:ay5′(x)+by5(x)=0,(5.23)p10p20:ay0′(x)+by0(x)+c-1=0,p11p20:ay1′(x)+by1(x)+(1-a)y0′(x)+(1-b)y0(x)-c=0,p10p21:ay2′(x)+by2(x)=0,p11p21:ay3′(x)+by3(x)+(1-a)y2′(x)+(1-b)y2(x)+y02(x)y0′(x)-4y02(x)=0,p12p20:ay4′(x)+by4(x)+(1-a)y1′(x)+(1-b)y1(x)=0,p10p22:ay5′(x)+by5(x)=0,(5.24)p10p20:ay0′(x)+by0(x)+c-1=0,p11p20:ay1′(x)+by1(x)+(1-a)y0′(x)+(1-b)y0(x)-c=0,p10p21:ay2′(x)+by2(x)=0,p11p21:ay3′(x)+by3(x)+(1-a)y2′(x)+(1-b)y2(x)+y02(x)y0′(x)+y03(x)=0,p12p20:ay4′(x)+by4(x)+(1-a)y1′(x)+(1-b)y1(x)=0,p10p22:ay5′(x)+by5(x)=0,(5.25)p10p20:ay0′(x)+by0(x)+c-1=0,p11p20:ay1′(x)+by1(x)+(1-a)y0′(x)+(1-b)y0(x)-c=0,p10p21:ay2′(x)+by2(x)=0,p11p21:ay3′(x)+by3(x)+(1-a)y2′(x)+(1-b)y2(x)-4y02(x)=0,p12p20:ay4′(x)+by4(x)+(1-a)y1′(x)+(1-b)y1(x)=0,p10p22:ay5′(x)+by5(x)=0.
respectively.In all cases, solution fory(x)=y0(x) and y(x)=y0(x)+y1(x) are(5.26)y(x)=236-(4415)exp⁡(-1317x),(5.27)y(x)=580153-(269)exp⁡(-1714x)-(3518)exp⁡(-1714x)x,
respectively.For theN1,N2,N3, and N4 cases, the results for y(x)=y0(x)+y1(x)+y2(x)+y3(x) are(5.28)y(x)=6857918104-(8696434529)exp⁡(-2419x)-(1199488)exp⁡(-2419x)x+(3448)exp⁡(-1225323x)-(104277)exp⁡(-4819x),(5.29)y(x)=25925768476-(316829130380)exp⁡(-1411x)-(55572160)exp⁡(-1411x)x-(919)exp⁡(-2811x)+(3176)exp⁡(-4211x),(5.30)y(x)=758584199745-(15797946956)exp⁡(-4136x)-(181140)exp⁡(-4136x)x-(5127)exp⁡(-4112x)+(4079)exp⁡(-4118x),(5.31)y(x)=38631020-(1322525)exp⁡(-2419x)-(226979238)exp⁡(-2419x)x-(1746)exp⁡(-4819x),
respectively.


## 5.1. Using Basic HPM

In order to analyse (5.1), it should be reordered(5.2)y′(x)+y2y′(x)+k2y3-(A+k1)y2+k2y-A=0.To exemplify, we assign the following valuesA=1,k1=3,k2=1, and the initial condition y(0)=0.9. Therefore, (5.2) is updated becoming(5.3)y′(x)+y2y′(x)+y3-4y2+y-1=0,wherey(0)=0.9.Equation (5.3) represents the transition from a state given by the initial condition y(0)=0.9 until a steady state for the calcium concentration. Now, the linear part is separated(5.4)L(y)=y′(x)+y-1,
and its nonlinear part(5.5)N(y)=y2(x)y′(x)+y3-4y2.Adjustment constants are added to (5.4), the result is(5.6)Lc(y)=ay′(x)+by+c-1,
here a,b, and c are adjustment constants.Now, we establish the homotopy equation(5.7)(1-p)Lc(y)+p(L(y)+N(y))=0,
where p is the homotopy parameter.We suppose that solution for (5.7) has the form(5.8)y(x)=y0(x)+py1(x)+p2y2(x)+⋯.Adjustingp=1, we can obtain an approximate solution(5.9)y(x)=y0(x)+y1(x)+y2(x)+⋯.Substituting (5.8) into (5.7) and equaling terms having potentials in the same order as p, it can be solved for y0(x),y1(x),y2(x), and so on (in order to fulfil initial conditions from y(0)=0.9, it follows that y0(0)=0.9,y1(0)=0,y2(0)=0, and so on). Therefore, the following differential equations are established:(5.10)ay0′(x)+by0(x)+c-1=0,ay1′(x)+by1(x)+(y02(x)-a+1)y0′(x)+y03(x)-4y02(x)+(1-b)y0(x)-c=0,⋮
where terms y2 and subsequent are ignored because the first two terms (y0(x) and y1(x)) contribute to obtain a good accurate solution. In a similar process presented in [20, 59], we use NonlinearFit (given a total of k-samples from the exact model, the NonlinearFit command finds values of the approximate model parameters such that the sum of the squared k-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:(5.11)y(x)=(-7397)exp⁡(-8633x)+(11536)exp⁡(-3319x)+5314+(-29355+3753x)exp⁡(-3338x)


## 5.2. Using the NDHPM

The NDHPMs establish that(5.12)(1-p)Lc(y)+p(L(y)+Ni(y,p)-f(r))=0,
where p is the homotopy parameter, Lc is established in (5.6), f(r)=0, and the nonlinear part Ni(y,p) (i=1,2,3,4) is one of the following representative cases:(5.13)N1(y,p)=(y2y′(x)+y3-4y2)p,N2(y,p)=y2y′(x)+y3p-4y2p,N3(y,p)=y2y′(x)p+y3-4y2p,N4(y,p)=y2y′(x)p+y3p-4y2.Consider that whenlimp→1Ni=N.Now, the exact basic steps for the HPM method are followed to obtain the following differential equation system for casesN1,N2,N3, and N4:(5.14)p0:ay0′(x)+by0(x)+c-1=0,p1:ay1′(x)+by1(x)+(-a+1)y0′(x)+(1-b)y0(x)-c=0,⋮p0:ay0′(x)+by0(x)+c-1=0,p1:ay1′(x)+by1(x)+(y02(x)-a+1)y0′(x)+(1-b)y0(x)-c=0,⋮p0:ay0′(x)+by0(x)+c-1=0,p1:ay1′(x)+by1(x)+(-a+1)y0′(x)+y03(x)+(1-b)y0(x)-c=0,⋮p0:ay0′(x)+by0(x)+c-1=0,p1:ay1′(x)+by1(x)+(-a+1)y0′(x)-4y02(x)+(1-b)y0(x)-c=0,⋮
with initial conditions y0(0)=0.9,y1(0)=0,y2(0)=0,…. This is done to fulfil initial conditions of (5.3).In this case, the termsy2 and subsequent are ignored because the first two terms (y0(x) and y1(x)) provide good accuracy to the solution. Then, we obtain the following accurate approximations for all cases N1,N2,N3 and N4:(5.15)y(x)=5314+(-269+1629x)exp⁡(-1529x),(5.16)y(x)=(11264)exp⁡(-4312x)-(11118)exp⁡(-4318x)+195+(-4516-3719x)exp⁡(-3731x),(5.17)y(x)=(554)exp⁡(-4311x)-(8571)exp⁡(-9938x)+195+(-7944-9429x)exp⁡(-4333x),(5.18)y(x)=-(1341)exp⁡(-3916x)+195+(-3112-2913x)exp⁡(-3932x),
respectively.


## 5.3. Using the Multiparameter HPM

The homotopy equation is established according to the basic guidelines for the multiparameter HPM:(5.19)(1-p1)Lc(y)+p1(L(y)+Ni(y,p2))=0,
where p1 is the first homotopy parameter, Lc  is defined in (5.6), and we propose the nonlinearoperator Ni(y,p2)(i=1,2,3,4) as the following representative case:(5.20)N1(y,p2)=(y2y′(x)+y3-4y2)p2,N2(y,p2)=y2y′(x)p2+y3p22-4y2p2,N3(y,p2)=y2y′(x)p2+y3p2-4y2p22,N4(y,p2)=y2y′(x)p22+y3p22-4y2p2.From (4.2) we obtain(5.21)y(x)=y0(x)+y1(x)p1+y2(x)p2+y3(x)p1p2+y4(x)p12+y5(x)p22+⋯.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 y0(0)=0.9,y1(0)=0,y2(0)=0,y3(0)=0,y4(0)=0, and y5(0)=0) for cases N1,N2,N3, and N4, we obtain(5.22)p10p20:ay0′(x)+by0(x)+c-1=0,p11p20:ay1′(x)+by1(x)+(1-a)y0′(x)+(1-b)y0(x)-c=0,p10p21:ay2′(x)+by2(x)=0,p11p21:ay3′(x)+by3(x)+(1-a)y2′(x)+(1-b)y2(x)+y02(x)y0′(x)-4y02(x)+y03(x)=0,p12p20:ay4′(x)+by4(x)+(1-a)y1′(x)+(1-b)y1(x)=0,p10p22:ay5′(x)+by5(x)=0,(5.23)p10p20:ay0′(x)+by0(x)+c-1=0,p11p20:ay1′(x)+by1(x)+(1-a)y0′(x)+(1-b)y0(x)-c=0,p10p21:ay2′(x)+by2(x)=0,p11p21:ay3′(x)+by3(x)+(1-a)y2′(x)+(1-b)y2(x)+y02(x)y0′(x)-4y02(x)=0,p12p20:ay4′(x)+by4(x)+(1-a)y1′(x)+(1-b)y1(x)=0,p10p22:ay5′(x)+by5(x)=0,(5.24)p10p20:ay0′(x)+by0(x)+c-1=0,p11p20:ay1′(x)+by1(x)+(1-a)y0′(x)+(1-b)y0(x)-c=0,p10p21:ay2′(x)+by2(x)=0,p11p21:ay3′(x)+by3(x)+(1-a)y2′(x)+(1-b)y2(x)+y02(x)y0′(x)+y03(x)=0,p12p20:ay4′(x)+by4(x)+(1-a)y1′(x)+(1-b)y1(x)=0,p10p22:ay5′(x)+by5(x)=0,(5.25)p10p20:ay0′(x)+by0(x)+c-1=0,p11p20:ay1′(x)+by1(x)+(1-a)y0′(x)+(1-b)y0(x)-c=0,p10p21:ay2′(x)+by2(x)=0,p11p21:ay3′(x)+by3(x)+(1-a)y2′(x)+(1-b)y2(x)-4y02(x)=0,p12p20:ay4′(x)+by4(x)+(1-a)y1′(x)+(1-b)y1(x)=0,p10p22:ay5′(x)+by5(x)=0.
respectively.In all cases, solution fory(x)=y0(x) and y(x)=y0(x)+y1(x) are(5.26)y(x)=236-(4415)exp⁡(-1317x),(5.27)y(x)=580153-(269)exp⁡(-1714x)-(3518)exp⁡(-1714x)x,
respectively.For theN1,N2,N3, and N4 cases, the results for y(x)=y0(x)+y1(x)+y2(x)+y3(x) are(5.28)y(x)=6857918104-(8696434529)exp⁡(-2419x)-(1199488)exp⁡(-2419x)x+(3448)exp⁡(-1225323x)-(104277)exp⁡(-4819x),(5.29)y(x)=25925768476-(316829130380)exp⁡(-1411x)-(55572160)exp⁡(-1411x)x-(919)exp⁡(-2811x)+(3176)exp⁡(-4211x),(5.30)y(x)=758584199745-(15797946956)exp⁡(-4136x)-(181140)exp⁡(-4136x)x-(5127)exp⁡(-4112x)+(4079)exp⁡(-4118x),(5.31)y(x)=38631020-(1322525)exp⁡(-2419x)-(226979238)exp⁡(-2419x)x-(1746)exp⁡(-4819x),
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.
Equationx=0x=2x=4x=6x=8A.A.R.EExact (3.12)0.900003.192383.703533.789423.80354Standard HPM  (5.11)0.900303.188493.710183.779593.785962.1E-3NDHPM (5.15)0.896833.151153.699573.804643.810054.6E-3NDHPM (5.16)0.895073.182853.710463.788753.798692.4E-3NDHPM (5.17)0.899953.182363.719513.791453.799181.8E-3NDHPM (5.18)0.899593.182003.712133.789353.798811.5E-3MHPM (5.26)0.900003.197793.695633.803493.826872.7E-3MHPM (5.27)0.901963.193303.707943.780873.789741.90E-3MHPM (5.28)0.900723.191413.709123.779243.787151.92E-3MHPM (5.29)0.899423.188993.707823.777483.785232.2E-3MHPM (5.30)0.900323.193083.708123.785783.796250.9E-3MHPM (5.31)0.899593.190693.708323.778433.786351.94E-3Figure 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 Table1). 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 Figure1.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.



---

*Source: 309123-2012-03-13.xml*