Arnoldi
methods for numerical continuation of stationary points
Zenaida Natividad Castillo Marrer=
o.[1], =
Gustavo
Adolfo Colmenares Pacheco.[2] &=
amp;
Víctor Oswaldo Cevallos Vique.[3]
Recibido:
24-04-2020 / Revisado: 15-05-2020 / Aceptado:
25-06-2020 / Publicado: 03-07-2020
Abstract. =
DOI: https://doi.org/10.33262/cienciadigital.v=
4i3.1385
This work presents an analysis of a continuation algorithm with an
embedded Arnoldi eigen solver for the numerical
computation of solutions of nonlinear systems G (x, α) =3D 0, where x =
is a
vector in Rn and α is a parameter which takes values in a g=
iven
interval. This technique allows us to detect and predict particular solutio=
ns
when computing the eigenvalues of the associate Jacobian matrix, and
simultaneously to get the solution of linear systems in each iteration. This
method could be applied to problems of electrical engineering, chemical
reactions or coating process. The main idea is to embed an Arnoldi
eigenvalue solver in a continuation algorithm to compute solutions of a
nonlinear system in order to get additional information of these solutions,
such as the stability or the periodicity. An example of the usability of th=
is
technique is presented, with preliminary results on models used in the
bibliography of this topic. The results are encouraged, and show the
reliability of this approach in the accurate detection of critical points.<=
/span>
Keywords: Path following, Nonlinear systems, Blo=
ck Arnoldi method, Stability Analysis, Critical points.<=
/span>
Resumen
En este trabajo se describe un método de continuación numérica, el cual uti=
liza
descomposiciones de Arnoldi para el cálculo de
soluciones de un sistema de la forma
Palabras claves: Continu=
ación
numérica, Sistemas no lineales, Métodos de Arnoldi,
Análisis de estabilidad, Puntos críticos.
Un número significativo de problemas que surgen en las áreas c=
ientíficas,
tales como dinámica poblacional, teoría de combustión, dinámica y mecánica =
de
fluidos, pueden ser modelados por sistemas de ecuaciones diferenciales
parciales (PDE’s), y luego discretizados en esp=
acio, dando
lugar a lo que se conoce como sistemas dinámicos no lineales, los cuales
dependen de uno o más parámetros de interés.
Una manera de analizar estos fenómenos físicos, una vez que son
modelados por PDE’s es analizando las caracterí=
sticas
de las soluciones de un sistema de ecuaciones no lineales de la forma
Algunos cambios en estos sistemas son detectados cuando se pierde
estabilidad en la continuación de las soluciones, o también cuando las
soluciones comienzan o dejan de ser periódicas. En este caso, se requiere c=
alcular
con suficiente precisión estos puntos o soluciones especiales, ya que esta
información puede ser usada para modificar los parámetros de operación con =
el
fin de mejorar la productividad o la calidad del proceso que se esté modela=
ndo.
Los métodos de continuación más usados en la actualidad son
aquellos basados en un esquema a dos pasos, del tipo predictor-corrector. E=
n el
primer paso (predictor) el método genera un nuevo valor para el parámetro <=
/span>
En este trabajo se muestra como diseñar e implementar un métod=
o de
continuación numérica con un esquema predictor-corrector, que sea competiti=
vo y
que además de generar el diagrama de bifurcación, nos permita visualizar los
puntos que representan soluciones estables o inestables o soluciones que
señalen cambios de comportamiento del proceso modelado. Para ello, combinar=
emos
el método de continuación con un método de proyección en espacios de Krylov del tipo Arnoldi, =
a fin de
resolver el sistema lineal que aparece en cada iteración del paso corrector=
y
al mismo tiempo obtener información sobre los autovalores y autovectores
de la matriz Jacobiana asociada.
Se plantea en este trabajo una propuesta, con experimentación =
preliminar
sobre el uso de la descomposición de Arnoldi de=
ntro
de un esquema de continuación, a fin de analizar la estabilidad de las
soluciones que se calculan en cada iteración que resuelve el problema
1) Análisis de los métodos de continuación y la forma general de sus implementaciones. 2) A partir de un algo=
ritmo
general para un método de continuación se estudiará la inserción del módulo
para el cálculo de algunos autovalores de interés, de tal forma que no se d=
egrade
la eficiencia y que se contribuya no solo al cálculo de la próxima solución
sino también a detectar pérdida de estabilidad en la solución, si la hubier=
a. 3) Diseño e implementación de la soluci=
ón
propuesta. 4) Pruebas experimentales en datos
sintéticos. 5) Organización y divulgación de result=
ados.
Muchas aplicaciones industriales son modeladas matemáticamente=
por
sistemas de ecuaciones en derivadas parciales o PDE's<=
/span>,
y luego discretizadas en espacio para construir soluciones temporales. Los
modelos más usados están basados en esquemas de diferencias finitas o eleme=
ntos
finitos, tal como explica Saad (2011). Cuando en el modelo existe un paráme=
tro
de importancia, cuya variación en un rango se desea analizar, esta discreti=
zación
da origen a un sistema de ecuaciones no lineales parametrizado de la forma:=
=
donde la función
Una metodología muy utilizada para encontrar soluciones que
satisfagan
El conjunto de soluciones o ramas obtenidas en un método de
continuación, pueden ser graficadas en 2D o 3D para establecer la relación
entre la variación del parámetro
El diagrama de bifurcación no es el gráfico de la solución de =
1) &n=
bsp;
Puntos donde dos o más ramas se intersecan BP (Branching
Points).
2) &n=
bsp;
Puntos donde la rama cambia de dirección TP (Turning
Points).
3)
La detenc=
ión
y cálculo eficiente de estos puntos es de vital importancia en aplicaciones
industriales, una vez que permite cambiar o ajustar el esquema de operación=
. A continuación,
detallamos el esquema de todo algoritmo de continuación.
Dado un p=
unto
estacionario inicial
La idea detrás del algoritmo predictor-corrector es seguir una
curva de soluciones, resolviendo sistemas no lineales parametrizados de la
forma dada en
Para entender este procedimiento, supongamos que ya tenemos una
solución inicial
Figura 1. Algoritmo
=
Fuente:=
Algoritmo
desarrollado en la investigación.
En el paso 1(a) de la i-ésima iter=
ación del
Algoritmo
En este caso, el
predictor se calcula tomando un punto sobre la recta tangente a
Donde
En ambos tipos de
continuación (orden cero y orden uno) el corrector puede usar un método
iterativo del tipo Newton o quasi-Newton. La Fi=
gura 2
muestra la diferencia entre el predictor de la continuación de orden cero y
aquel que determina la continuación de orden uno. En la misma se observa qu=
e el
predictor en la continuación de orden uno está más cerca de la nueva soluci=
ón
que el predictor de orden natural o cero, lo que sugiere que el método
iterativo convergerá más rápidamente a la solución
Figura 2. Predictor de orden cero y
Predictor de orden 1
=
Fuente:=
Elaboración
propia.
En cada paso del método iterativo que resuelve el funcional no
lineal se resuelve un sistema de ecuaciones lineales de la forma
Por otra parte, los problemas de interés en aplicaciones reales
conllevan a sistemas lineales con matrices de gran magnitud, generalmente
dispersas o ralas, lo que determina que métodos clásicos para matrices dens=
as
no son de utilidad. En las últimas
décadas se han desarrollados métodos iterativos basados en espacios de Krylov que además de contribuir a la solución del sis=
tema
Los métodos basados en Arnoldi se
fundamentan en la descomposición Hessenberg de =
la
matriz
Para obtener la descomposición de Hessenb=
erg
En este caso, dado un vector inicial
Figura 3. Algoritmo de Arnoldi
=
Fuente:=
Elaboración
propia.
Los pasos del (4) al (6) del algoritmo
=
llamada factorización de Arnoldi de
orden
La precisión y la velocidad de convergencia de este tipo de mé=
todo
mejora con el incremento de
También se
suelen hacer implementaciones en bloque del método de =
Arnoldi
con un tamaño de bloque que permita agrupar autovalores, siempre evitando
calcular la factorización completa, ya que en las aplicaciones de la
continuación numérica basta con conocer unos cuantos autovalores para saber=
si
estamos en presencia de un punto crítico o inestable, y lo que debemos aseg=
urar
es la obtención de una solución numéricamente estable del sistema de ecuaci=
ones
lineales que se resuelve en cada iteración del método de Newton. Ver por
ejemplo los trabajos de Calvetti =
&
Reichel <=
span
style=3D'font-size:12.0pt;line-height:115%;font-family:"Times New Roman",se=
rif;
color:black'>(2000) y Baglama et al. (1998).
Con las ideas mencionadas previamente se diseñó y programó en
Matlab un algoritmo basado en el método de Arnoldi y
ajustado para la continuación numérica de un conjunto de soluciones
provenientes del modelaje de una reacción química. La medida
Figura 4. Continuación con pseudo-longitud de arco
=
Fuente:
Algoritmo desarrollado en la investigación.
El algoritmo de la Figura 4 fue
implementado en Matlab y comparado con Cl-MatCont,
que es un módulo (toolbox) de Matlab, ver Bindel et al. (2008), que resuelve problemas de
continuación numérica. Se generaron los datos del problema de reacción quím=
ica
conocido como Brusselator, descrito en detalle =
por Prigogine y Lefever
(1974), el cual ha sido usado para pruebas
experimentales en estudios de continuación. La reacción química que da orig=
en a
este problema es oscilante con puntos de bifurcación.
Con la intención de comparar ambas
implementaciones en cuanto a uso de memoria y tiempo de cómputo para hallar=
la
solución, se consideraron mallas de discretización del problema seleccionad=
o, que
generaban matrices Jacobianas de tamaño variado, entre
Cuadro 1. Comparación Cl_Matcom(CM) y Continuación con Arnoldi
(CA)
Fuente: Elaboración propia.
Aunque las mallas definidas para este
problema pueden considerarse pequeñas, la generación de soluciones intermed=
ias
genera un consumo de memoria que llega a ser muy grande. Podemos observar e=
n la
tabla anterior que la implementación del algoritmo que usa un módulo de
Es importante señalar que, en términos=
de
precisión, tanto Cl-MatCont como la implementac=
ión PA
presentada como ejemplo, satisfacen el objetivo de la continuación numérica,
hallando diversas ramas de soluciones y detectando los puntos de interés, t=
al
como se muestra en la Figura 5.
=
Figura
5. Continuación
con pseudo-longitud de arco
Fuente: Elabora=
ción
propia.
Los resultados obteni=
dos
solo nos permiten mostrar el uso de la herramienta, y no son concluyentes en
términos generales; de hecho, otros autores como han logrado excelentes
resultados en otras aplicaciones usando otros métodos iterativos para
continuación, o diferentes métodos basados en espacios de Krylov,
ver por ejemplo los trabajos de Baglama et al. =
(1998)
y Calvetti & Reichel (=
2000).
El análisis de bifurcación es una herr=
amienta
importante en muchas aplicaciones que modelan procesos industriales, y una =
de
las técnicas más usadas para este análisis es la continuación numérica, que
consiste en seguir un conjunto de soluciones de una ecuación no lineal
Hemos presentado un ejemplo usando el
método de Arnoldi para el cálculo de autovalore=
s y
también para resolver estos sistemas lineales; sin embargo, pudiéramos usar
algún otro método de Krylov, ya que estos métod=
os
están en capacidad de resolver
El tamaño de la base de Krylov es un asunto que debe ser estudiado, ya que ex=
iste
un compromiso entre la precisión de la solución y el uso de recursos
computacionales como la memoria y el tiempo de procesamiento. Un algoritmo =
en
bloques basado en espacios de Krylov sería una =
buena
alternativa cuando se establece la característica que se desea evaluar de l=
as
soluciones; por ejemplo, si se desea evaluar la estabilidad, podría definir=
se
un bloque de tamaño
El presente documento estuvo dirigido =
a:
1)
Presentar la generali=
dad
de los métodos de continuación numérica, su utilidad e implementación básic=
a en
términos de un esquema predictor-corrector.
2)
Presentar las princip=
ales
características de los métodos de Krylov para el
cálculo de algunos autovalores y autovectores de
matrices de gran magnitud.
3)
Señalar la utilidad de
los métodos basados en espacios de Krylov, como=
los
de Arnoldi, cuando se usan en conjunto con un m=
étodo
de continuación numérica.
4)
Incluir algoritmos
básicos que pueden usarse como guía para diseñar un método de continuación.=
5)
Establecer una
comparación con un ejemplo práctico, para resaltar y soportar la propuesta =
de
uso de métodos tipo Arnoldi en la continuación
numérica de soluciones a problemas provenientes de aplicaciones en las cien=
cias
e ingeniería.
6)
Instruir en forma gen=
eral
con respecto a la algoritmia detrás de los métodos de continuación presenta=
dos.
Baglama, J., Calvetti, D.,
Bindel,
D., Friedman, M., Govaerts, W., Hughes, J., &am=
p; Kuznetsov, Y. A. (2008). CL MATCONTL: a continuation
toolbox for large equilibrium problems in MATLAB.
Calvetti,
D., & Reichel, L. (1999). A block-Lanczos m=
ethod
for large continuation problems. Numerical Algorithms, 21(1-4),
109-118.
Calvetti,
D., & Reichel, L. (2000). Iterative methods for large continuation
problems. Journal of computational and applied mathematics, 123(1=
-2),
217-240.
Castillo, Z. (2004). A new algorithm for
continuation and bifurcation analysis of large scale free surface flows.
Demmel=
, J.
W. (1997). Applied numerical linear algebra. Society for Industrial and
Applied Mathematics. Berkeley, California, first edition.
Golub, G. H., & Van Loan, C. F. (1996). Matrix
Computations. 3rd. edn. ed. The Johns Hopkins
University.
Govaerts=
, W.
J. (2000). Numerical methods for bifurcations of dynamical equilibria.
Society for Industrial and Applied Mathematics.
Prigogine, I., & Lefever, R.
(1974). Stability and self-organization in open systems (pp.
215-240). Wiley, New York.
Saad, Y. (2011). Numerical methods for large
eigenvalue problems: revised edition. Society for Industrial and Applied
Mathematics.
Seydel=
, R.
(2009). Practical bifurcation and stability analysis (Vol. 5).
Springer Science & Business Media.
Sorensen, D. C. (1997). Implicitly restarted Arnoldi/Lanczos methods f=
or large
scale eigenvalue calculations. In Parallel Numerical Algorithms (=
pp.
119-165). Springer, Dordrecht.
PARA CITAR EL ARTÍCULO INDEXADO.=
Castillo Marrero, Z. N., Colmenares Pacheco, G. A.,
& Cevallos Vique, V. O. (2020). Sobre el us=
o de
métodos de Arnoldi para la continuación numéric=
a de
puntos estacionarios. Ciencia Digital, 4(3), 378-390. https://doi.org/=
10.33262/cienciadigital.v4i3.1385
=
El artículo que se publica es de exclusiva responsabilid=
ad
de los autores y no necesariamente reflejan el pensamiento de la Revista Ciencia Digital.
El artículo queda en propiedad de la revista y, por tant=
o,
su publicación parcial y/o total en otro medio tiene que ser autorizado por=
el
director de la Revista Ciencia Dig=
ital.
[1] Escue=
la
Superior Politécnica de Chimborazo, Facultad de Ciencias, Escuela de
Estadística, Grupo CIDED. Riobamba, Ecuador, zenaida.castillo@espoch.edu.ec=
[2] Unive=
rsidad
de Investigación Experimental Yachay, Escuela de Matemáticas y Ciencias
Computacionales. Urcuquí, Ecuador, gcolmenares@yachaytech.edu.ec
[3] Escue=
la
Superior Politécnica de Chimborazo, Facultad de Administración de Empresas,
Escuela de Finanzas, Grupo CIDED, Riobamba, Ecuador,
victor.cevallos@espoch.edu.ec
www.cienciadigital.org