El presente artículo nace como un mecanismo de difusión para los métodos numéricos en la enseñanza superior para carreras de Ingeniería, Computación y Matemática. Su principal contribución no es la teoría matemática y física que envuelve los métodos numéricos para las Ecuaciones Diferenciales Parciales (EDP). Si no más bien, un documento que recopila todos aquellos aspectos importantes sobre una eficiente y clara implementación computacional de un método numérico de aproximación para EDP. Se espera que lo detallado aquí pueda ser utilizado como apoyo en los cursos de postgrado y por investigadores en sus respectivos trabajos. Es importante aclarar que, si bien, lo presentado aquí puede ser hallado en la literatura, al menos por el conocimiento de los autores, no se cuenta con un escrito que clarifique con gran detalle las técnicas a emplear a la hora de aproximar la solución de EDP por medio del método de diferencias finitas. Adicional a esto, no se cuenta con mucha bibliografía en el idioma español, que describa la implementación computacional del método de diferencias finitas.

En relación con la temática principal de este artículo, se resalta que el análisis y diseño de métodos numéricos para EDP ha sido un área de gran investigación en las últimas décadas, debido a que estas surgen de numerosos problemas de interés como la Ingeniería y la Física (ver, por ejemplo, Haberman (1998), Carslow y Jaeger (1959), y sus referencias). A pesar de lo enriquecedor que es la teoría matemática para la resolución de EDP (ver, por ejemplo, Evans (2010)), lo cierto es que esta no es flexible a la hora de cambiar datos en la ecuación diferencial, ni tampoco brinda, en la mayoría de los casos, la solución explícita de estas. Es debido a ello que la aplicación de métodos numéricos para obtener aproximaciones se convierte, no solo en una gran herramienta, sino en general, en la única para poder resolver un problema que involucre EDP. Aunque claro, llevar los problemas de la matemática a la computación, provoca nuevos inconvenientes. En particular, la dificultad de realizar la implementación de los métodos, lo cual será el punto principal de discusión en este texto.

Para establecer una acotación de la zona en la que se enfoca este trabajo, se debe tener en cuenta que la misma corresponde al uso del Método de Diferencias Finitas (MDF) para una EDP parabólica, conocida como la ecuación de calor, dado que esta modela el flujo de temperatura en una región que está bajo la influencia de una fuente de calor. Esa delimitación es de suma importancia, no solo por la variedad de EDP que se pueden abordar, sino por la cantidad de métodos numéricos existentes para una misma EDP. Para dar un ejemplo de esto, además del enfoque por MDF, el cual consiste en reemplazar cada operador diferencial dentro de la EDP por una estimación puntual del mismo (ver, Ciarlet (1995), LeVeque (2007), Morton y Mayers (2005), Strikwerda (2004), Thomas (1995), y sus referencias); también es posible obtener una aproximación al reemplazar el espacio vectorial donde se ubica la solución exacta por uno de dimensión finita, tal y como ocurre con los Métodos de Elementos Finitos (MEF) (ver Ciarlet (2002)). La diferencia más visual que se presenta entre los MDF y los MEF consiste en que los MEF pueden aproximar la solución de la EDP en dominio con geometrías más complejas que los MDF. A pesar de ello, cuando la EDP depende del tiempo (lo que es usual en los problemas de aplicación), independientemente de si se utiliza MDF o MEF para las variables espaciales, la dependencia del tiempo se discretiza por medio los MDF, dado que los MEF requieren de una relación entre puntos vecinos que resulta inconveniente para la derivada temporal. Esto se puede apreciar en investigaciones recientes como, por ejemplo, en (Guzmán, Shu y Sequeira, 2017).

Por otro lado, basta realizar una búsqueda rápida en MathSciNet para notar que muchas de las investigaciones actuales, que involucran resolver EDP, utilizan el paquete computacional MATLAB^® para sus implementaciones. La principal razón de esto obedece a la potencia y facilidad que brinda este programa a la hora de emplear métodos numéricos en un determinado problema. Sin embargo, ningún programa computacional cuenta con todos los métodos numéricos que existen, sobre todo los nuevos, por lo que se hace necesario que los interesados en resolver una EDP deban programar, lo cual no es una tarea fácil, en la mayoría de los casos en. Por tal razón, en este artículo se discuten los aspectos más relevantes a la hora de realizar una eficiente implementación computacional. Estos aspectos se exponen en general para cualquier lenguaje de programación, pero se ejemplificarán con ayuda de MATLAB^®.

Finalmente, respecto a la organización de este artículo, se inicia con el problema modelo a desarrollar, el cual corresponde a la ecuación de calor en dos dimensiones con condiciones de Dirichlet no homogéneas y cuya solución depende del tiempo. El dominio de definición de la ecuación será un rectángulo por lo que la aplicación del método de diferencias finitas es idónea en este caso. Luego, se detalla cómo la aproximación por esquemas de diferencias finitas en los operadores diferenciales permite obtener un sistema de ecuaciones lineales, cuya solución corresponde a la evaluación de esta buscada en ciertos puntos particulares. Se precisa, además, la construcción (o ensamblaje) eficiente de este sistema lineal para poder implementar la misma en cualquier lenguaje de programación. En particular se presenta ese ensamblaje con ayuda de MATLAB^®. Seguidamente, se consideran aspectos sobre los métodos a utilizar para resolver los sistemas lineales que se generan en cada paso de tiempo. Al final, se muestran algunos experimentos numéricos para validar las programaciones propuestas. Es decir, se comprueba el funcionamiento correcto del código generado. Más aún, se anexa un procedimiento, por medio de videos, para visualizar el comportamiento de la solución obtenida.

Método de diferencias finitas para la ecuación de calor

Por simplicidad en la descripción del modelo a trabajar, defina Ω como el cuadrado unitario en IR², definido como:

Asimismo, sobre Ω se considera la siguiente ecuación parabólica con condiciones de Dirichlet no homogéneas, conocida como la ecuación de calor:

donde corresponde a la derivada temporal, es el operador del laplaciano, el término fuente y T > 0. El objetivo de este escrito es hallar una aproximación para la solución del sistema (1), para lo cual se emplea el método de diferencias finitas (ver, por ejemplo, Morton y Mayers (2005), Strikwerda (2004) o LeVeque (2007)). Más precisamente, dado un entero M > 0, se considera una partición uniforme del intervalo [0, T], dada por:

Donde , se conoce como el paso en tiempo. Seguidamente, dado un entero N > 0, se construye una malla cartesiana (o cuadrícula) uniforme para Ω, la cual consiste de puntos , donde y , con y es el paso en espacio. Es importante aclarar que esta descomposición de Ω no tiene por qué se ser uniforme, sin embargo, se sigue de esta forma con el único fin de simplificar los cálculos posteriores.

Por otro lado, considere el punto el cual al ser evaluado en la ecuación diferencial se obtiene que:

donde, al definir las aproximaciones:

así como , en conjunto con el esquema de 5 puntos de esténcil para el laplaciano (ver, por ejemplo, Morton y Mayers (2005), Strikwerda (2004) o LeVeque (2007)):

se concluye que

Con respecto al espacio, este esquema relaciona a cinco puntos de la malla para estimar tal y como se muestra en la Figura 1. Mientras que, con respecto al tiempo, se tiene una ecuación diferencial ordinaria, cuya solución será aproximada por medio del método de Euler implícito (ver, por ejemplo, Strikwerda (2004)). La razón de emplear un esquema implícito es para generar un método que sea incondicionalmente estable cuando se avanza en cada paso de tiempo. Es decir, no hay condición CFL (Courant-Friedrichs-Lewy) en este caso.

Nota: Fuente propia de la investigación.

El siguiente paso consiste en emplear el método de Euler implícito (conocido también como “backward Euler”) para la derivada temporal. En efecto, aproxima el valor de la derivada respecto del tiempo, mediante el esquema de diferencias finitas:

,

donde, al sustituir este en (3), se obtiene el esquema de aproximación:

o bien,

Luego, despejando a la derecha los términos de la forma: (es decir, los de potencia n), se reescribe la igual previa de la forma:

Por lo tanto, denotando , se cumple que:

para todo y para todo . Es decir, la ecuación (4) establece una relación entre los valores de la solución de la ecuación de calor en dos aproximaciones de tiempo consecutivos. Por ende, la resolución del problema (1) requiere primero aproximar los valores de u en el primer tiempo t0 y luego emplear el resultado obtenido para alcanzar los valores en el tiempo siguiente. Sin embargo, nótese de la tercera ecuación de (1) que los valores de u en el primer tiempo t0 son conocidos a priori de manera exacta, es por ello, que a la ecuación se le conoce como la condición inicial. Más aún, siguiendo la notación previa, esta condición se puede reescribir de la forma:

Por último, de la condición de borde (segunda ecuación en (1)), se sabe que , con alguna función conocida y continua en la frontera de Ω. Así, empleando nuevamente la notación previa se tiene que:

Esquema de diferencias finitas

De la aproximación por diferencias finitas (4), en conjunto con la condición inicial (5) y la condición de borde (6), es posible apreciar que los puntos donde se conoce los valores de u, así como aquellos donde no se conoce (valores a determinar o incógnitas), son los ilustrados en la Figura 2.

Nota: Fuente propia de la investigación.

Ahora, se debe analizar el esquema (4) para los puntos dados por la condición de borde (6), teniendo en cuenta la relación presentada en la Figura 1. Para ello, nótese que:

Luego, para los puntos en los que se desea aproximar la solución (para un tiempo fijo), se pueden ordenar las incógnitas, en el nivel de filas según la Figura 2, en un vector columna en el nivel de bloques, de la forma:

Para todo . De esta manera es posible encontrar el vector , al resolver un sistema de ecuaciones lineales: , donde

Es una matriz de tamaño con I la matriz identidad de orden N, O la matriz nula de orden N, y T una matriz tridiagonal de orden N, definida por

Finalmente, el vector de la derecha del sistema viene dado por:

Para todo . En el caso de y , se tiene, respetivamente, que:

En resumen, es posible encontrar la solución en cada tiempo , con , mediante la resolución de un sistema lineal de orden . Observe que la matriz no depende de , es decir, es independiente del tiempo, lo que implica que en cada paso de tiempo se debe resolver un sistema lineal con la misma matriz de coeficientes. Esta última característica es de gran relevancia a la hora de optar por la eficiencia en la resolución del sistema lineal de interés.

Ensamblaje de la matriz de coeficientes

Tal y como se detalla previamente, la matriz de coeficientes no depende del tiempo, por lo que esta se puede ensamblar una única vez, a priori al recorrido de los pasos de tiempo . De esta forma, basta con ensamblar el vector de la derecha en cada , con , para luego resolver el sistema lineal. En el caso particular del ensamble de , es importante notar que:

Donde, se considera la matriz de orden , definida por:

Aquí, el símbolo ⊗ corresponde al producto tensorial entre matrices, definido como:

Para y , denominado el producto de Kronecker. Gracias a este, la construcción de la matriz se realiza de manera natural (por bloques), lo que además facilita las implementaciones con el paradigma de programación en paralelo (ver, por ejemplo, Quinn (2003)). Así, para ensamblar basta con ensamblar la matriz definida en (10) y junto con la matriz identidad y la matriz definida en (8), todas ellas de tamaño , y efectuar la suma de dos productos de Kronecker: .

Si se utiliza MATLAB® para llevar acabo la implementación del método propuesto en este escrito, entonces, dado que este posee estructuras para la manipulación de matrices esparcidas (o ralas), es recomendable almacenar la matriz en dichas estructuras con el fin de mejorar significativamente la eficiencia del método. En efecto, por medio de los comandos para matrices esparcidas de MATLAB^®, es posible ensamblar la matriz mediante la función:

Donde la rutina realiza l producto de Kronecker: .

Ensamblaje del vector de la derecha

A diferencia de la matriz , la cual se calcula una sola vez, el vector de la derecha , definido en (9), debe ser ensamblado en cada paso de tiempo , con . Más aún, este depende del término fuente , el dato de frontera , así como de la solución del tiempo anterio . Es importante recordar aquí que la solución en el primer tiempo , denominada , es definida a través de la condición inicial, es decir, .

A diferencia del ensamblaje de , donde se toma ventaja de su estructura e independencia respecto de n, el vector de la derecha se ensambla exactamente como fue definido en (9). Más precisamente, definiendo en MATLAB^® un tiempo fijo, con , así como la solución previa , se puede determinar el vector por medio de la siguiente función:

Donde x y y representan los vectores que almacenan la partición en el espacio, mientras que las funciones f y g, son implementaciones respectivas de la fuente y la información de borde , las cuales pueden ser implementadas en MATLAB®, respectivamente, de la siguiente manera:

Con ayuda del fragmento de código previo es posible ensamblar el vector , y, dada la matriz , basta con resolver el sistema lineal , para obtener una aproximación de la solución del problema (1) en el tiempo .

Evolución en el tiempo

Ahora que se conoce cómo ensamblar el sistema lineal para cada , se debe generar el algoritmo principal que recorre cada tiempo , desde , obteniendo estos sistemas, para luego proceder a resolverlos, por medio de un método directo o uno iterativo. Para iniciar, considere las descomposiciones del espacio y el tiempo, respetivamente, de la forma:

donde, se ignoran y , debido a que son valores para los cuales la solución es completamente conocida, según la Figura 2. Por otro lado, es necesario definir la solución inicial para el tiempo la cual es dada por la condición inicial , tal y como se muestra en el siguiente fragmento de código:

Donde la rutina u0 se puede implementar de la forma:

Finalmente, basta con recorrer cada paso en tiempo , , ensamblando el sistema lineal , y como se realiza en la siguiente rutina de MATLAB®:

Observación

Para validar este fragmento de código, basta con dar una función y determinar la función para la cual es solución del problema (1). Para luego verificar que la solución obtenida por el algoritmo anterior concuerda (bajo alguna medida del error) con la solución exacta u dada inicialmente. Por ejemplo, entre varias de las funciones u con las que se validó este algoritmo se cuenta con:

con la cual:

y claramente, , mientras que . Los resultados se pueden apreciar más adelante en la sección de experimentos numéricos.

Resolución del sistema lineal

Según el análisis de las secciones previas, con el fin de aproximar la solución del problema (1), se debe resolver el sistema lineal , en cada paso de tiempo , donde está definida en (7), mientras que está definido en (9). Luego de obtener dicho sistema, es necesario resolverlo por algún método directo o iterativo. En el fragmento de código anterior se emplea la instrucción u = A \ b, correspondiente al método de Eliminación Gaussiana, el cual es un método directo poco estable y poco eficiente. Aunque cabe mencionar que la implementación de MATLAB^® es de las más eficientes que existen, dado que además realiza los cálculos en paralelo (ver, por ejemplo, Hahn y Valentine (2016)).

A continuación, se detallan brevemente algunas opciones eficientes para resolver los sistemas lineales en cada paso de tiempo, tomando provecho principal de la estructura particular de la matriz de coeficientes y del hecho de que esta es invariante en cada paso. Para mayor información sobre los métodos de resolución enunciados aquí, se puede consultar literatura relacionada con álgebra lineal numérica, como, por ejemplo: Datta (2010), Demmel (1997) y Saad (2003).

Métodos directos

Además del método de Eliminación Gaussiana, es posible aplicar dos variantes adicionales, con el fin de acelerar la resolución de los sistemas lineales involucrados en el algoritmo previo. La primera variante se obtiene al observar que el sistema lineal tiene la particularidad de que la matriz es siempre la misma, para todo . Con ello, si se calcula inicialmente la factorización LU de la matriz y se usa esta en la solución del sistema, es posible reducir el tiempo de ejecución significativamente. Por otro lado, la segunda variante, consiste en notar que la matriz es simétrica definida positiva, por lo que en lugar de la factorización LU, se puede utilizar la factorización de Cholesky, la cual tiene como ventaja adicional la reducción del almacenamiento, dado que ahora se almacena en memoria solo una matriz triangular en lugar de dos.

En el caso de la factorización LU, esta puede ser calculada en MATLAB^® a través de la instrucción:

donde L es una matriz triangular inferior con unos en su diagonal y U es una matriz triangular superior, tales que A = LU. De esta forma, el sistema lineal Au = b, puede ser resuelto mediante el comando:

Luego, de manera similar, la instrucción de MATLAB^®:

Retorna una matriz triangular superior R, la cual satisface que A = R^tR. Es decir, determina la factorización de Cholesky, donde R:=L^t. Luego de esto, se puede proceder de manera análoga a la de la factorización LU.

Métodos iterativos

Los métodos iterativos para resolver sistemas lineales, aproximan la solución de un sistema , mediante el cálculo de algunos términos de una sucesión que converge a la solución de este. En relación con los métodos directos, los cuales obtienen la solución en una sola iteración, los iterativos son útiles para resolver sistemas que involucran un gran número de variables (del orden de miles), dado que estos no conllevan un alto costo computacional como los directos. Para poder iniciar un método iterativo, se requiere de una aproximación inicial , usualmente definida como el vector nulo, y entonces se calcula un conjunto de aproximaciones sucesivas , para .

Dado que la matriz de coeficientes en el método de diferencias finitas es simétrica definida positiva, entonces es posible emplear métodos iterativos clásicos como: Jacobi, Gauss-Seidel y SOR. En particular, el método de sobre-relajación o método SOR (Succesive Over Relaxation) corresponde en una aceleración del método de Gauss-Seidel en el que, para definir la iteración, se considera , para el cual se resuelve el sistema equivalente . Más precisamente, el método de SOR es un método de punto fijo cuya iteración viene dada por:

Donde es la parte triangular inferior de , la parte triangular superior y la diagonal de . Es importante tener en cuenta que según el Teorema de Ostrowski-Reich, presente por ejemplo en Saad (2003), la iteración de SOR converge siempre que . Más aún, la matriz de coeficientes generada por el método de diferencias finitas admite un parámetro óptimo, denominado , el cual viene dado por:

Donde y corresponde al radio espectral de una matriz, el cual en este caso satisface que:

Con y, el mínimo y máximo valor propio de , respectivamente.

Finalmente, otro método iterativo recomendado para matrices simétricas positivas definidas es el método de Gradiente Conjugado (CG, por sus siglas en inglés), que consiste en determinar la solución de un sistema lineal , mediante el cálculo de una familia de vectores -ortogonales y con ella se puede escribir la solución del sistema como una combinación lineal de esta. En el caso de MATLAB®, este cuenta con comandos predefinidos para ejecutar el método de Gradiente Conjugado Precondicionado (PCG, por sus siglas en inglés), de la forma:

Donde se desea resolver el sistema Ax = b, por medio de una tolerancia relativa al primer residuo tol > 0, con un máximo de iterMax iteraciones. El PCG se conoce como método precondicionado, dado que si se considera la siguiente descomposición para :

Entonces nótese que puede ser usado como un precondicionador, donde se conoce como la factorización incompleta de Cholesky, la cual no necesariamente existe, pero en el caso de hacerlo, se espera que esta ayude a mejorar el condicionamiento espectral de la matriz del sistema , al considerar el nuevo sistema equivalente:

Para obtener la factorización incompleta de Cholesky en MATLAB®, se puede utilizar la instrucción:

(en Windows) o (en Unix),

Donde, en el caso que L pueda ser creada y sea invertible, entonces se puede ejecutar el método de Gradiente Conjugado Precondicionado con la factorización incompleta de Cholesky de la siguiente manera:

Comparación: tiempo de ejecución y experimentos numéricos: Validación de la implementación

Memoria requerida

Considere el siguiente problema de valores iniciales y contorno, para la ecuación de calor sobre el dominio :

Es decir, se considera el problema modelo (1) con los datos: , , y .

Ahora, se procede a determinar la memoria requerida, así como el tiempo de ejecución del algoritmo final presentado en la sección anterior. Se considera y , así como cada uno de los métodos (directos e iterativos) de resolución de sistemas lineales descritos previamente.

En la Tabla 1, se determina el total de memoria creada por las implementaciones, donde no es posible conocer con exactitud la requerida por Eliminación Gaussiana dado que no se puede acceder a la programación MATLAB® de este método. A pesar de ello, se aprecia, tal y como se esperaba, que los métodos iterativos son los que menos memoria requieren. Por ende, cuando la cantidad de incógnitas sea significativamente grande, los resultados en la Tabla 1 sugieren utilizar métodos iterativos, para así no requerir más memoria que la utilizada para almacenar el sistema lineal generado por el método de diferencias finitas. En la Tabla 2, se muestra el tiempo de ejecución que se requiere previo al recorrido de los pasos temporales. Este tiempo considera la construcción de la matriz , así como aquellos procesos que requieren los métodos de resolución, como, por ejemplo, la determinación de factorizaciones. Observe que Eliminación Gaussiana fue el que menos tiempo necesitó, dado que solo se requiere ensamblar el sistema lineal para su ejecución. En la Tabla 3, se muestra el tiempo de ejecución exclusivo de la iteración que recorre cada paso de tiempo, donde se resuelven los sistemas lineales. Este es el tiempo más considerable, dado que depende de la cantidad de puntos considerados para la aproximación de la solución. Luego, nótese de estos resultados, que el método más eficiente corresponde a la factorización de Cholesky, lo cual es consecuencia de que previo a esta iteración ya se realizaron procedimientos que agilizan la resolución de los sistemas lineales. Más aún, es interesante observar que Eliminación Gaussiana y PCG poseen tiempos similares, a pesar de que el primero se realiza en paralelo. Finalmente, en la Tabla 4, se resumen el tiempo total de ejecución, donde nuevamente se aprecia que resolver los sistemas lineales con la factorización de Cholesky resulta ser el procedimiento más veloz, al menos, para los tamaños de problemas considerados. A pesar de ello, es importante observar que el método PCG obtuvo un buen rendimiento, lo que sugiere que utilizar este método no solo no requiere gran capacidad memoria comparado con los otros, sino que además es veloz.

Visualización del error en la aproximación

Además, del estudio de los recursos involucrados en el método propuesto, también es importante corroborar que se está aproximando la solución del problema (1) correctamente. Para ello, en esta sección se ilustra numéricamente el comportamiento del error:

Respecto al tiempo , siempre que , cuando . En este contexto, u representa la solución exacta al problema (1). En particular, se considera el problema (1) que queda completamente definido al considerar la solución exacta:

Para . Es decir, se considera , y .

Luego, considerando y , en la Gráfica 1 se muestra la disminución del error en el tiempo . Observe cómo el error decrece a medida que h disminuye.

Nota: fuente propia de la investigación. Finalmente, considerando y , en la Gráfica 2 se muestra el incremento del error , para . Observe cómo el error aumenta según se avanza en el paso en tiempo.

Nota: fuente propia de la investigación.

Visualización de la solución: Creación de películas en MATLAB ^®

Debido a que el dominio en el problema (1) está en , entonces la solución, respecto al espacio, es una función de tres variables. Sin embargo, la solución depende además del tiempo, por lo que en realidad sería una función de cuatro variables y, por ende, se debe graficar en , lo cual resulta claro que es visualmente imposible en este escrito, e incluso en la computadora. Por tal razón, se considera realizar una película que muestre el comportamiento de la solución en cada instante , para . Así, el objetivo de esta última sección es mostrar cómo crear una película en MATLAB® que permita visualizar la solución del problema (1). Para ello, considere la siguiente modificación del algoritmo presentado en una sección previa:

Nótese que esta nueva versión realiza las mismas funciones que la anterior, pero dejando el método de Cholesky como método de resolución por defecto. Además, retorna la malla en espacio y no la solución en el último tiempo, sino que almacena la solución del tiempo en la columna j de una matriz matU de orden . Esto puede verse como “el guion a ser filmado”.

Ahora, considere el siguiente conjunto de instrucciones que permiten realizar la creación de la película:

Con ayuda de este fragmento de código, se ejecuta el método y se crea la película para visualizar la solución obtenida. En efecto, considere en el problema (1) los valores:

En la Figura 3, se presentan algunas imágenes de la película que genera las instrucciones previas.

Figura 3. Tomas del video, para .

Nota: Fuente propia de la investigación.

Referencias

Fechas de Publicación

Histórico