ISSN Electrónico: 2215-3470
DOI: http://dx.doi.org/10.15359/ru.33-1.7		 UNICIENCIA Vol. 33, N° 1. Enero-Junio, 2019
URL: www.revistas.una.ac.cr/uniciencia
Correo electrónico: revistauniciencia@una.cr		 logos
uniciencia

Aspectos computacionales del método de diferencias finitas para la ecuación de calor dependiente del tiempo

Computational aspects of the finite difference method for the time-dependent heat equation

Aspectos computacionais do método de diferenças finitas para a equação de calor dependente do tempo

Filánder Sequeira-Chavarría

filander.sequeira@una.ac.cr

Escuela de Matemática

Universidad Nacional

Heredia, Costa Rica

Orcid: https://orcid.org/0000-0002-0593-3446

Melvin Ramírez-Bogantes

meramirez@itcr.ac.cr

Escuela de Matemática

Instituto Tecnológico de Costa Rica

Cartago, Costa Rica

Orcid: https://orcid.org/0000-0001-5516-0085

Recibido-Received: 1/abr/2018 • Corregido-Corrected: 4/jun/2018.

Aceptado-Accepted: 13/ago/2018 • Publicado-Published: 31/ene/2019.


Resumen

En este artículo se describe en detalle un algoritmo para la eficiente implementación computacional del método de diferencias finitas (MDF) en la ecuación de calor dependiente del tiempo, con condiciones de frontera de Dirichlet no homogéneas, en dos dimensiones. Para validar el método presentado aquí se utiliza el paquete computacional MATLAB®, sin embargo, los procesos se exponen independientes al lenguaje de programación. Finalmente se presentan resultados numéricos que validan el algoritmo propuesto.

Palabras claves: ecuación de calor; método de diferencias finitas; implementación computacional; MATLAB®

Abstract

In this paper we describe in detail an algorithm for the efficient computational implementation of the finite difference method (FDM) in the two-dimensional time-dependent heat equation with non-homogeneous Dirichlet boundary conditions. The MATLAB® software was used to validate the method mentioned here; however, the processes are presented independently from the programming language. Finally, numerical results are presented to validate the proposed algorithm.

Keywords: Heat equation; finite difference method; computational implementation; MATLAB®

Resumo

Neste artigo é descrito detalhadamente um algoritmo para a eficiente implementação computacional do método de diferenças finitas (MDF) na equação de calor dependente do tempo, com condições de fronteira de Dirichlet não homogêneas, em duas dimensões. Para validar o método apresentado aqui é utilizado o pacote computacional MATLAB®, porém, os processos são expostos independentes da linguagem de programação. Finalmente, são apresentados os resultados numéricos que validam o algoritmo proposto.

Palavras-chaves: Equação de calor; método de diferenças finitas; implementação computacional; MATLAB®


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 IR2, definido como:

Ω∶=[0,1] × [0,1] = {(x,y) ∈ IR2 / 0 ≤ 𝑥 ≤ 1, 0 ≤ y ≤1}

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

{

𝑢t = ∆𝑢 + 𝑓,

para (x, y)∈ Ω, 𝑡 ∈ ]0, 𝑇],

u (t, x, y) = g (t, x, y),

para (x, y)∈ 𝜕Ω, 𝑡 ∈ [0, 𝑇]

(1)

u (0, x, y) = u0 (x, y),

para (x, y)∈ Ω,

donde ut:= corresponde a la derivada temporal, ∆𝑢 := 𝑢xx+𝑢yy= + es el operador del laplaciano, f = f (t, x, y) 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:

tn:= nk, para todo n = 0,1,…, M,

donde k = ∆t = , 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 (N + 2)2 puntos (𝑥𝑖, 𝑦j ) , donde 𝑥𝑖:= 𝑖ℎ y 𝑦j := 𝑗ℎ, con 𝑖, 𝑗 = 0, 1,…, 𝑁+1 y h= 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 (tn , x𝑖 , yj) el cual al ser evaluado en la ecuación diferencial se obtiene que:

ut (tn ,xi ,yj) = uxx (tn ,xi ,yj) + uyy (tn , xi, yj) + f (tn ,xi ,yj),

donde, al definir las aproximaciones:

≈ u (𝑡n , 𝑥i, 𝑦j ) y ≈ ut (𝑡n , 𝑥i, 𝑦j),

así como := f (tn , 𝑥i , 𝑦j), 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)):

∆u (tn , xi , yj) ≈ (+ – ++),

(2)

se concluye que

= (+ – + + )+

(3)

Con respecto al espacio, este esquema relaciona a cinco puntos de la malla para estimar ∆𝑢 (𝑡n, 𝑥i, 𝑦j ) 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.

Figura 1. Los cinco puntos involucrados en la aproximación del laplaciano.

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: Un (es decir, los de potencia n), se reescribe la igual previa de la forma:

- - + (4+) - - = + .

Por lo tanto, denotando μ:=, se cumple que:

– – + (4+μ) – – = + μ,

(4)

para todo i, j =1 ,2,…, N y para todo n = 1, 2,…, M. 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 u(0, x, y) = u0(x, y) 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:

= u0 (xi , yj), para 𝑖, 𝑗 = 0, 1,2, . . . , 𝑁, 𝑁 + 1. (5)

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

=𝑔(𝑡n ,𝑥i ,𝑦j)=:𝑔, para 𝑖, 𝑗∈{0, 𝑁+1},
𝑡 = 0,1,…, 𝑀.

(6)


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.

Figura 2. Puntos dados (condición inicial y de borde) e incógnitas por calcular.

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:

– – + (4+ μ) – – = + μ

–𝑔 – + (4+ μ) – – = + μU

+ (4+ μ) – – = + μU+𝑔

– – +(4+ μ) – – =+

– –+(4+μ) – –𝑔=+

– – +(4+μ) – =++𝑔

– – +(4+μ) – – = +

– –𝑔 +(4+μ) – – = +

– +(4+μ) – – = + +𝑔

– – +(4+μ) – – = +

– – +(4+μ) – – = +

– – +(4+μ) – =++

Luego, para los N2 puntos en los que se desea aproximar la solución (para un tiempo tn 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:

, donde

,

Para todo j =1, 2,…, N. De esta manera es posible encontrar el vector un ∈ IRN2, , al resolver un sistema de ecuaciones lineales: Aun=bn , donde

(7)

Es una matriz de tamaño N2 × N2 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

(8)

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

, donde

(9)

Para todo j = 2, 3,…, N - 1. En el caso de j = 1 y j = N, se tiene, respetivamente, que:

Y

En resumen, es posible encontrar la solución un en cada tiempo tn, con n = 1, 2, … , M, mediante la resolución de un sistema lineal Aun=bn de orden N 2. Observe que la matriz A no depende de n, 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 A no depende del tiempo, por lo que esta se puede ensamblar una única vez, a priori al recorrido de los pasos de tiempo tn . De esta forma, basta con ensamblar el vector de la derecha bn en cada tn, con n=1, 2,…, M, para luego resolver el sistema lineal. En el caso particular del ensamble de A,es importante notar que:

+

= (I ⊗T)+(S ⊗I),

donde, se considera la matriz S de orden N, definida por:

S :

=

(10)

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

Para B∈IRmxn y C∈IRpxq, denominado el producto de Kronecker. Gracias a este, la construcción de la matriz A 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 A basta con ensamblar la matriz Sdefinida en (10) y junto con la matriz identidad I y la matriz T definida en (8), todas ellas de tamaño N × N, y efectuar la suma de dos productos de Kronecker: (I ⊗ T ) + (S ⊗ I ) .

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 A 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 A mediante la función:

function [ A ] = Laplacian2DMatrix( N, mu )

% Matriz identidad

I = speye( N );

% Matriz T

e = -1 * ones( N, 1 );

T = spdiags( [e, -(4 + mu)*e, e], [-1, 0, 1], N, N );

% Matriz S

S = spdiags( [e, e], [-1, 1], N, N );

% Matriz A

A = kron( I, T ) + kron( S, I );

end

Donde la rutina kron(B,C)realiza l producto de Kronecker: B ⊗ C.

Ensamblaje del vector de la derecha

A diferencia de la matriz A, la cual se calcula una sola vez, el vector de la derecha bn, definido en (9), debe ser ensamblado en cada paso de tiempo tn, con n = 1,2,…, M. Más aún, este depende del término fuente := f(tn, xi , yj), el dato de frontera :=𝑔(tn,xi ,yj), así como de la solución del tiempo anterior un-1. Es importante recordar aquí que la solución en el primer tiempo t0, denominada u0, es definida a través de la condición inicial, es decir, : = u0 (xi , yj) .

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

function [ b ] = Lapacian2DRhs( N, tn, x, y, h, f, g, mu, u )

b = zeros(N^2,1);

k = N^2 - N;

v = y(1); % para j = 1

w = y(N); % para j = N

for i = 1 : N

z = x(i);

b(i) = h^2 * f(tn, z, v) + g(tn, z, 0);

b(k + i) = h^2 * f(tn, z, w) + g(tn, z, 1);

end

b(1) = b(1) + g(tn, 0, v);

b(N) = b(N) + g(tn, 1, v);

b(k+1) = b(k+1) + g(tn, 0, w);

b(k+N) = b(k+N) + g(tn, 1, w);

k = N; % para j = 2, 3 ,..., N -1

for j = 2 : N-1

v = y(j);

for i = 1 : N

b(k + i) = h^2 * f(tn, x(i), v);

end

b(k+1) = b(k+1) + g(tn, 0, v);

k = k + N;

b(k) = b(k) + g(tn, 1, v);

end

b = b + mu * u; % contribución de la solución anterior

end

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 f y la información de borde 𝑔, las cuales pueden ser implementadas en MATLAB®, respectivamente, de la siguiente manera:

function [ z ] = f( t, x, y )

z = 0; % aquí el criterio de la función

return

function [ z ] = g( t, x, y )

z = 0; % aquí el criterio de la función

end

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

Evolución en el tiempo

Ahora que se conoce cómo ensamblar el sistema lineal para cada n = 1, 2,…, M, , se debe generar el algoritmo principal que recorre cada tiempo tn , desde n = 1 a n = M, 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:

% Malla uniforme para [0,1]x[0,1]

h = 1 / ( N + 1 );

x = h : h : ( 1 - h );

y = x;

% Malla uniforme para [0,T]

k = T / M;

t = k : k : T;

donde, se ignoran x0 = y0 = 0 xN+1 = yN+1 = 1 y t0 = 0, 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 u=u0 para el tiempo t0= 0 la cual es dada por la condición inicial u0(x, y), tal y como se muestra en el siguiente fragmento de código:

u = zeros( N^2, 1 );

for j = 1 : N

for i = 1 : N

u( (j-1)*N + i ) = u0( x(i), y(j) );

end

end

donde la rutina u0 se puede implementar de la forma:

function [ z ] = u0( x, y )

z = 0; % aquí el criterio de la función

end

Finalmente, basta con recorrer cada paso en tiempo tn, n = 1, 2,…M, ensamblando el sistema lineal Aun = bn , y como se realiza en la siguiente rutina de MATLAB®:

% Discretización de la ecuación de calor

% u_t = u_xx + u_yy + f, (x,y) en [0,1]x[0,1]

% u = g(t, x, y), (x,y) en la frontera

% u = u_0(x, y), para t = 0

% donde t pertenece a [0,T] y N es el número de puntos en

% espacio, mientras que M es el número de puntos en tiempo.

function [u] = HeatEquation2DScheme( T, N, M, f, g, u0 )

% Malla uniforme para [0,1]x[0,1]

h = 1 / (N + 1);

x = h : h : (1 - h);

y = x;

% Malla uniforme para [0,T]

k = T / M;

t = k : k : T;

% Solución en el tiempo t = 0

u = zeros(N^2, 1);

for j = 1 : N

for i = 1 : N

u( (j-1)*N + i ) = u0( x(i), y(j) );

end

end

% Obtener la matriz

mu = h^2 / k;

A = Laplacian2DMatrix( N, mu );

% Recorrer cada paso de tiempo

for n = 1 : M

b = Lapacian2DRhs( N, t(n), x, y, h, f, g, mu, u );

u = A \ b; % Resolver el sistema lineal

end

end


Observación:

Para validar este fragmento de código, basta con dar una función u (t, x, y) y determinar la función f (t, x, y) para la cual u 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:

u (t, x, y) : = sen (𝑡 + 𝜋 𝑥) sen (𝑡 + 𝜋 y),

con la cual:

f (t, x, y) : = cos (t + 𝜋𝑥) sen (𝑡 + 𝜋y) + sen (𝑡 + 𝜋𝑥) (cos (𝑡 + 𝜋y) + 2𝜋2 sen (𝑡 + 𝜋y)),

y claramente, g (t, x, y) = u (t, x, y), mientras que u0 (x, y) = u (0, x, y). 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 Aun=bn , en cada paso de tiempo tn , donde A ∈IRN2×N2 está definida en (7), mientras que b ∈ IRN2 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 Aun = bn tiene la particularidad de que la matriz A es siempre la misma, para todo n=1,…, M. Con ello, si se calcula inicialmente la factorización LU de la matriz A 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 A 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:

[L, U] = lu( A );

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:

u = U \ ( L \ b );

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

R = chol( A );

Retorna una matriz triangular superior R, la cual satisface que A = RtR. Es decir, determina la factorización de Cholesky, donde R:=Lt. 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 Ax=b,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 x0, usualmente definida como el vector nulo, y entonces se calcula un conjunto de aproximaciones sucesivas xk , para k = 1, 2,…,NúmeroMáximoIteraciones.

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 w∈ IE, w ≠ 0, para el cual se resuelve el sistema equivalente wAx = wb. Más precisamente, el método de SOR es un método de punto fijo cuya iteración viene dada por:

xk+1=(L+D)-1[b+(() D – U)xk],

donde L es la parte triangular inferior de A, U la parte triangular superior y D la diagonal de A. 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 0<w<2. Más aún, la matriz de coeficientes generada por el método de diferencias finitas admite un parámetro w óptimo, denominado wopt, el cual viene dado por:

donde y ρ (•) corresponde al radio espectral de una matriz, el cual en este caso satisface que:

con 𝜆mín (A) y 𝜆máx(A), el mínimo y máximo valor propio de A, 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 Ax=b, mediante el cálculo de una familia de vectores A-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:

x = pcg(A, b, tol, iterMax);

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 A:

A=LLt + R, con R ≠ O,

Entonces nótese que A ≈ LLt puede ser usado como un precondicionador, donde LLt 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 Ax = b, al considerar el nuevo sistema equivalente:

M-1Ax = M-1b, con M : = LLt.

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

L = ichol( A ); (en Windows) o L = cholinc( A ); (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:

x = pcg(A, b, tol, iterMax, L, L');


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 Ω:=[0,1]2:

Es decir, se considera el problema modelo (1) con los datos: f (t, x, y) : = 0, g (t, x, y) : = 0, u0 (x, y) : = xy (1 – x) (1 – y) y T : = 1.

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 N ∈{19, 99, 499} y M : = N + 1, 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 A, 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.

Tabla 1. Memoria requerida en megabytes (MB)

N

M

Eliminación Gaussiana

Factorización LU

Factorización Cholesky

SOR

PCG

19

20

0.0192

0.1241

0.1241

0.0379

0.0711

99

100

0.5227

15.3299

15.3299

1.0431

2.4909

499

500

13.2944

1909.2407

1909.2407

26.5774

72.7090

Nota: Fuente propia de la investigación.


Tabla 2. Tiempo de cómputo independiente de la iteración temporal (en min)

N

M

Eliminación Gaussiana

Factorización LU

Factorización Cholesky

SOR

PCG

19

20

0.00041

0.00008

0.00002

0.00288

0.00002

99

100

0.00050

0.00465

0.00091

0.01066

0.00102

499

500

0.00506

1.97844

0.18130

3.60424

0.10279

Nota: Fuente propia de la investigación.

Tabla 3. Tiempo de cómputo únicamente de la iteración temporal (en min)

N

M

Eliminación Gaussiana

Factorización LU

Factorización Cholesky

SOR

PCG

19

20

0.00027

0.00009

0.00008

0.00046

0.00107

99

100

0.03561

0.01298

0.01290

0.09668

0.03907

499

500

6.80456

4.07359

4.04021

34.71655

7.38384

Nota: Fuente propia de la investigación.

Tabla 4. Tiempo de cómputo total (en min)

N

M

Eliminación Gaussiana

Factorización LU

Factorización Cholesky

SOR

PCG

19

20

0.00068

0.00017

0.00009

0.00333

0.00109

99

100

0.03611

0.01763

0.01381

0.10734

0.04009

499

500

6.80962

6.05204

4.22151

38.32079

7.48663

Nota: Fuente propia de la investigación.

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 n = 0, 1, 2, …, M, 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:

u (t, x, y) = sen (t + 𝜋x) sen (t + 𝜋y),

Para T = 1. Es decir, se considera g (t, x, y) : = u (t, x, y) , u0 (x, y) : = u (0, x, y) y f (t, x, y) : = cos (t + 𝜋x) sen (t + 𝜋y) + (t + 𝜋x) (cos (t + 𝜋y) + 2 𝜋2 sen (t + 𝜋y))

Luego, considerando N ∈{20, 40, 60, …, 200} y M = 10N, en la Gráfica 1 se muestra la disminución del error eh,m en el tiempo tM = T = 1. Observe cómo el error decrece a medida que h disminuye.

Gráfica 1. Error eh,m para un tiempo fijo

Nota: fuente propia de la investigación.

Finalmente, considerando N = 100 y M = 1000, en la Gráfica 2 se muestra el incremento del error eh,n , para n = 0,1,…,M. Observe cómo el error aumenta según se avanza en el paso en tiempo.

Gráfica 2. Error eh,n para un h fijo.

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 IR2 , 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 IR4, 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 tn, para n=0,1,…,M. 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:

% Discretización de la ecuación de calor

% u_t = u_xx + u_yy + f, (x,y) en [0,1]x[0,1]

% u = g(t, x, y), (x,y) en la frontera

% u = u_0(x, y), para t = 0

% donde t pertenece a [0,T] y N es el número de puntos en

% espacio, mientras que M es el número de puntos en tiempo.

function [matU,x,y] = HeatEquation2DSchemeMovie(T, N, M, f, g, u0)

% Malla uniforme para [0,1]x[0,1]

h = 1 / (N + 1);

x = h : h : (1 - h);

y = x;

% Malla uniforme para [0,T]

k = T / M;

t = k : k : T;

% Solución en el tiempo t = 0

u = zeros(N^2, 1);

for j = 1 : N

for i = 1 : N

u( (j-1)*N + i ) = u0( x(i), y(j) );

end

end

matU = u; % Matriz donde la columna j contiene

% la solución en el tiempo t_j

% Obtener la matriz

mu = h^2 / k;

A = Laplacian2DMatrix( N, mu );

% Factorización de Cholesky

U = chol( A ); L = U';

% Recorrer cada paso de tiempo

for n = 1 : M

b = Lapacian2DRhs( N, t(n), x, y, h, f, g, mu, u );

u = U \ ( L \ b ); % Resolver el sistema lineal

matU = [ matU, u ]; % Almacenar la solución

end

end % Fin del programa

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 tj en la columna j de una matriz matU de orden N2 × (M+1). 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:

N = 99;

M = 300;

T = 10;

% Determinar la solución aproximada

[matU, x, y] = HeatEquation2DSchemeMovie( T, N, M, @f, @g, @u0 );

% Crear película

writerObj = VideoWriter('video.avi');

open( writerObj );

% Control de los ejes coordenados

axis tight

set(gca, 'nextplot', 'replacechildren');

set(gcf, 'Renderer', 'zbuffer');

% Fase de “filmación”

numframes = size(matU, 2);

[X, Y] = meshgrid(x,y);

for n = 1 : numframes

uh = matU(:, n);

U = reshape(uh, N, N);

Z = mesh(X, Y, U);

title(['Soluci\''on aproximada $t_{', num2str(n-1), '}$'],...

'interpreter','latex')

xlabel('$x$','interpreter’,'latex')

ylabel('$y$','interpreter','latex')

zlabel('$t$','interpreter','latex')

currentframe = getframe;

writeVideo(writerObj, currentframe);

end

close( writerObj );

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: g(t, x, y):= 0, u0 (x, y):= x(1–x) y (1–y) y f(t, x, y):= cos(t+πx) sen(t+ πx)(cos(t+ πy) + 2 π2sen(t+ πy)).

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

Figura 4. Tomas del video, para n=50, 100,…,300

Nota: Fuente propia de la investigación.

Referencias


Carslow, H. S., y Jaeger, J. C. (1959). Conduction of Heat in Solids. Segunda edición. Inglaterra: Oxford University Press.

Ciarlet, P. G. (1995). Introduction to numerical linear algebra and optimisation. Estados Unidos: Cambrige University Press.

Ciarlet, P. G. (2002). The finite Element Method for Elliptic Problems. Estados Unidos: SIAM. doi: https://doi.org/10.1137/1.9780898719208

Datta, B. N. (2010). Numerical Linear Algebra and Applications. Segunda Edición. Estados Unidos: SIAM. doi: https://doi.org/10.1137/1.9780898717655

Demmel, J. W. (1997). Applied Numerical Linear Algebra. Estados Unidos: SIAM. doi: https://doi.org/10.1137/1.9781611971446

Evans, L. C. (2010). Partial Differential Equations. Segunda edición. Estados Unidos: American Mathematical Society.

Guzmán, J., Shu, C-W., y Sequeira, F. A. (2017). H(div) conforming and DG methods for the incompressible Euler’s equations. IMA Journal of Numerical Analysis, 37(4), 1733-1771. doi: https://doi.org/10.1093/imanum/drw054

Haberman, R. (1998). Elementary Applied Partial Differential Equations, With Fourier Series and Boundary Value Problems. Tercera edición. Estados Unidos: Prentice-Hall.

Hahn, B. H., y Valentine D. T. (2016). Essential MATLAB for Engineers and Scientists. Sexta edición. Estados Unidos: Elsevier.

LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations, Steady-State and Time-Dependent Problems. Estados Unidos: SIAM. doi: https://doi.org/10.1137/1.9780898717839

Morton, K. W., y Mayers, D. F. (2005). Numerical Solution of Partial Differential Equations. Segunda edición. Inglaterra: Cambrige University Press. doi: https://doi.org/10.1017/CBO9780511812248

Quinn, M. J. (2003). Parallel Programming in C with MPI and OpenMP. Estados Unidos: McGraw-Hill.

Saad, Y. (2003). Iterative Methods for Sparse Linear Systems. Segunda Edición. Estados Unidos: SIAM. doi: https://doi.org/10.1137/1.9780898718003

Strikwerda, J. C. (2004). Finite Difference Schemes and Partial Differential Equations. Segunda edición. Estados Unidos: SIAM. doi: https://doi.org/10.1137/1.9780898717938

Thomas, J. W. (1995). Numerical Partial Differential Equations: Finite Difference Methods. Estados Unidos: Springer-Verlag. doi: https://doi.org/10.1007/978-1-4899-7278-1

Licencia Creative Commons

Aspectos computacionales del método de diferencias finitas para la ecuación de calor dependiente del tiempo (Filánder Sequeira-Chavarría y otros) por Revista Uniciencia se encuentra bajo una Licencia CreativeCommons Atribución-NoComercial-SinDerivadas 3.0 Unported.