See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/319407345

NEW ADD-IN FOR ENGINEERING SIMULATIONS IN SPREADSHEETS USING

HOMOTOPY CONTINUATION METHODS

Article  in  Revista mexicana de ingeniería química · December 2017

CITATIONS READS
0 194

3 authors, including:

Hugo Jiménez-Islas Rolando Barrera

Instituto Tecnológico de Celaya University of Antioquia
118 PUBLICATIONS   605 CITATIONS    40 PUBLICATIONS   89 CITATIONS   

SEE PROFILE SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Rheological modelling View project

Optimization of hyperspherical path tracking for homotopic continuation methods View project

All content following this page was uploaded by Hugo Jiménez-Islas on 13 October 2017.

The user has requested enhancement of the downloaded file.

 Vol. 16, No. 3 (2017) 1013-1029
Revista Mexicana de Ingeniería Química
CONTENIDO
NEW ADD-IN FOR ENGINEERING SIMULATIONS IN SPREADSHEETS USING
Volumen 8, número 3, 2009
HOMOTOPY / Volume 8, numberMETHODS
CONTINUATION 3, 2009

NUEVO COMPLEMENTO PARA SIMULACIONES DE INGENIERÍA EN HOJAS DE

CÁLCULO USANDO
213 Derivation MÉTODOS
and application DE CONTINUIDAD
of the Stefan-Maxwell equations HOMOTÓPICA
L. Dı́az-Montes1 , R. Barrera-Zapata1 * , H. Jiménez-Islas2
1 Grupo CERES, Departamento de ingenierı́a
(Desarrollo quı́mica,
y aplicación Facultad de
de las ecuaciones de Stefan-Maxwell)
ingenierı́a, Universidad de Antioquia UdeA, Calle 70 No.
Stephen Whitaker 52-21, Medellı́n, Colombia
2 Departamento de Ingenierı́a Bioquı́mica. Tecnológico Nacional de México, Instituto Tecnológico de Celaya, Ave. Tecnológico
y A. Garcı́a Cubas s/n. Celaya, Gto. CP 38010. México
Biotecnología / Received November 15, 2016; Accepted July 2, 2017
Biotechnology

Abstract 245 Modelado de la biodegradación en biorreactores de lodos de hidrocarburos totales del petróleo
Solving models represented by non-linear
intemperizados equations
en suelos is a common situation in engineering and can be a challenging task.
y sedimentos
Homotopy Continuation Methods have shown globally convergent behavior, capable of finding all the possible roots of algebraic
(Biodegradation
systems. Based on these methods, a new modeling
MicrosoftofOffice
sludgeExcel
bioreactors
add-inofwas
totaldeveloped,
petroleum hydrocarbons
SphereSolver,weathering
aiming toinenhance
soil the
ability to use this easily obtainable software
and sediments) as a simulation computing platform for engineering problems. SphereSolver was
applied over different typesS.A.of Medina-Moreno,
engineering models and the non-linear
S. Huerta-Ochoa, equations were solved
C.A. Lucho-Constantino, successfully each
L. Aguilera-Vázquez, time. Selected
A. Jiménez-
problems were solved with SphereSolver and compared with other solution tecniques, e.g., classical Newton’s method, Fixed
González y M. Gutiérrez-Rojas
Point homotopy or Affine Homotopy. It was observed that using SphereSolver the solution for selected problems was improved,
259 Crecimiento, sobrevivencia y adaptación de Bifidobacterium infantis a condiciones ácidas
i.e., the interval of convergence is extended massively and/or additional roots can be found. It can be concluded that SphereSolver
is a new and efficient tool (Growth,
that can survival
be usedand
within Microsoft
adaptation Excel for solving
of Bifidobacterium engineering
infantis problems where other software or
to acidic conditions)
traditional tools could fail or present limitations.
L. Mayorga-Reyes, P. Bustamante-Camilo, A. Gutiérrez-Nava, E. Barranco-Florido y A. Azaola-
Keywords: engineering non-linear models, homotopy continuation methods, excel add-in, hyperspherical path tracking.
Espinosa
Resumen
265 Statistical approach to optimization of ethanol fermentation by Saccharomyces cerevisiae in the
Resolver modelos matemáticos representados por ecuaciones no lineales es común en la ingenierı́a y puede ser una tarea difı́cil.
Los métodos de Continuidad Homotópica
presence hanzeolite
of Valfor® mostrado
NaAser globalmente convergentes y capaces de encontrar todas las posibles
raı́ces de sistemas algebraicos. Con base en este tipo de métodos, se presenta un nuevo complemento para Microsoft Office Excel,
(Optimización estadística de la fermentación etanólica de Saccharomyces cerevisiae en presencia de
SphereSolver, desarrollado con el objetivo de mejorar la capacidad de este software común, como una plataforma de computación
zeolita
de simulación para problemas deValfor® zeolite
ingenierı́a. NaA)
SphereSolver se aplicó sobre diferentes tipos de modelos de ingenierı́a y en todos
los casos, las ecuaciones noG. lineales se resolvieron
Inei-Shizukawa, con éxito. La G.
H. A. Velasco-Bedrán, solución con SphereSolver
F. Gutiérrez-López a los problemas seleccionados se
and H. Hernández-Sánchez
comparó con otras técnicas de solución, por ejemplo, el método clásico de Newton, Homotopı́a de punto fijo y Homotopı́a Afı́n.
Se observó que usandoIngeniería
SphereSolver el intervalo
de procesos convergencia se extiende masivamente y / o se pueden encontrar raı́ces
deengineering
/ Process
adicionales, es decir, se mejoró la solución. Se concluye que SphereSolver es una herramienta eficiente que se puede utilizar
271 Localización de una planta industrial: Revisión crítica y adecuación de los criterios empleados en
en Microsoft Excel para resolver problemas de ingenierı́a donde otros programas o herramientas tradicionales pueden fallar o
presentar limitaciones. esta decisión
Palabras clave: modelos de (Plant site selection:
ingenierı́a Criticalmétodos
no lineales, review and
deadequation
continuidadcriteria used in thiscomplemento
homotópica, decision) de excel, seguimiento
hiperesférico. J.R. Medina, R.L. Romero y G.A. Pérez

1 Introduction like methods are a classical approach used in these

cases, though they are locally convergent, i.e., they
require an initial approximation close enough to the
solution to find the value of the vector that satisfies a
The solution of systems of non-linear equations is a given convergence criterion (Davis, 1984). Homotopy
common practice in engineering, and transcendental Continuation Methods (HCMs) have emerged as an
functions are usually part of these models (Jiménez- alternative strategy, which uses the local convergence
Islas et al., 2014; Rodrı́guez and Niño, 2016). Newton-
* Corresponding author. E-mail: rolando.barrera@udea.edu.co
Tel. +57 +4 2198564

Publicado por la Academia Mexicana de Investigación y Docencia en Ingenierı́a Quı́mica A.C. 1013
 Dı́az-Montes et al./ Revista Mexicana de Ingenierı́a Quı́mica Vol. 16, No. 3 (2017) 1013-1029

properties of such methods in a gradual deformation parameter t that establishes specific conditions to be
scheme (Allgower and Georg, 1990). They have satisfied by the homotopic function, H(t, x), i.e., Eq.
been shown to be accurate, robust, and able to (1) (Cveticanin et al., 2012; Ward and King, 2012).
converge from almost any arbitrary starting point.
 E(x) if t = 0

These techniques have been successfully applied to H(t, x) = 

 F(x) if t = 1 (1)
solve polynomial systems of equations, constrained
and unconstrained optimization problems, and non-
Thus the standard convex homotopic function, Eq.
linear partial differential equations. Therefore, the
(2), is stated (Verschelde and Cools, 1994; Thompson,
applications of HCMs to solve engineering models
2006).
have been numerous, including phase separation and
H(t, x) = (1 − t)E(x) + tF(x) (2)
transport phenomena problems, among others (Jalali
et al,. 2008; Menglong and Zhenxin, 2008; Laiadi and Starting from t = t0 = 0, if a small increment ∆t is
Merzougui, 2012; Yang et al., 2012; Jimenez-Islas et applied over t0 , such as t1 = t0 + ∆t → t1 ≈ t0 →
al., 2013; Torres-Muñoz et al., 2015). H(t1 , x) ≈ H(0, x), it is expected that x( H, 1) ≈ xE ,
The solution of systems of non-linear equations is where H(t1 , xH,1 ) = 0. Once the value of xH1 is
a common practice in engineering, and transcendental obtained, t is slightly changed again, such as t2 =
functions are usually part of these models (Jiménez- t1 + ∆t → t2 ≈ t1 → H(t2 , x) ≈ H(t1 , x), and therefore
Islas et al., 2014; Rodrı́guez and Niño, 2016). Newton- x( H, 2) ≈ x( H, 1), where H(t2 , xH,2 ) = 0. Following this
like methods are a classical approach used in these pattern, an increase in t is made such as tk+1 = tk +∆t →
cases, though they are locally convergent, i.e., they tk+1 ≈ tk → H(tk+1 , x) ≈ H(tk , x), using xH,k as the
require an initial approximation close enough to the initial trial in a locally convergent numerical method
solution to find the value of the vector that satisfies a to find the vector xH,k+1 . Following this strategy, every
given convergence criterion (Davis, 1984). Homotopy time t = 1 is reached, a solution of the problem F(x) =
Continuation Methods (HCMs) have emerged as an 0 is found.
alternative strategy, which uses the local convergence Different types of HCMs can be obtained
properties of such methods in a gradual deformation according to the specific form of the auxiliary function
scheme (Allgower and Georg, 1990). They have E(x), the numerical method used to solve H(tk , xH,k ) =
been shown to be accurate, robust, and able to 0 for each discretized value tk , and the selected
converge from almost any arbitrary starting point. homotopic path tracking method (Menglong and
These techniques have been successfully applied to Zhenxin, 2008; Jimenez-Islas et al., 2013). Some
solve polynomial systems of equations, constrained typical auxiliary functions found in the literature
and unconstrained optimization problems, and non- include Fixed Point Homotopy (Eq. (3)), Affine
linear partial differential equations. Therefore, the Homotopy (Eq. (4)) and Newton or global Homotopy
applications of HCMs to solve engineering models (Eq. (5)).
have been numerous, including phase separation and E(x) = x − x0 (3)
transport phenomena problems, among others (Jalali
et al,. 2008; Menglong and Zhenxin, 2008; Laiadi and E(x) = A · (x − x0 ) (4)
Merzougui, 2012; Yang et al., 2012; Jimenez-Islas et E(x) = F(x) − F(x0 ) (5)
al., 2013; Torres-Muñoz et al., 2015).
Where x0 is a known value (initial estimate) and A is a
HCMs are recognized by their potential as a scaling matrix, e.g., Jacobian matrix.
reliable and efficient seeker of all the set of possible Regarding to the numerical method used to solve
roots of a multivariable function from one starting H(tk , xH,k ) = 0 for each discretized value tk , the
point (Rahimian et al. 2011). The general problem basic Predictor-Corrector approach is recommended
that faces HCMs is to find the solution xF of the (Wayburn and Seader, 1987). Taking the total
system of equations F(xF ) = 0, where the mapping derivative with respect to t of the homotopic function
F : <n → <n represents a set of equations, which given by Eq. (2), the Initial Value Problem, IVP (Eq.
in the general case includes transcendental functions. (6)), is obtained (Georg, 1981).
The basic idea of HCMs is fully described in the
literature and deals with the gradual deformation from dxH
= [t · JF (xH ) + (1 − t) · JE (xH )]−1 · [E(xH ) − F(xH )],
an auxiliary function, E(x), to the function to be dt
solved, F(x), where xE , such as E(xE ) = 0, is known. xH (t = 0) = xE
This deformation is developed using an artificial (6)

1014 www.rmiq.org
 Dı́az-Montes et al./ Revista Mexicana de Ingenierı́a Quı́mica Vol. 16, No. 3 (2017) 1013-1029

Where JF (xH ) and JE (xH ) represent the derivative derivative of each component of the homotopic
or Jacobian matrix of the original function and the function H(t, x) = 0 regarding to “p” will be described
auxiliary function, respectively. by Eq. (8).
In a basic Predictor-Corrector scheme, instead of
∂H1 dx1
dH1
+ · · · + ∂H 1 dxn ∂H1 dt ∂H1
∂xn d p + ∂t d p + ∂p
   
directly using the solution obtained at point k as the  dp   ∂x1 d p

initial guess for the problem for the next point k + 1, a ..
 
..
 
 =   = 0
. .

simple strategy such as Euler’s method is applied over    
∂Hn dx1
dHn
+ · · · + ∂H n dxn ∂Hn dt ∂Hn
∂xn d p + ∂t d p + ∂p
   
the IVP given by Eq. (6) using the new initial condition dp ∂x1 d p
xH (tk ) = xH,k , and this numerical result is used as the (8)
initial guess for the problem at the next point k + 1. where x = [x1 · · · xn ]T and H(t, x) = [H1 · · · Hn ]T .
By other side, among different homotopic path Considering that ∂H1 /∂p = · · · = ∂Hn /∂p = 0, and
tracking methods discussed in the literature, the using the compact notation of Eq. (9), Eq. (10) can be
circle tracking method, generalized as hyperspheres, obtained.
exhibited the best behavior (Yamamura, 1993;
dx dt ∂H
Jiménez-Islas, 1996). The hypersphere that is JHx =− (9)
represented by Eq. (7) must be solved for each dp d p ∂t
discretized point k following the algorithm described dx dt ∂H
= − [JHx ]−1 (10)
in Table 1. dp dp ∂t

kx − xH,k k2 + (t − tk )2 = R2 (7)
Fig. 1 shows how, for a two-dimensional xH , the
tangent line from the center of the sphere given by
Euler’s predictor intersects the surface of the sphere,
and this point feeds a locally convergent method that
finds the intersection of the path with the sphere. A
detailed analysis of this tracking method is presented
by Oliveros-Muñoz and Jiménez-Islas, (2013).
One known problem of HCMs appears when the
homotopic path presents turning points as described
by Watson (2001), where an inappropriate step
increase could skip roots of the homotopic path.
Such situation can be avoided by using the
arc length p as alternative tracking parameter
(Klopfenstein, 1961), wherein x and t are taken as Fig. 1. Euler’s predictor for hyperspheres as
functions of this new parameter. Hence, the total homotopic path tracking method.

Table 1. Algorithm steps for hyperspheres as homotopic path tracking method.

Step Description

1. Hypersphere setting Taking the point (tk , xH,k ) as the center of the sphere, with an established radius
R, Eq. (7) is set.
2. Euler’s predictor step Applying Euler’s method over Eq. (6) a tangent line is formed, and the
intersection point of this line and the hypersphere surface is obtained.
The point obtained from step 2 is used as an initial trial to find the simultaneous
3. Corrector step solution of H(t, x) = 0 and Eq. (7), using a locally convergent Newton-like method
as a corrector strategy. The solution will be the next point (tk+1 , xH,k+1 ) of the
homotopic path.
4. Loop The point (tk+1 , xH,k+1 ) is set as the center of a new hypersphere, and the tracking
strategy continues with step 1.

www.rmiq.org 1015
 Dı́az-Montes et al./ Revista Mexicana de Ingenierı́a Quı́mica Vol. 16, No. 3 (2017) 1013-1029

On the other hand, the approximate calculation range of applications for the HCMs methods in
of dt/dp can be performed using the determinant of solving engineering problems. In the next section, the
the homotopic Jacobian matrix (Oliveros-Muñoz and software is presented and then a problem is proposed
Jiménez-Islas, 2013), Eq. (11). as case study. For comparison proposes, the case
study is analyzed by using different approaches, i.e.,
dt
= (−1)n+1 · det(JHx ) (11) classical Newton´s method; Fixed Point Homotopy,
dp Affine Homotopy, and the SphereSolver. Finally, it is
In a Predictor-Corrector scheme, Euler’s method is illustrated the use of SphereSolver for solving different
applied over the problem at the center of the sphere. kinds of typical engineering problems.
Thus, for a given change in the new path tracking
parameter, ∆p, the initial estimate that feeds the
locally convergent method, t0 , and the vector, x0 , 2 Software presentation
are found using eqs. (12) and (13), respectively.
This approach permit to follow any n-dimensional
homotopic path only in terms of t(p). Microsoft Excel© by Microsoft© can be considered
a low-cost program with uses in several engineering
dt
t0 = t + ∆p (12) applications (Huddleston et al., 2004; Wong and Zhou,
dp 2004; Bhattacharjya, 2010; Le Roux et al., 2010).
dx Excel is widely available and easily adaptable to
x0 = x + ∆p (13) individual needs for solving engineering problems.
dp
Although several add-ins and templates have been
When this approach is combined with the hypersphere developed in Excel, as far as we know, there have been
radius R (Eq. (7)), the change in the arc length can no tools capable of finding several roots of system
be described by Eq. (14), where the sign for the square of non-linear equations using a globally convergent
root is set according to the positive or negative advance method. The newly developed add-in SphereSolver
of p (Oliveros-Muñoz and Jiménez-Islas, 2013). was designed with a friendly user interface, and with
v
u default settings suitable for users who are not familiar
R2
t
∆p = ± (14) with the numerical method. The size of the add-in is
 2
1 + k ddxp k2 + ddtp 143 kB.
The code was developed in Visual Basic for
According to the literature survey, it was found that Applications 7.1©. A description of the main
different computational tools have been developed algorithm is included as supporting material.
based on HCMs, most of them implemented The developed add-in uses Newton Homotopy;
in scientific computation packages as Matlab hyperspheres as the tracking device of the homotopic
(Thompson, 2006; Lee et al., 2008); while the path; and second-order finite central differences to
use of HCMs as an add-in for the widely known approximate the required derivatives.
software Microsoft Excel Office is scarce and limited To use the developed add-in, Microsoft Office
to Fixed Point and Affine Homotopies (Henao and Excel spreadsheets are used to set up the model that
Velásquez, 2004), without including hyperspheres as contains the system of equations to be solved. Specific
homotopic path tracking method, neither arc length p information shown in Figs. 2 and 3 is related with the
as alternative tracking parameter. case study described in the next section. As shown in
In this contribution, a new code for HCMs Fig. 2, at least two set of cells are required to capture
developed in Visual Basic for Applications 7.1©, the information of the model from the spreadsheet.
to be used as Microsoft Office Excel add-in The first set contains n cells, where the different values
(named SphereSolver) is presented. SphereSolver that the unknown variables adopt during the execution
uses Newton Homotopy (Eq. (5)) as auxiliary of the algorithm are downloaded (cell D7, Fig. 2). The
equation, hyperspheres as the tracking device of initial values for the HCM are recommended to be
the homotopic path; arc length p as alternative entered there. The second set contains other n cells,
tracking parameter; and second-order finite central where the formulas of the system of equations are
differences to approximate the required derivatives. introduced. They must depend directly or indirectly
The simplicity associated to the use of widely known on the first set of cells (cell F7, Fig. 2). In general,
software as Microsoft Office Excel improves the SphereSolver changes the values of the cells that

1016 www.rmiq.org
 Dı́az-Montes et al./ Revista Mexicana de Ingenierı́a Quı́mica Vol. 16, No. 3 (2017) 1013-1029

represent the unknown variables as the homotopic path satisfies the equations with in a convergence criterion.
is followed, and the solution will be the values that

Fig. 2. Example of a set up in a spreadsheet used to capture the information of a model (case study presented in the
next section).

Fig. 3. SphereSolver’s interface. a) Page to capture the information of the system of equations from an Excel Office
spreadsheet. b) Page to acquire basic numerical parameter values. c) Page to set advanced parameters for the HCM.

www.rmiq.org 1017
 Dı́az-Montes et al./ Revista Mexicana de Ingenierı́a Quı́mica Vol. 16, No. 3 (2017) 1013-1029

After installation, the add-in can be executed from where parameters a and b are functions of the critical
the Excel Office ribbon with the name SphereSolver. properties of the gas, and can be described by eqs. (16)
The interface consists of three pages: Problem (Fig. and (17), respectively (Murdock, 1993).
3a), Settings (Fig. 3b) and Advanced (Fig. 3c). In
Problem (Fig. 3a), using two different reference 0.42748023354R2 T c2.5
a= (16)
selection controls, the cells containing the functions Pc
and the variables are captured. In Settings (Fig. 3b),
0.086640349965RT c
some basic settings for the HCM are entered, i.e., b= (17)
tolerance of convergence (ε), the initial value of the Pc
radius (Rinit ) and the maximum number or iterations The case study consists in applying the Redlich-
in the correction step (Imax ). In this page the user Kwong model (Eq. (15)) to estimate the molar volume
has the option of taking one particular initial value v of n-butane with the conditions and parameters given
for all the variables, or taking the values initially in Table 2 (Smith et al. 2005).
entered in the spreadsheet. Finally in Advanced (Fig.
3c), advanced numerical parameters can be entered, 3.1 Use of the classical Newton’s method in
i.e., the minimum slope value of both sides of a return
and the maximum error at an intersection with the
the case study
hyperplane t = 1. The option to set the domain for Eq. (15) is written as a function of molar volume, Eq.
the homotopic parameter t is also available for users (18)
which are familiar with the HCMs and the homotopic RT a
F(v) = − √ −P=0 (18)
function behavior. In this study, the machine used v−b T v(v + b)
has an Intel(R) Core(TM) i7-4510U CPU @ 2.00GHz
A plot F(v) (Eq. 18) using condition from Table 1
60 GHz processor and 8.00 GB RAM installed.
reveals several asymptotes (Fig. 4).
Because the developed add-in is intended for a wide
To apply classical Newton’s method (Ramı́rez et
variety of users, the characteristics of the machine are
al. 2006; Bürgisser and Cucker, 2013) and find the root
intentionally easy to obtain.
of the function given by Eq. (18), the derivative with
respect to the unknown variable v is required, i.e., Eq.
(19).
3 Case of study RT a(2v + b)
F 0 (v) = − 2
+ √ (19)
(v − b) T [v(v + b)]2

Redlich-Kwong (Eq. (15)) represents one of the most Iterations with the classical Newton’s method are
classical and widely used equations of state in the given by vk+1 = vk − F(vk )/F 0 (vk ), in this case
calculation of thermodynamic properties required in represented by Eq. (20).
several chemical engineering applications (Murdock,  
1993). RT √ a
vk −b − − P
T vk (vk +b)
RT a vk+1 = vk −   (20)
P= a(2vk +b)
−
v − b T 1/2 v(v + b)
(15) −RT
2 + √
v −b)
k T [vk (vk +b)]
2

Table 2. Values for Redlich-Kwong cubic equation.

Parameter Name Units Value
P Pressure atm 50
T Temperature K 298.15
Pc Critical pressure atm 37.5*
Tc Critical temperature K 425.2*
R Universal gas constant L·atm/mol K 0.08205
* Taken from (Smith et al., 2005)

1018 www.rmiq.org
 Dı́az-Montes et al./ Revista Mexicana de Ingenierı́a Quı́mica Vol. 16, No. 3 (2017) 1013-1029

Fig. 4. Graph of the function F(v) = RT

v−b − √ a − P using conditions from Table 1.
T v(v+b)

Fig. 5. Estimative of vk+1 as a function of the estimate vk of Newton’s method (Eq. (18) with conditions from Table
1); a) Higher scale; b) Smaller scale.

www.rmiq.org 1019
 Dı́az-Montes et al./ Revista Mexicana de Ingenierı́a Quı́mica Vol. 16, No. 3 (2017) 1013-1029

Typical drawbacks of Newton’s method are related Thus, when x − x0 = F(x), t becomes ±∞. Hence,
with the adequate choice of the initial guess. Fig. 5 the horizontal asymptote in Fig. 6 is due to the type
shows the behavior of the first estimate as a function of of Homotopy itself, and not due to the nature of the
the initial value for Newton’s method in two different original function.
scales. The sequence reveals only one continuous For the given set of parameters (Table 1), the
interval containing the root where convergence horizontal asymptote occurs at v =0.109557347
can be achieved, i.e., from 0.0806047752495982 L/mol, which is 0.0195% away from the root. In other
to 0.127396259816616. In other words, the initial words, for Fixed Point Homotopy, following the path
estimate must go from -26.41% to 16.31% of in the real domain from the initial estimate will not
the solution. Thus, the interval of convergence of lead to the answer at t = 1.
Newton’s method for this case of study is tight.
Importantly, the common initial value for the molar
volume in a real model is the one obtained from 3.3 Use of the Affine Homotopy in the case
the ideal gas equation, and in this case it is
study
v=RT/P=0.48926415 L/mol, which falls outside the
interval of convergence. The homotopic path for the implicit equation
H(xH , t) = 0 using Affine Homotopy (Eqs. (2) and (4))
3.2 Use of the Fixed Point Homotopy in the over the case of study is shown in Fig. 7. It is observed
case study that the homotopic path presents similar behavior
to the Fixed Point Homotopy (Fig. 6). In a similar
Using the same case study, when the molar volume fashioned way, the horizontal asymptote occurs when
obtained from the ideal gas equation is used as the the denominator of the Eq. (22) equals 0, i.e., when
initial estimate (i.e., x0 = v0 = RT/P =0.48926415, F 0 (x)(x − x0 ) = F(x).
t = 0), and Fixed Point Homotopy (eqs. (2) and (3))
is applied; the homotopic path shown in Fig. 6 is F 0 (x0 )(x − x0 )
obtained. t= (22)
F 0 (x
0 )(x − x0 ) − F(x)
In Fig. 6, the portion of the homotopic path that
is above the horizontal asymptote comes from the For the case study this scenario is at v=
point (xH = xE = v0 =0.48926415,t = 0). Due to 0.111378181943983 L/mol, 1.68% away from the
the horizontal asymptote, there is no continuous path root. Although not shown in Fig. 7, if the path from
that connects the initial point to the desired root at the initial estimate is followed upwards, it leads to
t = 1. Setting the Homotopy function equal to 0, and a vertical asymptote at t = 1, and the other path
combining eqs. (2) and (3), Eq. (21) can be obtained. has a horizontal asymptote at the mentioned point.
(x − x0 ) This same behavior was observed for the Fixed Point
t= (21) Homotopy.
(x − x0 ) − F(x)

Fig. 6. Homotopic path for the study case of the Fig. 7. Homotopic path for the study case of the
Redlich-Kwong cubic equation under the conditions Redlich-Kwong cubic equation under the conditions
given in Table 1 using Fixed Point Homotopy. given in Table 1 using Affine Homotopy.

1020 www.rmiq.org
 Dı́az-Montes et al./ Revista Mexicana de Ingenierı́a Quı́mica Vol. 16, No. 3 (2017) 1013-1029

root in a continuous path. Thus, the new practical

continuous interval of convergence that contains the
root is from b =0.0806047752495982 to ∞. Further
results are shown in the use of SphereSolver, which
include Newton Homotopy with in the code.

3.5 Use of the SphereSolver in the case

study
The SphereSolver presented in Section 2 was also used
with the case study, using several initial values for the
HCM. Three of them are shown in Table 3, with their
corresponding performance within the application.
Fig. 8. Homotopic path for the study case of the The initial value of the first case in Table 3 was
Redlich-Kwong cubic equation under the conditions obtained using the ideal gas equation of state, which
given in Table 1 using Newton Homotopy. is common for this type of problem and, as previously
shown, is not useful in this case when Fixed Point
According to the previous analysis, for the Homotopy or Affine Homotopy are used. Even when
Fixed Point Homotopy or Affine Homotopy, there more initial values were used (with SphereSolver),
is no continuous path from the initial value they presented similar behavior.
v0 =0.48926415 to the solution of the problem at As analyzed in Section 3.1, the highest continuous
v =0.109535967857868. Again, as observed for Fixed interval where the classic Newton’s method converges
Point Homotopy, the horizontal asymptotes analyzed is from 0.0806047752495982 to 0.127396259816616,
do not rise from the nature of the function but from the i.e., the initial estimate must go from -26.41% to
type of homotopy. 16.31% of the solution. Using SphereSolver (which
applies Newton Homotopy), any value higher than
0.0806047752495982 can be used where a horizontal
3.4 Use of the Newton Homotopy in the asymptote exists due to the nature of the original
case study function and not the type of HCM, will go from the
Using the Newton Homotopy over the case study, initial value to the root of the equation at t = 1, i.e., the
i.e., combination of eqs. (2) and (5), a path from the interval of convergence is extended massively.
initial estimate, v0 =0.48926415, to the solution of the
problem exists, as shown in Figure 8.
This case study has shown the advantage of 4 Additional engineering
using the Newton Homotopy over the Fixed Point problems solved with
and Affine Homotopy. According to the results, the
convergence interval will be the same continuous SphereSolver
domain of the function F(v), given by Eq. (18),
where the root is contained. For other initial values,
a similar homotopic path could be found, but due Three additional engineering models that include non-
to the discontinuities created by the nature of the linear system of equations were selected to be solved
original function, it will not be possible to obtain the by using SphereSolver.

Table 3. Performance of SphereSolver over the case study with three different initial values. Initial sphere radius:
0.1. Tolerance: 1 × 10−6 . Max. Iterations for each step: 2. Max. error for t = 1: 5 × 10−7 .
Case Initial value Time [s] Number of Spheres Norm
1 0.48926 39.93 393 3.42 × 10−7
2 10 33.45 337 1.87 × 10−10
3 20 35.46 360 4.91 × 10−7

www.rmiq.org 1021
 Dı́az-Montes et al./ Revista Mexicana de Ingenierı́a Quı́mica Vol. 16, No. 3 (2017) 1013-1029

Table 4. Value of parameters for the adiabatic equilibrium temperature model, taken from Fogler (2016).
Parameter Name Units Value
HA0 Heat of formation for A at 298 K cal/mol -40,000
HB0 Heat of formation for B at 298 K cal/mol -60,000
C PA Heat capacity of A at 298 K cal/mol K 50
C PB Heat capacity of B at 298 K cal/mol K 50
R Universal gas constant cal/mol K 1.987
Ke Equilibrium constant at 298 K — 100,000

Table 5. Tests of SphereSolver solving Eq. (23) by using different initial guess. Minimum slope in a return: 0.005.
Max. error for t = 1 : 5 × 10−7 .
Case Initial value Initial sphere radius Tolerance Max. Iterations Time (s) Spheres Norm
for each step
1 100 10 1 × 10−3 2 140.84 517 7 × 10−8
2 10000 1000 1 × 10−6 5 29.55 113 4 × 10−7
3 100000 10000 1 × 10−3 5 26.18 125 8 × 10−4
4 30 100 1 × 10−3 5 22.09 93 4 × 10−7
5 30 10 1 × 10−3 5 22.6 94 1 × 10−7

Selected models include the adiabatic equilibrium A typical approach for solving the problem is by using
temperature in a chemical reactor, the Benzene Newton’s method, in such case, the derivative of Eq.
production in a plant performance evaluation, and an (23), (Eq. (24)), should be considered.
electronic circuit with two-tunnel exponential diodes.  0 0  
(HB −HA ) 1 1
Depending on the nature of the functions, different dF K e exp R T1 − T (HB0 − HA0 )
aspects are highlighted for each simulation, i.e., severe =
dt 
(H 0 −H 0 )   2
 RT 2
non-linearity, restricted convergence with traditional 1 + Ke exp B R A T11 − T1
methods, several sets of possible solutions, and so on.
C PA
− (24)
−(HB0 − HA0 )
4.1 Chemical reactor: Adiabatic equilibrium
Fig. 9 shows the behavior of the function
temperature convergence. For the Newton’s method, it is achieved
if the initial trial is within the interval 400.569432
The highest conversion that can be achieved in
< T < 522.389085, which is -13% to +13.46% away
reversible reactions is obtained at the equilibrium
from the root at T = 460.399642820735. In Fig. 9,
temperature. Fogler (2016) presents an example for
the dotted vertical lines are part of the convergence
an elementary solid-catalyzed liquid-phase reaction
interval for Newton’s method
A ↔ B, where the adiabatic equilibrium temperature
Alternatively, SphereSolver add-in was used to
T is determined when pure A is fed to the reactor at a
solve the non-linear equation under different initial
temperature of 3000 K. The resulting Eq. (23), which
guess. The convergence criteria were satisfied in each
deduction is presented by Fogler (2016), should be
occasion, and the performance of the method for five
solved to find the adiabatic equilibrium temperature.
different cases is reported in Table 5.
Supplied information is given in Table 4.

(HB0 −HA0 )  1  4.2 Plant performance evaluation: Benzene
Ke exp R T1 − T1
F(T ) =
production
(HB0 −HA0 )  1
 
1 + Ke exp R T1 − T1 The hydrodealkylation of toluene to benzene is an
C PA (T − T 0 ) important petrochemical process; an industrial process
− =0 (23) for this transformation is shown in Fig. 10 (Henao,
−(HB0 − HA0 ) 2010).

1022 www.rmiq.org
 Dı́az-Montes et al./ Revista Mexicana de Ingenierı́a Quı́mica Vol. 16, No. 3 (2017) 1013-1029

path in terms of the homotopic parameter t and the

unknown variables is not practical. As shown in the
Introduction section, the arc length p of the Homotopy
function could be used as a new parameter to represent
the path. In this manner, a graph of p vs. t describes
the hotomopic path accurately, and each intersection
with the horizontal line t = 1 represents a root of the
original function. A wide set of initial values were
tested, and at least one of the roots listed in Table 8 was
always found. In Fig. 11, different homotopic paths are
shown, where the same initial value, Xo , were used
for all unknown variables. Most of them, except for
Xo = 6, led to the three possible roots. For Xo = 6, only
Fig. 9. Graph of the function given by Eq. (23).
one root was found showing that the initial value of the
Values wherein the dotted vertical lines are part
homotopic path is still a relevant subject and must be
of the convergence interval for Newton’s method.
explored in future researches.

4.3 Electronic circuits: Circuit with two-

tunnel exponential diodes
Torres-Munoz et al. (2014), presents a circuit with
two tunnel diodes, one voltage source and a resistor
in series. Using the Kirchoff’s current law, with E=1,
R=20Ω, I p = 100×10−3 , V p = 50×10−3 , I0 = 1×10−9
and q/RT = 40, eqs. (25) - (27) are obtained.
1 v1
f1 = − + I v1 = 0 (25)
20 20
1 v1
f2 = − + + 2(v1 − v2 )e1−20v1 +20v2
20 20
Fig. 10. Hydrodealkylation toluene plant for producing
+ 1 × 10−9 e40v1 −40v2 = 0 (26)
benzene.
1−20v1 +20v2
f3 = 2(v1 − v2 )e + 1 × 10 e
−9 40v1 −40v2

The global reaction uses toluene and hydrogen as − 2v2 e1−20v2 −9 40v2
− 1 × 10 e =0 (27)
reactants to obtain bencene and methane as products
by the chemical reaction: C6 H5 CH3 +H2 →C6 Torres-Muñoz et al., (2014) reported 5 solutions to
H6 +CH4 . The system of equations that represents the this system of equations. Using SphereSolver with 0.5
process at T=620°C and P=3400 KPa is depicted in as the initial value for all three variables, 9 solutions
Table 6. Table 7 presents the meanings of the variables were found (Table 10). Therefore, it is demonstrated
used in the model. again that SphereSolver is a powerful tool for solving
Resulting non-linear system includes 44 equations different kind of engineering problems.
and 58 unknowns. To fix the degrees of freedom,
the information in Table 8 is used. Upon setting
these values, the system of equations becomes a Conclusions
performance evaluation. Henao (2010) reported one
solution of this problem. Using SphereSolver add-in,
three solutions were found (Table 9). The answer that SphereSolver, an add-in for use in Microsoft Office
satisfies the equilibrium conditions will be the one that Excel© to solve systems of non-linear equations,
represents the phenomena most accurately, but this was developed and applied to several engineering
question is beyond the scope of the current research. models. The Newton Homotopy Method was chosen
Taking into account that this system of equations over Affine and Fixed Point Homotopy as the main
is multivariable, the representation of the homotopic algorithm.

www.rmiq.org 1023
 Dı́az-Montes et al./ Revista Mexicana de Ingenierı́a Quı́mica Vol. 16, No. 3 (2017) 1013-1029

Table 6. Equations for benzene production model. Index 1, 2, 3, 4 refers to toluene; hydrogen; benzene; and
methane, respectively.

Unit Type of equation Equations

Mixer M-1001 Material balances C1 = A1 + B1 + H1 + L1
C2 = A2 + B2 + H2 + L2
C3 = A3 + B3 + H3 + L3
C4 = A4 + B4 + H4 + L4
Furnace F-1001 Material balances D1 = C 1
D2 = C 2
D3 = C 3
D4 = C 4
Reactor R-100 Material balances E 1 = D1 − r · V
(Continuous E 2 = D2 − r · V
stirred- E 3 = D3 + r · V
tank E 4 = D4 + r · V
reactor Reaction kinetics r = kc1,E c1/22,E
model) equations k = k0 exp(−E
 a /RT E ) 
PE E1
c1,E = RT E  E 1 +E 2 +E 3 +E 4 
PE E2
c2,E = RT E E 1 +E 2 +E 3 +E 4

Cooler C-1001 Material balances F1 = E1

F2 = E2
F3 = E3
F4 = E4
Vessel V-1001 Material balances G1 + J1 = F1
G2 + J2 = F2
G3 + J3 = F3
G4 + J4 = F4
Operation equations G1 = 0.0001 · F1
G2 = 0.99 · F2
G3 = 0.001 · F3
G4 = 0.95 · F4
Splitter S-1001 Material balances H1 + I1 = G1
H2 + I2 = G2
H3 + I3 = G3
H4 + I4 = G4
Operation equations H1 = α · G 1
H2 = α · G 2
H3 = α · G 3
H4 = α · G 4
Distillation Tower T-1001 Material balances K1 + L1 = J1
K2 + L2 = J2
K3 + L3 = J3
K4 + L4 = J4
Operation equations K1 = 0.0001J1
K2 = 0.9999 · J2
K2 = 0.999 · J3
K2 = 0.9999 · J4

1024 www.rmiq.org
 Dı́az-Montes et al./ Revista Mexicana de Ingenierı́a Quı́mica Vol. 16, No. 3 (2017) 1013-1029

Fig. 11. Discretized homotopic path t vs. p for differential initial conditions, Xo , and different scales.

Table 7. Notation for the benzene production model.

Name Definition
Xi Molar flow of substance i in stream X, for i = 1, ..4
V Reactor volume
ci,X Molar concentration of substance i in stream X
TE Stream E temperature
PE Stream E pressure
r Reaction rate
k Reaction rate coefficient
α Molar flow ratio of stream H to stream G
R Universal gas constant: 8.314 kJ/kmol K

Table 8. Parameter values for the benzene production model.

Parameter Units Value
k0 kmol1/2 /m3/2 h 2.268 × 1014
Ea kJ/K kmol 2.18 × 105
V m3 100
α - 0.7
TE K 893.15
PE kPa 3 400
A1 , A3 , A4 kmol/h 0
A2 kmol/h 130
B2 , B3 , B4 kmol/h 0
B1 kmol/h 100

www.rmiq.org 1025
 Dı́az-Montes et al./ Revista Mexicana de Ingenierı́a Quı́mica Vol. 16, No. 3 (2017) 1013-1029

Table 9. List of solutions for the benzene production model.

Variable Units Solution 1 Solution 2 Solution 3
C1 kmol/h 225.4 15 524.9 532 750.5
C2 kmol/h 197.76 202.25 354.04
C3 kmol/h 0.17016 0.16678 0.05234
C4 kmol/h 198.48 194.53 61.05
D1 kmol/h 225.4 15 524.9 532 750.5
D2 kmol/h 197.76 202.25 354.04
D3 kmol/h 0.17016 0.16678 0.05234
D4 kmol/h 198.48 194.53 61.05
E1 kmol/h 125.4 15 427.0 532 719.8
E2 kmol/h 97.77 104.25 323.29
E3 kmol/h 100.15 98.16 30.80
E4 kmol/h 298.46 292.53 91.80
F1 kmol/h 125.4 15 427.0 532 719.8
F2 kmol/h 97.77 104.25 323.29
F3 kmol/h 100.15 98.16 30.80
F4 kmol/h 298.46 292.53 91.80
G1 kmol/h 0.013 1.543 53.272
G2 kmol/h 96.80 103.21 320.05
G3 kmol/h 0.1002 0.0982 0.0308
G4 kmol/h 283.54 277.90 87.21
H1 kmol/h 0.0088 1.0799 37.2904
H2 kmol/h 67.76 72.25 224.04
H3 kmol/h 0.0701 0.0687 0.0216
H4 kmol/h 198.478 194.530 61.045
I1 kmol/h 0.0038 0.4628 15.9816
I2 kmol/h 29.039 30.963 96.016
I3 kmol/h 0.0300 0.0294 0.0092
I4 kmol/h 85.062 83.370 26.162
J1 kmol/h 125.4 15 425.4 532 666.5
J2 kmol/h 0.978 1.043 3.233
J3 kmol/h 100.05 98.06 30.77
J4 kmol/h 14.923 14.626 4.590
K1 kmol/h 0.013 1.543 53.267
K2 kmol/h 0.978 1.042 3.233
K3 kmol/h 99.954 97.965 30.743
K4 kmol/h 14.922 14.625 4.589
L1 kmol/h 125.4 15 423.9 532 613.2
L2 kmol/h 0.000098 0.000104 0.000323
L3 kmol/h 0.100054 0.098063 0.030773
L4 kmol/h 0.001492 0.001463 0.000459
c1,E kmol/m3 0.092370 0.443639 0.457490
c2,E kmol/m3 0.071992 0.002998 0.000278
r kmol/m3 h 0.999837 0.979947 0.307518
k kmol1/2 /m3/2 h 40.342 40.342 40.342

The proposed algorithm uses hyperspheres as a algorithm are Euler’s tangent plane as a predictor and
tracking method for the homotopic path, which has Newton’s method as the locally convergent method
been proven to be a reliable method, following the path for the corrector step at each discretized point of the
without losing it. The other two main elements of the homotopic path.

1026 www.rmiq.org
 Dı́az-Montes et al./ Revista Mexicana de Ingenierı́a Quı́mica Vol. 16, No. 3 (2017) 1013-1029

Table 10. List of solutions for the circuit with two-tunnel exponential diodes.
Solution # Iv1 v1 v2 Norm
1 -0.0282803 0.43439597 0.4285488 4.70 × 10−7
2 -0.0421756 0.15648531 0.14713161 2.99 × 10−7
3 -0.0488741 0.02252113 0.01126054 1.98 × 10−7
4 -0.0421755 0.15648584 0.00935372 1.99 × 10−7
5 -0.0336138 0.32772635 0.16386296 2.18 × 10−7
6 -0.0189197 0.62160289 0.41818366 8.27 × 10−7
7 -0.0099145 0.80171369 0.40085755 5.44 × 10−7
8 -0.0189198 0.62159913 0.20341635 6.88 × 10−7
9 -0.0282804 0.4343955 0.00584709 4.31 × 10−7

Based on the obtained results, SphereSolver is Pp. 283-294. Springer Science & Business
shown to be a powerful tool for developing simulations Media.
in Microsoft Excel office, finding the solutions of
Cveticanin, L., Kalami-Yazdi, M. and Askari,
the system of non-linear equations that represents the
H. (2012). Analytical approximations to the
corresponding model. By comparison on the use of
solutions for a generalized oscillator with
SphereSolver for solving the case study (section 3)
strong nonlinear terms. Journal of Engineering
and additional engineering problems (section 4) it is
Mathematics 77, 211-223.
concluded that SphereSolver improves the numerical
solution by increasing the interval range of the Davis, M. E. (1984). Numerical Methods and
possible initial values that lead to the solution of the Modeling for Chemical Engineers, John Wiley
problems and/or increasing the amount of found roots. & Sons, New York.
SphereSolver uses Microsoft Excel Office ©, thus
Fogler, H. S. (2016). Elements of Chemical Reaction
it offers a user-friendly spreadsheet environment,
Engineering, fifth ed., Pearson Education, New
where models represented by algebraic equations
Jersey.
(including transcendental functions) are introduced
easily, presenting the information and results in an Georg, K. (1981). On tracing an implicitly defined
organized and personalized set up. curve by quasi-Newton steps and calculating
bifurcation by local perturbations. SIAM
Acknowledgements Journal on Scientific & Statistical Computing
2, 35-50.
The authors are grateful for the financial support
provided by Tecnológico Nacional de México Henao, C. A. and Velásquez, J. A. (2004). Simulación
(TecNM) Grant 5037.13-P de Procesos Quı́micos Empleando MS Excel®
(I) -Aplicación de Métodos Homotópicos.
Ingenierı́a Quı́mica 26, 16-22.
References Henao, C. A. (2010). Simulación and Evaluación
de Procesos Quı́micos: Herramientas Básicas
para la Sı́ntesis de Procesos, 1a ed., Universidad
Allgower, E. L. and Georg, K. (1990). Numerical Pontificia Bolivariana, Medellı́n.
Continuation Methods: An Introduction,
Springer, Berlin. Huddleston, D., Alarcon, V. and Chen, W. (2004).
Water distribution network analysis using Excel.
Bhattacharjya, R. K. (2010). Solving groundwater Journal of Hydraulic Engineering 130, 1033-
flow inverse problem using spreadsheet solver. 1035.
Journal of Hydraulic Engineering 16, 472-477.
Jalali, F., Seader, J. D. and Khaleghi, S. (2008).
Bürgisser, P. and Cucker, F. (2013). Homotopy Global solution approaches in equilibrium and
Continuation and Newton’s Method. En: stability analysis using homotopy continuation
Condition: The Geometry of Numerical in the complex domain. Computers & Chemical
Algorithms, (P. Bürgisser, and F. Cucker, eds.), Engineering 32, 2333-2345.

www.rmiq.org 1027
 Dı́az-Montes et al./ Revista Mexicana de Ingenierı́a Quı́mica Vol. 16, No. 3 (2017) 1013-1029

Jiménez-Islas, H. (1996). SEHPE: Programa para la correction step in homotopic continuation

solución de sistemas de ecuaciones no lineales methods. Chemical Engineering Science 97,
mediante método homotópico con seguimiento 413-429.
hiperesférico. Avances en Ingenierı́a Quı́mica 6,
174-179. Rahimian, S. K., Jalali, F., Seader, J. D. and
White, R. E. (2011). A robust homotopy
Jimenez-Islas, H., Martinez-Gonzalez, G. M., continuation method for seeking all real roots
Navarrete-Bolanos, J. L., Botello-Alvarez, J. E., of unconstrained systems of nonlinear algebraic
and Oliveros-Muñoz, J. M. (2013). Nonlinear and transcendental equations. Industrial &
homotopic continuation methods: a chemical Engineering Chemistry Research 50, 8892-
engineering perspective review. Industrial & 8900.
Engineering Chemistry Research 52, 14729-
Ramı́rez, A. D., Rodriguez, D. and González, N. D.
14742.
(2006). Estudio experimental and numérico de
Jiménez-Islas, H., Calderon-Ramı́rez, M., Molina- la morfologı́a de mezclas PS/PEAD preparadas
Herrera, F., Martı́nez-Gonzalez, G., Navarrete- por extrusión. Revista Mexicana de Ingenierı́a
Bolanos, J. and Castrejon-González, E. (2014). Quı́mica 5, 67-78.
Algoritmo de bajo consumo de RAM para Rodrı́guez, R. and Niño, Z. (2016). Evaluación
resolver problemas de convección natural 3-d, de herramientas computacionales gratuitas para
usando colocación ortogonal. Revista Mexicana la simulación de procesos de combustión en
de Ingenierı́a Quı́mica 13, 259-277. motores de encendido por chispa. Revista
Klopfenstein, R. (1961). Zeros of nonlinear Mexicana de Ingenierı́a Quı́mica 15, 977-984.
functions. Journal of the ACM 8, 366-373. Smith, J. M., van Ness, H. G., and Abbott, M. M.
(2005). Introduction to Chemical Engineering
Laiadi, D. and Merzougui, A. (2012). Homotopy
Thermodynamics, 7th Ed., McGraw-Hill Book
method to predict liquid-liquid equilibria for
Co., New York.
ternary mixtures of (water + carboxylic acid +
organic solvent). Fluid Phase Equilibria 313, Thompson, S. L. (2006). A friendly Fortran DDE
114-120. solver. Applied Numerical Mathematics 56,
503-516.
Le Roux, J. P., Demirbilek, Z. and Flemming,
B. W. (2010). WAVECALC: an Excel-VBA Torres-Munoz, D., Vazquez-Leal, H., Hernandez-
spreadsheet to model the characteristics of fully Martinez, L. and Sarmiento-Reyes, A. (2014).
developed waves and their influence on bottom Improved spherical continuation algorithm with
sediments in different water depths. Geo-Marine application to the double-bounded homotopy
Letters 30, 549-560. (dbh). Computational and Applied Mathematics
33, 147-161.
Lee, T. L., Li, T. Y. and Tsai, C. H. (2008).
HOM4PS-2.0: a software package for solving Torres-Muñoz, D., Hernandez-Martinez, L.
polynomial systems by the polyhedral and Vázquez-Leal, H. (2015). Spherical
homotopy continuation method. Computing 83, Continuation Algorithm with Spheres of
109-133. Variable Radius to Trace Homotopy Curves.
International Journal of Computational and
Menglong, S. and Zhenxin, L. (2008). Modified Applied Mathematics 2, 1-13.
homotopy methods to solve fixed points of self-
mapping in a broader class of nonconvex sets. Verschelde, J. and Cools, R. (1994). Symmetric
Applied Numerical Mathematics 58, 236-248. homotopy construction. Journal of
Computational and Applied Mathematics 50,
Murdock, J. W. (1993). Fundamental Fluid 575-592.
Mechanics for the Practicing Engineer, first
ed., CRC Press, New York. Ward J. P. and King, J. R. (2012). Thin-film
modelling of biofilm growth and quorum
Oliveros-Muñoz, J. and Jiménez-Islas, H. (2013). sensing. Journal of Engineering Mathematics
Hyperspherical path tracking methodology as 73, 71-92.

1028 www.rmiq.org
 Dı́az-Montes et al./ Revista Mexicana de Ingenierı́a Quı́mica Vol. 16, No. 3 (2017) 1013-1029

Watson, L. (2001). Globally convergent homotopy Resources Management, (G. Sehlke, D. F.

methods. En: Encyclopedia of Optimization, (C. Hayes and D. K. Stevens, eds.), Pp. 1-8.
A. Floudas, P. M. Pardalos, eds.), Pp. 893-898, American Society of Civil Engineers, Salt Lake
Springer, US. City, UT, USA.
Wayburn, T. L. and Seader, J. D. (1987). Homotopy Yamamura, K. (1993). Simple algorithms for tracing
continuation methods for computer-aided solution curves. IEEE Transactions on Circuits
process design. Computers & Chemical and Systems I 40, 537-541.
Engineering 11, 7-25.
Yang, F., Chen, K. and Yu, B. (2012). Homotopy
Wong, T. S. and Zhou, M. (2004). Determination method for a mean curvature-based denoising
of critical and normal depths using excel. En: model. Applied Numerical Mathematics 62,
Critical Transitions in Water and Environmental 185-200.

www.rmiq.org 1029

View publication stats