Métodos iterativos de resolución de sistemas lineales basados en subespacios de Krylov. Aplicaciones a la resolución numérica de EDP y en matrices test.

FACULTAD DE MATEMÁTICAS TRABAJO FIN DE GRADO Métodos iterativos de resolución de sistemas lineales basados en subespacios de Krylov. Aplicaciones a la resolución numérica de EDP y en matrices test. Tutor: Samuele Rubino Marta Guijo Moreno Grado en Matemáticas Curso 2022/2023 Resumen Este trabajo se centra en los diferentes métodos iterativos para resolver sistemas lineales basados en subespacios de Krylov. El objetivo principal es realizar una comparativa entre los distintos métodos a la hora de resolver numéricamente un problema de EDP o problemas test con matrices grandes y dispersas. Para ello, introducimos en el capítulo 1 conceptos ya vistos en Cálculo Numérico II necesarios para poder definir los siguientes métodos. En el capítulo 2 definimos métodos de proyección y en el 3 nos centramos en los métodos iterativos de Krylov. Posteriormente, en el capítulo 4, consideramos un problema de EDP, su discretización e implementación en Matlab para la comparación de dichos métodos. Vamos a ver la diferencia en sistemas precondicionados y no precondicionados en el método GMRES considerando problemas test con matrices grandes y dispersas y por último, analizamos el resultado obtenido. Abstract The present work focuses on different iterative methods for solving linear systems based on Krylov subspaces. The main objective is to make a comparison among the different methods when solving a PDE problem numerically or test problems with large and sparse arrays. To do so, we introduce in Chapter 1 concepts already covered in the course Numerical Calculus II (Cálculo Numérico II), which are required to define the following methods. In Chapter 2, the projection methods are defined and in Chapter 3 we focus on Krylov’s iterative methods. Subsequently, in Chapter 4, we will consider a PDE problem, its discretization and implementation in MATLAB for the comparison of these methods. We will consider the difference in preconditioned systems and non-preconditioned systems in the GMRES method considering test problems with large and sparse arrays and finally, the result obtained will be analyzed. Índice general Introducción iii 1. Métodos iterativos básicos 1 1.1. Introducción........................................ 1 1.2. Métodos de Jacobi, Gauss-Seidel y relajación . . . . . . . . . . . . . . . . . . . . . . 4 1.2.1. MétododeJacobi................................. 5 1.2.2. Método de Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.2.3. Métododerelajación............................... 7 2. Métodos de proyección para la resolución de sistemas lineales 11 2.1. Definiciones básicas y algoritmos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.1.1. Métodos generales de proyección . . . . . . . . . . . . . . . . . . . . . . . . 12 2.1.2. Representación matricial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2. Teoríageneral ....................................... 14 2.2.1. Dos resultados de optimalidad . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2.2. Interpretación en términos de proyectores . . . . . . . . . . . . . . . . . . . 15 2.2.3. Límite de error general . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3. Procesos de proyección unidimensional . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.3.1. Descensomásrápido............................... 19 2.3.2. IteraciónMR ................................... 21 2.3.3. Norma residual del descenso más rápido . . . . . . . . . . . . . . . . . . . . 22 3. Métodos iterativos de Krylov 23 3.1. SubespaciodeKrylov.................................... 23 ii trabajo fin de grado 3.2. ElmétododeArnoldi ................................... 25 3.2.1. El método de Arnoldi para sistemas lineales . . . . . . . . . . . . . . . . . . 25 3.3. El método del residual mínimo generalizado (GMRES) . . . . . . . . . . . . . . . . . 26 3.3.1. Problemas de implementación práctica . . . . . . . . . . . . . . . . . . . . . 28 3.4. Algoritmo simétrico de Lanczos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.4.1. Elalgoritmo.................................... 31 3.4.2. Relación con los polinomios ortogonales . . . . . . . . . . . . . . . . . . . . 32 3.5. Algoritmo del gradiente conjugado . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.5.1. Derivaciónyteoría................................ 33 4. Aplicación a la resolución numérica de EDP y en matrices test 37 4.1. Aplicación a la discretización de EDP lineales elípticas . . . . . . . . . . . . . . . . 37 4.1.1. Ejemplo1 ..................................... 40 4.2. Aplicación en matrices test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.2.1. Precondicionamiento............................... 42 4.2.2. Ejemplo2 ..................................... 45 4.2.3. Ejemplo3 ..................................... 45 4.3. Conclusiones........................................ 46 Bibliografía 47 Introducción En este trabajo nos centramos en los métodos iterativos de resolución de sistemas lineales. Recordamos que los métodos directos de resolución de sistemas lineales (vistos en la asignatura Cálculo Numérico I) proporcionan la solución exacta 𝑥de un sistema lineal en un número finito de pasos. A diferencia de los métodos directos, los métodos iterativos conducen a la solución exacta después de un número infinito de pasos. Estos métodos iterativos pueden ser mejores que los métodos directos en algunas situaciones particulares, por ejemplo, para sistemas lineales grandes y, con frecuencia, muy eficientes en el caso de matrices huecas (matrices en las que la mayor parte de sus elementos son cero). En el capítulo 1 describimos los métodos iterativos de resolución de sistemas lineales vistos en Cálculo Numérico II. Los métodos iterativos son apropiados para sistemas lineales grandes y matrices huecas cuando se discretizan EDP (e.g., ecuación de Poisson, convección-reacción-difusión, etc.) mediante, por ejemplo, el método de diferencias finitas o el método de elementos finitos. Comenzamos el capítulo con conceptos sobre la convergencia de dichos métodos y las estimaciones a priori y a posteriori. Repasamos en particular los métodos de Jacobi, Gauss-Seidel y Relajación. A lo largo del capítulo 2 definimos métodos de proyección para la resolución de sistemas lineales con el fin de obtener una solución aproximada del problema planteado para un subespacio de ℝ𝑛. También exponemos dos resultados de optimalidad donde obtendremos la condición de Petrov-Galerkin. A continuación, desarrollamos varios procesos de proyección unidimensional como el algoritmo del descenso más rápido cuyo resultado es una consecuencia de la Desigualdad de Kantorovich; el algoritmo del residuo mínimo; y la norma residual del descenso más rápido. En el capítulo 3 describimos métodos iterativos de Krylov introduciendo los subespacios de Krylov que denotaremos 𝑚. A continuación, vemos un método de proyección ortogonal sobre 𝑚 llamado método de Arnoldi, para matrices no hermitianas. También describimos el método del residual mínimo generalizado (GMRES) en el que nos centramos resolviendo el problema de mínimos cuadrados transformando la matriz de Hessenberg en una de forma triangular superior. A continuación, describimos el algoritmo de Lanczos que se puede ver como una simplificación del método de Arnoldi cuando la matriz es simétrica. Y por último, describimos el algoritmo del gradiente conjugado que es una de las mejores técnicas iterativas conocidas para resolver sistemas lineales definidos positivos, simétricos y dispersos. 6trabajo fin de grado Observación 1.4. 1) Los pasos a seguir para el cálculo de 𝑢𝑘+1 𝑖son los siguientes. Si en un vector de memoria 𝑢tenemos guardado 𝑢𝑘, 1. 𝑟←𝑏−𝐴𝑢 (residuo en la etapa k), 2. 𝑢𝑖←𝑢𝑖+𝑟𝑖 𝑎𝑖𝑖 . 2) Para calcular cada componente 𝑢𝑘+1 𝑖se utiliza el vector 𝑢𝑘, con lo que se usan 2 vectores de dimensión n en cada iteración, es decir, 2n registros de memoria. 3) Las n componentes del vector 𝑥𝑘+1 son calculadas de manera similar (así, el método también se conoce como método de iteraciones simultáneas), lo que hace que el método sea fácilmente "paralelizable". 1.2.2 Método de Gauss-Seidel Para el método de Jacobi (1.11) se consideran los valores 𝑢𝑘 1, ..., 𝑢𝑘 𝑖−1 para calcular 𝑢𝑘+1 𝑖. Sin embargo, podemos considerar intuitivamente que las aproximaciones 𝑢𝑘+1 1, ..., 𝑢𝑘+1 𝑖−1 son mejores que 𝑢𝑘 1, ..., 𝑢𝑘 𝑖−1 para calcular 𝑢𝑘+1 𝑖, por lo que se puede definir el siguiente método de iteración de punto fijo en el que se van sustituyendo las antiguas aproximaciones 𝑢𝑘 1, ..., 𝑢𝑘 𝑖−1 por las nuevas 𝑢𝑘+1 1, ..., 𝑢𝑘+1 𝑖−1 conforme se van calculando. El método de Gauss-Seidel para la solución del sistema lineal de ecuaciones 𝐴𝑢 =𝑏viene dado por: 𝑢𝑘+1 𝑖=1 𝑎𝑖𝑖(𝑏𝑖− 𝑖−1 ∑ 𝑗=1 𝑎𝑖𝑗𝑢𝑘+1 𝑗− 𝑛 ∑ 𝑗=𝑖+1 𝑎𝑖𝑗𝑢𝑘 𝑗),1≤𝑖≤𝑛. (1.12) Con la notación anterior, introducimos el método invirtiendo la parte triangular inferior de 𝐴, es decir, tomando 𝑀=𝐷−𝐸y𝑁=𝐹. Con esta elección y, de nuevo, gracias a la hipótesis (1.8) tenemos 𝑑𝑒𝑡(𝐷−𝐸) = 𝑑𝑒𝑡𝐷 ≠0, es decir, se tiene (1.9). El método tiene la forma: {𝑢0∈𝕂𝑛𝑑𝑎𝑑𝑜 𝑢𝑘+1 = (𝐷−𝐸)−1𝐹 𝑢𝑘+ (𝐷−𝐸)−1𝑏, ∀𝑘≥0. Se llamará matriz de Gauss-Seidel a la matriz dada por: 1= (𝐷−𝐸)−1𝐹 . De nuevo, como consecuencia del Teorema 1.1, podemos escribir Corolario 1.5. Sea 𝐴∈(𝕂𝑛)una matriz satisfaciendo 𝑑𝑒𝑡𝐴 ≠0y (1.8), y sea 𝑏∈𝕂𝑛. Entonces, el método de Gauss-Seidel es globalmente convergente si y solo si 𝜌(1)<1. Obvsérvese que el método puede ser reescrito como (𝐷−𝐸)𝑢𝑘+1 =𝐹 𝑢𝑘+𝑏, ∀𝑘≥0. 1. métodos iterativos básicos 7 Así, el cálculo efectivo de la sucesión {𝑢𝑘}𝑘≥1⊂𝕂𝑛asociada al método se lleva a cabo del modo siguiente: 𝐷𝑢𝑘+1 =𝐸𝑢𝑘+1 +𝐹 𝑢𝑘+𝑏, que, matricialmente, se escribe como: ⎛⎜⎜⎜⎜⎝ 𝑎11 0⋯0 0𝑎22 ⋯0 ⋮ ⋮ ⋱ ⋮ 0 0 ⋯𝑎𝑛𝑛 ⎞⎟⎟⎟⎟⎠ ⎛⎜⎜⎜⎜⎝ 𝑢𝑘+1 1 𝑢𝑘+1 2 ⋮ 𝑢𝑘+1 𝑛 ⎞⎟⎟⎟⎟⎠ =⎛⎜⎜⎜⎜⎝ 𝑏1 𝑏2 ⋮ 𝑏𝑛 ⎞⎟⎟⎟⎟⎠ −⎛⎜⎜⎜⎜⎝ 0 0 ⋯0 𝑎21 0⋯0 ⋮ ⋮ ⋱ ⋮ 𝑎𝑛1𝑎𝑛2⋯0 ⎞⎟⎟⎟⎟⎠ ⎛⎜⎜⎜⎜⎝ 𝑢𝑘+1 1 𝑢𝑘+1 2 ⋮ 𝑢𝑘+1 𝑛 ⎞⎟⎟⎟⎟⎠ −⎛⎜⎜⎜⎜⎝ 0𝑎12 ⋯𝑎1𝑛 0 0 ⋯𝑎2𝑛 ⋮ ⋮ ⋱ ⋮ 0 0 ⋯0 ⎞⎟⎟⎟⎟⎠ ⎛⎜⎜⎜⎜⎝ 𝑢𝑘 1 𝑢𝑘 2 ⋮ 𝑢𝑘 𝑛 ⎞⎟⎟⎟⎟⎠ . Así, la primera componente se calcula como 𝑢𝑘+1 1=𝑢𝑘 1+1 𝑎11 [𝑏1− (𝑎11𝑢𝑘 1+𝑎𝑘 12 +... +𝑎1𝑛𝑢𝑘 𝑛)]. Fijemos ahora i: 2≤𝑖≤𝑛. Calculadas las componentes previas del vector 𝑢𝑘+1, es decir, hallados 𝑢𝑗+1 𝑗para j<i, determinamos 𝑢𝑘+1 𝑖como 𝑢𝑘+1 𝑖=1 𝑎𝑖𝑖 [𝑏𝑖−∑ 𝑗<𝑖 𝑎𝑖𝑗𝑢𝑘+1 𝑗−∑ 𝑗>𝑖 𝑎𝑖𝑗𝑢𝑘+1 𝑗]=𝑢𝑘 𝑖+1 𝑎𝑖𝑖 [𝑏𝑖−∑ 𝑗<𝑖 𝑎𝑖𝑗𝑢𝑘+1 𝑗−∑ 𝑗≥𝑖 𝑎𝑖𝑗𝑢𝑘+1 𝑗].(1.13) Observación 1.5. En el método de Gauss-Seidel, el valor 𝑢𝑘 𝑖no interviene en el cálculo de los 𝑢𝑘+1 𝑗, para j > i. Por tanto, podemos reemplazar las posiciones de memoria que almacenan los valores 𝑢𝑘 𝑖 por los nuevos valores 𝑢𝑘+1 𝑖. Eso hace que, para aplicar este método con un programa informático baste utilizar un único vector (n posiciones de memoria) que almacene los distintos vectores 𝑢𝑘que calcula el método. El cálculo de 𝑢𝑘+1 𝑖se hace del modo siguiente: Si en un vector de memoria u tenemos guardado 𝑢𝑘, se hace un ciclo descendente en las componentes, desde 𝑖= 1, ..., 𝑛, calculando 𝑢𝑖←𝑢𝑖+1 𝑎𝑖𝑖 (𝑏𝑖− (𝐴𝑢)𝑖). En este método, además de necesitar menos memoria (solo el vector 𝑢), se "invierte" más parte de la matriz 𝐴que en el de Jacobi. Es razonable pensar que será más rápido. Como desventaja respecto al método de Jacobi, ahora el cálculo de las diferentes componentes no se puede hacer "en paralelo". 1.2.3 Método de relajación A continuación, introducimos el método de relajación para la solución de 𝐴𝑢 =𝑏de la forma: 𝑢𝑘+1 𝑖=𝑢𝑘 𝑖+𝑤 𝑎𝑖𝑖(𝑏𝑖− 𝑛 ∑ 𝑗=1 𝑎𝑖𝑗𝑢𝑘+1 𝑗),1≤𝑖≤𝑛. (1.14) 8trabajo fin de grado Con la notación anterior consideramos la siguiente descomposición de 𝐷: 𝐷=1 𝑤𝐷+ (1 − 1 𝑤)𝐷, 𝑐𝑜𝑛 𝑤 ∈ℝ⧵{0}. De esta forma, podemos escribir 𝐴como: 𝐴=1 𝑤𝐷−𝐸−(1 − 𝑤 𝑤)𝐷−𝐹 y se pueden tomar 𝑀=1 𝑤𝐷−𝐸, 𝑁 =(1 − 𝑤 𝑤)𝐷+𝐹 , de modo que se ha pasado parte de la diagonal 𝐷de 𝐴a la matriz 𝑁. La matriz 𝑀es invertible gracias a la hipótesis (1.8). Por tanto, el sistema lineal (1.1) es equivalente a (1.10) y podemos plantear el correspondiente método iterativo. Este tiene la forma: ⎧ ⎪ ⎨ ⎪ ⎩ 𝑢0∈𝕂𝑛𝑑𝑎𝑑𝑜, 𝑢𝑘+1 =(1 𝑤𝐷−𝐸)−1[(1−𝑤 𝑤𝐷+𝐹)𝑢𝑘+𝑏],∀𝑘≥0 Se llamará matriz de relajación a la matriz dada por 𝑤=(1 𝑤𝐷−𝐸)−1(1 − 𝑤 𝑤𝐷+𝐹)= (𝐷−𝑤𝐸)−1[(1 − 𝑤)𝐷+𝑤𝐹 ]. Con esta notación y, como consecuencia del Teorema 1.1, se tiene: Corolario 1.6. Sean 𝐴∈(𝕂𝑛), una matriz satisfaciendo 𝑑𝑒𝑡𝐴 ≠0y (1.8), 𝑏∈𝕂𝑛y𝑤∈ℝ⧵{0}. Entonces, el método de relajación es globalmente convergente si y solo si 𝜌(𝑤)<1. El cálculo efectivo que se lleva a cabo es el siguiente: (1 𝑤𝐷−𝐸)𝑢𝑘+1 =(1 − 𝑤 𝑤𝐷+𝐹)𝑢𝑘+𝑏⇐⇒ (𝐷−𝑤𝐸)𝑢𝑘+1 = ((1 − 𝑤)𝐷+𝑤𝐹 )𝑢𝑘+𝑤𝑏. De manera equivalente 𝐷𝑢𝑘+1 =𝐷𝑢𝑘+𝑤[𝐸𝑢𝑘+1 − (𝐷−𝐹)𝑢𝑘+𝑏]. En notación vectorial, la igualdad anterior proporciona: ⎛⎜⎜⎜⎜⎝ 𝑎11 0⋯0 0𝑎22 ⋯0 ⋮ ⋮ ⋱ ⋮ 0 0 ⋯𝑎𝑛𝑛 ⎞⎟⎟⎟⎟⎠ ⎛⎜⎜⎜⎜⎝ 𝑢𝑘+1 1 𝑢𝑘+1 2 ⋮ 𝑢𝑘+1 𝑛 ⎞⎟⎟⎟⎟⎠ =⎛⎜⎜⎜⎜⎝ 𝑎11 0⋯0 0𝑎22 ⋯0 ⋮ ⋮ ⋱ ⋮ 0 0 ⋯𝑎𝑛𝑛 ⎞⎟⎟⎟⎟⎠ ⎛⎜⎜⎜⎜⎝ 𝑢𝑘 1 𝑢𝑘 2 ⋮ 𝑢𝑘 𝑛 ⎞⎟⎟⎟⎟⎠ +𝑤[⎛⎜⎜⎜⎜⎝ 𝑏1 𝑏2 ⋮ 𝑏𝑛 ⎞⎟⎟⎟⎟⎠ −⎛⎜⎜⎜⎜⎝ 0 0 ⋯0 𝑎21 0⋯0 ⋮ ⋮ ⋱ ⋮ 𝑎𝑛1𝑎𝑛2⋯0 ⎞⎟⎟⎟⎟⎠ ⎛⎜⎜⎜⎜⎝ 𝑢𝑘+1 1 𝑢𝑘+1 2 ⋮ 𝑢𝑘+1 𝑛 ⎞⎟⎟⎟⎟⎠ −⎛⎜⎜⎜⎜⎝ 𝑎11 𝑎12 ⋯𝑎1𝑛 0𝑎22 ⋯𝑎2𝑛 ⋮ ⋮ ⋱ ⋮ 0 0 ⋯𝑎𝑛𝑛 ⎞⎟⎟⎟⎟⎠ ⎛⎜⎜⎜⎜⎝ 𝑢𝑘 1 𝑢𝑘 2 ⋮ 𝑢𝑘 𝑛 ⎞⎟⎟⎟⎟⎠ ]. 1. métodos iterativos básicos 9 Así pues, la primera componente viene dada por: 𝑢𝑘+1 1=𝑢𝑘 1+𝑤 𝑎11 [𝑏1− (𝑎11𝑢𝑘 1+... +𝑎1𝑛𝑢𝑘 𝑛)]. y, una vez hallados 𝑢𝑘+1 𝑗para j < i, se determina 𝑢𝑘+1 𝑖=𝑢𝑘 𝑖+𝑤 𝑎𝑖𝑖 [𝑏𝑖− (𝑎𝑖1𝑢𝑘+1 1+... +𝑎𝑖,𝑖−1𝑢𝑘+1 𝑖−1 +𝑎𝑖𝑖𝑢𝑘 𝑖+𝑎𝑖,𝑖+1𝑢𝑘 𝑖+1 +... +𝑎𝑖𝑛𝑢𝑘 𝑛)]. Observación 1.6. 1) El método de relajación para 𝑤= 1 coincide con el de Gauss-Seidel. De ahí la notación usada para la matriz de Gauss-Seidel. El cálculo efectivo de una etapa es similar a Gauss- Seidel, para i = 1,...,n, se calcula 𝑢𝑖←𝑢𝑖+𝑤 𝑎𝑖𝑖 (𝑏𝑖− (𝐴𝑢)𝑖). 2) Aunque en principio el parámetro 𝑤podría ser un número real no nulo, para que el método converja, es necesario que 𝑤∈ (0,2). El método se llamará de sobrerrelajación si 𝑤∈ (1,2) y de subrrelajación si 𝑤∈ (0.1). 3) Si 𝜌(𝑤)es una función continua de 𝑤, entonces, el estudio del método consiste en a) Determinar un intervalo 𝐼 ⊂ ℝ⧵{0}, tal que se tenga 𝜌(𝑤)<1. b) Determinar 𝑤0∈𝐼tal que 𝜌(𝑤0) ≈ 𝑖𝑛𝑓𝑤∈𝐼𝜌(𝑤). 4) Para ciertos valores del parámetro de relajación se obtiene una convergencia más rápida que para 𝑤= 1 y por tanto un tiempo de cálculo menor que para el método de Gauss-Seidel. El número de operaciones es similar en ambos métodos. No obstante, hay que tener en cuenta el tiempo utilizado en la estimación preliminar del parámetro 𝑤0para comparar la eficacia de ambos métodos. Capítulo 2 Métodos de proyección para la resolución de sistemas lineales La mayoría de los métodos iterativos existentes para resolver grandes sistemas de ecuaciones lineales utilizan un proceso de proyección. Este proceso representa una forma canónica de extraer una aproximación de la solución de un sistema lineal a partir de un subespacio. Este capítulo describe estas técnicas en un marco muy general. El caso unidimensional está cubierto en detalle al final del capítulo, que proporciona una vista previa de los métodos de proyección más complejos que veremos en el capítulo 3. 2.1 Definiciones básicas y algoritmos Comenzamos con algunas definiciones y algoritmos necesarios para describir los métodos de proyección. Consideramos el sistema lineal 𝐴𝑥 =𝑏, (2.1) donde 𝐴es una matriz real 𝑛𝑥𝑛. La idea de los métodos de proyección es extraer una solución aproximada del problema anterior para un subespacio de ℝ𝕟. Si es este subespacio de candidatos aproximados, también llamado el subespacio buscado, y si 𝑚es su dimensión, entonces, las restricciones de 𝑚deben imponerse para poder extraer dicha aproximación. Un camino típico para describir estas restricciones es imponer a 𝑚(independiente) condiciones ortogonales. Especialmente, el vector residual 𝑏−𝐴𝑥 es restringido para ser ortogonal a 𝑚vectores linealmente independientes. Esto define otro subespacio de dimensión 𝑚, que llamaremos el subespacio de restricciones. Hay dos clases de métodos de proyección: ortogonal y oblicua. En la proyección ortogonal, el subespacio es el mismo que . En el método de proyección oblicua, es diferente de y puede no estar totalmente relacionado con esto. Esta distinción es importante y da lugar a diferentes tipos de algoritmos [7]. 12 trabajo fin de grado 2.1.1 Métodos generales de proyección Sea 𝐴una matriz real 𝑛𝑥𝑛 yydos subespacios m-dimensional de ℝ𝕟. Una técnica de proyección sobre el subespacio y ortogonal a es el proceso que busca una solución aproximada 𝑥 para (2.1) imponiendo las condiciones: 𝐻𝑎𝑙𝑙𝑎𝑟 𝑥 ∈𝑡𝑎𝑙 𝑞𝑢𝑒 𝑏 −𝐴 𝑥 ⊥ .(2.2) Si deseamos aprovechar el conocimiento de una suposición inicial 𝑥0de la solución, luego la aproximación debe buscarse en el espacio afín 𝑥0+en lugar del espacio vectorial homogéneo . Esto requiere una leve modificación en la formulación anterior, 𝐻𝑎𝑙𝑙𝑎𝑟 𝑥 ∈𝑥0+𝑡𝑎𝑙 𝑞𝑢𝑒 𝑏 −𝐴 𝑥 ⊥ .(2.3) Notar que, si escribimos 𝑥 de la forma 𝑥 =𝑥0+𝛿y el vector residual inicial 𝑟0es definido como 𝑟0=𝑏−𝐴𝑥0,(2.4) luego la ecuación anterior se convierte en 𝑏−𝐴(𝑥0+𝛿)⊥ó𝑟0−𝐴𝛿 ⊥ . En otras palabras, la solución aproximada puede ser definida como 𝑥 =𝑥0+𝛿, 𝛿 ∈(2.5) (𝑟0−𝐴𝛿, 𝑤)=0,∀𝑤∈(2.6) Esto es un paso básico de proyección, en la forma más general. La mayoría de las técnicas estándar utilizan una sucesión de tales proyecciones. Normalmente, un nuevo paso de proyección usa un nuevo par de subespacios yy una conjetura inicial 𝑋0igual a la aproximación más reciente obtenida para los pasos de proyección previos. Los métodos de proyección forman una estructura unificadora para muchos de los métodos conocidos en el cálculo científico. De hecho, practicamente todas las técnicas iterativas básicas vistas en el capítulo 1 pueden ser consideradas técnicas de proyección. Cuando una aproximación es definida a través de 𝑚grados de libertad (subespacio )ym restricciones () resulta un proceso de proyección. Los métodos de proyección ortogonales corresponden al caso particular en el que los dos subespacios yson idénticos. La diferencia es particularmente importante en el caso Hermitiano ya que se garantiza que el problema de proyección será Hermitiano en esta situación. Además, muchos resultados teóricos son ciertos para el caso ortogonal. Cuando =, las condiciones de Petrov- Galerkin, que definiremos más adelante (Proposición 2.3), son llamadas a menudo las condiciones de Galerkin. 2.1.2 Representación matricial Sea 𝑉= [𝑣1, ..., 𝑣𝑚], una matriz 𝑛𝑥𝑚 cuyos vectores columna forman una base de , y sea 𝑊= [𝑤1, ..., 𝑤𝑚], una matriz 𝑛𝑥𝑚 cuyos vectores columna forman una base de . Si la solución aproximada es escrita como 𝑥 =𝑥0+𝑉 𝑦, 2. métodos de proyección para la resolución de sistemas lineales 13 luego la condición de ortogonalidad conduce inmediatamente al siguiente sistema de ecuaciones para el vector y: 𝑊𝑇𝐴𝑉 𝑦 =𝑊𝑇𝑟0. Si se hace la suposición de que la matriz 𝑚𝑥𝑚,𝑊𝑇𝐴𝑉 es no singular, la expresión para la solución aproximada resulta ser: 𝑥 =𝑥0+𝑉(𝑊𝑇𝐴𝑉 )−1𝑊𝑇𝑟0(2.7) En algunos algoritmos, la matriz 𝑊𝑇𝐴𝑉 no es necesaria calcularla ya que la obtenemos como subproducto del algoritmo. Un pseudo-código de la técnica de proyección es representado por el siguiente algoritmo. ALGORITMO 2.1. Pseudo-código de la técnica de proyección. 1. Seleccionamos un par de subespacios y. 2. Elegimos las bases 𝑉= [𝑣1, ..., 𝑣𝑚]y𝑊= [𝑤1, ..., 𝑤𝑚]para y 3. 𝑟∶= 𝑏−𝐴𝑥 4. 𝑦∶= (𝑊𝑇𝐴𝑉 )−1𝑊𝑇𝑟 5. 𝑥∶= 𝑥+𝑉 𝑦 La solución aproximada está definida solo cuando la matriz 𝑊𝑇𝐴𝑉 es no singular, una propiedad que no se garantiza que sea verdadera incluso cuando A es no singular. Ejemplo 2.1. Consideramos la matriz 𝐴=(𝑂 𝐼 𝐼 𝐼) donde 𝐼es la matriz identidad 𝑚𝑥𝑚 y𝑂es la matriz cero 𝑚𝑥𝑚, y sea 𝑉=𝑊= [𝑒1, 𝑒2, ..., 𝑒𝑚]la base canónica. Aunque 𝐴es no singular, la matriz 𝑊𝑇𝐴𝑉 es precisamente el bloque 𝑂en la esquina superior izquierda de 𝐴y es entonces singular. Puede ser fácilmente verificado que 𝑊𝑇𝐴𝑉 es no singular si ningún vector del subespacio 𝐴 es ortogonal al subespacio . Hay dos casos particulares importantes donde la no singularidad de 𝑊𝑇𝐴𝑉 está garantizada. Estos son discutidos en la siguiente proposición. Proposición 2.1. Sean A, ysatisfaciendo cualquiera de las dos condiciones siguientes: (i) 𝐴es definida positiva y = (ii) 𝐴es no singular y =𝐴. Luego la matriz 𝐵=𝑊𝑇𝐴𝑉 es no singular para cualquier base 𝑉y𝑊de y, respectivamente. Demostración. Consideramos el primer caso (i). Sea 𝑉alguna base de y sea 𝑊alguna base de . De hecho, ya que yson iguales, 𝑊puede ser siempre expresada como 𝑊=𝑉 𝐺, donde G es una matriz no singular 𝑚𝑥𝑚. Luego 𝐵=𝑊𝑇𝐴𝑉 =𝐺𝑇𝑉𝑇𝐴𝑉 . 14 trabajo fin de grado Como 𝐴es definida positiva, también lo es 𝑉𝑇𝐴𝑉 , que nos demuestra que 𝐵es no singular. Consideramos ahora el caso (ii). Sea 𝑉alguna base de y sea 𝑊alguna base de . Ya que =𝐴, 𝑊puede ser expresado en este caso como 𝑊=𝐴𝑉 𝐺, donde 𝐺es una matriz 𝑚𝑥𝑚 no singular. Luego 𝐵=𝑊𝑇𝐴𝑉 =𝐺𝑇(𝐴𝑉 )𝑇𝐴𝑉 (2.8) Como 𝐴es no singular, la matriz 𝑛𝑥𝑚,𝐴𝑉 es de rango completo y, como consecuencia, (𝐴𝑉 )𝑇𝐴𝑉 es no singular. Esto, junto con (2.8), demuestra que 𝐵es no singular. Ahora consideramos el caso particular donde 𝐴es simétrica (real) y se utiliza una técnica de proyección ortogonal. En esta situación, las mismas bases pueden ser usadas para y, que son subespacios idénticos, y la matriz de proyección, que es 𝐵=𝑉𝑇𝐴𝑉 , es simétrica. Además, si la matriz 𝐴es simétrica definida positiva (SPD), entonces también lo es 𝐵. 2.2 Teoría general En esta sección daremos algunos resultados generales de teoría sin ser específicos sobre los subespacios yque se utilizan. El objetivo es aprender sobre la clase de la aproximación obtenida de un proceso general de proyección. Dos herramientas principales se utilizan para ello. La primera sirve para usar de manera óptima las propiedades de los métodos de proyección. La segunda consiste en interpretar el problema de proyección con la ayuda de los operadores de proyección en un intento de extraer límites residuales [9]. 2.2.1 Dos resultados de optimalidad En esta sección definimos dos resultados importantes de optimalidad considerarando el caso cuando 𝐴es (SPD). Proposición 2.2. Suponiendo que 𝐴es SPD y =, entonces un vector 𝑥 es el resultado de un método de proyección (ortogonal) sobre con el vector inicial 𝑥0si se minimiza la norma de 𝐴 del error sobre 𝑥0+, es decir, si 𝐸(𝑥) = 𝑚𝑖𝑛𝑥∈𝑥0+𝐸(𝑥), donde 𝐸(𝑥)=(𝐴(𝑥∗−𝑥), 𝑥∗−𝑥)1∕2. siendo 𝑥∗la solución exacta. Demostración. Para que 𝑥 sea el mínimo de 𝐸(𝑥), es necesario y suficiente que 𝑥∗−𝑥 sea 𝐴ortogonal para todo el subespacio . Esto produce (𝐴(𝑥∗−𝑥), 𝑣) = 0 ∀𝑣∈ 2. métodos de proyección para la resolución de sistemas lineales 15 o, equivalentemente, (𝑏−𝐴 𝑥, 𝑣) = 0 ∀𝑣∈ que es la condición de Galerkin definiendo un proceso de proyección ortogonal para la aproximación 𝑥. Ahora tomamos el caso donde es definido por =𝐴. Proposición 2.3. Sea 𝐴una matriz cuadrada arbitraria y suponiendo que =𝐴. Entonces, un vector 𝑥 es el resultado de un método de proyección (oblicuo) sobre ortogonalmente a con el vector inicial 𝑥0si se minimiza la norma-2 del vector residual 𝑏−𝐴𝑥 sobre 𝑥∈𝑥0+, es decir, si 𝑅(𝑥) = 𝑚𝑖𝑛𝑥∈𝑥0+𝑅(𝑥), donde 𝑅(𝑥) = ‖𝑏−𝐴𝑥‖2. Demostración. Para que 𝑥 sea el mínimo de 𝑅(𝑥), es necesario y suficiente que 𝑏−𝐴 𝑥 sea ortogonal a todos los vectores de la forma 𝑣=𝐴𝑦, donde 𝑦pertenece a , es decir, (𝑏−𝐴 𝑥, 𝑣) = 0 ∀𝑣∈𝐴, que es precisamente la condición de Petrov-Galerkin que define la solución aproximada 𝑥. 2.2.2 Interpretación en términos de proyectores Volvemos ahora a los dos casos particulares importantes señalados en la sección anterior, es decir, los casos =y=𝐴. En estos casos, el resultado del proceso de proyección puede ser interpretado fácilmente en cuanto a las acciones de las proyecciones ortogonales sobre el residual inicial o error inicial. Consideremos primero el segundo caso, ya que es más simple. Sea 𝑟𝑜el residuo inicial 𝑟0=𝑏−𝐴𝑥0y𝑟 =𝑏−𝐴 𝑥 el residuo obtenido después del proceso de proyección con =𝐴. Entonces 𝑟 =𝑏−𝐴(𝑥0+𝛿) = 𝑟0−𝐴𝛿. (2.9) Además, 𝛿es obtenido imponiendo la condición de que 𝑟0−𝐴𝛿 sea ortogonal a 𝐴. Por lo tanto, el vector 𝐴𝛿 es la proyección ortogonal del vector 𝑟0sobre el subespacio 𝐴. Proposición 2.4. Sea 𝑥 la solución aproximada obtenida de un proceso de proyección sobre ortogonal a =𝐴y sea 𝑟 =𝑏−𝐴 𝑥 el residuo asociado. Entonces 𝑟 = (𝐼−P)𝑟0,(2.10) donde Pdenota la proyección ortogonal sobre el subespacio 𝐴. Un resultado de la proposición 2.4 es que la norma-2 del vector residual obtenido después de un paso de proyección no excederá la norma-2 inicial del residuo, es decir, ‖𝑟‖≤‖𝑟0‖2. 22 trabajo fin de grado = (𝑟𝑘, 𝑟𝑘) − 𝛼𝑘(𝐴𝑟𝑘, 𝑟𝑘) =‖𝑟𝑘‖2 2(1 − (𝐴𝑟𝑘, 𝑟𝑘) (𝑟𝑘, 𝑟𝑘) (𝐴𝑟𝑘, 𝑟𝑘) (𝐴𝑟𝑘, 𝐴𝑟𝑘))(2.18) =‖𝑟𝑘‖2 2(1 − (𝐴𝑟𝑘, 𝑟𝑘)2 (𝑟𝑘, 𝑟𝑘)2‖𝑟𝑘‖2 2 ‖𝐴𝑟𝑘‖2 2). Aplicando la desigualdad de Kantorovich (2.13) se puede afirmar que (𝐴𝑥, 𝑥) (𝑥, 𝑥)≥𝜇 > 0(2.19) donde 𝜇=𝜆𝑚𝑖𝑛(𝐴+𝐴𝑇)∕2. El resultado deseado se obtiene usando la inecuación ‖𝐴𝑟𝑘‖2≤‖𝐴‖2‖𝑟𝑘‖2. 2.3.3 Norma residual del descenso más rápido En el algoritmo de la norma residual del descenso más rápido, la suposición de que 𝐴es definida positiva se relaja. De hecho, el único requisito es que 𝐴sea una matriz no singular. En cada paso, el algoritmo usa 𝑣=𝐴𝑇𝑟y𝑤=𝐴𝑣, dando la siguiente secuencia de operaciones: 𝑟←𝑏−𝐴𝑥, 𝑣 =𝐴𝑇𝑟 𝛼←‖𝑣‖2 2∕‖𝐴𝑣‖2 2(2.20) 𝑥←𝑥+𝛼𝑣 Sin embargo, un algoritmo basado en la secuencia anterior requeriría tres productos vector por matriz, que es tres veces más que los otros algoritmos vistos en esta sección. El número de operaciones de vector por matriz puede reducirse de manera diferente. Una manera está representada por esta variante. ALGORITMO 2.4. Norma residual del descenso más rápido. 1. Calcular 𝑟=𝑏−𝐴𝑥 2. 𝑣∶= 𝐴𝑇𝑟 3. Calcular 𝐴𝑣 y𝛼∶= ‖𝑣‖2 2∕‖𝐴𝑣‖2 2 4. 𝑥∶= 𝑥+𝛼𝑣 5. 𝑟∶= 𝑟−𝛼𝐴𝑣 En esta variante, cada paso minimiza 𝑓(𝑥) = ‖𝑏−𝐴𝑥‖2 2en la dirección −▿𝑓. Como resultado, esto es equivalente al algoritmo del descenso más rápido de la sección 2.3.1 aplicado a las ecuaciones normales 𝐴𝑇𝐴𝑥 =𝐴𝑇𝑏. Como 𝐴𝑇𝐴es definida positiva cuando 𝐴es no singular, entonces, acorde con el Teorema 2.9, el método convergerá siempre que 𝐴sea no singular. Capítulo 3 Métodos iterativos de Krylov En este capítulo desarrollamos métodos que son considerados actualmente como los métodos iterativos más importantes disponibles para resolver grandes sistemas lineales. Estos métodos se basan en procesos de proyección, ortogonales u oblicuos, sobre subespacios de Krylov, que son subespacios con vectores de la forma 𝑝(𝐴)𝑣, donde 𝑝es un polinomio. Este capítulo describe métodos relacionados con la ortogonalización de Arnoldi y métodos basados en la biortogonalización de Lanczos. 3.1 Subespacio de Krylov. Recordemos del capítulo anterior que un método de proyección general para resolver el sistema lineal 𝐴𝑥 =𝑏, (3.1) extrae una solución aproximada 𝑥𝑚de un subespacio afín 𝑥0+𝑚de dimensión 𝑚al imponer la condición de Petrov-Galerkin 𝑏−𝐴𝑥𝑚⊥𝑚, donde 𝑚es otro subespacio de dimensión 𝑚. Aquí, 𝑥0representa una estimación inicial arbitraria de la solución. En esta sección consideramos métodos de proyección sobre subespacios de Krylov [10], es decir, subespacios de la forma 𝑚=𝑠𝑝𝑎𝑛{𝑣, 𝐴𝑣, 𝐴2𝑣, ..., 𝐴𝑚−1𝑣},(3.2) que denotaremos por 𝑚. La dimensión del subespacio de aproximaciones aumenta en cada paso del proceso de aproximaciones. Algunas propiedades del subespacio de Krylov se pueden establecer. Una primera propiedad es que 𝑚es el subespacio de todos los vectores en ℝ𝑛que pueden escribirse como 𝑥=𝑝(𝐴)𝑣, donde 𝑝es un polinomio de grado no superior a m-1. Recordar que el polinomio mínimo de un vector 𝑣es el polinomio mónico distinto de cero 𝑝de menor grado tal que 𝑝(𝐴)𝑣= 0. Hemos mencionado que la dimensión de 𝑚es creciente. En efecto, la siguiente proposición determina la dimensión de 𝑚en general. 24 trabajo fin de grado Proposición 3.1. El subespacio de Krylov 𝑚es de dimensión 𝑚si y sólo si el grado 𝜇de 𝑣con respecto a 𝐴no es menor que 𝑚, es decir, 𝑑𝑖𝑚(𝑚) = 𝑚⇔𝑔𝑟𝑎𝑑𝑜(𝑣)≥𝑚. (3.3) Por lo tanto, 𝑑𝑖𝑚(𝑚) = 𝑚𝑖𝑛{𝑚, 𝑔𝑟𝑎𝑑𝑜(𝑣)}.(3.4) Demostración. Los vectores 𝑣, 𝐴𝑣, ..., 𝐴𝑚−1𝑣forman una base de 𝑚si y sólo si para cualquier conjunto de m escalares 𝛼𝑖, i = 0,...,m-1, donde al menos un 𝛼𝑖es distinto de cero, la combinación lineal 𝑚−1 ∑ 𝑖=0 𝛼𝑖𝐴𝑖𝑣es distinta de cero. Esto es equivalente a la condición de que el único polinomio de grado no superior a m-1 para el cual 𝑝(𝐴)𝑣= 0 es el polinomio cero. Dado un cierto subespacio 𝑋, recordemos que 𝐴|𝑋denota la restricción de 𝐴a𝑋. Si 𝑄es un proyector para 𝑋, la sección del operador 𝐴en 𝑋es el operador para 𝑋sobre sí mismo definido por 𝑄𝐴|𝑋. La siguiente proposición caracteriza el producto de polinomios de 𝐴con 𝑣términos de la sección de 𝐴en 𝑚. Proposición 3.2. Sea 𝑄𝑚algún proyector sobre 𝑚y sea 𝐴𝑚la sección de 𝐴a𝐾𝑚; que es, 𝐴𝑚= 𝑄𝑚𝐴|𝑚. Luego, para algún polinomio q de grado no superior a m-1, 𝑞(𝐴)𝑣=𝑞(𝐴𝑚)𝑣, y, para algún polinomio de grado no superior a m, 𝑄𝑚𝑞(𝐴)𝑣=𝑞(𝐴𝑚)𝑣. Demostración. Primero probamos que 𝑞(𝐴)𝑣=𝑞(𝐴𝑚)𝑣para algún polinomio q de grado no superior a m-1. Esto es suficiente para probar la propiedad para los polinomios mónicos 𝑞𝑖(𝑡)≡𝑡𝑖, i = 0,...,m-1. La demostración es por inducción. La propiedad es cierta para el polinomio 𝑞0(𝑡)≡1.Asumiendo que esto es cierto para 𝑞𝑖(𝑡)≡𝑡𝑖: 𝑞𝑖(𝐴)𝑣=𝑞𝑖(𝐴𝑚)𝑣. Multiplicando la ecuación anterior por 𝐴en ambos lados tenemos 𝑞𝑖+1(𝐴)𝑣=𝐴𝑞𝑖(𝐴𝑚)𝑣. Si 𝑖+ 1 ≤𝑚− 1 el vector en el lado izquierdo pertenece a 𝑚y, por lo tanto, si la ecuación anterior es multiplicada en ambos lados por 𝑄𝑚, luego 𝑞𝑖+1(𝐴)𝑣=𝑄𝑚𝐴𝑞𝑖(𝐴𝑚)𝑣. Mirando el lado derecho observamos que 𝑞𝑖(𝐴𝑚)𝑣pertenece a 𝑚. Por tanto, 𝑞𝑖+1(𝐴)𝑣=𝑄𝑚𝐴|𝑚𝑞𝑖(𝐴𝑚)𝑣=𝑞𝑖+1(𝐴𝑚)𝑣, que prueba que la propiedad es cierta para i+1, siempre que 𝑖+ 1 ≤𝑚− 1. Para el caso, 𝑖+ 1 = 𝑚, solo queda demostrar que 𝑄𝑚𝑞𝑚(𝐴)𝑣=𝑞𝑚(𝐴𝑚)𝑣, que se sigue de 𝑞𝑚−1(𝐴)𝑣=𝑞𝑚−1(𝐴𝑚)𝑣simplemente multiplicando ambos lados por 𝑄𝑚𝐴. 3. métodos iterativos de krylov 25 3.2 El método de Arnoldi El método de Arnoldi es un método de proyección ortogonal sobre 𝑚para matrices no hermitianas generales. El procedimiento de Arnoldi es un algoritmo para construir una base ortogonal del subespacio de Krylov 𝑚. En la aritmética exacta, una variante del algoritmo es la siguiente. ALGORITMO 3.1. Arnoldi 1. Elegir un vector 𝑣1tal que ‖𝑣1‖2= 1 2. Para 𝑗= 1,2, ..., 𝑚 3. Calcular ℎ𝑖𝑗 = (𝐴𝑣𝑗, 𝑣𝑖)para i = 1,2,...,j 4. Calcular 𝑤𝑗∶= 𝐴𝑣𝑗− 𝑗 ∑ 𝑖=1 ℎ𝑖𝑗𝑣𝑖 5. ℎ𝑗+1,𝑗 =‖𝑤𝑗‖2 6. Si ℎ𝑗+1,𝑗 = 0 entonces parar 7. 𝑣𝑗+1 =𝑤𝑗∕ℎ𝑗+1,𝑗 En cada paso, el algoritmo multiplica el vector de Arnoldi previo 𝑣𝑗por 𝐴y luego ortonormaliza el vector resultante 𝑤𝑗contra todos los anteriores 𝑣𝑖por un procedimiento Gram-Schmidt estándar. Esto parará si el vector 𝑤𝑗calculado en la linea 4 es nulo. Ahora demostraremos algunas propiedades simples del algoritmo. Proposición 3.3. Asumiendo que el algoritmo 3.1 no se detiene antes del m-ésimo paso. Entonces, los vectores 𝑣1, ..., 𝑣𝑚forman una base ortonormal del subespacio de Krylov 𝑚=𝑠𝑝𝑎𝑛{𝑣1, 𝐴𝑣1, ..., 𝐴𝑚−1𝑣1}. Demostración. Los vectores 𝑣𝑗, j = 1,2,...,m son ortonormales por construcción. 𝑚se obtiene del hecho de que cada vector 𝑣𝑗es de la forma 𝑞𝑗−1(𝐴)𝑣1, donde 𝑞𝑗−1 es un polinomio de grado j-1. Esto puede ser demostrado por inducción en j como sigue. El resultado es claramente cierto para j = 1, ya que 𝑣1=𝑞0(𝐴)𝑣1, con 𝑞0(𝑡)≡1. Asumir que el resultado es cierto para todos los enteros no superiores a j y considerar 𝑣𝑗+1. Tenemos ℎ𝑗+1,𝑗𝑣𝑗+1 =𝐴𝑣𝑗− 𝑗 ∑ 𝑖=1 ℎ𝑖𝑗𝑣𝑖=𝐴𝑞𝑗−1(𝐴)𝑣1− 𝑗 ∑ 𝑖=1 ℎ𝑖𝑗𝑞𝑖−1(𝐴)𝑣1,(3.5) que demuestra que 𝑣𝑗+1 puede ser expresada como 𝑞𝑗(𝐴)𝑣1, donde 𝑞𝑗es de grado j, y se concluye la prueba. 3.2.1 El método de Arnoldi para sistemas lineales Dada una suposición inicial 𝑥0al sistema lineal original 𝐴𝑥 =𝑏, consideramos ahora un método de proyección ortogonal, como se define en el capítulo anterior, que toma ==𝑚(𝐴, 𝑟0), con 𝑚(𝐴, 𝑟0) = 𝑠𝑝𝑎𝑛{𝑟0, 𝐴𝑟0, 𝐴2𝑟0, ..., 𝐴𝑚−1𝑟0},(3.6) 26 trabajo fin de grado en el cuál, 𝑟0=𝑏−𝐴𝑥0.Este método busca una solución aproximada 𝑥𝑚del subespacio afín 𝑥0+𝑚 de dimensión 𝑚imponiendo la condición de Galerkin 𝑏−𝐴𝑥𝑚⊥𝑚.(3.7) Si 𝑣1=𝑟0∕‖𝑟0‖2en el método de Arnoldi y si 𝛽=‖𝑟0‖2, entonces 𝑉𝑇 𝑚𝐴𝑉𝑚=𝐻𝑚, y 𝑉𝑇 𝑚𝑟0=𝑉𝑇 𝑚(𝛽𝑣1) = 𝛽𝑒1. Como consecuencia, la solución aproximada usando los subespacios anteriores m-dimensionales es dada por 𝑥𝑚=𝑥0+𝑉𝑚𝑦𝑚,(3.8) 𝑦𝑚=𝐻𝑚−1(𝛽𝑒1).(3.9) Un método basado en lo anterior y llamado método de ortogonalización completa es descrito a continuación. El método de Gram-Schmidt es usado en el procedimiento de Arnoldi. ALGORITMO 3.2. Ortogonalización completa 1. Calcular 𝑟0=𝑏−𝐴𝑥0,𝛽∶= ‖𝑟0‖2y𝑣1∶= 𝑟0∕𝛽 2. Definir la matriz m x m 𝐻𝑚= {ℎ𝑖𝑗}𝑖,𝑗=1,...,𝑚; Tomar 𝐻𝑚= 0 3. Para 𝑗= 1,2, ..., 𝑚 4. Calcular 𝑤𝑗∶= 𝐴𝑣𝑗 5. Para 𝑖= 1, ..., 𝑗 6. ℎ𝑖𝑗 = (𝑤𝑗, 𝑣𝑖) 7. 𝑤𝑗∶= 𝑤𝑗−ℎ𝑖𝑗𝑣𝑖 8. Calcular ℎ𝑗+1,𝑗 =‖𝑤𝑗‖2. Si ℎ𝑗+1,𝑗 = 0, tomar 𝑚∶= 𝑗e ir a 12 9. Calcular 𝑣𝑗+1 =𝑤𝑗∕ℎ𝑗+1,𝑗 10. Calcular 𝑦𝑚=𝐻−1 𝑚(𝛽𝑒1)y𝑥𝑚=𝑥0+𝑉𝑚𝑦𝑚. El algoritmo anterior depende de un parámetro 𝑚que es la dimensión del subespacio de Krylov. 3.3 El método del residual mínimo generalizado (GMRES) El método del residual mínimo generalizado es un método de proyección basado en tomar = 𝑚y=𝐴𝑚, donde 𝑚es el m-ésimo subespacio de Krylov, con 𝑣1=𝑟0∕‖𝑟0‖2. Como hemos vis- to en el capítulo anterior, tal técnica minimiza la norma residual sobre todos los vectores en 𝑥0+𝑚. 3. métodos iterativos de krylov 27 La implementación de un algoritmo basado en esta propuesta es similar al algoritmo de ortogonalización completa. Hay dos maneras de describir el algoritmo. La primera usa la propiedad de optimización y la relación 𝐴𝑉𝑚=𝑉𝑚+1  𝐻𝑚.Cualquier vector 𝑥en 𝑥0+𝑚puede ser escrito como 𝑥=𝑥0+𝑉𝑚𝑦, (3.10) donde 𝑦es un m-vector. Definiendo 𝐽(𝑦) = ‖𝑏−𝐴𝑥‖2=‖𝑏−𝐴(𝑥0+𝑉𝑚𝑦)‖2,(3.11) obtenemos 𝑏−𝐴𝑥 =𝑏−𝐴(𝑥0+𝑉𝑚𝑦) =𝑟0−𝐴𝑉𝑚𝑦 =𝛽𝑣1−𝑉𝑚+1  𝐻𝑚𝑦 =𝑉𝑚+1(𝛽𝑒1− 𝐻𝑚𝑦). Como los vectores columna de 𝑉𝑚+1 son ortonormales, entonces 𝐽(𝑦)≡‖𝑏−𝐴(𝑥0+𝑉𝑚𝑦)‖2=‖𝛽𝑒1− 𝐻𝑚‖2.(3.12) La aproximación de GMRES es el único vector de 𝑥0+𝑚que minimiza (3.11). Por (3.10) y (3.12), esta aproximación puede ser obtenida de forma muy simple como 𝑥𝑚=𝑥0+𝑉𝑚𝑦𝑚, donde 𝑦𝑚minimiza la función 𝐽(𝑦) = ‖𝛽𝑒1− 𝐻𝑚𝑦‖2, es decir, 𝑥𝑚=𝑥0+𝑉𝑚𝑦𝑚, 𝑑𝑜𝑛𝑑𝑒 (3.13) 𝑦𝑚=𝑎𝑟𝑔𝑚𝑖𝑛𝑦‖𝛽𝑒1− 𝐻𝑚𝑦‖2(3.14) El minimizador 𝑦𝑚es fácil de calcular ya que requiere la solución de un problema de mínimos cuadrados (m+1) x m, donde m es típicamente pequeño. Esto nos devuelve el siguiente algoritmo: ALGORITMO 3.3: GMRES 1. Calcular 𝑟0=𝑏−𝐴𝑥0,𝛽∶= ‖𝑟0‖2, y 𝑣1∶= 𝑟0∕𝛽 2. Para j = 1,2,...,m 3. Calcular 𝑤𝑗∶= 𝐴𝑣𝑗 4. Para i = 1,...,j 5. ℎ𝑖𝑗 ∶= (𝑤𝑗, 𝑣𝑖) 6. 𝑤𝑗∶= 𝑤𝑗−ℎ𝑖𝑗𝑣𝑖 7. ℎ𝑗+1,𝑗 =‖𝑤𝑗‖2. Si ℎ𝑗+1,𝑗 = 0, tomar 𝑚∶= 𝑗e ir a 11 8. 𝑣𝑗+1 =𝑤𝑗∕ℎ𝑗+1,𝑗 9. Definir la matriz (m+1) x m de Hessengerg  𝐻𝑚= {ℎ𝑖𝑗}1≤𝑖≤𝑚+1,1≤𝑗≤𝑚 10. Calcular 𝑦𝑚, el minimizador de ‖𝛽𝑒1− 𝐻𝑚𝑦‖2, y 𝑥𝑚=𝑥0+𝑉𝑚𝑦𝑚 28 trabajo fin de grado 3.3.1 Problemas de implementación práctica Una clara dificultad de usar el algoritmo 3.3 es que no proporciona una solución aproximada 𝑥𝑚 explícita en cada paso. Como resultado, no es fácil determinar cuando parar. Un remedio es calcular la solución aproximada 𝑥𝑚en intervalos regulares y comprobar la convergencia mediante una prueba sobre el residuo, por ejemplo. Sin embargo, hay una solución más elegante, que está relacionada con la forma en la que el problema de mínimos cuadrados (3.14) está resuelto. Una técnica común para resolver el problema de mínimos cuadrados (3.14) es tránsformar la matriz de Hessenberg en una de forma triangular superior usando notaciones planas. Definimos las matrices de notación dadas, Ω𝑖= ⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ 1 ⋱ 1 𝑐𝑖𝑠𝑖 −𝑠𝑖𝑐𝑖 1 ⋱ 1 ⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ (3.15) donde 𝑐2 𝑖+𝑠2 𝑖= 1. Si los m pasos de la iteración GMRES son realizados, entonces estas matrices tienen dimensión (𝑚+ 1)𝑥(𝑚+ 1). Se multiplica la matriz de Hessenberg  𝐻𝑚y el correspondiente lado derecho 𝑔0≡𝛽𝑒1por una secuencia de tales matrices desde la izquierda. Los coeficientes 𝑠𝑖,𝑐𝑖son seleccionadas para eliminar ℎ𝑖+1 en cada iteración. Por lo tanto, si 𝑚= 5, tendriamos  𝐻5=⎛⎜⎜⎜⎜⎜⎜⎝ ℎ11 ℎ12 ℎ13 ℎ14 ℎ15 ℎ21 ℎ22 ℎ23 ℎ24 ℎ25 ℎ32 ℎ33 ℎ34 ℎ35 ℎ43 ℎ44 ℎ45 ℎ54 ℎ55 ℎ65 ⎞⎟⎟⎟⎟⎟⎟⎠ , 𝑔0=⎛⎜⎜⎜⎜⎜⎜⎝ 𝛽 0 0 0 0 0 ⎞⎟⎟⎟⎟⎟⎟⎠ Luego, se multiplica  𝐻5por Ω1=⎛⎜⎜⎜⎜⎜⎝ 𝑐1𝑠1 −𝑠1𝑐1 1 1 1 ⎞⎟⎟⎟⎟⎟⎠ con 𝑠1=ℎ21 √ℎ2 11 +ℎ2 21 , 𝑐1=ℎ11 √ℎ2 11 +ℎ2 21 , 3. métodos iterativos de krylov 29 para obtener la matriz y el lado derecho  𝐻(1) 5=⎛⎜⎜⎜⎜⎜⎜⎜⎝ ℎ(1) 11 ℎ(1) 12 ℎ(1) 13 ℎ(1) 14 ℎ(1) 15 ℎ(1) 22 ℎ(1) 23 ℎ(1) 24 ℎ(1) 25 ℎ32 ℎ33 ℎ34 ℎ35 ℎ43 ℎ44 ℎ45 ℎ54 ℎ55 ℎ65 ⎞⎟⎟⎟⎟⎟⎟⎟⎠ , 𝑔1=⎛⎜⎜⎜⎜⎜⎜⎝ 𝑐1𝛽 −𝑠1𝛽 0 0 0 0 ⎞⎟⎟⎟⎟⎟⎟⎠ .(3.16) Ahora podemos multiplicar la matriz anterior y el lado derecho de nuevo por una matriz de notación Ω2para eliminar ℎ32. Esto se consigue tomando 𝑠2=ℎ32 √(ℎ(1) 22 )2+ℎ2 32 , 𝑐2=ℎ(1) 22 √(ℎ(1) 22 )2+ℎ2 32 . Este proceso de eliminación continúa hasta que la m-ésima notación es aplicada, que transforma el problema en uno que involucra la matriz y el lado derecho  𝐻(5) 5=⎛⎜⎜⎜⎜⎜⎜⎜⎝ ℎ(5) 11 ℎ(5) 12 ℎ(5) 13 ℎ(5) 14 ℎ(5) 15 ℎ(5) 22 ℎ(5) 23 ℎ(5) 24 ℎ(5) 25 ℎ(5) 33 ℎ(5) 34 ℎ(5) 35 ℎ(5) 44 ℎ(5) 45 ℎ(5) 55 0 ⎞⎟⎟⎟⎟⎟⎟⎟⎠ , 𝑔5=⎛⎜⎜⎜⎜⎜⎝ 𝛾1 𝛾2 𝛾3 ⋮ 𝛾6 ⎞⎟⎟⎟⎟⎟⎠ .(3.17) Generalmente, los escalares 𝑐𝑖y𝑠𝑖de la i-ésima notación Ω𝑖son definidas como 𝑠𝑖=ℎ𝑖+1,𝑖 √(ℎ(𝑖−1) 𝑖𝑖 )2+ℎ2 𝑖+1,𝑖 , 𝑐𝑖=ℎ(𝑖−1) 𝑖𝑖 √(ℎ(𝑖−1) 𝑖𝑖 )2+ℎ2 𝑖+1,𝑖 .(3.18) Definir 𝑄𝑚como el producto de matrices Ω𝑖: 𝑄𝑚= Ω𝑚Ω𝑚−1...Ω1,(3.19) y tomar  𝑅𝑚= 𝐻(𝑚) 𝑚=𝑄𝑚 𝐻𝑚,(3.20) 𝑔𝑚=𝑄𝑚(𝛽𝑒1)=(𝛾1, ..., 𝛾𝑚+1)𝑇.(3.21) Como 𝑄𝑚es unitario 𝑚𝑖𝑛‖𝛽𝑒1− 𝐻𝑚𝑦‖2=𝑚𝑖𝑛‖𝑔𝑚− 𝑅𝑚𝑦‖2. La solución del problema de mínimos cuadrados anterior es obtenida simplemente resolviendo el sistema triangular resultante de la eliminación de la última fila de la matriz  𝑅𝑚y del lado derecho 𝑔𝑚en (3.17). Proposición 3.4. Sea 𝑚≤𝑛yΩ𝑖,𝑖= 1, ..., 𝑚, las matrices de notación usadas para transformar  𝐻𝑚en una de forma triangular superior. Denotar por  𝑅𝑚,𝑔𝑚= (𝛾1, ..., 𝛾𝑚+1)𝑇la matriz resultante y el lado derecho, definidos por (3.20) y (3.21), y por 𝑅𝑚,𝑔𝑚la matriz triangular superior 𝑚𝑥𝑚 y el vector m-dimensional obtenido de  𝑅𝑚,𝑔𝑚por eliminar su última fila y componente, respectivamente. Entonces tenemos las siguientes propiedades: 30 trabajo fin de grado 1. El rango de 𝐴𝑉𝑚es igual al rango de 𝑅𝑚. En particular, si 𝑟𝑚𝑚 = 0, entonces A debe ser singular. 2. El vector 𝑦𝑚que minimiza ‖𝛽𝑒1− 𝐻𝑚𝑦‖2es dado por 𝑦𝑚=𝑅−1 𝑚𝑔𝑚. 3. El vector residual en el paso m satisface 𝑏−𝐴𝑥𝑚=𝑉𝑚+1(𝛽𝑒1− 𝐻𝑚𝑦𝑚) = 𝑉𝑚+1𝑄𝑇 𝑚(𝛾𝑚+1𝑒𝑚+1),(3.22) y como resultado, ‖𝑏−𝐴𝑥𝑚‖2=|𝛾𝑚+1|.(3.23) Hasta ahora hemos descrito un proceso para calcular la solución de mínimos cuadrados 𝑦𝑚(3.14). Esto junto con las notaciones del plano sirve para resolver el sistema lineal (3.9) del método de ortogonalización completa. Es posible implementar el proceso anterior de una manera progresiva, es decir, en cada paso del algoritmo GMRES. Este enfoque permitirá obtener la norma residual en cada paso, sin operaciones adicionales de aritmética. Para ilustrar esto, asumir que la primera 𝑚rotación ya ha sido aplicada, es decir, comenzar en (3.17). Ahora, la norma residual es válida para 𝑥5y el criterio de parada puede ser aplicado. Un paso más en el algoritmo de Arnoldi debe ser ejecutado para obtener 𝐴𝑣6y la sexta columna de  𝐻6. Esta columna se adjunta a  𝑅5, que ha sido aumentada por una fila de ceros para que coincida la dimensión. A continuación, hemos obtenido la siguiente matriz y lado derecho  𝐻(5) 6= ⎛⎜⎜⎜⎜⎜⎜⎜⎜⎝ ℎ11 ℎ12 ℎ13 ℎ14 ℎ15 ℎ16 ℎ22 ℎ23 ℎ24 ℎ25 ℎ26 ℎ33 ℎ34 ℎ35 ℎ36 ℎ44 ℎ45 ℎ46 ℎ55 ℎ56 0ℎ66 0ℎ76 ⎞⎟⎟⎟⎟⎟⎟⎟⎟⎠ , 𝑔(5) 6=⎛⎜⎜⎜⎜⎜⎜⎝ 𝛾1 𝛾2 𝛾3 ⋮ 𝛾6 0 ⎞⎟⎟⎟⎟⎟⎟⎠ .(3.24) El algoritmo ahora continúa de la misma forma que antes. Necesitamos multiplicar la matriz por una matriz de notación Ω6(ahora de dimensión 7 x 7), con 𝑠6=ℎ76 √(ℎ66)2+ℎ2 76 , 𝑐6=ℎ(5) 66 √(ℎ(5) 66 )2+ℎ2 76 ,(3.25) para obtener la matriz y lado derecho  𝑅6= ⎛⎜⎜⎜⎜⎜⎜⎜⎜⎝ 𝑟11 𝑟12 𝑟13 𝑟14 𝑟15 𝑟16 𝑟22 𝑟23 𝑟24 𝑟25 𝑟26 𝑟33 𝑟34 𝑟35 𝑟36 𝑟44 𝑟45 𝑟46 𝑟55 𝑟56 𝑟66 0 ⎞⎟⎟⎟⎟⎟⎟⎟⎟⎠ , 𝑔6=⎛⎜⎜⎜⎜⎜⎜⎝ 𝛾1 𝛾2 𝛾3 ⋮ 𝑐6𝛾6 −𝑠6𝛾6 ⎞⎟⎟⎟⎟⎟⎟⎠ .(3.26) 3. métodos iterativos de krylov 31 Si la norma residual dada por |𝛾𝑚+1|es suficientemente pequeña, el proceso debe ser parado. Las últimas filas de  𝑅𝑚y𝑔𝑚son eliminadas y el sistema triangular superior resultante es resuelto para obtener 𝑦𝑚. Entonces, la solución aproximada 𝑥𝑚=𝑥0+𝑉𝑚𝑦𝑚es calculada. Anotar para (3.26) que la siguiente relación conveniente para 𝛾𝑗+1 sea: 𝛾𝑗+1 = −𝑠𝑗𝛾𝑗.(3.27) En particular, si 𝑠𝑗= 0, entonces la norma residual debe ser igual a cero, que significa que la solución es exacta en el paso j. 3.4 Algoritmo simétrico de Lanczos Este algoritmo se puede ver como una simplificación del método de Arnoldi para el caso particular cuando la matriz es simétrica. Cuando 𝐴es simétrica, entonces, la matriz de Hessenberg 𝐻𝑚 se convierte en una tridiagonal simétrica. Esto lleva a tres términos recurrentes en el proceso de Arnoldi y recurrencia a corto plazo para soluciones algorítmicas como la ortogonalización completa y GMRES. 3.4.1 El algoritmo Para introducir el algoritmo de Lanczos comenzamos haciendo la observación escrita en el siguiente teorema. Teorema 3.5. Asumimos que el método de Arnoldi es aplicado a una matriz real simétrica 𝐴. Entonces, los coeficientes ℎ𝑖𝑗 generados por el algoritmo son tales que ℎ𝑖𝑗 = 0 𝑝𝑎𝑟𝑎 1≤𝑖<𝑗− 1 (3.28) ℎ𝑗,𝑗+1 =ℎ𝑗+1,𝑗, 𝑗 = 1, ..., 𝑚. (3.29) En otras palabras, la matriz 𝐻𝑚obtenida del proceso de Arnoldi es tridiagonal y simétrica. Demostración. La prueba es una consecuencia inmediata de que 𝐻𝑚=𝑉𝑇 𝑚𝐴𝑉𝑚sea una matriz simétrica que es también una matriz de Hessenberg por construcción. Por lo tanto, debe ser una matriz tridiagonal simétrica. La notación estándar usada para describir el algoritmo de Lanczos es obtenida imponiendo 𝛼𝑗≡ℎ𝑖𝑗, 𝛽𝑗≡ℎ𝑗−1,𝑗 y, si 𝑇𝑚denota la matriz 𝐻𝑚resultante, esta es de la forma 𝑇𝑚=⎛⎜⎜⎜⎜⎜⎝ 𝛼1𝛽2 𝛽2𝛼2𝛽3 ⋯ 𝛽𝑚−1 𝛼𝑚−1 𝛽𝑚 𝛽𝑚𝛼𝑚 ⎞⎟⎟⎟⎟⎟⎠ .(3.30) 38 trabajo fin de grado 𝑢′(𝑥𝑗) = 𝑢(𝑥𝑗+1) − 𝑢(𝑥𝑗−1) 2ℎ−1 6𝑢′′′(𝜁1)⋅ℎ2, 𝜁1∈ (𝑥𝑗−1, 𝑥𝑗+1),(4.3) 𝑢′′(𝑥𝑗) = 𝑢(𝑥𝑗+1)−2𝑢(𝑥𝑗) + 𝑢(𝑥𝑗−1) ℎ2−1 12𝑢(4)(𝜁2)⋅ℎ2, 𝜁2∈ (𝑥𝑗−1, 𝑥𝑗+1),(4.4) Usando (4.3) y (4.4) llegamos a: {𝑢0=𝛼, 𝑢𝑁+1 =𝛽 𝑢𝑗−1−2𝑢𝑗+𝑢𝑗+1 ℎ2=𝑎(𝑥𝑗)𝑢𝑗+1−𝑢𝑗−1 2ℎ+𝑟(𝑥𝑗)𝑢𝑗+𝑓(𝑥𝑗), 𝑗 = 1, ..., 𝑁. (4.5) En el siguientes ejemplo consideramos la ecuación de Poisson para n = 2,3. Para ello vamos a tomar 𝑑= 1,𝑎=𝑟= 0, con lo que obtenemos, − 𝑛 ∑ 𝑖=1 𝑢𝑥𝑖𝑥𝑖=𝑓(𝑥).(4.6) 𝑛= 2 Tenemos el problema {−𝑢𝑥1𝑥1−𝑢𝑥2𝑥2=𝑓, (𝑥1, 𝑥2) ∈ Ω = (0,1)2 𝑢=𝑔, (𝑥1, 𝑥2) ∈ Γ = 𝜕Ω(4.7) y si consideramos la malla Ω(2) ℎ= {(𝑖ℎ, 𝑗ℎ)}𝑁+1 𝑖,𝑗=0, con ℎ=1 𝑁+1 , y diferencias centrales de orden dos: 𝑢𝑥1𝑥1(𝑥1, 𝑥2) ≃ 𝑢(𝑥1+ℎ, 𝑥2)−2𝑢(𝑥1, 𝑥2) + 𝑢(𝑥1−ℎ, 𝑥2) ℎ2, 𝑢𝑥2𝑥2(𝑥1, 𝑥2) ≃ 𝑢(𝑥1, 𝑥2+ℎ)−2𝑢(𝑥1, 𝑥2) + 𝑢(𝑥1, 𝑥2−ℎ) ℎ2, se obtiene el esquema en diferencias finitas: (−𝑈𝑖+1,𝑗 + 2𝑈𝑖,𝑗 −𝑈𝑖−1,𝑗) ℎ2+(−𝑈𝑖,𝑗+1 + 2𝑈𝑖,𝑗 −𝑈𝑖,𝑗−1) ℎ2=𝑓𝑖,𝑗 1≤𝑖, 𝑗 ≤𝑁, (4.8) que corresponde a un sistema lineal de dimensión 𝑁2 𝐴(2) ⋅𝑈=𝑏, (4.9) donde 𝑏∈ℝ𝑁2contiene la discretización del dato 𝑓(𝑥1, 𝑥2)y las condiciones de frontera 𝑔(𝑥1, 𝑥2), 𝑈= (𝑈11, ..., 𝑈𝑁1, 𝑈12, ..., 𝑈𝑁2, ..., 𝑈1𝑁, ..., 𝑈𝑁𝑁 )𝑇∈ℝ𝑁2,𝑈𝑖𝑗 ≃𝑈(𝑖ℎ, 𝑗ℎ)y𝐴(2) es la matriz de dimensión 𝑁2𝑥𝑁2, 4. aplicación a la resolución numérica de edp y en matrices test 39 A(2) = ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ 4 −1 −1 −1 ⋱ ⋱ ⋱ ⋱ ⋱ −1 ⋱ −1 4 −1 −1 4 −1 ⋱ ⋱−1 ⋱ ⋱ ⋱ ⋱ ⋱ ⋱ −1 ⋱ −1 −1 4 ⋱ ⋱ ⋱ −1 ⋱ ⋱ ⋱ ⋱ ⋱ ⋱ ⋱ ⋱ −1 −1 4 −1 ⋱−1 ⋱ ⋱ ⋱ ⋱ ⋱ −1 −1 −1 4 ⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ (4.10) 𝑛= 3 Tenemos el problema {−𝑢𝑥1𝑥1−𝑢𝑥2𝑥2−𝑢𝑥3𝑥3=𝑓, (𝑥1, 𝑥2, 𝑥3) ∈ Ω = (0,1)3, 𝑢=𝑔, (𝑥1, 𝑥2, 𝑥3) ∈ Γ = 𝜕Ω,(4.11) 40 trabajo fin de grado y consideramos la malla Ω(3) ℎ= {(𝑖ℎ, 𝑗ℎ, 𝑘ℎ)}𝑁+1 𝑖,𝑗,𝑘=1, con ℎ=1 𝑁+1 , se obtiene el esquema en diferencias finitas: (−𝑈𝑖+1,𝑗,𝑘+2𝑈𝑖,𝑗,𝑘−𝑈𝑖−1,𝑗,𝑘)+(−𝑈𝑖,𝑗+1,𝑘+2𝑈𝑖,𝑗,𝑘−𝑈𝑖,𝑗−1,𝑘)+(−𝑈𝑖,𝑗,𝑘+1+2𝑈𝑖,𝑗,𝑘−𝑈𝑖,𝑗,𝑘−1) = 𝑓𝑖,𝑗,𝑘,1≤𝑖, 𝑗, 𝑘 ≤𝑁, (4.12) que corresponde a un sistema lineal de dimensión 𝑁3 𝐴(3) ⋅𝑈=𝑏, (4.13) donde 𝑏∈ℝ𝑁3contiene la discretización del dato 𝑓(𝑥1, 𝑥2, 𝑥3)y las condiciones de frontera 𝑔(𝑥1, 𝑥2, 𝑥3), 𝑈= (𝑈𝑖.𝑗,𝑘)𝑁 𝑖,𝑗,𝑘=1 ∈ℝ𝑁3,𝑈𝑖𝑗 ≃𝑈(𝑖ℎ, 𝑗ℎ, 𝑘ℎ)y𝐴(3) es la matriz de dimensión 𝑁3𝑥𝑁3, 𝐴(3) = ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ 𝐴(2) −𝐼0⋯0 0 −𝐼 𝐴(2) −𝐼⋯0 0 ⋱ ⋱ ⋱ ⋱⋱⋱ ⋱ ⋱ −𝐼 0 0 −𝐼 𝐴(2) ⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ (4.14) A continuación, describimos un ejemplo en dimensiones 𝑛= 2 y𝑛= 3 donde obtenemos que la solución de (4.6) es un polinomio de grado menor o igual que dos. El vector 𝑈se usará para medir el error de las iteraciones provistas por los métodos vistos en los capítulos anteriores. La función 𝑓 de la ecuación de Poisson la obtenemos a partir de la solución exacta que tomemos y el vector 𝑏se obtiene discretizando dicha 𝑓junto con las condiciones de frontera. 4.1.1 Ejemplo 1 En este ejemplo consideramos la ecuación de Poisson con condiciones de frontera de tipo Dirichlet homogéneas en [0,1]𝑛, n = 2,3, con solución exacta 𝑢=𝑢(𝑥1, ..., 𝑥𝑛)=4𝑛⋅ 𝑛 ∏ 𝑖=1 𝑥𝑖(1 − 𝑥𝑖). Para las dimensiones n = 2 y n = 3, la discretización de (4.2) mediante diferencias centrales de orden 2 produce los sistemas lineales (4.9) y (4.13) 𝐴(𝑛)⋅𝑈=𝑏, de dimensión 𝑁𝑛, con 𝐴(𝑛)simétrica y definida positiva. Tomamos en cada caso N = 100. 4. aplicación a la resolución numérica de edp y en matrices test 41 •Para n = 2, la solución exacta es 𝑢(𝑥1, 𝑥2) = 42⋅𝑥1(1 − 𝑥1)⋅𝑥2(1 − 𝑥2), de donde obtenemos 𝑓(𝑥1, 𝑥2) = 32⋅[𝑥1(1−𝑥1)+𝑥2(1−𝑥2)]. Y obtenemos la matriz 𝑏discretizando la 𝑓. A continuación, mostramos el comportamiento de la norma del error en relación al número de iteraciones. Hemos comparado los métodos de GMRES (algoritmo 3.3), del Gradiente conjugado (algoritmo 3.7) y del Residual mínimo (algoritmo 2.3). Figura 4.1: Velocidad de convergencia en los métodos GMRES, GC y MINRES aplicados a la ecuación de Poisson 2D con solución exacta 𝑢(𝑥1, 𝑥2)=42⋅𝑥1(1 − 𝑥1)⋅𝑥2(1 − 𝑥2). Observamos en la Figura 4.3 que la convergencia de GC y MINRES acaban muy cerca en apenas 150 iteraciones, mientras que GMRES necesita más iteraciones para converger. •Para n = 3, la solución exacta es 𝑢(𝑥1, 𝑥2, 𝑥3) = 43⋅𝑥1(1 − 𝑥1)⋅𝑥2(1 − 𝑥2)⋅𝑥3(1 − 𝑥3), de donde obtenemos la 𝑓(𝑥1, 𝑥2, 𝑥3) = 192 ⋅[𝑥1(1 − 𝑥1) + 𝑥2(1 − 𝑥2) + 𝑥3(1 − 𝑥3)]. Y obtenemos la matriz 𝑏 discretizando el dato f. A continuación, mostramos el comportamiento del error para los métodos GC y GMRES, y sus variantes con reinicialización GMRES(m) tras m iteraciones, con m = 10,50,100. Observamos nuevamente que el decaimiento del error con las iteraciones es similar para GC y GMRES con m = 50,100. Mientras que para GMRES con m = 10 la convergencia necesita el máximo de iteraciones. 42 trabajo fin de grado Figura 4.2: Error frente al número de iteraciones de los métodos GC y GMRES(m), m = 10,50,100, aplicados a la ecuación de Poisson 3D. 4.2 Aplicación en matrices test En esta sección vamos a centrarnos en el método GMRES precondicionado y vamos a observar la convergencia frente al número de iteraciones en dos ejemplos, uno con reinicialización y otro sin reinicialización usando problemas test con matrices grandes y dispersas [4]. 4.2.1 Precondicionamiento El precondicionamiento se basa en intentar encontrar una matriz 𝑀tal que 𝑀−1𝐴tenga mejores propiedades espectrales de forma que se tenga convergencia rápida cuando se aplica al sistema modificado 𝑀−1𝐴=𝑀−1𝑏. Para el caso 𝑀=𝐴, todos los métodos alcanzarían la solución exacta en un único paso. Por tanto, buscamos una aproximación 𝑀de 𝐴de manera que el correspondiente método de Krylov aplicado a 𝑀−1𝐴=𝑀−1𝑏solo necesite una pocas iteraciones para alcanzar una buena aproximación de la solución del sistema. La matriz 𝑀en este contexto recibe el nombre de precondicionador para la matriz 𝐴. El mismo precondicionador puede ser implementado de diferentes formas. Esto no cambia los autovalores de la matriz precondicionada pero sí los autovectores, lo cual puede repercutir en el comportamiento de la convergencia. Existen tres implementaciones diferentes: 1. Precondicionamiento por la izquierda. Aplicamos el método iterativo a 𝑀−1𝐴=𝑀−1𝑏. Notemos que la simetría de 𝐴y𝑀no implica la simetría de 𝑀−1𝐴. Sin embargo, si 𝑀es simétrica definida positiva entonces [𝑥, 𝑦]=(𝑥, 𝑀𝑦)define un producto interno adecuado. 4. aplicación a la resolución numérica de edp y en matrices test 43 Es fácil verificar que 𝑀−1𝐴es simétrica respecto a este nuevo producto interno, y por tanto, podemos usar métodos como MINRES y GC. Si usamos un método que minimice la norma del residuo (como GMRES o MINRES) debemos tener en cuenta que estamos minimizando el residuo precondicionado 𝑀−1(𝑏−𝐴𝑥𝑘), que podría ser bastante diferente del residuo 𝑏−𝐴𝑥𝑘. 2. Precondicionamiento por la derecha. Aplicamos el método iterativo a 𝐴𝑀−1𝑦=𝑏 𝑐𝑜𝑛 𝑥 =𝑀−1𝑦. Esta forma tampoco deja un producto simétrico aunque 𝐴y𝑀sean simétricas. Este precondicionamiento tiene la ventaja de que sólo afecta al operador y no al lado derecho del sistema. 3. Precondicionamiento por ambos lados. Para un precondicionador 𝑀con 𝑀=𝑀1𝑀2, el método iterativo puede ser aplicado a 𝑀−1 1𝐴𝑀−1 2𝑧=𝑀−1 1𝑏 𝑐𝑜𝑛 𝑥 =𝑀−1 2𝑥. Esta forma podría ser útil para aquellos precondicionadores que vienen en forma factorizada o para obtener un operador simétrico cuando 𝑀no pueda ser usada para la definición de producto interno (es decir, 𝑀no es simétrica definida positiva). Observación 4.1. La elección del precondicionador 𝑀puede influir bastante en la eficiencia del correspondiente método precondicionado y, por lo general, 𝑀ha de elegirse adaptado al sistema lineal que se desea resolver. No obstante, en este trabajo consideraremos un tipo de precondicionador general basado en las denominadas factorizaciones matriciales incompletas. Consideramos el sistema lineal disperso 𝐴𝑥 =𝑏, con 𝐴∈ℝ𝑁𝑥𝑁 una matriz 𝑛𝑥𝑛. Con frecuencia este es resuelto utilizando eliminación Gaussiana que es equivalente a factorizar la matriz 𝐴como 𝐴=𝐿𝑈. Con lo cuál, elegimos 𝑀=𝐿⋅𝑈, donde 𝐿y𝑈forman una descomposición 𝐿𝑈 de la matriz 𝐴, con 𝐿 triangular inferior y 𝑈triangular superior. Sin embargo, tal descomposición puede ser muy costosa si la dimensión de la matriz 𝐴es suficientemente grande. Una primera estrategia para relajar la exigencia computacional de la descomposición 𝐿𝑈 consiste en fijar un subconjunto de índices 𝑆 ⊆ {1, .., 𝑁}𝑥{1, .., 𝑁}, y seguir los pasos usuales de los algoritmos para computar una descomposición 𝐿𝑈 pero sólo en aquellos índices (𝑖, 𝑗) ∈ 𝑆. En particular, una usual opción a considerar es tomar 𝑆= {(𝑖, 𝑗) ∈ {1, .., 𝑁}2|𝑎𝑖,𝑗 ≠0},esto es, realizar los pasos de la descomposición 𝐿𝑈 pero respetando las componentes nulas de la matriz 𝐴original. Esta elección da lugar a la factorización incompleta 𝐿𝑈 con "nivel 0 de relleno" denotado como 𝐼𝐿𝑈(0). Observar que si 𝐴es una matriz sparse, es decir, con una elevada proporción de elementos nulos, el cálculo de factorizaciones 𝐿𝑈 hace que los factores 𝐿y𝑈tengan un mayor número de componentes no nulas, lo cual conlleva un mayor gasto computacional en la resolución de sistemas 𝐿⋅𝑈⋅𝑦=𝐶por sustitución hacia atrás y hacia adelante. Por ello, la elección de un precondicionador 𝑀=𝐿⋅𝑈a través de una factorización incompleta suele ser ventajosa si la dimensión N de 𝐴es elevada, a pesar de que 𝐴−𝐿𝑈 ≠0. Consideramos en particular el siguiente pseudocódigo para la factorización 𝐼𝐿𝑈(0) [2]. 44 trabajo fin de grado ALGORITMO 4.1: ILU(0) (por el método de Doolittle, esto es, 𝐿con unos en la diagonal). 1. Para 𝑟=1∶𝑁− 1 2. Para 𝑖=𝑟+ 1 ∶ 𝑁 3. Si (𝑖, 𝑟) ∈ 𝑆, entonces 𝑎𝑖,𝑟 =𝑎𝑖,𝑟⟋𝑎𝑟,𝑟. 4. Para 𝑗=𝑟+ 1 ∶ 𝑁 5. Si (𝑖, 𝑗) ∈ 𝑆y(𝑟, 𝑗) ∈ 𝑆, entonces 𝑎𝑖,𝑗 =𝑎𝑖,𝑗 −𝑎𝑖,𝑟 ⋅𝑎𝑟,𝑗 𝑎𝑟,𝑟 .(4.15) Observación 4.2: (Factorizaciones incompletas modificadas). En la factorización 𝐼𝐿𝑈(0), si (𝑖, 𝑗) ∉ 𝑆, entonces el valor 𝑎𝑖,𝑟 ⋅𝑎𝑟,𝑗∕𝑎𝑟,𝑟 no se sustrae de 𝑎𝑖,𝑗, como correspondería a una factorización 𝐿𝑈 completa. Una opción alternativa que suele ser ventajosa en la práctica consiste en no despreciar completamente el valor 𝑎𝑖,𝑟 ⋅𝑎𝑟,𝑗∕𝑎𝑟,𝑟 cuando (𝑖, 𝑗) ∉ 𝑆y sustraerlo de la componente diagonal de la fila correspondiente, eso es, 𝑎𝑖,𝑖. Esto da lugar a la descomposición imcompleta modificada 𝐿𝑈,𝑀𝐼𝐿𝑈(0). ALGORITMO 4.2: MILU(0) 1. Si (𝑖, 𝑗) ∈ 𝑆, entonces 2. Si (𝑖, 𝑗) ∈ 𝑆, entonces 𝑎𝑖,𝑗 =𝑎𝑖,𝑗 −𝑎𝑖,𝑟⋅𝑎𝑟,𝑗 𝑎𝑟,𝑟 ; en caso contrario, 𝑎𝑖,𝑖 =𝑎𝑖,𝑖 −𝑎𝑖,𝑟 ⋅𝑎𝑟,𝑗 𝑎𝑟,𝑟 (4.16) A continuación, nos centraremos en el método GMRES precondicionado. Para ello, primero desarrollamos su algoritmo: ALGORITMO 4.1: MÉTODO GMRES PRECONDICIONADO (PGMRES): Consideramos el sistema lineal 𝐴𝑥 =𝑏donde 𝐴∈ℝ𝑁𝑥𝑁 es una matriz regular, 𝑏∈ℝ𝑁un vector no nulo y 𝑀una matriz regular. Dado 𝑥0∈ℝ𝑁y𝑟0=𝑏−𝐴𝑥0, 1. Definimos 𝑟0∶= 𝑀−1𝑟0,𝑣1=𝑟0 ‖𝑟0‖2 y𝑐=‖𝑟0‖2𝑒1∈ℝ𝑁. 2. Sean 𝑣𝑛=𝑀−1𝐴𝑣𝑛yℎ𝑖𝑛 ∶= ⟨𝑣𝑛, 𝑣𝑖⟩2para 1≤𝑖≤𝑛. Calcular 𝑣𝑛+1 =𝑣𝑛− 𝑛 ∑ 𝑖=1 ℎ𝑖𝑛𝑣𝑖y ℎ𝑛+1,𝑛 =‖𝑣𝑛+1‖2. 3. Continuar con los pasos 7-10 del Algoritmo del Método GMRES. Para los siguientes ejemplos vamos a tomar una matriz dispersa no simétrica real de dimensión 479 por 479 predefinida en MATLAB [4]. Esta matriz puede surgir de la discretización de ciertas EDP (e.g., ecuación de convección-difusión-reacción). 4. aplicación a la resolución numérica de edp y en matrices test 45 4.2.2 Ejemplo 2 Vamos a observar el efecto de usar una matriz precondicionada y no precondicionada sin reinicialización en el método GMRES para resolver un sistema lineal. Para el caso precondicionado, como 𝐴es no simétrica, usamos el método de precondicionamiento ILU descrito anteriormente para generar la matriz precondicionada 𝑀=𝐿𝑈. A continuación, resolvemos el sistema precondicionado 𝑀−1𝐴𝑥 =𝑀−1𝑏donde 𝐿y𝑈son argumentos de entrada de GMRES. Figura 4.3: Error frente al número de iteraciones del método GMRES sin reinicialización para los casos precondicionado y no precondicionado. Para el método no precondicionado se ha realizado el número máximo de iteraciones sin conseguirse la convergencia. Mientras que el método precondicionado ha convergido con la tolerancia deseada en un número menor de iteraciones que las máximas establecidas. El uso del método de precondicionamiento ILU produce un error menor que la tolerancia descrita en la sexta iteración. 4.2.3 Ejemplo 3 Veamos ahora el efecto de usar una matriz precondicionada con distintas reinicialización en el método GMRES(m), m = 3,4,5. La siguiente gráfica muestra el historial de convergencia de cada método GMRES reiniciado. Observamos que gmres(3) converge en la iteración 5 por encima de la tolerancia, por otro lado, para el caso gmres(4) y gmres(5), la convergencia se aproxima a la iteración 3. Podemos observar que a mayor reinicialización converge más rápido en menos iteraciones. 46 trabajo fin de grado Figura 4.4: Error frente al número de iteraciones del método GMRES con reinicialización para m = 3,4 y 5. 4.3 Conclusiones El comportamiento de los métodos de proyección presentados en este trabajo queda recogido en los ejemplos pertenecientes al último capítulo. Se intuye que el método del Gradiente Conjugado y del Residual mínimo resultan más eficientes que el método GMRES para una ecuación de Poisson con condiciones de frontera tipo Dirichlet homogéneas, para una matriz simétrica definida positiva, tomando n = 2. Para el caso en el que n = 3, tomando la misma ecuación y condiciones que para n = 2 y con matriz simétrica, concluimos que el método GMRES con reinicialización mayor nos lleva a una convergencia muy parecida a la obtenida por el método del gradiente conjugado. A continuación, para el caso de una matriz test no simétrica que puede surgir, por ejemplo, de la discretización de una EDP de la forma convección-difusión-reacción, estudiamos el error frente al número de iteraciones para el método GMRES precondicionado. En el ejemplo 2, en el cuál vemos la diferencia entre dicho método precondicionado y no precondicionado sin reinicialización, obtenemos que el método no precondicionado no converge mientras que el precondicionado si lo hace y además en pocas iteraciones. Para finalizar, en el ejemplo 3, donde vemos el comportamiento del método gmres precondicionado con reinicialización, obtenemos que mientras mayor reinicialización, menos iteraciones son necesarias para que converga dicho método. Por lo tanto, el objetivo de nuestro trabajo queda expuesto en estas gráficas como resultado del comportamiento de los diferentes métodos de proyección basados en subespacios de Krylov cuando se tratan de distintos sistemas lineales mostrando las condiciones más favorables a cada uno. Bibliografía [1] Aubanell, A., Benseny, A. Delshams, A., Útiles básicos de Cálculo Numérico. Labor (1993). [2] Axelsson, 0. Iterative solution methods. Cambridge University Press (1996). [3] Ayres, F. Ecuaciones diferenciales. Mc Graw-Hill (1991). [4] Barrett, R., Berry, M., Chan, T., Demmen, J., Donat, J., Dongarra, J., Eijkhout,V., Pozo, R., Romine, C., Van der Vorst, H. Templates for the solution of linear systems: building blocks for iterative systems. SIAM (1994) [5] Ciarlet, P.G., Introducción à l’analyse numérique matricielle et à l’optimization. Masson (1990). [6] Cullum, J., Willoughby, R. Lanczos Algoritms for Large Symmetric Eigenvalue Computations. Birkhäuser (1984). [7] Hackbusch, W. Iterative Solution of Large Sparse Systems of Equations. Springer (2016). [8] Meurant, G., Duintjer Tebbens, J. Krylov methods for nonsymmetric linear systems: from theory to computations. Springer (2020). [9] Quarteroni, A., Sacco, R., Saleri, F. Numerical Mathematics. Springer (2007). [10] Saad, Y. Iterative methods for sparse linear systems. SIAM (2003).