Electromagnetismo Numerico

Resolución numérica de campos electromagnéticos Julián Barquı́n ii Prefacio El propósito de este escrito es servir de apoyo al laboratorio de métodos numéricos de la asignatura de Electrotecnia I de la Escuela Técnica Superior de Ingenierı́a Industrial de la Universidad Pontificia Comillas. Por otra parte, se ha intentado incluir no solamente el material que se imparte en esta asignatura, sino también material adicional que pueda ser de utilidad en cursos superiores. Se ha diferenciado estos dos tipos de materiales señalándoles como “básico” y “avanzado”. En este texto no se trata de los aspectos fı́sicos de las ecuaciones tratadas, lo que es, por supuesto, la preocupación principal en la asignatura de Electrotecnia I. En cambio, si que se incluye un capı́tulo (el segundo) sobre la resolución numérica de sistemas lineales. Aunque este tema debiera serle familiar a los alumnos al finalizar el segundo curso, no lo es todavı́a cuando comienza el laboratorio de métodos numéricos, por lo que unas ligeras nociones preliminares son deseables. Por otra parte, creo que es útil que tengan una pequeña descripción de los algoritmos útiles en la resolución de las ecuaciones de campo. Tras este capı́tulo, el tercero trata de la resolución de la ecuación de Poisson mediante diferencias finitas, y el cuarto mediante elementos finitos. En el quinto se estudian otras ecuaciones elı́pticas de interés en electromagnetismo. El sexto capı́tulo aborda el tratamiento de la ecuación de difusión, y el último de su resolución en régimen senoidal. No se han tratado ecuaciones hiperbólicas. La razón es que estas ecuaciones no se tratan en la asignatura con el detenimiento de las elı́pticas o las parabólicas, tanto por el (muy relativo) menor interés técnico, como porque se estudian con menos intensidad en segundo curso, debido, entre otras cosas, a que los alumnos todavı́a no han visto análisis de Fourier. Madrid, octubre de 1996. iii iv Índice general 1 Introducción 1 2 Resolución de sistemas lineales 2.1 Introducción . . . . . . . . . . . . . . . . . . . . . 2.2 Método directo: la factorización LU . . . . . . . . 2.2.1 Factorización LU de matrices llenas . . . 2.2.2 Factorización LU de matrices ralas . . . . 2.3 Métodos iterativos . . . . . . . . . . . . . . . . . 2.3.1 Los métodos de Jacobi y Gauss-Seidel . . 2.3.2 El método de la sobrerrelajación sucesiva 2.4 Métodos semidirectos: el gradiente conjugado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 3 5 5 9 11 11 17 19 3 Ecuaciones de Laplace y Poisson: diferencias finitas 3.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Preliminares matemáticos . . . . . . . . . . . . . . . . . . . 3.3 Construcción de la malla y planteamiento de las ecuaciones 3.4 El principio del máximo . . . . . . . . . . . . . . . . . . . . 3.5 Errores y convergencia . . . . . . . . . . . . . . . . . . . . . 3.5.1 Anexo: prueba del teorema . . . . . . . . . . . . . . 3.6 Problemas con los contornos . . . . . . . . . . . . . . . . . . 3.6.1 Condiciones de von Neumann no nulas . . . . . . . . 3.6.2 Tratamiento de contornos curvos . . . . . . . . . . . 3.7 Diferencias finitas en otros sistemas coordenados . . . . . . 3.7.1 Diferencias finitas en coordenadas polares . . . . . . 3.7.2 Diferencias finitas con simetrı́a axial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 27 28 29 33 34 36 38 38 39 40 41 41 4 Ecuaciones de Laplace y Poisson: elementos finitos 4.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . 4.2 Forma variacional de la ecuación de Poisson . . . . . 4.3 Aproximación funcional . . . . . . . . . . . . . . . . 4.3.1 Un ejemplo . . . . . . . . . . . . . . . . . . . 4.4 Nomenclatura y propiedades . . . . . . . . . . . . . . 4.5 El espacio de funciones admisibles . . . . . . . . . . 4.6 Elementos finitos triangulares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 45 45 48 51 54 55 58 v . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi ÍNDICE GENERAL 4.7 4.8 4.9 Otros elementos . . . . . . . . . . . . . . . . . . . . . . . . . . . . Errores y convergencia . . . . . . . . . . . . . . . . . . . . . . . . Elementos finitos y diferencias finitas . . . . . . . . . . . . . . . . 64 66 69 5 Otros problemas estáticos 5.1 Introducción . . . . . . . . . . . . . . . . . . . . . . 5.2 Corrientes estacionarias . . . . . . . . . . . . . . . 5.3 Magnetostática en medios lineales . . . . . . . . . 5.4 Magnetostática no lineal . . . . . . . . . . . . . . . 5.4.1 Planteamiento de la forma variacional . . . 5.4.2 Resolución numérica de la forma variacional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 71 71 72 76 76 78 6 La ecuación de difusión 6.1 Introducción . . . . . . . . . . . . . . . . . 6.2 La ecuación de difusión en 1 dimensión . . 6.3 Resolución de las ecuaciones . . . . . . . . 6.3.1 Método explı́cito clásico . . . . . . 6.3.2 Método implı́cito clásico . . . . . . 6.3.3 Método de Crank-Nicholson . . . . 6.4 Estabilidad de la solución . . . . . . . . . 6.5 Convergencia de los algoritmos . . . . . . 6.6 Elementos finitos y la ecuación de difusión . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 83 84 86 87 88 89 91 94 97 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Problemas armónicos 101 7.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 7.2 Notación compleja . . . . . . . . . . . . . . . . . . . . . . . . . . 102 7.3 Difusión magnética e histéresis . . . . . . . . . . . . . . . . . . . 104 Capı́tulo 1 Introducción Hace ya más de 200 años, Leonard Euler, de quien dicen que jamás dejaba de calcular1 , decidió abordar el problema del flujo del agua. El conocı́a las ecuaciones de Newton, y sabı́a como aplicarlas a partı́culas o sólidos rı́gidos. En estos casos, el estado del sistema fı́sico está descrito por un número finito de cantidades: las posiciones y las velocidades de las partı́culas, o la posición y velocidad del centro de masas junto con la orientación del cuerpo y su velocidad angular. No ası́ con el agua: para definir su estado se necesita saber la velocidad de cada partı́cula. Matemáticamente, su estado se puede representar mediante el campo de velocidades v(x, y, z). Euler logró derivar una ecuación que relacionaba el cambio de la velocidad con el tiempo con su cambio en el espacio, una ecuación en derivadas parciales, o una ecuación de campo. 100 años después, James Maxwell logró dar forma matemática a las ideas de Michel Farady sobre las leyes del campo electromagnético. De nuevo se encontró con ecuaciones en derivadas parciales, que en teorı́a permitirı́an calcular los campos eléctrico E(x, y, z) y magnético B(x, y, z). Este fue, sin duda alguna2 , el hecho más importante del siglo XIX, el comienzo de la era eléctrica. Durante los años siguientes, el final del siglo XIX y el siglo XX, los ingenieros aplicarı́an estas ecuaciones en un número creciente de dispositivos. Sin embargo, es curioso de que, salvo en raras ocasiones, no las resolvı́an: las simplificaban. La razón, por supuesto, es que no se sabı́a resolverlas, salvo en casos de una extrema idealización. Hoy la situación es distinta. La razón es que la aparición de los ordenadores (lejanos descendientes de Farady y Maxwell, por otra parte) permiten la resolución de grandes sistemas de ecuaciones algebraicas. Este escrito trata, por tanto, de la aproximación de ecuaciones en derivadas parciales por ecuaciones algebraicas. En electromagnetismo, casi todas las ecuaciones de interés tienen un mismo aspecto: son ecuaciones de segundo orden. Por ejemplo, en dos dimensiones 1 De 2Y hecho, murió rico. que digan lo que quieran. 1 2 CAPÍTULO 1. INTRODUCCIÓN toman la forma: a ∂2φ ∂2φ ∂2φ ∂φ ∂φ + b + c +d +e + fφ + g = 0 2 2 ∂x ∂x∂y ∂y ∂x ∂y (1.1) donde a, b, c, d, e, f y g son funciones conocidas de x e y, y φ es la función (el campo) a resolver. El caracter cualitativo de la solución depende, sin embargo, de los coeficientes de los términos de segundo orden: a, b y c. En concreto: • b2 − 4ac < 0. Se llaman ecuaciones elı́pticas. Describen situaciones de equilibrio, como el potencial eléctrico en problemas electrostáticos, o las tensiones mecánicas en una viga. • b2 −4ac = 0. Son ecuaciones parabólicas. Describen fenómenos de difusión, tales como la conducción del calor, o el efecto pelicular en conductores. • b2 − 4ac > 0. Se trata de las ecuaciones hiperbólicas. Describen campos que se propagan a una velocidad dada, tales como el sonido o la luz. Es curioso constatar como comportamientos tan diferentes dependen de modificaciones aparentemente tan pequeñas en la ecuación. De todas formas, aquı́ se trataran tan solo los dos primeros tipos de ecuaciones. Capı́tulo 2 Resolución de sistemas lineales ▽ Básico 2.1 Introducción La mayor parte de este escrito trata sobre la forma en que las ecuaciones de campo son aproximadas por sistemas de ecuaciones algebraicas. Estos sistemas son muy a menudo sistemas lineales. Cuando no son sistemas lineales, su resolución suele requerir la de varios sistemas lineales. Por todo ello, las técnicas numéricas de resolución de estos sistemas son de una importancia crucial. El tratamiento sistemático de estos procedimientos hace muy conveniente, de hecho casi imprescindible, el que estos sistemas lineales se planteen de forma matricial. Ası́ pues, se puede reescribir el párrafo anterior diciendo que este capı́tulo trata de la resolución del sistema Ax = b (2.1) A, la matriz del sistema, es una matriz cuadrada que depende de la ecuación a resolver y, naturalmente, del método escogido. El término independiente b depende del valor de las condiciones de contorno y de la densidad de carga (o de fuerzas volumétricas, o de generación de calor) dentro del dominio. x es, naturalmente, la incógnita, los valores del campo. Una caracterı́stica muy importante de la matriz A es que casi todos sus elementos son nulos. Por ejemplo, la figura 2.1 muestra la malla de elementos finitos utilizada para analizar el flujo alrededor del ala de un avión y la matriz A que tiene asociada. Se ha dibujado un punto en la posición de todos los elementos de la matriz cuyo valor es distinto de cero. Obsérvese que exigua minorı́a constituyen. Una matriz con esta caracterı́stica es llamada rala o cuasivacı́a1 . Esta caracterı́stica ofrece la posibilidad de tratar matrices de mucha mayor dimensión de 1 En inglés, sparse. Cierta similitud fonética parece haber llevado en ocasiones a traducirlo (erróneamente) como dispersa. 3 4 CAPÍTULO 2. RESOLUCIÓN DE SISTEMAS LINEALES The adjacency matrix. 0 1 500 0.9 0.8 1000 0.7 1500 0.6 2000 0.5 2500 0.4 3000 0.3 0.2 3500 0.1 4000 0 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 500 1000 1500 2000 2500 nz = 28831 3000 Figura 2.1: Ala de avión, su malla y su matriz asociadas. lo que serı́a de otra forma posible. Por ejemplo, es posible almacenar solamente los elementos que sean distintos de cero, junto con sus posiciones. En el caso de la matriz del ejemplo, de dimensión 1700 y con 3500 elementos no nulos, no es necesario almacenar 1700 × 1700 = 3, 500, 000 posiciones de memoria, sino que bastan con 3500 × 3 = 10,500 posiciones (1 para el valor del elemento y dos para su posición en la matriz). De hecho, no es solamente la necesidad de memoria la que se puede reducir de una forma muy marcada, sino también el tiempo de cómputo requerido. Son ası́ las técnicas de resolución de sistemas ralos las que hacen posibles el tratamiento de problemas que no sean de dimensión muy reducida. Existen tres tipos de métodos fundamentales para la resolución de sistemas lineales: los directos, los iterativos y los semidirectos. Los métodos directos resuelven las ecuaciones en un número finito de pasos. Los métodos iterativos van mejorando una estimación inicial de la solución (a menudo x = 0), aunque sin llegar jamás a la solución exacta. Tienen, sin embargo, menores requerimientos de memoria que los métodos directos. Los semidirectos pueden, en principio, resolver también el sistema en un número finito de pasos, pero en su planteamiento y análisis son más parecidos a los iterativos. El resto del capı́tulo explicará algún método, o algunos métodos, pertenecientes a cada categorı́a. Son, de todas formas, los más usados en la resolución numérica de campos. 3500 4000 2.2. MÉTODO DIRECTO: LA FACTORIZACIÓN LU 2.2 5 Método directo: la factorización LU △ Básico Un método directo es un método que alcanza la solución tras un número finito ▽ Avanzado de pasos. De todos ellos, el más empleado en el problema que aquı́ se trata es una generalización del método de eliminación de Gauss conocida como factorización LU. En el resto de la sección se explica la factorización de Gauss en el contexto de las matrices llenas. Después se exponen las modificaciones necesarias para su aplicación a matrices ralas. 2.2.1 Factorización LU de matrices llenas La eliminación de Gauss (y la factorización LU) se basan en una observación aparentemente trivial: es muy fácil resolver un sistema lineal cuando la matriz A es triangular. En efecto, si se tiene un sistema como      2 0 0 0 3 x1  3 1 0 0   x2   −1       (2.2)  −4 0 3 0   x3  =  5  0 x4 7 1 0 −4 se puede resolver primero x1 de la primera ecuación (la primera fila), conocido x1 se resuelve x2 de la segunda, conocidos x1 y x2 se puede resolver x3 de la tercera, y ası́ sucesivamente. Evidentemente, se puede proceder ası́ porque la matriz del sistema tiene elementos no nulos solamente por debajo de la diagonal principal (los elementos aij , con i ≥ j). Si la matriz fuera triangular superior (elementos no nulos solamente por encima de la diagonal superior) se podrı́a proceder de la misma forma, salvo que se tendrı́a que comenzar por la última componente de x, en vez de por la primera (¡compruébese!). Esta observación no dejarı́a de ser trivial sino fuera por que toda matriz A no singular se puede escribir como el producto de una matriz triangular inferior por una superior: A = LU (2.3) Es costumbre denotar al factor triangular inferior por L (del inglés lower) y al superior por U (de upper). La factorización LU es, esencialmente, el procedimiento de cálculo de estos dos factores a partir de la matriz A. Supóngase entonces que se ha logrado calcular L y U . Ahora la resolución del sistema Ax = b es inmediata. En efecto, este sistema es equivalente a LU x = b (2.4) Consideremos también el sistema Ly = b (2.5) Como este es un sistema triangular inferior, es fácil de resolver. La solución es y = L−1 b. Además, de la ecuación (2.4) 6 CAPÍTULO 2. RESOLUCIÓN DE SISTEMAS LINEALES U x = L−1 b = y (2.6) ¡Pero esto vuelve a ser una ecuación con matriz triangular!. Como y es conocida (se ha resuelto (2.5)), es ahora fácil calcular x. En resumen, se han resuelto sucesivamente los sistemas Ly = b Ux = y (2.7) (2.8) El primer sistema es triangular inferior, luego las incógnitas se resuelven iendo de la primera a la última. El segundo es triangular superior, en el que las incógnitas se resuelven de la última a la primera. Es por esto que el procedimiento completo se conoce por la substitución hacia delante y hacia atrás. El problema que queda es, naturalmente, el cálculo de los factores L y U . Es esta la tarea, con mucho, más complicada y costosa del procedimiento. La técnica consiste en ir eliminado sucesivamente los elementos de A que se encuentren por debajo de la diagonal. Quizá la forma más sencilla de explicarlo es mediante un ejemplo. Supóngase ası́ que se desea factorizar la matriz   3 1 0 1  1 4 2 0   A= (2.9)  0 2 5 2  1 0 2 6 Considérese el elemento (2,1) (segunda fila, primera columna). Este elemento (un 1) está por debajo de la diagonal principal e impide, por tanto, que la matriz A sea diagonal superior. Para eliminar este elemento, se puede premultiplicar por una triangular inferior:  1  −l2,1   0 0  3 0 0 0  1 1 0 0   0 1 0  0 0 0 1 1 El término l2,1 matriz, es decir  0 1 2 −l2,1   5 2  2 6 (2.10) se escoge de manera que anule al término (2, 1) de la nueva 1 4 2 0 0 2 5 2   3 1  1 − l2,1 3 0  = 0 2   1 6 1 − 3l2,1 = 0 ⇒ l2,1 = 1 4 − l2,1 2 0 1 3 (2.11) Luego,  1  −1  3  0 0  0 0 0 3  1 1 0 0   0 1 0  0 1 0 0 1 1 4 2 0 0 2 5 2   3 1  0 0  = 2   0 6 1 1 11 3 2 0  0 1 2 − 13   5 2  2 6 (2.12) 7 2.2. MÉTODO DIRECTO: LA FACTORIZACIÓN LU Ahora, obsérvese que  1  −1  3  0 0 −1  0 0 0 1  1 1 0 0   = 3  0 0 1 0  0 0 1 0  0 0 0 1 0 0   0 1 0  0 0 1 (2.13) En realidad, es en general cierto para todas las matrices triangulares inferiores formadas por una diagonal de 1’s y un solo elemento adicional debajo de la diagonal, que su inversa es ella misma salvo por el elemento adicional que cambia signo. Se sigue ası́ de (2.12) que:  3  1 A=  0 1 1 4 2 0 0 2 5 2   1 1  1 0  = 3 2   0 6 0  3 0 0 0  0 1 0 0   0 1 0  0 1 0 0 1 1 11 3 2 0  0 1 2 − 13   = L1 A1 (2.14) 5 2  2 6 La matriz A1 tiene todavı́a un elemento en la primera columna por debajo de la diagonal distinto de 0: el 1 de la cuarta fila. Procediendo de manera análoga a la anterior se tiene:    1 0 1  0 2 − 13  = 5 2   0 1 2 6 3  3 1 0 0 0  0 11 1 0 0  3  0 1 0  0 2 0 0 1 0 − 31   0 1 1  0 2 − 13  = 5 2   0 0 2 17 3 0 1  0 1 11 2 − 31  3  = L2 A2 5 2  2 0 2 17 3 (2.15) Nótese que la matriz L2 no afecta más que a la cuarta fila de A2 . Por tanto, no puede cambiar los ceros de las filas segunda y tercera a algo no nulo. Esto es, naturalmente, lo que se pretende. En todo caso, ya no hay elementos a eliminar en la primera columna de A2 . En la segunda columna el elemento (3, 2) es distinto de cero (2). Ası́, se tiene: 3  0 A1 =   0 1  1 3 1  0 11 3 A2 =   0 2 0 − 31 6 11 0  3 1 0 0  0 11 0 0  3  1 0  0 0 0 1 0 − 13 0 2 43 11 2  1 − 31   24  = L3 A3 11 17 3 (2.16) Hay que notar ahora que la matriz L3 no afecta a la primera columna de A3 . Es por esto que se ha comenzado por eliminar la primera columna, y por lo que una vez que se elimine la segunda se comenzará por la tercera. Las matrices L’s que afectan a cada columna no pueden modificar las columnas anteriores. En A3 el elemento (4, 2) es distinto de cero. Por tanto, se elimina: 8 CAPÍTULO 2. RESOLUCIÓN DE SISTEMAS LINEALES  3 1  0 11 3 A3 =   0 0 0 − 31 0 2 43 11 2   1 1  0 − 31    24  =  0 11 17 0 3  3 0 0  0 0 0   1 0  0 0 1 0 0 1 0 1 11 1  1 − 13   24  = L4 A4 0 2 11 3 43 11 24 11 0 0 11 62 11 (2.17) Finalmente, solo queda eliminar el término (4, 3) de A4 . De aquı́ resulta:  3  0 A4 =   0 0 1 11 3 0 0 0 2 43 11 24 11   1 1  0 − 13    24  =  0 11 62 0 11 0 1 0 0 0 0 1 24 43  3 0  0 0   0  0 1 0 1 11 3 0 0  1 − 31   24  = L5 U 0 2 43 11 11 2048 473 0 (2.18) U es una matriz triangular superior. Ahora, de todo lo anterior se sigue: A = L1 A1 = L1 L2 A2 = L1 L2 L3 A3 = L1 L2 L3 L4 A4 = L1 L2 L3 L4 L5 U (2.19) Por otra parte, el producto de matrices triangulares inferiores es triangular inferior. Luego L1 L2 L3 L4 L5 es una triangular inferior que se llamará L. Es, además, muy fácil de calcular. En efecto   1 0 0 0  1 1 0 0  3  (2.20) L=  0 6 1 0  11 1 24 1 1 3 11 43 Es decir, el efecto de cada Li es añadir su término fuera de la diagonal. Esto es válido en general, para matrices Li como las que resultan del proceso de factorización (una diagonal de 1’s más un término fuera de diagonal). Nótese que para calcular cada matriz Li hay que resolver la ecuación −li,(k,l) ai−1,(l,l) + ai−1,(k,l) = 0 (2.21) siendo (k, l) el término a eliminar. Para resolver esta ecuación en (li,(k,l) es preciso que el término de la diagonal ai−1,(l,l) 6= 0. Este término es llamado el pivote. Caso de ser nulo, se puede intentar permutar la fila l-ésima por otra situada más abajo. Esto equivale a reordenar las variables del sistema lineal. Si no hubiera ninguna permutación que resolviera el problema (que todas dieran pivotes nulos) es que la matriz es singular y no hay, por tanto, solución única. Por otra parte, en el ejemplo que se estudia, se puede comprobar que  3  0 U =  0 0 1 11 3 0 0 0 2 43 11 0   3 1  0 − 31    24  =  0 11 2048 0 473 0 11 3 0 0 0 0 43 11 0 0 0 0 2048 473  1  0   0 0 1 3 1 0 0 0 6 11 1 0 1 3 1 11 24 43 1    = DLT  (2.22) 9 2.2. MÉTODO DIRECTO: LA FACTORIZACIÓN LU Luego A = LU = LDLT (2.23) Esta expresión es conocida como la factorización triple LDLT . Es posible llevarla a cabo si la matriz A es simétrica (¡pruébese!). Esta simetrı́a se presenta a menudo en problemas de campos, y se puede explotar para mejorar el procedimiento de factorización LU arriba explicado. 2.2.2 Factorización LU de matrices ralas La factorización LU expuesta arriba tiene el inconveniente de que las operaciones se realizan sobre matrices llenas, es decir, sin tener en cuenta que la mayor parte de los elementos de la matriz son nulos. Esencialmente, la factorización LU con matrices ralas procede de la misma forma que la factorización con matrices llenas. Las diferencias que existen son las siguientes: 1. Solamente se almacenan los elementos no nulos de las matrices, junto con su localización. Ası́, por ejemplo, en Matlab, la matriz:  2  1   0 0 1 2 1 0 0 1 2 1  0 0   1  2 se almacena en tres vectores, un primero que almacena los valores no nulos, y otros dos que almacenan las posiciones (fila y columna) donde se encuentran: ( ( ( 2 1 1 1 1 2 1 2 1 2 2 2 1 2 3 1 3 2 2 3 3 1 1 3 4 4 3 2 4 4 ) ) ) Existen otras formas de codificar la posición más eficientes que la anterior, aunque más complejas. 2. Existe una lógica que determina las operaciones y sumas que es preciso realizar, saltándose todas las operaciones que consistan en multiplicar por cero, o sumar cero. Esto es posible porque se saben los elementos que son nulos (los que no estén codificados en los vectores anteriores). En cualquier caso, es claro que interesa conseguir que los factores LU que resulten tengan tantos elementos nulos como sea posible. En este sentido, el orden que adopten las variables y ecuaciones es de crucial importancia. Por ejemplo, considérese la matriz 10 CAPÍTULO 2. RESOLUCIÓN DE SISTEMAS LINEALES  20  15 A=  6 4  15 6 4 6 0 0   0 4 0  0 0 10 Al ser simétrica, admite una descomposición LDLT . Es fácil comprobar que el factor L es   1 0 0 0  0,75 1 0 0     0,30 0,8571 1 0  0,20 0,5714 0,2264 1 Es decir, aún cuando la matriz original tenı́a bastantes ceros, el factor L tiene todos los elementos posibles (los de debajo de la diagonal) no nulos. En cambio, si se factoriza la matriz   10 0 0 4  0 4 0 6   B=  0 0 6 15  4 6 15 20 se obtiene el factor L  1  0   0 0,4 0 1 0 1,5  0 0 0 0   1 0  2,5 1 Nótese que ahora el factor L tiene tantos ceros como la matriz B. Por otra parte, las matrices A y B representan al mismo sistema lineal. En efecto, lo único que se ha hecho es cambiar el orden de ecuaciones (filas) y variables (columnas), de manera que la que antes era la primera es ahora la última, y viceversa. Por tanto, es importante ordenar adecuadamente las ecuaciones y las variables al formar la matriz del sistema lineal a resolver. “A grosso modo”, el mejor sitio para los ceros es al principio de las filas, antes de los elementos no nulos. Esto es porque al hacer la eliminación de Gauss, estos elementos ya no requieren ser eliminados. Existen una serie de algoritmos que buscan el orden óptimo en el que colocar ecuaciones y variables. Matemáticamente, este es un problema para el que todavı́a no se ha encontrado solución, aunque la experiencia demuestra que algunos de los algoritmos propuestos son muy eficaces. En el ejemplo anterior se ve como lo más eficaz es colocar la fila con más elementos no nulos al final de la matriz. Esta es la idea subyacente al algoritmo del grado mı́nimo, quizá el más popular de los algoritmos de reordenación. Por ejemplo, la figura 2.2 muestra la reordenación que este algoritmo da para la matriz de la figura 2.1. Obsérvense las estructuras en forma de flecha que se forman, que recuerdan lo obtenido en el ejemplo. 11 2.3. MÉTODOS ITERATIVOS Minimum degree 0 500 1000 1500 2000 2500 3000 3500 4000 0 500 1000 1500 2000 2500 nz = 28831 3000 3500 4000 Figura 2.2: Reordenación de la matriz. 2.3 Métodos iterativos Los métodos iterativos resuelven el sistema lineal Ax = b (2.24) calculando una serie {x0 , x1 , . . . , xk , . . .} que aproxima cada vez mejor la solución lı́m xk = x = A−1 b k→∞ (2.25) Los dos métodos iterativos de uso más frecuente en el cálculo de campos son el de Jacobi y el de Gauss-Seidel, que se estudian en la siguiente sección. Después, se expone una mejora a estos métodos, importante para obtener la convergencia de la serie a una velocidad razonable, conocida como sobrerelajación sucesiva. 2.3.1 Los métodos de Jacobi y Gauss-Seidel Escrı́base de nuevo la ecuación (2.24): Ax − b = 0 (2.26) Sumando un término P x a cada término, siendo P una matriz que se definirá más tarde, de la misma dimensión que A: (A + P )x − b = P x (2.27) 12 CAPÍTULO 2. RESOLUCIÓN DE SISTEMAS LINEALES La serie {xk } se construye dando una estimación inicial {x0 } (muy a menudo un vector con todas las componentes nulas), y calculando, a partir de aquı́, cada xk a partir de la estimación anterior xk−1 . La regla para pasar de una estimación a la siguiente es: P xk = (A + P )xk−1 − b (2.28) En principio, no hay garantı́a alguna de que la sucesión ası́ generada converja a la solución del sistema (2.24), salvo que la matriz P haya sido adecuadamente escogida. Pero si la sucesión converge, en el lı́mite se cumple la ecuación (2.27), que es equivalente a (2.24). Para determinar en que condiciones la sucesión generada por la ley (2.28) converge, es conveniente analizar el comportamiento del error ek = xk − xk−1 (2.29) Si lı́mk→∞ xk = x, entonces lı́mk→∞ ek = 0. Ahora bien, de (2.28) ek = P −1 (A + P )xk−1 − b − P −1 (A + P )xk−2 + b = P −1 (A + P )ek−1 = Gek−1 (2.30) con la matriz de error G = P −1 (A + P ) Se sigue entonces: = = = = ek Gek−1 G2 ek−2 G3 ek−3 G k e0 (2.31) Ası́ pues, lı́mk→∞ ek = 0 implica que lı́mk→∞ Gk = 0. Para ver que significa esta condición, supóngase que G tiene n autovalores distintos λi con autovectores asociados vi : Gvi = λi vi (2.32) Como los autovectores forman una base, es posible escribir e0 como una combinación lineal de los mismos: e0 = n X (2.33) ci vi i=1 Entonces Ge1 = Ge0 = n X i=1 Gci vi = n X i=1 ci λi vi (2.34) 13 2.3. MÉTODOS ITERATIVOS Iterando el procedimiento, se obtiene: Gek = n X ci λki vi (2.35) i=1 El sumatorio solamente tenderá a cero si todos los λki tienden a cero, es decir, si kλi k < 1 ∀i. Definiendo el radio espectral de la matriz G como ρ(G) = máx kλi k (2.36) i esta condición se puede reescribir como que ρ(G) < 1. En resumen, se ha de escoger una matriz P tal que 1. ρ(P −1 (A + P )) < 1 2. Los sistemas de la forma P xk = (A + P )xk−1 − b, donde la incógnita sea xk , sean fáciles de resolver. Desde el punto de vista de la primera condición, la mejor P serı́a −A. En efecto, en este caso G = −A−1 (A − A) = 0, y por tanto el error se anuları́a en la primera iteración (e1 = Ge0 = 0). El inconveniente es que, naturalmente, la matriz A no da un sistema de fácil solución: si lo fuera se resolverı́a directamente el sistema Ax = b. De todas formas, el párrafo anterior sugiere que la matriz P debiera ser parecida a −A. Los métodos de Jacobi y Gauss-Seidel se basan en tomar sus partes diagonal o triangular. Concretando, descompóngase A en A=N +D+S donde N es la parte inferior, D la diagonal y  2 1 0  1 2 1 A=  0 1 2 0 0 1 Entonces N D  0  1 =   0 0  2  0 =   0 0 (2.37) S la superior. Por ejemplo, si  0 0   (2.38) 1  2 0 0 1 0 0 0 0 1 0 2 0 0 0 0 2 0  0 0   0  0  0 0   0  2 14 CAPÍTULO 2. RESOLUCIÓN DE SISTEMAS LINEALES S Entonces, se tiene que  0  0  =  0 0 1 0 0 0 0 1 0 0  0 0   1  0 (2.39) • Jacobi hace P = −D. Por tanto, la iteración de Jacobi es −Dxk = (N + S)xk−1 − b (2.40) • Gauss-Seidel hace P = −(D + N ). Luego, la iteración es −(D + N )xk = Sxk−1 − b (2.41) Ambas fórmulas son fáciles de aplicar, pues la primera implica la resolución de un sistema cuya matriz es diagonal, y la segunda la de una matriz triangular inferior. Lo que es preciso es probar, además, que convergen, es decir que ρ(D−1 (N + S)) < 1 y que ρ((D + N )−1 S) < 12 . En general estas desigualdades son falsas. Hay, sin embargo, un caso de gran importancia en que se verifican: cuando la matriz A es de diagonal dominante. Esto significa que los valores absolutos elementos de la diagonal son mayores que la suma de los valores de los elementos fuera de la diagonal. Es decir X | Ai,i |≥ | Ai,j | (2.42) j6=i De hecho, las matrices que aparecen en el cálculo de la ecuación de Poisson, y en otras, verifican esta desigualdad. Se probará ahora esta afirmación para los dos métodos. Para ello es conveniente introducir otras medidas del “tamaño” de una matriz, análoga al radio espectrales: las normas matriciales. Para ello, considérese primero las normas (“longitudes” ) de un vector x. Una norma de un vector de n dimensiones es una aplicación de ℜn en ℜ que verifica 1. kxk > 0 si x 6= 0, kxk = 0 si x = 0 2. kαxk =| α | kxk, siendo α un escalar. 3. kx + yk ≤ kxk + kyk Existen varias de estas normas, tales como • La norma euclı́dea kxk2 = 2 Es p x(1)2 + . . . + x(n)2 (2.43) evidente que ρ(G) = ρ(−G), puesto que el radio espectral se define por los módulos de los autovalores, que se limitan a cambiar de signo cuando la matriz cambia signo. 15 2.3. MÉTODOS ITERATIVOS • La norma-0 kxk0 =| x(1) | + . . . + | x(n) | (2.44) kxk∞ = máx | x(k) | (2.45) • La norma-∞ k A cada una de estas normas vectoriales se les puede asociar una norma matricial. Por ejemplo, para la norma-∞ kAk∞ = máx kAxk∞ kxk=1 (2.46) y de forma análoga para las demás normas. Nótese que, a partir de esta definición, si kxk∞ = 1 kAxk∞ ≤ kAk∞ kxk∞ (2.47) Si kxk∞ 6= 1, esta fórmula es también válida. En efecto, se ha de verificar para el vector y = kxk1 ∞ x kAyk∞ ≤ kAk∞ kyk∞ (2.48) ya que kyk∞ = 1. Multiplicando ahora por el escalar kxk∞ se obtiene la expresión buscada. Es ahora posible demostrar que una condición suficiente (aunque no necesaria) para que un algoritmo iterativo converja es que la norma de la matriz G = P −1 (A + P ) sea menor que 1. En efecto, escribiendo la ecuación de error, y teniendo en cuenta la desigualdad anterior kek k∞ = kGek−1 k∞ ≤ kGk∞ kek−1 k∞ ≤ kGkk∞ ke0 k∞ (2.49) Si kGk∞ < 1, entonces lı́mk→∞ kek k = 0. Si kGk∞ = 1, se tiene al menos la garantı́a de que el error no puede crecer, y de hecho normalmente disminuirá. Si kGk∞ > 1, la serie podrı́a ser, pese a todo, convergente. Aunque el razonamiento anterior es válido para cualquier norma, ocurre que la norma-∞ es particularmente fácil de calcular. En efecto, si G tiene dimensión n, se tiene que kGk∞ = máx 1≤i≤n n X | Gij | (2.50) j=1 Para probar este resultado, partimos de la definición Pde norma (2.46). Teniendo en cuenta que la componente i del vector Ge es j Gij ej , se tiene que 16 CAPÍTULO 2. RESOLUCIÓN DE SISTEMAS LINEALES kGk∞ = máx | i X Gij ej | (2.51) j supuesto, claro está, que kek∞ = 1. Pero esta igualdad significa que todas las componentes de e están comprendidas entre 1 y -1, y que una al menos alcanza este extremo. Teniendo estoPen cuenta, es fácil ver que la manera de maximizar la componente i, que es j Gij ej , es considerar un vector ej cuyas componentes Pvalgan 1 si Gij > 0, y -1 si Gij < 0. En este caso, se tiene que P G e = ij j j | Gij |, a partir de lo cual se deduce (2.50). j Es ahora posible probar la convergencia de los métodos de Jacobi y GaussSeidel para matrices diagonalmente dominantes. • Jacobi. La matriz de la iteración es G = D−1 (N + S) (2.52) es decir Gij = Aij , i 6= j, Gii = 0 Aii (2.53) Entonces kGk∞ = máx i X | Aij | ≤1 | Aii | j (2.54) donde la última desigualdad se sigue del hecho de que A es diagonalmente dominante. • Gauss-Seidel. La matriz de iteración es G = (D + N )−1 S (2.55) Sea la ecuación y = Gx. Esta ecuación es equivalente a (D + N )y = Sx (2.56) Defı́nanse las matrices N̂ = D−1 N y Ŝ = D−1 S. Los elementos de estas matrices son: Aij Aii N̂ij =  j<i 0 j≥i Ŝij =  j>i 0 j≤i Aij Aii (2.57) 17 2.3. MÉTODOS ITERATIVOS Se tiene entonces que y = Ŝx − N̂ y (2.58) Sea entonces k la componente para la cual kyk∞ = yk . De la ecuación k-ésima del sistema de ecuaciones anterior se sigue que: yk = kyk∞ ≤ sk kxk∞ + nk kyk∞ (2.59) donde si = n i−1 X X | Aij | | Aij | ni = | Aii | | Aii | j=i+1 j=1 (2.60) Ası́ se tiene que kyk∞ ≤ sk kxk∞ 1 − nk (2.61) Por otra parte, recuérdese la definición de norma kGk∞ = máx kGxk kxk=1 (2.62) Como y = Gx, es claro que kGk∞ ≤ sk 1 − nk (2.63) Y como, al ser A diagonalmente dominante, sk +nk ≤ 1, luego sk ≤ 1−nk , y por tanto kGk∞ ≤ 1 2.3.2 (2.64) El método de la sobrerrelajación sucesiva Se ha visto en la sección anterior que los métodos de Jacobi y Gauss-Seidel son convergentes cuando se cumple que la matriz A es diagonalmente dominante. Sin embargo, las matrices que aparecen tı́picamente en los problemas de campos presentan están dominadas “por los pelos” (se cumple el = pero no el < en la definición de dominancia), por lo que cabe temer que la convergencia sea lenta. Este temor se confirma en la práctica. Existe una forma sencilla de modificar estos métodos que a menudo mejora substancialmente su velocidad de convergencia. Por ejemplo, consideremos la iteración de Jacobi 18 CAPÍTULO 2. RESOLUCIÓN DE SISTEMAS LINEALES xk = −D−1 (N + S)xk−1 + D−1 b = xk−1 + rk−1 (2.65) siendo por tanto rk−1 = −xk−1 − D−1 (N + S)xk−1 + D−1 b = −D−1 Axk−1 + D−1 b (2.66) El residuo rk−1 es la cantidad en que corregimos xk−1 para intentar obtener la solución. Los residuos son cada vez más pequeños, y tienden a cero conforme nos acercamos a la solución. El método de la sobrerrelajación sucesiva de Jacobi viene dado por la ecuación xk = xk−1 + ωrk−1 (2.67) siendo ω un número que, como se verá más adelante, está comprendido entre 0 y 2. Por supuesto, si ω = 1 se recupera el método de Jacobi. El método de Jacobi es también conocido como método de relajación de Jacobi. Como normalmente ω > 1, al método propuesto se le llama método de sobrerrelajación de Jacobi. Existe también un método de sobrerrelajación de Gauss- Seidel, dado por la misma fórmula, pero donde rk−1 es el residuo que resulta de la iteración de Gauss-Seidel rk−1 = −(D + N )−1 Axk−1 + (D + N )−1 b (2.68) Es esta segunda versión la que es utilizada con más frecuencia. La iteración sobrerrelajada de Gauss-Seidel se puede escribir, substituyendo la expresión del residuo en (2.67), como: xk = (I − ω(D + N )−1 A)xk−1 + ω(D + N )−1 b (2.69) En estos algoritmos el error ek = ωrk . Por otra parte, la matriz de error para este algoritmo es Gω = (I − ω(D + N )−1 A) Utilizando las matrices N̂ = D Gω = = = = = = −1 N y Ŝ = D −1 (2.70) S, se tiene que I − ω(D + N )−1 A I − ω(I + N̂ )−1 D−1 A I − ω(I + N̂ )−1 D−1 (I + N + S) I − ω(I + N̂ )−1 (I + N̂ + Ŝ) (I + N̂ )−1 ((I + N̂ ) − ω(I + N̂ + Ŝ) (I + N̂ )−1 ((1 − ω)(I + N̂ ) − ω Ŝ) (2.71) 19 MÉTODOS SEMIDIRECTOS Por otra parte, recuérdese que el determinante de una matriz es, por una parte, el producto de sus autovalores; y por otra el producto de los elementos de la diagonal. N y S son matrices inferior y superior, con elementos nulos en la diagonal. Por tanto det(Gω ) = Πni=1 λG i = det(I + N̂ )−1 det((1 − ω)(I + N̂ ) − ω Ŝ) = (1 − ω)n (2.72) Por consiguiente ρ(G) = máx | λG i |≥| 1 − ω | i (2.73) Y como es condición necesaria para que el algoritmo converja que el radio espectral sea menor que 1, se sigue que 0 ≤ ω ≤ 2. Se puede demostrar también que para matrices simétricas y positivas definidas, el método converge para todo ω entre 0 y 2. Por otra parte, existen algunas matrices para los que los valores óptimos de ω son conocidos. En el caso de las matrices que resultan al aplicar el método de diferencias finitas a la ecuación de Laplace, un valor de ω ligeramente por encima de 1 (1.1 por ejemplo) suele dar buenos resultados. 2.4 Métodos semidirectos: el gradiente conjugado Los métodos semidirectos permiten, en principio, como los directos, resolver el sistema lineal Ax = b en un número finito de pasos. Sin embargo, tanto en su implantación como en su comportamiento práctico son más similares a los iterativos. De entre todos ellos, el que tiene más predicamento es el método del gradiente conjugado. Tiene sin embargo una limitación: solamente se aplica a matrices A que sean simétricas definidas positivas. Sin embargo, este es muy a menudo el caso en problemas de campos. Recuérdese que A es positiva definida cuando se cumple que xT Ax > 0 ∀x (2.74) De la ecuación anterior se sigue que la función F (x) = 1 T x Ax − xT b 2 (2.75) tiene un mı́nimo para x∗ = A−1 b (2.76) 20 CAPÍTULO 2. RESOLUCIÓN DE SISTEMAS LINEALES En efecto, escrı́base x = x∗ + ∆x. Entonces F (x) = = = = = 1 ∗ (x + ∆x)T A(x∗ + ∆x) − (x∗ + ∆x)T b 2 1 ∗T 1 x Ax∗ + ∆xT Ax∗ + ∆xT A∆x − x∗T b − ∆xT b 2 2 1 1 ∗T ∗ ∗T x Ax − x b + +∆xT (Ax∗ − b) + ∆xT A∆x 2 2 1 ∗T 1 ∗ ∗T T x Ax − x b + ∆x A∆x 2 2 1 ∗ T F (x ) + ∆x A∆x ≥ F (x∗ ) (2.77) 2 Se ha empleado el hecho de que A es simétrica para pasar de la primera a la segunda ecuación, y el hecho de que es positiva definida para poder escribir la desigualdad en la última. Se ha reducido, pues, el problema de resolver un sistema lineal al de encontrar el mı́nimo de una función. Una algoritmo popular para encontrar el mı́nimo es el algoritmo del gradiente. Reuérdese que el gradiente de una función indica la dirección de máxima pendiente. Parece pues razonable escoger esta dirección para minimizar F . En el caso presente, el gradiente de F en el punto x es ∇F = Ax − b = g (2.78) xk+1 = xk − αk gk (2.79) y un posible algoritmo serı́a La cantidad αk se calcula de forma que el decremento sea tan grande como sea posible. Pero F (xk+1 ) = F (xk − αk gk ) 1 T T (xk − αk gk ) A (xk − αk gk ) − (xk − αk gk ) b = 2

1 2 T αk gk Agk − αk gkT Axk − gkT b + = 2 1 T x Axk − xT b (2.80) 2 k Como xk y gk ya estan fijos, el αk óptimo se encuentra buscando el mı́nimo de esta expresón respecto a αk e igualando a 0. Luego
gkT Axk − gkT b gkT (Axk − b) gkT gk


= = (2.81) αk = gkT Agk gkT Agk gkT Agk Ası́ pues, el algoritmo tomarı́a la forma 21 MÉTODOS SEMIDIRECTOS gk αk xk+1 = Axk − b gkT gkT
= gkT Agk = xk − αk gk (2.82) siendo necesaria, por supuesto, alguna estimación inicial x0 (por ejemplo, x0 = 0). El algoritmo anterior (el algoritmo del gradiente) converge, aunque, en general, de una manera más bien lenta. Una forma de acelerar la convergencia es modificarlo de la siguiente manera: gk βk dk αk xk+1 = Axk − b dTk−1 Agk = dTk−1 Adk−1 = gk − βk dk−1 dTk gk
= dTk Adk = xk − αk dk (2.83) La diferencia básica es que, en lugar de actualizar xk según el gradiente gk , se baja a lo largo de una dirección modificada dk , egún se observa en la última ecuación. Esta dirección se calcula a partir del gradiente y de la dirección que se utilizó en la iteración anterior: dk = gk − βk dk−1 . El número βk se ha calulado de forma que dTk−1 Adk = 0 (¡comprúebese!). El valor αk se calcula según el mismo procediento que ya se empleó en el algoritmo del gradiente. Es de importancia básica para comprender la eficacia del algoritmo el que se verifica que giT gj = 0 ∀i 6= j (2.84) Para demostrar este hecho es preciso obtener unos resultados previos: • gk = gk−1 − αk−1 Adk−1 Esta fórmula se obtiene reescribiendo la definición de gk gk = Axk − b = A [xk−1 − αk−1 dk−1 ] − b = [Axk−1 − b] − αk−1 Adk−1 = gk−1 − αk−1 Adk−1 (2.85) 22 CAPÍTULO 2. RESOLUCIÓN DE SISTEMAS LINEALES • dTk−1 gk = 0 Esta condición viene exigida por el hecho de que el punto xk se calcula como el punto de la recta (el valor de α) x = xk−1 − αdk−1 que hace mı́nimo a F . Pero en este punto mı́nimo, el gradiente de F (es decir, gk ) debe ser perpendicular a la recta, cuyo vector director es dk−1 . Con estos resultados, el algoritmo puede resumirse en cuatro ecuaciones: gk dk dTk Adk−1 gkT dk−1 = gk−1 − αk−1 Adk−1 = gk − βk dk−1 = 0 = 0 (2.86) (2.87) (2.88) (2.89) De hecho, a partir de estas fórmulas se pueden despejar los valores de αk y βk . Se empezará por demostrar que gkT gk−1 = 0. Naturalmente esto es lo mismo que gkT gk+1 = 0. Para ello, primero hay que ver que gkT dk−1 = 0 T = (gk+1 + αk Adk ) dk−1 T = gk+1 dk−1 = gkT dk−2 (2.90) donde se ha usado el hecho de que AT = A. Pero, por otra parte gkT dk−1 = 0 = gkT (gk−1 − βk−2 dk−2 ) = gkT gk−1 (2.91) Con esto, ya se ha probado que el gradiente obtenido en una iteración es perpendicular a los que se obtienen en la iteración anterior y posterior. Pero además Agk = Adk + βk Adk−1 βk 1 = (gk − gk+1 ) + (gk−1 − gk ) αk αk−1   βk βk 1 1 gk + = − gk+1 + − gk−1 αk αk αk−1 αk−1 = ak+1 gk+1 + bk gk + ck−1 gk−1 (2.92) 23 MÉTODOS SEMIDIRECTOS Es esta fórmula de 3 términos, junto con la simetrı́a de A, lo que permite asegurar la ortogonalidad de los vectores gk . Se empezará por comprobar que T gk−2 gk = 0. Para ello, se necesita probar una fórmula recursiva para el producto T T gk gk . Premultiplicando dk = gk − βk dk−1 por gk−1 : T gk−1 dk T = −βk gk−1 dk−1 T = −βk gk−1 (gk−1 − βk−1 dk−2 ) T = −βk gk−1 gk−1 (2.93) Por otra parte T gk−1 dk = T (gk + αk−1 Adk−1 ) dk = gkT dk = gkT (gk − βk dk−1 ) = gkT gk (2.94) Por lo tanto, se concluye que T gkT gk = −βk gk−1 gk−1 (2.95) T Premultiplicando ahora la fórmula de tres términos (2.92) por gk−1 , se obtiene T gk−1 Agk = − βk T 1 T g gk+1 + g gk−1 αk k−1 αk−1 k−1 (2.96) Es claro que también es válida la fórmula Agk−1 = − 1 αk−1 gk +  1 αk−1 − βk−1 αk−2  gk−1 + βk−1 gk−2 αk−1 (2.97) Premultiplicando por gkT : gkT Agk−1 = − 1 αk−1 gkT gk + βk−1 T g gk−2 αk−1 k (2.98) T Como la matriz A es simétrica, gk−1 Agk = gkT Agk−1 , y por tanto − βk T 1 βk−1 T 1 T gk−1 gk+1 + gk−1 gk−1 = − gkT gk + g gk−2 αk αk−1 αk−1 αk−1 k (2.99) Y teniendo en cuenta (2.95): − 1 T βk−1 T g gk+1 = + g gk−2 αk k−1 αk−1 k (2.100) 24 CAPÍTULO 2. RESOLUCIÓN DE SISTEMAS LINEALES Ahora bien, considérese la primera iteración del algoritmo. Se parte de x1 = 0. Entonces g1 = −b. Ahora, al ir a escoger la dirección d1 , no existe d0 de alguna iteración anterior. Por tanto, es natural hacer d1 = g1 = −b. Con esta selección se tiene que g3T g1 = g3T d1 = 0 (2.101) como se demuestra en (2.90). Pero entonces, se puede usar la expresión (2.100) para demostrar, por inducción, que T gk+2 gk = 0 ∀k (2.102) Se probará ahora que giT gj = 0 ∀i, j. Para ello supóngase que giT gj = 0 si i, j ≤ k (es una prueba por inducción). De hecho, ya sabemos que, por lo menos, se cumplirá para las primeras iteraciones: i, j = 1, 2, 3 o k = 3. Entonces, si i≤k−2 giT (Agk ) = (Agi )T gk = (ai+1 gi+1 + bi gi + ci−1 gi−1 )T gk = 0 (2.103) T T gk = 0. Pero, de la gk = giT gk = gi+1 puesto que los productos escalares gi−1 fórmula de tres términos ak+1 gk+1 = Agk − bk gk − ck−1 gk−1 Premultiplicando por giT (2.104) i < k − 1, y teniendo en cuenta (2.103) se sigue que giT gk+1 = 0 i < k − 1 (2.105) Como ya se habı́an probado los casos i = k − 1 e i = k, queda demostrado que los vectores gk son ortogonales. Ahora bien, si la dimensión de la matriz A es n, solamente puede haber n vectores ortogonales no nulos. Esto significa que, a más tardar para i = n + 1, se ha tenido que cumplir que gi = 0. Pero como g = Ax − b, esto es lo mismo que decir que el algoritmo converge en, a lo más, n + 1 iteraciones. De hecho a menudo la situación es todavı́a más favorable. En efecto, obsérvese que xk+1 = xk − αk dk = x1 − k X i=1 αi di = − k X αi di (2.106) i=1 ya que x1 = 0. Ahora bien, el algoritmo tiende a calcular los mayores di en las primeras iteraciones, dejando los más pequeños para el final. Por lo tanto, a menudo se llegan a precisiones razonables con muchas menos iteraciones que n. La restricción de que x1 = 0 no es tan seria como parece. Por ejemplo, supóngase que se tiene una solución aproximada x̃ desde la que se desea iniciar elalgoritmo. Entonces, se puede escribir la solución exacta x como: x = x̃ + ∆x (2.107) MÉTODOS SEMIDIRECTOS 25 La ecuación Ax = b es equivalente entonces a A (x̃ + ∆x) = b ⇒ A∆x = b − Ax̃ = b̃ siendo ya razonable resolver este sistema por ∆x1 = 0. (2.108) 26 CAPÍTULO 2. RESOLUCIÓN DE SISTEMAS LINEALES Capı́tulo 3 Ecuaciones de Laplace y Poisson: diferencias finitas ▽ Básico 3.1 Introducción Las ecuaciones de Laplace y Poisson aparecen con frecuencia en el estudio de situaciones estáticas, de situaciones de equilibrio. Por ejemplo, la ecuación que rige la distribución del campo electrostático en el vacı́o es ǫ0 ∆V = −ρ (3.1) donde V es el potencial eléctrico y ρ la densidad de carga. La ecuación anterior, escrita de una forma menos compacta, serı́a: ǫ0 ( ∂2 ∂2 ∂2 + 2 + 2 )V (x, y, z) = −ρ(x, y, z) 2 ∂x ∂y ∂z (3.2) Esta ecuación es la ecuación de Poisson en tres dimensiones. Frecuentemente se requiere su resolución en dos, o incluso una dimensión. Esto significa que en la ecuación anterior no se considera la coordenada z, o la y y la z. Por ejemplo, la ecuación de Poisson en 2 dimensiones es: ǫ0 ( ∂2 ∂2 + 2 )V (x, y) = −ρ(x, y) 2 ∂x ∂y (3.3) Sin embargo, se utilizará el operador laplaciano ∆ para indicar también el ∂2 ∂2 operador diferencial en dos dimensiones ∂x 2 + ∂y 2 . De hecho, la dimensión del operador denotado por ∆ será determinada por el contexto. La ecuación de Poisson no sirve solamente para calcular el campo electrostático. De hecho, es preciso resolver la misma ecuación para determinar el campo gravitacional en la teorı́a de Newton, la distribución estacionaria de temperaturas en los cuerpos, la velocidad con que fluyen los lı́quidos ideales, o las tensiones que aparecen en una barra sometida a torsión. 27 28 DIFERENCIAS FINITAS V ∆V = f ✏✏ ❜❜ ∂V ∂n Figura 3.1: Dominio de la ecuación de Poisson Un caso particular de gran interés es cuando la función ρ es nula. En tal caso se tiene la ecuación de Laplace: ∆V = 0 (3.4) En este capı́tulo se expondrá como resolver estas ecuaciones mediante el método de las diferencias finitas. El próximo abordará su resolución mediante el método de los elementos finitos. 3.2 Preliminares matemáticos Antes de proceder a explicar cualquier procedimiento numérico, es preciso que las ecuaciones a resolver hayan sido planteadas de una forma matemáticamente correcta. En el caso de la ecuación de Poisson, el problema a resolver es determinar el potencial V en el interior de un determinado dominio D. En la frontera de dicho dominio se da, o el propio potencial V o su derivada respecto a la normal de la superficie ∂V ∂n (ver figure 3.1) En el caso de problemas electrostáticos, el campo eléctrico es el gradiente del potencial (E = −∇V ), por lo que la derivada normal a la superficie no es más que la componente normal del campo (En = − ∂V ∂n ). Es esta, de hecho, una condición natural en problemas electrostáticos. Por ejemplo, si se está intentando calcular la resistencia de una determinada pieza, sabemos que no puede salir corriente por la superficie no conectada a electrodos (ver figura 3.2). Ası́, jn = 0. Pero como, por la ley de Ohm, j = σE, se tiene también que En = 0. O dicho de otra forma ∂V ∂n = 0. 29 CONSTRUCCIÓN DE LA MALLA V=0 I V=10 j=0 I Figura 3.2: Una resistencia En otro tipo de problemas (térmicos, hidraúlicos, ...) esta condición suele tener también un sentido fı́sico bastante claro. Además, y en cualquier caso, aparece también cuando se consideran ejes de simetrı́a en la configuración de estudio, como se verá más adelante. De cualquier manera, se tiene el siguiente problema matemático: ∆V = f sujeto a V =g ∂V ∂n = h en D en ∂D1 en ∂D2 (3.5) (3.6) (3.7) donde f (x, y, z) = − ǫρ0 (definida en todo D), g(x, y, z) (definida sólo en ∂D1 ) y h(x, y, z) (definida sólo en ∂D2 ) son funciones conocidas, y V es la función incógnita, a determinar en D. Las condiciones V = g en ∂D1 y ∂V ∂n = h en ∂D2 son colectivamente conocidas como condiciones de contorno. El primer tipo (V = g) es conocido como condición de Dirichlet, y el segundo como condición de von Neumann. Nótese que entre la dos han de cubrir todo el contorno, pero que su dominio de definición no se puede superponer (∂D1 ∩ ∂D2 = ∅). 3.3 Construcción de la malla y planteamiento de las ecuaciones En el método de las diferencias finitas se intenta calcular el valor de la función V en una serie finita de puntos, situados en las intersecciones de una red o malla (ver figura 3.3). Por regla general, esta es una red cuadrada (en 2 dimensiones) o cúbica (en 3 dimensiones), aunque otros tipos de red (rectangular, hexagonal, ...) son también posibles y se utilizan ocasionalmente. 30 DIFERENCIAS FINITAS Figura 3.3: Dominio mallado En el resto de esta sección se estudiarán problemas bidimensionales en redes cuadradas. La generalización a otras dimensiones y redes es sencilla, y se estudiarán algunos casos en otras secciones. La ecuación de Poisson es una ecuación diferencial. Sin embargo, de lo único que se dispone es del valor de V en una serie discreta de puntos, por lo que es preciso aproximar las derivadas por diferencias. De ahı́ el nombre del método. Por ejemplo, considérese el punto P de la figura 3.4. Como la malla tiene anchura h, el punto P se encuentra situado en las coordenadas (ih, jh). Para simplificar la notación, se escribirá: V (ih, jh) = Vi,j (3.8) El valor de V en el punto (i + 1, j) puede relacionarse con el que tiene en el punto (i, j) mediante la serie de Taylor: Vi+1,j = Vi,j + 1 ∂2V 1 ∂3V 1 ∂4V ∂V 2 3 |i,j h + | h + | h + |ξ ,j h4 (3.9) i,j i,j ∂x 2 ∂x2 6 ∂x3 24 ∂x4 1 Donde x(i) < ξ1 < x(i + 1). Análogamente, se pueden relacionar el valor de V en el punto (i − 1, j) con el de (i, j): Vi−1,j = Vi,j − ∂V 1 ∂2V 1 ∂3V 1 ∂4V 2 3 |i,j h + | h − | h + |ξ ,j h4 (3.10) i,j i,j ∂x 2 ∂x2 6 ∂x3 24 ∂x4 2 Sumando las ecuaciones (3.9) y (3.10) se obtiene: Vi+1,j + Vi−1,j = ∂2V | h2 + 2Vi,j + 2 i,j ∂x   ∂4V 1 ∂4V 4 | + | ξ ,j ξ ,j h 24 ∂x4 1 ∂x4 2 (3.11) 31 CONSTRUCCIÓN DE LA MALLA i, j + 1 P jh i − 1, j i, j i + 1, j i, j − 1 h O h ih Figura 3.4: La malla, de cerca = 2Vi,j + ∂2V |i,j h2 + O(h4 ) ∂x2 (3.12) donde O(h4 ) denota términos de grado cuarto o mayor en h. Suponiendo que estos términos son despreciables (lo que sucederá si h es lo bastante pequeño) se tiene: ∂2V Vi+1,j + Vi−1,j − 2Vi,j |i,j = 2 ∂x h2 De forma similar se obtiene que: ∂2V Vi,j+1 + Vi,j−1 − 2Vi,j |i,j = ∂y 2 h2 (3.13) (3.14) Luego ( ∂2V ∂2V Vi+1,j + Vi−1,j + Vi,j+1 + Vi,j−1 − 4Vi,j + )|i,j = 2 2 ∂x ∂y h2 (3.15) Ahora bien, el término izquierdo de esta ecuación es precisamente el operador diferencial ∆ (el laplaciano), luego el término derecho es la aproximación buscada que solo involucra el valor de V en los nodos de la red. De esta manera, la ecuación de Poisson en el punto (i, j) se puede aproximar mediante la fórmula en diferencias finitas: Vi+1,j + Vi−1,j + Vi,j+1 + Vi,j−1 − 4Vi,j = fi,j h2 (3.16) 32 DIFERENCIAS FINITAS Figura 3.5: Problema electrostático Ası́ se tiene ya una ecuación para todos los puntos internos a la malla. Los nodos que estén en el contorno, o próximos a él, requieren un tratamiento especial. La situación se simplifica notablemente si la malla es escogida de tal forma que se adapte al contorno. Esto es lo que se supondrá en esta sección, dejando el caso más general para más adelante. Por ejemplo, supóngase que se desea calcular el potencial electrostático en el espacio entre los dos conductores que se muestran en la figura 3.5. El conductor interno está a un potencial V = 100 y el externo a V = 0. En el espacio intermedio se verifica la ecuación de Laplace: ∆V = 0 (3.17) Es también claro que, debido a la simetrı́a de la configuración, basta con resolver un octavo. En concreto, se tienen las siguientes ecuaciones: 2Vb + 0 + 100 − 4Va Va + Vc + 0 + 100 − 4Vb Vb + Vd + Ve + 0 − 4Vc 2Vc + 2 ∗ 0 − 4Vd 2Vc + 2 ∗ 100 − 4Ve = = = = = 0 0 0 0 0 punto a punto b punto c punto d punto e (3.18) (3.19) (3.20) (3.21) (3.22) Nótese que los nodos de la red donde se dan condiciones de Dirichlet (los conductores) no son incógnitas, puesto que, por definición, se conocen allı́ los valores de los potenciales. Hay además dos nodos internos (b y c) y otros tres en el contorno (a, d y e). En este caso, debido a la simetrı́a de la configuración, se puede suponer que estos tres nodos son internos, puesto que sabemos como se relacionan los potenciales de los nodos que están a ambos lados de la lı́nea de simetrı́a (ası́, por 3.4. EL PRINCIPIO DEL MÁXIMO 33 ejemplo, en el caso del nodo a los valores de las tensiones a derecha e izquierda son iguales, que es por lo que ambos nodos se han denotado como b). Pero hay otra forma de considerar este problema. Notemos que sobre la lı́nea de simetrı́a se ha de verificar la anulación de la derivada normal: ∂V = 0 en linea simetria (3.23) ∂n En efecto, siendo las cargas simétricas a ambos lados de la lı́nea de simetrı́a, no puede haber campo eléctrico normal a la lı́nea. Ası́, en particular, se ha de verificar ∂V |a = 0 (3.24) ∂x puesto que en el punto a no puede haber componente del campo Ex , ya que ambas mitades “ tiran ” lo mismo. Por lo tanto, en general una lı́nea de simetrı́a implica la anulación de la derivada normal. Pero es esta una condición de contorno de von Neumann. Por lo tanto, las condiciones de contorno nulas de von Neumann ∂V ∂n = 0 pueden tratarse añadiendo nodos adicionales por fuera del dominio simétricos a los nodos internos. El caso de condiciones no nulas se verá en una sección posterior. Las ecuaciones (3.18-3.22) se pueden escribir también en forma matricial:      −4 2 0 0 0 Va −100  1 −4 1     0 0     Vb   −100   0   Vc  =  0  1 −4 1 1 (3.25)       0 0 2 −4 0   Vd   0  0 0 2 0 −4 Ve −200 O, en notación más compacta AV = b (3.26) En esta ecuación tanto la matriz A como el vector b son conocidos, por lo que puede resolverse para determinar V. Nótese que los valores concretos de las condiciones de contorno están ı́ntegramente incluidos en el término independiente b. Ası́, si se cambian el valor de los potenciales (por ejemplo, haciendo que el potencial del conductor interior valga 17, y el del exterior -313), lo único que serı́a preciso modificar serı́a el vector b, pero no la matriz A. Existen diversos métodos para resolver la ecuación lineal (3.26). Nótese que las ecuaciones que resultan del método de diferencias finitas (!y de otros muchos!) tienen matrices A con un gran número de ceros, o sea, son matrices ralas o cuasivacı́as. Su resolución se ha abordado en el capı́tulo anterior. 3.4 El principio del máximo Considérese de nuevo la ecuación (3.16), válida para puntos en el interior de la malla: △ Básico ▽ Avanzado 34 DIFERENCIAS FINITAS Vi+1,j + Vi−1,j + Vi,j+1 + Vi,j−1 − 4Vi,j = fi,j h2 (3.27) En el caso de la ecuación de Laplace, el término de la derecha fi,j es nulo. Ası́ pues, se puede escribir: Vi,j = 1 (Vi+1,j + Vi−1,j + Vi,j+1 + Vi,j−1 ) 4 (3.28) Es decir, el valor del potencial en el centro es la media de los cuatro que tiene alrededor. Esto significa que el valor central Vi,j no puede ser ni el máximo ni el mı́nimo, ya que la media se encuentra entre estos valores. Como esta conclusión es válida para todos los puntos en el interior del dominio, se sigue que el valor máximo y el mı́nimo de V se tienen que dar en puntos del contorno. Este resultado es conocido como el principio del máximo. El párrafo anterior se refiere a los valores de V que resultan de resolver las ecuaciones en diferencias. Ahora bien, la misma conclusión se puede establecer para los valores de V solución de la ecuación diferencial. En efecto, conforme va disminuyendo la anchura de la malla h la solución en diferencias va tendiendo a la solución de la ecuación diferencial (se da una prueba en la sección siguiente). Por tanto, como la solución de la ecuación diferencial es el lı́mite de la ecuación en diferencias, y para todas estas se verifica el principio del máximo, se sigue que también se cumple para la ecuación diferencial. Existe una extensión del principio del máximo que será de utilidad más adelante. Sea la ecuación de Poisson ∆V = f (3.29) con f > 0. Entonces se tiene, de la fórmula de diferencias finitas que Vi,j 1 = 4  Vi+1,j + Vi−1,j 1 + Vi,j+1 + Vi,j−1 − 2 fi,j h  (3.30) Es decir, V en el nodo central ha de ser menor que en al menos uno de los que le rodean. Por tanto, el máximo de V ha de alcanzarse en el contorno (aunque nada puede decirse del mı́nimo en este caso). 3.5 Errores y convergencia En la sección anterior se ha derivado una ecuación algebraica con la que se espera aproximar la ecuación diferencial que se desea resolver. Sin embargo, serı́a interesante contar con alguna idea del error cometido, y de como este error disminuye cuando la anchura de la malla h disminuye también. Se va a realizar este análisis para un problema particular: el problema puro de Dirichlet en dos dimensiones: 35 3.5. ERRORES Y CONVERGENCIA  ∂2 ∂2 + ∂x2 ∂y 2  V =f sujeto a V = g en D (3.31) en ∂D (3.32) siendo D el rectángulo definido por 0 < x < a, 0 < y < b. Las ideas básicas de este análisis son, sin embargo, generalizables a otros casos; y las conclusiones que se obtengan también. Para resolver las ecuaciones (3.31,3.32) por diferencias finitas se crearı́a una malla, que definiremos por la intersección de las lı́neas verticales xi = ih, i = 1 . . . m y las lı́neas horizontales yj = jh, j = 1 . . . n. Se supondrá además que mh = a y nh = b, de forma que la malla se adapta perfectamente al dominio D. Denotemos además por Dh el conjunto de nodos interiores al dominio y por ∂Dh el conjunto de nodos en la frontera del dominio. La aproximación a la ecuación de Poisson dentro del dominio es 1 (vi+1,j + vi−1,j + vi,j+1 + vi,j−1 − 4vi,j ) = fi,j (3.33) h2 En esta sección se denotará por V la solución de la ecuación diferencial y por v la solución de la ecuación en diferencias. El problema en diferencias se puede escribir por tanto como Lvi,j = fi,j sujeto a vi,j = gi,j en Dh en ∂Dh (3.34) (3.35) siendo L es operador en diferencias 1 (vi+1,j + vi−1,j + vi,j+1 + vi,j−1 − 4vi,j ) (3.36) h2 Considérese ahora el error ei,j en los nodos de la malla entre la solución de la ecuación diferencial Vi,j y la de la ecuación en diferencias vi,j : Lvi,j = Lei,j = LVi,j − Lvi,j = LVi,j − fi,j = Ti,j (3.37) (3.38) (3.39) El error de truncación Ti,j puede estimarse a partir del desarrollo de Taylor de V . En efecto, de la fórmula de Taylor (3.11) se sigue que: Ti,j 1 2 h = 24  ∂4V ∂4V ∂4V ∂4V + | + | + | |i,η2 ξ ,j ξ ,j i,η 1 2 1 ∂x4 ∂x4 ∂y 4 ∂y 4  (3.40) Denótese por D̄ la unión de D y ∂D. Entonces, defı́nase una constante M mediante la ecuación 36 DIFERENCIAS FINITAS  ∂4V ∂4V M = máx máx | 4 |, máx | 4 | ∂x ∂y D̄ D̄  (3.41) Ası́ pues, puede escribirse máx |Ti,j | ≤ Dh 1 2 h M 6 (3.42) Y por tanto máx |Lei,j | = máx |Ti,j | ≤ Dh Dh 1 2 h M 6 (3.43) Ahora lo que se necesita es algún resultado que relacione ei,j con Lei,j . En el anexo de esta sección se prueba que en el rectángulo D definido por 0 < x < a y 0 < y < b se tiene para cualquier función u definida en la malla 1 máx |u| ≤ máx |u| + (a2 + b2 ) máx |Lu| Dh ∂Dh 4 D̄h (3.44) Aplicando este teorema el error ei,j se tiene que 1 máx |ei,j | ≤ máx |ei,j | + (a2 + b2 ) máx |Lei,j | Dh ∂Dh 4 D̄h (3.45) Pero, en la frontera ∂Dh , ei,j = 0 ya que se están cumpliendo las condiciones de contorno. Luego máx |ei,j | ≤ D̄h 1 2 (a + b2 )h2 M 24 (3.46) Este resultado muestra que, si la solución tiene derivada cuarta, el error se reduce conforme la anchura de la malla se reduce. Prueba además que el error es de orden h2 , información de utilidad para varios métodos de mejora de la solución. En el caso más general, se obtienen fórmulas semejantes, con errores proporcionales a h2 . El paso más complicado suele ser el de buscar una relación entre ei,j y Lei,j . Sin embargo, existen algunas pocas situaciones peligrosas donde es imposible encontrar cotas de error como la anterior. La más importante es si hay esquinas con ángulos internos mayores de 180 grados. En estos casos, hay que desconfiar de la solución cerca de estos puntos1 . 3.5.1 Anexo: prueba del teorema Si u es cualquier función definida en el conjunto de los nodos D̄h = Dh ∪ ∂Dh de la región rectangular 0 ≤ x ≤ a, 0 ≤ y ≤ b, se verifica 1 máx |u| ≤ máx |u| + (a2 + b2 ) máx |Lu| Dh ∂Dh 4 D̄h 1 Cosa (3.47) lógica, por otra parte. Recuérdese que en una punta el campo tiende a infinito 37 3.5. ERRORES Y CONVERGENCIA siendo Lui,j = ui+1,j + ui−1,j + ui,j+1 + ui,j−1 − 4ui,j (3.48) Prueba Defı́nase la función φi,j mediante la ecuación: φi,j = Claramente
1 2 i + j 2 h2 , (i, j) ∈ D̄h 4
1 2 a + b2 ∀(i, j) ∈ D̄h 4 Además, en todos los puntos (i, j) ∈ Dh 0 ≤ φi,j ≤ Lφi,j = = 1 (i + 1)2 + j 2 + (i − 1)2 + j 2 + i2 4
+(j + 1)2 + i2 + (j − 1)2 − 4i2 − 4j 2 1 (3.49) (3.50) (3.51) (3.52) Defı́nanse ahora las funciones w+ y w− mediante w+ = u + N φ y w− = u − N φ (3.53) N = máx |Lui,j | (3.54) donde Dh Operando en las ecuaciones (3.53) con L y usando (3.52) se obtiene que: ± Lwi,j = ±Lui,j + N ≥ 0 ∀(i, j) ∈ Dh (3.55) donde la desigualdad se sigue de la definición de N . Pero, por la extensión del principio de máximo expuesta al final de la sección anterior, se sigue de esta última inecuación que w± alcanza su máximo en el contorno. Ası́ pues ± máx wi,j D̄h ± ≤ máx wi,j (3.56) = máx(±vi,j + N φi,j ) (3.57) ≤ máx |vi,j | + máx N |φi,j | (3.58) ∂Dh ∂Dh ∂Dh ≤ máx |vi,j | + ∂Dh ∂Dh
1 2 a + b2 N 4 (3.59) donde, para dar el último paso, se utiliza la expresión (3.50). Por otra parte, ± . Luego, por (3.53) como w± = u ± N φ, y N φ ≥ 0, se sigue que ±ui,j ≤ wi,j 38 DIFERENCIAS FINITAS Figura 3.6: Condiciones de von Neumann no nulas máx ±ui,j ≤ máx |vi,j | + D̄h Es decir ∂Dh
1 2 a + b2 N 4 1 máx |ui,j | ≤ máx |ui,j | + (a2 + b2 ) máx |Lui,j | Dh ∂Dh 4 D̄h △ Avanzado ▽ Básico 3.6 (3.60) (3.61) Problemas con los contornos En esta sección se estudian algunos problemas que se pueden presentar en los contornos. Especı́ficamente, se estudia el tratamiento de condiciones de von Neumann no nulas y qué hacer cuando la malla no se adapta bien al contorno. Existe, por supuesto, otros problemas distintos a estos que se tratan. Las fórmulas a aplicar en estos casos suelen ser bastante complicadas, aunque, después de leer esta sección, el estudiante interesado no debiera tener dificultades en establecerlas. 3.6.1 Condiciones de von Neumann no nulas Considérese la figura 3.6. Se desea encontrar una fórmula de diferencias finitas que aproxime la condición de Dirichlet ∂V = g(y) en x = x0 (3.62) ∂x Para obtener la ecuación en el punto O existen varias posibilidades. Una de ellas es introducir un punto “fantasma” F de forma que la ecuación en O se reduzca a la fórmula habitual 3.6. PROBLEMAS CON LOS CONTORNOS 39 1 (vF + vA + vB + vC − 4vO ) = fO (3.63) h2 Ahora se tiene una variable más: el potencial en el nudo “fantasma” F . Ası́ pues, se necesita una ecuación adicional para este nudo, que debiera estar relacionada con la condición de von Neumann, que no se ha usado todavı́a. Para ello, considérese de nuevo los desarrollos de Taylor expuestos en la sección 3.3: Vi+1,j = Vi,j + 1 ∂2V 1 ∂3V 1 ∂4V ∂V 2 3 |i,j h + | h + | h + |ξ ,j h4 (3.64) i,j i,j ∂x 2 ∂x2 6 ∂x3 24 ∂x4 1 Donde x(i) < ξ1 < x(i + 1). Análogamente, Vi−1,j = Vi,j − ∂V 1 ∂2V 1 ∂3V 1 ∂4V |i,j h + |i,j h2 − |i,j h3 + |ξ ,j h4 (3.65) 2 3 ∂x 2 ∂x 6 ∂x 24 ∂x4 2 Si las dos ecuaciones anteriores se restan se obtiene: ∂V |i,j h + O(h3 ) (3.66) ∂x siendo O(h3 ) términos de orden cúbico o superior. Despreciando estos términos, se obtiene: Vi+1,j − Vi−1,j = 2 1 ∂V |i,j = (Vi+1,j − Vi−1,j ) ∂x h (3.67) o, en nuestro caso 1 ∂V |0 = gO = (VA − VF ) ∂x h que es la ecuación adicional buscada. 3.6.2 (3.68) Tratamiento de contornos curvos Considérse una malla cuadrada pero con la frontera curva como en la figura 3.7. Se desea determinar que ecuación en diferencias finitas hay que aplicar al nodo O. Aplicando de nuevo el teorema de Taylor: VA V3 ∂2V ∂V |O + h2 θ12 2 |O + O(h3 ) ∂x ∂x 2 ∂V ∂ V = VO − h |O + h2 2 |O + O(h3 ) ∂x ∂x = VO + hθ1 (3.70) ∂V ∂x |O se sigue que:   1 2 2 2 ∂2V | = V + V − V O A 3 O ∂x2 h2 θ1 (1 + θ1 ) (1 + θ1 ) θ1 Eliminando (3.69) (3.71) 40 DIFERENCIAS FINITAS Figura 3.7: Contorno curvo Realizando cuentas similares en la dirección y, se obtiene la aproximación a la ecuación de Poisson: 2VA 2VB 2V3 2V4 + + + = −2 θ1 (1 + θ1 ) θ2 (1 + θ2 ) (1 + θ1 ) (1 + θ2 ) 3.7  1 1 + θ1 θ2  VO = h2 fO (3.72) Diferencias finitas en otros sistemas coordenados En este capı́tulo se ha venido tratando el análisis de problemas bidimensionales en coordenadas cartesianas. Hay, sin embargo, circunstancias en que es preciso tratar problemas de más dimensiones o en sistemas coordenados no cartesianos. El establecimiento de las ecuaciones tridimensionales en coordenadas cartesianas es bastante trivial: lo único que hay que hacer es repetir lo que se ha hecho añadiendo la coordenada z. En particular, con una malla cúbica se obtiene la aproximación a la ecuación de Laplace: Vi+1,j,k + Vi−1,j,k + Vi,j+1,k + Vi,j−1,k + Vi,j,k+1 + Vi,j,k+1 − 6Vi,j,k = 0 (3.73) en que de nuevo se ve que el potencial central es la media de los adyacentes. En el caso de coordenadas no cartesianas el análisis es también similar. Se parte de la ecuación diferencial, cuyas derivadas se aproximan por diferencias finitas a partir de desarrollos de Taylor en las coordenadas utilizadas. Los dos siguientes apartados realizan estas aproximaciones, de forma somera, para dos casos particulares de cierta relevancia. Otros casos se tratarı́an de forma similar. OTROS SISTEMAS COORDENADOS 41 Figura 3.8: Diferencias finitas en polares 3.7.1 Diferencias finitas en coordenadas polares En configuraciones con contornos circulares, el uso de coordenadas polares puede estar indicado. Considérese, por ejemplo, la ecuación de Laplace: 1 ∂V 1 ∂2V ∂2V + + =0 (3.74) ∂r2 r ∂r r2 ∂θ2 Defı́nanse los nodos de la malla en el plano r - θ por las intersecciones de los cı́rculos r = ih, (i = 1, 2, . . .) y las lı́neas rectas θ = jl, (j = 0, 1, 2, . . .) (ver la figura 3.8) La ecuación de Laplace puede aproximarse entonces por: 1 (Vi,j+1 − 2Vi,j + Vi,j−1 ) (Vi+1,j − 2Vi,j + Vi−1,j ) 1 (Vi+1,j − Vi−1,j ) + + =0 h2 ih 2h (ih)2 l2 (3.75) Reordenando términos     1 1 Vi,j Vi−1,j + 1 + Vi+1,j − 2 1 + 2i (il)2 1 1 + Vi,j−1 + Vi,j+1 = 0 (il)2 (il)2  1− 1 2i  (3.76) Obsérvese que también en este caso el potencial central es la media de los que hay alrededor. Por otra parte, el planteamiento de las condiciones de contorno es similar al del caso cartesiano. 3.7.2 Diferencias finitas con simetrı́a axial Existen configuraciones que presentan simetrı́a respecto a un eje. En estos casos, el uso de coordenadas cilı́ndricas está indicado. Por ejemplo, la ecuación de 42 DIFERENCIAS FINITAS z l h r Figura 3.9: Diferencias finitas con simetrı́a axial Laplace en coordenadas cilı́ndricas es: 1 ∂V ∂2V ∂2V + + =0 (3.77) ∂r2 r ∂r ∂z 2 donde se han anulado las derivadas respecto a θ, debido a la simetrı́a axial de la configuración. Considérese la malla en el plano r - z por las intersecciones de las lı́neas r = ih, (i = 0, 1, 2, . . .) y las lı́neas rectas z = jl, (j = 0, 1, 2, . . .) (ver la figura 3.9) Una aproximación en diferencias finitas a esta ecuación viene dada por: (Vi+1,j − 2Vi,j + Vi−1,j ) 1 (Vi+1,j − Vi−1,j ) (Vi,j+1 − 2Vi,j + Vi,j−1 ) + + =0 h2 ih 2h l2 (3.78) Reordenando términos       1 1 1 1 1 1 1− Vi−1,j + 2 1 + Vi+1,j − 2 + 2 Vi,j h2 2i h 2i h2 l 1 1 + 2 Vi,j−1 + 2 Vi,j+1 = 0 (3.79) l l Esta ecuación no es válida en los puntos del eje, donde i = 0. No obstante, se pueda aplicar a estos nodos la ecuación de Laplace en cartesianas, puesto que la malla alrededor de estos nodos en 3 dimensiones en nada se diferencia de una malla cartesiana (ver figura 3.10). Luego se tiene 43 OTROS SISTEMAS COORDENADOS Figura 3.10: La malla en el eje 4 Vi+1,j − h2 para los puntos del eje.  4 2 + 2 h2 l  Vi,j + 1 1 Vi,j−1 + 2 Vi,j+1 = 0 l2 l (3.80) 44 DIFERENCIAS FINITAS Capı́tulo 4 Ecuaciones de Laplace y Poisson: elementos finitos ▽ Básico 4.1 Introducción En la sección anterior se estudió una manera de aproximar ecuaciones diferenciales en derivadas parciales mediante un sistema de ecuaciones algebraicas: el método de las diferencias finitas. Existe otro método importante para efectuar la discretización de las ecuaciones en diferencias finitas: el método de los elementos finitos. El propósito de este capı́tulo es exponer sus fundamentos. No obstante, antes de comenzar, dos comentarios parecen oportunos. El primero es que estos dos métodos, aunque sean los más importantes, no son los únicos. El segundo es que ambos tienen ventajas e inconvenientes, y que cada uno está indicado en diferentes circunstancias. Más sobre esto será dicho al final del presente capı́tulo. 4.2 Forma variacional de la ecuación de Poisson La base del método de los elementos finitos se basa en la posibilidad de replantear el problema de Poisson: ∇ (ǫ∇V ) = f sujeto a V =g ǫ ∂V ∂n = h en D en ∂D1 en ∂D2 (4.1) como un problema de optimización. En efecto, considérese el siguiente problema: mı́n F V = 1 2 ZZZ 2 ǫ (∇V ) dv − ZZZ D D 45 f V dv − ZZ ∂D2 hV ds 46 ELEMENTOS FINITOS sujeto a V = g en ∂D1 (4.2) Es este un problema de optimización funcional: se trata de encontrar la función V (x, y, z) que cumpla la restricción y minimice el funcional1 F (V ). Es decir, para cada V que cumpla ser igual a g en la frontera ∂D1 , se puede calcular un número real F (ya que f y h son datos). Después se escoge la V que de el F más pequeño. El resto de esta sección se empleará en demostrar que la función buscada V verifica la ecuación de Poisson. Sea V ∗ la solución de la ecuación de Poisson (4.1). Cualquier función V en la ecuación (4.2) puede escribirse como V = V ∗ + δV (4.3) Por supuesto, tanto V como V ∗ tienen que cumplir la restricción V = V ∗ = g en ∂D1 (4.4) V ∗ + δV = V ∗ = gen ∂D1 (4.5) δV = 0 en ∂D1 (4.6) Por tanto Luego Escrı́base el funcional F (V ) como la suma de tres términos F (V ) F1 (V ) F2 (V ) F3 (V ) = F1 (V ) + F2 (V ) + F3 (V ) ZZZ 1 2 ǫ (∇V ) dv = 2 D ZZZ = f V dv D ZZ = hV ds (4.7) (4.8) (4.9) ∂D2 (4.10) Se estudian ahora cada uno de los tres términos • F1 (V ). Se tiene que F1 (V ) = = 1 Un ZZZ 1 2 ǫ (∇V ) dv 2 ZZZD 1 2 ǫ (∇(V ∗ + δV )) dv 2 D funcional es una ”función”que asigna a una función un número. (4.11) (4.12) 4.2. FORMA VARIACIONAL DE LA ECUACIÓN DE POISSON ZZZ ZZZ 1 1 2 ∗ 2 = ǫ (∇V ) dv + ǫ (∇δV ) dv + 2 2 D D ZZZ (ǫ∇V ∗ ) (∇δV ) dv D ZZZ ZZZ 1 1 2 2 = ǫ (∇V ∗ ) dv + ǫ (∇δV ) dv + 2 2 D D ZZZ ZZZ ∗ δV ∇ (ǫ∇V ∗ ) dv ∇ (δV ǫ∇V ) dv − D D ZZZ ZZZ 1 1 2 2 = ǫ (∇V ∗ ) dv + ǫ (∇δV ) dv + 2 2 D ZZZ D ZZ ∗ δV f dv (δV ǫ∇V ) ds − D ∂D ZZZ ZZZ 1 1 2 2 = ǫ (∇V ∗ ) dv + ǫ (∇δV ) dv − 2 2 D D ZZ ZZZ ∂V ∗ ds − δV ǫ δV f dv ∂n ∂D2 D ZZZ ZZZ 1 1 2 2 ǫ (∇V ∗ ) dv + ǫ (∇δV ) dv − = 2 2 D ZZ D ZZZ δV f dv δV h − 47 (4.13) (4.14) (4.15) (4.16) (4.17) D ∂D2 Para realizar este desarrollo se ha aplicado para pasar de (4.13) a (4.14) la fórmula de integración por partes2 ; para pasar de (4.14) a (4.15), el teorema de Gauss y el hecho de que V ∗ verifica la ecuación de Poisson; para pasar de (4.15) a (4.16) la condición de contorno de que δV se ha de anular en ∂D1 ; y por último, para pasar de (4.16) a (4.17), la condición de contorno que verifica V ∗ en ∂D2 . • F2 (V ). En este caso F2 (V ) = = ZZZ (4.18) f V dv ZZZD f V ∗ dv + ZZ hV ds D ZZZ f δV dv (4.19) D • Finalmente F3 (V ) F3 (V ) = (4.20) ∂D2 = ZZ ∂D2 hV ∗ ds + ZZ hδV ds (4.21) ∂D2 2 O sea, ∇ (A∇B) = ∇(A)∇(B) + A∆B. Para probar esta fórmula, basta con escribir los operadores diferenciales en, por ejemplo, cartesianas y operar. 48 ELEMENTOS FINITOS Sumando ahora F1 , F2 y F3 : F (V ) = F1 (V ) + F2 (V ) + F3 (V ) ZZZ ZZZ 1 1 2 2 = ǫ (∇V ∗ ) dv + ǫ (∇δV ) dv − 2 2 D D ZZZ ZZ δV f dv + δV h − D ∂D ZZZ ZZZ 2 f δV dv + f V ∗ dv + D ZZ ZZ D hV ∗ ds + hδV ds ∂D2 ∂D2 ZZZ 1 2 = ǫ (∇V ∗ ) dv + 2 ZZ ZZZ D hV ∗ ds + f V ∗ dv + ∂D2 D ZZZ 1 2 ǫ (∇δV ) dv 2 D ZZZ 1 2 ∗ ǫ (∇δV ) dv = F (V ) + 2 D (4.22) (4.23) (4.24) Se observa que F tiene un mı́nimo para V = V ∗ , ya que la segunda integral es positiva para cualquier δV 6= 0 puesto que ǫ > 0. Pero esto es equivalente a decir que la solución a la ecuación de Poisson (4.1) es la misma que la del problema (4.2). Al llegar aquı́, se hace necesario introducir alguna nomenclatura. Las ecuaciones (4.1) y (4.2) son, como se ha visto, esencialmente la misma. Ası́, la ecuación (4.1) es conocida como la forma fuerte, y la (4.2) como la forma débil de la ecuación. Hay también dos tipos de condiciones iniciales. Unas (en este caso, las de Dirichlet) tienen que ser introducidas como una restricción adicional en la forma débil. Estas condiciones de contorno son conocidas como condiciones esenciales. Existen otras (en el caso presente, las de von Neumann) que se introducen en la forma débil dentro del funcional F (V ) a minimizar. Estas condiciones se conocen como naturales. En general, cuando se consigue reformular un problema de campos como la minimización de determinado funcional, se dice que se ha obtenido la forma variacional de problema. En nuestro caso particular, se puede comprobar que el funcional F tiene dimensiones de energı́a. De hecho, lo único que se ha hecho es establecer el principio de trabajos virtuales para este problema particular. 4.3 Aproximación funcional La forma débil (4.2) es la base para el desarrollo del método de los elementos finitos. Sin embargo, esta expresión da la solución como el mı́nimo sobre el 49 4.3. APROXIMACIÓN FUNCIONAL conjunto de todas las funciones3 . Este es, por supuesto, un conjunto infinito. Lo que se hace entonces es aproximar la función V como una suma de funciones conocidas, preestablecidas: X ai Vi (x, y, z) (4.25) V (x, y, z) = V0 (x, y, z) + i En esta expresión, las funciones V0 y Vi son conocidas, y lo que se va a intentar es encontrar los coeficientes, los números reales ai , que mejor aproximen la función solución V ∗ . La función V0 se diferencia de las funciones Vi en lo siguiente: V0 verifica las condiciones de contorno esenciales: V0 = g en ∂D1 (4.26) mientras que las Vi verifican condiciones nulas en el mismo contorno Vi = 0 en ∂D1 (4.27) Es claro, por tanto, que V verifica las condiciones esenciales. Substituyendo la aproximación en el funcional a minimizar se obtiene: F (V ) = 1 2 ZZZ 2 ǫ (∇V ) dv − D ZZZ f V dv − D ZZ hV ds !2 ZZZ X 1 ai Vi ) dv − ǫ ∇(V0 + = 2 D i ZZZ ZZ X X ai Vi )h ds ai Vi )f dv − (V0 + (V0 + D = 1 2 ∂D2 i ZZZ 2 ǫ (∇V0 ) dv + D X ai i (4.28) ∂D2 ZZZ (4.29) i ǫ (∇V0 ∇Vi ) dv + D ZZZ 1 XX ǫ∇Vi ∇Vj dv − ai aj 2 i j D ZZZ X ZZZ ai Vi f dv − V0 f dv − D ZZ V0 h ds − ∂D2 = D i 1 XX ai aj 2 i j X i ZZZ ai ZZ Vi h ds (4.30) ∂D2 ǫ∇Vi ∇Vj dv + D 3 O más bien, las funciones derivables que sean además integrables en D, para que la forma débil tenga sentido. Aunque la caracterización rigurosa del conjunto de funciones admisibles es importante para demostrar la consistencia matemática del esquema, aquı́ se adoptará un punto de vista menos riguroso, y se supondrá que todas las funciones consideradas se pueden derivar e integrar sin problemas. 50 ELEMENTOS FINITOS X ai i ZZZ ǫ∇V0 ∇Vi dv − ZZZ Vi f dv − D D ZZ ∂D2  Vi h ds +  ZZZ  ZZZ ZZ 1 2 ǫ (∇V0 ) dv − V0 f dv − V0 h ds 2 D D ∂D2   ! X 1 X X = ai aj Ki,j  + ai bi + c 2 i i j (4.31) (4.32) donde, en la última lı́nea, se han definido las cantidades: = ZZZ ǫ∇Vi ∇Vj dv (4.33)  bi = ǫ (∇V0 ∇Vi ) dv − Vi f dv − Vi h ds (4.34) D D ∂D2   ZZZ ZZ ZZZ 1 2 V0 h ds (4.35) V0 f dv − ǫ (∇V0 ) dv − c = 2 ∂D2 D D Ki,j D ZZZ ZZZ ZZ En forma matricial, la ecuación (4.32) puede escribirse como F = 1 T a Ka + bT a + c 2 (4.36) donde es claro que a es el vector columna de componentes ai , b el de componentes bi , K la matriz cuyos elementos son Ki,j y el superı́ndice T denota la transposición. F es una función de a, ya que las cantidades K, b y c solamente dependen del problema y de las funciones V0 y Vi , que son conocidas. Por otra parte, sabemos que la solución de la ecuación de Poisson es la función que minimiza F . Por tanto, los coeficientes ai buscados serán aquellos que minimicen F , es decir, que se ha de cumplir: ∂F = 0 ∀i ∂ai (4.37) Ka + b = 0 (4.38) a+ = −K −1 b (4.39) o, en forma matricial Luego los a+ i buscados serán y la solución V + = V0 + X i a+ i Vi (4.40) 51 4.3. APROXIMACIÓN FUNCIONAL La solución aproximada V + no es la solución exacta V ∗ , ya que esta última es el mı́nimo sobre todas las funciones que cumplen las condiciones de contorno esenciales, mientras que V + es el mı́nimo sobre un subconjunto de estas: las que se pueden generar con V0 y las Vi . Sin embargo, si el conjunto de funciones Vi es lo bastante grande, la solución aproximada V + no será muy distinta de la solución V ∗ . 4.3.1 Un ejemplo Sea el problema de Poisson   d2 V dx2 = x(1 − x) x ∈ [0, 1] V (0) = 0  V (1) = 1 (4.41) Este es un problema en una sola dimensión. El dominio D es el segmento [0, 1], y la “superficie” que lo encierra los dos puntos x = 0 y x = 1, donde se nos dan condiciones de contorno esenciales. La teorı́a explicada anteriormente se aplica también a este caso, sin más que substituir las integrales de volumen por integrales de lı́nea, y las integrales de superficie por los valores de las funciones en los puntos extremos del intervalo4 . De todas formas, para el problema (4.41) es fácil encontrar una solución analı́tica integrando la ecuación diferencial. Se obtiene ası́: V (x) = − 1 4 1 3 11 x + x + x 12 6 12 (4.42) Se va ahora a obtener una solución aproximada por el método variacional. Se ha de desarrollar V como V = V0 + X ai Vi (4.43) i V0 ha de cumplir las condiciones de contorno esenciales. Una posibilidad sencilla es: V0 (x) = x (4.44) En efecto V0 (0) = 0 y V0 (1) = 1, que son las condiciones esenciales. Por otra parte, las Vi han de tener condiciones esenciales nulas. Se escogen ası́: Vi (x) = sin (iπx) (4.45) ya que Vi (0) = Vi (1) = 0. Ahora se puede proceder al cálculo de la matriz K y del vector b. Recuérdese que se obtuvo que 4 De todas formas, el lector interesado puede rehacer los desarrollos anteriores para el caso monodimensional, más sencillo. 52 ELEMENTOS FINITOS Ki,j = ZZZ ǫ∇Vi ∇Vj dv (4.46) D Luego 1 = Z = 1 δij π 2 i2 2 dVi dVj dx 0 dx dx    Z 1 d d sin(πix) sin(πjx) dx = dx dx 0 Z 1 πi sin(πix)πj sin(πjx)dx = Ki,j 0 (4.47) donde δij denota el sı́mbolo de Kronecker (δij = 1 si i = j, deltaij = 0 si i 6= j). Por otra parte, se obtuvo antes que bi = ZZZ ǫ (∇V0 ∇Vi ) dv − ZZZ Vi f dv − D D ZZ ∂D2  Vi h ds (4.48) En este ejemplo, se tiene pues bi   Z 1 d d x(1 − x) sin(πix) dx x sin(πix) dx − = dx dx 0 0 Z 1 Z 1 x(1 − x) sin(πix) dx πi cos(πix) dx − = 1  Z 1 Z 0 0 = − x(1 − x) sin(πix) dx 0 = −
2 1 − (−1)i (πi)3 (4.49) Ahora la ecuación Ka+ = −b (4.50) da X Ki,j a+ j = −bi (4.51) j O sea X1 j 2 δij π 2 i2 a+ j =
2 1 − (−1)i 3 (πi) (4.52) 53 4.3. APROXIMACIÓN FUNCIONAL 1 0.8 0.6 0.4 0.2 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -6 1 x 10 0.5 0 -0.5 -1 0 Figura 4.1: Solución analı́tica y errores de la aproximación De donde
4 1 − (−1)i 5 (πi) a+ i = (4.53) Por tanto, la solución aproximada es V = V0 + X a+ i Vi i = x+ X i
4 1 − (−1)i sin (iπx) (πi)5 (4.54) En la figura 4.1 se muestra la solución analı́tica exacta y el error entre esta y la aproximada cuando se toman 8 Vi (i = 1, . . . , 8). Se ve que el error es menor de una millonésima. En este ejemplo ha sido fácil el cálculo de los coeficientes ai . De hecho, hay en general dos dificultades para este cálculo: la evaluación de las integrales necesarias y la resolución del sistema Ka+ = b (4.55) En nuestro caso, las integrales eran sencillas y, sobre todo, la resolución de la ecuación anterior trivial, debido al factor δij en Ki,j o, dicho de otro modo, al hecho de que la matriz K era diagonal. 54 ELEMENTOS FINITOS Funciones que llevan a matrices K diagonales se llaman ortogonales. Una gran parte de los métodos “ clásicos” de resolución de ecuaciones de campo están relacionados con la determinación de estas funciones. El encontrarlas analı́ticamente es, en general, una tarea nada trivial y, a menudo, imposible. Solamente en configuraciones especiales (cilindros, esferas, ...) es ello posible5 . Por contra, el método de los elementos finitos busca matrices K no diagonales, pero si “sencillas”, utilizando funciones también sencillas. 4.4 Nomenclatura y propiedades La matriz K es conocida como matriz de rigidez, mientras que el vector b recibe el nombre de vector de fuerzas. Los nombres tienen su origen en problemas mecánicos, de elasticidad, que fueron para los que se desarrolló originariamente el método. A partir de la definición de la matriz de rigidez ZZZ ǫ∇Vi ∇Vj dv (4.56) Ki,j = D se observa que es simétrica (Ki,j = Kj,i ). Además, ha de ser positiva definida. Recuérdese que una matriz K es positiva definida si verifica que xT Kx > 0 ∀x (4.57) Supongamos entonces que K no es positiva definida. Entonces, existe algún y tal que yT Ky ≤ 0 (4.58) Hay pues dos posibilidades: que se verifique el signo < o el =. Supóngase primero que yT Ky < 0 (4.59) Considérese entonces el valor del funcional a minimizar F (αy) = α2 yT Ky + αyT b + c (4.60) Es claro que F puede hacerse tan negativo como se quiera sin más que escoger un α lo suficientemente grande. Algo similar pasa si es válido el signo igual. Entonces F (αy) = αyT b + c (4.61) de manera que F puede hacerse tan negativo como se quiera escogiendo α ≪ 0 si bT y > 0, o α ≫ 0 en caso contrario. 5 Para los cilindros se han de emplear funciones de Bessel, y para las esferas armónicos esféricos. 4.5. EL ESPACIO DE FUNCIONES ADMISIBLES 55 En cualquier caso, F puede hacerse tan negativo como se quiera. Pero este F está acotado inferiormente, puesto que no es más que una aproximación al funcional que se obtiene cuando se consideran todas las funciones posibles, que adquiere su valor mı́nimo cuando se le da como argumento la solución de la ecuación de Poisson. Esta contradicción prueba que K es definida positiva. En particular, como K es positiva definida, K −1 existe siempre. 4.5 El espacio de funciones admisibles Poco se ha dicho hasta aquı́ de cuales deben ser las caracterı́sticas que han de cumplir las funciones Vi . Por una parte, se sabe que la solución V ∗ de la forma fuerte (4.1) debe ser dos veces derivable6 , por lo que al menos debieran incluirse como posibles Vi las funciones con dos derivadas bien definidas. Pero por otra parte, en la forma débil (4.2) solamente aparecer derivadas primeras, lo que sugiere que el espacio de posibles funciones Vi debiera ampliarse hasta las funciones que tienen primera derivada, aún cuando no tengan bien definida la segunda. No existe aquı́ el propósito de definir con rigor la naturaleza del espacio de funciones Vi admisibles. Sin embargo, si se justificará la posibilidad de usar una clase particular de funciones: las funciones continuas que además tienen derivada segunda a trozos7 . La idea es dividir el dominio de definición de Vi en varias piezas. Por ejemplo, en 1 dimensión en intervalos, en 2 en un teselado de polı́gonos, en tres en un conjunto de poliedros. Ahora, dentro de cada intervalo, polı́gono o poliedro, Vi tiene que ser 2 veces derivable, aunque puede haber discontinuidades en la segunda derivada al pasar de un intervalo, polı́gono o poliedro a otro. La función Vi tiene que ser, de todas formas, continua en todo el dominio del problema. Por ejemplo, la figura 4.2 muestra una de estas funciones para el caso unidimensional:  0 x≤0 f (x) = (4.62) x x>0 Esta función es, en efecto, continua y tiene definidas las derivadas primera y △ Básico segunda en todas partes, salvo en x = 0, es decir, tiene definidas las dos primeras ▽ Avanzado derivadas a trozos. Esta función f podrı́a ser utilizada como una posible función Vi para calcular, por ejemplo, los elementos de la matriz de rigidez (4.33). Sin embargo, no tiene derivada en x = 0, con lo que el procedimiento puede parecer sospechoso8 . Sin embargo, se puede transformas f en una función f˜ que tiene todas las derivadas que se requieran en todos los puntos, y es por tanto ciertamente 6 Supuesto que ǫ es una vez derivable en inglés 8 Nótese que se demostró la equivalencia entre las formas fuerte y la débil usando ciertas operaciones diferenciales (integración por partes, etc.) que requieren que la función V tenga el número adecuado de derivadas en todos los puntos. 7 “piecewise” 56 ELEMENTOS FINITOS f x δ Figura 4.2: Una función admisible admisible, uniendo los dos tramos de f mediante un empalme suave (ver figura 4.2) de anchura δ. El punto importante es que δ puede ser tan pequeño como se desee. Pero, los elementos de K o b calculados con f o f˜ difieren en cantidades de orden δ, y por tanto son a todos los efectos equivalentes. Por ejemplo, supóngase que D es el intervalo [−1, 1], y que se desea calcular la integral que aparece en la definición de la matriz de rigidez Z I(f ) = 1 −1  df dx 2 dx = 1 2 (4.63) Si se calcula esta integral con f˜ se obtiene: I(f˜) = Z 1 Z −δ −1 = −1 = df˜ dx !2 df˜ dx !2
1 1 − δ2 + 2 dx dx + Z δ −δ Z δ −δ df˜ dx Por otra parte, en el intervalo [−δ, δ] !2 df˜ dx !2 dx dx + Z δ 1 df˜ dx !2 dx (4.64) 57 4.5. EL ESPACIO DE FUNCIONES ADMISIBLES f δ x Figura 4.3: Una función no admisible 0≤ df˜ ≤1 dx (4.65) !2 (4.66) Por tanto 0≤ Z δ −δ df˜ dx dx ≤ 2δ Luego Finalmente

1 1 1 − δ 2 ≤ I(f˜) ≤ 1 + 4δ + δ 2 2 2 1 |I(f ) − I(f˜)| ≤ 2δ + δ 2 2 (4.67) (4.68) Como δ puede ser tan pequeño como se quiera, es claro que ambas integrales son esencialmente iguales, y por lo tanto los elementos correspondientes de la matriz de rigidez. Por otra parte, el argumento es de un carácter general, válido para todos los elementos de K y b, por lo que queda demostrado que la clase considerada de funciones definidas a trozos es admisible. El argumento que se acaba de indicar no es, sin embargo, válido con funciones discontinuas, que no son admisibles. Por ejemplo, considérese la función escalón de la figura 4.3: 58 ELEMENTOS FINITOS f (x) =  0 x≤0 1 x>0 (4.69) Esta función puede aproximarse mediante la función admisible  x ≤ −δ  0 1 1 + x −δ < x≤δ f˜(x) =  2 2δ 1 x>δ (4.70) f˜ es una función admisible, ya que es continua y tiene derivada segunda a trozos. Por otra parte, es claro que aproxima a f como δ va tendiendo a cero. Sin embargo, ahora se tiene I(f˜) = Z 1 Z −δ Z δ −1 = −1 = −δ = df˜ dx !2 dx !2 Z δ df˜ dx + dx −δ !2 df˜ dx dx df˜ dx !2 dx + Z 1 2δ δ 1 df˜ dx !2 dx (4.71) Ahora la integral I, y por tanto el elemento correspondiente de la matriz de rigidez K (y, en general, todos aquellos elementos de K y b donde intervenga f ) no están definidos, tienden a infinito conforme δ tiende a cero. Por tanto, f no es una función Vi admisible. △ Avanzado ▽ Básico 4.6 Elementos finitos triangulares Ha llegado ya el momento de escoger unas funciones Vi concretas. Estas funciones deberán dar una matriz de rigidez K que sea “sencilla”. Esto significa, básicamente, que deberá tener muchos ceros, es decir, deberá ser una matriz rala. Por otra parte, las propias funciones Vi deberán ser también sencillas, para poder calcular fácilmente los valores Ki,j y bi , que resultan de hacer integrales con estas funciones. Una forma sencilla de cumplir estos requisitos en dos dimensiones, aunque muy a menudo eficaz, es la utilización de los llamados elementos triangulares. En esta sección se estudiarán estos elementos. Ante todo es necesario triangularizar el dominio donde queremos resolver la ecuación, es decir, dividirlo en triángulos tal como muestra la figura 4.4, a la derecha. Se requiere para las triangularizaciones admisibles que cada triángulo no tenga más nodos de la red que sus vértices. 4.6. ELEMENTOS FINITOS TRIANGULARES 59 Figura 4.4: Triangularizaciones admisibles e inadmisibles. Figura 4.5: Función Vi Las funciones Vi son pirámides. Cada pirámide tiene su vértice encima de un nodo de la red, mientras que su base está formada por todos los triángulos que coinciden en el vértice. La altura de la pirámide es 1. La figura 4.5 muestra una de estas Vi . Por supuesto, puede considerarse que Vi fuera del polı́gono que forma la base de la pirámide vale 0, con lo que Vi es una función continua con derivada segunda a trozos en todo el dominio D. El potencial se aproxima mediante una serie V = m X ai Vi (4.72) i=1 En este caso, las ai tienen además un significado sencillo: es el valor de V en el vértice i-ésimo. En efecto, en este vértice todas las V ’s salvo la i-ésima valen 0. Es por ello, que las ai son llamadas variables nodales. De entre todos los vértices, existirán un subconjunto que se encontrarán en 60 ELEMENTOS FINITOS el contorno donde se aplican condiciones de contorno esenciales. En estos nodos el valor de V , y por lo tanto el valor de los ai correspondientes es conocido. Dividánse entonces los nodos en dos conjuntos: uno formado por n nodos donde no existen condiciones de contorno esenciales, y por tanto las ai son incógnita; y otro formado por los m − n nodos restantes en el contorno esencial. Entonces: = V = n X i=1 n X ai Vi + m X ai Vi (4.73) i=n+1 ai Vi + V0 (4.74) i=1 V0 es conocido, puesto que se conocen las ai correspondientes. Se podrı́a, a partir de aquı́, calcular la matriz K y el vector b. Sin embargo, es más conveniente hacer estas operaciones de una forma menos directa. Considérese, entonces, cada uno de los triángulos que forman la red. Cada uno de ellos se llamará un elemento, y se denotará al i-ésimo elemento Ti . Si hay en total N elementos, se tiene que ZZZ D △ Básico ▽ Avanzado . . . dv = N ZZZ X i=1 . . . dv (4.75) Ti Dentro de cada elemento, las pirámides son planos. Esto permite calcular las integrales de forma muy eficaz. Especı́ficamente: sean (x, y)1 , (x, y)2 y (x, y)3 las coordenadas de los tres vértices del elemento Ti . Sea Vj1 la pirámide cuyo vértice esta en el punto 1. Entonces, la función Vj1 restringida al elemento Ti es Vj1 |Ti = Tj1 = 1 (x1 − x2 )(y1 − y3 ) − (x1 − x3 )(y1 − y2 ) [(y2 − y3 )x − (x2 − x3 )y+  1 ((y2 + y3 )(x2 − x3 ) − (x2 + x3 )(y2 − y3 )) 2 (4.76) (4.77) Esta ecuación se ha deducido de la siguiente manera: dado que es un plano tiene que ser de la forma Tj1 = ax + by + c. Por otra parte, Tj1 vale 1 en el nodo 1( el vértice de la pirámide ), y 0 en los nodos 2 y 3, que están en la base de la pirámide. Estas tres condiciones permiten despejar a, b y c, resultando la ecuación anterior. De manera análoga se calculan Tj2 y Tj3 . Cada función Tj es llamada función de forma. Con su ayuda es ya fácil calcular la matriz de rigidez y el vector de fuerzas. En efecto, aplicando las ecuaciones (4.56 y (4.75): 4.6. ELEMENTOS FINITOS TRIANGULARES Ki,j = ZZ 61 ǫ∇Vi ∇Vj ds D = N ZZ X k=1 ǫ∇Vi ∇Vj ds (4.78) Tk (4.79) Dentro de cada elemento solamente hace falta considerar 3 Vi ’s: las tres pirámides con vértice en cada uno de los tres vértices del elemento. En efecto, el resto de las pirámides valen 0 en ese elemento. Entonces, la matriz de rigidez elemental ZZ k ǫ∇Vi ∇Vj ds (4.80) Ki,j = Tk es, esencialmente, una matriz de 3 × 3 (o, si se quiere, una matriz n × n pero que solo tiene elementos no nulos en la submatriz correspondiente a los tres vértices del elemento). Más aún, las funciones Vi restringidas al elemento son las funciones de forma. Ası́ pues, se tiene por ejemplo que ZZ k Ki1,i2 = ǫ∇Ti1 ∇Ti2 ds (4.81) Tk k k , , Ki1,i3 Ki1,i1 etc. y de forma análoga los elementos Las funciones de forma son funciones lineales, por lo que sus derivadas son constantes. Ası́ pues, la integración es simplemente una constante por el área del triángulo. Especı́ficamente: k Ki1,i2   ∂Ti1 ∂Ti2 ∂Ti1 ∂Ti2 ǫ = ds + ∂x ∂x ∂y ∂y T ZZ k 1 × ǫ = Tk (x1 − x2 )(y1 − y3 ) − (x1 − x3 )(y1 − y2 ) 1 × (x2 − x3 )(y2 − y1 ) − (x2 − x1 )(y2 − y3 ) ((y2 − y3 )(y3 − y1 ) + (x2 − x3 )(x3 − x1 )) ds 1 × = (x1 − x2 )(y1 − y3 ) − (x1 − x3 )(y1 − y2 ) 1 × (x2 − x3 )(y2 − y1 ) − (x2 − x1 )(y2 − y3 ) ((y − y3 )(y3 − y1 ) + (x2 − x3 )(x3 − x1 )) × ZZ 2 ǫ ds ZZ Tk = 1 × (x1 − x2 )(y1 − y3 ) − (x1 − x3 )(y1 − y2 ) (4.82) 62 ELEMENTOS FINITOS 1 × (x2 − x3 )(y2 − y1 ) − (x2 − x1 )(y2 − y3 ) ((y2 − y3 )(y3 − y1 ) + (x2 − x3 )(x3 − x1 )) × 1 ǫ ((x1 − x2 )(y1 − y3 ) − (x1 − x3 )(y1 − y2 )) 2 ǫ (y2 − y3 )(y3 − y1 ) + (x2 − x3 )(x3 − x1 ) = 2 (x2 − x3 )(y2 − y1 ) − (x2 − x1 )(y2 − y3 ) (4.83) Se ha supuesto que ǫ es constante en todo el elemento, lo que suele ser habitual. Si no lo fuera, habrı́a que calcular la integral teniendo en cuenta esta dependencia. Se tiene ya calculada la matriz de rigidez elemental K k para cada elemento. Es ahora preciso juntar estas matrices elementales a fin de obtener la matriz de rigidez total K. Esto se realiza empleando la expresión (4.79): Ki,j = N ZZ X ǫ∇Vi ∇Vj ds = Tk k=1 N X k Ki,j (4.84) k=1 Recuérdese que cada matriz K k se puede interpretar como una matriz que solo tiene elementos no nulos en las posiciones que corresponden a las variables nodales del elemento k. Supóngase, entonces, que se tiene una lista de elementos tal como: Elemento 1 Elemento 2 ... Elemento N Vértices 3,1,2 Vértices 3,4,2 ... Vértices 4,5,n Al combinar las matrices K k correspondientes a estos elementos se obtendrá una matriz de rigidez de la forma:  K1 + K2 + KN 1 1 1 1+2 1 1+2 2      =    ..  . .. . 1 1+2 2 1+2 2 2 2+N N .. .. . . N N N .. . N  ...  ...   ...  ... N   ... N   ..  .. . .  ... N (4.85) donde los números indican las posiciones donde aparecen términos provenientes de cada K k . El ensamblaje del vector de fuerzas se realiza de una forma similar. Recuérdese que 63 4.6. ELEMENTOS FINITOS TRIANGULARES bi ZZ = ǫ (∇V0 ∇Vi ) ds − D N X Z Vi h dl  Z Vi h dl ∂D2 ZZ ǫ (∇V0 ∇Vi ) ds − Vi f ds − Tk Tk1 k=1 = Vi f ds − D N ZZ X = ZZ (4.86) ∂Tk2 bki  (4.87) (4.88) k=1 En esta expresión se han denotado por Tk2 los elementos que tienen algún lado en un contorno natural y por ∂Tk2 el contorno natural del elemento k. Puede haber elementos tengan ambos tipos de contornos. Conocidas las funciones V0 , f y h, pueden evaluarse todas las integrales por algún procedimiento numérico. En el caso de tener una malla lo suficientemente fina, las funciones h y f se pueden considerar constantes, que salen por tanto fuera de la integral. Esto simplifica mucho el cálculo del vector de fuerzas, puesto que queda tan solo unas integrales sobre las funciones Vi asociadas al elemento, que son sencillas de calcular. En particular ZZ 1 (4.89) Vi ds = superficie elemento 2 Tk y 1 Vi dl longitud lado 2 ∂Tk2 Z (4.90) Por otra parte, la función V0 en la integral (4.87) puede ser escrita como V0 = = (ai1 Vi1 + ai2 Vi2 )|Tk (aj1 Ti1 + ai2 Ti2 ) (4.91) (4.92) donde i1 y i2 son los dos vértices que limitan el segmento del elemento k que cae sobre el contorno esencial, y Ti1 y Ti2 las respectivas funciones de forma. Las constantes ai1 y ai2 son el valor de V en los nodos i1 e i2, que es conocido. Es claro que la única Vi que entra en la integral (4.87) es la que corresponde al tercer vértice del elemento: Vj3 . Por tanto ZZ ǫ (∇V0 ∇Vj3 ) ds = (4.93) Tk1 aj1 ZZ Tk1 ǫ (∇Tj1 ∇Tj3 ) ds + aj2 ZZ ǫ (∇Tj1 ∇Tj3 ) ds Tk1 Las dos integrales de la segunda lı́nea se evalúan como en (4.83). (4.94) 64 ELEMENTOS FINITOS Figura 4.6: Elemento triangular cuadrático Con esto, se han calculado ya la matriz K y el vector b, y queda por resolver la ecuación Ka = b 4.7 (4.95) Otros elementos Los elementos que se acaban de estudiar son elementos triangulares lineales. Se pueden generar otra serie de elementos para problemas planos substituyendo los triángulos por otros polı́gonos o las funciones lineales por polinomios de mayor grado. Una opción es, por ejemplo, seguir utilizando triángulos pero utilizar como funciones de forma, en vez de planos, funciones cuadráticas. Un elemento cuadrático general tiene la expresión a + bx + cy + dx2 + exy + f y 2 (4.96) por lo que existen 6 coeficientes a identificar (a, b, c, d, e, f ). En el caso de los elementos lineales, solamente existı́an tres coeficientes. También habı́a tres variables nodales por elemento: los valores de V en cada uno de los vértices. De hecho, cada función de forma valı́a 1 en “su” vértice y 0 en los otros dos. Esto se puede generalizar al caso cuadrático, pero ahora es preciso escoger 6 vértices. Lo adecuado es escoger los tres vértices y los tres puntos medios de los lados, tal como se muestra en la figura 4.6. Para cada nodo, es posible conseguir un polinomio cuadrático que valga 1 en el nodo considerado y 0 en los demás, ya que esto nos da 6 ecuaciones para las 6 incógnitas a, b, c, d, e, f . Las funciones de forma son estos polinomios. En el caso del elemento triangular lineal, las funciones de forma estaban coordinadas porque cada variable nodal (cada vértice) está compartida por tres 65 4.7. OTROS ELEMENTOS Figura 4.7: Elemento cuadrilátero bilineal elementos. Juntando las tres funciones de forma que confluı́an en cada vértice se obtenı́an funciones Vi en forma de pirámide. En el caso cuadrático las cosas son muy similares. Cada variable nodal está compartida por dos o tres elementos, lo que junta las funciones de forma para crear Vi que tienen forma de “pirámides” con las caras alabeadas. Al construir, de la misma manera que antes, las matrices de rigidez elementales se obtienen matrices de 6 × 6, de forma análoga a como antes se obtenı́an de 3 × 3. El ensamblamiento de estas matrices es análogo al estudiado antes. La teorı́a es ahora generalizable a polinomios de mayor grado. Se hace preciso, eso sı́, introducir un mayor número de variables nodales. Por otra parte, las matrices elementales son también de dimensión cada vez mayor. Esto significa que la matriz de rigidez total tiene cada vez un número de ceros menor. Por lo tanto, elementos de un orden elevado tienen una utilidad muy limitada, al no dar lugar a matrices lo suficientemente ralas. Otra posibilidad es cambiar la forma del elemento. En dos dimensiones, además de triángulos, tan solo se usan, normalmente, rectángulos. La posibilidad más sencilla es usar como variables nodales los cuatro vértices. Esto significa que no se pueden usar funciones de forma lineales, que tienen tres coeficientes y por tanto se quedan “cortas”, ni cuadráticas, que con seis se “pasan”. La solución es usar polinomios bilineales a + bx + cy + dxy (4.97) Estos polinomios se llaman bilineales porque son lineales en x e y separadamente. De nuevo, es posible encontrar para cada vértice un polinomio bilineal que valga 1 en el vértice considerado y 0 en el resto. Se pueden ahora generalizar la teorı́a a elementos bicuadráticos, bicúbicos, etc, en rectángulos añadiendo el número adecuado de variables nodales. Es también posible combinar elementos triángulo y rectángulo. 66 ELEMENTOS FINITOS Nótese P que se está utilizando siempre como coeficientes ai del desarrollo V = V0 + i ai Vi el valor de V en determinados nodos. Es importante hacerlo ası́, ya que esto facilita muchı́simo el asegurarse de que las funciones Vi construidas a partir de las funciones de forma de cada elemento son admisibles9 . En tres dimensiones la teorı́a es muy parecida. Los elementos triángulo pasan a ser tetraedros, pudiendo ser también lineales, cuadráticos, cúbicos, etc. Los rectángulos se transforman en paralelepı́pedos bilineales, bicuadráticos, bicúbicos, etc. Mención especial merecen los elementos usados en problemas con simetrı́a axial, conocidos como elementos axisimétricos. En estos casos es conveniente utilizar coordenadas cilı́ndricas, con el eje de simetrı́a como coordenada z. Además existen las coordenadas r y φ, siendo la derivada respecto a φ nula debido a la simetrı́a axial. En estos casos, el problema queda reducido a un problema en dos dimensiones en el plano z − r. Puede utilizarse cualquier elemento plano, pero modificando adecuadamente las integrales que definen K y b. Por ejemplo, se tiene que Ki,j ZZZ ǫ∇Vi ∇Vj dv    ZZZ ∂ ∂ ∂ ∂ 1 (rVi ) + (zVi ) (rVj ) + (zVj ) dr dz rdφ = ǫ 2 r ∂r ∂z ∂r ∂z D    ZZ3 1 ∂ ∂ ∂ ∂ = 2π ǫ (rVi ) + (zVi ) (rVj ) + (zVj ) dr dz ∂r ∂z ∂r ∂z D2 r = D La última lı́nea es la integral que se hace en el plano r−z, pudiéndose utilizar las funciones Vi vista en el caso bidimensional. 4.8 Errores y convergencia Sea V ∗ la solución de la ecuación de Poisson. Se sabe que V ∗ minimiza el funcional F que se define como F (V ) = sujeto a V = g ZZ ZZZ ZZZ 1 2 hV ds f V dv − ǫ (∇V ) dv − 2 ∂D2 D D en ∂D1 (4.98) Por tanto, al ser V ∗ el minimizador, se ha de cumplir para cualquier función V que cumpla las condiciones de contorno esenciales F (V ∗ ) ≤ F (V ) ⇒ F (V ) − F (V ∗ ) > 0 (4.99) 9 No obstante, es también frecuente, para algunos elementos de orden elevado, usar los valores de las derivadas de V en ciertos nodos como coeficientes del desarrollo. La razón última de hacerlo ası́ suele ser similar a la dada. 67 4.8. ERRORES Y CONVERGENCIA La cantidad E = F (V )−F (V ∗ ) mide ası́ la bondad de V como aproximación de la solución. En el caso de la ecuación de Poisson E admite una expresión particularmente sencilla. En efecto F (V ) − F (V ∗ ) = = = = = = ZZZ ZZZ ZZ 1 2 ǫ (∇V ) dv − f V dv − hV ds − 2 D D ∂D2 ZZZ ZZZ ZZ 1 2 ǫ (∇V ∗ ) dv − f V ∗ dv − hV ∗ ds (4.100) 2 D D ∂D2 ZZZ ZZZ ZZ 1 2 ǫ (∇(V − V ∗ )) dv − f (V − V ∗ ) dv − h(V − V ∗ ) ds + 2 D D ∂D2 ZZZ ǫ∇V ∇V ∗ dv (4.101) D ZZ ZZZ ZZZ 1 2 h(V − V ∗ ) ds + f (V − V ∗ ) dv − ǫ (∇(V − V ∗ )) dv − 2 ∂D2 D D ZZZ ZZZ ǫV ∆V ∗ dv (4.102) ∇(V ∇V ∗ ) dv − D D ZZ ZZZ ZZZ 1 2 h(V − V ∗ ) ds + f (V − V ∗ ) dv − ǫ (∇(V − V ∗ )) dv − 2 ∂D2 D D ZZZ ZZ V f dv (4.103) V ∇V ∗ ds + D ∂D ZZ ZZZ ZZZ 1 2 h(V − V ∗ ) ds + f V ∗ dv − ǫ (∇(V − V ∗ )) dv + 2 ∂D2 D D ZZ ZZ ∗ V h ds (4.104) V ∇V ds + ∂D2 ∂D1 ZZZ ZZZ 1 2 ǫ (∇(V − V ∗ )) dv + f V ∗ dv + 2 D ZZ ZZ D ∗ ∗ g∇V ds (4.105) hV ds + ∂D1 ∂D2 Ası́ pues 1 F (V ) − F (V ) = E(V ) = 2 ∗ ZZZ 2 ǫ (∇(V − V ∗ )) dv + C siendo la constante C (no depende de V ) ZZ ZZ ZZZ hV ∗ ds + f V ∗ dv + C= D ∂D2 (4.106) D g∇V ∗ ds (4.107) ∂D1 Por otra parte E(V ∗ ) = F (V ∗ ) − F (V ∗ ) = 0. Como ZZZ 1 2 ǫ (∇(V ∗ − V ∗ )) dv = 0 2 D (4.108) 68 ELEMENTOS FINITOS se sigue que C = 0. Se ve entonces que el error E no es más que como el error cuadrático, total R o medio(si se divide por D ǫ dv) a la solución. Por otra parte, la cantidad E(V ) tiene dimensiones de energı́a. Es por ello que se suele llamar la energı́a de tensión10 . Considérese ahora la aproximación a la solución obtenida con elementos finitos V + . Se puede escribir como X ai Vi (4.109) V + = V0 + i y se considera también otra combinación de las Vi X bi V i V a = V0 + (4.110) F (V + ) ≤ F (V a ) ∀V a (4.111) i El método de los elementos finitos minimiza F (v). Luego Se sigue entonces que F (V + ) − F (V ∗ ) ≤ F (V a ) − F (V ∗ ) ⇒ E(V + ) ≤ E(V a ) (4.112) Es decir, el método de los elementos finitos calcula la combinación de funciones Vi que tiene la menor diferencia cuadrática media a la solución. En este sentido, es un método “óptimo”. Ahora ya es muy fácil demostrar que la solución aproximada de elementos finitos converge a la solución del problema conforme se consideran mallas más finas. En efecto, considérese que ahora V a no es una aproximación arbitraria sino la aproximación que cumple que es igual a V ∗ en los nodos de la malla. De todas formas, se seguirá verificando que E(V + ) ≤ E(V a ), luego en un sentido medio V + aproxima mejor que V a a la solución V ∗ , aunque en los nodos lo haga peor. Conforme la malla se va haciendo más tupida, V a va aproximando cada vez mejor a V ∗ . Considérese entonces una serie de mallas cada vez más finas. Sea h la mayor anchura del mayor elemento de la malla, y V +h y V ah las soluciones de elementos finitos y la solución igual a V ∗ en los nodos para la malla de anchura h, respectivamente. Claramente, cuando h tienda a 0, V ah tenderá, de alguna forma, a V ∗ , y por tanto lı́m E(V ah ) = 0 (4.113) lı́m E(V +h ) = 0 (4.114) h→0 Y como E(V +h ) ≤ E(V ah ) h→0 Es decir, la solución aproximada V +h tiende a la solución V . 10 En inglés “strain”. 4.9. ELEMENTOS FINITOS Y DIFERENCIAS FINITAS 4.9 Elementos finitos y diferencias finitas 69 △ Avanzado El método de los elementos finitos se ha impuesto, para un gran número de ▽ Básico problemas de resolución de campos, como el método de resolución por excelencia. De hecho, en ecuaciones para las que es posible encontrar una forma variacional, suele ser el método más recomendable. Esto es debido a que permite utilizar mallas muy irregulares, en geometrı́as con diferentes materiales, sin que esto suponga esfuerzos adicionales a la hora de establecer el sistema de ecuaciones algebraico y sin preocuparse por problemas de convergencia. Por contra, en diferencias finitas, el uso de mallas irregulares y distintos materiales lleva a ecuaciones en diferencias complicadas, no siendo siempre fácil establecer errores y propiedades de convergencia para las soluciones obtenidas. Sin embargo, el método de diferencias finitas tiene a su favor su mayor generalidad, ya que parte directamente de la ecuación diferencial. En cambio, no es siempre posible encontrar una forma variacional para una ecuación diferencial dada. En términos generales, los problemas estáticos, en los que no interviene el tiempo, suelen tener forma variacional. Sin embargo, los problemas dinámicos en los que se pretende estudiar la evolución de un cierto campo inicial, no suelen tenerla. Es en estos problemas donde el método de las diferencias finitas tiene más predicamento. Por último, es preciso comentar que existen métodos que combinan ideas procedentes de ambas técnicas, ası́ como otros métodos basados en diferentes principios. 70 ELEMENTOS FINITOS Capı́tulo 5 Otros problemas estáticos ▽ Básico 5.1 Introducción El objeto de este capı́tulo es exponer la forma de resolver otros problemas estáticos del electromagnetismo. En concreto, se abordará la resolución por elementos finitos de los problemas de corrientes estacionarias y de magnestostática en medios lineales y no lineales. 5.2 Corrientes estacionarias La ecuación de continuidad para la corriente reza: ∂ρ (5.1) ∂t En situaciones estacionarias, la densidad de carga permanece constante, y por tanto ∇j = ∇j = 0 (5.2) La ley de Ohm establece que la corriente es proporcional al campo eléctrico: j = σE. Por otra parte, el campo eléctrico a menudo se puede considerar irrotacional: E = −∇V . Por tanto, se obtiene la ecuación: ∇ (σ∇V ) = 0 (5.3) Esta ecuación es la ecuación de Laplace que se vio en el capı́tulo anterior, substituyendo ǫ por σ. Mas aún, también las condiciones de contorno son del mismo tipo. En efecto, un conductor tiene algunos de sus lados a potenciales conocidos (condiciones de Dirichlet o esenciales) y otros aislados, por los que no se escapa corriente. Matemáticamente esto significa que, en este último caso, la componente normal de la densidad de corriente es nula: jn = 0 o, como j = −σ∇V , ∂V ∂n = 0 (condiciones de von Neumann o naturales). 71 72 CAPÍTULO 5. OTROS PROBLEMAS ESTÁTICOS Si se desea resolver este problema por elementos finitos, podemos escribir directamente el funcional de optimización a partir de los resultados del capı́tulo anterior: F = = ZZZ 1 2 σ (∇V ) dv 2 ZZZD 1 jE dv 2 D (5.4) Es decir, la distribución de corrientes que realmente se da es aquella que es compatible con las condiciones de contorno y que minimiza el calor producido. 5.3 Magnetostática en medios lineales Las ecuaciones de campo para campos magnetostáticos son ∇B = 0 ∇×H = j (5.5) (5.6) La anulación de la divergencia de B permite escribir que B=∇×A (5.7) En medios lineales, existe una relación entre B y H del tipo: B = µH (5.8) donde µ es una constante que depende del material, pero es independiente del campo magnético al que está sometido. Es preciso, además, especificar las condiciones de contorno. En la mayor parte de los casos solamente hay que considerar dos tipos de condiciones de contorno: • El campo B no cruza el contorno, es decir, Bn = 0, siendo Bn la componente normal al contorno. • El campo B cruza perpendicularmente al contorno. Esta condición se puede escribir Bt = 0, siendo Bt las dos componentes de B tangentes al contorno. Con las condiciones de contorno que se acaban de prescribir, las ecuaciones pueden derivarse de un principio variacional, de forma análoga a la ecuación de Poisson. En efecto, considérese el problema 5.3. MAGNETOSTÁTICA EN MEDIOS LINEALES ∇B = ∇×H = B = sujeto a Bn = 0 en ∂Dn Bt = 0 en ∂Dt 73 0 j µH (5.9) En este problema son datos la densidad de corriente j y las condiciones de contorno. La solución puede obtenerse como la solución del problema: mı́n F = sujeto a At = 0 ZZZ ZZZ 1 jA dv HB dv − 2 D D en ∂Dn (5.10) Nótese que el funcional de optimización es, de forma análoga al caso electrostático, la energı́a del campo menos la energı́a asociada a las corrientes. Por otra parte, la condición de que B sea perpendicular al contorno es una condición △ Básico natural, es decir, no requiere la introducción de restricciones en el problema de ▽ Avanzado optimización. La forma de demostrar este resultado es análoga a la que se empleó en el caso electrostático. No obstante, ahora hay que emplear el potencial vector A en vez del potencial escalar V . Sea entonces A∗ la solución de (5.10). Consideremos otras posibles soluciones A = A∗ + δA (5.11) Estas posibles soluciones verifican la condición de contorno de (5.10), es decir At = 0 en ∂Dn ⇒ δAt = 0 en ∂Dn (5.12) Nótese que esta condición garantiza que Bn = 0. En efecto, supóngase que el eje normal al contorno coincide con el eje z. Entonces Bn = Bz . Pero Bz = ∂Ax ∂Ay − ∂x ∂y (5.13) Como las componentes tangenciales son, en este caso, x e y, que son constantes (iguales a 0) en el contorno (es decir, el plano x − y), se sigue que las derivadas, y por tanto Bz , se anulan. El caso general, en función de las componentes tangenciales y normales, se sigue de este resultado. Volviendo al problema variacional, se procede a evaluar F para cualquier A admisible, es decir F (A) = F (A∗ + δA) F es la suma de dos términos (5.14) 74 CAPÍTULO 5. OTROS PROBLEMAS ESTÁTICOS ZZZ ZZZ 1 F = jA dv HB dv − 2 D D Se comienza por el primer sumando 1 HB dv 2 RRR D 1 2 1 B dv 2 D µ R R R 2 1 1 (∇ × A) 2 D µ RRR 1 2 RRR 1 2 D 1 D µ ∗ 2 RRR ∗ (5.15) = = = dv 2 (∇ × (A + δA)) dv = R R R 1 1 ∗ µ (∇ × A ) dv + D µ (∇ × A ) (∇ × δA) dv+ R R R 2 1 1 (∇ × δA) dv 2 D µ (5.16) El segundo sumando puede transformarse aplicando la fórmula del análisis vectorial, válida para dos campos a y b cualesquiera: b (∇ × a) = a (∇ × b) + ∇ (a × b) (5.17) Aplicando esta expresión al segundo sumando de (5.16): 1 D µ (∇ × A∗ ) (∇ × δA) dv i i h RRR RRR 1 1 ∗ ∗ (∇ × A ) dv − (∇ × A ) dv ∇ δA × δA∇ × µ µ D D RRR RR δA∇ × H dv − ∂D (δA × H) ds RRR D RR δA∇ × H dv − ∂Dt (δA × H) ds D RRR h = = = (5.18) El último paso se ha dado porque en ∂Dn las componentes tangenciales de δA son nulas, por lo que δA es paralelo a ds, y por tanto δA × H ds = (ds × δA) H = 0 Por otra parte, se tiene que ZZZ ZZZ ZZZ jδA dv jA∗ dv + jA dv = D D (5.19) (5.20) D Ası́ pues, se obtiene finalmente F = ZZZ ZZZ 1 1 2 jA∗ dv + (∇ × A∗ ) dv − 2 D D µ ZZZ ZZZ δA∇ × H dv − jδA dv − D ZZ D (δA × H) ds + ∂D ZZZt 1 1 2 (∇ × δA) dv 2 µ D (5.21) 75 5.4. MAGNETOSTÁTICA NO LINEAL Ahora bien, puesto que el mı́nimo de F se alcanza para A = A∗ , los términos que son lineales en δA se han de anular. En efecto, tanto δA como −δA son funciones admisibles, puesto que verifican las condiciones de contorno (5.12). Pero si los términos que multiplican δA fuesen distintos de cero, tendrı́a signos distintos para δA y −δA. Escogiendo el signo negativo, se tendrı́a un F que serı́a menor que el F en A∗ , lo que contradice que el mı́nimo de F se alcance en A∗1 . Ası́ pues se ha de verificar: ZZZ jδA dv = δA∇ × H dv − D D ZZ (δA × H) ds = ZZZ 0 0 (5.22) ∂Dt Estas ecuaciones se han de cumplir para cualquier δA que cumpla las condiciones de contorno (5.12). La primera ecuación implica que ∇×H=j (5.23) como se puede ver tomando funciones δA que sean no nulas en pequeños dominios interiores a D. Por otra parte, la segunda ecuación se puede escribir como ZZ (H × ds) δA = 0 (5.24) ∂Dt En ∂Dt , δA puede, en principio, valer cualquier cosa. Por tanto, la única forma de que se anule la integral es que H × ds = 0, o lo que es lo mismo, que H no tenga más componente que la normal al contorno. Es importante notar que el problema variacional (5.10) está escrita de forma vectorial. Ası́ pues, es posible ahora escoger las coordenadas que más convengan para resolver el problema. En casos en dos dimensiones (por ejemplo, cuando se estudian configuraciones muy largas o con simetrı́a cilı́ndrica) se ha de utilizar un potencial A que tenga tan solo una componente no nula. En este caso, la similaridad con la ecuación de Poisson es particularmente marcada. En cualquier caso, la resolución del problema por elementos finitos se plantea de la misma forma que en el capı́tulo anterior. Las tres componentes de A se escriben como una suma finitas de las funciones de forma estudiadas en el capı́tulo anterior, y se calculan la matriz análoga a la matriz de rigidez K y el vector análogo al vector de fuerzas b, resolviendo posteriormente el sistema lineal resultante. 76 CAPÍTULO 5. OTROS PROBLEMAS ESTÁTICOS 5.4 △ Avanzado Magnetostática no lineal ▽ Básico 5.4.1 Planteamiento de la forma variacional Prácticamente la totalidad de los aparatos que utilizan campos magnéticos tienen partes de hierro. El hierro es empleado porque amplifica enormemente el campo magnético. Sin embargo, desde el punto de vista analı́tico, complica grandemente el problema debido a la existencia de fenómenos de histéresis y saturación. Muy a menudo, los efectos de la histéresis pueden despreciarse. En efecto, son muchos los aparatos que utilizan hierros con ciclos muy estrechos, a fin de minimizar pérdidas. Incluso en casos en que la histéresis no sea despreciable, esta aproximación da una primera solución útil. En este caso, existe una dependencia funcional entre B y H. En efecto, puede escribirse que B = B(H) (5.25) donde B(H) es una función no lineal de H que tiende a saturarse, es decir, que verifica lı́m H→∞ ∂B(H) = µ0 ∂H (5.26) mientras que, muy a menudo, ∂B(H) ∂H |H=0 tiene un valor mucho más alto (por ejemplo, 1500µ0 ). La presencia de histéresis harı́a imposible determinar esta relación funcional, ya que a cada valor de H le corresponderı́a más de uno de B. Puede verificarse que B(0) 6= 0. Los materiales con esta caracterı́stica son imanes permanentes, que crean campo en ausencia de corriente. En lo que sigue, se supondrá que la dependencia funcional entre B y H puede describirse mediante una relación del tipo: H = H0 + B̂(H) B B (5.27) El término H0 es la fuerza coercitiva, que existe tan solo en imanes permanentes, y resulta de la imanación permanente del imán. El segundo término expresa que, en la parte del campo no atribuible a la imanación permanente (toda, salvo en imanes permanentes), los campos H y B son paralelos. Esta relación es falsa si a un imán se le aplica un campo externo muy intenso, ya que en este caso B y H acabarán siendo paralelos. Sin embargo, es aplicable a la mayor parte de los casos de interés, en los que en el interior del imán el primer término suele ser claramente mayor que el segundo. Es fundamental entonces el hecho de que las ecuaciones de campo: 1 Los términos cuadráticos en δA pueden despreciarse eligiendo δA lo bastante pequeño. 77 5.4. MAGNETOSTÁTICA NO LINEAL ∇B = ∇×H = B = sujeto a Bn = 0 en ∂Dn Bt = 0 en ∂Dt 0 j B(H) (5.28) son equivalentes al problema variacional: mı́n F Z ZZZ = B H dB 0 D ! dv − ZZZ jA; dv D sujeto a At = 0 en ∂Dn (5.29) El funcional es idéntico al del problema lineal, salvo por el término cuadrático en el campo. De hecho, la integral Z B (5.30) H dB 0 ha de interpretarse a partir de la substitución de (5.27): Z B H dB = H0 B + B (5.31) H(B) dB 0 0 Nótese que en el caso lineal H = Z Z B H(B) dB = Z B 0 0 1 µ B, H0 = 0, y por tanto 1 2 1 1 1 B dB = B = BH = BH µ 2µ 2 2 (5.32) △ Básico con lo que el funcional se reduce al del apartado anterior. ▽ Avanzado La demostración de la igualdad entre las formas fuerte (5.28) y débil (5.29) de las ecuaciones de campo se demuestra de una manera muy similar a la del caso lineal. Sea ası́ A∗ y B∗ la solución del problema (5.29), y se consideran perturbaciones δA, y por tanto δB = ∇ × δA al mismo, que cumplan las condiciones de contorno esenciales del problema (5.29). Naturalmente, se ha de cumplir que F (A∗ ) ≤ F (A∗ + δA) ∀δA (5.33) Ahora bien, se tiene que; ZZZ D Z 0 ∗ B +δB H dB ! F (A∗ + δA) − F (A∗ ) ZZZ j (A∗ + δA) dv − dv − D = 78 CAPÍTULO 5. OTROS PROBLEMAS ESTÁTICOS Z ZZZ ZZZ D Z dv − ! dv − H dB 0 D ! B∗ B ∗ +δB H dB B∗ ZZZ ZZZ jA∗ dv = ZZZ jδA dv = ZZZ jδA dv D D HδB dv − (5.34) D D donde el último paso se ha podido dar debido a la suposición de que δB es pequeño. Entonces F (A∗ + δA) − F (A∗ ) = ZZZ H∇ × δA dv − ZZZ jδA dv (5.35) D D Aplicando las fórmulas (5.17): F (A∗ + δA) − F (A∗ ) = ZZZ jδA dv = H∇ × δA dv − D D ZZZ ZZZ ZZZ jδA dv = ∇ (HδA) dv − δA∇ × H dv − D D ZZZ ZZ ZZZD δA∇ × H dv − HδA dv − jδA dv ZZZ D (5.36) D ∂D Esta expresión coincide con la que se obtenı́a en el caso lineal (ver ecuación (5.18)). A partir de este punto, el resto del análisis es idéntico: la anulación de los términos lineales en δA (que son los únicos que existen cuando δA es lo bastante pequeño), exigida por la condición de que F (A∗ ) es mı́nimo, lleva a las ecuaciones (5.28). 5.4.2 Resolución numérica de la forma variacional Volviendo de nuevo a las ecuaciones (5.29), que se vuelven a reproducir a continuación: mı́n F = ZZZ D sujeto a At = 0 en ∂Dn Z 0 B H dB ! dv − ZZZ jA dv D (5.37) es preciso establecer un proceso numérico para su resolución. De forma análoga a como se hacı́a en el capı́tulo, o en la sección, anteriores se desarrolla el potencial vector en una suma finita de funciones conocidas: 79 5.4. MAGNETOSTÁTICA NO LINEAL Ak = Ak0 + N X aki Aki (5.38) i=1 donde k = x, y, z, o k = r, θ, φ, u otras posibilidades según el sistema de coordenadas utilizado. Las funciones Aki son funciones conocidas, y las incógnitas son los coeficientes aki . Por lo tanto F (A) ≈ F (aki ) (5.39) Encontrar el mı́nimo de F es resolver el sistema de ecuaciones ∂F = 0 ∀aki ∂aki (5.40) En el caso lineal (la ecuación de Poisson o las ecuaciones de la magnetostática lineal), el funcional F es una una función cuadrática de los coeficientes aki . Por lo tanto, la ecuación anterior es una ecuación lineal, cuya resolución ya se ha abordado. En el caso no lineal, el funcional F no es una función cuadrática de los aki . En efecto, en el término ! Z ZZZ B H dB dv (5.41) 0 D B es una función lineal de aki , puesto que B depende linealmente de A (B = ∇ × A) y A depende linealmente de aki ; pero H no es función lineal de aki , ya que la función H(B) no es lineal en este caso. Por tanto, la ecuación (5.40) no es lineal. Existen varios métodos para resolver problemas no lineales. En los programas de elementos finitos el más utilizado es el método de Newton-Raphson. Este método, para resolver el sistema de ecuaciones f (x) = 0 (5.42) plantea la iteración l+1 x l =x −  ∂f ∂x −1 f (xl ) (5.43) donde xl es la estimación de la solución en la iteración l-ésima. Supóngase entonces que, para simplificar la notación2 , existe una sola componente del potencial A, como sucede, por ejemplo, en problemas cartesianos de 2 dimensiones o en problemas con simetrı́a cilı́ndrica. Entonces A = A0 k + N X ai Ai k i=1 2 Aunque el caso general se trata de la misma manera. (5.44) 80 CAPÍTULO 5. OTROS PROBLEMAS ESTÁTICOS para el caso cartesiano con el campo B con solamente componentes x e y. k es el vector unitario en la dirección z. En el caso de simetrı́a cilı́ndrica habrı́a que substituir k por uφ . Aplicando el algoritmo de Newton-Raphson a la ecuación (5.40) se obtiene que al+1 = ali − i X  ∂ 2 F −1  ∂F  | l ∂ai ∂aj ∂aj ai =ai j (5.45) Ası́ pues, es preciso evaluar la matriz ∂2F ∂ai ∂aj Kij = (5.46) y el vector ∂F (5.47) ∂aj El cálculo de estas cantidades se realiza de una manera muy similar a la vista en el capı́tulo anterior. Por ejemplo, para calcular la matriz Kij es preciso calcular dos términos: b̃j = Kij = = ∂2 ∂ai ∂aj (ZZZ Z B ∂2 ∂ai ∂aj (ZZZ Z B D D ! dv − ! ) H(al ) dB 0 H(al ) dB 0 ZZZ j A0 + D dv N X i=1 ai Ai ! ; dv ) (5.48) La derivada segunda se puede hacer en dos pasos. Se calculará primero la ∂ . Sea derivada primera ∂a i ! Z ZZZ B H(al ) dB FB = D (5.49) dv 0 Entonces, por la definición de derivada, si se considera un incremento pequeño de las cantidades ai , el incremento de FB se calculará por la fórmula: δFB = X ∂FB i ∂ai (5.50) δai Nótese que δai Ai es simplemente δA. Comparando con (5.36) ZZZ ZZ δA∇ × H dv − HδA D ∂D ZZZ δA∇ × H dv D ZZZ δai Ai k∇ × H dv D = = = X ∂FB i ∂ai δai (5.51) 5.4. MAGNETOSTÁTICA NO LINEAL 81 El paso de la primera a la segunda lı́nea se ha podido dar porque las perturbaciones consideradas δA tienen condiciones de contorno nulas (recuérdese que A0 lleva toda la información sobre las condiciones de contorno). Como los δai son arbitrarios, se sigue que ZZZ ∂FB = Ai k∇ × H dv (5.52) ∂ai D Ası́ pues, la segunda derivada se calcula como ZZZ ∂ 2 FB ∂H Ai k∇ × = dv ∂ai ∂aj ∂a j D (5.53) Pero, escribiendo H = H0 + H(B) B B (5.54) se obtiene ∂H = (5.55) ∂aj   H(B) ∂ B = (5.56) ∂aj B   ∂ H(B) ∂B H(B) B + = (5.57) ∂aj B B ∂aj   H(B) H(B) ∂B ∂B ∂ + (5.58) B ∂aj ∂B B B ∂aj P Ahora, como B = ∇ × A = ∇ × (A0 + i ai Ai )k, las derivadas de B y de B con respecto a ai son fáciles de evaluar. Conociendo B (es decir, el Bl que resultare de la estimación   l en el algoritmo de Newton-Raphson), también se H(B) ∂ puede calcular ∂B . Por lo tanto, la expresión anterior es calculable. B Supóngase además que las Ai son funciones piramidales como las del capı́tulo ∂H anterior. Es entonces fácil de comprobar que ∂a es distinto de cero solamente j en los elementos que se juntan en el vértice j. Por tanto, la matriz Kij , que se calcula mediante la fórmula (5.53), solamente tendrá elementos no nulos si i y j son vértices vecinos, o el mismo vértice. Es decir, será tan cuasivacı́a como la de la ecuación de Poisson. Esto es muy importante, puesto que en el algoritmo de Newton-Raphson hace falta invertir esta matriz. El cálculo del vector b̃, más sencillo, se deja al lector interesado. 82 CAPÍTULO 5. OTROS PROBLEMAS ESTÁTICOS Capı́tulo 6 La ecuación de difusión ▽ Básico 6.1 Introducción La ecuación de difusión es una de las más sencillas ecuaciones de sistemas que evolucionan en el tiempo. Esta ecuación es: ∂u ∂t sujeto a = u∈ u(x, 0) = f (x) u(x, t) = g(x, t) ∂u(x,t) ∂n ∇(κ∇u) D × [0, T ] x ∈ ∂D1 = h(x, t) x ∈ ∂D2 (6.1) Es decir, la ecuación de difusión iguala la laplaciana de una determinada función en un cierto dominio a la derivada parcial de dicha función respecto al tiempo. Como en el caso de la ecuación de Poisson, existen también condiciones de contorno de Dirichlet y Neumann, que en este caso pueden también depender del tiempo. Además, se da el valor de la función para el instante inicial t = 0. El problema es calcular la función u(x, t) para todos los instantes hasta un determinado tiempo máximo T . Se puede apreciar que esta es una ecuación aparentemente muy similar a la de Poisson. Sin embargo, sus propiedades y comportamiento son completamente diferentes, y los métodos de resolución numérica, aunque aparentemente similares, tienen también propiedades muy diferentes a los estudiados hasta ahora. En este capı́tulo se estudiará con especial atención la ecuación de difusión con una sola dimensión espacial. Esta ecuación ya presenta la mayor parte de los nuevos problemas y fenómenos. Más adelante, se estudiará su generalización a varias dimensiones. 83 84 CAPÍTULO 6. LA ECUACIÓN DE DIFUSIÓN 6.2 La ecuación de difusión en 1 dimensión En 1 dimensión, la ecuación de difusión se reduce a   ∂u ∂ ∂u κ = ∂t ∂x ∂x (6.2) más las condiciones de contorno. El dominio D es ahora un intervalo, que se puede tomar de la forma D = [0, L]. La constante κ puede en principio depender de la posición. En esta sección se considerará el caso en que κ sea constante. Por tanto, la ecuación anterior se reduce a ∂2u ∂u =κ 2 (6.3) ∂t ∂x Esta ecuación puede ser simplificada aún más definiendo las nuevas variables x̃ = t̃ = x L κt L2 (6.4) (6.5) Si se hace esta substitución se obtiene que (6.3) se transforma en: ∂u ∂2u = x ∈ [0, 1] (6.6) ∂ x̃2 ∂ t̃ Por lo tanto, se estudiará, sin pérdida de generalidad, la ecuación (6.6) (aunque se omitirá la tilde sobre las variable). El procedimiento que se ha aplicado es conocido como normalización, y es conveniente llevarlo a cabo, aún en los casos en que no se simplifique la ecuación, para que el ordenador maneje números en torno a la unidad. Considérense ahora las condiciones de contorno. Tras la normalización, el contorno de D son los puntos 0 y 1. Por tanto, las condiciones de contorno serán el valor u(0) o el valor ∂u ∂x |x=0 para el punto x = 0, y similarmente para x = 1. Para concretar, supóngase que se impone una condición de Dirichlet en x = 0 y una de Neumann en x = 1. Entonces, el problema es ∂u ∂t sujeto a = u(x, t) ∈ u(0, t) ∂u(x,t) ∂x |x=1 u(x, 0) ∂2u ∂x2 [0, 1] × [0, T ] = g(t) = h(t) = f (x) (6.7) Para resolver numéricamente el problema, el primer paso es discretizar el espacio. Ası́ pues, en vez de considerar los valores de u para todos los puntos de 85 6.2. LA ECUACIÓN DE DIFUSIÓN EN 1 DIMENSIÓN 1 2 3 4 N-1 N Figura 6.1: Malla unidimensional. x, se va a considerar en los puntos de una malla de anchura h, tal como muestra la figura 6.1. En los puntos internos de la malla, la laplaciana se puede aproximar por diferencias finitas: ui+1 + ui−1 − 2ui ∂2u ∂ui = |x=xi ≈ ∂t ∂x2 h2 (6.8) Cada ui es, naturalmente, una función del tiempo. La condición de contorno de Dirichlet es u0 (t) = g(t) (6.9) La condición de contorno de Neumann se puede aproximar introduciendo un nodo auxiliar imaginario en x = 1 + h1 . Este será el nodo N + 1. Entonces, en el último nodo N se aproxima la laplaciana mediante uN +1 + uN −1 − 2uN ∂2u ∂uN |x=1 ≈ = 2 ∂t ∂x h2 (6.10) Por otra parte, se tiene que la condición de Neumann ∂u uN +1 − uN −1 |x=1 = h(t) ≈ ∂x 2h (6.11) uN +1 = uN −1 + 2hh(t) (6.12) Luego Substituyendo esta expresión en (6.10) ∂2u 2uN −1 − 2uN 2h(t) ∂uN = |x=1 ≈ + ∂t ∂x2 h2 h Las ecuaciones (6.8,6.9,6.13) pueden combinarse en el sistema (6.13) 86 CAPÍTULO 6. LA ECUACIÓN DE DIFUSIÓN  u1 u2 u3 .. .   d    dt    uN −1 uN        1    = 2  h      −2 1 0 1 −2 1 0 1 −2 .. .. .. . . . 0 0 0 0 0 0 ... ... ... .. . 0 0 0 .. . 0 0 0 .. .  u1 u2 u3 .. .        . . . −2 1   uN −1 . . . 2 −2 uN           +       g(t) h2 0 0 .. . 0 2h(t) h          (6.14) En forma compacta d u(t) = Au(t) + b(t) dt 6.3 (6.15) Resolución de las ecuaciones En general, la ecuación de difusión lleva a ecuaciones diferenciales lineales de la forma anterior, donde la matriz A es la discretización de la ecuación de Laplace en el dominio considerado y el vector b contiene información sobre las condiciones de contorno. La solución se puede obtener de la expresión: u(t2 ) = exp [A(t2 − t1 )] u(t1 ) + Z t2 exp [A(t2 − s)] b(s) ds (6.16) t1 donde el exponencial de una matriz se define por la serie exp [Aτ ] = I + Aτ + 1 1 2 2 A τ + A3 τ 3 + . . . 2! 3! (6.17) Se puede comprobar que (6.16) es la solución mediante su substitución directa en (6.15), substituyendo las exponenciales por sus series definitorias. En principio, la ecuación (6.16) permite la resolución directa del problema, sin más que hacer t1 = 0 y substituir t2 por el tiempo (o tiempos) para el cual se desea la solución. Sin embargo, el calculo del exponencial de una matriz grande (que surgen si el mallado del espacio es estrecho) es una tarea numéricamente compleja. La solución es mallar también el tiempo, de forma que se evalúa la solución en los instantes tj = j ∗ l, siendo l el intervalo temporal. Si l es lo bastante pequeño, puede sustituirse la exponencial por una aproximación de cálculo más sencillo. Esto permite obtener una fórmula sencilla para pasar del instante tj al tj+1 . Algunas posibilidades son: • Método explı́cito clásico exp [Al] ≈ I + Al (6.18) 87 6.3. RESOLUCIÓN DE LAS ECUACIONES • Método implı́cito clásico exp [Al] ≈ (I − Al)−1 (6.19) • Método de Crank-Nicholson exp [Al] ≈ I + 12 Al I − 12 Al (6.20) Se analizan ahora estas posibilidades. 6.3.1 Método explı́cito clásico La aproximación empleada (6.18) resulta de cortar la serie definitoria de la exponencial (6.17) a partir del segundo término. Por tanto, aproximará bien la exponencial si l es lo bastante pequeño. Aplicando la fórmula (6.16) entre los instantes tj y tj+1 , y utilizando la aproximación (6.18) uj+1 = exp [A(tj+1 − tj )] uj + Z tj+1 exp [A(tj+1 − s)] b(s) ds tj = exp [Al] uj + Z tj+1 exp [A(tj+1 − s)] b(s) ds tj l bj+1 + bj [I + As] ds 2  0 2 l bj+1 + bj ≈ [I + Al] uj + l + A 2 2 ≈ [I + Al] uj + cj ≈ [I + Al] uj + Z (6.21) La gran ventaja de este método es que permite calcular la solución uj+1 a partir de la solución uj sin hacer nada más que unas multiplicaciones. En este sentido, podrı́a parecer que presenta una carga computacional baja. Sin embargo, se verá que esto es más apariencia que realidad, ya que se requiere un l muy pequeño para que el método funcione bien. Antes de estudiarlo, es conveniente mirar con más atención la estructura de la ecuación (6.21). Por ejemplo, para el caso monodimensional, la actualización del valor de u en la posición espacial i interior al dominio es 1 (ui+1,j − 2ui,j + ui−1,j ) l (6.22) h2 No aparece ningún término proveniente de cj porque este vector tiene componentes nulas para los nodos interiores. En efecto, b solamente tiene componentes no nulas en los nodos de frontera, y c, que resulta de una composición de b y Ab, solamente puede tenerlos en los nodos de frontera y en sus vecinos. La ecuación anterior se puede reescribir como ui,j+1 = ui,j + 88 CAPÍTULO 6. LA ECUACIÓN DE DIFUSIÓN 1 ui,j+1 − ui,j = 2 (ui+1,j − 2ui,j + ui−1,j ) l h (6.23) 2 ∂ u Comparando con la ecuación de difusión ∂u ∂t = ∂x2 , se comprueba que se ha discretizado la ecuación mediante las aproximaciones ∂u |i,j ∂t ∂2u |i,j ∂x2 ui,j+1 − ui,j l 1 (ui+1,j − 2ui,j + ui−1,j ) h2 = = (6.24) (6.25) La aproximación de las derivadas respecto al espacio es la misma que se utilizó al estudiar la ecuación de Poisson, y se puede justificar de la misma forma. En cambio, la aproximación de la derivada respecto al tiempo es una diferencia retrasada, ya que estima el valor de la derivada a partir del valor de la función en un instante posterior. 6.3.2 Método implı́cito clásico La aproximación de la exponencial utilizada en este método puede desarrollarse en serie (I − Al)−1 = I + Al + A2 l2 + A3 l3 + . . . (6.26) en el supuesto de que la serie del lado derecho converja. La fórmula puede demostrarse, en este caso, sin más que multiplicar la serie por I − Al. Se puede comprobar que la convergencia se produce si l es lo bastante pequeño1 . Comparando esta expresión con la serie que define la exponencial se puede apreciar que en este caso se comete un error principal por exceso de 21 A2 l2 , comparado con el error por defecto de la misma cantidad que se tenı́a en el método explı́cito. Utilizando esta aproximación, se obtiene que uj+1 = exp [A(tj+1 − tj )] uj + Z tj+1 exp [A(tj+1 − s)] b(s) ds tj = exp [Al] uj + Z tj+1 exp [A(tj+1 − s)] b(s) ds tj −1 ≈ [I − Al] uj + Z l −1 ds [I + As] ds bj+1 + bj 2 0 −1 ≈ [I − Al] uj + Z 0 1 Especı́ficamente, bj+1 + bj 2 [I − As] l debe ser lo bastante pequeño como para que todos los autovalores de Al tengan un módulo inferior a 1. 89 6.3. RESOLUCIÓN DE LAS ECUACIONES −1 ≈ [I − Al] −1 ≈ [I − Al]   l2 bj+1 + bj uj + Il + A 2 2 uj + cj (6.27) El método parece más trabajoso que el método implı́cito, puesto que hay que se observa la inversa de una matriz. La manera eficiente de tratar este problema es reescribir la ecuación como [I − Al] uj+1 = uj + [I − Al] cj (6.28) y resolver el sistema lineal para calcular uj+1 a partir de los valores conocidos de uj y cj . Como la matriz que aparece es siempre la misma en todos los pasos del problema, de j = 1 hasta el último, lo más eficaz es calcular al principio del algoritmo la factorización LU de I − Al, y utilizarla en adelante. Esta factorización inicial es la mayor diferencia, en términos de coste computacional, entre el método explı́cito y el implı́cito. Sin embargo, a menudo queda compensada por el hecho de que en el método implı́cito se puede utilizar pasos l mucho mayores que en el explı́cito, por razones que se verán más adelante. La ecuación (6.28) se puede escribir, para nodos del interior, como ui,j+1 − 1 (ui+1,j+1 − 2ui,j+1 + ui−1,j+1 ) l = ui,j h2 (6.29) Reordenando términos 1 ui,j+1 − ui,j = 2 (ui+1,j+1 − 2ui,j+1 + ui−1,j+1 ) l h (6.30) Comparando con la ecuación diferencial, se comprueba que se han realizado las aproximaciones ∂u |i,j+1 ∂t ∂2u |i,j+1 ∂x2 = = ui,j+1 − ui,j l 1 (ui+1,j+1 − 2ui,j+1 + ui−1,j+1 ) h2 (6.31) (6.32) Es decir, se están utilizando derivadas temporales adelantadas, en comparación con las retrasadas utilizadas en el método explı́cito. 6.3.3 Método de Crank-Nicholson En este caso, la exponencial queda aproximada mediante la matriz  1 I − Al 2 −1  1 I + Al 2  Desarrollando en serie la matriz inversa, se obtiene (6.33) 90 CAPÍTULO 6. LA ECUACIÓN DE DIFUSIÓN −1   1 1 I − Al I + Al = 2 2    1 1 1 1 I + Al = I + Al + A2 l2 + A3 l3 + . . . 2 4 8 2 1 1 I + Al + A2 l2 + A3 l3 + . . . 2 4  (6.34) Comparando con la serie que define la exponencial, se observa que el error 1 A3 l3 . Es un término de orden cúbico, en vez de cuadrático principal es de 12 como en los casos anteriores, por lo que, en principio, el método debiera ser más exacto. Procediendo como en los casos anteriores, se encuentra que uj+1 = exp [A(tj+1 − tj )] uj + Z tj+1 exp [A(tj+1 − s)] b(s) ds tj = exp [Al] uj + Z tj+1 exp [A(tj+1 − s)] b(s) ds tj  −1   1 1 ≈ I − Al I + Al uj + 2 2 −1   Z l 1 1 bj+1 + bj I − As I + As ds 2 2 2 0  −1   1 1 ≈ I − Al I + Al uj + 2 2  Z l bj+1 + bj 1 I + As + As2 ds 2 2 0  −1     1 1 1 2 1 3 bj+1 + bj ≈ I − Al I + Al uj + l + Al + Al 2 2 2 6 2  −1   1 1 ≈ I − Al I + Al uj + cj (6.35) 2 2 Como en el anterior método implı́cito, aparece la inversa de una matriz. Este problema se trata de forma semejante, calculado la factorización LU de I − 12 Al al principio del procedimiento. La ecuación anterior se puede escribir también como     1 1 I − Al uj+1 = I + Al uj + cj 2 2 Para puntos del interior, esto es lo mismo que (6.36) 91 6.4. ESTABILIDAD DE LA SOLUCIÓN ui,j+1 − 1 (ui+1,j+1 − 2ui,j+1 + ui−1,j+1 ) l h2 1 ui,j + 2 (ui+1,j − 2ui,j + ui−1,j ) l h = (6.37) Reordenando términos ui,j+1 − ui,j = l   1 1 1 (ui+1,j+1 − 2ui,j+1 + ui−1,j+1 ) + 2 (ui+1,j − 2ui,j + ui−1,j ) 2 2 h h (6.38) Se puede considerar que se han hecho las aproximaciones ∂u 1 | ∂t i,j+ 2 ∂2u 1 | ∂x2 i,j+ 2 = = ui,j+1 − ui,j  l 1 1 (ui+1,j+1 − 2ui,j+1 + ui−1,j+1 ) + 2 h2  1 (ui+1,j − 2ui,j + ui−1,j ) h2 (6.39) (6.40) Es decir, se han empleado diferencias centradas para aproximar las derivadas, tanto en el tiempo como en el espacio. 6.4 Estabilidad de la solución Todos los métodos expuestos se pueden escribir como uj+1 = Guj + cj (6.41) La matriz G depende del método escogido, siendo I + Al para el método explı́cito, (I − Al)−1 para el implı́cito clásico y (I − 21 Al)−1 (I + Al) para el Crank-Nicholson. Esta ecuación se resuelve empezando desde j = 0 y calculando para valores sucesivos de j. Esta ecuación será inútil si no es robusta frente a los errores. En efecto, es inevitable que se introduzcan errores al resolver la ecuación, aunque solamente sea debido a errores de redondeo. Los errores producidos en el paso j causan errores en el (se propagan al) paso j + 1, y de allı́ a los pasos sucesivos. Si estos errores propagados no disminuyen, sino que crecen conforme avanza el algoritmo, la solución que se obtenga será de escasa utilidad. Un algoritmo que tenga la propiedad de hacer que el efecto de los errores disminuya al avanzar la solución se llama estable. Supongamos entonces que en el instante j la solución obtenida ũj tiene un error ej . O sea, 92 CAPÍTULO 6. LA ECUACIÓN DE DIFUSIÓN ũj = uj + ej (6.42) Si se aplica el algoritmo (6.41) se obtiene, suponiendo que no haya errores adicionales que ũj+1 = uj+1 + ej+1 = Gũj = Guj + Gej + cj (6.43) Comparando con la ecuación original (6.41) se obtiene la ecuación del error ej+1 = Gej (6.44) Es decir, la ecuación del error es la misma que la de la solución, salvo que desaparece el término independiente. En general se tiene que ej+k = Gk ej (6.45) Ası́ pues, si el error se ha de anular cuando k → ∞, se debe verificar que lı́m Gk → 0 k→∞ (6.46) Esta ecuación se cumple si todos los autovalores de G tienen módulo menor que 1. Por otra parte, los autovalores de G son función de los autovalores de A, ya que G es función de A. Los autovalores de A son las soluciones de la ecuación Avi = λA i vi (6.47) Para el método explı́cito clásico, G = I + Al
Gvi = (I + Al) vi = 1 + λA i vi (6.48) λEC = 1 + λA i i l (6.49) Ası́ pues, los autovalores de G para este caso son Análogamente, se encuentra para el método implı́cito clásico λIC = i 1 1 − λA i l (6.50) 1 + 21 λA i l 1 A 1 − 2 λi l (6.51) Para el de Crank-Nicholson λCN = i Se puede demostrar que los autovalores de A están, para el caso monodimensional, entre 0 y − h42 . Ası́ pues, la condición de que el módulo de los autovalores de G sea menor que 1 es 6.4. ESTABILIDAD DE LA SOLUCIÓN • Método explı́cito clásico: l h2 < 93 1 2 • Método implı́cito clásico: vale cualquier l y h. Se dice que el método es incondicionalmente estable. • Método de Crank-Nicholson: es la misma situación que antes. El método es incondicionalmente estable. Se observa que en el método explı́cito el uso de mallas espaciales pequeñas (h pequeñas) obliga al uso de pasos temporales muy cortos (h2 extremadamente △ Básico pequeñas). Es por esto por lo que este método no es muy usado, a pesar de su sencillez. ▽ Avanzado La demostración anterior, aunque muy general, no es quizá muy transparente. Para analizar con más precisión la situación puede utilizarse el llamado método de Neumann. Este método parte de imponer un error de estructura particular. En concreto ei,0 (k) = sin (πki) (6.52) fi,0 (k) = cos (πki) (6.53) o La constante k es un número real. Si todos los senos y cosenos definidos de esta manera dan lugar a un error que decrece, también decrecerá un error inicial arbitrario. En efecto, un error arbitrario se puede escribir como suma de senos y cosenos. Para analizar la estabilidad del método explı́cito se escribe la ecuación asociada al i-ésimo nodo: ei,j+1 (k) = ei,j (k) + 1 (ei+1,j (k) − 2ei,j (k) + ei−1,j (k)) l h2 (6.54) 1 (fi+1,j (k) − 2fi,j (k) + fi−1,j (k)) l h2 (6.55) Análogamente se tiene fi,j+1 (k) = fi,j (k) + El sı́mbolo  denota la unidad imaginaria. Sea entonces gi,j (k) = fi,j (k) + ei,j (k) (6.56) Entonces, gi,j+1 (k) = gi,j (k) + Por otra parte 1 (gi+1,j (k) − 2gi,j (k) + gi−1,j (k)) l h2 (6.57) 94 CAPÍTULO 6. LA ECUACIÓN DE DIFUSIÓN gi,0 (k) = fi,0 (k) + ei,0 (k) = cos (πki) +  sin (πki) = exp (πki) Por otra parte, aplicando (6.57) gi,1 (k) 1 (gi+1,0 (k) − 2gi,0 (k) + gi−1,0 (k)) l h2 1 exp (πki) + 2 (exp (πk(i + 1)) − 2 exp (πki) + exp (πk(i − 1))) l   h l 2l l exp (πki) 1 + 2 exp (πk) − 2 + 2 exp (−πk) h h h   2l 2l exp (πki) 1 − 2 + 2 cos (πk) h h exp (πki) G(k) G(k)gi,0 (6.58) = gi,0 (k) + = = = = = Ası́ pues, los senos y cosenos con un k dados resultan multiplicados por el factor (real) G(k). Si los errores van a disminuir, se requiere que | G(k) |< 1. El k más desfavorable es aquel para el que cos (πk) = −1. En este caso 4l l 1 > −1 → 2 < (6.59) h2 h 2 que coincide con la condición deducida anteriormente. Incidentalmente, el cálculo demuestra que los errores que crecen más rápido son los senos y cosenos con k’s para los cuales cos (πk) = −1. Estas funciones cambian signo de un nodo i al vecino (i + 1 o i − 1). La figura 6.2 muestra la evolución del algoritmo para pasos temporales por encima y por debajo del valor crı́tico. En concreto, se simuló el caso monodimensional con 10 nodos y h = 1. El contorno de la izquierda era de Dirichlet, con potencial nulo; y el de derecha de Neumann, con derivada nula. La simulación estable, mostrada a la izquierda, se hizo con un paso temporal l = 0,4. La de la derecha, inestable, con l = 0,6. Se observa, en el caso inestable, como crecen los errores del tipo indicado en el párrafo anterior. Los algoritmos implı́cito clásico y de Crank-Nicholson se analizan de forma similar. G(k) = 1 − 6.5 Convergencia de los algoritmos En las secciones precedentes se ha analizado de que manera es posible discretizar la ecuación diferencial de forma que sea resoluble en un ordenador. Es de esperar 95 6.5. CONVERGENCIA DE LOS ALGORITMOS 14 14 12 12 10 10 8 8 6 6 4 4 2 2 0 0 5 10 0 0 5 10 Figura 6.2: Algoritmo estable e inestable que conforme se disminuya la anchura de la malla h y el paso temporal de integración l se vayan obteniendo soluciones que aproximen cada vez mejor la de la ecuación diferencial. Además de los algoritmos expuestos arribas, existen muchos más que presentan sus propias ventajas e inconvenientes, y que muy a menudo se presentan como una ecuación en diferencias sobre la malla. Un requerimiento para que la solución aproximada converja a la exacta es que el algoritmo sea estable. Es natural, además, exigir que sea consistente, es decir, que aproxime la ecuación diferencial correcta cuando se reduce la malla. Un ejemplo aclarará este concepto. Para aproximar la ecuación de difusión ∂u ∂ 2 u − =0 ∂t ∂x2 (6.60) se propone la ecuación en diferencias ui,j+1 − ui,j−1 ui+1,j − 2 [θui,j+1 + (1 − θ)ui,j−1 ] + ui−1,j =0 − 2l h2 (6.61) con 0 < θ < 1. En esta aproximación, el primer término puede interpretarse como una aproximación a la derivada respecto al tiempo, y el segundo como una a la segunda derivada respecto al espacio, en el punto (i, j). En efecto, ambos términos tienden a la derivada que aproximan cuando l o h, respectivamente, tiende a cero. 96 CAPÍTULO 6. LA ECUACIÓN DE DIFUSIÓN Es fácil comprobar, de hecho, que para una función arbitraria u(x, t) ui,j+1 − ui,j−1 |i,j 2l ui+1,j − 2 [θui,j+1 + (1 − θ)ui,j−1 ] + ui−1,j h2 ∂u l2 ∂ 3 u |i,j + ∂t 6 ∂t3 4 +O(l ) (6.62) 2 2 4 ∂ u h ∂ u = + − ∂x2 12 ∂x4 2l ∂u (2θ − 1) 2 + h ∂t l2 ∂ 2 u + h2 ∂t2 l3 (6.63) O(h4 , 2 ) h = Restando una ecuación de la otra se obtiene ui+1,j ui,j+1 − ui,j−1 − 2l − 2 [θui,j+1 + (1 − θ)ui,j−1 ] + ui−1,j = h2 ∂u ∂ 2 u − + ∂t ∂x2 2 3 l ∂ u h2 ∂ 4 u − + 6 ∂t3 12 ∂x4 2l ∂u l2 ∂ 2 u (2θ − 1) 2 + 2 2 + h ∂t h ∂t 3 l O(h4 , 2 , l4 ) (6.64) h Para aproximar con más precisión la ecuación diferencial hace falta que tanto h como l tiendan a cero. Además, por razones de estabilidad, conviene hacer que el cociente entre l y h2 permanezca constante. Por tanto, se supondrá que en todas las mallas k = rh2 , siendo r alguna constante real que asegure estabilidad. Entonces, cuando la anchura de la malla tiende a cero, la ecuación anterior muestra que el lı́mite de la ecuación en diferencias es: ∂u ∂u ∂ 2 u − 2 + 2(2θ − 1)r (6.65) ∂t ∂t ∂t Esta ecuación solamente coincide con la ecuación de difusión original si θ = 1 . En este caso se dice que la ecuación en diferencias es consistente. En los 2 demás casos, en el que el lı́mite no es igual a la ecuación diferencial, se dice que la ecuación es inconsistente. Es importante resaltar que el concepto de consistencia requiere alguna especificación de cómo se reduce la anchura de la malla. Además, el que se tenga 97 6.6. ELEMENTOS FINITOS Y LA ECUACIÓN DE DIFUSIÓN una fórmula que aproxime correctamente cada derivada por separado no es suficiente para asegurar la consistencia, sino que es preciso examinar la ecuación completa. Estabilidad y consistencia son, por tanto, dos condiciones necesarias para asegurar que la ecuación en diferencias converja a la ecuación diferencial conforme se reduce la anchura de malla. Existe un resultado, conocido como el teorema de Lax, que garantiza que en el caso de ecuaciones lineales es también suficiente. 6.6 Elementos finitos y la ecuación de difusión La ecuación (6.1) se puede también resolver mediante el método de elementos finitos. Sin embargo, dado que no es posible derivar esta ecuación de un principio variacional, es preciso obtener la forma débil en la que se basa el método de los elementos finitos de alguna otra manera. Recuérdese que la ecuación de difusión tiene la forma: ∂u ∂t sujeto a ∇(κ∇u) = u∈ u(x, 0) = f (x) u(x, t) = g(x, t) ∂u(x,t) ∂n D × [0, T ] x ∈ ∂D1 = h(x, t) x ∈ ∂D2 (6.66) Escrı́base la función u como una serie u = v0 + X (6.67) ai vi i donde, como es habitual en el método de elementos finitos, v0 cumple las condiciones de contorno esenciales, mientras que las vi se anulan en el contorno esencial ∂D1 . De la ecuación (6.66) se sigue que  ZZZ  ∂u − ∇(κ∇u) φ dv = 0 (6.68) ∂t D para cualquier función φ. En particular, para las funciones vi : ZZZ   ∂u − ∇(κ∇u) vi dv ∂t D i  h   P ZZZ ∂ v0 + j aj vj X  aj vj ) vi dv − ∇(κ∇ v0 + ∂t D j = = 98 CAPÍTULO 6. LA ECUACIÓN DE DIFUSIÓN ZZZ D − ZZZ X ∂v0 vi dv ∂t j ∇(κ∇v0 )vi dv − ZZZ X D ∂aj vj vi dv ∂t aj ∇(κ∇vj )vi dv = j 0 ∀i (6.69) Es costumbre tomar las funciones vi independientes del tiempo, aunque la función v0 puede ser variable si las condiciones de contorno esenciales lo son. P En este caso, en el desarrollo u = v0 + i ai vi los coeficientes ai son variables, y la ecuación anterior se simplifica a:  ∂u − ∇(κ∇u) vi dv ∂t D ZZZ X daj ZZZ ∂v0 vi dv + vj vi dv dt D ∂t D j ZZZ X ZZZ aj ∇(κ∇vj )vi dv − ∇(κ∇v0 )vi dv − ZZZ  = = j X daj j dt Mj,i − ci + X Ij,i = 0 ∀i (6.70) j La matriz definida por Mi,j = ZZZ (6.71) vj vi dv D es conocido como matriz de masas. En cuanto a la matriz Ij,i puede escribirse de un modo más familiar utilizando la fórmula de Green: ZZZ ∇ (vi κ∇vj ) ZZ ∂vj vi κ n ∂D ZZ ∂vj vi κ n ∂D2 ZZZ ∇(κ∇vj )vi ZZZ dv − κ∇vj ∇vi ZZZ ds − κ∇vj ∇vi ZZZ ds − κ∇vj ∇vi dv = dv = dv = dv = dj,i + Kj,i (6.72) donde Ki,j es la matriz de rigidez introducida en el capı́tulo 4. Substituyendo en la ecuación base  X daj Mi,j − Ki,j aj − (di,j aj + ci ) = 0 (6.73) dt j 6.6. ELEMENTOS FINITOS Y LA ECUACIÓN DE DIFUSIÓN 99 Por otra parte X di,j aj + ci = j X j aj ZZ ZZZ ∂vj ds + ∇(κ∇v0 )vi dv ∂n P ZZ j aj ∂vj + ∂v0 ds+ vi κ ∂n ∂D2 ZZ ∂u vi κ ds+ ∂n ∂D2 ZZ vi h ds vi κ ∂D2 = = = = bi (6.74) ∂D2 nótese que este es el coeficiente que aparece como término independiente en el capı́tulo 4. Finalmente, se obtiene X Mi,j j o, en forma matricial X daj Ki,j aj + bi =− dt j (6.75) da = −Ka + b (6.76) dt Esta ecuación diferencial puede integrarse para obtener la solución aproximada a la ecuación original. Las condiciones iniciales vienen dadas por la aproximación a la condición inicial u(x, 0) = f (x). La matriz de rigidez K es positiva definida, según se demostró en el capı́tulo 4. También lo es la matriz de masas M . Por tanto, la matriz M −1 K que aparece en el sistema M da = −M −1 Ka + M −1 b (6.77) dt es positiva definida. Por tanto, el algoritmo de Crank-Nicholson, aplicado a esta ecuación debe ser estable. Es por ello que es una opción muy popular para resolver este sistema. La matriz del algoritmo es G  −1 1 −1 1 −1 = I− M A I+ M A 2 2   −1 1 1 = M− A M+ A 2 2  Por lo tanto, solamente es preciso factorizar la matriz M − 12 A. (6.78) (6.79) 100 CAPÍTULO 6. LA ECUACIÓN DE DIFUSIÓN Capı́tulo 7 Problemas armónicos 7.1 Introducción Es frecuente que en problemas dependientes del tiempo las condiciones de contorno y las densidades de carga o de corriente tengan una dependencia senoidal (o cosenoidal) del tiempo. Por ejemplo, la ecuación de difusión toma a menudo la forma particular: ∂u ∂t sujeto a = u∈ u(x, 0) = f (x) u(x, t) = g0 (x) cos (ωt + φg (x)) ∂u(x,t) ∂n ∇(κ∇u) D × [0, T ] x ∈ ∂D1 = h0 (x) cos (ωt + φh (x)) x ∈ ∂D2 (7.1) La constante ω es la frecuencia del problema. Nótese que aunque se han utilizado funciones coseno para describir la dependencia temporal de las condiciones de contorno, podrı́an haberse utilizado igualmente bien funciones seno, sin más que cambiar las fases φ∗ (x) por φ∗ (x) + pi 2. La solución de esta ecuación es la suma de un término transitorio, que se va amortiguando conforme el tiempo final T tiende a infinito, y un término permanente, que presenta una dependencia cosenoidal del tiempo. Una de las maneras más sencillas de demostrar este hecho es utilizar la ecuación de difusión discretizada (6.15): d u(t) = Au(t) + b(t) dt (7.2) Recuérdese que la matriz A representa el operador laplaciano, y es constante; mientras que el vector b(t) lleva la información sobre las condiciones de contorno. Como éstas son cosenoidales, es fácil ver que 101 102 CAPÍTULO 7. PROBLEMAS ARMÓNICOS bi (t) = b0i cos (ωt + φi ) Por otra parte, la solución de (7.2) es Z t exp [A(t − s)] b(s) ds u(t) = exp [At] u(0) + (7.3) (7.4) 0 El primer término tiende a cero cuando t tiende a infinito. En efecto, como se comentó en la sección 6.4, los autovalores de A son negativos, por lo que que la matriz exp [At] tiende a cero cuando t tiende a infinito. Esto significa que el efecto de las condiciones iniciales, que están representadas en el vector u(0), tiende a desaparecer conforme el tiempo crece, siendo el estado al cabo de un tiempo largo independiente de la manera en que se originó, lo que cabe esperar por razones fı́sicas. También significa que la solución u(t) permanece acotada: si los autovalores de A fueran positivos la exponencial crecerı́a sin lı́mite, lo que es fı́sicamente absurdo. En cualquier caso, se obtiene para tiempos grandes: Z t u(t) = exp [A(t − s)] b(s) ds (7.5) 0 7.2 Notación compleja La integral anterior puede resolverse de una forma particularmente sencilla utilizando números complejos. En efecto, bi (t) = b0i cos (ωt + φi ) = ℜ [b0i exp  (ωt + φi )] = ℜ(b̃0i exp ωt) (7.6) (7.7) (7.8) En este capı́tulo se seguirá la convección de que  es la unidad imaginaria, z̃ un número complejo, z̃ ∗ su complejo conjugado y ℜ e ℑ las partes real e imaginaria respectivamente. Entonces b(t) = ℜ(b̃0 eωt )  1 b̃0 eωt + b̃∗0 e−ωt = 2 Utilizando este resultado para calcular la integral u(t) = Z t exp [A(t − s)] b(s) ds 0 = Z 0 t exp [A(t − s)]  1 b̃0 eωs + b̃∗0 e−ωs ds 2 (7.9) (7.10) (7.11) 103 7.2. NOTACIÓN COMPLEJA = = t  1 b̃0 exp [A(t − s) + ωs] + b̃∗0 exp [A(t − s) − ωs] ds 0 2  −1  −1  ωt b̃0 [A − ω] e − eAt + 2  −1  −ωt (7.12) b̃∗0 [A + ω] e − eAt Z Como eAt tiende a cero cuando t tiende a infinito, para tiempos grandes esta expresión se puede simplificar a: u(t) siendo  −1  −1 −1 b̃0 [A − ω] eωt + b̃∗0 [A + ω] e−ωt 2
1 = ũ0 eωt + ũ∗0 e−ωt 2
= ℜ u0 eωt = −1 (7.13) (7.14) (7.15) ũ0 = −b̃0 [A − ω] (7.16) ω ũ0 = Aũ0 + b̃0 (7.17) De otra forma Comparando con la ecuación original (7.2), d u(t) = Au(t) + b(t) dt (7.18) se observa que para pasar de la una a la otra se han hecho los siguientes cambios: 1. Se han substituido la función del tiempo b(t) por el vector de números complejos b̃0 . Estos números complejos se calculan mediante las ecuaciones (7.6,7.7,7.8). Análogamente con la función u(t) y el vector ũ0 . 2. Se ha substituido la derivada respecto al tiempo por el producto por ω. △ Básico Estas reglas son generales para obtener la respuesta permanente en régimen senoidal de ecuaciones diferenciales lineales. ▽ Avanzado Por ejemplo, en el caso de que se prefiriera resolver la ecuación de difusión mediante elementos finitos, habrı́a que pasar de da(t) = −Ka(t) + b(t) dt (7.19) ωM ã0 = −K ã0 + b̃0 (7.20) M a 104 7.3 CAPÍTULO 7. PROBLEMAS ARMÓNICOS Difusión magnética e histéresis Un problema frecuente es el cálculo de campos magnéticos y eléctricos en circunstancias donde las corrientes de desplazamiento son despreciables. La cuarta ecuación de Maxwell es entonces: ∇×H=j (7.21) Aplicando la ley de Ohm ∇ × H = σE (7.22)   ∂A (7.23) = σ ∇V − ∂t ∂A (7.24) = js − σ ∂t El término js son las corrientes impuestas por el potencial V , que suele ser dato del problema. Por otra parte, como B = µH,   1 ∂A ∇ × A = js − σ (7.25) ∇× µ ∂t En medios ferromagnéticos, el número µ no es una constante, sino que depende del propio campo B, debido a la existencia de fenómenos de saturación e histéresis. Esto causa que, aún cuando el campo H tenga una dependencia senoidal del tiempo, el campo B (y por tanto A) no lo tenga, haciendo imposible la aplicación de los métodos de la sección anterior, que se basa en que todas las magnitudes tengan dependencia senoidal. Existe, no obstante, una forma aproximada de tratar el problema de la histéresis. Esta consiste en la substitución del ciclo real de histéresis por un ciclo aproximado elı́ptico (ver figura 7.1). Esta aproximación es equivalente a la utilización de una constante µ compleja. En efecto, supóngase que H presenta una dependencia senoidal del tiempo H(x, t) = H0 (x) cos (ωt + φH (x)) (7.26) De acuerdo con los métodos de la sección anterior, se puede asociar a este campo, el campo complejo no dependiente del tiempo: H̃0 = H0 eφH (7.27) Se puede calcular ahora el campo B̃0 multiplicando por la constante compleja µ̃: B̃0 = = = = µ̃H̃0 µeφµ H0 eφH µH0 e(φµ +φH ) B0 eφB (7.28) (7.29) (7.30) (7.31) 105 7.3. DIFUSIÓN MAGNÉTICA E HISTÉRESIS B H Figura 7.1: Ciclos de histéresis real y aproximado Este campo complejo tiene asociado el campo real dependiente del tiempo B(x, t) = B0 (x) cos (ωt + φB (x)) (7.32) Ahora, a partir de (7.26) y (7.32) B(x, t) = B0 cos (ωt + φB ) = µH0 cos (ωt + (φµ + φH )) = µH0 cos (ωt + φH ) cos φµ − µH0 sin (ωt + φH ) sin φµ s  2 H H = µH0 sin φµ (7.33) cos φµ − µH0 1 − H0 H0 Luego µH0 sin φµ µ2 H02 sin2 φµ s 2 H = µH cos φµ − B H0  2 ! H 1− = µ2 H 2 cos2 φµ + B 2 − 2µHB cos φµ H0 1−  Ya que al ser los vectores B y H son paralelos, BH = BH. Ası́ pues: 106 CAPÍTULO 7. PROBLEMAS ARMÓNICOS µ2 H02 sin2 φµ − µ2 H 2 sin2 φµ µ2 H 2 + B 2 − 2µHB cos φµ = µ2 H 2 cos2 φµ + B 2 − 2µHB cos φµ = µ2 H02 sin2 φµ (7.34) que es claramente la ecuación de una elipse en el pano B − H. Pero esto es lo que se querı́a demostrar. Nótese que con este método se ha modelado la histéresis, pero no la saturación. En efecto, al aumentar H̃0 , también aumenta proporcionalmente B̃0 = µ̃H̃0 , en vez de aumentar a un ritmo menor al saturarse el material. Con estas cualificaciones, la ecuación (7.25) se transforma en   1 ∇ × Ã0 = j̃s0 − σ Ã0 ∇× (7.35) µ̃ Ahora, no queda más que considerar las condiciones de contorno en Ã0 , discretizar las derivadas respecto al espacio, y resolver la ecuación lineal compleja resultante. Bibliografı́a La bibliografı́a existente sobre la resolución de ecuaciones de campo es sumamente extensa, y ciertamente no existe aquı́ el propósito de presentarla de una forma exhaustiva1 , sino más bien indicar algunos textos en los que el lector interesado puede ampliar conocimientos. Un libro excelente de introducción a la resolución moderna (como contrapuesta a la de preguerra) de sistemas lineales es: • “Linear Algebra and its Applications”, G. Strang, Saunders College Publishing. Y para los más exigentes, la referencia es: • “Matrix Computations”, G. H. Golub and C. F. van Loan, John Hopkins University Press. Una buena introducción a las diferencias finitas, con un gran número de ejemplos trabajados: • “Numerical Solution of Partial Differetial Equations: Finite Difference Methods”, G. D. Smith, Oxford University Press. Uno de los autores de la primera referencia también tiene escrita una introducción muy clara a la teorı́a matemática de los elementos finitos: • “An Analysis of the Finite Element Method”, G. Strang and G. J. Fix, Pretince-Hall. Más completo (o más pesado), tanto desde el aspecto matemático (cubre elementos y diferencias finitas, y otras cosas) como del algorı́tmico, es: • “Numerical Approximation of Partial Differential Equations”, A. Quarteroni and A. Valli, Springer-Verlag. Desde el punto de vista de las aplicaciones a problemas de ingenierı́a, es entretenido: • “What Every Engineer should know about Finite Element Analysis”, J. R. Brauer, Marcel Dekker. 1 Tarea para la que, por otra parte, el autor no se siente capacitado. 107