2019
242
07122018
29062019
Resumen: El problema de los valores propios cuadráticos ha captado el interés de muchos investigadores debido a sus numerosas aplicaciones en áreas tales como mecánica de fluidos, Acústica e Ingeniería, entre otras. En este artículo se presenta un algoritmo híbrido para resolver este problema, el cual combina un algoritmo para el cálculo de soluciones de la ecuación cuadrática matricial con uno tradicional para el cálculo de valores propios. Se seleccionaron doce problemas de aplicación, en diversas áreas, para hacer un análisis comparativo del desempeño del algoritmo propuesto con el tradicional método de Newton-Schur. Como conclusión se determinó que la propuesta algorítmica es competitiva frente al método de Newton-Schur y más eficiente.
Palabras clave:convergencia q-cuadráticaconvergencia q-cuadrática,función cuadrática matricialfunción cuadrática matricial,Método cuasi-NewtonMétodo cuasi-Newton,método de Newtonmétodo de Newton,problema de valores propios cuadráticoproblema de valores propios cuadrático.
Abstract: The quadratic eigenvalue problem has attracted the interest of many researchers due to its numerous applications in areas such as fluid mechanics, acoustics and engineering, among others. In this paper we present a hybrid algorithm to solve this problem, which combines an algorithm for the calculation of solutions of the matrix quadratic equation with a traditional one for the calculation of eigenvalues. Twelve application problems were selected, in different areas, to make a comparative analysis of the performance of the proposed algorithm with the traditional Newton method. In conclusion, it was determined that the algorithmic proposal is competitive with Newton's method and more efficient.
Keywords: Newton method, quadratic matrix function, quasi-Newton method, q-quadratic convergence, the quadratic eigenvalue problem.
Ciencias Básicas
Un algoritmo híbrido para resolver el problema de los valores propios cuadrático
A hybrid algorithm to solve the quadratic eigenvalue problem
Recepción: 07 Diciembre 2018
Aprobación: 29 Junio 2019
En una gran variedad de aplicaciones es necesario encontrar, si es posible, soluciones o de la llamada ecuación cuadrática matricial Ec.1.
(1)donde 
y Q es una función Fréchet diferenciable.
Entre esas aplicaciones están problemas estocásticos [1], análisis dinámico de mecánica estructural y sistemas acústicos; simulación de circuitos eléctricos; mecánica de fluidos; procesamiento de señales; problemas de vibración y modelamiento de sistemas mecánicos micro electrónicos [2] [3] [4].
En gran parte, el estudio y los desarrollos teóricos de la ecuación (1) han sido motivados por un problema que surge con frecuencia en el análisis de sistemas estructurales y en problemas de vibración [5] [6] [7] conocido como el problema de los valores propios cuadrático, el cual consiste en encontrar números complejos λ y vectores no nulos z tales que Ec.2.
(2)en este caso, z es llamado un vector propio de la matriz 
correspondiente al valor propio λ. La matriz 
, cuyas componentes son polinomios en λ de grado menor o igual que dos, es conocida como una matriz λ [8] [9] [10].
Una idea útil para resolver el problema de los valores propios cuadrático es la siguiente [11] [12]. Si la matriz X es un solvente de (1) entonces 
con lo cual 
Así mediante una manipulación algebraica, la matriz 
puede factorizarse como en (3),
(3)A partir de (3), se tiene que los valores propios de 
son los valores propios del solvente junto con los del problema de valores propios generalizado 
[13].
Desde un punto de vista algorítmico, la idea expuesta en el párrafo anterior plantea la necesidad de disponer de un algoritmo híbrido que combine un algoritmo que calcule una solución de (1), con un algoritmo que resuelva el problema de valores propios tanto estándar como generalizado.
Como se evidencia en la factorización (3), un punto clave en la solución del problema de los valores propios cuadrático es obtener una solución de la ecuación cuadrática matricial. Respecto a la solución numérica de esta ecuación, se destaca el método de Newton que ha tenido un papel central en la solución de problemas matriciales no lineales y, en particular, en la solución de la ecuación matricial. Su iteración básica para resolver (1), es la siguiente,
donde 
denota la derivada Fréchet de Q en en XK la dirección de SK .
Observe que, en cada iteración, se debe resolver un caso particular de la ecuación conocida como la ecuación 
de Silvester generalizada [14] [11] para encontrar la matriz 
La forma usual de resolver este tipo de ecuación es a través de descomposición generalizada de Schur [13], la cual es muy costosa computacionalmente [11]. Esta versión del método de Newton para resolver (1), se conoce como el método de Newton-Schur.
Una alternativa que disminuye el costo computacional y el tiempo de ejecución del método de Newton-Schur fue presentada en [12] y consiste en un método cuasi-Newton cuya iteración básica es dada por
Una globalización de este algoritmo fue presentada recientemente en [15], en donde se demuestra que la estrategia de globalización usada no interfiere en la convergencia del algoritmo cuasi-Newton.
Teniendo en cuenta la factorización de la matriz Q(λ) dada en (3) y motivados no solo por las numerosas aplicaciones del problema de valores propios cuadrático, sino por el buen desempeño del algoritmo cuasi-Newton global propuesto recientemente en [15] para encontrar solventes de la ecuación cuadrática matricial, se presenta en este artículo un algoritmo híbrido para resolver dicho problema y se presenta un estudio numérico detallado del mismo, en el cual se consideran doce problemas de aplicación que surgen en diferentes contextos de la ingeniería.
Organizamos este artículo en la siguiente forma. En la Sección 2, se presenta el algoritmo híbrido para resolver la ecuación (1) y con el propósito de comparar su desempeño numérico, se presenta una versión híbrida del método de Newton-Schur. En la Sección 3, se hace una breve descripción de cada uno de los doce problemas de aplicación considerados. En la Sección 4, se analiza numéricamente el comportamiento del algoritmo propuesto. Finalmente, en la Sección 6, se presentan algunas conclusiones del trabajo realizado.
En esta sección, se presenta un algoritmo híbrido para resolver el problema de los valores propios cuadráticos, en el cual combina el algoritmo tipo cuasi-Newton globalizado propuesto en [15] para resolver la ecuación cuadrática matricial con un algoritmo de Matlab que calcula valores propios. El algoritmo híbrido propuesto (Algoritmo 1), es el siguiente.
Con el propósito de comparar el desempeño global de nuestra propuesta algorítmica, se implementa también el algoritmo global de Newton-Schur híbrido (Algoritmo 2) basado en la propuesta que hacen los autores en [11].
Para el análisis comparativo del desempeño numérico del algoritmo propuesto, se consideran doce problemas de aplicación. Estos se describen brevemente a continuación.
En [16], los autores presentan un estudio de auto estabilidad de una bicicleta y para ello, toman como referencia un modelo que consta de cuatro partes rígidas, lateralmente simétricas e idealmente articuladas: dos ruedas, un marco y un ensamblaje frontal para el cual, muestran que el modelo de ecuaciones linealizadas que describe el movimiento de esta clase de bicicleta puede escribirse como el siguiente sistema de ecuaciones diferenciales,
(4)donde las componentes constantes de las matrices M,C1 ,Ko y K2 son definidas en términos de 25 parámetros de diseño [15]. M es una matriz de masa, simétrica, la cual da la energía cinética del sistema de la bicicleta, la matriz 
es lineal en la velocidad de avance y almacena los pares de torsión debido a las velocidades de dirección y de inclinación; la matriz de 
es la suma de dos partes: una parte simétrica, 
, independiente de la velocidad, proporcional a la aceleración gravitacional, y una parte no simétrica 
, cuadrática en la velocidad de avance.
La ecuación diferencial homogénea asociada a (4) conduce a la correspondiente ecuación característica
(5)cuyas soluciones λ, los valores propios de 
dan información sobre el tipo de movimiento. En efecto, valores propios con parte real positiva corresponden a movimientos inestables; con parte real negativa, a movimientos estables asintóticamente y valores propios imaginarios corresponden a movimientos oscilatorios. En particular, las matrices M,C1,Ko y k2 son las siguientes,
En [17], los autores presentan un modelo matemático para sistemas ecológicos que incorpora la variabilidad del ambiente en procesos de nacimiento y muerte, lo cual hace que dicho modelo sea apropiado para el estudio de la evolución de muchas poblaciones. Al incorporar una variable auxiliar, el modelo se convierte en uno de casi nacimiento y muerte.
El modelo Bilby, es un ejemplo del modelo descrito en el párrafo anterior, donde la población considerada está formada por pequeños marsupiales australianos conocidos como bilbies, los cuales son una especie en peligro de extinción. En el modelo se tiene en cuenta que los patrones de reproducción de los bilbies dependen de la calidad de las estaciones precedentes; una medida de esta calidad es la cantidad de lluvia que hubo durante esa estación; así, el modelo bilby incorpora una variable que se específica, si una estación es de buena o de mala calidad.
Algunas de la propiedades cuadrática claves del modelo bilby son capturadas en la solución de la ecuación matricial (1) y el correspondiente polinomio cuadrática matricial, 
con
donde,
Los vectores b y d son vectores de probabilidad y es un vector de unos.
En este problema se considera una viga delgada con soportes simples en ambos extremos, y amortiguada en el punto medio. La ecuación del movimiento que da el desplazamiento transversal de la viga es una ecuación diferencial en derivadas parciales, con condiciones de frontera, cuya discretización mediante el método de elementos finitos es un problema de valores propios cuadráticos, con
donde, la matriz de masa resultante A y la matriz de rigidez C son simétricas y definidas positivas, por construcción; en tanto que B es simétrica y semidefinida positiva. Esto implica que todos los valores propios de 
pertenecen al semiplano izquierdo (cerrado) lo cual garantiza que el problema de la viga es estable (débilmente) [21].
Encontrar los valores propios de la matriz 
donde 
(matrices definidas positivas) y 
(matriz semidefinida positiva). En particular, 
y 
En este problema se considera un polinomio matricial cuadrático 
el cual es hiperbólico, es decir, las matrices 
y 
son hermitianas,
es definida positiva y, para todo 
se satisface,
Estos polinomios desempeñan un papel central en sistemas de vibración [7]. En [19] los autores generan polinomios de este tipo a partir de un conjunto de valores y vectores propios reales 
, para lo cual, definen las matrices,
con 
no singular y para alguna matriz ortogonal Con esta información se tiene que la matriz simétrica 
es hiperbólica, donde
y tiene valores y vectores propios 
.
En particular tomamos n=3 ; los valores propios son 
y los vectores propios asociados son
En este problema se considera un circuito de prueba que representa un circuito equivalente de onda completa para una tira de metal pequeña el cual, está discretizado en dos. El modelo matemático de este problema conduce a la siguiente ecuación diferencial con retardo 
donde
Las ecuaciones diferenciales de retardo ocurren en diferentes campos incluyendo teoría de circuitos. En efecto, circuitos que incluyen elementos de retardo han incrementado su importancia debido a que aumentan en el rendimiento del sistema electrónico respectivo. Una vez modelado el problema un aspecto clave es el estudio de la estabilidad del mismo. En dicho estudio usando el método propuesto en [21] [22] conduce a un problema de valores propios cuadráticos 
con
Una cámara omnidireccional es la que tiene un ángulo de visión mayor que 180 grados. Cámaras con tal ángulo de visión son ventajosas en aplicaciones como vigilancia, seguimiento, estructuras en movimiento y navegación; la tecnología de estas cámaras es muy usada en algunos campos como la robótica. Actualmente, a raíz de la aparición de nuevos sistemas de reproducción de realidad virtual se está introduciendo en el sector del entretenimiento audiovisual.
En [23] los autores proponen un modelo de cámara omnidireccional, el cual asigna un punto de la escena a su imagen digitalizada por una matriz de proyección en perspectiva. El número de parámetros en este modelo permite formular el proceso de calibración de la cámara omnidireccional como un problema de valores propios cuadrático con matriz 
de tamaño 9x9, donde la matriz 
tiene una columna no nula, 
tiene 5 columnas no nulas y de rango 5, y 
es una matriz de rango completo,
Este problema surge del análisis de estabilidad de una ecuación parcial de retardo diferencial (PDDE) [21], [22]. Su discretización da lugar a un sistema de retardo de tiempo
donde,
es tridiagonal y 
son diagonales con
Aquí,
y 
son parámetros reales y 
es el número de puntos de cuadrícula interiores espaciados uniformemente en la discretización del PDDE. Preguntando por los retrasos 
;
tal que el sistema de retardo es estable conduce al problema de valor propio cuadrático 
de dimensión 
donde
donde i es la unidad imaginaria y 
es un parámetro. Los valores tomados son:
En este problema, el comportamiento dinámico de una planta de energía nuclear simplificada en un sistema de 8 grados de libertad es descrito por un problema de valores propios cuadrático donde la matriz de masa M y la matriz de amortiguamiento D son reales, simétricas y la matriz de rigidez K tiene la forma 
con K0una matriz real simétrica. El parámetro 
describe la amortiguación histerética (cuando un elemento estructural es sometido a inversiones en el sentido de la carga aplicada debido a que el material del elemento se encuentra en rango inelástico o no lineal) del problema.
En este problema se considera una modelo de vía férrea soportada elásticamente (por “traviesas”) en puntos espaciados uniformemente. El modelo admite mecanismos de amortiguación externos e internos. El largo del riel es mucho más grande que h, y las deformaciones elásticas del riel deben ser investigadas.
El modelo está formulado haciendo una discretización adecuada de un problema continuo comparable gobernado por una ecuación diferencial en derivadas parciales. En esta discretización se muestra que las oscilaciones de la vía férrea están descritas por los valores propios cuadrático cuya matriz 
está dada por
donde A es una matriz circulante cuya primera fila es 
En este problema, se considera un sistema masa-resorte amortiguador [10], donde la masa i-ésima de peso 
está conectada a su 
vecino por un resorte y un amortiguador con constantes 
y 
respectivamente. La masa i-ésima también está conectada a tierra por un resorte y un amortiguador con constantes 
y 
respectivamente (ver Fig. 1).
El estudio de la vibración de este sistema conduce al problema de valores propios cuadrático 
Para una escogencia por defecto de los parámetros, la matrices 
están dadas por,
En este problema se considera un modelo del ala de un avión, la cual está colocada en un pequeño ángulo de incidencia con una corriente de aire. Los tres grados de libertad considerados son, la flexibilidad del ala, giro del ala y el movimiento del alerón relativo al del ala. Los modos de movimiento y las columnas modales, son las soluciones al problema de valores propios cuadrático 
En el caso particular en el que la velocidad del viento es 12 pies por segundo, las matrices de coeficientes dinámicas son las siguientes.
Para entender y predecir el nivel de sonido en un sistema acústico es de gran importancia conocer los valores propios del sistema, conocidos en este contexto como frecuencias de respuesta. En particular, modificando estas frecuencias de respuesta se pueden evitar resonancias indeseables en una cavidad.
Dado un revestimiento es posible predecir las frecuencias de respuesta modelando el comportamiento del sistema por medio de la ecuación de la onda acústica harmónica en el tiempo
para una presión acústica p en un dominio acotado con condiciones de frontera y de impedancia apropiadas [26]. En esta ecuación, f es la frecuencia, c es la velocidad del sonido en el medio y es la impedancia, posiblemente compleja.
Los valores propios del sistema discretizado [27] pueden ser calculados para una impedancia dada resolviendo el problema de valores propios cuadrático
Los valores propios de 
son las frecuencias de resonancia del sistema.
Sobre el dominio 
con un tamaño de malla h y las matrices de coeficientes 
de tamaño 
con 
están definidas por,
donde 
denota el producto Kronecker,
representa la impedancia y 
y 
son matrices definidas, respectivamente, por
En esta sección, se presenta un análisis numérico comparativo del algoritmo híbrido propuesto comparado con el Algoritmo 2. Para ello, consideramos los doce problemas de aplicación descritos en la sección anterior.
Para los códigos de los algoritmos y de las funciones que definen los problemas usamos MATLAB; los experimentos numéricos se desarrollaron en un computador INTEL (R) Core (TM) i5-3450 de 2.8 GHz.
La Tabla I, contiene los problemas considerados, su tamaño y las matrices iniciales (aleatorias), en cada caso.
En esta Tabla, i representa 
rand es un número aleatorio entre 0 y 1.
Para efectos de la comparación numérica se procedió a realizar 100 pruebas numéricas para cada algoritmo con la misma matriz inicial; se encontró el número de iteraciones, se promedió el tiempo computacional obtenido y el número de éxitos con cada algoritmo.
La Fig. 1, muestra el tiempo de ejecución de cada algoritmo para cada problema considerado. Se observa que en todos los casos el algoritmo híbrido propuesto (Algoritmo 1) tiene un menor costo computacional con respecto al método Newton-Schur híbrido (Algoritmo 2). En particular para los problemas P1 y P2 el comportamiento es casi similar, pero podemos apreciar un excelente comportamiento del Algoritmo 1, en los demás casos, en particular, para los problemas P3 y P10.
La Fig. 2, ilustra el comportamiento de los dos algoritmos respecto al número de iteraciones.
En general, se puede apreciar que el Algoritmo 2 tiene un mejor comportamiento que el Algoritmo, en cuanto al número de iteraciones se refiere. Vale la pena mencionar que, para el problema P7 el número de iteraciones con el Algoritmo 1 es muy alto comparado con el Algoritmo 2, pero su costo computacional es menor (ver Fig. 1).
La Fig. 3 relaciona el porcentaje de éxitos para cada problema.
Este gráfico permite apreciar el buen comportamiento de los dos algoritmos, ya que el nivel de éxito fue prácticamente del 100%.
Finalmente, se presentan las gráficas de la ubicación, en el plano complejo, de los valores propios obtenidos por el algoritmo híbrido propuesto, para cada uno de los problemas Fig.(4.5.6.7.8.9.10.11.12,13,14,15).
Dadas las numerosas aplicaciones en las que surge el problema de los valores propios cuadrático y la necesidad de disponer de algoritmos eficientes que lo resuelvan, presentamos un nuevo algoritmo para su solución. Este es un algoritmo híbrido que combina un algoritmo para el cálculo de soluciones de la ecuación cuadrática matricial con un algoritmo tradicional para el cálculo de valores propios.
Las pruebas numéricas realizadas muestran que nuestra propuesta algorítmica es competitiva frente al tradicional método de Newton-Schur y más eficiente, haciendo de la propuesta híbrida una nueva herramienta computacional para resolver el problema de valores propios cuadrático.
Profesora Titular del Departamento de Matemáticas de la Universidad del Cauca, integrante y coordinadora del Grupo de Investigación OPTIMIZACIÓN. Licenciada en matemáticas de la Universidad del Cauca; Magíster en Matemáticas de la Universidad del Valle; Doctora en Matemática Aplicada de la Universidad de Campinas, SP, Brasil. Ha sido coordinadora de la Maestría en Ciencias Matemáticas de la Universidad del Cauca. Ha coordinado proyectos de investigación en optimización, álgebra lineal numérica y redes neuronales, y ha dirigido trabajos de grado a nivel de pregrado y maestría.
El autor recibió el título Licenciado en Matemáticas en el año 2000, es Especialista en Matemática Aplicada 2001, Magíster en Ciencias Matemáticas 2013 y actualmente estudiante de Doctorado en Ciencias Matemáticas en la Universidad del Cauca. Pertenece al grupo de Optimización, grupo de investigación de la Universidad del Cauca. Sus áreas de interés son el Análisis Numérico, la Optimización Matemática y el Álgebra Lineal Numérica.
Profesor titular del Departamento de Matemáticas de la Universidad del Cauca y miembro del grupo de Optimización de la Universidad del Cauca, Licenciado en Matemáticas de La Universidad de Sucre, Especialista en Matemáticas de la Universidad de Córdoba y Magíster en Ciencias Matemáticas de la Universidad del Cauca, ha trabajado en el grupo interesado en los temas: sistemas de ecuaciones no lineales y el problema de complementariedad no Lineal.
En esta Tabla, i representa rand es un número aleatorio entre 0 y 1.
