Precondicionamiento de sistemas de ecuaciones de matrices variables en la modelización de campos de viento
Abstract
Programa de Doctorado: Sistemas Inteligentes y Aplicaciones Numéricas en Ingeniería
Full text
Instituto Universitario de Sistemas Inteligentes y Aplicaciones Numéricas en Ingeniería Tesis Doctoral PRECONDICIONAMIENTO DE SISTEMAS DE ECUACIONES DE MATRICES VARIABLES EN LA MODELIZACIÓN DE CAMPOS DE VIENTO Héctor Sarmiento Almeida Las Palmas de Gran Canaria, Abril de 2010 UNIVERSIDAD DE LAS PALMAS DE GRAN CANARIA
A mis nietos, Tom´as y Sara
Agradecimientos A Antonio Su´arez Sarmiento y Gustavo Montero Garc´ıa, directores de esta tesis, sin cuyo aliento, asesoramiento y tutela no hubiera sido posible la realizaci´on de este trabajo. A Eduardo Rodr´ıguez Barrera que, pacientemente, me ha ayudado a resolver los imprevistos y dudas inform´aticas que fueron surgiendo durante el desarrollo del tema. A Mar´ıa Dolores Garc´ıa Le´on y Elizabeth Fl´orez V´azquez, por su inestimable colaboraci´on, gracias a las cuales, las ideas y experimentos necesarios para el desarrollo y presentaci´on de esta tesis, se han hecho realidad. A todos los compa˜neros del Departamento de Matem´aticas que de una forma u otra me han animado y ayudado, en todo momento, a la terminaci´on de esta tesis. A mi familia, en especial, por su amable paciencia durante los largos tiempos de ausencia que ha supuesto para ellos mi dedicaci´on a este tema. A todos, muchas gracias.
Resumen En la formulaci´on matem´atica de los modelos de campos de viento surgen grandes sistemas de ecuaciones lineales, caracterizados por tener matrices variables, al depender estas de un cierto par´ametro, , tal que: Ax=bdonde Aes una matriz sim´etrica del tipo A=M+ N, siendo M y N dos matrices, tipo sparse, diferentes, Sim´etricas Definidas Positivas (SDP) y el par´ametro ≥0. Los m´etodos basados en los subespacios de Krylov constituyen la mejor alternativa para la resoluci´on de los sistemas de ecuaciones sparse. En el caso particular de sistemas cuya matriz es SDP, el m´etodo del Gradiente Conjugado, es el que presenta los mejores resultados En esta tesis se trata de extender el uso del algoritmo del Gradiente Conjugado Precondicionado a los sistemas de ecuaciones de matrices variables, estudiando los Precondicionadores m´as adecuados para mejorar su convergencia. El presente trabajo est´a estructurado en tres partes.En una primera parte se describen los distintos tipos de modelizaciones de campos de viento y el proceso de generaci´on de sus sistemas de ecuaciones lineales de matrices variables. En la segunda parte se presenta el estado del arte de los m´etodos iterativos para la resoluci´on de sistemas lineales basados en los subespacios de Krylov, analizando la influencia de su Precondicionamiento y Reordenaci´on. Y en la tercera parte, se adapta la construcci´on de Precondicionadores al caso de los sistemas de ecuaciones de matrices variables. Se ilustra su eficacia mediante numerosos experimentos num´ericos y se destaca la importancia de las t´ecnicas presentadas en la aplicaci´on de Algoritmos Gen´eticos, para la selecci´on de los par´ametros ´optimos del modelo de viento. Finalmente, se extraen las conclusiones oportunas y se exponen las posibles l´ıneas futuras.
´ Indice general 1. INTRODUCCI´ ON 1 1.1. ANTECEDENTES.......................... 1 1.2. OBJETIVOS ............................. 3 1.3. METODOLOG´ IA........................... 4 2. CAMPOS DE VIENTO 8 2.1. PRELIMINARES........................... 8 2.2. MODELOS DE CAMPOS DE VIENTO . . . . . . . . . . . . . . 10 2.3. BREVES NOCIONES DE C´ ALCULO VARIACIONAL . . . . . . 13 2.3.1. FUNCIONAL: SU DEFINICI´ ON .............. 13 2.3.2. C´ ALCULO VARIACIONAL . . . . . . . . . . . . . . . . . 14 2.3.3. ECUACIONES DE EULER . . . . . . . . . . . . . . . . . 14 2.3.4. PROBLEMAS VARIACIONALES CON LIGADURAS . . 16 2.4. MODELO DE MASA CONSISTENTE . . . . . . . . . . . . . . . 17 2.5. CONSTRUCCI´ ON DEL CAMPO INICIAL . . . . . . . . . . . . 21 2.5.1. INTERPOLACI´ ON HORIZONTAL . . . . . . . . . . . . . 21 2.5.2. EXTRAPOLACI´ ON VERTICAL . . . . . . . . . . . . . . 22 3. DISCRETIZACI´ ON MEDIANTE ELEMENTOS FINITOS 30 3.1. GENERALIDADES ......................... 30 3.2. MALLAS ADAPTATIVAS . . . . . . . . . . . . . . . . . . . . . 31 3.3. GENERACI´ ON DE MATRICES VARIABLES . . . . . . . . . . . 37 4. ESTIMACI´ ON DE PAR´ AMETROS 39 4.1. CONSIDERACIONES PREVIAS . . . . . . . . . . . . . . . . . . 39 4.2. ALGORITMOS GENETICOS . . . . . . . . . . . . . . . . . . . . 42
´ Indice general iv 5. M´ ETODOS ITERATIVOS BASADOS EN SUBESPACIOS DE KRYLOV 47 5.1. PRELIMINARES........................... 47 5.2. SUBESPACIOS DE KRYLOV . . . . . . . . . . . . . . . . . . . . 48 5.3. M´ ETODO DEL GRADIENTE . . . . . . . . . . . . . . . . . . . 50 5.4. M´ ETODO DEL GRADIENTE CONJUGADO . . . . . . . . . . . 53 5.5. OTROS M´ ETODOS DE KRYLOV . . . . . . . . . . . . . . . . . 58 5.5.1. M´ ETODOS DE ORTOGONALIZACI´ ON ......... 59 5.5.2. M´ ETODOS DE BIORTOGONALIZACI´ ON........ 60 5.5.3. M´ ETODOS BASADOS EN LA ECUACI´ ON NORMAL . 64 6. PRECONDICIONAMIENTO 66 6.1. CONSIDERACIONES PREVIAS . . . . . . . . . . . . . . . . . . 66 6.2. CONDICIONAMIENTO DE UN SISTEMA . . . . . . . . . . . . 67 6.3. T´ ECNICAS DE PRECONDICIONAMIENTO . . . . . . . . . . . 69 6.4. M´ ETODO DEL GRADIENTE CONJUGADO PRECONDICIONADO................................. 72 6.5. PRECONDICIONADORES EXPL´ ICITOS E IMPL´ ICITOS . . . 75 6.6. PRECONDICIONADORES EXPL´ ICITOS............. 76 6.6.1. PRECONDICIONADOR AINV . . . . . . . . . . . . . . . 76 6.6.2. PRECONDICIONADOR SAINV . . . . . . . . . . . . . . 79 6.7. PRECONDICIONADORES IMPL´ ICITOS............. 81 6.7.1. POR COMPARACI´ ON CON EL M´ ETODO DE RICHARDSON.............................. 81 6.7.2. POR FACTORIZACIONES INCOMPLETAS . . . . . . . 84 7. REORDENACI´ ON 89 7.1. PRELIMINARES........................... 89 7.2. ALGORITMO DE CUTHILL-McKEE INVERSO (RCM) . . . . . 91 7.3. ALGORITMO DEL M´ INIMO VECINO (MN) . . . . . . . . . . . 91 7.4. ALGORITMO DE GEORGE . . . . . . . . . . . . . . . . . . . . 92 7.5. ALGORITMO MULTICOLORING (MC) . . . . . . . . . . . . . 93
´ Indice general v 8. PRECONDICIONAMIENTO DE SISTEMAS DE MATRIZ VARIABLE 95 8.1. PROPUESTA DE ESTRATEGIA . . . . . . . . . . . . . . . . . . 95 8.2. ADAPTACI´ ON DEL PRECONDICIONADOR SAINV . . . . . . 97 8.3. ADAPTACI´ ON DE LA FACTORIZACI´ ON DE CHOLESKY . . 99 9. EXPERIMENTOS NUM´ ERICOS 102 9.1. APLICACIONES TEST . . . . . . . . . . . . . . . . . . . . . . . 102 9.1.1. PRELIMINARES . . . . . . . . . . . . . . . . . . . . . . . 102 9.1.2. EJEMPLO1 ......................... 104 9.1.3. EJEMPLO2 ......................... 107 9.1.4. EJEMPLO3 ......................... 110 9.2. ELECCI´ ON DEL PAR´ AMETRO ´ OPTIMO ............ 116 9.2.1. PRELIMINARES . . . . . . . . . . . . . . . . . . . . . . . 116 9.2.2. EJEMPLO4 ......................... 118 9.2.3. EJEMPLO5 ......................... 119 9.2.4. AN´ ALISIS DE RESULTADOS . . . . . . . . . . . . . . . 120 10.CONCLUSIONES Y LINEAS FUTURAS 122 10.1.CONCLUSIONES .......................... 122 10.2. LINEAS FUTUTRAS . . . . . . . . . . . . . . . . . . . . . . . . 125 BIBLIOGRAF´ IA ............................. 127
METODOLOG´ IA 4 precondicionadores, tanto expl´ıcitos como impl´ıcitos. Exponer las t´ecnicas de reordenaci´on y su efecto en la resoluci´on de sistemas precondicionados. Implementar los precondicionadores SAINV y los obtenidos como consecuencia de la factorizaci´on incompleta de Cholesky, para su aplicaci´on a los sistemas de ecuaciones lineales de matrices variables. Realizar un estudio comparativo del efecto que producen los precondicionadores, anteriormente citados, sobre la convergencia del algoritmo del Gradiente Conjugado, al aplicarlos a los sistemas de ecuaciones lineales de matrices variables, obtenidos en la modelizaci´on pr´actica de diversos campos de viento de la isla de Gran Canaria. Comprobar la influencia de la reordenaci´on en la soluci´on de varios de esos sistemas lineales de coeficientes variables, obtenidos en la modelizaciones pr´acticas mencionadas. Realizar, asimismo, un estudio comparativo de la convergencia del algoritmo del Gradiente Conjugado, con la utilizaci´on de los precondicionadores propuestos, en la selecci´on, mediante Algoritmos Gen´eticos, de los valores param´etricos ´optimos, para diferentes campos de viento. 1.3. METODOLOG´ IA El presente trabajo se inicia con la Introducci´on, Cap´ıtulo 1, donde se presentan los antecedentes, los objetivos propuestos en el desarrollo de estas tesis, y la metodolog´ıa seguida a lo largo de la misma. El contenido b´asico, se estructura en tres grandes bloques. Un primer bloque, que abarca desde el cap´ıtulo 2 hasta el cap´ıtulo 4, en el que se expone la problem´atica de la modelizaci´on de los campos de viento, su discretizaci´on por MEF y loa procesos de estimaci´on de los par´ametros que intervienen en su formulaci´on. Un segundo bloque, que va del cap´ıtulo 5 al 7, en el que se da una visi´on general del estado del arte de los m´etodos iterativos, para la resoluci´on de grandes
METODOLOG´ IA 5 sistemas de ecuaciones lineales, basados en los subespacios de Krylov, as´ı como de las t´ecnicas de precondicionamiento y reordenaci´on de dichos sistemas. Y un tercer bloque, constituido por el cap´ıtulo 8, donde se trata de hacer una nueva aportaci´on, extendiendo las t´ecnicas de precondicionamiento a los sistemas de ecuaciones lineales de matrices variables. Por ´ultimo, en el cap´ıtulo 9, se presentan los resultados de diversos experimentos num´ericos, donde se aplican los precondicionadores propuestos a los sistemas de matrices variables, se utilizan diferentes reordenaciones y se comprueba su influencia sobre la convergencia del algoritmo del Gradiente Conjugado. El primer bloque, dedicado a la modelizaci´on de campos de viento est´a constituido por: El Cap´ıtulo 2, donde se exponen el estado del arte de la modelizaci´on de los campos de viento, as´ı como, unas breves nociones de C´alculo Variacional, dada su aplicaci´on en la formulaci´on matem´atica de los modelos en cuesti´on. Se presta especial atenci´on al Modelo de Masa Consistente, por considerarse, actualmente, como el m´as eficiente, llegando a establecer la ecuaci´on diferencial el´ıptica que lo define. Asimismo, se analiza la construcci´on del campo inicial, por su trascendencia en la consecuci´on de un resultado final correcto. El Cap´ıtulo 3, en el que se afronta la discretizaci´on, por el M´etodo de Elementos Finitos, de la ecuaci´on el´ıptica, mencionada anteriormente, correspondiente a los Modelos de Masa Consistente. Se expone la problem´atica de la generaci´on de mallas adaptativas, para conseguir discretizar con m´as eficacia, sobre todo en los terrenos de orograf´ıa compleja. Llegando finalmente a la generaci´on del sistema de ecuaciones lineales del modelo matem´atico, que resulta ser de matriz Sim´etrica Definida Positiva (SDP), de coeficientes variables, dependientes de un par´ametro, denominado par´ametro de estabilidad del modelo. El Cap´ıtulo 4, dedicado al an´alisis y valoraci´on, no solo, de los par´ametros que intervienen an la construcci´on del campo inicial, sino tambi´en, del par´ametro de estabilidad, por su gran protagonismo en el sistema de ecuaciones
METODOLOG´ IA 6 lineales obtenido. El segundo bloque, correspondiente a los m´etodos del C´alculo Num´erico para la resoluci´on de grandes sistemas de ecuaciones lineales, est´a formado por: El Cap´ıtulo 5, en el que se exponen los fundamentos de los m´etodos iterativos basados en los subespacios de Krylov, por considerarse los m´as adecuados para la resoluci´on de grandes sistemas lineales. Prestando especial atenci´on al algoritmo del Gradiente Conjugado, por tratarse del m´etodo m´as eficaz para resolver los sistemas de matrices SDP, que, como se ha indicado, es el tipo de sistema que se obtiene en la modelizaci´on de campos de viento de Masa Consistente. El Cap´ıtulo 6, dedicado por entero al Precondicionamiento de sistemas, por ser una herramienta que mejora sensiblemente la convergencia de los m´etodos de Krylov. Adem´as de exponer el estado del arte de esta t´ecnica, se establece el algoritmo del Gradiente Conjugado Precondicionado, pues es el m´etodo que se utilizar´a m´as adelante para afrontar los experimentos num´ericos y, tambi´en, se describen los precondicionadores expl´ıcitos AINV y SAINV, as´ı como, los impl´ıcitos, tanto los obtenidos por comparaci´on con el m´etodo de Richardson, como los que se fundamentan en factorizaciones incompletas de la matriz inicial del sistema. El Cap´ıtulo 7, en el que se estudian los m´etodos m´as pr´acticos de Reordenaci´on de sistemas, como son: el algoritmo de Cuthill-McKee Inverso, el del M´ınimo Vecino y el Multicoloring, que basados en la teor´ıa de grafos, proporcionan matrices con ancho de banda o perfil menor, lo que incide notablemente en una mayor simplificaci´on, a la hora de construir un precondicionador m´as eficaz. El tercer bloque, donde se hace la aportaci´on novedosa de esta tesis, est´a constituido por: El Cap´ıtulo 8, en el que, a partir del conocimiento de los principales precondicionadores aplicables a los sistemas de matrices SDP, como son los SAINV y los basados en la factorizaci´on incompleta de Cholesky (ICHOL),
METODOLOG´ IA 7 se implementa su adaptaci´on a los sistemas de ecuaciones lineales de matrices variables, que, como ya se ha indicado, son los que se obtienen en la modelizaci´on matem´atica de los campos de viento. Los resultados de los experimentos num´ericos, realizados para confrontar las propuestas aportadas en esta tesis, se recogen en El Cap´ıtulo 9, que a su vez se distribuye en dos secciones. Una, dedicada a las aplicaciones test, donde se han utilizado ambos tipos de precondicionadores, SAINV e ICHOL, sobre sistemas de ecuaciones de matrices variables obtenidos en tres modelos de campos de viento distintos, conseguidos con tres discretizaciones diferentes, sobre una regi´on de la isla de Gran Canaria, comprobando en la pr´actica, la eficacia de los distintos precondicionadores, as´ı como, la influencia de las diferentes t´ecnicas de Reordenaci´on. Otra, orientada a la optimizaci´on del par´ametro de estabilidad, donde se han aplicado los precondicionadores tipo ICHOL, dado que, con los resultados obtenidos en los experimentos anteriores, han demostrado ser los m´as eficaces para estos modelos. Se han utilizado sobre los sistemas de ecuaciones de matrices variables correspondientes a la modelizaci´on de dos campos de viento diferentes, usando siempre el algoritmo del Gradiente Conjugado Precondicionado. En cada uno de estos ejemplos se recogen los tiempos de computaci´on, para dos gamas distintas de valores del par´ametro, pues estos resultados influir´an notablemente a la hora de seleccionar el precondicionador m´as adecuado para usar en la selecci´on del valor ´optimo del par´ametro, mediante Algoritmos Gen´eticos, dado la repetitividad del proceso. Las conclusiones obtenidas y las futuras lineas de investigaci´on, con las que finaliza este trabajo, se recogen en el Cap´ıtulo 10.
Cap´ıtulo 2 CAMPOS DE VIENTO 2.1. PRELIMINARES El viento no es otra cosa que el aire en movimiento, entendiendo por aire la masa de gases que constituyen nuestra atm´osfera terrestre. Hace unos cuatro mil seiscientos millones de a˜nos el Sistema Solar se condens´o a partir de una nube de gas y polvo interestelar, la Nebulosa Solar. Las atm´osferas de la Tierra, Venus y Marte se formaron a partir de materia vol´atil que escap´o de cada planeta. La primitiva atm´osfera de la Tierra estaba compuesta por di´oxido de carbono (CO2), nitr´ogeno (N2) y vapor de agua (H2O), con trazas de hidr´ogeno (H2), una mezcla muy similar a la emitida hoy en d´ıa por los volcanes. La aparici´on del ox´ıgeno (O2) como componente de la atm´osfera fue el resultado de su producci´on por la actividad fotosint´etica. Se estima que el nivel actual de O2 se alcanz´o hace aproximadamente cuatrocientos millones de a˜nos y se mantiene gracias a un balance entre su producci´on por fotos´ıntesis y su desaparici´on por la respiraci´on de los seres vivos y gradual descenso del carbono org´anico. La atm´osfera actual est´a compuesta principalmente por los gases N2(78 %), O2(21 %), Ar(1 %) y una proporci´on muy variable de vapor de agua, que puede alcanzar hasta un 3 %. Desde la m´as remota antig¨ uedad, el hombre se di´o cuenta que el viento pod´ıa ser aprovechado como fuente de energ´ıa, as´ı los egipcios navegaban ya a vela en el a˜no 4500 a.C. M´as tarde el aprovechamiento de la energ´ıa e´olica continu´o con la aparici´on de los molinos, se tienen datos de ellos desde el siglo II a.C. Los
PRELIMINARES 9 m´as antiguos eran de eje vertical, pero hacia el siglo VIII aparecieron en Europa, procedentes del Este, los grandes molinos de eje horizontal con cuatro aspas. A partir de los siglos XII y XIII se generaliza el uso de los molinos de viento para la molienda de granos y para la elevaci´on de agua, actividades que se mantienen hasta bien entrado el siglo XIX. La llegada de la revoluci´on industrial, con la utilizaci´on masiva del vapor, la electricidad y los combustibles f´osiles como fuente de energ´ıa, interrumpe su desarrollo. Sin embargo, en la segunda mitad del siglo XIX, aparece el popular molino americano multipala, utilizado para el bombeo de agua, practicamente en todo el mundo, y cuyas caracter´ısticas habr´ıan de sentar las bases para el dise˜no de los modernos aerogeneradores. Durante el siglo XX, la superficie del planeta se ha ido cubriendo paulatinamente de m´as y m´as parques e´olicos, que van surgiendo como alternativa viable de las centrales t´ermicas. La sociedad ha ido adquiriendo conciencia de los problemas medio ambientales y valora cada vez m´as el uso de las energ´ıas renovables. Esta creciente inquietud social ha adquirido una gran importancia desde el punto de vista pol´ıtico (no hay formaci´on pol´ıtica que se sustraiga a los problemas ecol´ogicos y no los incluya en su programa electoral) y econ´omico (las empresas invierten cada d´ıa m´as en estudios de impacto ambiental y en publicidad para alardear de sus valores ecol´ogicos, sean reales o no), lo que ha producido en los ´ultimos a˜nos un crecimiento notable de la producci´on de energ´ıa el´ectrica de origen e´olico. La Conferencia de Madrid, marzo de 1994, consider´o viable que las energ´ıas renovables contribuyeran con un 15 % a la demanda total de energ´ıa primaria en la CE, antes de 2010. Espa˜na ocupa un lugar destacado en el panorama e´olico comunitario, con el quinto puesto por potencia e´olica instalada, detr´as de Dinamarca, Alemania, Reino Unido y Holanda. Las empresas del sector necesitan herramientas cada vez m´as sofisticadas que les permita hacer frente a las demandas de un mercado cada vez m´as competitivo y exigente. Por otra parte el desarrollo industrial ha tra´ıdo como consecuencia el vertido masivo a la atm´osfera de sustancias contaminantes. Es cada vez m´as evidente que la contaminaci´on atmosf´erica tiene graves repercusiones que provocan la alteraci´on de las condiciones medioambientales del planeta. Las consecuencias de esta contaminaci´on van desde la lluvia ´acida, hasta el aumento de los trastornos respi-
MODELOS DE CAMPOS DE VIENTO 10 ratorios y al´ergicos de la poblaci´on, pasando por el preocupante efecto invernadero de graves consecuencias a largo plazo. 2.2. MODELOS DE CAMPOS DE VIENTO Los modelos de campos de viento son herramientas que permiten afrontar diversos problemas relacionados con el impacto del viento en nuestro entorno, tales como, el estudio de sus efectos sobre una determinada estructura (especialmente puentes y edificios de gran altura), la dispersi´on de contaminantes en la atm´osfera, la propagaci´on de incendios o el estudio del emplazamiento y rendimiento de parques e´olicos. Concretamente en este ´ultimo segmento, con los modelos de campo de viento se pueden afrontar diversos problemas surgidos en el seno de las empresas dedicadas a la explotaci´on de este tipo de parques, tales como la evaluaci´on de la potencia producida por un aerogenerador, en funci´on de su situaci´on, y su comparaci´on con las curvas suministradas por el fabricante; el estudio de la ubicaci´on ´optima de la red de estaciones de medida previa a la instalaci´on del parque y la distribuci´on m´as eficaz de los aeorgeneradores. Los modelos de campos de viento son, pues, herramientas cada vez m´as importantes para afrontar con eficacia una amplia gama de problemas de inter´es social, pol´ıtico y econ´omico, y cada vez se exige m´as de ellos. Inicialmente los modelos meteorol´ogicos se dividen en dos grandes grupos: los Modelos F´ısicos y los Modelos Matem´aticos. Los primeros usan t´uneles de viento sobre reproducciones a peque˜na escala del terreno en estudio. En los Modelos Matem´aticos, que ser´an los objetos de esta tesis, se emplean t´ecnicas algebraicas y de c´alculo para resolver ecuaciones metereol´ogicas. A su vez los Modelos Matem´aticos se dividen en Anal´ıticos y Num´ericos; los primeros, por la gran complejidad de las ecuaciones que describen la atm´osfera, hacen muy dif´ıcil la resoluci´on exacta en dominios irregulares, mientras que los segundos ofrecen mejores perspectivas. Por ello el objetivo de esta tesis va a ser tratar de aportar herramientas de C´alculo Num´erico que permitan la modelizaci´on matem´atica de campos de viento con la mayor eficacia posible. Seg´un la extensi´on del dominio a estudiar, los Modelos de Viento se pueden
MODELOS DE CAMPOS DE VIENTO 11 clasificar en Modelos de Macroescala, cuando el ´area de estudio abarca desde un continente hasta el globo terr´aqueo completo. Modelos de Mesoescala, cuando se refieren a extensiones que van desde unos pocos kil´ometros hasta alrededor de cien, y de Microescala, para regiones locales que tienen como m´aximo alrededor de un kil´ometro. Los Modelos Matem´aticos de Campos de Viento tambi´en se pueden clasificar en Modelos de Pron´ostico o Din´amicos y Modelos de Diagn´ostico o Cinem´aticos. Los Modelos de Pron´ostico se basan en la soluci´on de ecuaciones hidrodin´amicas y termodin´amicas que depende del tiempo (llamadas tambi´en ecuaciones primitivas porque derivan directamente de los principios de conservaci´on) modificadas para su aplicaci´on a la atm´osfera. Sin embargo la soluci´on del conjunto completo de ecuaciones sigue siendo una tarea costosa. Adem´as, cuanto m´as elaborado es el modelo, m´as fiables deben ser los datos de entrada para aprovechar las ventajas ofrecidas, y con frecuencia estos datos no suelen estar disponibles. Autores como Lalas et al. [58] incluyen en los Modelos Din´amicos algunos c´odigos que introducen aproximaciones a las ecuaciones primitivas, a la vez que desprecian su dependencia del tiempo. Estos c´odigos, denominados JH, se basan en una propuesta realizada por Jackson y Hunt [51] y son ampliamente utilizados. Como ejemplo podemos citar el modelo empleado por Troen y Petersenen [110] en la confecci´on del Atlas Europeo de Viento. Los Modelos de Diagn´ostico deben su nombre a que no se utilizan para realizar previsiones a trav´es de la integraci´on de las relaciones conservativas. Eliminan directamente de sus ecuaciones la dependencia del tiempo y por esta raz´on se les llama tambi´en Cinem´aticos. Estos modelos generan un campo de viento que satisface algunas restricciones f´ısicas. Si la ´unica restricci´on que se les impone es que cumplan la ecuaci´on de continuidad, lo que supone la conservaci´on de la masa, el modelo se denomina de Masa Consistente. Los Modelos de Diagn´ostico no requieren muchos datos de entrada y son f´aciles de usar, por lo que resultan atractivos desde el punto de vista pr´actico. Autores como Pennel [81] han comprobado que en algunos casos los Modelos de Masa Consistente mejorados, tales como NOABL y COMPLEX, superaron los resultados de Modelos Din´amicos m´as complicados y costosos. Sin embargo hay que tener en cuenta que los Modelos de Diagn´ostico no
MODELOS DE CAMPOS DE VIENTO 12 consideran los efectos t´ermicos ni los debidos a cambios de gradientes de presi´on. Por ello, flujos tales como las brisa marinas, vientos en pendiente y otros tales como los de separaci´on a favor del viento, no pueden simularse con esta modelos, a no ser que se incorporen en los datos de viento inicial, a partir de observaciones realizadas en lugares apropiados a tal efecto [54, 76]. Los Modelos de Diagn´ostico est´an dise˜nados espec´ıficamente para predecir los efectos de la orograf´ıa sobre el flujo medio de viento considerado de manera estacionaria, esto es, flujos promediados en intervalos de tiempo entre 10 minutos y 1 hora. NOABL [82] es un modelo meteorol´ogico que proporciona una representaci´on precisa del terreno gracias a una transformaci´on de la coordenada vertical en la que la coordenada m´as baja es conforme a la superficie del terreno. Posteriormente, diversos autores [55, 56, 108] realizaron algunas modificaciones en la inicializaci´on del mismo, con el fin de que el modelo considerara el efecto de la rugosidad del terreno sobre el perfil del viento, de forma que el modelo dispusiera de perfiles m´as realistas que los del c´odigo original y tuviera en cuenta el cambio debido a la fuerza de Coriolis en la direcci´on del viento, en la capa l´ımite atmosf´erica. Los c´odigos resultantes fueron bautizados como NOABL* [58] y EOLOS [108], aunque actualmente es m´as conocido por AIOLOS, por su referencia a la etimolog´ıa griega. M´as tarde se introdujeron modificaciones que describen con m´as precisi´on los perfiles de viento en diferentes condiciones de estabilidad. Este nuevo c´odigo [86] se llam´o WINDS(Wind-field Interpolation by Non Divergent Squemes). Ambos modelos, AIOLOS y WINDS, usan datos de estaciones situadas en tierra y, opcionalmente, datos observados tanto de perfiles verticales como de viento geostr´ofico. Tambi´en utilizan coordenadas conformes al terreno. Ambos se basan en la minimizaci´on de los cuadrados de las diferencias entre las velocidades de un viento inicial, obtenido por interpolaci´on de los datos observados, y el viento a ajustar, sujeto a la restricci´on de que el campo de viento ajustado ha de tener divergencia nula [99]. Adem´as de los ya citados, existe toda una gama de Modelos de Diagn´ostico que la comunidad cient´ıfica ha venido utilizando en problemas relacionados con la meteorolog´ıa y/o con la contaminaci´on atmosf´erica. Dentro de los m´as conocidos
MODELO DE MASA CONSISTENTE 19 Utilizando el grupo de ecuaciones (2.3), que hemos visto anteriormente para las Ecuaciones de Euler, resulta que: F0 p1= Φ; F0 q1= 0; F0 r1= 0 F0 p2= 0; F0 q2= Φ; F0 r2= 0 F0 p3= 0; F0 q3= 0; F0 r3= Φ Y sustituyendo y derivando convenientemente tendremos: 2(u−uo)α2 1−∂Φ ∂x = 0 2(v−v0)α2 1−∂Φ ∂y = 0 2(ω−ω0)α2 2−∂Φ ∂z = 0 u=1 2α2 1 ∂Φ ∂x +u0 v=1 2α2 1 ∂Φ ∂y +v0 ω=1 2α2 2 ∂Φ ∂z +ω0 Que pueden resumirse como: ~u =~v0+T~ ∇Φ (2.10) Expresi´on en la que T= (Th, Th, Tv), se puede considerar como un tensor diagonal de transmisi´on, tal que: T=diagh1 2α2 1 ,1 2α2 1 ,1 2α2 2i(2.11) constante para un dominio dado, Ω. As´ı el campo soluci´on, ~u, se obtendr´a a partir del campo inicial, ~v0, y de los valores de Φ, para cuyo c´alculo se recurre a la ecuaci´on de condici´on ~ ∇·~u = 0,o sea,∂u ∂x +∂v ∂y +∂ω ∂z = 0 resultando la ecuaci´on diferencial siguiente: 1 2α2 1 ∂2Φ ∂x2+∂u0 ∂x +1 2α2 1 ∂2Φ ∂y2+∂v0 ∂y +1 2α2 2 ∂2Φ ∂z2+∂ω0 ∂z = 0. Teniendo en cuenta que los par´ametros α1yα2, se suelen considerar constantes en todo el dominio Ω, la ecuaci´on anterior se puede simplificar multiplicando todos los t´erminos por 2 α2 1y llamando a α2 1 α2 2 =Tv Th =α2=. (2.12) se obtiene ∂2Φ ∂x2+∂2Φ ∂y2+∂2Φ ∂z2=−1 Th∂u0 ∂x +∂v0 ∂y +∂ω0 ∂z (2.13)
MODELO DE MASA CONSISTENTE 20 ecuaci´on el´ıptica en Φ, donde recibe el nombre de par´ametro de estabilidad del modelo. Teniendo en cuenta las expresiones (2.6) y (2.10), esta ecuaci´on diferencial, en derivadas parciales, estar´a sujeta a una condici´on tipo Neuman en las fronteras impermeables (terreno y frontera superior); ya que ~n ·T~ ∇Φ = −~n ·~v0en Γb(2.14) que se completa con la condici´on de Dirichlet, nula en las fronteras permeables (fronteras verticales del dominio): Φ = 0 en Γa(2.15) Obs´ervese que en la frontera superior, al ser el campo inicial ~v0horizontal, la condici´on (2.14) se transforma en ~n ·T~ ∇Φ = 0 (2.16) Por tanto, desde el punto de vista matem´atico, la parte esencial de la construcci´on de un modelo de viento de Masa Consistente, se convierte en la resoluci´on de la ecuaci´on diferencial en derivadas parciales (??), con las condiciones de contorno (2.14) y (2.15). Obs´ervese que el par´ametro de estabilidad del modelo, , va a estar presente en la soluci´on del problema, de ah´ı que los modelos de Masa Consistente sean criticados por su alta dependencia de par´ametros. La elecci´on acertada de sus valores es de gran importancia para la fiabilidad de los resultados. Teniendo en cuenta (2.11) y (2.12): =Tv Th (2.17) resultando entonces que para 1, predomina el ajuste del flujo en la direcci´on vertical, es decir, el aire tiende a sobrepasar las barreras del terreno, m´as que a pasar horizontalmente alrededor de ellas. Mientras que para 1, el ajuste del flujo ocurre primeramente en el plano horizontal, por tanto el aire pasar´a alrededor de las barreras del terreno m´as que sobre ellas. En particular, → ∞ significa ajuste vertical puro, mientras que →0 significa ajuste horizontal puro [7].
CONSTRUCCI ´ ON DEL CAMPO INICIAL 21 2.5. CONSTRUCCI´ ON DEL CAMPO INICIAL Para la construcci´on del campo inicial partimos de los valores de la velocidad del viento y de su direcci´on obtenidos en las estaciones de medida. Los datos de viento se toman de estaciones ubicadas en el dominio de estudio. Cada estaci´on de medida proporciona la velocidad (en m/s) y direcci´on del viento a una altura zssobre el nivel del terreno (t´ıpicamente 10 metros). La direcci´on del viento viene dada en grados sexagesimales medidos en sentido horario y tomando como referencia la direcci´on norte. As´ı el norte se corresponde a 0 grados, el sur a 180 grados, el este a 90 grados y el oeste a 270 grados. A efectos de c´alculo en el modelo, es necesario obtener el ´angulo medido en sentido antihorario, tomando como referencia el semieje positivo horizontal. Por otro lado, como las estaciones miden el viento en intervalos discretos de tiempo, en general es necesario interpolar las medidas para calcular el viento en un instante concreto. El campo inicial ~v0se construye en dos etapas: 2.5.1. INTERPOLACI´ ON HORIZONTAL En primer lugar, se calcula mediante interpolaci´on horizontal el valor de ~v0en los puntos del dominio situados a la misma altura zs(sobre el terreno) que las estaciones de medida. La t´ecnica m´as com´un de interpolaci´on se formula en t´erminos de la inversa de la distancia al cuadrado entre el punto y la estaci´on de medida [116]. Sin embargo, otros autores usan simplemente la altitud de los puntos de medida [79]. Aqu´ı se propone una f´ormula que tiene en cuenta ambas consideraciones, ~v0(ze) = β N P n=1 ~vn d2 n N P n=1 1 d2 n + (1 −β) N P n=1 ~vn |∆hn| N P n=1 1 |∆hn| (2.18) El valor de ~vncorresponde a la velocidad observada en la estaci´on n,Nes el n´umero de estaciones utilizadas en la interpolaci´on, dnes la distancia horizontal desde la estaci´on nal punto donde estamos calculando la velocidad del viento,
CONSTRUCCI ´ ON DEL CAMPO INICIAL 22 |∆hn|es la diferencia de altura entre la estaci´on ny el punto en estudio, y βes un par´ametro de peso que toma valores entre 0 y 1. Cuando β→1 aumenta la importancia de la distancia horizontal desde cada punto a las estaciones de medida. Esta aproximaci´on se emplea en problemas con una orograf´ıa regular o en an´alisis bidimensionales. De manera an´aloga, si β→0 es entonces la diferencia de altura entre cada punto y las estaciones de medida la que resulta determinante, en detrimento de la distancia horizontal. Esta segunda aproximaci´on es la que se usa cuando la orograf´ıa del terreno es irregular. En la pr´actica, las regiones geogr´aficas estudiadas suelen combinar zonas de orograf´ıa irregular con otras de orograf´ıa mas regular, por lo que tomar un valor intermedio para βsuele ser lo m´as apropiado. 2.5.2. EXTRAPOLACI´ ON VERTICAL Con la informaci´on obtenida en el paso anterior se realiza una extrapolaci´on vertical para definir el campo de velocidades en la totalidad del dominio. El viento se desarrolla, en primer lugar, como consecuencia de diferencias espaciales en la presi´on atmosf´erica. Estas diferencias de presi´on normalmente son causadas por una diferente absorci´on de la radiaci´on solar. En un plano horizontal, el viento fluye de las zonas de alta presi´on a zonas de baja presi´on y verticalmente de zonas de baja presi´on a zonas de alta presi´on. La velocidad del viento es proporcional al cambio de presi´on por unidad de distancia o gradiente de presi´on. Las zonas con presiones similares se representan en los mapas meteorol´ogicos unidas mediante l´ıneas imaginarias denominadas isobaras. Cuanto m´as juntas est´an unas isobaras, mayor ser´a la fuerza del viento. Un segundo factor que afecta el movimiento del aire es la fuerza de Coriolis, debida a la rotaci´on terrestre. El par´ametro f= 2Θ sen φlse denomina par´ametro de Coriolis, siendo Θ = 7,292 ×10−5s−1la velocidad de rotaci´on de la Tierra y φlla latitud. Se considera positiva en el hemisferio norte, nula en el ecuador y negativa en el hemisferio sur. En tercer lugar puede aparecer una aceleraci´on centr´ıpeta, cuando el viento gira en torno a un centro. Por ´ultimo, aparece la fricci´on debida al desplazamiento del aire. Los vientos influenciados por el gradiente de presi´on y la fuerza de Coriolis
CONSTRUCCI ´ ON DEL CAMPO INICIAL 23 se denominan vientos geostr´oficos. Estratificaci´on atmosf´erica: En este modelo se considera una divisi´on de la capa m´as baja de la atm´osfera en distintas subcapas, en las que la extrapolaci´on vertical de las velocidades de viento se realiza de forma diferente, como puede observarse en la figura 2.1. ~v0(z) = ~v∗ klog z z0−Φm Zs Z0 Zsl hm Zpbl Capa de Mezcla Viento Geostr´ofico ~v0(z) = ρ(z)~v0(zsl) + [1 −ρ(z)]~vg Figura 2.1: Perfil vertical de viento definido sobre cada capa de la estratificaci´on atmosf´erica. As´ı, la capa l´ımite planetaria est´a situada a una altitud zpbl sobre el nivel del terreno, y es la capa de la atm´osfera, situada por debajo de la atm´osfera libre, que est´a afectada directamente por la fricci´on de la superficie de la tierra (conocida tambi´en como capa l´ımite atmosf´erica). La altitud de la capa l´ımite planetaria zpbl sobre el terreno se ha tomado tal que la direcci´on e intensidad del viento es constante a partir de esa altura [101]: zpbl =γ|~v∗| f(2.19) siendo γun par´ametro comprendido entre 0,15 y 0,45 que depende de la estabilidad de la atm´osfera y ~v∗la velocidad de fricci´on que ser´a definida m´as adelante a partir
CONSTRUCCI ´ ON DEL CAMPO INICIAL 24 de los valores obtenidos en la interpolaci´on horizontal. La capa de mezcla, tambi´en llamada capa l´ımite convectiva, es la capa l´ımite atmosf´erica sujeta a fen´omenos convectivos causados por el calor superficial. El aire est´a bien mezclado, es decir, el viento y el potencial de temperatura son pr´acticamente constantes con la altura. La altitud de la capa de mezcla hmse considerar´a igual a zpbl para condiciones neutras e inestables. En condiciones estables se aproxima por hm=γ0s|~v∗|L f(2.20) donde usualmente se toma el par´ametro γ0= 0,4 [117] y Les la longitud de Monin-Obukov, que se calcula a trav´es de la f´ormula de Liu [85], 1 L=azb 0(2.21) con ayb, definidas por la clase de estabilidad de Pasquill (Ver Tabla 2.1): Clase de estabilidad de Pasquill a b A (Extremadamente Inestable) -0.08750 -0.1029 B (Moderadamente Inestable) -0.03849 -0.1714 C (Ligeramente Inestable) -0.00807 -0.3049 D (Neutra) 0.00000 0.0000 E (Ligeramente Estable) 0.00807 -0.3049 F (Moderadamente Estable) 0.03849 -0.1714 Tabla 2.1: Coeficientes aybpara el c´alculo de la longitud de Monin Obukov seg´un la clase de estabilidad de Pasquill. La capa superficial, localizada a una altura zsl sobre la superficie, es la capa baja, dentro de la capa l´ımite planetaria, inmediatamente adyacente a la capa de la superficie de la tierra, en la que la fuerza de arrastre de fricci´on es dominante. Conocido el valor de la altura de la capa de mezcla hm, la altitud de la capa superficial se suele fijar en [117] zsl =hm 10 (2.22)
CONSTRUCCI ´ ON DEL CAMPO INICIAL 25 Estabilidad atmosf´erica: El concepto de estabilidad atmosf´erica est´a relacionado tanto con la turbulencia atmosf´erica como con el gradiente vertical de temperatura y las situaciones de inversi´on t´ermica. La estabilidad atmosf´erica nos proporciona una medida cualitativa de las variaciones de la densidad del aire, debidas a los cambios de presi´on y temperatura y que influyen en determinados movimientos atmosf´ericos. Las condiciones atmosf´ericas pueden clasificarse como: Estable: Si una masa de aire sube se encontrar´a rodeada de aire m´as caliente y, por tanto, menos denso que ella, lo que la har´a bajar; y si baja, se encontrar´a rodeada de aire m´as fr´ıo (m´as denso), y tender´a a subir. Esta tendencia que tiene el aire de permanecer en la misma capa es lo que se denomina estabilidad de la estratificaci´on atmosf´erica. Inestable: En condiciones inestables la temperatura potencial disminuye con la altura, increment´andose los movimientos verticales, es decir si el aire sube se encontrar´a rodeado de aire m´as fr´ıo y denso que ´el, y tender´a a seguir subiendo; y si baja se encontrar´a con aire m´as caliente y ligero, y tender´a a seguir bajando. Neutra: Si un volumen de aire (despu´es de un desplazamiento vertical en una capa atmosf´erica sin mezclar con el aire circundante) experimenta una fuerza neta vertical nula, los movimientos ascensionales no se ver´an perturbados por el gradiente t´ermico, entonces la capa atmosf´erica se asume neutralmente estratificada. Bajo tales condiciones, dicho volumen ni tiende a volver a su posici´on original (estratificaci´on estable) ni acelera alej´andose de ella (estratificaci´on inestable). La estabilidad atmosf´erica puede ser caracterizada mediante la tabla definida por Pasquill.(Ver Tabla 2.2 ) Perfil vertical de velocidades de viento: Como se muestra en la Figura 2.1, se considera un perfil logar´ıtmico-lineal [57] en la capa l´ımite planetaria, que tiene en cuenta la interpolaci´on horizontal [70],
CONSTRUCCI ´ ON DEL CAMPO INICIAL 26 Clase de estabilidad de Pasquill Insolaci´on Noche Velocidad del Cubierto viento en la ´o ≥4/8≤3/8 superficie (m/s) Fuerte Moderada Ligera nubes nubes <2 A A-B B - - 2-3 A-B B C E F 3-5 B B-C C D E 5-6 C C-D D D D >6 C D D D D Para A-B, tomar la media de los valores de A y B, etc. Tabla 2.2: Clases de estabilidad de Pasquill seg´un la velocidad del viento en la superficie y la insolaci´on. Insolaci´on fuerte corresponde al mediod´ıa soleado de mitad de verano en Inglaterra; insolaci´on ligera a condiciones similares en mitad del invierno. La noche se refiere al periodo que va desde una hora antes de ponerse el sol hasta una hora despu´es de salir. La clase neutra D deber´ıa ser usada tambi´en, a pesar de la velocidad del viento, para cielos cubiertos durante el d´ıa o la noche, y para cualquier condici´on del cielo durante las horas precedente y siguiente de la noche definida anteriormente. el efecto de la rugosidad en la intensidad y direcci´on del viento, y la estabilidad del aire (neutra, estable o inestable) seg´un la clasificaci´on de Pasquill. En la capa superficial se construye un perfil logar´ıtmico de velocidades de viento definido por, ~v0(z) = ~v∗ klog z z0−Φmz0< z ≤zsl (2.23) donde ~v0es la velocidad del viento, k≃0,4 es la constante de von Karman y zes la altura sobre el terreno del punto estudiado. El t´ermino ~v∗representa la velocidad de fricci´on. En el flujo turbulento atmosf´erico las fuerzas que se oponen al movimiento est´an caracterizadas por la acci´on que ejercen las rugosidades o asperezas propias de la orograf´ıa del terreno. La velocidad de fricci´on se obtiene en cada punto a partir de las medidas interpoladas a la altura de las estaciones
CONSTRUCCI ´ ON DEL CAMPO INICIAL 27 (interpolaci´on horizontal), ~v∗=k ~v0(ze) ln ze z0−Φm(ze)(2.24) Asimismo, z0corresponde a la longitud de rugosidad de la zona. El concepto de longitud de rugosidad viene a definir una altura por encima del terreno diferente de z= 0, donde, en teor´ıa de la capa superficial, la velocidad del viento es cero. El valor de z0depende de las caracter´ısticas del terreno. Una forma de estimarla es mediante valores est´andar para diferentes tipos de terreno [64]; ver figura 2.2. Otros autores la definen como z0=e 30 , donde ees la altura media de los obst´aculos existentes en la zona de estudio. Por ´ultimo, Φm, es una funci´on que depende de la estabilidad del aire [117]: Φm= 0 (neutra) Φm=−5z L(estable) Φm= log "θ2 m+ 1 2θm+ 1 22#−2 arctan θm+π 2(inestable) donde θm= (1 −16 z L)1/4(2.25) El viento geostr´ofico es una buena aproximaci´on al viento real con flujo uniforme en la alta atm´osfera (atm´osfera libre), donde la fricci´on y aceleraciones no son importantes. La forma general de la expresi´on usada para calcular el viento geostr´ofico es la bien conocida ley de resistencia geostr´ofica (geostrophic drag low) [85] |~ Vg|=|~v∗| kslog |~v∗| fz0−A2 +B2(2.26) Los valores de los coeficientes A y B tienden a ser ∼1,8 y ∼1,5 respectivamente, que son los valores aceptables para condiciones neutras de estabilidad atmosf´erica. El viento en la superficie se supone que gira un ´angulo φgcon respecto a ~ Vg, dado por la relaci´on φg= sin−1 −B|~v∗| k|~ Vg|!(2.27) En nuestro modelo, desde zsl hasta zpbl se realiza una interpolaci´on lineal en ρ(z) con el viento geostr´ofico ~vg ~v0(z) = ρ(z)~v0(zsl) + [1 −ρ(z)]~vgcon zsl < z ≤zpbl (2.28)
CONSTRUCCI ´ ON DEL CAMPO INICIAL 28 Z (m) 0 ~ ~ 2 3 4 5 7 8 9 6 10 2 3 4 5 7 8 9 6 10 2 3 4 5 7 8 9 6 10 2 3 4 5 7 8 9 6 10 2 3 4 5 7 8 9 6 10 2 3 4 5 7 8 9 6 10 2 3 4 5 7 8 9 6 10 2 3 4 5 7 8 9 6 10 2 3 4 5 7 8 9 6 10 2 3 4 5 7 8 9 6 10 2 3 4 5 7 8 9 6 10 2 3 4 5 7 8 9 6 10 2 3 4 5 7 8 9 6 10 2 3 4 5 7 8 9 6 10 2 3 4 5 7 8 9 6 10 10 −5 −4 −3 −2 −1 2 6 5 4 3 7 8 9 1 2 3 4 5 Centros de ciudades con edificios muy altos Centros de grandes poblaciones, ciudades Centros de pequeñas poblaciones Alrededores de poblaciones Bosques Areas escarpadas o montañosas Region de nivel medio de bosque Muchos arboles, setos, pocos edificios Muchos setos Pocos arboles, verano Tierras de cultivo Hierba alta( 60 cms) campos cultivados Aeropuertos Llanuras de hierbas medianas Arboles aislados Pocos arboles, invierno Hierba sin cortar Hierba cortada ( 3 cm) Superficie natural nevada (tierras de cultivo) Desierto (llanura) Planicie cubierta de nieve Terreno ondulado Mar abierto en calma Hielo, planicie (llanura) enlodada. Grandes extensiones de agua Viento mar adentro en zonas costeras Figura 2.2: Longitud de rugosidad: valores aproximados de z0para distintos tipos de terreno definidos por McRae (1982)
MALLAS ADAPTATIVAS 34 Despu´es de estos antecedentes nos proponemos describir a continuaci´on el proceso de creaci´on de una malla de tetraedros que respete la topograf´ıa de una regi´on rectangular con una precisi´on determinada, disponiendo ´unicamente de la informaci´on digitalizada del terreno. Este problema posee cierta dificultad debido a la fuerte irregularidad de la superficie del terreno. Por otra parte, deseamos que la malla est´e adaptada, es decir, que exista una densidad de nodos mayor donde sea necesario para definir las caracter´ısticas geom´etricas de nuestro dominio. La malla generada podr´a utilizarse como malla base para la simulaci´on num´erica de procesos naturales en el dominio; por ejemplo, ajuste de campos de viento [116, 70], propagaci´on de fuego [68], contaminaci´on atmosf´erica [115], etc. Estos fen´omenos tienen su mayor efecto en las zonas pr´oximas al terreno, de ah´ı que tambi´en sea deseable que la densidad de nodos aumente al acercarnos a ´este. Sobre esta malla base, adaptada a las caracter´ısticas geom´etricas del dominio, se podr´an aplicar posteriormente algoritmos de refinamiento y desrefinamiento de tetraedros para mejorar la soluci´on num´erica del problema [61, 62, 45, 44]. Estos algoritmos tendr´an un especial inter´es en los problemas evolutivos. Nuestro dominio est´a limitado en su parte inferior por el terreno y en su parte superior por un plano horizontal situado a una altura en la que las magnitudes objeto del estudio puedan ser consideradas estables. Las paredes laterales est´an formadas por cuatro planos verticales, paralelos dos a dos. Las ideas b´asicas para la construcci´on de la malla inicial combinan, por un lado, la utilizaci´on de un algoritmo de refinamiento y desrefinamiento para dominios bidimensionales y, por otro lado, un algoritmo de generaci´on de mallas de tetraedros basado en la triangulaci´on de Delaunay. Es bien conocido que para construir una triangulaci´on de Delaunay es necesario definir una nube de puntos en el dominio y su frontera. Estos nodos ser´an precisamente los v´ertices de los tetraedros que conforman la malla. La generaci´on de puntos en nuestro dominio se realizar´a sobre diferentes capas, reales o ficticias, definidas desde el terreno hasta la frontera superior del dominio. En concreto, se construye una triangulaci´on con una distribuci´on uniforme de puntos en el plano superior del dominio. Esta malla bidimensional puede ser obtenida a partir de la
MALLAS ADAPTATIVAS 35 realizaci´on de un cierto n´umero de refinamientos globales sobre una malla simple o, por ejemplo, puede tambi´en construirse realizando una triangulaci´on de Delaunay sobre la distribuci´on uniforme de puntos establecida. Consideraremos la malla obtenida como el nivel m´as bajo de la secuencia que define la distribuci´on de los puntos en el resto de las capas. Sobre esta malla regular aplicamos a continuaci´on el algoritmo de refinamiento y desrefinamiento, [26, 83], para definir la distribuci´on de los nodos de la capa correspondiente a la superficie del terreno. Para ello, en primer lugar se construye una funci´on que interpola las cotas obtenidas a partir de una digitalizaci´on de la topograf´ıa de la zona rectangular estudiada. En segundo lugar, realizamos una serie de refinamientos globales sobre la malla uniforme hasta conseguir una malla regular capaz de captar la variaci´on topogr´afica del terreno. El m´aximo grado de discretizaci´on viene definido por el nivel de detalle de la digitalizaci´on. Posteriormente, se realizar´a un desrefinamiento sobre estos ´ultimos niveles de malla utilizando como par´ametro de desrefinamiento el m´aximo error de cotas permitido entre la superficie real del terreno y la superficie definida mediante la interpolaci´on a trozos obtenida con la malla bidimensional resultante. Una vez que se ha definido la distribuci´on de nodos sobre el terreno y sobre el plano superior del dominio, comenzamos a distribuir los nodos situados entre ambas capas. Esta distribuci´on se puede realizar mediante diferentes estrategias, en las que interviene una funci´on de espaciado vertical . La caracter´ıstica fundamental de esta funci´on es que el grado de discretizaci´on obtenido sobre la vertical debe disminuir con la altura, o a lo sumo mantenerse constante. Esta nube de puntos ser´a utilizada por nuestro mallador tridimensional basado en la triangulaci´on de Delaunay. Para evitar posibles problemas de conformidad con la superficie del terreno, se propone construir la malla de tetraedros con la ayuda de un paralelep´ıpedo auxiliar. Sobre su cara inferior se sit´uan todos los nodos distribuidos sobre el terreno, proyectados sobre un plano horizontal situado a la altura definida por la cota m´ınima de la regi´on de estudio, y sobre su cara superior se sit´uan los puntos distribuidos en el plano superior del dominio a su altura real. Esto conlleva una transformaci´on de coordenadas, atendiendo a la funci´on de espaciado sobre cada vertical, para situar el resto de puntos en el
MALLAS ADAPTATIVAS 36 paralelep´ıpedo auxiliar. Estos detalles nos asegurar´an que la distancia m´axima entre dos puntos consecutivos sobre la misma vertical del dominio real ser´a siempre igual o inferior que la correspondiente distancia establecida en el paralelep´ıpedo auxiliar. Se define la nube de puntos en el dominio real y se analiza la transformaci´on entre el dominio real y el paralelep´ıpedo auxiliar en el que se construye la malla mediante una versi´on del m´etodo de triangulaci´on de Delaunay [24]. Proponemos cuatro estrategias diferentes para determinar el n´umero de puntos generados sobre la vertical de cada nodo de la malla bidimensional adaptada a la superficie del terreno, y analizamos las caracter´ısticas fundamentales de cada una de ellas. Las dos primeras estrategias generan puntos sobre capas definidas entre el terreno y la frontera superior del dominio. En estos dos casos, el n´umero de capas reales que se desea crear ser´a introducido como dato. En la primera estrategia el grado de concentraci´on de las capas hacia el terreno se impone, mientras que en la segunda se obtiene autom´aticamente en funci´on del tama˜no de los elementos existentes en la malla bidimensional adaptada a la superficie del terreno. En las dos ´ultimas estrategias las capas generadas ser´an virtuales, es decir, no se define un n´umero concreto de superficies interiores al dominio sobre las que se sit´uan los puntos. Por ello, diremos que en estas dos ´ultimas estrategias el n´umero de capas es variable, y ser´a calculado autom´aticamente en funci´on de los tama˜nos de los elementos existentes en la malla bidimensional que define el terreno, o, tambi´en, en la correspondiente a la frontera superior del dominio. En concreto, la tercera estrategia concentrar´a los puntos hacia el terreno en funci´on del tama˜no de los elementos definidos sobre ´el. En cambio, la ´ultima estrategia determina autom´aticamente, para cada nodo del terreno, una funci´on de espaciado vertical con el objeto de respetar las distancias desde el primer punto generado hasta el terreno, y desde el ´ultimo punto generado hasta la frontera superior, en funci´on de los tama˜nos de los elementos existentes sobre ambas superficies. Una vez que se ha construido la triangulaci´on de Delaunay de la nube de puntos en el paralelep´ıpedo, procedemos a situar los puntos en sus posiciones reales manteniendo la topolog´ıa de la malla. Hay que tener en cuenta que este proceso de compresi´on de la malla puede dar lugar a cruces de tetraedros que habr´a que
GENERACI ´ ON DE MATRICES VARIABLES 37 deshacer posteriormente. Asimismo, ser´a aconsejable aplicar una etapa de suavizado para mejorar la calidad de los elementos de la malla resultante. 3.3. GENERACI´ ON DE MATRICES VARIABLES Para la discretizaci´on mediante elementos finitos de la formulaci´on cl´asica del problema el´ıptico que aparece en los modelos de Campos de Viento de Masa Consistente, dada por la expresi´on (13), con las condiciones de contorno (14) y (15) se ha utilizado una malla de tetraedros, generada mediante las t´ecnicas adaptativas ya descritas. Esto conduce a un conjunto de matrices elementales de dimensi´on 4 ×4 asociadas al elemento Ωe, siendo ˆ ψila funci´on de forma correspondiente a su i-´esimo nodo, i= 1,2,3,4, definidos en el elemento de referencia ˆ Ωey|J|el jacobiano de la transformaci´on de Ωeaˆ Ωe, {Ae}ij =Zˆ Ωe{(∂ˆ ψi ∂ξ ∂ξ ∂x +∂ˆ ψi ∂η ∂η ∂x +∂ˆ ψi ∂ϕ ∂ϕ ∂x )(∂ˆ ψj ∂ξ ∂ξ ∂x +∂ˆ ψj ∂η ∂η ∂x +∂ˆ ψj ∂ϕ ∂ϕ ∂x )+ +(∂ˆ ψi ∂ξ ∂ξ ∂y +∂ˆ ψi ∂η ∂η ∂y +∂ˆ ψi ∂ϕ ∂ϕ ∂y )(∂ˆ ψj ∂ξ ∂ξ ∂y +∂ˆ ψj ∂η ∂η ∂y +∂ˆ ψj ∂ϕ ∂ϕ ∂y )+ (3.1) +(∂ˆ ψi ∂ξ ∂ξ ∂z +∂ˆ ψi ∂η ∂η ∂z +∂ˆ ψi ∂ϕ ∂ϕ ∂z )(∂ˆ ψj ∂ξ ∂ξ ∂z +∂ˆ ψj ∂η ∂η ∂z +∂ˆ ψj ∂ϕ ∂ϕ ∂z )}·|J|dξ dη dϕ y de vectores elementales de 4 ×1, {be}i=Zˆ Ωe−1 Th{u0(∂ˆ ψi ∂ξ ∂ξ ∂x +∂ˆ ψi ∂η ∂η ∂x +∂ˆ ψi ∂ϕ ∂ϕ ∂x )+ +v0(∂ˆ ψi ∂ξ ∂ξ ∂y +∂ˆ ψi ∂η ∂η ∂y +∂ˆ ψi ∂ϕ ∂ϕ ∂y )+ (3.2) +w0(∂ˆ ψi ∂ξ ∂ξ ∂z +∂ˆ ψi ∂η ∂η ∂z +∂ˆ ψi ∂ϕ ∂ϕ ∂z )}·|J|dξ dη dϕ N´otese que la matriz elemental puede escribirse como {Ae}ij ={Me}ij +{Ne}ij (3.3)
GENERACI ´ ON DE MATRICES VARIABLES 38 El ensamblaje de tales matrices elementales conduce a un sistema lineal de la forma: Ax=b(3.4) donde Aes una matriz sim´etrica variable del tipo A=M+ N (3.5) siendo M y N dos matrices, tipo ”sparse”, diferentes, Sim´etricas Definidas Positivas,(SDP), constantes para un nivel de discretizaci´on dado y el llamado par´ametro de estabilidad del modelo, ya mencionado; debiendo resolverse el sistema para cada valor diferente de . Precisamente el objetivo principal de esta tesis consiste en plantear propuestas que hagan posible solucionar por m´etodos iterativos, no demasiado costosos, este tipo de sistemas de ecuaciones lineales de matrices variables.
Cap´ıtulo 4 ESTIMACI´ ON DE PAR´ AMETROS 4.1. CONSIDERACIONES PREVIAS La eficiencia de los modelos de masa consistente para ajuste de campos de viento depende en gran medida de ciertos par´ametros que aparecen en las distintas etapas del proceso, especialmente de algunos de los que intervienen en la construcci´on del campo de viento inicial y de los m´odulos de precisi´on de Gauss. En general, los valores de estos par´ametros se toman usando una serie de reglas emp´ıricas. Se puede plantear su estimaci´on de manera autom´atica, tal que las velocidades observadas en las estaciones de medida sean regeneradas de la forma m´as exacta posible por el modelo, dando lugar por tanto a un problema inverso. Existen diversos m´etodos de resoluci´on de problemas inversos relacionados con la estimaci´on de par´ametros. De entre ellos, se han elegido los algoritmos gen´eticos, por ser una herramienta robusta y flexible, que puede ser competitiva ya que los c´alculos pueden paralelizarse. Para el c´alculo autom´atico de ciertos par´ametros del modelo de viento se plantea el siguiente problema inverso: de las Nestaciones de medida disponibles se toman Nrcomo referencia; el resto, se utilizan para el c´alculo del viento. El viento as´ı obtenido se compara con el medido en las Nrestaciones de referencia. Para la estimaci´on de los par´ametros del modelo se procede a minimizar la diferencia entre los resultados obtenidos y las medidas observadas en las estaciones de referencia.
CONSIDERACIONES PREVIAS 40 Esta t´ecnica supone una mejora sustancial sobre la propuesta de Barnard [8] para la estimaci´on de uno solo de los par´ametros por un procedimiento de ensayo y error. E. Rodr´ıguez, en su Tesis: Modelizaci´on y simulaci´on num´erica de campos de viento mediante elementos finitos adaptativos en 3-D, [89], propone extenderlo a un total de cuatro de los par´ametros que intervienen en el modelo y automatizar su c´alculo; de tal forma que la funci´on a minimizar sea: F(α, β, γ, γ0) = 1 Nr Nr X n=1 |~vn−~v(xn, yn, zn)| |~vn| donde ~v(xn, yn, zn) es la velocidad del viento obtenida por el modelo en la posici´on de la estaci´on n, y Nres el n´umero de estaciones de referencia, 1 ≤Nr≤N, siendo Nel n´umero total de estaciones de medida disponibles. En primer lugar se considera el par´ametro de estabilidad α=α1 α2 =rTv Th =√ ya establecido como f´ormula 2.12, que se deriva del funcional 2.8, y cuyo m´ınimo no var´ıa si se divide por α2 2. Hay que se˜nalar que para α >> 1 predomina el ajuste de viento en la direcci´on vertical, mientras que para α << 1 el ajuste tiene lugar predominantemente sobre el plano horizontal. Por lo tanto la elecci´on de α, ´o , determina que el viento tienda a rodear los obst´aculos o a sobrepasarlos. Diversos experimentos num´ericos han hecho patente que el comportamiento de los modelos de masa consistente depende sensiblemente de la elecci´on de los valores de , por lo que se presta particular atenci´on a este problema. Diversos autores han estudiado c´omo parametrizar la estabilidad debido a que la dificultad en la determinaci´on de los valores de αhan limitado el uso de modelos de masa consistente en terrenos de orograf´ıa compleja. En [103, 54, 14], los autores proponen tomar α= 10−2, o sea, proporcional a la magnitud de w/u. Otros, como Ross [91] y Moussiopoulos [76], relacionan αcon el n´umero de Froude, mientras que Geai, [38], Lalas, [58], y Tombrou, [108], proponen que el par´ametro αvar´ıe en la direcci´on vertical. Finalmente, Barnard et al., [8], proponen un procedimiento para obtener αen cada simulaci´on del campo de viento. La idea es usar Nvelocidades de viento observadas para obtener el campo de viento y usar las restantes Nrcomo referencia. Entonces se realizar´ıan diversas simulaciones con distintos valores de , lo que supondr´ıa
CONSIDERACIONES PREVIAS 41 resolver la ecuaci´on (3.4) Ax=b varias veces, confirmando la idea de que es muy interesante poder resolver este tipo de ecuaciones, de matrices variables, de la forma m´as eficiente posible. El valor que m´as acerque el viento estimado al observado en las estaciones de referencia es el que se toma como valor del par´ametro de estabilidad. Este m´etodo proporciona valores de que s´olo son v´alidos para cada caso particular y, por tanto, no proporciona valores v´alidos a priori para otras simulaciones. E. Rodriguez, en su Tesis, [89], estudia una versi´on del m´etodo propuesto por Barnard et al., [8], utilizando algoritmos gen´eticos como herramienta de optimizaci´on que permite una selecci´on autom´atica de . El segundo par´ametro a estimar es el coeficiente de peso β(0 ≤β≤1) de la ecuaci´on 2.18, correspondiente a la interpolaci´on horizontal de las medidas de viento observadas. Cuando β→1 adquiere m´as importancia la distancia horizontal de cada punto a las estaciones de medida, mientras que para β→0 se da m´as peso a la distancia vertical entre cada punto y las estaciones [70]. En general, para terrenos complejos se utiliza la segunda aproximaci´on [79]. En orograf´ıas m´as llanas o en an´alisis horizontales en 2-D, se utiliza la primera. En aplicaciones m´as realistas existir´an zonas con orograf´ıa compleja y zonas de orograf´ıa m´as regular, lo que sugiere el uso de valores intermedios de β. El siguiente par´ametro objeto de estimaci´on es γ, que aparece en la ecuaci´on 2.19 y est´a relacionado con la capa l´ımite planetaria en la estratificaci´on atmosf´erica. Existen diferentes autores que proponen distintos rangos para este par´ametro. Panofsky y Dutton, [80], proponen el intervalo [0.15, 0.25]. Sin embargo, en [85] se utiliza directamente el valor γ= 0,3 en el c´odigo de su programa WINDS, mientras que para Baas, [7], γha de estar dentro del intervalo [0.3, 0.4]. En nuestras simulaciones el espacio de b´usqueda de γincluye todas estas posibilidades. Finalmente, tambi´en resulta de inter´es obtener estimaciones de los valores del par´ametro γ0, que interviene en el c´alculo de la altura de la capa de mezcla en el caso de condiciones atmosf´ericas estables, v´ease la f´ormula (2.20). Garrat propone directamente γ0= 0,4. Tambi´en en el c´odigo de WINDS el valor de γ0est´a en
ALGORITMOS GENETICOS 42 torno a 0,4. As´ı, hemos definido el espacio de b´usqueda para el valor de γ0en el entorno de 0,4. Desde el punto de vista de esta tesis, donde nos planteamos el precondicionamiento de los sistemas de ecuaciones variables de los modelos de campos de viento (3.4): Ax=b conviene aclarar que los par´ametros β,γyγ0solo intervienen en la interpolaci´on del viento inicial y por tanto sus valores afectan, ´unicamente, al c´alculo del vector del segundo miembro b, mientras que el par´ametro de estabilidad α, ´o , afecta al c´alculo del viento resultante, ya que con ´el cambia la estructura de la matriz Apuesto que viene expresada como (3.5) A=M+ N A la hora de recurrir al precondicionamiento del sistema para mejorar su resoluci´on por m´etodos iterativos los valores asignados al par´ametro (α) tienen un gran protagonismo, de ah´ı que se le dedique especial atenci´on a su estimaci´on ´optima. Para ello se propone la utilizaci´on de Algoritmos Gen´eticos, como herramienta robusta, para su selecci´on. En el Cap´ıtulo 9.2 se presentan los resultados de numerosos experimentos num´ericos, con los que se comprueba el comportamiento de los distintos precondicionadores propuestos, para una amplia gama de valores de α() 4.2. ALGORITMOS GENETICOS Los algoritmos gen´eticos (en los sucesivo AG) son herramientas de optimizaci´on basadas en el mecanismo de evoluci´on natural. Producen intentos sucesivos que tienen una probabilidad cada vez mayor de alcanzar el ´optimo global. Los aspectos m´as importantes de los AG son la construcci´on de una poblaci´on inicial, la evaluaci´on de cada individuo a trav´es de la funci´on de aptitud o funci´on objetivo, la selecci´on de los padres de la siguiente generaci´on, el cruce de esos padres para crear los hijos y la mutaci´on, que incrementa la diversidad. En la figura 4.1 puede verse una representaci´on esquem´atica del funcionamiento de los AG. Se parte de una poblaci´on inicial a la que se somete a prueba a trav´es
SUBESPACIOS DE KRYLOV 48 En los m´etodos iterativos , a partir de un vector inicial, x0, se genera una secuencia de vectores (xi), que converge a la soluci´on buscada. A pesar de la convergencia relativamente lenta, el car´acter sparse de Ahace posible efectuar un elevado n´umero de iteraciones sin un trabajo excesivo. Por otra parte, los errores de redondeo, importantes en los m´etodos directos, influyen, por lo general, en la velocidad de convergencia, pero no en la aproximaci´on final. Adem´as de los m´etodos cl´asicos de este tipo, (Jacobi, Gauss-Seidel, SOR y SSOR) que se ajustan a estas precisiones, se han desarrollado otros que presentan mayores ventajas respecto a los mismos, sobre todo en estos sistemas sparse. Entre ellos han adquirido ´ultimamente especial relevancia los algoritmos basados en los Subespacios de Krylov [94]. La r´apida evoluci´on experimentada por los sistemas inform´aticos ha contribuido a facilitar su implementaci´on [41]. 5.2. SUBESPACIOS DE KRYLOV En los m´etodos iterativos, los sucesivos valores de la soluci´on aproximada en la resoluci´on de A x =b, vienen dados por la relaci´on de recurrencia xi+1 =xi+B−1(b−A xi) (5.1) ´o bien B xi+1 =C xi+b, siendo C =B−A(5.2) Diferentes elecciones para las matrices ByC, en funci´on de la matriz A, conducen a los m´etodos cl´asicos de relajaci´on (Jacobi, Gauss-Seidel, SOR). Pero estas expresiones tambi´en se pueden escribir en funci´on del vector residuo, ri=b−A xi con lo cual xi+1 =xi+B−1r1 De esta forma, eligiendo una aproximaci´on inicial, x0, los valores de las sucesivas iteraciones se podr´ıan calcular por las respectivas expresiones: x1=x0+B−1r0
SUBESPACIOS DE KRYLOV 49 x2=x1+B−1r1=x0+B−1r0+B−1(b−A x1) = =x0+B−1r0+B−1(b−A xo−A B−1r0) = x0+B−1r0+B−1(r0−A B−1r0) = =x0+ 2B−1r0−B−1A B−1r0 x3=x2+B−1r2 x4=x3+B−1r3 ......................... ......................... xi=xi−1+B−1ri−1 En el caso de que B=Iquedar´ıa: x1=x0+r0 x2=x0+ 2 r0−A r0 x3=x2+r2=x0+ 2 r0−A r0+ (b−A x2) = =x0+ 2 r0−A r0+ (b−A x0+ 2 A r0−A2r0) = =x0+ 2 r0−A r0+ (r0+ 2 A r0−A2r0) = =x0+ 3 r0+A r0−A2r0 x4=x3+r3 ......................... ......................... xi=xi−1+ri−1 Es decir, la i-´esima iteraci´on de la soluci´on aproximada se puede expresar como la suma de la aproximaci´on inicial y una combinaci´on lineal de ivectores: xi=x0+C.L.{r0, Ar0, A2r0, ......, Ai−1r0} Por tanto xi=x0+ [r0, Ar0, A2r0, ......, Ai−1r0] (5.3) El subespacio Ki(A;r0), de base [r0, Ar0, A2r0, ......, Ai−1r0], es llamado subespacio de Krylov de dimensi´on i, correspondiente a la matriz Ay residuo inicial
M´ ETODO DEL GRADIENTE 50 r0 5.3. M´ ETODO DEL GRADIENTE En los sistemas de ecuaciones lineales A x =b, cuya matriz de coeficientes es Sim´etrica y Definida Positiva (SDP), el m´etodo del Gradiente ´o del M´aximo Descenso es una buena herramienta para su resoluci´on. Sea A∈ <n×nuna matriz SDP y b∈ <n. Se considera como soluci´on ´optima x∈ <ndel sistema A x =b, aquella que minimiza la funci´on de error: E(x) = A e(x), e(x)∈ < /E(x) : <n→ < (5.4) siendo e(x), el error de una soluci´on, = x−¯x, ∈ <n, en la que ¯xes la soluci´on exacta del sistema y r(x), el residuo, = b−A x =A¯x−A x =A(¯x−x). Minimizar E(x) = A(x−¯x), x−¯x=A x−A¯x, x−¯x=A x, x−2A¯x, x+ A¯x, ¯x, implica, puesto que A¯x, ¯xes constante, obtener el m´ınimo para la funci´on A x, x−2A¯x, x= 2J(x), siendo J(x) = 1 2A x, x−b, x/J(x) : <n→ < (5.5) No cabe duda que encontrar una soluci´on x, lo m´as cercana posible a ¯xpasa por minimizar J(x). La condici´on necesaria para que una funci´on de varias variables, diferenciable, alcance un valor m´ınimo en x, es que su gradiente sea cero: grad J(x) = J0(x) = 1 2A x, x0 −b, x0 =1 2{[(A x)0]Tx+(x0)TA x}−(x0)Tb= =1 2(ATx+A x)−b=1 2(AT+A)x−b= 0 Si adem´as Aes sim´etrica: J0(x) = A x −b(5.6) La naturaleza del punto cr´ıtico depender´a del signo de J00 (x) = A Si Aes definida positiva, no cabe duda que la soluci´on del gradiente nulo se corresponder´a con un m´ınimo del error, y, en este caso, con la soluci´on exacta. La ventaja de utilizar J(x), es que el problema de optimizar x∈ <nes reemplazado por un problema unidimensional que se puede describir como sigue:
M´ ETODO DEL GRADIENTE 51 I. Evaluar Jcon una aproximaci´on inicial x0. II. Determinar una direcci´on a partir de x0que produzca una descenso en J. III. Calcular cuanto debemos movernos en esa direcci´on para obtener una mejor soluci´on x1. IV. Volver al paso I, reemplazando x0por x1. Proceso de iteraci´on, que escribiremos como: xi+1 =xi+αipi(5.7) donde pies un vector direcci´on y αi, un escalar, que determinaremos minimizando J(xi) a lo largo de pi: J(xi+αipi) = minαJ(xi+α pi) con xi, pi∈ <nfijos: f(α) = J(xi+α pi) = 1 2(xi+α pi)TA(xi+α pi)−bT(xi+α pi) = =1 2α2pT iA pi+αpT i(A xi−b) + 1 2xT i(A xi−2b) El punto cr´ıtico de la funci´on parab´olica f(α) se puede determinar haciendo f0(α) = 0: f0(α) = αpT iA pi+pT i(A xi−b) = 0 y teniendo en cuenta que f00 (α) = pT iA pi, siendo Auna matriz definida positiva, se corresponder´a con un m´ınimo, por lo cual αi=αopt(xipi) = pT i(b−A xi) pT iA pi =ri, pi A pi, pi(5.8) El conocimiento del factor αinos permite tambi´en establecer, con facilidad, una formulaci´on iterativa para el vector residuo: As´ı, partiendo de la expresi´on 5.7, podemos transformarla con las operaciones siguientes: A xi+1 =A xi+αiA pi b−A xi+1 =b−A xi−αiA pi ri+1 =ri−αiA pi(5.9)
M´ ETODO DEL GRADIENTE 52 Formulaci´on que nos va a permitir demostrar que cualquier vector residuo resulta siempre ortogonal a la direcci´on de descenso anterior. En efecto, recurriendo al producto escalar de ambos vectores, tendremos: ri+1, pi=ri, pi−αiA pi, pi= =ri, pi−ri, pi A pi, piA pi, pi= 0 resultado nulo que confirma su ortogonalidad. Sabemos que el gradiente de una funci´on, en un punto, define la direcci´on de la derivada direccional m´axima y adem´as tiene el sentido en que aumenta la funci´on; por tanto, parece l´ogico, que para minimizar la funci´on J(x) tomemos como direcci´on de descenso, pi, la del vector gradiente, pero en sentido contrario, o sea pi=−∇J(xi) = −J0(xi) y que seg´un 5.6, se convertir´a en pi=b−A xi=ri por lo que la f´ormula de iteraci´on se convertir´a en: xi+1 =xi+αiri(5.10) y la formulaci´on iterativa para el vector residuo se transformar´a en: ri+1 =b−A xi+1 =b−A(xi+αiri) = ri−αiA ri(5.11) dando lugar al algoritmo siguiente [94]: ALGORITMO DEL GRADIENTE Valor inicial: x0 r0=b−A x0 para i= 0,1,2, .... hasta la convergencia, hacer: αi=ri,ri A ri,ri xi+1 =xi+αiri ri+1 =ri−αiA ri
M´ ETODO DEL GRADIENTE CONJUGADO 53 5.4. M´ ETODO DEL GRADIENTE CONJUGADO Hemos visto en el m´etodo anterior, que la condici´on J(x)≤J(x+α p),∀α∈ < supone conseguir el valor ´optimo de xen la direcci´on p, lo que implica que cada nuevo vector residuo es ortogonal a la direcci´on de descenso anterior y, a su vez, dicha direcci´on de descenso viene dada por un vector opuesto al gradiente de J(x), que coincide con el vector residuo, con lo que resultar´a: ri+1⊥piy pi≡ri⇒ri+1⊥ri La dificultad del M´etodo del Gradiente reside en que esta relaci´on de ortogonalidad no es transitiva, es decir si bien ri+1⊥riyri+2⊥ri+1, esto no supone que ri+2⊥ri y, por consiguiente, con las sucesivas iteraciones se vaya perdiendo la condici´on de optimizaci´on de x. Para mantener esta condici´on, se debe trabajar con unas direcciones de descenso, que a diferencia del M´etodo del Gradiente, conserven los requisitos de ortogonalidad con respecto a todas las direcciones anteriores y que por tanto todas las direcciones sean conjugadas. Seg´un hemos visto anteriormente, si x0es un valor ´optimo de xen la direcci´on presulta que: x0=x+α p ⇒r0⊥p Hallemos ahora un nuevo valor x00 en una nueva direcci´on q: x00 =x0+α q el nuevo residuo ser´a r00 =b−Ax00 =b−A(x0+α q) = b−Ax0−α A q =r0−α A q para que el nuevo valor r00 siga siendo ´optimo respecto a la direcci´on de pdeber´a cumplirse que r00 ⊥p, es decir: r0−α A q, p= 0
M´ ETODO DEL GRADIENTE CONJUGADO 54 o sea r0, p−αA q, p= 0 como r0⊥p⇒r0, p= 0 ⇒A q, p= 0 y, por tanto, se dice que los vectores pyq, que cumplen dicha condici´on, son A conjugados; como adem´as la matriz Aes sim´etrica y definida positiva, podemos afirmar que pyqson A-ortogonales. En lo sucesivo usaremos direcciones de descenso p0, p1, p2, ......pique sean Aortogonales dos a dos, es decir, tal que: pT i·A pj= 0,∀i6=j . Adem´as, teniendo en cuenta que la matriz Aes sim´etrica y definida positiva, podemos comprobar que el conjunto de los vectores p0, p1, p2, ......pi,A-ortogonales, resultan ser linealmente independientes. En efecto, si existe una colecci´on de coeficientes escalares αk, tal que: α0p0+α1p1+α2p2+...... +αipi= 0 se cumplir´a tambi´en que: α0A p0+α1A p1+α2A p2+...... +αiA pi= 0 que multiplicando escalarmente por p0, se convierte en α0pT 0A p0+α1pT 0A p1+α2pT 0A p2+...... +αipT 0A pi= 0 como por hip´otesis todos los pT i·A pj,∀i6=j, son nulos, resultar´a que α0pT 0A p0= 0 y al ser Adefinida positiva implicar´a que: α0= 0 De forma similar se puede demostrar para cualquier αk/ k = 1,2, .....i ⇒ ∀k∈ {0,1,2, ....i}, αk= 0 y por tanto todos los vectores A-ortogonales ser´an linealmente independientes.
M´ ETODO DEL GRADIENTE CONJUGADO 55 El M´etodo del Gradiente Conjugado [50, 52, 1, 84] viene a ser una variante del M´etodo del Gradiente, en el que las sucesivas direcciones de descenso se generan como versiones conjugadas de los gradientes que se van obteniendo seg´un progresa el m´etodo. Cada nueva direcci´on,pi+1, se obtiene en el plano formado por las direcciones ortogonales ri+1 ypisiguiendo la expresi´on pi+1 =ri+1 +βipi(5.12) determin´andose el escalar βide tal forma que pi+1 ypisean A-ortogonales, es decir: pT i+1 A pi= 0 por tanto (ri+1 +βipi)TA pi= (rT i+1 +βipT i)A pi= 0 resultando βi=−rT i+1 A pi pT iA pi =−A pi, ri+1 A pi, pi(5.13) lo que nos va a garantizar que los residuos sucesivos sean ortogonales, o sea que ri+1, ri= 0 . En efecto: ri+1 =ri−αiA pi por tanto ri+1, ri=ri−αiA pi, ri=ri, ri−αiA pi, ri= =ri, ri−αiA pi, pi−βi−1pi−1=ri, ri−αiA pi, pi+αiβi−1A pi, pi−1 pero teniendo en cuenta que piypi−1son A-ortogonales, y el valor de αi, resultar´a que
M´ ETODO DEL GRADIENTE CONJUGADO 56 ri+1, ri=ri, ri−ri, pi =ri, ri−ri, ri+βi−1pi−1 =ri, ri−ri, ri−βi−1ri, pi−1 = 0 Como consecuencia de todo lo anterior, se puede afirmar que este M´etodo del Gradiente Conjugado tiene dos propiedades esenciales: A) Cada direcci´on de descenso es A-ortogonal a todas las direcciones anteriores, con lo cual no se pierde la condici´on de ´optimo de cada nuevo vector, x, hallado: Partiendo de cada producto A pi+1, pi= 0, se puede comprobar f´acilmente que A pi+1, pk= 0, para 0≤k≤i. B) Prescindiendo de los errores de redondeo, el M´etodo converge a lo sumo en niteraciones [112], siendo nel orden de la matriz cuadrada A: En efecto, tomando como base que ri+1, pi= 0, se puede demostrar con facilidad que ri+1, pk= 0, para 0≤k≤i. o sea, que cada vector residuo es ortogonal a todas las direcciones de descenso anteriores. Ahora bien, para k < n −1, puede ocurrir que rk= 0, con lo cual b−A xk= 0, habr´ıamos encontrado ya la soluci´on correcta y el M´etodo converger´ıa en k iteraciones. Pero mientras rk6= 0, habr´a un residuo no nulo, cada vez m´as peque˜no. As´ı llegar´ıamos hasta rn, que ser´ıa ortogonal a p0, p1, p2, ......pn−1, que son nvectores linealmente independientes y A-ortogonales , con lo cual rntendr´ıa que ser nulo y por tanto b−A xn= 0, con lo que habr´ıamos llegado a la soluci´on exacta, y por tanto a la ´ultima y n-sima iteraci´on. Precisamente los nvectores, p0, p1, p2, ......pn−1, linealmente independientes, definen un subespacio de ndimensiones, A-ortogonales, que resulta ser un Subespacio de Krylov, correspondiente a la matriz A=<n×n, de vector inicial p0, que se hace corresponder con el residuo inicial r0.
M´ ETODO DEL GRADIENTE CONJUGADO 57 Hay que se˜nalar, sin embargo, que en la pr´actica, la aparici´on de errores de redondeo, hace que las direcciones de descenso no sean exactamente A-ortogonales y, por tanto, que el Gradiente Conjugado se comporte como un m´etodo iterativo cualesquiera [4]. El esquema principal del algoritmo del Gradiente Conjugado consiste en obtener xi+1 =xi+αipi sabiendo, por la expresi´on (5.8), que αi=ri, pi A pi, pi es el escalar ´optimo, que minimiza la funci´on error, y tomando cada direcci´on de descenso, seg´un (5.12), como pi+1 =ri+1 +βipi siendo βi=−A pi, ri+1 A pi, pi que nos permite mantener la ortogonalidad entre todas las direcciones ´optimas, pi. Teniendo, adem´as, en cuenta la expresi´on (5.9), de iteraci´on de los residuos: ri+1 =ri−αiA pi Las condiciones de ortogonalidad del proceso, entre vectores residuo y direcciones de descenso, nos van a permitir expresar αiyβide forma m´as sencilla, de tal manera que hagan este algoritmo m´as eficiente: As´ı αi=rT ipi pT iA pi =rT i(ri+βi−1pi−1) pT iA pi =rT iri+βi−1rT ipi−1 pT iA pi =rT iri pT iA pi =ri, ri A pi, pi
OTROS M´ ETODOS DE KRYLOV 64 problema de m´ınimos cuadrados que plantean estos m´etodos; es decir, haciendo una factorizaci´on LU, en lugar de la factorizaci´on QR tradicional [36]. 5.5.3. M´ ETODOS BASADOS EN LA ECUACI´ ON NORMAL La resoluci´on del sistema de ecuaciones A x =b, donde la matriz Aes no sim´etrica, es equivalente a resolver el sistema ATA x =ATbde matriz ATA, sim´etrica definida positiva, al que se puede aplicar el algoritmo del Gradiente Conjugado. La ecuaci´on ATA x =ATb recibe el nombre de Ecuaci´on Normal. Al igual que el CG, los m´etodos desarrollados a partir de la Ecuaci´on Normal cumplen las dos condiciones fundamentales de minimizaci´on de la norma residual y optimizaci´on del coste computacional, si embargo presentan el inconveniente de que el condicionamiento del nuevo sistema es el cuadrado del sistema inicial: K(ATA) = K(A)2, lo cual, para sistemas mal condicionados, puede resultar desastroso y adem´as en cada iteraci´on aparecen dos productos matriz por vector correspondientes a las matrices AyATaumentando el coste computacional. Resultan as´ı m´etodos como: El M´etodo CGN(M´etodo del Gradiente Conjugado para la Ecuaci´on Normal). Este m´etodo construye una sucesi´on de vectores: xk=x0+hATr0,(ATA)ATr0,(ATA)2ATr0, ....., (ATA)k−1ATr0i con residuo m´ınimo en cada paso, sin efectuar el c´alculo expl´ıcito del producto ATA. El M´etodo LSQR(Least-Square QR). Propuesto por Paige y Saunders, en 1982 [78]. Este m´etodo intenta corregir el posible empeoramiento del n´umero de condici´on de los sistemas al aplicar el m´etodo de la Ecuaci´on Normal.
OTROS M´ ETODOS DE KRYLOV 65 La idea b´asica del LSQR es hallar la soluci´on del sistema sim´etrico: I A AT−λ2I r x = b 0 minimizando, A λI x− b 0 2 donde λes un n´umero real arbitrario.
Cap´ıtulo 6 PRECONDICIONAMIENTO 6.1. CONSIDERACIONES PREVIAS Aunque los m´etodos iterativos basados en los subespacios de Krylov est´an, te´oricamente, bien fundamentados, casi todos ellos adolecen de lentitud en la convergencia, principalmente aquellos que se usan para resolver problemas que surgen de temas como din´amica de fluidos o simulaciones de mecanismos electr´onicos. Precondicionar un sistema es un paso clave para el ´exito de los m´etodos de Krylov que se usan en estas aplicaciones [74]. Es ampliamente reconocido que la falta de robustez es una de las debilidades de los m´etodos iterativos, este inconveniente ha impedido la amplia aceptaci´on de estos m´etodos en aplicaciones industriales, a pesar de su intr´ınseco atractivo para resolver grandes sistemas de ecuaciones lineales. Tanto la eficacia, como la robustez de las t´ecnicas iterativas pueden mejorarse con el uso de los precondicionadores. Precondicionar es simplemente un medio de transformar un sistema lineal original en otro que tenga la misma soluci´on, pero que sea m´as f´acil de resolver por m´etodos iterativos. En general la fiabilidad de estos m´etodos dependen mucho m´as de la calidad del precondicionador, que del m´etodo de Krylov elegido. Encontrar un buen precondicionador para resolver un sistema lineal sparse est´a considerado, con frecuencia, como una combinaci´on de arte y ciencia [94]. algunos m´etodos de precondicionamiento funcionan sorpresivamente bien, a pesar de sus escasas expectativas te´oricas. N´otese que en principio no hay virtualmente l´ımites para elegir opciones que
CONDICIONAMIENTO DE UN SISTEMA 67 permitan obtener buenos precondicionadores. Por ejemplo, los precondicionadores pueden derivarse del conocimiento de los problemas f´ısicos de los que surge el sistema lineal o pueden construirse a partir de la matriz de coeficientes del sistema original. En l´ıneas generales un precondicionador es cualquier forma expl´ıcita o impl´ıcita de modificaci´on de un sistema lineal original que lo haga m´as f´acil de resolver por un m´etodo iterativo dado. El sistema resultante debe poderse resolver por un m´etodo basado en los subespacios de Krylov y deben requerirse menos pasos para su convergencia, que si se aplicara el mismo m´etodo al sistema original (Aunque esto no pueda garantizarse te´oricamente, sino confirmarse con la experiencia). En definitiva, para resolver ciertos problemas es indispensable la implementaci´on de un precondicionador adecuado para asegurar la convergencia del m´etodo de Krylov elegido. En esta secci´on introduciremos ideas generales sobre precondicionamiento de sistemas y relacionaremos algunos de los precondicionadores m´as utilizados para resolver sistemas de matriz sim´etrica definida positiva, que son los que surgen en los modelos de campos de viento. 6.2. CONDICIONAMIENTO DE UN SISTEMA Decimos que un sistema de ecuaciones, A x =b, est´a bien o mal condicionado, cuando peque˜nas variaciones en sus coeficientes, o en sus t´erminos independientes, producen una peque˜na ´o gran variaci´on en la soluci´on del mismo. Con objeto de dar una medida del buen o mal condicionamiento de un sistema se introduce la noci´on de n´umero de condici´on (relativo a una norma matricial dada) [15, 32]. Limit´andonos al caso m´as sencillo, si se produce una perturbaci´on, ∆b, en el vector columna de los t´erminos independientes, vamos a analizar la perturbaci´on, ∆x, que se producir´a en la soluci´on exacta, x, del sistema. Empleando una norma matricial arbitraria y una norma vectorial cualquiera, compatible con ella, se tendr´a: A x =b A(x+ ∆x) = b+ ∆b ⇒A∆x= ∆b
CONDICIONAMIENTO DE UN SISTEMA 68 ∆x=A−1∆b b=A x ⇒ k∆xk≤k A−1k k ∆bk kbk≤k Ak k xk ⇒ k∆xk k bk≤k Ak k A−1k k ∆bk k xk ⇒ ⇒k∆xk kxk≤k Ak k A−1kk∆bk kbk Al factor kAk k A−1k=K(A), se le denomina n´umero de condicionamiento del sistema, relativo a la norma matricial k·k; obteni´endose as´ı la siguiente relaci´on del error relativo de la soluci´on, en funci´on del error relativo en los t´erminos independientes: k∆xk kxk≤K(A)k∆bk kbk(6.1) Evidentemente, cuanto menor sea el valor num´erico, K(A), menor ser´a la variaci´on de la soluci´on ante las posibles fluctuaciones en los valores de los elementos de b. Dicho n´umero ser´a siempre mayor o igual que la unidad: 1≤k Ik=kA·A−1k≤k Ak k A−1k=K(A) (6.2) Cuanto m´as pr´oximo a la unidad sea K(A), mejor condicionado estar´a el sistema: l´ım K(A)→1k∆xk/kxk k∆bk/kbk= 1 Conviene resaltar que el n´umero de condici´on, K(A), del sistema A x =b, es una cantidad intr´ınseca a su matriz de coeficientes, A. En otras palabras, el condicionamiento del sistema es independiente del vector b, de t´erminos independientes. De forma similar, si denominamos ∆Aa las variaciones introducidas en los coeficientes de la matriz, no cabe duda que el sistema se transformar´a en (A+ ∆A)(x+ ∆x) = b A partir del cual, como en el caso anterior, se puede comprobar f´acilmente que: k∆xk kx+ ∆xk≤k Ak k A−1kk∆Ak kAk lo que nos indica que el error relativo de la nueva soluci´on es tambi´en funci´on del error relativo producido por las variaciones de la matriz de coeficientes y, ambos, est´an relacionados por el mismo n´umero K(A), ya definido anteriormente.
T´ ECNICAS DE PRECONDICIONAMIENTO 69 Por tanto, cuanto m´as cercano a 1 est´e el n´umero de condicionamiento, K(A), tanto menor ser´a la variaci´on de la soluci´on del sistema, ante las posibles fluctuaciones en los valores, tanto de los elementos de A, como de by por tanto se dir´a que el sistema est´a mejor condicionado. Usando la norma espectral k · k2, K(A) = µM µm≥1, donde µMyµmson, respectivamente, los valores singulares m´aximo y m´ınimo de la matriz del sistema, valores que se pueden calcular por µi=pλi(AAT), siendo λi(AAT) los correspondientes valores propios de AAT. Cuando Aes sim´etrica,µi(A)≡λi(A), y el n´umero de condicionamiento ser´ıa K(A) = λM λm≥1 Con esta definici´on, dado que la bondad del condicionamiento de una matriz viene dada, en principio, por la proximidad de K(A) a la unidad, en el caso que los valores singulares extremos coincidieran, el n´umero de condicionamiento adquirir´ıa este valor ´optimo. 6.3. T´ ECNICAS DE PRECONDICIONAMIENTO En numerosas ocasiones, la matriz y el vector segundo miembro del sistema se calculan de forma aproximada, pudiendo existir ciertas diferencias con los valores num´ericos que reflejen exactamente el problema. En estos casos, un mal condicionamiento del sistema afectar´ıa negativamente a la convergencia. Se hace necesario as´ı, mejorar este condicionamiento utilizando adecuadas t´ecnicas de precondicionamiento. T´ecnicas que, en general, consisten en transformar el sistema en otro de id´entica soluci´on, pero con menor K(A). Para ello multiplicaremos la expresi´on A x =bpor una matriz M, llamada matriz de precondicionamiento: M A x =M b tal que K(M A)< K(A).
T´ ECNICAS DE PRECONDICIONAMIENTO 70 El menor valor de K(M A) corresponder´ıa a M=A−1, puesto que quedar´ıa K(A−1A) = 1, que es, obviamente, el caso ideal y el sistema converger´ıa en una sola iteraci´on, pero el coste computacional del c´alculo de A−1equivaldr´ıa a resolver el sistema por un m´etodo directo. Esta circunstancia sugiere para Muna matriz lo m´as pr´oxima posible a A−1, sin que su determinaci´on suponga un elevado coste. Generalmente, se opta por considerar como matriz de precondicionamiento a M−1, y obtener Mcomo aproximaci´on de A, escribiendo, M−1A x =M−1b Los ordenadores, con m´ultiples procesadores en paralelo, ofrecen una gran versatilidad [92]. Grote y Simon [47] proponen obtener Mcomo aproximaci´on de A−1, minimizando una norma que reduce el c´alculo a npeque˜nos problemas independientes de m´ınimos cuadrados. Asimismo, en [102] y [77] se plantea efectuar esta aproximaci´on por una expresi´on polin´omica P(A). El campo de posibles precondicionadores aplicables en ordenadores que no tengan estas caracter´ısticas es tambi´en muy amplio. Por otro lado, en los algoritmos precondicionados de los distintos m´etodos figurar´an productos de matriz inversa por vector, que no deben exigir excesivo trabajo adicional, por ello, la matriz Mdebe ser f´acilmente invertible. Por ejemplo una matriz diagonal, o una matriz factorizada adecuadamente para efectuar esos productos por procesos de remonte, sin necesidad de calcular M−1. Dependiendo de la forma de plantear el producto de la inversa de la matriz de precondicionamiento por la matriz del sistema, y aprovechando la descomposici´on en factores de aquella, se distinguen los siguientes casos: -a) Precondicionamiento por la izquierda: M−1A x =M−1b M−1A=˜ A M−1b=˜ b ˜ A x =˜ b -b) Precondicionamiento por la derecha: A M−1M x =b A M−1=˜ A M x = ˜x ˜ A˜x=b -c) Precondicionamiento por ambos lados:
T´ ECNICAS DE PRECONDICIONAMIENTO 71 Expresando Mfactorizada como M=M1M2, M−1 1AM−1 2M2x=M−1 1b M−1 1AM−1 2=˜ A M2x= ˜x M−1 1b=˜ b ˜ A˜x=˜ b Las caracter´ısticas de cada problema y, en definitiva, de la matriz A, del precondicionador utilizado Me, incluso, de la tolerancia exigida para el criterio de parada adoptado, hacen m´as eficiente una forma u otra de precondicionamiento, sin que pueda establecerse a priori bases de elecci´on que nos inclinen por una de ellas. El precondicionamiento de un sistema tiene por finalidad mejorar la convergencia del m´etodo aplicado respecto a la convergencia del sistema sin precondicionar. En la resoluci´on de sistemas sim´etricos, la raz´on de convergencia del Gradiente Conjugado kx−xikA≤2 pK(A)−1 pK(A)+1!i kx−x0kA depende del n´umero de condicionamiento K(A), funci´on de los autovalores mayor y menor de la matriz del sistema. Sin embargo, en la pr´actica, despu´es de cierto n´umero de iteraciones la convergencia se hace superlineal, como si el n´umero de condicionamiento inicial fuese sustituido por otro menor, de tal forma que la raz´on de convergencia depende, adem´as, de la distribuci´on total de los valores propios de A. Las t´ecnicas de precondicionamiento, tienen por objeto transformar el sistema original A x =b, en otro, con otra nueva matriz ˜ A, con unos nuevos autovalores que conduzcan a una convergencia m´as r´apida, bien disminuyendo el n´umero de condicionamiento K(A), bien mejorando la distribuci´on de los autovalores m´as peque˜nos del espectro [48, 113]. Las distintas iteraciones xi, del sistema precondicionado verificar´an: kx−xik˜ A≤2 qK(˜ A)−1 qK(˜ A)+1 i kx−x0k˜ A con xi∈x0+Ki(˜ A; ˜r0). En los sistemas no sim´etricos, es complicado probar que el sistema precondicionado resultante posee un espectro de autovalores que mejore la convergencia
M´ ETODO DEL GRADIENTE CONJUGADO PRECONDICIONADO 72 respecto al sistema original. Pero, por analog´ıa, y, a´un sin demostraciones matem´aticas que lo confirmen, se podr´ıa esperar un comportamiento similar. Esta conclusi´on, corroborada num´ericamente en todas las aplicaciones, permite tratar t´ecnicas de precondicionamiento en los m´etodos tipo doble-gradiente, al igual que se utiliza en el Gradiente Conjugado, con las salvedades correspondientes que contemplen la no simetr´ıa del sistema. 6.4. M´ ETODO DEL GRADIENTE CONJUGADO PRECONDICIONADO Dado que en la Modelizaci´on de los Campos de Viento el tipo de matrices que aparecen en sus Sistemas de Ecuaciones, son Sim´etricas Definidas Positivas (SDP), el mejor m´etodo iterativo para resolverlas es el del Gradiente Conjugado (GC), cuyo algoritmo ya ha sido expuesto anteriormente (Secci´on 4.4). Teniendo en cuenta, adem´as, que mediante un Precondicionamiento adecuado se consigue mejorar la convergencia del proceso iterativo, es por ello que adoptamos, como m´etodo m´as adecuado para la resoluci´on de los Sistemas de Ecuaciones originados en los Modelos de Campos de Viento, el del Gradiente Conjugado Precondicionado (GCP) [74]. De hecho, todos los ejemplos realizados, cuyos resultados se recogen en el apartado de Experimentos Num´ericos, han sido afrontados con este procedimiento, por considerarlo el m´as id´oneo, dado que el tipo de matrices sobre las que se aplica son SDP. En los casos de Precondicionamiento expuestos anteriormente, tanto por la izquierda, como por la derecha, aunque la matriz original del sistema, A, sea Sim´etrica Definida Positiva, las nuevas matrices precondicionadas, ˜ A, tanto ˜ A=M−1A como ˜ A=A M−1 en general, no tienen por que continuar siendo Sim´etricas, de ah´ı la necesidad de introducir estrategias que al Precondicionar contin´uen conservando la Simetr´ıa, para que el m´etodo del Gradiente Conjugado no pierda, sino que mejore, su
PRECONDICIONADORES EXPL´ ICITOS 79 N´otese la simplificaci´on que supone que, para el c´alculo de los productos internos, s´olo se requiere el uso de las filas de A, lo cual hace el procedimiento bastante atractivo, puesto que no es necesario almacenar expl´ıcitamente la matriz Acompleta. Una vez calculados ZyD, la soluci´on del sistema A x =b, puede computarse como x∗=A−1b=Z D−1ZTb= n X i=1 zT ib pizi Aunque Asea una matriz sparse, el coste de esta algoritmo, aplicado tal como se ha descrito, lo hace inviable, ya que Ztiende a ser una matriz densa. Para evitar este efecto fill-in, conforme se van obteniendo los vectores zi, se efect´ua la simplificaci´on de despreciar aquellas entradas inferiores en valor absoluto a una cierta tolerancia prefijada 0 ≤δ≤1. Llamando ˜ Za la matriz triangular obtenida: A−1≈˜ Z˜ D−1˜ ZT que ser´ıa la incompleta factorizada de A. 6.6.2. PRECONDICIONADOR SAINV Aunque no es frecuente, el algoritmo AINV puede conducir (para matrices A,SDP) a valores de pi(entradas de la diagonal principal) nulos o negativos. En el primer caso, en cuanto se obtuviera un picero, no podr´ıa proseguir el proceso,puesto que los pison denominadores en las f´ormulas de recurrencia que se utilizan en el algoritmo AINV para calcular las zj i, y, si lo que ocurre, es la aparici´on de un pi<0,como A−1≈Z D−1ZT, dar´ıa lugar a una aproximada inversa que ya no es Definida Positiva Evidentemente, esto no ocurre en el caso te´orico de calcular exactamente los valores de zi, ya que siendo A, Sim´etrica Definida Positiva, la expresi´on (6.9) resultar´a pi=zT iA zi>0.
PRECONDICIONADORES EXPL´ ICITOS 80 La raz´on de este breakdown es la siguiente: Los pivotes pi, entradas de la diagonal principal, D, se obtienen, seg´un el algoritmo AINV, por la relaci´on pi=p(i−1) i=aT iz(i−1) i= i−1 X l=1 ailzl−1 li +aii (1 ≤i≤n), si hacemos cero algunas de las entradas de los vectores zi, algunos de los productos ailzli desaparecer´ıan y, en el caso general de A, SDP, estos sumandos que se eliminan pueden resultar positivos con lo que realmente estamos haciendo, al prescindir de ellos, es disminuir el valor te´oricamente positivo de los pi, pudiendo llegar a hacerse nulos o negativos. Llamando ˜zia los vectores modificados (al eliminar las entradas inferiores a δ), para algunas matrices puede ocurrir que el c´alculo de los correspondientes pivotes sea: ˜pi=aT i˜zi˜zT iA˜zi con lo cual se va perdiendo la ortogonalidad de los vectores ˜zi, aumentando as´ı las probabilidades del breakdown. Para robustecer el proceso, [12], el algoritmo SAINV propone que los ˜pise computen usando la expresi´on ˜pi= ˜zT iA˜zi, procedimiento algo m´as costoso que el del AINV, pero m´as seguro y que garantiza la no ruptura del mismo. As´ı el c´alculo de los ˜pise puede expresar como ˜pi= ˜vT i˜zidonde ˜vT i= ˜zT iA. La diferencia de coste depender´a de cuanto m´as denso resulte el vector ˜vien comparaci´on con aT i. Los experimentos realizados demuestran que el resultado del nuevo procedimiento da lugar a un precondicionador de m´as alta calidad y coste muy parecido al anterior, con lo cual compensa definitivamente el uso de una aproximaci´on algo m´as costosa. Por todo ello, el nuevo algoritmo de la factorizaci´on inversa estabilizada (SAINV) puede escribirse de la siguiente forma:
PRECONDICIONADORES IMPL´ ICITOS 81 ALGORITMO SAINV (1) Hacer z(0) i=ei(1 ≤i≤n) (2) Para i= 1,2, ..., n hacer vi=A zi−1 i (3) para j=i, i + 1, ..., n hacer p(i−1) j=vT iz(i−1) j fin. si i=nir a (4) para j=i+ 1, ..., n hacer z(i) j=z(i−1) j− p(i−1) j p(i−1) i!z(i−1) i fin. fin. (4) Hacer zi=z(i−1) iypi=p(i−1) i,para 1 ≤i≤n. Volver a Z= [z1, z2, ...., zn] y D=diag(p1, p2, ..., pn) Obviamente los algoritmos AINV y SAINV [30] son matem´aticamente equivalentes; sin embargo, con el estabilizado se consigue una aproximada inversa m´as fiable. Aplic´andolo a cualquier matriz SDP se obtiene un proceso sin rupturas. 6.7. PRECONDICIONADORES IMPL´ ICITOS 6.7.1. POR COMPARACI´ ON CON EL M´ ETODO DE RICHARDSON El m´etodo de Richardson, m´etodo iterativo muy simple, de relaci´on de recurrencia xi+1 =xi+α(b−A xi)
PRECONDICIONADORES IMPL´ ICITOS 82 con α > 0, permite, por comparaci´on con otros m´etodos iterativos, definir cierta matrices utilizables como precondicionadores. Para ello, contrastaremos la f´ormula de recurrencia para la soluci´on que resulta de aplicar el m´etodo Richardson al sistema precondicionado por una matriz gen´erica M, con la f´ormula correspondiente que se obtiene aplicando los m´etodos, te´oricamente superiores, de Jacobi, SOR y SSOR al sistema sin precondicionar. -Precondicionador de Jacobi ´o Diagonal Aplicando el m´etodo de Richardson al sistema precondicionado M−1A x =M−1b queda, para el c´alculo de los sucesivos valores de la soluci´on: xi+1 =xi+α(M−1b−M−1A xi) y, multiplicando por la matriz de precondicionamiento, resulta: M xi+1 =M xi+α(b−A xi).(6.10) Por otro lado, considerando la descomposici´on de la matriz Ade la forma A=D−E−F, (siendo Dla matriz diagonal formada por los elementos de la diagonal principal de AyE,Fmatrices triangulares), y utilizando el m´etodo de Jacobi para la resoluci´on del sistema A x =b, resulta, xi+1 =D−1(E+F)xi+D−1b multiplicando por la matriz diagonal, y operando: D xi+1 =D xi+ (b−A xi).(6.11) Comparando las expresiones de recurrencia finales de ambos m´etodos, (6.10) y (6.11), se observa que el m´etodo de Jacobi aplicado al sistema sin precondicionar, equivale al de Richardson, con α= 1, menos robusto y m´as simple, cuando se aplica al sistema precondicionado con la matriz diagonal D=diag(A). Resulta as´ı un precondicionador elemental, que se conoce como precondicionador Diagonal, f´acil de implementar y con matriz inversa que se determina con
PRECONDICIONADORES IMPL´ ICITOS 83 muy bajo coste computacional. -Precondicionador SOR Aplicando el m´etodo SOR al sistema A x =b, y con la misma descomposici´on anterior A=D−E−F, queda para la f´ormula de recurrencia de la soluci´on: xi+1 = (D−w E)−1[(1 −w)D+w F]xi+w(D−w E)−1b, donde wes el llamado par´ametro de relajaci´on y, operando convenientemente resulta: (D−w E)xi+1 = (D−w E)xi+w(b−A xi).(6.12) Comparando de nuevo con el m´etodo de Richardson aplicado al sistema precondicionado, (6.10), definir´ıamos, en esta ocasi´on, la matriz de precondicionamiento como: M= (D−w E). -Precondicionador SSOR Aplicando ahora el m´etodo SSOR al sistema sin precondicionar, se obtiene para la soluci´on: xi+1 =D w−F−11−w wD+ED w−E−11−w wD+Fxi+ +D w−F−12−w wDD w−E−1 b operando, para expresar esta relaci´on de forma que se pueda comparar con la (6.10), queda 1 w(2 −w)(D−wE)D−1(D−wF)xi+1 =1 w(2 −w)(D−wE)D−1(D−wF)xi+(b−A xi) (6.13) con lo que resulta como matriz de precondicionamiento 1 w(2 −w)(D−wE)D−1(D−wF). En el caso de aplicar este precondicionador a sistemas sim´etricos, ya que en estos casos, (D−wF) = (D−wE)T, se puede expresar como un producto de dos
PRECONDICIONADORES IMPL´ ICITOS 84 matrices triangulares transpuestas, M="(D−wE)D−1/2 pw(2 −w)#"(D−wE)D−1/2 pw(2 −w)#T Para sistemas no sim´etricos, se puede escribir como producto de dos matrices triangulares, inferior y superior, respectivamente M= (I−wED−1)D−wF w(2 −w) 6.7.2. POR FACTORIZACIONES INCOMPLETAS La propiedad que en principio, debe verificar una matriz de precondicionamiento M, de ser una aproximaci´on, m´as o menos cercana, de la matriz de coeficientes del sistema, sugiere que un procedimiento para obtenerla sea el descomponer Ade la forma A≈A1A2, de tal manera, que no suponga un excesivo esfuerzo computacional, y adoptar para la misma M=A1A2 Adem´as de otras factorizaciones posibles destacaremos dos de ellas: la basada en la factorizaci´on en dos matrices triangulares LU, usualmente utilizada en la resoluci´on de sistemas por m´etodos directos, y la factorizaci´on incompleta de Cholesky . -Precondicionador ILU(0) Resulta de descomponer Aen dos matrices triangulares, inferior y superior, respectivamente, LyU, A≈LU =M cuyos elementos, mij, sean tales que: mij = 0 si aij = 0 (A−LU)ij = 0 si aij 6= 0
PRECONDICIONADORES IMPL´ ICITOS 85 es decir, que los elementos nulos de la matriz del sistema, contin´uen siendo nulos en las posiciones respectivas de las matrices triangulares. Si no realiz´aramos esta simplificaci´on, el coste computacional se incrementar´ıa y equivaldr´ıa a resolver el sistema por un m´etodo directo. -Precondicionador ILU(n) Para las matrices tridiagonales o pentadiagonales que resultan de la discretizaci´on de problemas con ecuaciones en derivadas parciales el´ıpticas, se pueden practicar otros niveles de factorizaci´on, consistentes en rellenar alguna diagonal de las matrices factores de la descomposici´on, que en ILU(0) ser´ıan nulas. La aproximaci´on ser´ıa mayor a costa de incrementar el trabajo computacional de la descomposici´on. -Factorizaci´on incompleta de Cholesky Est´a especialmente indicado para factorizar matrices Sim´etricas Definidas Positivas, (SDP).Su objetivo consiste en descomponer Aen tres matrices, una matriz central diagonal y dos matrices laterales: una triangular inferior y su transpuesta. Sea A= (aij) una matriz n×n, SDP. En orden a realizar su factorizaci´on podemos considerarla formada por A= (aij) = a11 fT 1 f1A2 o sea, compuesta por: su primer elemento, a11; una matriz columna, (n−1)×1, f1, formada por los elemtos de su primera columna menos el primero; su transpuesta, la matriz fila, 1 ×(n−1), fT 1y la matriz A2, (n−1) ×(n−1), SDP, resultado de eliminar la primera fila y la primera columna de la matriz inicial A. A partir de esta estructura, como primer paso de su factorizaci´on, f´acilmente puede descomponerse en un producto de tres matrices, de la siguiente forma: A= (aij) = a11 fT 1 f1A2 = a11 0 f1I a−1 11 0 0C2 a11 fT 1 0I =L1Z1LT 1 Siendo Ila matriz unidad y C2una matriz (n−1) ×(n−1), SDP, obtenida a partir de A2, tal que: C2=A2−1 a11 f1fT 1= (aij)(2),∀i, j ≥2
PRECONDICIONADORES IMPL´ ICITOS 86 en la que, por tanto, su primer elemento ser´a: (a22)(2) =a22 −a2 12 a11 . Esta nueva matriz C2, as´ı obtenida, a su vez puede estructurarse de la misma forma que se hizo con la matriz inicial A: C2= (aij)(2) = a(2) 22 fT 2 f2A3 Siendo f2una matriz columna (n−2)×1 y A3una matriz SDP, (n−2)×(n−2).En orden a evitar el efecto fill-in, que aparecer´a al calcular los elementos (ai2)(2), para i≥3, correspondientes a la matriz columna f2y su transpuesta, realizaremos una aproximaci´on, tomando en su lugar una nueva matriz columna, l2, que tenga las mismas entradas que f2, pero manteniendo nulas aquellas que lo son en la matriz inicial A, y as´ı resultar´a: C2= (aij)(2) = a(2) 22 fT 2 f2A3 ≈ a(2) 22 lT 2 l2A3 nueva matriz, igual de sparse que la inicial, que volvemos a factorizar de la misma forma anterior, en una matriz central y dos triangulares a ambos lados: C2= (aij)(2) ≈ a(2) 22 lT 2 l2A3 = a(2) 22 0 l2I a(2)−1 22 0 0C3 a(2) 22 lT 2 0I siendo ahora C3= (aij)(3) =A3−1 a(2) 22 l2lT 2 una matriz SDP, (n−3) ×(n−3),∀i, j ≥3. En dicha matriz su primera entrada ahora ser´a: a(3) 33 =a(2) 33 −a(2)2 23 a(2) 22 . Con las factorizaciones realizadas hasta ahora, teniendo en cuenta que f1=l1, la matriz inicial Ase podr´a descomponer de la siguiente forma: A≈ a11 0 l1I 1 0 0 0a(2) 22 0 0l2I a−1 11 0 0 0a(2)−1 22 0 0 0 C3 1 0 0 0a(2) 22 lT 2 0 0 I a11 lT 1 0I =
PRECONDICIONADORES IMPL´ ICITOS 87 =L1L2Z2LT 2LT 1 producto de una matriz central y dos matrices triangulares a cada lado. Continuando con la factorizaci´on de C3, de forma similar a la anterior, puede descomponerse la matriz inicial como: A≈ a11 0 l1I 1 0 0 0a(2) 22 0 0l2I 1 0 0 0 0 1 0 0 0 0 a(3) 33 0 0 0 l3I a−1 11 0 0 0 0a(2) 22 −10 0 0 0 a(3) 33 −10 0 0 0 C4 1 0 0 0 0 1 0 0 0 0 a(3) 33 lT 3 0 0 0 I 1 0 0 0a(2) 22 lT 2 0 0 I a11 lT 1 0I =L1L2L3Z3(L1L2L3)T Y procediendo de la misma forma, al llegar a la n-sima factorizaci´on resultar´a: A≈L1L2L3.....LnZn(L1L2L3.....Ln)T=L D LT tal que L= a11 000··· 0 a12 a(2) 22 0 0 ··· 0 a13 a(2) 23 a(3) 33 0··· 0 a14 a(2) 24 a(3) 33 a(4) 44 ··· 0 · · · · · · · · · · · · · · · · · · · · · · · · a1na(2) 2na(3) 3na(4) 44 ···a(n) nn y
PRECONDICIONADORES IMPL´ ICITOS 88 D= a−1 11 000··· 0 0a(2) 22 −10 0 ··· 0 0 0 a(3) 33 −10··· 0 0 0 0 a(4) 44 −1··· 0 · · · · · · · · · · · · · · · · · · · · · · · · 0 0 0 0 ···a(n) nn −1 Con lo cual quedar´a la matriz Aaproximadamente igual a un producto de dos matrices triangulares laterales, LyLT, y una matriz central diagonal. Los precondicionadores impl´ıcitos m´as usuales son el DIAGONAL y el ILU(0). El DIAGONAL, es con mucho, el que menos esfuerzo computacional exige. Se calcula directamente tomando los elementos de la diagonal principal de Ay los productos matriz inversa por vector se efect´uan, asimismo, de forma inmediata. Sin embargo su aplicaci´on se ve reducida, puesto que en muchos sistemas mal condicionados, representativos de problemas f´ısicos con capas l´ımites, singularidades o condiciones de contorno especiales, no mejoran sustancialmente la convergencia. Los precondicionadores ILU(0), exigen m´as esfuerzo computacional inicial para su construcci´on y los productos de matriz inversa por vector se realizan por procesos de remonte, pero su aplicaci´on da lugar a buenos resultados en sistemas que con el precondicionador DIAGONAL no convergen.
ALGORITMO MULTICOLORING (MC) 94 En general el n´umero de colores necesarios no exceder´a al del m´aximo grado de cada nodo +1. Una vez asignados los colores a todos los nodos del grafo asociado a la matriz, se reordena ´esta reagrupando todos los v´ertices, es decir, todos los elementos de la diagonal principal, del mismo color; as´ı se consigue una nueva estructura de la matriz reordenada por bloques, en la que los bloques diagonales ser´an precisamente matrices diagonales y el n´umero de bloques coincidir´a con el n´umero de colores. El resto de las entradas configuraran dos matrices triangulares sparse situadas a ambos lados de los bloques diagonales. Si bien la t´ecnica del multicoloring resulta muy barata, al utilizar los precondicionadores ILU(0) puede ocurrir, seg´un Saad, que el n´umero de iteraciones, necesarias para alcanzar la convergencia, probablemente resulte mucho m´as alto precondicionando la matriz con el reordenamiento multicolor, que precondicionando la matriz original directamente.
Cap´ıtulo 8 PRECONDICIONAMIENTO DE SISTEMAS DE ECUACIONES LINEALES DE MATRIZ VARIABLE 8.1. PROPUESTA DE ESTRATEGIA El objetivo principal de esta tesis consiste en extender las t´ecnicas de Precondicionamiento, as´ı como las de Reordenaci´on, a los sistemas de ecuaciones lineales de matrices variables [105, 96], que surgen de la modelizaci´on de campos de viento [75], que como ya se ha visto anteriormente, son del tipo: Aεxε=bε donde εrepresenta el par´ametro de estabilidad del modelo y Aε=M+ε N siendo MyNmatrices constantes para un nivel de discretizaci´on dado y, adem´as, Sim´etricas Definidas Positivas (SDP), por lo que, para su resoluci´on, utilizaremos siempre el algoritmo del Gradiente Conjugado Precondicionado. Para precondicionar estos sistemas se podr´ıa recurrir, en principio, a dos estrategias extremas:
PROPUESTA DE ESTRATEGIA 96 a) Por un lado, se podr´ıa construir un ´unico precondicionador para un cierto valor de ε,εo, y lo aplicar´ıamos para la resoluci´on de los distintos sistemas que surgen para cada valor de , lo que conducir´a a convergencias cada vez m´as lentas a medida que los valores de εse vayan alejando del valor inicial. b) Por otro lado, y´endonos al extremo opuesto, usar´ıamos un precondicionador diferente para cada sistema, o sea para cada valor de ε, lo que resultar´ıa muy costoso. En esta tesis se propone una soluci´on intermedia, que consiste en construir un precondicionador, que pueda ser actualizado f´acilmente para cada valor de ε, y cuya aplicaci´on al algoritmo del Gradiente Conjugado de lugar a velocidades de convergencia comprendidas entre las conseguidas al aplicar las estrategias extremas mencionadas. Con lo que en definitiva se conseguir´a mejorar el grado de eficacia del algoritmo a utilizar en la resoluci´on del sistema. Este tipo de soluci´on intermedia ya ha sido propuesto por Benzi [11] y Meurant [65], independientemente, usando cada uno de ellos un modelo de Precondicionador diferente que, a un bajo coste computacional, se adaptan f´acilmente para cada valor del par´ametro, ε; estrategia que conduce a conseguir un grado de convergencia, con el algoritmo del Gradiente Conjugado Precondicionado, intermedio entre las convergencias alcanzables usando las dos opciones extremas. Benzi desarrolla su estudio utilizando un Precondicionador Expl´ıcito, mediante la construcci´on de Inversas Aproximadas, usando el algoritmo SAINV, para el caso especial de matrices variables, Aε, del tipo Aε=M+ε I siendo Ila matriz unitaria. Sin embargo Meurant lo hace para la matriz Aε=M+ε D siendo D una matriz diagonal y considerando un Precondicionador Impl´ıcito, a partir de una Factorizaci´on Incompleta de Cholesky de la matriz M.
ADAPTACI ´ ON DEL PRECONDICIONADOR SAINV 97 8.2. ADAPTACI´ ON DEL PRECONDICIONADOR SAINV En este estudio, proponemos seguir un camino paralelo al propuesto por Benzi [97] para el caso especial de matrices variables del tipo: Aε=M+ε I siendo Muna matriz SDP e Ila matriz unitaria; pero extendi´endolo al caso m´as gen´erico de Aε=M+ε N Con el algoritmo SAINV, ya expuesto en el Cap´ıtulo 6, (6.6.2), se puede construir una inversa aproximada factorizada de una matriz A, SDP, a partir de la obtenci´on por congruencia de una forma diagonal de la misma: ZTA Z =D=diag(d1, d2, ....., dn) mediante la matriz triangular superior Z= [z1, z2, ....., zn] conseguida por un proceso de A-conjugaci´on de Grand-Schmidt, a partir del conjunto de vectores unitarios linealmente independientes {e1, e2, ....., en} ∈ <n, donde dj=zT jA zj>0,1≤j≤n. Si realizamos el proceso de c´alculo de los vectores zide forma incompleta, descartando, en cada caso, las entradas respectivas menores de una cierta tolerancia escogida, δ, tal que: 0 < δ < 1, se obtiene una matriz sparse aproximada de Z, que denominamos ˜ Z, con la que se podr´a construir la matriz aproximada inversa de A: A−1≈˜ Z˜ D−1˜ ZT. Pues bien, aplicando dicho algoritmo, se puede obtener una aproximada inversa de M: M−1≈˜ Z˜ D−1˜ ZT=P−1 y, a partir de ah´ı, se considera un precondicionador para Aε=M+ε N de la forma: P−1 ε=˜ Z(˜ D+ε E)−1˜ ZT
ADAPTACI ´ ON DEL PRECONDICIONADOR SAINV 98 donde Eser´ıa una matriz gen´erica, sim´etrica, a determinar, f´acilmente computable, que haga a ( ˜ D+ε E) Sim´etrica Definida Positiva y tal que los productos P−1 εpor vector, que figuran en el Gradiente Conjugado, no supongan un coste elevado. A efectos de definir Ey, dando por supuesto que la inversa exacta de Mser´ıa M−1=Z D−1ZT, se establece la diferencia Pε−Aε=Z−T(D+ε E)Z−1−(M+ε N) = ε(Z−TEZ−1−N). Si se tomara E=ZTNZ, resultar´ıa: Pε−Aε= 0, con lo cual se conseguir´ıa el precondicionador ideal P−1 ε=A−1 ε. Evidentemente, esto no es viable, dado que no se dispone de la matriz Z, sino de su aproximaci´on ˜ Z, pero ello sugiere como mejor expresi´on para la matriz E: E=˜ ZTN˜ Z y as´ı, de esta forma, Ecumplir´ıa con las condiciones necesarias mencionadas anteriormente. En lugar de iniciar el proceso con la obtenci´on de la aproximada inversa de M, que se corresponde con la aproximada inversa de Aε=M+ε N, para ε= 0, se puede obtener inicialmente una aproximada inversa de Aε0=M+ε0N. Y as´ı las sucesivas matrices, para los distintos valores de εse escribir´ıan: Aε=M+ε N =Aε0−ε0N+ε N =Aε0+ ∆ε N donde ∆ε=ε−ε0. Dado que εes siempre positivo, la matriz Aε, obviamente, es definida positiva, aunque ∆εsea negativo. De todo ello resulta que la aproximada inversa ser´ıa: A−1 ε0≈˜ Z˜ D−1˜ ZT=P−1 ε0 y el precondicionador para la matriz variable: P−1 ε=˜ Z(˜ D+ ∆ε E)−1˜ ZT
ADAPTACI ´ ON DE LA FACTORIZACI ´ ON DE CHOLESKY 99 siendo por supuesto E=˜ ZTN˜ Z, al igual que antes y con las mismas prestaciones. Una opci´on para construir E, con los requisitos previstos, es tomar una aproximaci´on de ˜ Zque nombraremos como ˜ Zk, que se obtiene extrayendo solamente su diagonal principal, si k= 1, y adem´as sus k−1 diagonales superiores, si k > 1, y considerando para Nla aproximaci´on Nh, extrayendo su diagonal principal, si h= 1, y las h−1 diagonales secundarias para h > 1. Con lo cual denominaremos Eh,k =˜ ZT kNh˜ Zk En la pr´actica, a efectos de no incrementar el coste por iteraci´on del Gradiente Conjugado, resulta ´util considerar las parejas h= 1 y k= 2 ´o h= 2 y k= 1, que dan lugar, respectivamente, a las matrices E1,2´o E2,1, tridiagonales. Incluso se puede conseguir una mayor simplificaci´on considerando h=k= 1, resultando as´ı la matriz E1,1, diagonal. 8.3. ADAPTACI´ ON DE LA FACTORIZACI´ ON DE CHOLESKY Aqu´ı optamos por generalizar la factorizaci´on incompleta de Cholesky, propuesta por Meurant [98] para el caso de matrices Aε=M+εD, siendo Duna matriz diagonal, al caso m´as general de las matrices Aε=M+ε N, siendo MyN dos matrices sim´etricas definidas positivas n×n. As´ı podremos escribir Aεcomo sigue: Aε= (mij) + ε(nij) = m11 +εn11 (f1M+εf1N)T f1M+εf1NM2+εN2 donde f1M, f1Nrepresentan matrices columnas ((n−1,1) y M2, N2matrices de orden n−1. Factorizando Aε, Aε= m11 +εn11 0 l1M+εl1NI (m11 +εn11)−10 0C2 m11 +εn11 (l1M+εl1N)T 0 I con lo que nos queda, Aε=L1Z1LT 1(8.1)
ADAPTACI ´ ON DE LA FACTORIZACI ´ ON DE CHOLESKY 100 siendo l1M=f1Myl1N=f1N. Identificando, t´ermino a t´ermino, se obtiene para la matriz C2: C2=M2+εN2−1 m11 +εn11 (l1M+εl1N) (l1M+εl1N)T(8.2) Si, a efectos de construir el precondicionador, tomamos como primera aproximaci´on s´olo los elementos de la diagonal de N, ser´ıa l1N= 0, quedando C2=εD2+M2−1 m11 +εn11 l1MlT 1M y la aproximaci´on de orden cero, C2=εD2+M2−1 m11 l1MlT 1M con lo cual, las entradas de C2se obtendr´ıan f´acilmente, a˜nadiendo εD2a la matriz que resulta en la factorizaci´on de M. Otra aproximaci´on consiste en considerar, en (8.2), todas las entradas de N2 y despreciar los productos εl1N, quedando C2, de forma similar, C2=εN2+M2−1 m11 l1MlT 1M Continuando con esta aproximaci´on, C2=εN2+ m(2) 22 fT 2M f2MM3 = m(2) 22 +εn22 (f2M+εl2N)T f2M+εl2NM3+εN3 Haciendo ceros en f2Mlas correspondientes entradas nulas de M, para evitar el efecto fill-in, obtenemos l2M. Factorizando C2: C2≈ m(2) 22 +εn22 0 l2M+εl2NI m(2) 22 +εn22−10 0C3 m(2) 22 +εn22 (l2M+εl2N)T 0 I donde, C3=M3+εN3−1 m(2) 22 +εn22 (l2M+εl2N) (l2M+εl2N)T. Utilizando las mismas simplificaciones anteriores, se consigue C3=M3+εN3−1 m(2) 22 l2MlT 2M= m33 +εn33 (f3M+εl3N)T f3M+εl3NM4+εN4
ADAPTACI ´ ON DE LA FACTORIZACI ´ ON DE CHOLESKY 101 resultando as´ı, para C3, la misma ley de formaci´on que se obtuvo para C2. Con estos criterios de formaci´on de las matrices Ci, la factorizaci´on aproximada de Aε,queda: Aε≈L1Z1LT 1=L1L2Z2LT 2LT 1= (L1L2···Ln)Zn(L1L2···Ln)T(8.3) siendo Znla matriz diagonal (m11 +εn11)−1. m(2) 22 +εn22−1. m(3) 33 +εn33−1. . . . . . .m(n) nn +εnnn−1 Las entradas diagonales de la matriz triangular inferior (L1L2···Ln) ser´an m(i) ii + εnii. Y las respectivas columnas inferiores a los elementos diagonales vendr´an definidas por matrices ljM +εljN de orden (n−j)×1.
APLICACIONES TEST 108 εICHOL(Aε0) ICHOLDICHOLNFull-ICHOL 0noIter. t(s) - - - - - - 184 6.21 10−6noIter. t(s) 184 5.96 184 5.99 184 6.00 184 6.21 10−5noIter. t(s) 184 5.96 184 5.99 184 6.00 184 6.23 10−4noIter. t(s) 184 5.98 184 5.99 184 6.01 184 6.21 10−3noIter. t(s) 181 5.89 181 5.90 181 5.91 181 6.12 10−2noIter. t(s) 170 5.55 170 5.56 170 5.56 169 5.73 10−1noIter. t(s) 148 4.81 135 4.43 131 4.29 126 4.35 1noIter. t(s) 232 7.50 149 4.88 105 3.46 78 2.79 10 noIter. t(s) 454 14.66 303 9.84 145 4.74 76 2.73 102noIter. t(s) 995 32.09 675 22.03 261 8.46 >5000 – 103noIter. t(s) 1452 46.58 965 31.29 354 11.50 >5000 – 104noIter. t(s) 1583 50.73 1049 33.98 384 12.45 >5000 – 105noIter. t(s) 1604 51.40 1059 34.25 388 12.57 >5000 – 106noIter. t(s) 1605 51.43 1060 34.29 388 12.58 >5000 – Tabla 9.4: Ejemplo 2, 43.954 ecuaciones:N´umero de iteraciones y tiempo de computaci´on (en s.) del Gradiente Conjugado con diferentes Precondicionadores por Factorizaci´on Incompleta de Cholesky
APLICACIONES TEST 109 Observando los resultados de la Tabla 9.4, con los Precondicionadores por Factorizaci´on Incompleta, podemos concluir que para peque˜nos valores de εno es necesario adaptar la factorizaci´on inicial puesto que con ella se consigue alcanzar la convergencia a muy bajo coste. Sin embargo, con valores altos de ε, el ICHOLN tiene el mejor comportamiento. Conviene tener en cuenta que para los valores de εcomprendidos entre 1 y 10 la re-computarizaci´on de la factorizaci´on incompleta resulta m´as aconsejable.
APLICACIONES TEST 110 9.1.4. EJEMPLO 3 En las Tablas 9.5 y 9.7 se presentan los resultados conseguidos para el sistema de 98.999 ecuaciones, obtenido con otro refinamiento posterior, utilizando los mismos tipos de Precondicionadores. En esta caso, tanto con los Precondicionadores SAINV como ICHOL, en todas sus variantes, la convergencia fue extremadamente lenta (m´as de 5000 iteraciones) para valores de ε≥104, por lo que a partir de ε= 103ya no se recogen los resultados en ambas tablas. εFull-SAINV SAINV11 SAINV12 SAINV21 SAINV(Aε0) 0Iter. t(seg.) 278 3541.68 - - - - - - - - 10−6Iter. t(seg.) 278 3540.47 278 25.83 278 26.70 278 30.21 278 19.04 10−5Iter. t(seg.) 278 3541.55 278 25.78 278 26.71 278 30.24 279 19.05 10−4Iter. t(seg.) 280 3566.12 279 25.81 278 26.73 279 30.30 279 19.08 10−3Iter. t(seg.) 277 3540.38 278 25.82 276 26.61 278 30.17 276 18.37 10−2Iter. t(seg.) 258 3538.34 257 24.32 258 25.31 257 28.48 256 17.39 10−1Iter. t(seg.) 214 3543.85 227 22.34 229 23.21 224 25.79 307 20.82 100Iter. t(seg.) 193 3524.28 283 26.14 291 27.64 275 30.00 653 44.06 101Iter. t(seg.) 257 3591.23 591 47.11 589 48.93 581 54.65 1724 116.03 102Iter. t(seg.) 548 3724.73 1724 124.27 1670 126.10 1659 127.71 >5000 – 103Iter. t(seg.) 1290 3785.34 4234 297.10 4002 304.06 4128 305.60 >5000 – Tabla 9.5: Ejemplo 3, 98.999 ecuaciones: N´umero de iteraciones y tiempo de computaci´on (en segundos) del Gradiente Conjugado para los distintos Precondicionadores SAINV Una vez m´as se comprueba que para valores peque˜nos de εbasta con un precondicionador ´unico, el SAINV(Aε0); es especialmente notorio que este tipo de precondicionador empieza a fallar a partir de ε= 102. Sin embargo para valores de ε≥1 el SAINV11 presenta los mejores resultados.
APLICACIONES TEST 111 εOrden Inicial MN (69,43s) RCM (0,70s) MC (0,45s) 0Iter. t(s) 278 3541.68 264 3092.93 273 2757.27 263 4626.74 10−6Iter. t(s) 278 25.83 265 16.35 274 16.25 263 20.13 10−5Iter. t(s) 278 25.78 264 16.29 274 16.22 263 20.15 10−4Iter. t(s) 279 25.81 264 16.27 273 16.18 263 20.15 10−3Iter. t(s) 278 25,82 262 16.16 272 16.11 260 19.93 10−2Iter. t(s) 257 24.32 236 14.61 247 14.66 240 28.40 10−1Iter. t(s) 227 22.34 208 12.89 212 12.58 219 16.81 1Iter. t(s) 283 26.14 261 16.15 265 15.71 270 20.69 10 Iter. t(s) 591 47.11 520 32.04 549 32.32 553 42.24 102Iter. t(s) 1724 124.27 1512 92.96 1589 92.23 1623 123.67 103Iter. t(s) 4234 297,10 3701 227.79 3756 220.51 3924 294.04 Tabla 9.6: Ejemplo 3, 98.999 ecuaciones: N´umero de iteraciones y tiempo de computaci´on (en segundos) del Gradiente Conjugado con el Precondicionador SAINV11 para diferentes Reordenaciones Dado que con el precondicionador SAINV11 se consigui´o un mejor comportamiento, fue elegido para aplicarlo tambi´en sobre los Sistemas Reordenados [28, 106] utilizando los Algoritmos citados en el Cap´ıtulo 7: M´ınimo Vecino (MN), CuthillMcKee Inverso (RCM) y Multicoloring (MC). Los resultados se recogen en la Tabla 9.6
APLICACIONES TEST 112 Como puede observarse, reordenando con los Algoritmos MN y RCM, mejora la convergencia del Gradiente Conjugado Precondicionado. Sin embargo el Algoritmo MC no es tan eficiente, como ya era de prever, de acuerdo con los comentarios expuestos en el apartado 6.5 sobre su falta de eficacia. (a) Orden Inicial (b) M´ınimo Vecino (c) Cuthill-McKee Inverso (d) Multicoloring Figura 9.1: Patrones de ‘sparsidad’ de las matrices, con su orden inicial y una vez reordenadas
APLICACIONES TEST 113 En la figura 9.1 se muestran los patrones de sparsidad de la matriz inicial del sistema y de sus reordenaciones, seg´un los Algoritmos indicados. La imagen correspondiente a la matriz con el reordenamiento Cuthill-McKee Inverso es la que m´as se asemeja a una matriz diagonal y en efecto es la reordenaci´on que produce mejores resultados. Sin embargo la simple observaci´on del patr´on obtenido con el reordanamiento Multicoloring ya nos anticipa que los resultados a conseguir no van a ser mucho mejores que los del Gradiente Conjugado Precondicionado aplicado directamente sobre la matriz con su orden inicial, sobre todo para valores elevados de ε. En la Tabla 9.7 se puede comprobar como en este ejemplo, para 98.999 ecuaciones, se obtienen resultados similares a los conseguidos con los precondicionadores ICHOL en los sistemas de 43.954 ecuaciones: en general, para peque˜nos valores de εno es necesario adaptar la factorizaci´on inicial puesto que con ella se consigue alcanzar la convergencia a muy bajo coste. Sin embargo, con valores altos de ε, el ICHOLNtiene el mejor comportamiento. Teniendo en cuenta que para los valores de εcomprendidos entre 1 y 10 la re-computarizaci´on de la factorizaci´on incompleta resulta m´as efectiva. De forma similar a como se hizo con los precondicionadores SAINV, aqu´ı se ha elegido el precondicionador ICHOLN, por su mejor comportamiento, para aplicarlo sobre los sistemas reordenados y comprobar as´ı la eficacia de la Reordenaci´on en los Sistemas de Ecuaciones Variables. Los resultados se recogen en la Tabla 9.8. En este caso, para valores de εcomprendidos entre 0 y 10, se ha conseguido mejorar el coste de las iteraciones con la reordenaci´on del M´ınimo Vecino y a partir de valores de ε≥102los mejores resultados se han conseguido con el reordenamiento de Cuthill-McKee Inverso. Pero teniendo en cuenta que el tiempo de implantaci´on del MN(69,43s.) es muy superior al del RCM(0,7s.), en definitiva, la resoluci´on se abarata mucho m´as con el Cuthill-McKee Inverso. Por lo que respecta al Multicoloring, para ning´un valos de ε, se ha conseguido mejorar la eficacia conseguida con el ICHOLNaplicado directamente a la matriz
APLICACIONES TEST 114 εICHOL(Aε0) ICHOLDICHOLNFull-ICHOL 0noIter. t(s) - - - - - - 201 16.81 10−6noIter. t(s) 201 16.14 201 16.16 201 16.19 201 16.82 10−5noIter. t(s) 201 16.15 201 16.16 201 16.19 201 16.83 10−4noIter. t(s) 201 16.15 201 16.16 201 16.19 201 16.83 10−3noIter. t(s) 201 16.14 201 16.16 200 16.11 200 16.76 10−2noIter. t(s) 188 15.22 191 15.35 189 15.24 189 15.87 10−1noIter. t(s) 225 18.08 157 12.65 155 12.52 151 12.85 1noIter. t(s) 483 38.63 211 16.94 148 11.97 132 11.33 10 noIter. t(s) 1350 107.71 540 43.14 259 20.91 236 19.64 102noIter. t(s) 3973 317.16 1466 116.86 593 47.62 >5000 – 103noIter. t(s) >5000 – 3468 277.11 1269 101.60 >5000 – Tabla 9.7: Ejemplo 3: 98.999 ecuaciones. N´umero de iteraciones y coste computacional (en s.) del Gradiente Conjugado con diferentes Precondicionadores por Factorizaci´on Incompleta de Cholesky inicial, lo cual revela la ineficacia de la reodenaci´on MC con los precondicionadores ICHOL, confirm´andose as´ı la opini´on de Yousef Saad en su obra Iterative Methods [94], donde comenta que sobre todo con los precondicionadores ILU(0), puede ocurrir una p´erdida de eficacia, puesto que el n´umero de iteraciones para alcanzar la convergencia puede aumentar, resultando m´as alto que si se precondicionara directamente la matriz original.
APLICACIONES TEST 115 εInitial Ordering MN (69.43 s) RCM (0.70 s) MC (0.45 s) 0noIter. t(s) 201 16.81 158 12.86 175 13.84 223 24.01 10−6noIter. t(s) 201 16,19 159 12.53 175 13.46 223 22.51 10−5noIter. t(s) 201 16.19 159 12.54 175 13.45 222 22.41 10−4noIter. t(s) 201 16.19 158 12.47 176 13.53 223 22.51 10−3noIter. t(s) 200 16.11 157 12.37 174 13.38 220 22.21 10−2noIter. t(s) 189 15.24 143 11.31 160 12.30 207 20.90 10−1noIter. t(s) 155 12.52 116 9.19 131 10.12 170 17.21 1noIter. t(s) 148 11.97 123 9.73 129 9.97 160 16.20 10 noIter. t(s) 259 20,91 227 17.88 240 18.38 277 27.91 102noIter. t(s) 593 47.62 533 41.61 536 40.73 632 63.43 103noIter. t(s) 1269 101.60 1212 94.45 1177 89.23 1395 139.79 Tabla 9.8: Ejemplo 3, 98.999 ecuaciones: N´umero de iteraciones y tiempo de computaci´on (en segundos) del Gradiente Conjugado con el Precondicionador ICHOLNpara diferentes Reordenaciones
ELECCI ´ ON DEL PAR ´ AMETRO ´ OPTIMO 116 9.2. ELECCI´ ON DEL PAR´ AMETRO ´ OPTIMO 9.2.1. PRELIMINARES Como ya se indicaba en el Cap´ıtulo 4, ESTIMACI´ ON DE PAR´ AMETROS, el ajuste eficiente de un modelo de campo de viento depende, en gran medida, de la adecuada valoraci´on de los par´ametros que aparecen en las distintas etapas del proceso, especialmente de aquel que interviene de forma b´asica en la formulaci´on del sistema (M+ε N)xε=bε puesto que afecta de forma directa al c´alculo del viento resultante. Dicho par´ametro, ε, de estabilidad del modelo, est´a estrechamente relacionado con los m´odulos de precisi´on de Gauss, a trav´es de las relaciones (2.12) y (2.17): ε=α2=Tv Th =α2 1 α2 2 E. Rodr´ıguez en su Tesis Modelizaci´on y simulaci´on num´erica de campos de viento mediante elementos finitos adaptativos en 3D [89], propone la utilizaci´on de Algoritmos Gen´eticos, como m´etodos de optimizaci´on basados en un mecanismo de evoluci´on natural, para realizar la selecci´on autom´atica del par´ametro ε(α), pues se trata de una herramienta robusta, flexible, competitiva y cuyos c´alculos pueden paralelizarse. Ver Cap´ıtulo 4(4.2). Uno de los aspectos m´as importantes de los Algoritmos Gen´eticos es la construcci´on de una poblaci´on inicial y la posterior evaluaci´on, de cada individuo de la misma, de acuerdo con los resultados de su aplicaci´on a una funci´on objetivo. En el caso de la modelizaci´on de campos de viento se adopta como funci´on objetivo la minimizaci´on de las diferencias entre el viento resultante, calculado mediante la resoluci´on del sistema lineal Aεxε=bε, y el observado, en las estaciones de referencia. Ello lleva como consecuencia que para cada individuo de la poblaci´on inicial hay que resolver dicho sistema. Luego, a tenor de los resultados conseguidos, se proceder´a a crear una nueva poblaci´on, mediante operadores de Selecci´on, Cruce y Mutaci´on, que nuevamente hay que someterla a la evaluaci´on de la funci´on objetivo, lo cual supone volver a resolver el sistema tantas veces como valores seleccionados. A partir de estos
ELECCI ´ ON DEL PAR ´ AMETRO ´ OPTIMO 117 valores se vuelve a generar otra nueva poblaci´on, con los mismos criterios, que se vuelve a evaluar, y as´ı sucesivamente, hasta conseguir un individuo que cumpla con el criterio de parada, que ser´a el valor ´optimo del par´ametro a considerar para ese Modelo. Ello implica que en el proceso de selecci´on hay que resolver el sistema Aεxε= bε, tantas veces como posibles valores de ε(α) a estimar, multiplicado por el n´umero de iteraciones a realizar, en cada paso sucesivo del Algoritmo Gen´etico; por lo que se hace necesario disponer de m´etodos eficaces que permitan la implementaci´on r´apida de un precondicionador para cada valor del par´ametro. De ah´ı la importancia de conocer el coste computacional que supone el utilizar distintos tipos de precondicionadores. Como conclusi´on de los tests realizados anteriormente se estableci´o, que los precondicionadores basados en la Factorizaci´on Incompleta de Cholesky, parecen ser la herramienta m´as eficaz para mejorar la convergencia del m´etodo del Gradiente Conjugado, en la resoluci´on de los Sistemas Lineales de Ecuaciones Variables, por tanto estos han sido los tipos de precondicionadores elegidos para comprobar su comportamiento a la hora de evaluar una poblaci´on inicial de par´ametros, ε(α), de cara a la elecci´on de su valor ´optimo mediante Algoritmos Gen´eticos [29]. En este apartado presentamos los resultados obtenidos al resolver los Sistemas Lineales de Ecuaciones Variables, correspondientes a la modelizaci´on de dos campos de viento diferentes, usando siempre el m´etodo del Gradiente Conjugado Precondicionado, el m´as eficaz cuando hay matrices Sim´etricas Definidas Positivas, como es el caso, y utilizando los precondicionadores basados en la Factorizaci´on Incompleta de Cholesky, ya estudiados anteriormente, en sus variantes ICHOL(Aε0), ICHOLD,ICHOLNyFull −ICHOL. Todos los experimentos han sido realizados en un equipo XENON Precisi´on 530 con Fortran de Doble Precisi´on. Los procesos de iteraci´on siempre se han iniciado a partir de un vector nulo y finalizados si ri r0 2≤10−10 o si el n´umero de iteraciones se hac´ıa superior a 10.000.
CONCLUSIONES 123 Positivas, SDP, se justifica el considerar como mejor m´etodo iterativo para su resoluci´on, por su probada eficacia, el Gradiente Conjugado Precondicionado, GCP, lo que conduce a tener que afrontar, como novedoso, el Precondicionamiento de Matrices Variables. Dado que para cada valor del par´ametro, ε, se origina una matriz diferente, el Precondicionamiento de las mismas debe conseguirse de forma r´apida y eficaz, para una amplia gama de valores de su par´ametro. Se propone como estrategia, para conseguirlo, la construcci´on de un ´unico Precondicionador, f´acilmente adaptable a cada valor diferente del par´ametro, o sea un Precondicionador Variable . Esto se plantea a trav´es de dos modelos de Precondionamiento diferentes: uno Expl´ıcito, construyendo Inversas Aproximadas (SAINV), y otro Impl´ıcito, fundamentado en la Factorizaci´on Incompleta de Cholesky (ICHOL). La amplia gama de experimentos num´ericos realizados, con ambos tipos de Precondicionadores, nos permite obtener las siguientes conclusiones: Por lo que al Precondicionamiento SAINV se refiere, puede afirmarse que para peque˜nos valores del par´ametro, ε, no parece rentable la construcci´on de un Precondicionador Variable, bastar´ıa con disponer de un ´unico Precondicionador, elaborado a partir de un valor inicial de ε, al que se ha denominado SAINV(Aε0). Sin embargo, para valores de ε≥1, los llamados SAINV11, SAINV12 SAINV21 presentan notables ventajas sobre el anterior y obviamente resultan m´as econ´omicos que el hecho de construir un Precondicionador diferente para cada valor distinto del par´ametro (nombrado como Full-SAINV). Conviene destacar que, entre todos ellos, el que proporciona los mejores resultados es el SAINV11. En lo que respecta a los Precondicionadores ICHOL, ocurre algo similar. En general, para peque˜nos valores de ε, bastar´ıa con construir un ´unico Precondicionador, para un determinado valor inicial del par´ametro, el denominado ICHOL(Aε0); pero, para valores altos de ε, el ICHOLDy el ICHOLNtienen el mejor comportamiento. En especial el ICHOLNpresenta los mejores resultados de todos a partir de valores de ε≥102, para los que ni el construir un Precondicionador diferente para cada ε(α) (Full-ICHOL), dar´ıa resultado, puesto que as´ı no se alcanza la convergencia. Esto se debe a que, para valores altos de ε, al realizar la factorizaci´on
CONCLUSIONES 124 incompleta de la matriz Aε, se pierde su positividad, con lo cual la aplicaci´on del Gradiente Conjugado resulta inestable. Sin embargo esto no ocurre con los nuevos ICHOL, el D y el N. En los alrededores del valor unitario del par´ametro, los experimentos realizados no permiten obtener una conclusi´on definitiva acerca de cual ser´ıa la mejor estrategia. Probablemente la re-computarizaci´on para cada valor de εsea la elecci´on m´as fiable. Los Precondicionadores basados en al algoritmo SAINV no son tan eficientes como los que se fundamentan en la Factorizaci´on Incompleta de Cholesky, diferencia que se acent´ua, a favor de estos ´ultimos, a medida que crece el n´umero de ecuaciones del Sistema. Por ello, puede concluirse que, de todos los experimentados, el Precondicionador denominado ICHOLNes el que conduce a los mejores resultados. Teniendo en cuenta lo anterior, tambi´en puede afirmarse, que el Precondicionador ICHOLNconstituye la mejor opci´on, a la hora de optar por seleccionar, de forma autom´atica, el valor ´optimo de un par´ametro, para un Campo de Viento determinado, puesto que, hace posible utilizar para ello la potente herramienta de los Algoritmos Gen´eticos, al permitir afrontar con rapidez la gran cantidad de veces que se debe resolver el Sistema, para una amplia gama de valores del par´ametro ε. Por lo que respecta a la influencia de la Reordenaci´on, previa al Precondicionamiento, podemos concluir que, para grandes sistemas de ecuaciones, un orden adecuado mejora la eficacia del m´etodo del Gradiente Conjugado Precondicionado, ya que con ello se producen Precondicionadores con mejores cualidades que permiten reducir el n´umero de pasos para alcanzar la convergencia. Por tanto, con una adecuada Reordenaci´on puede mejorarse la eficacia de los Precondicionamientos tanto con Inversas Aproximadas, como con Factorizaciones Incompletas de Cholesky. Analizando conjuntamente, no s´olo el nuevo coste de las iteraciones una vez Reordenado el Sistema, sino a˜nadiendo tambi´en los tiempos requeridos para su implantaci´on, se comprueba que el algoritmo de Reordenaci´on m´as efectivo es el Cuthill-McKee Inverso y en particular asociado con la Factorizaci´on Incompleta
LINEAS FUTUTRAS 125 de Cholesky. Sin embargo el algoritmo Multicoloring no aporta ninguna mejora a la ejecuci´on de los m´etodos iterativos precondicionados. Queda pues patente, a modo de resumen, que el uso de los Precondicionadores basados en la Factorizaci´on Incompleta de Cholesky es una herramienta eficaz para mejorar la convergencia del algoritmo del Gradiente Conjugado Precondicionado, en el caso de los Sistemas Lineales de Ecuaciones Variables, con matrices Sim´etricas Definidas Positivas; en especial el que se ha denominado como ICHOLN, debido a su inferior coste computacional. Al menos hay una amplia gama de valores de εpara los cuales con dichos Precondicionadores se consiguen convergencias m´as r´apidas, que con el uso de los Precondicionadores obtenidos expl´ıcitamente con la Aproximada Inversa (SAINV) y, por supuesto, mejorando los resultados que se puedan conseguir con un re-computarizaci´on completa para cada valor distinto del par´ametro. 10.2. LINEAS FUTUTRAS Con este trabajo se abren varias l´ıneas futuras que admiten ser estudiadas en profundidad: En primer lugar, ser´ıa interesante comprobar el comportamiento del algoritmo del Gradiente Conjugado Precondicionado utilizando tambi´en Precondicionadores Variables, tipos SAINV e ICHOL, similares a los aqu´ı propuestos, pero aplic´andolos sobre Sistemas Lineales de Matrices Variables Sim´etricas originados en la resoluci´on de otros tipos de problemas diferentes a la Modelizaci´on de Campos de Viento. As´ı como, tambi´en ser´ıa de inter´es, estudiar los resultados a conseguir con dicho algoritmo al aplicarlo a Matrices Variables Sim´etricas en general, pero utilizando Precondicionadores Variables distintos de los SAINV e ICHOL aqu´ı descritos.
LINEAS FUTUTRAS 126 Otras posibles l´ıneas de inter´es, ser´ıan las que surjan al afrontar el estudio de la eficacia de nuevos Precondicionadores Variables, distintos a los SAINV e ICHOL, apropiados para Matrices Variables No Sim´etricas, que se originen al abordar modelizaciones diferentes a las de los Campos de Viento y tener que recurrir a algoritmos distintos al Gradiente Conjugado, pero basados en los Subespacios de Krylov, como pueden ser otros m´etodos de ortogonalizaci´on, como el GMRES o m´etodos de biortogonalizaci´on, como el Bi-CGSTAB y los QMR, TFQMR y QMRGCSTAB. Finalmente, debe se˜nalarse, que incluso ser´a interesante ampliar los estudios indicados utilizando otras T´ecnicas de Reordenaci´on, no solo basadas exclusivamente en la posici´on de los elementos de la matriz, sino tambi´en, que tuviesen en cuenta la influencia de los valores num´ericos de dichos elementos.
Bibliograf´ıa [1] L. Adams. m-Step preconditioned conjugate gradient methods. SIAM J. Sci. Stat. Comput. 6,2, 453–463 (1985). [2] P. Almeida. “Resoluci´on directa de sistemas sparse por grafos.” Tesis Doctoral, Universidad de Las Palmas de Gran Canaria (1989). [3] W.E. Arnoldi. The Principle of Minimized Iteration in the Solution of the Matrix Eingenvalue Problem. Quart. Appl. Math. 9, 17–29 (1951). [4] S.F. Ashby, T.A. Manteuffel y P.E. Saylor. A taxonomy for conjugate gradient methods. SIAM J. Numer. Anal. 27, 1542–1568 (1990). [5] O. Axelsson. A Restarted Version of a Generalized Preconditioned Conjugate Gradient Method. Comunications in Applied Numerical Methods 4, 521–530 (1988). [6] O. Axelsson. “Iterative Solution Methods.” Cambridge University Press (1996). [7] A.F. de Baas. “Modelling of Atmospheric Flow Fields.”, cap´ıtulo Scaling Parameters and their Estimation, p´aginas 87–102. World Sci. Singapore (1996). [8] J.C. Barnard, H.L. Wegley y T.R. Hiester. Improving the performance of mass consistent numerical models using optimization. J. Climate. Appl. Meteorol. 26, 675–686 (1987). [9] R. Barret, M. Berry, T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine y H.A. Van der Vorst.
Bibliograf´ıa 128 “Templates for the solution of linear systems: Building Blokcs for Iterative Methods.” SIAM, Philadelphia (1994). [10] M. Benzi. Preconditioning Techniques for Large Linear Systems: A Survey. Journal of Computational Physics 182, 418–477 (2002). [11] M. Benzi y D. Bertaccini. Approximate inverse preconditioning for shifted linear systems. BIT Num. Math. 43, 231–244 (2003). [12] M. Benzi, J.K. Cullum y M. Tuma. Robust approximate inverse preconditioning for the conjugate gradient mathod. SIAM J. Sci. Comput. 22, 1318–1332 (2000). [13] R. Boubel, D. Fox, D. Turner y A. Stern. “Fundamentals of Air Pollution”. Academic Press, San Diego (1994). [14] J.A. Businger y S.P.S. Arya. Heights of the mixed layer in the stably stratified planetary boundary layer. Adv. Geophys 18A, 73–92 (1974). [15] T.F. Chan, E. Gallopoulos, V. Simonsini, T. Szeto y C.H. Tong. A Quasi-Minimal residual variant of the BI-CGSTAB algorithm for nonsymmetric systems. SIAM J. Sci. Comput. 15,2, 338–347 (1994). [16] C. Conde y G. Winter. “M´etodos y algoritmos b´asicos de ´algebra num´erica.” Editorial Revert´e, Barcelona (1990). [17] E. H. Cuthill y J.M. Mckee. Reducing the Bandwidth of Sparse Symmetric Matrices. En “Proc. 24th National Conference of the Association for Computing Machinery”, p´aginas 157–172. Brondon Press, New Jersey, U.S.A. (1969). [18] C.G. Davis, SS. Bunker y J.P. Mutschlecner. Atmospheric Transport Models for Complex Terrain. J. Climate. Appl. Meteorol. 23(2), 235– 238 (1984). [19] M. Dikerson. A mass-consistent atmosferic flux model for regions with complex terrain. J. Appl. Meteor. 17, 241–253 (1978).
Bibliograf´ıa 129 [20] Q.V. Dinh, V. Mantel, J. Periaux y B. Stoufslet. “Contribution to problems T4 and T6 finit element GMRES and conjugate gradient solvers.” Informe T´ecnico. Dassault Aviation. (1993). [21] S. Douglas y R. Kessler. “User’s guide to the Diagnostic Wind Model (Version 1.0)”. System Applications, Inc., San Rafael. California (2009). [22] L.C. Dutto. The Effect of Ordering on Preconditioned GMRES Algorithm. Int. Jour. Num, Meth. Eng. 36, 457–497 (1993). [23] L. Elsgoltz. “Ecuaciones Diferenciales y C´alculo Variacional”. Editorial MIR. Mosc´u (1969). [24] J.M. Escobar y R. Montenegro. Several aspects of three-dimensional Delaunay triangulation. Advances in Engineering Software 1/2(27), 27–39 (1996). [25] J.M. Escobar, E. Rodr´ ıguez, R. Montenegro, G. Montero y J.M. Gonz´ alez-Yuste. Simultaneous untangling and smoothing of tetrahedral meshes. Comput. Methods Appl. Mech. Engrg. 192, 2775–2787 (2003). [26] L. Ferragut, R. Montenegro y A. Plaza. Efficient refinement/derefinement algorithm of nested meshes to solve evolution problems. Comm. Num. Meth. Eng. 10, 403–412 (1994). [27] R. Fletcher. Conjugate Gradient Methods for Indefite Systems. Lectures Notes in Math. 506, 73–89 (1976). [28] E. Fl´ orez, M.D. Garcia, A. Su´ arez y H. Sarmiento. The Effect of Ordering on the Convergence of the Conjugate Gradient Method for Solving Preconditioned Shifted Linear Systems. En B.H.V. Topping, G. Montero y R. Montenegro, editores,“Proceeding of The Fifth International Conference on Engineering Computational Techonology. Las Palmas de G. C.”, p´aginas 191–192. Civil-Comp Press, Stirlingshire, U.K. (2006).
Bibliograf´ıa 130 [29] E. Fl´ orez, H. Sarmiento, M.D. Garcia, A. Su´ arez y G. Montero. Incomplete factorisation for preconditioning shifted linear systems arising from a parameter estimation problem in wind modelling. En B.H.V. Topping, editor, “Proceedings of the Sixth Int. Conference on Engineering Computational Techonology. Atenas”. Civil-Comp Press, Stirlingshire, U.K. (2008). [30] E. Fl´ orez V´ azquez. “Construcci´on de inversas aproximadas tipo ”sparse”basada en la proyecci´on ortogonal de Frobenius para el precondicionamiento de sistemas de ecuaciones no sim´etricos.” Tesis Doctoral, Universidad de Las Palmas de G. C. (2003). [31] R.W. Freund. A transpose-free quasi-minimal residual algorithm for nonHermitian linear systems. SIAM J. Sci. Comput. 14, 470–482 (1993). [32] R.W. Freund y N.M. Nachtigal. Qmr: a quasi-minimal residual method for non-Hermitian linear systems. Numerische Math. 60, 315–339 (1991). [33] R.W. Freund y N.M. Nachtigal. An implementation of the QMR method based on coupled two-term recurrences. SIAM J. Sci. Comp. 15,2, 313–337 (1994). [34] M. Gal´ an. “Avances en el M´etodo de Residuo M´ınimo Generalizado (algoritmo GMRES), su Desarrollo en ANSI-C con Algoritmos de Paralelizaci´on y Vectorizaci´on, y sus Aplicaciones al M´etodo de los Elementos Finitos.” Tesis Doctoral, Universidad de Las Palmas de Gran Canaria (1994). [35] M. Gal´ an, G. Montero y G. Winter. A direct solver for the least square problem arising from GMRES(k). Com. Num. Meth. Eng. 10, 743– 749 (1994). [36] M.D. Garc´ ıa, E. Fl´ orez, A. Su´ arez, L. Gonz´ alez y G. Montero. New implementation of QMR-type algorithms. Computers and Structures 83, 2414–2422 (2005).
Bibliograf´ıa 131 [37] M.D. Garc´ ıa Le´ on. “Estrategias para la resoluci´on de grandes sistemas de ecuaciones lineales. M´etodos de Cuasi-M´ınimo Residuo Modificados.” Tesis Doctoral, Universidad de Las Palmas de G. C. (2003). [38] P. Geai. “Methode d’interpolation et de reconstitution tridimensionelle d’un champ de vent: le code d’analyse objetive MINERVE”. InformeT´ecnico. Electricit`e de France (1985). [39] P. Geai. “Reconstitution tridimensionnelle d’un champ de vent dans un domaine a’topographie complexe a partir de meseures in situ.” InformeT´ecnico. DER/HE/34-87.05, EDF, Chatou. France (1987). [40] I.M. Gelfand y S.V. Fomin. “Calculus of Variations”. Prentice-Hall, INC (1963). [41] A. George. Computer Implementation of the Finite Element Method. Report Stan CS p´aginas 71–208 (1971). [42] A. George y J.W. Liu. The Evolution of the Minimum Degree Ordering Algorithms. SIAM Rev. 31, 1–19 (1989). [43] G.H. Golub y G.A. Meurant. “R´esolution num´erique des grands syst`emes lin´eaires.” Editions Eyrolles, Par´ıs. (1983). [44] J.M. Gonz´ alez-Yuste. “Un Algoritmo de Refinamiento/Derefinamiento Local par Mallas de Tetraedros.” Tesis Doctoral, Universidad de Las Palmas de Gran Canaria (2004). [45] J.M. Gonz´ alez-Yuste, R. Montenegro, J. Escobar, G. Montero y E. Rodr´ ıguez. Local refinement of 3-D triangulations using objectoriented methods. Adv. in Eng. Softw. (2003). [46] A. Greenbaum. “Iterative Methods for Solving Linear Systems”. SIAM, Philadelphia (1997). [47] M. Grote y H. Simon. “Parallel Preconditioning and Approximate Inverses on the Connection Machine.” Informe T´ecnico. NASA Contract No. NAS2-12961 (1986).