SciELO - Scientific Electronic Library Online

 
vol.32Variation of zooplankton in two urban lakes located in recreational parks in the state of Morelos, MexicoDeveloping productive projects in rural and indigenous areas of the state of Oaxaca: Student experiences author indexsubject indexsearch form
Home Pagealphabetic serial listing  

Services on Demand

Journal

Article

Indicators

Related links

  • Have no similar articlesSimilars in SciELO

Share


Acta universitaria

On-line version ISSN 2007-9621Print version ISSN 0188-6266

Acta univ vol.32  México  2022  Epub Aug 28, 2023

https://doi.org/10.15174/au.2022.3358 

Artículos

Enseñanza del método de continuación homotópica con seguimiento hiperesférico para estudiantes de ingeniería

Teaching the homotopic continuation method with hyperspherical tracking for engineering students

Miriam Lucero Quemada-Villagómez1 

María de la Luz López-González1 

Juan Manuel Oliveros-Muñoz2 

Hugo Jiménez-Islas3  * 

1Departamento de Ingeniería Bioquímica, Tecnológico Nacional de México.

2Departamento de Ingenierías, Arquitectura y Diseño. Universidad Iberoamericana Torreón. Coahuila, México.

3Departamento de Ingeniería Bioquímica, Tecnológico Nacional de México, Celaya, Guanajuato, México. Tel. 461 61 17575 ext. 5469.


Resumen

La continuación homotópica es un método numérico de convergencia global que permite encontrar las soluciones a sistemas de ecuaciones algebraicas no lineales, sin embargo, no es muy utilizado como un método estandarizado debido a su aparente complejidad. En este trabajo se da a conocer, de una manera educativa, el método de continuación homotópica con seguimiento hiperesférico aplicado en la resolución de sistemas de ecuaciones algebraicos no lineales reportados en la literatura, los cuales están presentes en diversas áreas de la ingeniería.

Palabras clave: Ecuaciones algebraicas; homotopía; continuación homotópica; seguimiento hiperesférico

Abstract

Homotopic continuation is a global convergence numerical method that allows finding solutions to systems of nonlinear algebraic equations; however, it is not widely used as a standardized method due to its apparent complexity. This work presents, in an educational way, the homotopic continuation method with hyperspherical tracking applied in the resolution of nonlinear algebraic equations systems reported in the literature, which are present in various areas of engineering.

Keywords: Algebraic equations; homotopy; homotopic continuation; hyperspheric tracking

Introducción

La determinación de los vectores solución de un sistema de ecuaciones es uno de los problemas más frecuentemente encontrados en todas las áreas de la ingeniería, ya que se presentan como la solución de una gran variedad de problemas físicos. El método de Newton-Raphson, tradicionalmente, se ha utilizado para resolver sistemas de ecuaciones no lineales. Este algoritmo requiere de una buena aproximación inicial para que exista convergencia y, aunque en muchos casos se obtiene una solución, con frecuencia aparecen situaciones donde el método de Newton falla, o bien, se presentan multiplicidad de soluciones (Jiménez-Islas et al., 2013). La continuación homotópica es un método utilizado para encontrar las raíces de sistemas de ecuaciones algebraicas no lineales. El método consiste, principalmente, en perturbar el sistema de ecuaciones con la introducción de una ecuación homotópica y su respectivo parámetro homotópico. El seguimiento homotópico, sin embargo, es un método clásico y poco explorado. En 1988 los primeros investigadores utilizaron el método de continuación homotópica para encontrar las raíces de diversos sistemas de ecuaciones algebraicas no lineales con múltiples soluciones (Kuno & Seader, 1988).

Posteriormente, el seguimiento homotópico fue desarrollado de manera limitada debido a su naturaleza iterativa y a las capacidades computacionales de la época. Fue hasta 2011 que varios autores utilizaron el método de continuación homotópica con la ayuda de herramientas (Toolbox) pertenecientes al software MATLAB®, aplicadas en la resolución de casos de estudio unidimensionales (Rahimian et al., 2011a). Poco después, en ese mismo año, se analizaron sistemas de ecuaciones algebraicas no lineales y trascendentes (Rahimian et al., 2011b).

Sin embargo, el desarrollo del seguimiento homotópico no está limitado al uso de Toolbox pertenecientes a algún software. En 2013 se desarrolló un algoritmo computacional en código FORTRAN® que utiliza hiperesferas para mejorar el seguimiento predictivo-correctivo del camino homotópico (Oliveros-Muñoz & Jiménez-Islas, 2013). Dicho método demostró ser más eficaz que otros métodos al encontrar más soluciones a sistemas con multiplicidad de soluciones previamente reportados en la literatura.

El desarrollo del algoritmo de seguimiento homotópico hiperesférico (SHH) en el lenguaje de programación FORTRAN® destacó las características del método para poder modificarse y aplicarse de acuerdo con las necesidades de cada usuario. Esto lo demostraron diversos autores, quienes analizaron el método junto con distintos métodos homotópicos y vectores de predicción, a fin de utilizar el seguimiento homotópico para tratar problemas de ingeniería química (Jiménez-Islas et al., 2013).

Posteriormente, también se utilizó el método de SHH para encontrar los puntos de operación de circuitos integrados cuya modelación es un sistema de ecuaciones algebraicas no lineales (Torres-Muñoz et al., 2016). Por otro lado, otros autores (Diaz-Arango et al. 2018) optaron por optimizar aún más el método de SHH para emplearlo en el desarrollo de algoritmos de planeación de ruta para robots móviles. Este nuevo enfoque de optimización fue utilizado ese mismo año en aplicaciones de diseño de procesos sostenibles, como la producción de dimetil-éter (Asadi & Jalali, 2018).

En otras áreas de la ingeniería se han desarrollado algoritmos para aproximar las soluciones verdaderas a problemas de campo eléctrico cuyas ecuaciones no lineales son transformadas a ecuaciones de homotopía (Tang et al., 2019). Ese mismo año, por medio del método de continuación homotópica se resolvió el problema cinemático de avance de un manipulador paralelo 3-RRS (Gallardo-Alvarado, 2019).

Se puede concluir entonces que el método de SHH es un método numérico global, aplicable a diversas áreas de la ingeniería, el cual, una vez conocida su metodología, puede ser adaptable y optimizable de acuerdo con las necesidades del usuario. Por lo tanto, es importante conocer cómo funciona el algoritmo y los diferentes aspectos que intervienen en su ejecución, explicado de una manera sencilla, desde detalles fundamentales, algebraicos, numéricos y geométricos hasta representaciones abstractas de algunas estrategias colaterales de solución.

Continuación Homotópica

Partiendo de lo anterior, la continuación homotópica es un método numérico de convergencia global, es decir, que es capaz de encontrar múltiples soluciones a sistemas de ecuaciones algebraicas no lineales mediante una aproximación inicial. El método de continuación homotópica tiene como premisa el resolver sistemas de ecuaciones algebraicas no lineales a través de la deformación continua de una función G(x) hasta llegar a convertirse en otra función continua f(x) perteneciente al sistema de ecuaciones no lineales. Este mapeo se lleva a cabo con la ayuda del parámetro homotópico τ para todas las ecuaciones del sistema (Diaz-Arango et al., 2018). De esta manera el sistema se vuelve un sistema homotópico expresado como:

Hx,τ=τfx+1-τGx (1)

donde Hfx=Hx,τ:Rn+1Rn, x Rn τ0,1.

Dada la ecuación 1, y tomando en cuenta las condiciones impuestas para Hx,τ, es conveniente definir a la ecuación Gx de manera que tenga una solución propuesta o conocida; de este modo, se cumplen las condiciones propias de la continuación homotópica:

Si τ=0 y H(x, τ)= Gx, entonces la solución propuesta es obtenida de manera muy sencilla (dada su definición).

En contraste, si τ=1 y H(x, τ)=f(x), entonces la solución del sistema original es obtenida (Oliveros-Muñoz & Jiménez-Islas, 2013).

Con respecto a estos tres primeros conceptos (ecuación 1, condición 1 y condición 2), resulta útil visualizar que el proceso de variación, llamado en lo sucesivo SHH, consiste en deformar ese sistema intencionalmente construido para que sea fácil resolver Gx hasta que se convierta en f(x), todo a través de la variación paulatina y sistemática de τ, desde 0 hasta 1. Un ejemplo geométrico sencillo para integrar este proceso de manera intuitiva es considerar la clásica deformación homotópica necesaria para obtener un toroide a partir de una taza (Figura 1).

Fuente: Elaboración propia.

Figura 1 Interpretación esquemática de la deformación de una taza para obtener un toroide al llevar τ de 0 a 1. 

Los sistemas algebraicos asociados al toroide y la taza son bastante complejos; el ejemplo es usado únicamente para establecer el concepto intuitivo de deformación geométrica asociada a la ecuación homotópica (1).

Una de las consecuencias de la condición descrita para el dominio de las variables de la ecuación 1 (Rn+1Rn, x Rn τ0,1) es la noción de hiperespacio, que denota la generalización de las propiedades del espacio euclídeo (donde operan las geometrías clásica y analítica) de tres a cuatro o más dimensiones. Téngase en consideración que todo sistema con tres o más variables (y ecuaciones) cumple con x Rn (en este caso n sería 3), y al considerar el parámetro homotópico τ0,1, se tiene entonces Rn+1Rn (que en este ejemplo es R4R3). Esta aclaración es importante porque el SHH hace uso de propiedades euclídeas de las figuras regulares, como la definición R3 de la esfera en sistemas con tres o más ecuaciones.

En sistemas algebraicos, las formas más conocidas y sencillas para Gx (Diaz-Arango et al., 2018; Oliveros-Muñoz & Jiménez-Islas, 2013) son denominadas:

  • a) Homotopía de punto fijo, si Gx=x-x0,

  • b) Homotopía de Newton, si Gx=F(x)-Fx0.

Aunque existen más formas homotópicas Gx disponibles en la literatura científica, dada la intención didáctica del presente trabajo, solo se analizará la homotopía de Newton, la cual se sustituye en la ecuación 1 para así obtener:

Hx,τ=fx-1-τfx0 (2)

Materiales y métodos

Para poder efectuar el seguimiento del camino homotópico (H-1) por medio de hiperesferas, es necesario considerar los cuatro puntos base descritos a continuación:

1) Se parte de un hiperpunto inicial x0=x10, x20, x30, , xn0t, el cual se convertirá en el centro de la primer hiperesfera de radio r fijado previamente.

2) A partir de ese punto central se genera un vector tangente al camino homotópico que intercepte la superficie de la hiperesfera en el hiperpunto xp=x1p, x2p, x3p, , xnpt (predictor).

3) Tomando este punto como aproximación, se busca la intersección de la hiperesfera con el camino homotópico a través de la resolución del sistema de ecuaciones que conforman ambos elementos con el método de Newton-Raphson. De esta manera se obtiene el hiperpunto xc=x1c, x2c, x3c, , xnct (corregido).

4) El hiperpunto xc será el centro para la generación de la siguiente hiperesfera y el punto inicial para el siguiente vector tangente, con el fin de encontrar un nuevo punto xp que conducirá a un nuevo xc, y así sucesivamente.

De acuerdo con los puntos anteriores, se puede inferir que la trayectoria del camino homotópico está delimitado por una sucesión de hiperesferas (Oliveros-Muñoz & Jiménez-Islas, 2013), tal como se muestra en la Figura 2, donde Vp corresponde a los vectores de predicción y H-1P es el camino homotópico que conforman; Vc corresponde a los vectores corregidos por medio de la corrección hiperesférica y H-1C es el camino homotópico formado por estos; CN corresponde a los puntos de corrección hiperesférica realizada por el método de Newton-Raphson y H-1 sería el camino homotópico cuya trayectoria es desconocida a priori.

Fuente: Elaboración propia.

Figura 2 Diagrama ilustrativo de seguimiento homotópico hiperesférico. 

Por otro lado, la Figura 3 muestra un diagrama de flujo del seguimiento homotópico con corrección hiperesférica detallado de la manera más simple posible, con el fin de ilustrar paso a paso el desarrollo de la metodología que se explicará más adelante.

Fuente: Elaboración propia.

Figura 3 Diagrama de flujo, paso a paso, de algoritmo de seguimiento homotópico hiperesférico (SHH). 

Deformación homotópica y Parametrización del sistema

Para generar los vectores de predicción se necesitan las pendientes de dichos vectores, estas se obtienen de las ecuaciones que componen el camino homotópico. Por otra parte, se mencionó anteriormente que el sistema de ecuaciones algebraicas no lineales se deforma por medio de la adición de una ecuación G(x) y su respectivo parámetro homotópico τ. Dicha deformación da como resultado la ecuación 1, de la cual se puede deducir:

H=fx,τ=H(x,τ) (3)

Si se considera que H(x, τ) es también función de un parámetro p, con el cual se estima la distancia recorrida a lo largo del espacio homotópico, entonces todas las variables implicadas (xRn y τ 0,1) tienen funcionalidad de p.

Hp=Hxp,τp=h[x1p,x2p,x3p,,xnp,τp] (4)

De esta manera se parametriza el sistema.

Definición de matriz Jacobiana

Una vez que se tiene el sistema homotópico, se procede a agrupar el vector de x y el parámetro homotópico τ en un vector auxiliar y.

Hyp=0 (5)

ypH-1

Diferenciando la ecuación anterior, y empleando la regla de la cadena, se obtiene:

i=1n+1Hyi  yi˙=0 (6)

donde Hyi es una columna del Jacobiano H' del sistema homotópico, por lo que se tiene:

H1x1H1x2H1x3H1xnH1τH2x1H2x2H2x3H2xnH2τH3x1Hnx1H3x2Hnx2H3x3Hnx3H3xnHnxnH3τHnτx1px2px3pxnpτp=0 (7)

De esta manera se obtiene un sistema de n ecuaciones con n+1 variables. Al tratarse de un sistema homogéneo donde el número de variables es mayor al número de ecuaciones, existen soluciones diferentes a la trivial.

Cálculo de pendientes aplicando la Regla de Cramer

El sistema anterior (7) se puede resumir en notación vectorial como:

H'xp,τpy˙=0 (8)

H'yy˙=0

Desarrollando algebraicamente la ecuación anterior para un ejemplo genérico de dos incógnitas se aplica la regla de Cramer, haciendo:

x1˙=x1p

x2˙=x2p

τ˙=τp (9)

De manera que, usando movimientos de determinantes, se puede hacer que xi˙ quede en función de τ˙ y pueda expresarse de la siguiente forma:

x1˙=-1-1τ˙D_1D_3

x2˙=-1-1τ˙D_2D_3 (10)

donde D_i pertenece al determinante del Jacobiano H' en el cual se ha eliminado la columna i. Como se tienen dos ecuaciones de tres incógnitas, se propone:

τ˙=(-1)D_3 (11)

Estas ecuaciones se transforman como se muestra a continuación:

x1˙=(-1)D_1

x2˙=D_2

τ˙=(-1)D_3 (12)

Estas se pueden generalizar como:

yi˙=-1iD_i=xipτp (13)

La expresión anterior representa los valores de las pendientes del camino homotópico, necesarias para calcular los vectores predictores y, por lo tanto, la trayectoria del camino homotópico H-1.

Cálculo de vector predictor

Para realizar el cálculo de los vectores de predicción se toman en cuenta las pendientes del camino homotópico mencionadas en el punto anterior y el algoritmo de Euler, que en este caso corresponde a:

xip=xi0+pxip (14)

τp=τ0+pτp (15)

donde xip y τp son los valores de predicción obtenidos por el método de Euler que se encontrarán en la superficie de la hiperesfera de radio r. Por otro lado, para calcular el avance horizontal p se hace uso de las ecuaciones 14 y 15 en conjunto con la ecuación de la hiperesfera:

i=1nxi-xi02+τ-τ02-r2=0 (16)

Por lo tanto, p se calcula con la siguiente fórmula:

p=±r2i=1nxip2+τp212 (17)

donde el signo dependerá de si el seguimiento se realiza a la derecha o a la izquierda, es decir, en sentido positivo o negativo. De esta manera es posible calcular los vectores de predicción para x y τ de acuerdo con las ecuaciones 14 y 15.

Para calcular el punto p correspondiente al avance del camino homotópico se utilizará la ecuación:

p=p0+p (18)

Corrección hiperesférica

Una vez calculados los vectores de predicción, estos se vuelven el valor inicial para empezar con la corrección hiperesférica, dicho procedimiento tiene como objetivo el encontrar la intersección entre la hiperesfera y el camino homotópico.

Este procedimiento se realiza empleando el método de Newton-Raphson para resolver las ecuaciones homotópicas del sistema, descritas como la ecuación 2 en conjunto con la ecuación de la hiperesfera, ecuación 16.

Recordemos que el método de Newton-Raphson requiere una matriz de derivadas parciales, denominado Jacobiano (J), y un vector con la función evaluada en un punto inicial (B) para comenzar la corrección. Esto puede verse de forma más clara con la expresión siguiente:

JX=B

J-1B=X (19)

donde los valores del vector X se sumarán a los valores del vector predictor para obtener un vector corrector:

Xp+X=Xc (20)

Este proceso también es iterativo, se recomienda un rango de tres a cinco iteraciones para obtener un buen resultado.

Por lo tanto, para un sistema de dos ecuaciones, el Jacobiano (J) quedaría expresado de la siguiente forma:

J=H1x1H1x2H1τH2x1H2x2H2τcx1cx2cτ (21)

donde c es la ecuación de la hiperesfera derivada parcialmente en cada uno de los componentes del sistema. De manera similar, el vector denominado B, que contiene las funciones homotópicas evaluadas en los valores del vector de predicción y la función de la hiperesfera, quedará expresado como:

B=-H1-H2-c (22)

De esta manera el punto corregido se convertirá en el centro de la siguiente hiperesfera y, a su vez, en el punto inicial del siguiente vector de predicción.

Detección y refinamiento de cruces con τ=1

De acuerdo con los puntos anteriores, para que el camino homotópico pueda trazarse, es necesario un alto número de iteraciones. Estas iteraciones tienen que monitorearse para detectar el momento en el que el camino homotópico cruce el hiperplano τ=1, ya que, de acuerdo con la teoría, este será el punto en el que se encuentre el vector solución del sistema.

Sin embargo, es casi imposible que el cruce del camino se localice en el punto exacto de τ=1. Por lo tanto, cuando se detecte el paso de un número inferior (τ<1) a un número superior a (τ>1), o viceversa, es cuando hay que realizar el refinamiento de cruces con τ=1.

Este proceso se llevará a cabo con el mismo método que la corrección hiperesférica, es decir, se empleará el método de Newton-Raphson. En esta ocasión se contará solamente con las funciones del sistema de ecuaciones algebraicas no lineales originales y el último vector corregido que haya sobrepasado a τ=1, esto aplica tanto en el sentido +τ como en la dirección contraria -τ.

Para un sistema de dos ecuaciones se expresaría como:

JR-1BR=XR=f1x1f1x2f2x1f2x2-1-f1-f2 (23)

donde el subíndice R indica que se trata de elementos usados en el refinamiento de cruces con τ=1. Finalmente, se recomienda un número mínimo de cinco iteraciones para obtener el vector solución lo más exacto posible.

El pseudo-código del algoritmo de SHH puede encontrarse en el Apéndice B de este mismo artículo.

Ejemplo práctico

Deformación homotópica y parametrización del sistema:

En un ejemplo práctico considere el sistema de dos ecuaciones:

f1=x12+x22-25

f2=x1+x2-7 (24)

con un vector inicial x10,x20,τ0,r=[1,1,0,0.1].

De acuerdo con la metodología del seguimiento homotópico, se debe deformar el sistema introduciendo la ecuación de homotopía de Newton y reescribirse de acuerdo con la ecuación 2 como:

H1(x1,x2,τ)=[x12+x22-25]-(1-τ)[x102+x202-25]

H2(x1,x2,τ)=[x1+x2-7]-(1-τ)[x10+x20-7] (25)

donde xi_0 representan los valores del vector inicial.

NOTA: Estos valores xi_0 permanecen invariables mientras se trate de una función deformada, es decir, una función homotópica (Hi).

Definición de matriz Jacobiana

Tomando en cuenta lo anterior, se procede a obtener el Jacobiano (H') del sistema derivando parcialmente las ecuaciones homotópicas. De esta manera se obtiene:

H'yy˙=H1x1H1x2H1τH2x2H2x2H2τ-1x1px2pτp=0 (26)

Sustituyendo el Jacobiano H' del sistema se obtiene:

2x12x2x1_02+x2_02-2511x1_0+x2_0-7-1 (27)

Sustituyendo los valores iniciales:

2(1)2(1)12+12-25111+1-7-1=22-2311-5-1

Cálculo de pendientes aplicando la Regla de Cramer

En este paso se tiene por objetivo obtener los determinantes de cada uno de los elementos del sistema, (xi,τ). Este proceso se puede resumir a la ecuación 13, la cual expresa el cálculo de los determinantes de matrices a las cuales se les ha retirado la columna i perteneciente a la variable en turno. Esto es:

x1˙=2-231-5detx1˙=13

x2˙ =2-231-5 detx2˙=13

τ˙ =2211detτ˙=0

Aplicando la ecuación 13 para el cambio de signo, los determinantes quedan expresadas como:

yi˙=-13130 (28)

Cálculo de vector predictor

El vector predictor se calculará usando las pendientes que hemos calculado en el punto anterior, es decir, y˙i:

yi˙=x1px2pτp=-13130

Las ecuaciones necesarias para calcular el vector de predicción son las ecuaciones 14 y 15, las cuales requieren un elemento p que se calcula con la ecuación 17. Sustituyendo los valores de xip, τp y r en la formula, la ecuación se reduce a:

p=±0.12-132+132+0212=±0.0005439

Posteriormente, se calcula el vector de predicción con las ecuaciones 14 y 15. Esto da como resultado:

x1p=x10+px1p=1+0.0005439-13=0.9929

x2p=x20+px2p=1+0.000543913=1.0070

τp=τ0+pτp=0+0.00054390=0

Por último, se calcula el avance de p.

p=p0+p=0+0.0005439=0.0005439

Corrección hiperesférica

Para llevar a cabo la corrección hiperesférica se calculará un Jacobiano (J) de acuerdo con la ecuación 21, detallada previamente. En dicho Jacobiano se sustituirán los valores pertenecientes al vector predictor y al vector inicial de la siguiente manera:

J=2x1_p2x2_px1_02+x2_02-2511x1_0+x2_0-72(x10-x1p)2(x20-x2p)2(τ0-τp) (29)

donde xi_0 permanecerá constante, mientras que xi0 cambiará en cada iteración de la corrección con el punto anterior. Por otro lado, el vector B, que contiene las funciones homotópicas evaluadas en el punto predictor y la ecuación de la hiperesfera, se definirá de acuerdo con la forma de la ecuación 22.

B=-[x1_p2+x2_p2-25]-(1-τp)[x102+x202-25]-[x1_p+x2_p-7]-(1-τp)[x10+x20-7]-x10-x1p2+x20-x2p2+τ0-τp2-r2

Una vez definida la matriz del Jacobiano (J) y el vector B, se sustituyen los valores de manera que xi0, xi_0, τ0, xi_p, xip y τp contengan los valores del vector inicial y el vector de predicción.

J=1.9858578642.014142136-2311-50.014142136-0.0141421360

B=-1E-0400.0099

El siguiente paso es obtener la matriz inversa de J, es decir J-1:

J-1=-0.1923076920.88461538535.16303137-0.1923076920.884615385-35.54764675-0.0769230770.153846154-0.076923077

Dicha matriz se multiplica por el vector B para obtener un vector X que se sumará al vector de predicción para obtener un nuevo vector X.

X=0.348133241-0.351902472-0.000753846+0.9929289321.0070710680=1.3410621730.655168596-0.000753846

En la siguiente iteración el vector X ocupará el lugar del vector de predicción hasta que se cumplan las iteraciones del proceso de corrección hiperesférica, se recomienda que sea de tres a cinco iteraciones. Finalmente, el último vector X será el vector corregido Xc que pasará a ser el punto inicial xi0 para la siguiente iteración global del proceso de continuación homotópica (no confundir con xi_0).

Detección y refinamiento de cruces con τ=1

El proceso iterativo de continuación homotópica anteriormente explicado fue realizado durante cuatro ciclos cuando se registró un cruce en τ=1 donde τ cambia de 0.9636 a 1.0540. Dado que no se trata del punto exacto τ=1, se procede a refinar la solución utilizando el último punto corregido xi_c como punto inicial para el método de Newton-Raphson descrito por la ecuación 23, donde la matriz Jacobiana correspondiente al proceso de refinado JR es igual a:

JR=f1x1f1x2f2x1f2x2=2x1_c2x2_c11=2(3.3062)2(3.9637)11=6.612567.9275611

y el vector BR corresponde a:

BR=-f1-f2=-x1_c2+x2_c2-25-x1_c+x2_c-7=-1.64305-0.27006

Posteriormente, se calculó la inversa de JR y se multiplicó por el vector BR para obtener el vector XR de acuerdo con la ecuación 23 que se mencionó anteriormente.

JR-1=-0.7604536.0285440.760453-5.0285445

JR-1BR=XR=-0.7604536.0285440.760453-5.0285445-1.64305-0.27006=-0.3786130.108551

Dicho vector XR se sumó al vector inicial para obtener el nuevo vector XR.

XR=-0.3786130.108551+3.30623.9637=2.9276664.072333

El nuevo vector XR pasará a ser el vector inicial para la segunda iteración en el proceso de refinamiento de cruces en τ=1. Este proceso se repitió hasta que los valores del vector BR, que contiene las funciones del sistema original evaluadas con el vector inicial en turno, se volvieron prácticamente 0 (Es decir, los valores están por debajo de la tolerancia fijada de antemano, por ejemplo 10-5.) De esta manera es como se obtuvo el primer vector solución del sistema, para este caso particular es de:

xi=34

En este caso de estudio se realizó un seguimiento de camino homotópico corto con 12 hiperesferas, las cuales fueron suficientes para encontrar los dos vectores solución conocidos para este sistema. Dichos vectores se muestran en la Tabla 1.

Tabla 1 Vectores solución del sistema. 

x1 x2
3 4
4 3

Fuente: Elaboración propia.

Todos los vectores corregidos obtenidos durante el seguimiento hiperesférico se muestran en la Tabla 2, donde los cruces en τ=1 están resaltados en color gris oscuro. Se registra un total de tres cruces que fueron refinados con el método de Newton-Raphson, sin embargo, los dos primeros cruces resultaron ser el mismo vector solución, por lo que se omitió en la Tabla 1. En el Apéndice A, la Figura A1 muestra las gráficas 3D de las superficies H1(x1,x2,τ) y H2(x1,x2,τ) en un momento inicial (τ=0), y la Figura A2 muestra las seis primeras deformaciones y su influencia sobre las intersecciones de las funciones homotópicas con H=0. En este caso, la deformación homotópica se representa por el cruce simultaneo de las dos superficies con el plano H=0, y la deformación homotópica implica la expansión de ambas superficies. En sistemas con más ecuaciones y variables, este tipo de visualización geométrica se dificulta; sin embargo, dada la generalización de las propiedades euclídeas del espacio R3 a Rn, numéricamente pasa algo que tiene propiedades geométricas similares.

Tabla 2 Vectores corregidos obtenidos durante el seguimiento del camino homotópico. 

x1 x2 τ p
1.0 1.0 0.0 0.0
2.45000 0.36404 0.16281 0.00054
3.81500 1.10509 0.58402 0.00105
4.20730 2.61113 0.96369 0.00155
3.30628 3.96378 1.05401 0.00206
1.71130 4.14636 0.77153 0.00257
0.55608 3.12339 0.33589 0.00307
0.51537 1.54971 0.01302 0.00357
1.71810 0.44395 0.03241 0.00409
3.28956 0.63835 0.38558 0.00459
4.18549 1.89260 0.81562 0.00509
3.85895 3.45418 1.06263 0.00560

Fuente: Elaboración propia.

Finalmente, el trazo del camino homotópico se muestra en la Figura 4, donde se graficó p en el eje horizontal y τ en el eje vertical. El camino homotópico se muestra en color azul y una línea ilustrativa de τ=1 en color rojo, esta línea sirve para detectar a simple vista los tres cruces que se mencionaron anteriormente.

Fuente: Elaboración propia.

Figura 4 Gráfico del camino homotópico para el caso ejemplo de dos ecuaciones. 

Resultados

El caso ejemplo fue trabajado y graficado en Excel® debido a la simplicidad del sistema, los casos de estudio expuestos a continuación se comprobaron por medio de un algoritmo de seguimiento homotópico hiperesférico en FORTRAN® (Oliveros-Muñoz & Jiménez-Islas, 2013). El equipo utilizado para llevar a cabo la resolución de los sistemas fue una computadora con un procesador Intel(R) Core (TM) i7-7700HQ @2.80GHz con 8GB de RAM.

Caso de estudio 1

Consideremos el caso de estudio utilizado por Malinen & Tanskanen (2010). Se trata de un sistema de dos ecuaciones algebraicas no lineales, ecuación 30, cuyas soluciones se encuentran situadas entre -4x1,x24. Para este sistema se encuentran reportados un total de nueve vectores solución encontrados con el método de homotopía con un límite en su parámetro homotópico.

f1=4x13+4x1x2+2x22-42x1-14

f2=4x23+4x1x2+2x12-26x2-22 (30)

El programa de seguimiento homotópico hiperesférico (SHH) fue capaz de encontrar el mismo número de vectores solución (Vn), los cuales se muestran en la Tabla 3. Este resultado se obtuvo con un vector inicial x10,x20,τ0,r0=1,1,0,0.1, el tiempo de cómputo utilizado fue de 16 s.

Tabla 3 Vectores solución para el caso de estudio 1. 

Vector solución V1 V2 V3 V4 V5 V6 V7 V8 V9
X1 0.0866 -2.8051 -3.0730 -3.7793 -0.2708 -0.1279 3.5844 3.3851 3
X2 2.8842 3.1313 -0.0813 -3.2831 -0.9230 -1.9537 -1.8481 0.0738 2

Nota: Los datos obtenidos fueron comparados con los obtenidos por (Malinen & Tanskanen, 2010).

Fuente: Elaboración propia.

Gracias a la parametrización del sistema homotópico con respecto a p, es posible graficar el camino homotópico, tal como se mostró en el caso ejemplo. Así pues, para este caso de estudio el camino homotópico se muestra en la Figura 5, donde p recorrió el espacio de p tanto en sentido negativo como en positivo. El sentido negativo se encuentra marcado de color rojo, mientras que el sentido positivo se muestra de color azul. La línea de color verde representa τ=1, de manera que se puedan detectar a simple vista los nueve cruces que representan los nueve vectores solución.

Fuente: Elaboración propia.

Figura 5 Camino homotópico del caso de estudio 1. 

De acuerdo con la definición del método de SHH como un método de convergencia global, que tiene como premisa encontrar la mayor cantidad de vectores solución con una sola aproximación inicial, los resultados presentados anteriormente son los esperados comparados con otros métodos de convergencia local.

A continuación, se muestra un breve análisis de este caso de estudio realizado por un método numérico basado en el concepto de continuación homotópica (BCH) (Martín, 2013). Dicho método deforma el sistema original en un sistema homotópico de tal forma que se transforma en un problema de valor inicial; posteriormente, se utiliza el algoritmo de Runge-Kutta de 4° orden para hacer la búsqueda del vector solución en el intervalo [0<τ<1] (Apéndice C). Para complementar la evaluación del caso de estudio, también se utilizó el método de Newton-Raphson para sistemas de ecuaciones no lineales.

Los resultados de los dos métodos anteriores, en conjunto con los resultados del método de SHH, se muestran en la Tabla 4. Los tres métodos fueron evaluados con diferentes vectores iniciales en un equipo computacional con procesador Intel® Xeon® E5-2620 V3 @2.4GHz con 64GB de RAM. Los resultados incluyen los vectores solución encontrados por cada método, numerados como se muestra en la Tabla 3, así como el tiempo de cómputo utilizado durante la evaluación y la proximidad del vector solución obtenido con la solución en fx=0.

Tabla 4 Análisis de caso de estudio 1 con métodos de convergencia local y global. 

Vector
Inicial
(x10,x20)
Método Solución
Encontrada (Vi)
Tiempo de
cómputo
Tolerancia con
respecto a
f(x)=0
(-4,-4) BCH V4 0.62 s 6.90E-06
(-4,-4) Newton-Raphson V4 0.1 s 1.20E-14
(-4,-4) SHH V4, V3, V2, V1, V5, V6 6 s -5.70E-10
(0,0) BCH V5 0.62 s 2.30E-06
(0,0) Newton-Raphson V5 0.1 s 4.00E-15
(0,0) SHH V1, V2, V3, V4, V5, V6, V7, V8, V9 24 s -2.20E-08
(2,2) BCH V9 0.64 s 4.90E-06
(2,2) Newton-Raphson V7 0.1 s -1.00E-15
(2,2) SHH V1, V2, V3, V9 6 s 3.00E-08

Fuente: Elaboración propia.

Caso de estudio 2

El siguiente caso de estudio es un sistema de tres ecuaciones algebraicas no lineales con cuatro vectores solución reportados en la literatura (Kuno & Seader, 1988; Rahimian et al., 2011b). Dichos vectores fueron encontrados por los autores utilizando un método de homotopía global de punto fijo.

F1=(x1-x22)(x1-Sen(x2))

F2=Cosx2-x1x2-Cosx1

F3=(x2)(x2-1)+x32  (31)

En el método de SHH se utilizó un vector inicial de: x10,x20,x30,τ0,r0=[0,0,0,0,0.1]. El camino homotópico se muestra en la Figura 6; sin embargo, debido a la cercanía de los cruces, es imposible detectarlos a simple vista, por lo que se recurrió a varios acercamientos.

Fuente: Elaboración propia.

Figura 6 Camino homotópico del caso de estudio 2. 

La Figura 7a muestra un acercamiento a la sección del camino homotópico en sentido negativo que atraviesa la línea transversal de τ=1, mientras que en la Figura 7b se muestra un acercamiento a la sección del sentido positivo que atraviesa la misma línea. Cada una de estas líneas representan dos cruces, uno de forma ascendente y otro de forma descendente.

Fuente: Elaboración propia.

Figura 7 a) Sección del camino homotópico en sentido negativo; b) Sección del camino homotópico en sentido positivo. 

Finalmente, la Figura 8 muestra un acercamiento mucho más cercano a los cruces en la sección en sentido positivo del camino homotópico. Esto es para demostrar cómo el camino homotópico cruza ascendentemente τ=1 para luego descender generando otro cruce. Tomando en cuenta esto, hay un total de ocho cruces del camino homotópico con τ=1, cuatro más que los reportados en la literatura.

Fuente: Elaboración propia.

Figura 8 Sección en sentido positivo del camino homotópico que muestra dos cruces en τ=1. 

Los vectores solución correspondientes a cada uno de los cruces antes mencionados se muestran en la Tabla 5, señalando cuáles fueron los vectores solución adicionales encontrados. Este proceso llevó un tiempo de cómputo de 25 s.

Tabla 5 Vectores solución encontrados y reportados para el caso de estudio 2. 

x1 x2 x3 Vector solución reportado
0.70710678 0.78539816 0.41054584
0.67919407 0.82413231 0.38070756
0.64171437 0.80107077 0.39919468
0.69481969 0.76816915 0.42200155
0.70710678 0.78539816 -0.4105458 No
0.67919407 0.82413231 -0.3807075 No
0.64171437 0.80107077 -0.3991946 No
0.69481969 0.76816915 -0.4220015 No

Nota: Los datos obtenidos se compararon con los obtenidos por (Kuno & Seader, 1988; Rahimian et al., 2011b).

Fuente: Elaboración propia.

Caso de estudio 3

En un caso de estudio aplicado en el área de ingeniería química, considere el evaporador de triple efecto mostrado en la Figura 9. Dicho evaporador tiene por objetivo realizar la concentración de una solución de azúcar del 10% al 50% en peso. La alimentación de la primera etapa usa un vapor saturado a 280 °F y produce líquido a 100 °F a un flujo másico de 50000 lb/h. La tercera etapa debe operarse a una presión absoluta, de manera que el punto de ebullición de la solución sea de 125 °F. Descarte el efecto del incremento del punto de ebullición en las variaciones del calor específico y el calor latente con respecto a la temperatura y la concentración. Determine el área de transferencia para cada etapa (se consideran iguales), las temperaturas (T1,T2), los flujos (L1,L2,L3), las componentes (x1,x2) y el flujo V0 del vapor de alimentación. Tome Cp=1 BTU/lb°F y un calor latente λ=1000 BTU/lb con coeficientes de transferencia de calor U1=500 BTU/ft2h °F, U2=300 BTU/ft2h °F y U3=200 BTU/ft2h °F (Oliveros-Muñoz & Jiménez-Islas, 2013) .

Fuente: Oliveros-Muñoz & Jiménez-Islas (2013).

Figura 9 Esquema de un evaporador de triple efecto. 

De un balance de materia global se obtiene que L3 = 10000 lb/h, por lo que, planteando balances macroscópicos de materia, energía y la transferencia de calor regida por la ley de enfriamiento de Newton y sustituyendo los valores conocidos, resulta el siguiente sistema de ecuaciones no lineales:

f1=50000(100-x1)+1000x2-1000(50000-x4)

f2=500x3(250-x1)-1000x2

f3=300x3*(x1-x5)-1000(50000-x4)

f4=x4x1-x5+1000(50000-x4)-1000(x4-x6)

f5=x6(x5-125)+1000(x4-x6)-1000(x6-10000)

f6=200x3(x5-125)-1000(x4-x6) (33)

en donde:

x 1 = T 1 =

Temperatura de ebullición del primer efecto, °F

x 2 = V 0 =

Suministro de vapor, lb/h

x 3 = A =

Área de transferencia de calor en cada efecto, ft 2

x4 = L1 =

Flujo másico de concentrado del primer efecto, lb/h

x 5 = T 2 =

Temperatura de ebullición del segundo efecto, °F

x6 = L2 =

Flujo másico de concentrado del segundo efecto, lb/h

Para realizar el SHH se utilizó el vector inicial x10,x20,x30,x40,x50,x60,τ0,r0=0,0,10,10,100,100,0,0.1 con un tiempo de cómputo de 20 s.

El camino homotópico se muestra en la Figura 10. Se puede observar un total de cinco cruces, sin embargo, cuando se refinan, se muestra que solo tres de ellos son vectores solución que no se repiten. Dichos vectores solución se muestran en la Tabla 6. Para una mejor visualización del camino homotópico se consideró un factor de multiplicación en p de 1X1025. El vector solución factible termodinámicamente es V3, que no es fácil de encontrar con otros métodos como el de Newton-Raphson.

Fuente: Elaboración propia.

Figura 10 Camino homotópico para el caso de estudio 3. 

Tabla 6 Vectores solución encontrados en el caso de estudio 3. 

Vector solución V1 V2 V3
x1 -2051.32 -1232.79 218.5346
x2 -73 749.84 -34 398.43 17 888.59
x3 -64.0936 -46.397 1137.03
x4 16 183.98 17 759.04 38 038.14
x5 -292.6381 1083.52 183.467
x6 10 830.40 26 653.57 24 742.38

Nota: Los datos obtenidos se compararon con los obtenidos por Oliveros-Muñoz & Jiménez-Islas (2013).

Fuente: Elaboración propia.

Al igual que el caso de estudio 1, este caso de estudio se evaluó con el método BCH y Newton-Raphson multivariable en igualdad de condiciones, los resultados se muestran en la Tabla 7. Dichos resultados muestran los vectores iniciales utilizados, los vectores solución encontrados por cada método, numerados de acuerdo con la Tabla 6, mientras que la precisión de cada vector solución con respecto a fx=0 se muestra en la parte inferior de la Tabla.

Tabla 7 Análisis de caso de estudio 3 con métodos de convergencia local y global. 

Vector inicial Método Solución encontrada (Vi) Tiempo de cómputo
(0.1, 0.1, 1, 1, 100, 100) BCH V2 2.3 s
Newton-Raphson V1 6.4 s
SHH V1 , V2 , V3 15 s
(0.3, 0.3, 3, 3, 300, 300) BCH V2 2.3 s
Newton-Raphson V1 0.1 s
SHH V1 , V2 , V3 10 s
(0.4, 0.4, 4, 4, 400, 400) BCH V2 2.3 s
Newton-Raphson V1 0.1 s
SHH V1 , V2 , V3 12 s
(0.6, 0.6, 6, 6, 600, 600) BCH V2 2.3 s
Newton-Raphson V2 0.1 s
SHH V1 , V2 , V3 8 s
(0.7, 0.7, 7, 7, 700, 700) BCH V2 2.3 s
Newton-Raphson V3 0.1 s
SHH V1 , V2 , V3 11 s
i=1nabsfix1,x2,,xn0.1  y i=1nabsfix1,x2,,xn1×10-7

Fuente: Elaboración propia.

Caso de estudio 4

El siguiente caso describe el problema cinemático que presenta un manipulador paralelo 3-RRS. Este sistema consta de dos plataformas triangulares, la plataforma inferior es fija mientras que la superior es móvil; ambas están unidas por tres brazos con tres articulaciones, uno en cada una de las plataformas y otro en medio del brazo de unión, el sistema se muestra en la Figura 11.

Fuente: (Gallardo-Alvarado, 2019).

Figura 11 Esquema de manipulador paralelo 3-RRS. 

El problema consiste en encontrar las posiciones alcanzables de la plataforma móvil que cumple con un conjunto prescrito de coordenadas generalizadas, para este caso son q1=120°,  q2=135° y q3=95°. Dichas coordenadas están relacionadas a la plataforma fija del manipulador paralelo, por lo tanto, el sistema se encuentra descrito por el siguiente sistema de ecuaciones algebraicas no lineales (Gallardo-Alvarado, 2019):

F1=x12-1.37x1+x22+x32-1.212x3+0.416

F2=x2

F3=x12-2x1x4+x42+x22-2x2x5+x52+x32-2x3x6+x62-0.25

F4=-0.866x4-0.5x5

F5=x42+0.831x4+x52-1.439x5+x62-0.989x6+0.513

F6=x42-2x4x7+x72+x52-2x5x8+x82+x62-2x6x9+x92-0.25

F7=0.866x7-0.5x8

F8=x72+0.397x7+x82+0.688x8+x92-1.394x9+0.221

F9=x12-2x1x7+x72+x22-2x2x8+x82+x32-2x3x9+x92-0.25 (34)

Para resolver el sistema por el método de SHH se utilizó un vector inicial xi0,τ0,r0=0.1,0,2.5,0.5,1.5,1,0.5,-0.1,1.5,0,0.1. Los vectores solución encontrados tomaron un tiempo de cómputo de 1.5156 s, fueron un total de seis vectores, los cuales se presentan en la Tabla 8 y coinciden con el número reportado en la literatura.

Tabla 8 Vectores solución encontrados para el caso de estudio 4. 

Vector
solución
x1 x2 x3 x4 x5 x6 x7 x8 x9
V1 0.05224 0.00000 0.74768 -0.20478 0.35468 0.98882 -0.02069 -0.03584 1.24103
V2 0.29758 0.00000 0.08603 -0.10100 0.17493 0.33206 -0.13938 -0.24140 0.05784
V3 0.29085 0.00000 0.09111 -0.14152 0.24512 0.14557 -0.14256 -0.24692 0.05669
V4 0.27404 0.00000 0.10443 -0.15146 0.26233 0.11618 0.11607 0.20103 0.53412
V5 0.05145 0.00000 0.46788 -0.15643 0.27094 0.10267 -0.13088 -0.22668 0.06122
V6 0.18985 0.00000 1.02468 -0.18599 0.32213 0.95415 0.12074 0.20912 0.57580

Nota: Los datos obtenidos se compararon con los obtenidos por (Gallardo-Alvarado, 2019).

Fuente: Elaboración propia.

En la Figura 12a se muestra el camino homotópico del problema. Al igual que en problemas anteriores, se presenta el camino en sentido positivo en color azul y en sentido negativo en color rojo, además de una línea τ=1 donde se puede ver una gran cantidad de cruces. La Figura 12b muestra un acercamiento en una sección del camino homotópico en sentido positivo, de esta manera se puede ver con más claridad que, a pesar de tratarse de una gran cantidad de cruces, todo se reduce a una repetición periódica de seis cruces, cuya longitud de onda va en aumento conforme el camino se desplaza a la derecha o a la izquierda. Dicho de otra manera, se trata de los mismos vectores solución presentados con anterioridad y repetidos con frecuencia uno detrás de otro.

Fuente: Elaboración propia.

Figura 12 a) Camino homotópico del caso de estudio 4; b) Acercamiento a una sección del camino homotópico en sentido positivo. 

Al igual que en los casos de estudio 1 y 3, este caso de estudio se evaluó con el método BCH y Newton-Raphson multivariable en igualdad de condiciones, los resultados se muestran en la Tabla 9. En dicha Tabla se muestra el vector inicial, los vectores solución encontrados por cada método, numerados de acuerdo con la Tabla 8, mientras que la precisión de cada vector solución con respecto a fx=0 se muestra en la parte inferior de la Tabla.

Tabla 9 Análisis de caso de estudio 4 con métodos de convergencia local y global. 

Vector inicial Método Solución encontrada
(Vi)
Tiempo de cómputo
(-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5) BCH No converge ---
Newton-Raphson V1 0.1 s
SHH V1 , V2 , V3 , V4 , V5 , V6 15 s
(-0.4,-0.4,-0.4,-0.4,-0.4,-0.4,-0.4,-0.4,-0.4) BCH No converge ---
Newton-Raphson V3 0.1 s
SHH V1 , V2 , V3 , V4 , V5 , V6 17 s
(0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2) BCH No converge ---
Newton-Raphson V6 0.1 s
SHH V2 , V3 17 s
(0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5) BCH No converge ---
Newton-Raphson V2 0.1 s
SHH V2 , V3 , V5 , V6 16 s
i=1nabsfix1,x2,,xn1×10-7

Fuente: Elaboración propia

Caso de estudio 5

A continuación, se presenta una reacción de esterificación catalizada por ácido:

ftalato de monobutilo + n-butanol -> ftalato de dibutilo + agua

la cual se lleva a cabo continuamente en dos reactores de tanque de agitado continuo (CSTR) conectados en serie (Figura 13) (Seader et al., 1990).

Fuente: (Seader et al., 1990).

Figura 13 Diagrama de flujo de producción de ftalato de dibutilo. 

El equilibrio de vapor / líquido se establece en el reactor con ventilación de vapor. El primer CSTR es alimentado con una cantidad conocida de ácido, y una cantidad adicional de ácido ingresa al segundo CSTR para lograr una concentración de ácido específica en su efluente. Las ecuaciones de modelado consisten en cuatro ecuaciones de equilibrio de fase y tres ecuaciones de equilibrio de materia, la versión polinómica simplificada de dichas ecuaciones se muestra a continuación (Seader et al., 1990):

f1=5.5(x4-x2)(x2+x3)-x2(x1-x2-x3)

f2=3.5(0.888x1-x3-x4)(x2+x3)-x3(x1-x2-x3)

f3=5.5(0.0986x1-0.011.098x1-x2-9x4-x5-x6+x7-x2+x5)(x5+x6)-x5(1.098x1-x2-9x4-x5-x6+x7)

f4=3.5(0.986x1-9x4-x6-0.0986x1+0.011.098x1-x2-9x4-x5-x6+x7)(x5+x6)-x6(1.098x1-x2-9x4-x5-x6+x7)

f5 = 2056.40.0986x1-x42-x4x1-x2-x32

f6=0.61177-(0.0986x1-x4)+0.01(1.098x1-x2-9x4-x5-x6+x7)

f7 = 0.04(74.12(0.986x1-10x4)+222.24(0.0986x1-x4)+18(x4-x2)+278.84x4+98.090.0136x1)/98.01-0.0136x1-x7 (35)

Este sistema cuenta con un total de 44 soluciones, de las cuales 35 son soluciones complejas, de las nueve soluciones restantes solo una satisface los requerimientos físicos del sistema de acuerdo con la literatura (Gupta, 1995; Seader et al., 1990).

Para resolver el sistema por el método de SHH se utilizó un vector inicial xi0,τ0,r0=50, 5, 25, 5, 0.5, 4, 0.05, 0, 0.5 con el cual se encontraron cuatro vectores solución, incluyendo el vector que cumple con los requerimientos del sistema (V1) (Tabla 10). El tiempo de cómputo utilizado fue de 23 s.

Tabla 10 Vectores solución encontrados para el caso de estudio 5. 

Vector
solución
V1 V2 V3 V4
x1 80.539627 84.319221 25.133557 23.065198
x2 6.973011 7.283761 0.000000 0.000000
x3 61.123339 64.007810 0.000000 0.000000
x4 7.204681 7.525764 1.745845 1.625707
x5 0.547728 0.000000 0.000000 2.046727
x6 3.653712 0.000000 0.000000 5.123147
x7 0.059712 0.067266 0.171321 0.150819

Nota: Los datos obtenidos se compararon con los obtenidos por (Seader et al., 1990) (Gupta, 1995).

Fuente: Elaboración propia.

En la Figura 14 se puede observar el camino homotópico trazado para este caso de estudio. El camino en sentido positivo se mantiene en color azul mientras el negativo es de color rojo. En la línea verde perteneciente a τ=1 se puede observar un total de cuatro cruces correspondientes a los vectores solución antes mencionados.

Fuente: Elaboración propia.

Figura 14 Camino homotópico para el caso de estudio 5. 

Caso de estudio 6

El siguiente caso de estudio consiste en calcular el punto de burbuja (x5=T) de una mezcla binaria no ideal: 20%(mol) isobutanol-80% agua. En dicha mezcla el punto de burbuja forma dos fases líquidas y se puede encontrar resolviendo el sistema de ecuaciones descrito a continuación (Shacham et al., 1998):

f1=x1-0.2/(x6+1-x6k11/k12)

      f2=x2-x1k11/k12

      f3=x3-0.8/(x6+1-x6k21/k22)

      f4=x4-x3k21/k22

      f5=x1(1-k11)+x3(1-k21)

      f6=(x1-x2)+(x3-x4) (36)

donde A=1.7 y B=0.7 se utilizan para calcular los siguientes parámetros:

      γ11=10(Ax32/(Ax1B+x32))

      γ21=10(Bx12/(x1+Bx3A2))

      γ12=10(Ax42/(Ax2B+x42))

      γ22=10(Bx22/(x2+Bx4A2))

      q1=10(7.62231-1417.9/(191.15+x5))

      q2=10(8.10765-1750.29/(235+x5))

      k11=γ11q1/760

      k21=γ21q2/760

      k12=γ12q1/760

      k22=γ22q2/760

El vector inicial para resolver este sistema por medio de SHH fue xi0,τ0,r0=0, 1, 1, 4,100,0.8,  0, 0.1, con el cual se encontró el vector solución del sistema en un tiempo de cómputo de 27 s (Shacham et al., 1998) (Tabla 11).

Tabla 11 Vector solución encontrado para el caso de estudio 6. 

Vector solución x1 x2 x3 x4 x5=T x6
0.02270 0.68675 0.97730 0.31325 88.53784 0.73300

Nota: Los datos obtenidos se compararon con los obtenidos por (Shacham et al., 1998).

Fuente: Elaboración propia.

En la Figura 15 se observa el camino homotópico trazado para este caso de estudio con el vector inicial antes mencionado. Dicho camino homotópico presenta un cruce en (=1 únicamente en la sección del camino homotópico en sentido negativo en color rojo, mientras que la sección en color azul representando el camino homotópico en sentido positivo tiende a -( en τ.

Fuente: Elaboración propia.

Figura 15 Camino homotópico para el caso de estudio 6. 

Caso de estudio 7

El siguiente caso de estudio contempla un análisis de transferencia de calor y masa para una reacción catalítica en una partícula plana de un catalizador poroso (Lin et al., 2008). Suponiendo una geometría plana para la partícula, y que la transferencia de calor por conducción es despreciable comparada con la transferencia de calor por convección, se obtiene la siguiente ecuación diferencial como resultado directo del balance de materia y energía:

y''=λyeγβ1-y1+β1-y            t[0,1] (37)

donde γ=20 es una energía de activación adimensional y β=0.4 es un parámetro adimensional que describe la variación del calor. Las condiciones de frontera son y'(0)=0 y y(1)=1.

La ecuación 37 se discretiza por diferencias finitas, dando como resultado el siguiente sistema de ecuaciones algebraicas no lineales:

f1 y1,y2,,yn =-y3+4y2-y12Δx=0

fi y1,y2,,yn =yi+1-2yi+yi-1Δx2-λyi eγβ1-yi 1+β1-yi =0

para i=2,3,,n-1

fn (y1,y2,,yn )=yn-1=0 (38)

De acuerdo con la literatura, para diferentes valores de λ se obtiene diferente número de soluciones. Utilizando el método de SHH con un número de nodos de 21 y una aproximación inicial de 1.0, se encontraron los resultados mostrados en la Tabla 12.

Tabla 12 Soluciones encontradas para el caso de estudio 7 con diferentes valores de λ. 

λ Solución encontrada
(y0=Vi)
Tiempo de cómputo
0.05 V1= 0.9703436507687250 1.0 s
0.10 V1= 0.0641951468811831
V2=0.9226651927924107
V3= 0.5059615086540787
3.0 s
0.15 V1=0.01647925314865238 2.0 s

Nota: Los datos obtenidos se compararon con los obtenidos por (Lin et al., 2008).

Fuente: Elaboración propia.

En la Figura 16a se muestran las soluciones obtenidas con diferentes aproximaciones iniciales para una malla con 51 nodos con λ=0.1 y una aproximación inicial de 1.0. En la Figura 16b se muestra el camino homotópico trazado para este caso de estudio con el sentido positivo en color azul, el sentido negativo en color rojo y una línea en color verde que representa a τ=1 a fin de visualizar el número de cruces que representan las soluciones del sistema.

Fuente: Elaboración propia.

Figura 16 a) Gráfico de soluciones para el caso de estudio 7; b) Gráfico del camino homotópico para el caso de estudio 7. 

Finalmente, este caso de estudio, al igual que casos de estudio anteriores, se evaluó con los métodos BCH y Newton-Raphson multivariable con las mismas condiciones iniciales. Los resultados para λ=0.1 se muestran en la Tabla 13. Las soluciones encontradas se muestran numeradas de acuerdo con la Tabla 12, mientras que la precisión de cada solución con respecto a fx=0 se muestra en la parte inferior de la Tabla.

Tabla 13 Análisis de caso de estudio 7 con métodos de convergencia local y global. 

Aproximación inicial y0i=1nodos Método Solución encontrada (Vi)
[11 nodos; 15 nodos; 21 nodos]
Tiempo de cómputo
[11 nodos; 15 nodos; 21 nodos]
-1.0 BCH [V1 ; V1 ; V1 ] [2.3 s; 4.0 s; 8.3 s]
Newton-Raphson [V1 ; V1 ; V1 ] [0.0 s; 0.0 s; 0.0 s]
SHH [V1 , V2 , V3 ; V1 , V2 , V3 ; V1 , V2, V3 ] [3.0 s; 3.0 s; 5.0 s]
0.0 BCH [V1 , V1 , V1 ] [2.3 s; 4.2 s; 8.4 s]
Newton-Raphson [V1 ; V1 ; V1 ] [0.0 s; 0.0 s; 0.0 s]
SHH [V1 , V2 , V3 ; V1 , V2 , V3 ; V1 , V2, V3 ] [5.0 s; 4.0 s; 7.0 s]
2.0 BCH [V2 , V2 , V2 ] [2.3 s; 4.2 s; 8.3 s]
Newton-Raphson [V2 , V2 , V2 ] [0.0 s; 0.0 s; 0.0 s]
SHH [V1 , V2 , V3 ; V1 , V2 , V3 ; V1 , V2 , V3 ] [5.0 s; 4.0 s; 4.0 s]
3.0 BCH [V2 , V2 , V2 ] [2.3 s; 4.2 s; 8.3 s]
Newton-Raphson [V2 , V2 , V2 ] [0.0 s; 0.0 s; 0.0 s]
SHH [V1 , V2 , V3 ; V1 , V2 , V3 ; V1 , V2 , V3 ] [3.0 s; 5.0 s; 5.0 s]
i=1nabsfix1,x2,,xn0.1 y i=1nabsfix1,x2,,xn1×10-7

Fuente: Elaboración propia.

Simulador de Seguimiento Homotópico Circular (SHC)

Se desarrolló un simulador en MATLAB® de seguimiento homotópico circular; es decir, a diferencia del seguimiento hiperesférico que es n-dimensional, este simulador solo puede resolver una única ecuación algebraica no lineal. Esto con el objetivo de comprender los principios del seguimiento homotópico con corrección a base de una geometría circular.

Para este método, la ecuación algebraica se deforma con la ecuación de homotopía de Newton y se escribe de acuerdo con la ecuación 2. Posteriormente, se resuelve en conjunto con la ecuación del círculo:

Cx,τ=x-x02-τ-τ02-R2=0 (39)

De la ecuación 2 se despeja τ y se sustituye en la ecuación anterior, obteniéndose:

fx0-fxfx0-τ02+x-x02+R2=0 (40)

La metodología de seguimiento es la siguiente:

  • 1) A partir del punto τ0=0 y x=x0, que será el centro del círculo, se genera un círculo con una R fijada previamente.

  • 2) Se genera un vector tangente que pase por (τ0, x0) y se prolonga hasta que intersecte el círculo en el punto (τa, xa).

  • 3) A partir de ese punto, se evalúa el valor x1 usando el método iterativo Newton-Raphson descrito previamente. Esto es, propiamente, un proceso de corrección.

  • 4) Calcular τ1 a través de la ecuación 40.

  • 5) El punto (τ1, x1), el cual es la intersección entre H(x, τ) y C(x, τ), será a su vez el centro del siguiente círculo. Este proceso se muestra en la Figura 17a y se repetirá sucesivamente hasta trazar el camino homotópico, tal como se muestra en la representación gráfica de la Figura 17b con un radio constante.

  • 6) Cuando el camino cruce la región donde τ=1, se debe refinar la solución encontrada con el uso del método de Newton-Raphson (Oliveros-Muñoz & Jiménez-Islas, 2013) .

Fuente: Elaboración propia

Figura 17 a) Representación gráfica de método homotópico circular; b) Representación gráfica del seguimiento homotópico circular con radio constante. 

Para la primera prueba del simulador se utilizó la función polinómica f(x)=x3-1.473x2-5.738x+6.763. En la Figura 18 se muestran los gráficos que arrojó el simulador. En dichos gráficos se puede observar la búsqueda del camino homotópico para la función estipulada, en sentido positivo y negativo, así como el camino homotópico completo, excluyendo el trazado de círculos para este último. En esta prueba se utilizó un radio constante de 0.5 y una aproximación inicial de 0. El tiempo de cálculo de la computadora fue de 9.971 s, y se encontraron todas las raíces de la función. Los resultados se muestran en la Tabla 14.

Fuente: Elaboración propia.

Figura 18 Trazo del camino homotópico realizado por el simulador para la función a resolver. 

Tabla 14 Raíces de la función encontradas por el simulador. 

x 1.09997
-2.30007
2.67306

Fuente: Elaboración propia.

Discusión

Observaciones

Caso de estudio 1: El SHH iguala al método de homotopía con límite en su parámetro homotópico, encontrando las nueve soluciones reportadas en la literatura. En cuanto al análisis con diferentes métodos de convergencia global y local, el SHH demostró una notable ventaja en cuanto al número de vectores solución encontrados con la misma aproximación inicial utilizada en todos los métodos y con un tiempo de cómputo similar.

Caso de estudio 2: El SHH supera el método utilizado por la literatura para resolver el caso de estudio presentado, dado que es capaz de encontrar el doble de vectores solución con un solo vector inicial.

Caso de estudio 3: En el trazo del camino homotópico p presenta una disminución considerable en su magnitud, por lo que es necesario utilizar un factor de multiplicación grande para poder tener una mejor visualización al momento de graficar el camino homotópico. En este caso en particular, los vectores iniciales no pueden tener el mismo valor para cada una de las variables, esto se debe al problema de aplicación del que surgen las ecuaciones algebraicas no lineales, se observa que las variables cambian considerablemente de magnitud.

Caso de estudio 4: El SHH es capaz de encontrar todas las soluciones encontradas en la literatura con un solo vector inicial en un corto periodo de tiempo, tomando en cuenta que se trata de un sistema de nueve ecuaciones algebraicas no lineales. En el análisis con diferentes métodos de convergencia, el método de SHH demostró que es capaz de encontrar múltiples soluciones, mientras que con el método BCH no pudo obtenerse ningún vector solución con ninguna de las aproximaciones iniciales propuestas.

Caso de estudio 5: Los métodos utilizados para resolver el problema en la literatura son capaces de encontrar una solución por cada vector inicial utilizado, sin embargo, el SHH fue capaz de encontrar cuatro soluciones en una sola corrida incluyendo la única solución que es satisfactoria para el problema que representa el sistema de ecuaciones.

Caso de estudio 6: El SHH presenta dificultades para encontrar la solución con los vectores iniciales propuestos en la literatura, sin embargo, fue capaz de encontrar la solución en la mayoría de los vectores aleatorios utilizados para evaluarlo.

Caso de estudio 7: Problemas que incluyen ecuaciones diferenciales ordinarias también puede ser resueltos con el método de SHH discretizando las ecuaciones por métodos como el de diferencias finitas. Para este caso en particular se utilizó un valor de λ que aseguraba la aparición de al menos tres soluciones y se buscó un número de nodos óptimo para obtener la mejor aproximación de la solución. Por otro lado, con el propósito de una mejor visualización en el gráfico, se utilizó un número mucho mayor de nodos. En cuanto al análisis por diferentes métodos de convergencia, el método de SHH ocupó un poco de más tiempo de cómputo, pero se obtuvo las tres soluciones en cada ocasión, al contrario que los métodos BCH y Newton-Raphson.

Conclusiones

En comparación con otros métodos de resolución de ecuaciones algebraicas no lineales, el método de SHH es capaz de encontrar múltiples soluciones partiendo de un vector inicial arbitrario, lo cual es una clara ventaja contra métodos de convergencia local como el método de Newton-Raphson y el método basado en continuación homotópica (BCH) (Martín, 2013). En algunos casos de estudio, el algoritmo ha encontrado más vectores solución que los reportados en la literatura. Finalmente, gracias al parámetro p, es posible graficar el camino homotópico de una manera simple para facilitar la visualización del número de soluciones encontradas por el método de SHH. Los resultados obtenidos a la fecha son promisorios, pero hace falta mayor investigación para satisfacer las dos premisas básicas de la continuación o seguimiento homotópico que son: 1) El método funciona desde cualquier aproximación inicial y 2) El camino homotópico contiene todas las soluciones (al menos las soluciones del plano real) del sistema algebraico analizado, lo que sigue siendo un desafío para continuar las futuras investigaciones sobre el tema.

Conflictos de interés

Los autores manifiestan que no existen conflictos de interés.

Agradecimientos

Los autores agradecen al Tecnológico Nacional de México el soporte financiero de la investigación a través del proyecto 7607-20. P.

Referencias

Asadi, J., & Jalali, F. (2018). Optimization of dimethyl ether production process based on sustainability criteria using a homotopy continuation method. Computers & Chemical Engineering, 115, 161-178. doi: https://doi.org/10.1016/j.compchemeng.2018.03.014 [ Links ]

Diaz-Arango, G., Vázquez-Leal, H., Hernandez-Martinez, L., Sanz, M. T., & Sandoval-Hernandez, M. (2018). Homotopy path planning for terrestrial robots using spherical algorithm. IEEE Transactions on Automation Science and Engineering, 15(2), 567-585. doi: https://doi.org/10.1109/TASE.2016.2638208 [ Links ]

Gallardo-Alvarado, J. (2019). An application of the Newton-Homotopy continuation method for solving the forward kinematic problem of the 3-RRS parallel manipulator. Mathematical Problems in Engineering, 2019. doi: https://doi.org/10.1155/2019/3123808 [ Links ]

Gupta, Y. P. (1995). Bracketing method for online solution of low-dimensional nonlinear algebraic equations. Industrial & Engineering Chemistry Research, 34(2), 536-544. doi: https://doi.org/10.1021/ie00041a014 [ Links ]

Jiménez-Islas, H., Martínez-González, G. M., Navarrete-Bolaños, J. L., Botello-Álvarez, J. E., & Oliveros-Muñoz, J. M. (2013). Nonlinear homotopic continuation methods: A chemical engineering perspective review. Industrial and Engineering Chemistry Research, 52(42), 14729-14742. doi: https://doi.org/10.1021/ie402418e [ Links ]

Kuno, M., & Seader, J. D. (1988). Computing all real solutions to systems of nonlinear equations with a global fixed-point homotopy. Industrial and Engineering Chemistry Research, 27(7), 1320-1329. doi: https://doi.org/10.1021/ie00079a037 [ Links ]

Lin, Y., Enszer, J. A., & Stadtherr, M. A. (2008). Enclosing all solutions of two-point boundary value problems for ODEs. Computers & Chemical Engineering , 32(8), 1714-1725. doi: https://doi.org/10.1016/j.compchemeng.2007.08.013 [ Links ]

Malinen, I., & Tanskanen, J. (2010). Homotopy parameter bounding in increasing the robustness of homotopy continuation methods in multiplicity studies. Computers and Chemical Engineering, 34(11), 1761-1774. doi: https://doi.org/10.1016/j.compchemeng.2010.03.013 [ Links ]

Martín, C. (2013). Homotopía y continuación numérica en sistemas no lineales. Matemática, 11(1), 13-23. http://www.revistas.espol.edu.ec/index.php/matematica/article/view/506Links ]

Oliveros-Muñoz, J. M., & Jiménez-Islas, H. (2013). Hyperspherical path tracking methodology as correction step in homotopic continuation methods. Chemical Engineering Science, 97, 413-429. doi: https://doi.org/10.1016/j.ces.2013.03.053 [ Links ]

Rahimian, S. K., Jalali, F., Seader, J. D., & White, R. E. (2011a). A new homotopy for seeking all real roots of a nonlinear equation. Computers and Chemical Engineering, 35(3), 403-411. doi: https://doi.org/10.1016/j.compchemeng.2010.04.007 [ Links ]

Rahimian, S. K., Jalali, F., Seader, J. D., & White, R. E. (2011b). A Robust Homotopy Continuation Method for Seeking All Real Roots of Unconstrained Systems of Nonlinear Algebraic and Transcendental Equations. Industrial & Engineering Chemistry Research, 50(15), 8892-8900. doi: https://doi.org/10.1021/ie101966b [ Links ]

Seader, J. D., Kuno, M., Lin, W. J., Johnson, S. A., Unsworth, K., & Wiskin, J. W. (1990). Mapped continuation methods for computing all solutions to general systems of nonlinear equations. Computers and Chemical Engineering, 14(1), 71-85. doi: https://doi.org/10.1016/0098-1354(90)87006-B [ Links ]

Shacham, M., Brauner, N., & Pozin, M. (1998). Comparing software for interactive solution of systems of nonlinear algebraic equations. Computers and Chemical Engineering, 22(1-2), 323-331. doi: https://doi.org/10.1016/s0098-1354(97)88454-2 [ Links ]

Tang, L., Ruan, J., & Chen, R. (2019). A homotopy continuation method for nonlinear electric field analysis. IEEE Transactions on Dielectrics and Electrical Insulation, 26(4), 1125-1133. doi: https://doi.org/10.1109/TDEI.2019.007882 [ Links ]

Torres-Muñoz, D., Hernandez-Martinez, L., & Vázquez-Leal, H. (2016). Spherical continuation algorithm with spheres of variable radius to trace homotopy curves. International Journal of Applied and Computational Mathematics, 2, 421-433. doi: https://doi.org/10.1007/s40819-015-0067-1 [ Links ]

Cómo citar: Quemada-Villagómez, M. L., López-González, M. L., Oliveros-Muñoz, J. M. & Jiménez-Islas, H. (2022). Enseñanza del método de continuación homotópica con seguimiento hiperesférico para estudiantes de ingeniería. Acta Universitaria 32, e3358. doi. http://doi.org/10.15174.au.2022.3358

Apéndices Apéndice A

Fuente: Elaboración propia.

Figura A1 Superficies de H1(x1, x2) y H2(x1, x2) para la aproximación inicial x0=1,1t

Fuente: Elaboración propia.

Figura A2 Vista 2D (x1x2) de la deformación que sufren las superficies de H1(x1, x2) y H2(x1, x2) y sus intersecciones con H= 0, para la aproximación inicial x0=1,1t

Apéndice B

Pseudo Código de seguimiento homotópico hiperesférico

Paso 1. Definir valores iniciales de las variables x1ini,x2ini,,xnini, radio (r), valor inicial del parámetro homotópico (τ0), valor inicial del parámetro de avance (p0), número máximo de generación de esferas (itermáximasseguimiento), número máximo de iteraciones en la corrección de Newton ( itermáximascorreción), tolerancia para identificar cruces del camino homotópico (toleranciacruce).

Paso 2. Definir el número de ecuaciones (n) y el sistema de ecuaciones algebraicas no lineales.

f1x1,x2,, xnf2x1,x2,, xnfnx1,x2,, xn=000

Paso 3. Iniciar k=1, x10,x20,,xn0=x1ini,x2ini,,xnini

Paso 4. Calcular la deformación homotópica con homotopía de Newton.

H1x1,x2,, xn,τ H2x1,x2,, xn,τHnx1,x2,, xn,τ=f1x1,x2,, xn-1+τf1x1ini,x2ini,, xninif2x1,x2,, xn-1+τf2x1ini,x2ini,, xninifnx1,x2,, xn-1+τfnx1ini,x2ini,, xnini

Paso 5. Parametrización del sistema homotópico, xixip y ττp.

Paso 6. Definir un sistema matricial, HMn×n+1R, con la derivada parcial de Hixip;τp con respecto al parámetro de avance p.

H=H1x1x1pH2x1x1pH1x2x2pH1xnxnpH2x2x2pH2xnxnpH1ττpH2ττpHnx1x1pHnx2x2pHnxnxnpHnττp

Paso 7. El vector predictor, x1,x2,,xn,τ,pkpredictor, se calcula al resolver el sistema H con la regla de Cramer.

Paso 7.1. Se generaliza la pendiente para el método de Euler como,

xi=xip=-1idetH-i, donde i=1, 2, ,n

τ=τp=-1n+1detH-n+1

Paso 7.2. Calcular el avance del parámetro de avance p,

p±r2i=1nxi2+τ2

Paso 7.3. Entonces el vector predictor es,

xipredictor=xi0±pxip, para i=1,2,,n

τpredictor=τ0±pτp

ppredictor=p0±p

Paso 8. Realizar la corrección del vector predictor (j=1, 2, ,itermáximascorreción).

Paso 8.1. Definir el sistema de ecuaciones para la corrección incluyendo la ecuación de la hiperesfera.

f1x1,x2,, xn-1+τf1x1ini,x2ini,, xninif2x1,x2,, xn-1+τf2x1ini,x2ini,, xninifnx1,x2,, xn-1+τfnx1ini,x2ini,, xninii=1nx1-xi02+τ-τ02-r2  =0000

Paso 8.2. Generar el vector corrector x1,x2,,xn,τ,pkcorrector resolviendo con Newton-Raphson el sistema definido en el Paso 8.1.

Paso 9. Guardar camino homotópico generado: x1,x2, , xn,τ,pkcorrector.

Paso 10. If absτk-1toleranciacruce then

Paso 10.1. Refinar el cruce del camino homotópico con el método de Newton-Raphson.

Paso 10.2. Guardar la solución del sistema, x1,x2, , xnk.

Else

Continuar en el paso 11

End if

Paso 11. k=k+1, x10,x20, , xn0,τ0,p0=x1,x2, , xn,τ,pkcorrector  

Paso 12. If k=itermáximasseguimiento then

Paso 12.1. Terminar el seguimiento homotópico hiperesférico.

Else

Continuar en el paso 4

End if

Apéndice C

Comparativa del SHH y el BCH

El algoritmo basado en continuación homotópica denominado BCH es descrito por Martín (2013), partiendo de que la curva X=ψτ es una solución del sistema homotópico HX,τ=0 para todo 0τ1, si y solo si X=ψτ es la solución del problema de valor inicial descrito en la ecuación (C1) con aproximación inicial de HX0,0=0 y τ[0, 1].

X'τ=-HXX, τ-1HτX,τ (C1)

El problema de valor inicial considera el sistema homotópico HX,τ=τFX+1-τGX donde GX=FX-FX0. Por lo tanto, el sistema es HX,τ=τFX+1-τFX-FX0=0 y es simplificado como FX=1-τFX0=0 con su derivada respecto al parámetro homotópico (τ), F'XX'τ=-FX0. Lo anterior permite generalizar el problema de valor inicial para n ecuaciones en la ecuación (C2):

X'τ=-F'X-1FX0 (C2)

donde F'X es la matriz jacobiana del sistema de n ecuaciones, FX0 es la función evaluada con la aproximación inicial (X0) y X'τ es la derivada de cada variable independiente con respecto al parámetro homotópico τ. Entonces, el problema descrito en la ecuación (C2) es resuelto con el método de integración numérica Runge-Kutta de 4° orden, partiendo de la aproximación inicial X0 desde τ=0 hasta τ=1, donde se encuentra la solución del sistema X. La precisión en la solución encontrada dependerá del número de iteraciones (n=1/h), las cuales son definidas por el tamaño de paso (h) utilizado durante la integración numérica (Martín, 2013).

El algoritmo BCH de manera similar al seguimiento homotópico hiperesférico/circular (SHH/SHC) determina la solución al sistema cuando el parámetro homotópico τ es igual a 1. Sin embargo, el intervalo con respecto al parámetro homotópico, [0τ1], limita de forma sustancial la búsqueda de soluciones y el algoritmo BCH dependerá de la aproximación inicial (X0) para encontrar las soluciones del sistema (Xi); es decir, es necesario modificar la aproximación inicial en cada corrida del algoritmo para obtener soluciones diferentes, funcionamiento parecido al método de convergencia local de Newton-Raphson, donde una aproximación inicial siempre converge a la misma solución del sistema. Aun así, se debe considerar la principal diferencia del algoritmo BHC con el método de convergencia local Newton-Raphson; a diferencia de este, el algoritmo BCH no necesitará una aproximación inicial lo suficientemente cercana a la solución para alcanzar convergencia.

A continuación, se muestran los resultados de la comparación del algoritmo SHH con el algoritmo BCH propuesto por Martín(2013). Para fines de esta comparación, se replicó el algoritmo basado en continuación homotópica (BCH) en FORTRAN(, y de manera adicional se incluye la evaluación con el método de Newton-Raphson (NR). Los casos de estudio que se presentarán a continuación fueron evaluados en las mismas condiciones en un equipo computacional con procesador Intel® Xeon® E5-2620 V3 @2.4GHz con 64GB de RAM.

Caso de estudio A

En el primer caso se tomó la ecuación (C3) junto a las aproximaciones iniciales reportadas por Martín (2013), en las cuales se asegura convergencia a una solución de f(x).

fx=x3-exsenx (C3)

En la Tabla C1 se muestran los resultados de tres métodos numéricos empleados para encontrar las soluciones de la ecuación (C3). En este caso, se aprecia que el algoritmo SHC encuentra un mayor número de soluciones a partir de una sola aproximación inicial y se observa la capacidad del algoritmo de encontrar soluciones adicionales a las reportadas en la literatura, a diferencia de los algoritmos BCH y NR, donde solo se observa una solución por cada aproximación inicial.

Tabla C1 Comparación del SHC, BCH y NR para el caso de estudio A.  

Método 𝒙𝟎 Soluciones encontradas(𝒙𝒊, 𝒇(𝒙𝒊) = 𝒙𝟑 −𝒆𝒙𝒔𝒆𝒏(𝒙)) Tiempo de
cómputo (s)
NR 1.0 x3 = −0.683126494476553 , f(x3) = 0.000000000000000 0.10 s
1.5 x2 = 1.811257124852834, f(x2) = 0.000000000000001 0.10 s
-1.0 x3 = −0.683126494476553, f(x3) = 0.000000000000001 0.10 s
10.0 x4 = 9.353816905054659, f(x4) = 0.000000000003865 0.10 s
BCH 1.0 x1 = 0.000000000002075, f(x1) = 0.000000000002075 0.10 s
1.5 x2 = 1.811257124852842, f(x2) = 0.000000000000038 0.10 s
-1.0 x3 = −0.683126494476544, f(x3) = 0.000000000000012 0.10 s
10.0 x4 = 9.353816905054657, f(x4) = 0.000000000015348 0.10 s
SHC 1.0 x1 = 0.000000000000000, f(x1) = 0.000000000000000
x2 = 1.811257124852834, f(x2) = 0.000000000000001
x3 = −0.683126494476553, f(x3) = 0.000000000000001
x4 = 9.353816905054659, f(x4) = 0.000000000003865
x5 = 6.669246061151054, f(x5) = 0.000000000000001
1.45 s
1.5 x1 = 0.000000000000000, f(x1) = 0.000000000000000
x2 = 1.811257124852834, f(x2) = 0.000000000000001
x3 = −0.683126494476553, f(x3) = 0.000000000000001
x4 = 9.353816905054659, f(x4) = 0.000000000003865
x5 = 6.669246061151054, f(x5) = 0.000000000000001
1.46 s
-1.0 x1 = 0.000000000000000, f(x1) = 0.000000000000000
x2 = 1.811257124852834, f(x2) = 0.000000000000001
x3 = −0.683126494476553, f(x3) = 0.0000000000000001
x5 = 6.669246061151054, f(x5) = 0.000000000000001
1.62 s
10.0 x1 = 0.000000000000000, f(x1) = 0.000000000000000
x2 = 1.811257124852834, f(x2) = −0.000000000000001
x3 = −0.683126494476553, f(x3) = 0.000000000000001
x4 = 9.353816905054659, f(x4) = 0.000000000003865
x5 = 6.669246061151054, f(x5) = 0.000000000000001
x6 = 12.573254781521957, f(x6) = 0.000000000242835
x7 = 15.707378904286850, f(x7) = 0.000000004949015
1.51 s

Nota. Los datos obtenidos se compararon con los obtenidos por Martín (2013).

Fuente: Elaboración propia.

En la Figura C1 se observa la interpretación gráfica del camino homotópico generado por los algoritmos BCH y SHC. En primera instancia, la Figura C1-a muestra el camino homotópico generado con la aproximación inicial de 1.0 con el algoritmo SHC y se indica la Sección del camino donde se observan los cruces del camino homotópico con τ=1. La Figura C1-b muestra la Sección 1 del camino homotópico ampliada, donde se presentan al menos cinco cruces que representan soluciones diferentes de f(x). En tanto la Figura C1-c muestra el funcionamiento del algoritmo BCH, donde cada aproximación inicial converge solo a una solucion del sistema.

Figura C1 Representación gráfica de los algoritmos BCH y SHC. 

Caso de estudio B

A continuación, se presenta un sistema de dos ecuaciones no lineales.

f1x,y=5x2-y2=0

f2x,y=4y-senx-cosy=0 (C4)

En la Tabla C2 se muestran los resultados de los algoritmos SHH, BCH y NR con el sistema descrito en la ecuación (C4), empleando las mismas aproximaciones iniciales sugeridas por Martín (2013). El valor de las funciones con la solución se reporta como Exi,yi=j=1nabsfjxi,yi, donde j es el número de ecuaciones del sistema y los subíndices i indican el número de solución encontrada. En este caso, se aprecia que el SHH localiza las dos soluciones reportadas con una sola aproximación inicial en todas las pruebas realizadas.

Tabla C2 Comparación de los algoritmos SHH, BCH y NR para el caso de estudio B. 

Método (𝐱𝟎,𝐲𝟎) (𝐱𝐢, 𝐲𝐢, 𝐄(𝐱𝐢, 𝐲𝐢)) Tiempo de
cómputo (s)
NR (1, 1) x1 = 0.121241911480502
y1 = 0.271105155792415
E(x1, y1) = 0.000000000000001
0.1 s
(-3, -3) x2 = −0.098163444944170
y2 = 0.219500135800723
E(x2, y2) = 0.000000000000001
0.1 s
(2, 0) x1 = 0.121241911480502
y1 = 0.271105155792415
E(x1, y1) = 0.000000000000001
0.1 s
(-2, -1) x2 = −0.098163444944170
y2 = 0.219500135800723
E(x2, y2) = 0.000000000000001
0.1 s
BCH (1, 1) x1 = 0.121242254781049
y1 = 0.271105284876040
E(x1, y1) = 0.00000056356444
1.1 s
(-3, -3) x2 = −0.098166046641315
y2 = 0.219499328329200
E(x2, y2) = 0.0000003724963196
1.2 s
(2, 0) x1 = 0.121243429542644
y1 = 0.271105478716969
E(x1, y1) = 0.0000001794187867
1.2 s
(-2, -1) x2 = −0.098164854603032
y2 = 0.219499750906795
E(x2, y2) = 0.00000031773258715
1.1 s
SHH (1, 1) x1 = 0.1212419114796510
y1 = 0.27110515579220072
E(x1, y1) = 0.000000057209067483
x2 = −0.09816344494457177
y2 = 0.2195001358003615
E(x2, y2) = 0.000000002643687533
2 s
(-3, -3) x2 = −0.09816344526613148
y2 = 0.2195001357245536
E(x2, y2) = 0.000000000004846054
x1 = 0.1212419120849824
y1 = 0.271105155935472
E(x1, y1) = 0.000000000299723871
4 s
(2, 0) x1 = 0.1212419116733243
y1 = 0.2711051558379828
E(x1, y1) = 0.000000000001206840
x2 = −0.09816344498488395
y2 = 0.219500135791123
E(x2, y2) = 0.000000045184565936
7 s
(-2, -1) x2 = −0.0987634451807615
y2 = 0.2195001357447709
E(x2, y2) = 0.000000070768595067
x1 = 0.1212419114798090
y1 = 0.2711051557922
E(x1, y1) = 0.000000000001167844
5 s

Nota. Los datos obtenidos se compararon con los obtenidos por (Martín, 2013).

Fuente: Elaboración propia.

Caso de estudio C

Finalmente, se replica un caso de tres ecuaciones algebraicas no lineales.

f1x,y,z=x-22+y-12+x y-3=0

f2x,y,z=x ex+y+x y-3=0

f3x,y,z=senx+y+x-y-z+1=0 (C5)

En la Tabla C3 se observan las dos soluciones encontradas con los algoritmos BCH, SHH y NR para el sistema de ecuaciones (ecuación C5). El valor de las funciones con la solución se obtiene con Exi,yi,zi=j=1nabsfjxi,yi,zi, donde j es el número de ecuaciones del sistema y los subíndices i indican el número de solución encontrada. Al igual que en el caso anterior, a partir de una sola aproximación inicial el algoritmo SHH encuentra las soluciones reportadas en la literatura.

Tabla C3 Comparación de los algoritmos SHH, BCH y NR para el caso de estudio C. 

Método (𝐱𝟎, 𝐲𝟎, 𝐳𝟎) [𝐱𝐢, 𝐲𝐢, 𝐳𝐢, 𝐄(𝐱𝐢, 𝐲𝐢, 𝐳𝐢)]𝐓𝐓 Tiempo de cómputo (s)
NR (1, 1, 1) 𝑥1 = 0.388522261069885
𝑦1 = 1.034550195014107
𝑧1 = 1.341346675089112
𝐸(𝑥1, 𝑦1, 𝑧1) = 0.0000000000000001
0.1 s
(4, -2, 7) 𝑥2 = 3.039247388237472
𝑦2 = −1.610482827797872
𝑧2 = 6.013364389595542
𝐸(𝑥2, 𝑦2, 𝑧2) = 0.0000000000000001
0.1 s
BCH (1, 1, 1) 𝑥1 = 0.3885222610700
𝑦1 = 1.0345501950139
𝑧1 = 1.341346675089)
𝐸(𝑥1, 𝑦1, 𝑧1) = 0.00000000000001
0.1 s
(4, -2, 7) 𝑥2 = 3.039247388237
𝑦2 = −1.610482827797
𝑧2 = 6.013364389595
𝐸(𝑥2, 𝑦2, 𝑧2) = 0.00000000000001
0.1 s
SHH (1, 1, 1) 𝑥1 = 0.3885222605894432
𝑦1 = 1.034550202288378
z1 = 1.341346668659938
𝐸(𝑥1, 𝑦1, 𝑧1 0.0000000000000001
𝑥2 = 3.039247401546912
𝑦2 = −1.610482840131396
𝑧2 = 6.013364415908800
𝐸(𝑥2, 𝑦2, 𝑧2) = 0.0000000000000001
4 s
(4, -2, 7) 𝑥1 = 0.3885222605894432
𝑦1 = 1.034550202288378
𝑧1 = 1.341346668659938
𝐸(𝑥1, 𝑦1, 𝑧1) = 0.0000000000000001
𝑥2 = 3.039247401546912
𝑦2 = −1.610482840131396
𝑧2 = 6.013364415908800
𝐸(𝑥2, 𝑦2, 𝑧2) = 0.0000000000000001
8 s

Nota. Los datos obtenidos se compararon con los obtenidos por Martín (2013).

Fuente: Elaboración propia.

De acuerdo con las observaciones de los tres casos de estudio previamente presentados, se pueden concluir varios puntos importantes:

  • 1. El SHH/SHC son algoritmos de continuación homotópica que se destacan por su capacidad para encontrar múltiples soluciones con una sola aproximación inicial.

  • 2. El método de SHH facilita la representación gráfica del camino homotópico con el cual es posible visualizar el número de soluciones que se pueden encontrar para cualquier caso de estudio.

  • 3. Una vez que se conoce el método de SHH, es posible adecuarlo para poder realizar un mejor seguimiento del camino homotópico y, de esta manera, poder encontrar más soluciones, o bien, optimizar el tiempo de cómputo.

Recibido: 26 de Octubre de 2021; Aprobado: 23 de Febrero de 2022; Publicado: 06 de Abril de 2022

*Autor de correspondencia hugo.jimenez@itcelaya.edu.mx

Creative Commons License Este es un artículo publicado en acceso abierto bajo una licencia Creative Commons