Universidad Central de Venezuela

Facultad de Ciencias
Escuela de Matemática
Postgrado de Matemática

Un método mimético de diferencias finitas para la

ecuación no estática de difusión

Trabajo Especial de Grado presentado ante la Ilustre

Universidad Central de Venezuela por la Lic. Iliana
A. Mannarino S., para optar al Tı́tulo de Magister
Scientiarium en Matemática.

Tutores: Dr. J. M. Guevara-Jordan.

Dra. Y. Quintana.

Caracas - Venezuela
Octubre del 2007
A la memoria de mis tios Filiberto y Alberto.

2
Agradecimientos

Antes que nada quisiera agradecer a Dios todopoderoso por haberme dado la pa-
ciencia y voluntad para llegar hasta el final de este camino.
A mis padres y a mi hermana por estar siempre a mi lado. Este es mi otro regalo para
ustedes!!!
A Tony por darme ánimo sin permitir desvanecerme ante las dificultades. Gracias
amor!!!!
A Belkys por apoyarme en los buenos y malos momentos.
A mis tutores por creer en mi y darme la oportunidad de finalizar esta etapa de mi
formación académica.
A la profesora Zenaida Castillo por brindarme sus valiosos conocimientos en todo mo-
mento.
A la UCV por seguir siendo, a pesar de todo, la casa que vence la sombra!!!
Iliana

3
Contenido

Agradecimientos 3

Lista de tablas 4

Lista de figuras 4

Resumen 5

Introducción 6

1. La ecuación no estática de difusión 9

1.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.2. La ecuación no estática de difusión . . . . . . . . . . . . . . . . . . . . 9

2. Esquema mimético tipo Crank-Nicolson 13

2.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2. Definición del método mimético . . . . . . . . . . . . . . . . . . . . . . 13
2.3. Discretización mimética unidimensional de los operadores gradiente, di-
vergencia y derivada direccional en el borde . . . . . . . . . . . . . . . 14
2.4. Discretización mimética bidimensional de los operadores gradiente, di-
vergencia y derivada direccional en el borde . . . . . . . . . . . . . . . 16

4
 2.5. Resolución de la ecuación no estática de difusión por el esquema de
Crank-Nicolson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.6. Esquema mimético tipo Crank-Nicolson . . . . . . . . . . . . . . . . . . 22
2.7. Ecuación no estática de difusión bidimensional discretizada . . . . . . . 23
2.7.1. Discretización mimética de las ecuaciones bidimensionales . . . 24

3. Análisis de convergencia 32
3.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.2. Preliminares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.3. Convergencia del esquema mimético . . . . . . . . . . . . . . . . . . . . 33
3.3.1. Discretización mimética de las ecuaciones unidimensionales . . . 34
3.3.2. Consistencia del esquema mimético tipo Crank-Nicolson . . . . 39
3.3.3. Estabilidad del esquema mimético tipo Crank-Nicolson . . . . . 45
3.3.4. Estabilidad Numérica . . . . . . . . . . . . . . . . . . . . . . . . 51

4. Resultados numéricos 53
4.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.2. Método de diferencias finitas . . . . . . . . . . . . . . . . . . . . . . . . 53
4.3. Análisis 1D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.4. Análisis 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.5. Conclusiones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

Bibliografı́a 62

1
Lista de figuras

1.1. Volumen de control para deducir la ecuación no estática de difusión . . 10

2.1. Malla mimética unidimensional . . . . . . . . . . . . . . . . . . . . . . 14

2.2. Malla mimética bidimensional de orden 5 × 5 . . . . . . . . . . . . . . . 17
2.3. Estructura matricial del operador Divergencia en 2D . . . . . . . . . . 18
2.4. Matriz del operador gradiente en 2D . . . . . . . . . . . . . . . . . . . 18
2.5. Matriz del operador derivada direccional en el borde B en 2D . . . . . . 19
2.6. Molécula computacional del esquema explı́cito . . . . . . . . . . . . . . 20
2.7. Molécula computacional del esquema implı́cito . . . . . . . . . . . . . . 21
2.8. Molécula computacional del esquema de Crank-Nicolson 1D . . . . . . 22
2.9. Molécula computacional del esquema de Crank-Nicolson 2D . . . . . . 22
2.10. Nodos de la malla mimética bidimensional de orden 5 × 5 . . . . . . . . 24
2.11. Molécula computacional correspondiente al nodo N1 . . . . . . . . . . . 25
2.12. Molécula computacional correspondiente al nodo N2 . . . . . . . . . . . 26
2.13. Molécula computacional correspondiente al nodo N3 . . . . . . . . . . . 27
2.14. Molécula computacional correspondiente al nodo N4 . . . . . . . . . . . 28
2.15. Molécula computacional correspondiente al nodo N5 . . . . . . . . . . . 29
2.16. Molécula computacional correspondiente al nodo N6 . . . . . . . . . . . 30
2.17. Molécula computacional correspondiente al nodo N7 . . . . . . . . . . . 31
2.18. Estructura del sistema bidimensional . . . . . . . . . . . . . . . . . . . 31

2
3.1. Nodo ubicado en el extremo izquierdo. . . . . . . . . . . . . . . . . . . 35
3.2. Nodo correspondiente a la segunda ecuación. . . . . . . . . . . . . . . . 35
3.3. Nodo correspondiente a la tercera ecuación. . . . . . . . . . . . . . . . 36
3.4. Nodos correspondientes a las ecuaciones estándard. . . . . . . . . . . . 36
3.5. Nodo correspondiente a la antepenúltima ecuación. . . . . . . . . . . . 37
3.6. Nodo correspondiente a la penúltima ecuación. . . . . . . . . . . . . . . 37
3.7. Nodo correspondiente a la última ecuación. . . . . . . . . . . . . . . . . 38
3.8. Estructura del sistema unidimensional . . . . . . . . . . . . . . . . . . 38
3.9. Radio espectral para un sistema de órden 50 × 50 . . . . . . . . . . . . 51
3.10. Radio espectral para un sistema de órden 5 × 5 . . . . . . . . . . . . . 52

4.1. Malla referencial para el esquema de diferencias finitas 1D . . . . . . . 54

4.2. Malla referencial para el esquema de diferencias finitas 2D . . . . . . . 54
4.3. Comparación de los errores numéricos en norma infinito 1D . . . . . . . 57
4.4. Comparación de los errores numéricos en norma infinito 2D . . . . . . . 60

3
Lista de tablas

3.1. Órdenes de los errores de truncación 1D . . . . . . . . . . . . . . . . . 45

4.1. Errores en norma infinito cometidos por ambos métodos . . . . . . . . 56

4.2. Órdenes de convergencia 1D . . . . . . . . . . . . . . . . . . . . . . . . 56
4.3. Errores en norma infinito cometidos por ambos métodos . . . . . . . . 59
4.4. Órdenes de convergencia 2D . . . . . . . . . . . . . . . . . . . . . . . . 59

4
Resumen

Se presenta un nuevo esquema mimético de diferencias finitas tipo Crank-Nicolson

para resolver la ecuación no estática de difusión, el cual extiende las formulaciones
miméticas para ecuaciones elı́pticas [16] al caso parabólico. Se demuestra que el nuevo
esquema distingue perfectamente las condiciones de borde a nivel de malla, es consi-
stente con la ecuación de difusión, incondicionalmente estable y por lo tanto conver-
gente. Un estudio numérico y comparativo de las tasas de convergencia evidencia que
el nuevo método mimético produce sistemáticamente menores errores que el esquema
de diferencias finitas tradicional, basado en “puntos fantasmas” para las condiciones
de borde. Los resultados indican que la extensión propuesta del método mimético al
caso no estático de la ecuación de difusión representa una mejora sobre los esquemas
de diferencias finitas tradicionales, sin ser más compleja su implementación.

5
Introducción

Muchos procesos de distinta naturaleza son descritos mediante el planteamien-

to de sistemas de ecuaciones en derivadas parciales. Especı́ficamente, la ecuación no
estática de difusión juega un papel fundamental en la descripción matemática de diver-
sos fenómenos fı́sicos, biológicos, económicos y tecnológicos. En Fı́sica, dicha ecuación
es utilizada para modelar la evolución de la temperatura en un cuerpo mediante el
proceso de conducción de calor. Mientras que en Biologı́a la misma ecuación permite
describir la evolución de la población de seres vivos de distintas especies en función de
su distribución espacial. También se puede mencionar el reciente modelo matemático
de los productos financieros llamados “opciones”, el cual viene dado por la llamada
ecuación de Black-Scholes que representa un tipo especial de ecuación no estática de
difusión. Todos estos ejemplos parciales de aplicaciones evidencian la versatilidad de
esta ecuación.
A pesar de todas estas aplicaciones la ecuación no estática de difusión carece de
solución analı́tica en la gran mayorı́a de las situaciones de interés. Esta dificultad puede
tener diversos orı́genes: desde la presencia de fronteras irregulares hasta la existencia
de términos no lineales en las ecuaciones mismas. Con el fin de evitar tal dificultad es
necesario presentar el problema de manera puramente algebraica. Mediante el proceso
de discretización, el conjunto infinito de números que representa las funciones incógni-
tas en el contı́nuo es reemplazado por un número finito de parámetros incógnita y

6
este proceso requiere alguna forma de aproximación. Es por ello que se han propuesto
variados métodos numéricos [3, 11] para aproximar la solución sin que, hasta el momen-
to, exista alguno que pueda ser calificado como el mejor. En consecuencia, el estudio,
análisis e implementaciones de nuevos esquemas numéricos son tópicos de investigación
de interés actual.
Entre los métodos numéricos utilizados para resolver la ecuación no estática de
difusión se encuentran los esquemas de diferencias finitas tradicionales [3, 17, 18]. Es
bien conocido que dichos esquemas presentan deficiencias como por ejemplo el uso
de puntos fantasmas o imaginarios en la discretización de las condiciones de borde y
las inconsistencias existentes entre las discretizaciones de la ecuación diferencial con
las condiciones de borde, siendo ellas la causa de que muchos esquemas de diferencias
finitas tradicionales no sean completamente conservativos. Esto es una severa limitación
para cualquier método numérico. A pesar de ello, los esquemas tradicionales de dife-
rencias finitas siguen siendo los métodos numéricos más populares por su simplicidad
de implementación y sus bajos tiempos de ejecución.
En los últimos diez años se han venido presentando una serie de nuevos métodos
numéricos que se basan en diferencias finitas pero que eliminan las deficiencias de los
esquemas tradicionales. Entre esos nuevos esquemas numéricos se encuentran los lla-
mados métodos miméticos [5, 7, 8], desarrollados con la finalidad de obtener mejores
softwares para la simulación de diversos fenómenos. Entre sus ventajas se pueden men-
cionar: el tratamiento riguroso de las discretizaciones de las condiciones de borde sin
hacer uso de puntos imaginarios, la compatibilidad de las aproximaciones para las
condiciones de borde, las ecuaciones diferenciales son completamente conservativas y
satisfacen versiones discretas de las principales identidades entre los distintos oper-
adores diferenciales del cálculo diferencial.
Recientemente, Castillo-Grone [1, 2] desarrollaron un nuevo tipo de discretización
mimética para los operadores gradiente y divergencia que representa una mejora en el

7
orden de aproximación de las condiciones de borde. En [4, 20] se realizaron diversos
estudios comparativos entre los esquemas miméticos basados en el trabajo de Castillo-
Grone, los métodos de bajo orden en el borde y los métodos de diferencias finitas
tradicionales encontrándose que los primeros produjeron las mejores aproximaciones
a la solución en todos los casos analizados. En dichos trabajos se resolvió la ecuación
estática de difusión. El estudio teórico de convergencia del esquema mimético para la
ecuación estática de difusión ha sido realizado recientemente en [4]. Debido a las ven-
tajas del uso del método mimético en el caso estático, en este trabajo de investigación
se propone la combinación del esquema mimético antes mencionado con el esquema
de Crank-Nicolson [18] para el caso no estático. Parte de los resultados numéricos
obtenidos se pueden encontrar en [12, 13].
El contenido de este trabajo está organizado de la siguiente manera: en el capı́tulo
uno se presenta la deducción de la ecuación no estática de difusión a partir de las
leyes fı́sicas de conservación de la energı́a y de conducción de calor de Fourier, en
el segundo capı́tulo se describen el método mimético y el esquema de resolución de
Crank-Nicolson. Además se define el nuevo esquema mimético tipo Crank-Nicolson
y se dan las discretizaciones de las ecuaciones bidimensionales. Luego, en el capı́tulo
tres se realiza un análisis de convergencia a través del estudio de la consistencia y la
estabilidad del nuevo esquema para el caso unidimensional. Finalmente, en el cuarto
capı́tulo se muestra el estudio numérico-comparativo entre el nuevo esquema mimético
y el método de diferencias finitas tradicional, luego se exponen las conclusiones más
relevantes obtenidas a lo largo del trabajo investigativo.

8
Capı́tulo 1

La ecuación no estática de difusión

1.1. Introducción
Este capı́tulo presenta la deducción de la ecuación no estática de difusión a partir
de las leyes fı́sicas de conservación de la energı́a y de conducción de calor de Fourier.

1.2. La ecuación no estática de difusión

Se considera un sólido tridimensional en el cual el calor fluye libremente. Sea
u (x, y, z, t) la temperatura medida en el punto x = (x, y, z) en el instante t.
Sea q (x, t) una función vectorial que designa a la densidad de conductividad de
calor y representa la razón de flujo de calor por unidad de tiempo. La ley de Fourier
sobre conducción del calor [9] establece que:

q = −k∇u, (1.1)

siendo k la constante de conductividad térmica del material. La aparición del signo

menos en la igualdad anterior se debe a que el calor fluye desde las regiones más
calientes a las más frı́as (k > 0) y el vector gradiente señala la dirección de máximo

9
 Figura 1.1: Volumen de control para deducir la ecuación no estática de difusión

crecimiento. Por lo tanto q señala a la dirección de máximo crecimiento y q es la

razón del flujo de calor en esa dirección.
Se deduce de lo anterior que si v es cualquier vector, entonces q · v es el flujo de calor
en la dirección de v .
Durante un breve intervalo de tiempo (t, t + Δt) el calor fluye a través del material y
puede también generarse por fuentes internas f (x, t). Por consiguiente, la cantidad de
calor que entra en cualquier región V del material en el intervalo (t, t + Δt) viene dada
por:

   
Q= − q · n dS + f (x, t) dV Δt, (1.2)
∂V V

siendo n el vector normal exterior de la superficie ∂V . El signo menos de la integral

aparece porque q ·n dS es la densidad de flujo de calor hacia el exterior de la superficie.
Sin embargo, este calor Q tiene el efecto de elevar la temperatura en ut Δt. Es decir,
puede escribirse como:

 
Q= cρut dV Δt, (1.3)
V

10
donde c es la capacidad por unidad de masa y ρ es la densidad de masa del material.
Aplicando el teorema de la divergencia, la integral sobre el borde de la superficie se
transforma en:

 
q · n dS = ∇ · q dV. (1.4)
∂V V

Sustituyendo la expresión (1.4) en la ecuación (1.2) se obtiene:

   
Q= − ∇ · q dV + f (x, t) dV Δt. (1.5)
V V

Luego, las ecuaciones (1.3) y (1.5) son iguales:

     
cρut dV = − ∇ · q dV + f (x, t) dV . (1.6)
V V V

Aplicando la ley de conducción de calor de Fourier [9] en la ecuación (1.6) se obtiene:

cρut = ∇ · (k∇u) + f, (1.7)

que es la forma general de la ecuación del calor. En la mayorı́a de los problemas k es

independiente de x y la ecuación (1.7) se transforma en

ut = ∇ · (K∇u) + F, (1.8)

donde
k f
K= ρ, F = .
c cρ
k
El coeficiente c
ρ es una constante llamada coeficiente de difusividad termal del
material y ∇u es el operador gradiente. Las funciones F (x, t) y u (x, t) representarán
el término no homogéneo de la ecuación y la solución del problema respectivamente.
Luego, sustituyendo estas expresiones en la ecuación (1.8) se obtiene la ecuación no
estática de difusión:

11
 ∂u
(x, t) = ∇ · (K∇u(x, t)) + F (x, t). (1.9)
∂t
A la ecuación (1.9) se le deben imponer condiciones iniciales u(x, 0) = u0 y condi-
ciones de frontera para obtener un problema bien planteado matemáticamente y ase-
gurar la existencia, unicidad y estabilidad de su solución [10]. Existen varios tipos de
condiciones de frontera dependiendo del estado de la temperatura que se encuentra en
el borde de la región de integración [14], como por ejemplo:

Condiciones de Dirichlet: corresponden a tener una temperatura prefijada en el

contorno, en este caso:
u|V = u0 .

Condiciones de Neumann: imponen un flujo calorı́fico u1 a través de la superficie

V , esto es:
∂u
−K |V = u 1 .
∂n

Condiciones de Robin o Mixtas: corresponden con preestablecer el intercambio

de calor en la frontera de acuerdo con la ley de Newton y toman la forma:

∂u
K + s(u − u0 )|V = 0,
∂n

donde s es el coeficiente de transferencia de calor y u0 es la temperatura del

medio.

Las que se utilizarán en este trabajo serán las condiciones de borde de tipo Robin no
triviales por ser las más generales.

12
Capı́tulo 2

Esquema mimético tipo

Crank-Nicolson

2.1. Introducción
Las soluciones de la ecuación no estática de difusión generalmente no se pueden
obtener mediante métodos analı́ticos. Por tal motivo se recurre a la utilización de
métodos numéricos que aproximan de forma sistemática las soluciones de la mencionada
ecuación. En este capı́tulo se presenta una descripción del método numérico mimético
desarrollado por Castillo-Grone-Yasuda [1, 2] el cual será combinado con el esquema
de resolución de Crank-Nicolson [18] para resolver la ecuación no estática de difusión
en una y dos dimensiones.

2.2. Definición del método mimético

Un método numérico es mimético si produce las expresiones discretas de los ope-
 ∂u 
radores gradiente (∇), divergencia (∇·) y derivada direccional en el borde ∂ n
, las

13
cuales serán denotadas por G, D y B respectivamente y si satisface la identidad de
Green generalizada:
Dv, f Q + v, Gf P = Bv, f Q , (2.1)

donde P y Q son sus pesos y el producto interno es definido por:

x, yA = y t Ax ;

con x e y vectores en Rn y A ∈ Rn×n una matriz definida positiva.

2.3. Discretización mimética unidimensional de los

operadores gradiente, divergencia y derivada
direccional en el borde
La discretización unidimensional de los operadores mediante el método mimético
requiere de una malla que difiere sustancialmente en los bordes de aquellas mallas
empleadas por los esquemas de diferencias finitas tradicionales. La siguiente figura
representa la malla mimética unidimensional para un intervalo cerrado [x0 , xN ], en
nuestro caso [0, 1]. En ella los puntos xi = ih para i = 0, ..., N definen la posición
de los lados de los bloques representados por las lı́neas verticales, el tamaño de cada
1
bloque viene dado por h = N
y los centros de cada celda vienen dados por la expresión
(xi +xi+1 )
2
y son representados por los puntos negros internos.

Figura 2.1: Malla mimética unidimensional

14
 Como se puede observar en la figura 2.1 existen nodos superpuestos en los bordes
externos; esta es una caracterı́stica propia de la malla mimética. También están re-
presentadas las posiciones donde se evalúan la solución de la ecuación no estática de
difusión u, el gradiente G y la divergencia D respectivamente. La forma discreta del
operador gradiente G para un punto cualquiera xi , ubicado en el interior de la malla,
viene dada por la siguiente expresión:

ui+ 1 − ui− 1
(Gu)i = 2 2
, 1≤i≤N , (2.2)
h
y para los puntos extremos x0 y xN viene dada por:

8 3 1
(Gu)0 = − u0 + u 1 − u 3 , (2.3)
3h h 2 3h 2

8 3 1
(Gu)N = uN − uN − 1 + uN − 3 , (2.4)
3h h 2 3h 2

respectivamente.
Las expresiones anteriores se obtienen al realizar los desarrollos de Taylor de se-
gundo orden de la función u y pueden ser representadas de manera matricial como
sigue:
⎛ ⎞
u0
⎛ ⎞ ⎜ ⎟
⎜ ⎟
− 83 3 − 13 0 ··· 0 ⎜ u1 ⎟
⎜ ⎟ ⎜ 2 ⎟
⎜ ⎟ ⎜ ⎟
⎜ 0 −1 1 0 ··· 0 ⎟ ⎜ u3 ⎟
⎜ ⎟ ⎜ 2 ⎟
⎜ .. .. .. .. ⎟ ⎜ .. ⎟
Gu = h1 ⎜ . . . . ⎟ ·⎜ . ⎟
⎜ ⎟ ⎜ ⎟
⎜ ⎟ ⎜ ⎟
⎜ 0 ··· 0 −1 1 0 ⎟ ⎜ uN − 1 ⎟
⎝ ⎠ ⎜ 2 ⎟
⎜ ⎟
0 ··· 0 1
−3 8 ⎜ uN + 1 ⎟
3 3 (N +1)×(N +2) ⎝ 2 ⎠
uN +1 .
De manera similar, la forma discreta del operador divergencia D viene dada por:
vi+1 − vi
Di+ 1 v = , (2.5)
2 h
15
y su expresión matricial es:
⎛ ⎞ ⎛ ⎞
0 0 0 0 ··· 0 v
⎜ ⎟ ⎜ 0 ⎟
⎜ ⎟ ⎜ ⎟
⎜ −1 1 0 0 ··· 0 ⎟ ⎜ v1 ⎟
⎜ ⎟ ⎜ ⎟
⎜ ⎟ ⎜ ⎟
⎜ 0 −1 1 0 ··· 0 ⎟ ⎜ v2 ⎟
⎜ ⎟ ⎜ ⎟
⎜ . . ⎟ ⎜ . ⎟
Dv = h1 ⎜ .. . . . . . . . . . . . . .. ⎟ · ⎜ .. ⎟
⎜ ⎟ ⎜ ⎟
⎜ ⎟ ⎜ ⎟
⎜ 0 · · · 0 −1 1 0 ⎟ ⎜ vN −2 ⎟
⎜ ⎟ ⎜ ⎟
⎜ ⎟ ⎜ ⎟
⎜ 0 ··· 0 0 −1 1 ⎟ ⎜ vN −1 ⎟
⎝ ⎠ ⎝ ⎠
0 ··· 0 0 0 0 vN .
(N +2)×(N +1)

Obsérvese que la discretización del operador divergencia D es mucho más sencilla

que la del operador gradiente G. Este hecho se debe a que la discretización de D
está asociada sólo a los nodos que se encuentran en el centro de cada bloque. Nótese
también que la matriz D se ha completado con dos filas de ceros al principio y al final
para evitar problemas de dimensiones al componer los operadores D y G.

2.4. Discretización mimética bidimensional de los

operadores gradiente, divergencia y derivada
direccional en el borde
Al igual que el caso unidimensional, para hacer la discretización de los operadores
anteriormente mencionados, se necesita de la malla mimética bidimensional la cual se
obtiene haciendo el producto cartesiano de las mallas miméticas unidimensionales y
quitando los puntos correspondientes a las esquinas. En la figura que sigue se muestra
una malla de orden 5 × 5.

16
 Figura 2.2: Malla mimética bidimensional de orden 5 × 5

En la malla de la izquierda de la figura 2.2 se definen los operadores D, B y la

solución u, mientras que en la malla de la derecha de la misma figura se muestran los
puntos donde se evalúan las componentes Gx y Gy del operador G.
Los resultados obtenidos en el caso unidimensional son fundamentales para el estu-
dio del método en dos dimensiones. Las discretizaciones miméticas de los operadores
gradiente G y divergencia D para el caso bidimensional son simples extensiones uni-
dimensionales. En efecto, la expresión analı́tica para el operador divergencia D viene
dada por:
ui+1,j − ui,j vi,j+1 − vi,j
D(u, v)i,j = + , (2.6)
h h
y su estructura matricial se muestra en la siguiente figura:

17
 Figura 2.3: Estructura matricial del operador Divergencia en 2D

La estructura de la matriz asociada al operador gradiente se puede apreciar en la

siguiente figura:

Figura 2.4: Matriz del operador gradiente en 2D

Mientras que la estructura de la matriz asociada al operador derivada direccional

en el borde B es la siguiente :

18
 Figura 2.5: Matriz del operador derivada direccional en el borde B en 2D

donde cada recuadro de la figura 2.5 representa las siguientes secciones de sub-
matrices de B:
⎛ ⎞⎛ ⎞⎛ ⎞⎛ ⎞
−1 −1 − 18 1
1
⎜ ⎟⎜ ⎟⎜ 8 ⎟⎜ ⎟
⎜ .. ⎟⎜ 1 ⎟⎜ ⎟⎜ .. ⎟
⎜ . ⎟⎜ − 18 ⎟ ⎜ 1
− 18 ⎟ ⎜ . ⎟
⎝ ⎠⎝ 8 ⎠⎝ 8 ⎠⎝ ⎠
−1 − 18 1
8
1 1

Hay que destacar que estas matrices dependen fuertemente de la manera en que
se enumeran los nodos de la malla bidimensional. En nuestro caso en particular se
consideró la numeración natural que se muestra en la figura 2.2.
Después de haber mostrado las versiones miméticas de todos los operadores, se puede
generar la expresión discretizada de la ecuación no estática de difusión (1.9) como
sigue:
∂u
− ∇ · (K∇u) ≈ T u − D (KG) u, (2.7)
∂t
y de las condiciones de borde:

∂u
αu + βK ≈ α(A)u + βK (BG) u, (2.8)
∂n
donde K es una matriz cuadrada y diagonal que representa los coeficientes de tran-
sferencia de calor, T es la matriz del tiempo y α y β son escalares distintos de cero.

19
2.5. Resolución de la ecuación no estática de di-
fusión por el esquema de Crank-Nicolson
Para encontrar una discretización completa de la ecuación no estática de difusión
se necesita discretizar también la variable temporal. Para ello existen varios esquemas
entre los cuales se destacan: el esquema explı́cito y el esquema implı́cito [18]. El primero
permite que las soluciones en el paso n + 12 puedan obtenerse a partir de las soluciones
en el paso n utilizando esta formulación:
n+ 12
uj − unj unj+1 − 2unj + unj−1
Δt
= + Fjn . (2.9)
2
(Δx)2
Δt
A la cantidad r = (Δx)2
se le llama número de Courant. La molécula computacional,
que indica la ubicación de los nodos de la malla en el espacio y en el tiempo para el
esquema explı́cito, se muestra en la siguiente figura:

Figura 2.6: Molécula computacional del esquema explı́cito

En el esquema implı́cito, en cambio, la formulación analı́tica para obtener la dis-

cretización del tiempo viene dada por la expresión siguiente:
n+ 21
un+1
j − uj un+1 n+1
j+1 − 2uj + un+1
j−1
Δt
= 2
+ Fjn+1 . (2.10)
2
(Δx)

20
 n+ 12
Nótese que el cálculo de los valores de un+1
j a partir de los ui no se puede realizar
directamente sino resolviendo primero un sistema de ecuaciones lineales. Su molécula
computacional es la representada en la figura que sigue:

Figura 2.7: Molécula computacional del esquema implı́cito

Otra posibilidad para aproximar la derivada temporal con mayor precisión es la

combinación de los esquemas (2.9) y (2.10) obteniéndose la fórmula analı́tica del es-
quema de Crank-Nicolson:

  
un+1
j − unj n+1 n+1 n+1
1 uj+1 − 2uj + uj−1 1 unj+1 − 2unj + unj−1 1  n+1 n

= + + Fj + Fj .
Δt 2 (Δx)2 2 (Δx)2 2
(2.11)
La molécula computacional para el caso unidimensional se muestra en la siguiente
figura:

21
 Figura 2.8: Molécula computacional del esquema de Crank-Nicolson 1D

Mientras que en la proxima figura se puede apreciar la molécula computacional

para el caso bidimensional.

Figura 2.9: Molécula computacional del esquema de Crank-Nicolson 2D

2.6. Esquema mimético tipo Crank-Nicolson

Después de haber presentado las descripciones del método mimético desarrollado
por Castillo-Grone-Yasuda [1, 2] y del esquema de resolución de Crank-Nicolson [18],

22
se propone la combinación de ambos esquemas para resolver la ecuación no estática de
difusión en una y dos dimensiones; este método será llamado Esquema Mimético tipo
Crank-Nicolson. Al aplicar el esquema (2.11) a las ecuaciones (2.7) y (2.8) se obtiene
la aproximación mimética de la ecuación no estática de difusión:
1 1 1
[T ](un+1 − un ) = [D](KG)] un + [[D](KG)] un+1 + (F n+1 + F n ) , (2.12)
2 2 2
y la discretización mimética de la condición de borde tipo Robin:
1 1 1
[−α[A] − β[BG]]un + [−α[A] − β[BG]]un+1 + (f n+1 + f n ) = 0 , (2.13)
2 2 2
donde f representa el término no homogéneo. Utilizando el principio de superposición
se acoplan las expresiones (2.12) y (2.13) formándose el siguiente sistema:
1 1
[T ](un+1 −un ) = [D(KG)−α[A]−β[BG]]un + [D(KG)−α[A]−β[BG]]un+1 + (2.14)
2 2

1 1
+ (F n+1 + F n ) + (f n+1 + f n ).
2 2
La ecuación (2.14) se transforma en la expresión:
1 1
[T − ([DKG] − α[A] − β[BG])]un+1 = [T + ([DKG] − α[A] − β[BG])]un + (2.15)
2 2
1 1
+ (F n+1 + F n ) + (f n+1 + f n ),
2 2
la cual representa un sistema lineal cuya resolución numérica produce soluciones apro-
ximadas para la ecuación no estática de difusión.

2.7. Ecuación no estática de difusión bidimensional

discretizada
La ecuación no estática de difusión en dos dimensiones posee la siguiente forma:
 2 
∂u ∂ u ∂2u
− + = F (x, y, t) (2.16)
∂t ∂x2 ∂y 2

23
sujeta a las condiciones de Robin siguientes:
 
∂u
α1 u − (x, 0, t) = f (x, 0, t), (2.17)
∂y
 
∂u
α1 u + (x, 1, t) = f (x, 1, t), (2.18)
∂y
 
∂u
α2 u + (1, y, t) = f (1, y, t), (2.19)
∂x
 
∂u
α2 u − (0, y, t) = f (0, y, t), (2.20)
∂x
y condición inicial:
u(x, y, 0) = u0 (x, y). (2.21)

Los coeficientes α1 , α2 son escalares distintos de cero. A continuación se planteará la

discretización para el problema de contorno (2.16)-(2.21).

2.7.1. Discretización mimética de las ecuaciones bidimensionales

En esta sección se presentarán las ecuaciones discretizadas que conforman el sistema

lineal asociado a la numeración mostrada en la figura 2.2.

Figura 2.10: Nodos de la malla mimética bidimensional de orden 5 × 5

24
 La figura 2.10 representa una malla bidimensional correspondiente a la región Ω =
[0, 1] × [0, 1] y será tomada como referencia a lo largo de esta sección.
La primera ecuación viene dada por la expresión:
   
1 4 n+1 3 n+1 1 n+1 1 4
α+ um i+ 1 ,j− 1 − um i+ 1 ,j+ 1 + um i+ 1 ,j+ 3 = − α − umn(i+ 1 ,j− 1 )
2 3h ( 2 2 ) 2h ( 2 2 ) 6h ( 2 2 ) 2 3h 2 2

 
3 n 1 n 1 n+1 n
+ um(i+ 1 ,j+ 1 ) − um(i+ 1 ,j+ 3 ) + f + f(i+ 1 ,j )
2h 2 2 6h 2 2 2 (i+ 12 ,j ) 2

y corresponde a la discretización de la condición de borde tipo Robin en cada uno de

los lados del cuadrado unitario. Los tres términos de esta ecuación denotan los efectos
del operador gradiente sobre el borde. El vector um
 representa la aproximación a la
solución u. En la figura que sigue se muestra la molécula computacional centrada en el
nodo N1=(xi+ 1 , yj− 1 ):
2 2

Figura 2.11: Molécula computacional correspondiente al nodo N1

Esta ecuación:
 
3 9 1 3 3
− 2 umn+1 + + umn+1 1 − umn+1
− umn+1
2h (i,j+ 21 ) 2h2 Δt 1
(i+ 2 ,j+ 2 ) 4h2 1 3
(i+ 2 ,j+ 2 ) 2h2 (i+ 12 ,j )

 
3 n+1 3 n 9 1 3
− 2 um i+ 3 ,j+ 1 = 2 um(i,j+ 1 ) + − 2 + umn(i+ 1 ,j+ 1 ) + 2 umn(i+ 1 ,j+ 3 ) +
4h ( 2 2 ) 2h 2 2h Δt 2 2 4h 2 2

25
  
3 n 3 n 1 n+1 n
+ 2 um(i+ 1 ,j ) + 2 um(i+ 3 ,j+ 1 ) + F i+ 1 ,j+ 1 + F(i+ 1 ,j+ 1 )
2h 2 4h 2 2 2 ( 2 2) 2 2
 
representa la discretización del nodo N2= xi+ 1 , yj+ 1 el cual puede ser visualizado en
2 2

la siguiente figura:

Figura 2.12: Molécula computacional correspondiente al nodo N2

La ecuación que sigue representa la discretización mimética del nodo N3 de coor-

 
denadas xi+ 3 , yj+ 1 :
2 2

 
1 3 10 1 1
umn+1 1 − umn+1 + + umn+1 − umn+1
6h 2 ( i,j+ 2 ) 4h 2 (i+ 32 ,j+ 23 ) 3h2 Δt (i+ 32 ,j+ 21 ) 2h2 (i+ 52 ,j+ 21 )

3 3 1 3
− 2
umn+1 3 − 2
umn+1
1 1 = − 2
umn n
( 2 ) 4h2 um(i+ 23 ,j+ 23 ) +
i,j+ 1 +
2h (i+ 2 ,j ) 4h (i+ 2 ,j+ 2 ) 6h

 
10 1 1 3 3
+ − 2+ umn(i+ 3 ,j+ 1 ) + 2 umn(i+ 5 ,j+ 1 ) + 2 umn(i+ 3 ,j ) + 2 umn+1 +
3h Δt 2 2 2h 2 2 2h 2 4h (i+ 12 ,j+ 21 )

 
1 n+1 n
+ F i+ 3 ,j+ 1 + F(i+ 3 ,j+ 1 )
2 ( 2 2) 2 2

26
el cual es mostrado en la siguiente figura:

Figura 2.13: Molécula computacional correspondiente al nodo N3

1
En la ecuación, el coeficiente que posee el término Δt está asociado al nodo N3 de
 
coordenadas xi+ 3 , yj+ 1 , mientras que los otros términos se refieren a los nodos que
2 2

se encuentran alrededor del mismo incluyendo el nodo que se encuentra en el borde

izquierdo.

Siguiendo la numeración de los nodos mostrados en la figura 2.10, la próxima

ecuación:
 
1 13 1 1 3
− 2 umn+1 + + umn+1 1 − umn+1
− umn+1
2h (i+ 72 ,j+ 21 ) 4h2 Δt 5
(i+ 2 ,j+ 2 ) 2h2 3 1
(i+ 2 ,j+ 2 ) 2h2 (i+ 52 ,j )

 
3 n+1 1 n 13 1 1
− 2 um i+ 5 ,j+ 3 = 2 um(i+ 7 ,j+ 1 ) + − 2 + umn(i+ 5 ,j+ 1 ) + 2 umn(i+ 3 ,j+ 1 ) +
4h ( 2 2 ) 2h 2 2 4h Δt 2 2 2h 2 2

 
3 n 3 n 1 n+1 n
+ 2 um i+ 5 ,j + 2 um i+ 5 ,j+ 3 + F i+ 5 ,j+ 1 + F i+ 5 ,j+ 1
2h ( 2 ) 4h ( 2 2) 2 ( 2 2) ( 2 2)
 
le correspondiente al nodo N4= xi+ 5 , yj+ 1 el cual se repite N − 4 veces dependiendo
2 2

del tamaño de la malla y se puede apreciar en la figura que sigue:

27
 Figura 2.14: Molécula computacional correspondiente al nodo N4

Esta ecuación:
 
1 3 13 1 1
− 2 umn+1 3 − umn+1 + + umn+1 − umn+1 +
6h (i,j+ 2 ) 4h 2 (i+ 21 ,j+ 23 ) 6h2 Δt (i+ 32 ,j+ 23 ) 2h2 (i+ 52 ,j+ 23 )

1 3 1 1 3
+ 2
umn+1
i+ 3
,j
− 2 umn+1
i+ 3
,j+ 1 − 2
umn+1
i+ 3
,j+ 5 = 2 umn(i,j+ 3 ) + 2 umn(i+ 1 ,j+ 3 ) +
6h ( 2 ) 4h ( 2 2 ) 2h ( 2 2 ) 6h 2 4h 2 2

 
13 1 1 1 3
+ − 2+ umn(i+ 3 ,j+ 3 ) + 2 umn(i+ 5 ,j+ 3 ) − 2 umn(i+ 3 ,j ) + 2 umn(i+ 3 ,j+ 1 ) +
6h Δt 2 2 2h 2 2 6h 2 4h 2 2

 
1 1
+ 2 umn(i+ 3 ,j+ 5 ) + F n+1 + F(ni+ 3 ,j+ 3 )
2h 2 2 2 (i+ 32 ,j+ 23 ) 2 2

es la discretización correspondiente al nodo N5 cuya coordenadas cartesianas son:

 
xi+ 3 , yj+ 3 y es obtenida utilizando el esquema mimético tipo Crank-Nicolson de-
2 2

scrito en la sección anterior. Nótese que esta formulación es diferente a la que se ob-
tiene con el método de diferencias finitas tradicional pues en ella intervienen los nodos
ubicados en los bordes izquierdo e inferior sin hacer uso de puntos auxiliares. La figura
siguiente ilustra esta situación:

28
 Figura 2.15: Molécula computacional correspondiente al nodo N5

En la próxima ecuación:

 
1 1 25 1 1
− 2 umn+1 7 3 − umn+1 + + umn+1 3 + umn+1
2h (i+ 2 ,j+ 2 ) 2h2 (i+ 32 ,j+ 23 ) 12h 2 Δt 5
(i+ 2 ,j+ 2 ) 6h2 (i+ 52 ,j )

3 1 1 1
− 2
umn+1
i+ 5
,j+ 1 + 2
umn+1
i+ 5
,j+ 5 = 2 umn(i+ 7 ,j+ 3 ) + 2 umn(i+ 3 ,j+ 3 ) +
4h ( 2 2 ) 6h ( 2 2 ) 2h 2 2 2h 2 2

 
25 1 1 3
+ − 2
+ umn(i+ 5 ,j+ 3 ) − 2 umn(i+ 5 ,j ) + 2 umn(i+ 5 ,j+ 1 ) −
12h Δt 2 2 6h 2 4h 2 2

 
1 1
− 2 umn(i+ 5 ,j+ 5 ) + F n+1 + F(ni+ 5 ,j+ 3 )
6h 2 2 2 (i+ 52 ,j+ 23 ) 2 2
 
se describe la discretización de los N − 4 nodos centrales N6= xi+ 5 , yj+ 3 ubicados
2 2

en la tercera fila como se ve en la siguiente figura:

29
 Figura 2.16: Molécula computacional correspondiente al nodo N6

Obsérvese que al igual que el nodo N5 también ésta no es la tı́pica ecuación que se
obtiene con el método de diferencias finitas ya que en ella se incluyen los nodos de los
bordes izquierdo e inferior sin usar puntos imaginarios.

Finalmente la ecuación:

 
1 2 1 1 1
− 2 umn+1 + + umn+1 5 − umn+1
− umn+1
2h (i+ 52 ,j+ 23 ) h2 Δt 5
(i+ 2 ,j+ 2 ) 2h2 3 5
(i+ 2 ,j+ 2 ) 2h2 (i+ 72 ,j+ 23 )

 
1 n+1 1 n 2 1 1
− 2 um i+ 5 ,j+ 7 = 2 um(i+ 5 ,j+ 3 ) + − 2 + umn(i+ 5 ,j+ 5 ) + 2 umn(i+ 3 ,j+ 5 ) +
2h ( 2 2 ) 2h 2 2 h Δt 2 2 2h 2 2

 
1 1 1
+ 2 umn(i+ 7 ,j+ 3 ) + 2 umn(i+ 5 ,j+ 7 ) + F n+1 + F(ni+ 5 ,j+ 5 )
2h 2 2 2h 2 2 2 (i+ 52 ,j+ 25 ) 2 2

representa la discretización tı́pica de diferencias finitas tradicionales de los nodos cen-

 
trales N7 de coordenadas xi+ 5 , yj+ 5 mostrados en la siguiente figura:
2 2

30
 Figura 2.17: Molécula computacional correspondiente al nodo N7

Una vez aplicado el esquema mimético de Crank-Nicolson a todos los nodos siguien-
do la numeración preestablecida en la figura 2.2, se obtiene este sistema de ecuaciones
lineales que al ser resuelto numéricamente produce soluciones aproximadas para la
ecuación no estática de difusión en dos dimensiones. En la siguiente figura se muestra
la estructura de dicho sistema:

Figura 2.18: Estructura del sistema bidimensional

Se puede observar que la estructura de la matriz asociada al sistema bidimensional

es esparcida y estructuralmente simétrica respecto a su diagonal principal.

31
Capı́tulo 3

Análisis de convergencia

3.1. Introducción
En este capı́tulo se expondrá un análisis de convergencia para la ecuación no estática
de difusión. La presentación de dicho resultado teórico tiene el siguiente orden: inicial-
mente se describirán los conceptos de consistencia, estabilidad y convergencia de un
método numérico, luego, se darán a conocer las ecuaciones discretizadas obtenidas por
medio de la aplicación del esquema mimético tipo Crank-Nicolson en una dimensión,
también se presentarán las pruebas de consitencia y estabilidad del nuevo esquema
solamente para el caso de problemas unidimensionales y finalmente haciendo uso el
teorema de Lax [10] se demostrará la convergencia del esquema propuesto.

3.2. Preliminares
A continuación se presentan las definiciones de unos conceptos básicos necesarios
para hacer el análisis de convergencia del método que se está proponiendo.

Definición 3.1 Si el esquema numérico se asemeja a la ecuación en derivadas par-

32
ciales cuando los parámetros h y Δt tienden a cero, entonces se dice que el esquema
es consistente.

Definición 3.2 Si la solución de la ecuación diferencial discreta permanece acotada

para cualquier intervalo de tiempo (t, t + Δt), entonces el esquema es estable.

Definición 3.3 Si en un punto (x, t) cualquiera del dominio Ω se cumple que la solu-
ción aproximada en ese punto tiende a la solución exacta cuando h y Δt tienden a
cero, entonces se dice que el esquema numérico es convergente.

Estos tres conceptos [19] juegan un rol importantı́simo en la teorı́a de las aproxima-
ciones en diferencias finitas y están relacionados entre sı́ a través del siguiente:

Teorema 3.1 (Lax) Dado un problema de valores iniciales bien planteado matemática-
mente y una aproximación en diferencias finitas a él que sea consistente y estable,
entonces es convergente.

Este teorema es de gran importancia pues relaciona el objetivo final de toda apro-
ximación en diferencias finitas, es decir, la convergencia a la solución exacta, con una
propiedad que es mucho más conveniente para demostrar: la estabilidad.

3.3. Convergencia del esquema mimético

Se discretizará la ecuación no estática de difusión en una dimensión

∂u ∂ 2 u
− 2 = F (x, t) (3.1)
∂t ∂x

sujeta a las condiciones de Robin siguientes:

∂u
α1 u(0) − (0) = f (0) en x=0 (3.2)
∂x

33
 ∂u
α2 u(1) − (1) = f (1) en x=1
∂x
y condición inicial:
u(x, 0) = u(0) (3.3)

Los coeficientes α1 , α2 son escalares distintos de cero.

3.3.1. Discretización mimética de las ecuaciones unidimen-

sionales

Para hacer el estudio de consistencia y estabilidad del método, se necesita presen-

tar las discretizaciones de las ecuaciones asociadas a los nodos de la malla mimética
uniforme mostrada en la figura 2.1 obtenidas de la aplicación del esquema mimético
de Crank-Nicolson (2.15). La aproximación mimética de la solución u es dada por el
vector:
 t
 = um0 , um 1 , . . . , umN − 1 , umN .
um
2 2

La primera ecuación viene dada por la siguiente expresión:

   
1 4 n+1 3 n+1 1 n+1 1 4
α1 + um0 − um 1 + um 3 = − α1 − umn0 + (3.4)
2 3h 2h 2 6h 2 2 3h

3 1 1  n+1 
umn1 − umn3 + f0 + f0n
2h 2 6h 2 2

y representa la discretización de segundo orden de la condición de borde tipo Robin

evaluada en el extremo izquierdo x0 del intervalo que se muestra en la siguiente figura:

34
 Figura 3.1: Nodo ubicado en el extremo izquierdo.

La segunda ecuación del sistema es:

     
1 4 n+1 1 2 1 n+1 1 2
− − 2 um0 + + + um 1 + − − umn+1
3 = (3.5)
6 3h 4h h2 Δt 2 12h 3h2 2

     
1 4 1 2 1 1 2 1  n+1 
+ 2 umn0 + − − 2+ n
um 1 + + 2 n
um 3 + n
F1 + F1
6 3h 4h h Δt 2 12h 3h 2 2 2 2

y representa la discretización de la ecuación no estática de difusión centrada en el

nodo x 1 el cual se puede observar en la siguiente figura:

2

Figura 3.2: Nodo correspondiente a la segunda ecuación.

La tercera ecuación viene dada por:

   
1 n+1 1 1 n+1 1 1 1 1
um0 + − − 2 um 1 + + 2+ umn+1
3 − 2 umn+1
5 =
6h 4h 2h 2 12h h Δt 2 2h 2

   
1 1 1 1 1 1
− umn0 + + 2 n
um 1 + − − 2+ umn3 + (3.6)
6h 4h 2h 2 12h h Δt 2

1 n 1  n+1 n

+ um 5 + F 3 + F 3
2h2 2 2 2 2

y es la discretización centrada en el nodo x 3 el cual se muestra en la siguiente figura:

2

35
 Figura 3.3: Nodo correspondiente a la tercera ecuación.

Las ecuaciones siguientes:

 
1 n+1 1 1 1 1
− 2 umi− 1 + 2
+ umn+1
i+ 1 −
2
umn+1
i+ 3 =
2
umni− 1 + (3.7)
2h 2 h Δt 2 2h 2 2h 2

 
1 1 1 1  n+1 
+ − 2+ umni+ 1 + 2 umni+ 3 + n
Fi+ 1 + Fi+ 1
h Δt 2 2h 2 2 2 2

corresponden a la discretización estándar de los nodos centrales que se pueden apreciar

en la próxima figura:

Figura 3.4: Nodos correspondientes a las ecuaciones estándard.

La antepenúltima ecuación:
   
1 n+1 1 1 1 1 1 1 n+1
− 2 uN − 5 + + 2+ uN − 3 + − − 2 un+1
n+1
1 + u =
2h 2 12h h Δt 2 4h 2h N− 2 6h N

   
1 n 1 1 1 n 1 1
u 5+ − − + uN − 3 + + unN − 1 − (3.8)
2h2 N − 2 12h h2 Δt 2 4h 2h2 2

36
 1 n 1  n+1 n

− uN + F 3 + FN − 3
6h 2 N− 2 2

es la discretización centrada en el nodo xN − 3 mostrado en la figura:

2

Figura 3.5: Nodo correspondiente a la antepenúltima ecuación.

La penúltima ecuación es dada por:

     
1 2 n+1 1 2 1 1 4
− − umN − 3 + + + umN − 1 + − − 2 umn+1
n+1
N = (3.9)
12h 3h2 2 4h h2 Δt 2 6h 3h

     
1 2 1 2 1 1 4 1  n+1 
= + 2 umnN − 3 + − − 2+ n
umN − 1 + + n n
umN + FN − 1 + FN − 1
12h 3h 2 4h h Δt 2 6h 3h2 2 2 2

y representa la discretización de la ecuación no estática de difusión centrada en el

nodo xN − 1 que se muestra en la próxima figura:

2

Figura 3.6: Nodo correspondiente a la penúltima ecuación.

Finalmente la ecuación:
   
1 4 n+1 3 n+1 1 n+1 1 4
α2 + umN − umN − 1 + umN − 3 = − α2 − umnN + (3.10)
2 3h 2h 2 6h 2 2 3h

37
 3 1 1  n+1 
+ umnN − 1 − umnN − 3 + f1 + f1n
2h 2 6h 2 2
aproxima la condición de borde tipo Robin en el extremo derecho del intervalo y
está centrada en el nodo xN que se muestra a continuación:

Figura 3.7: Nodo correspondiente a la última ecuación.

La aplicación del esquema mimético tipo Crank-Nicolson a los nodos de la malla

unidimensional produce un sistema de ecuaciones lineales cuya estructura se muestra
a continuación:

Figura 3.8: Estructura del sistema unidimensional

Obsérvese que la matriz asociada a este sistema es casi tridiagonal lo cual hace que
el costo computacional para su resolución sea menor.

38
3.3.2. Consistencia del esquema mimético tipo Crank-Nicolson

El primer paso de la prueba consiste en la sustitución de las condiciones de borde

dadas por (3.4) y (3.10) en las ecuaciones (3.5), (3.6) y (3.8), (3.9) respectivamente. A
partir de la ecuación:
 
1 4  n+1  3 1 3
α1 + um0 + umn0 − umn+1
1 + umn+1
3 − umn1 + (3.11)
2 3h 2h 2 6h 2 2h 2

1 1  n+1 
+ umn3 = f0 + f0n
6h 2 2
sea
 
1 4   3 1 3 1
α1 + umn+1
0 + umn0 − umn+1
1 + umn+1
3 − umn1 + umn3 = γ0 (3.12)
2 3h 2h 2 6h 2 2h 2 6h 2

 
Despejando el término umn+1
0 + umn0 de la ecuación (3.12) se obtiene:
 
 n+1 n
 6h 3 n+1 1 n+1
um0 + um0 = γ0 + um 1 − um 3 (3.13)
3hα1 + 8 2h 2 6h 2

 
6h 3 n 1 n
+ um 1 − um 3
3hα1 + 8 2h 2 6h 2

Sustituyendo la expresión (3.13) en las ecuaciones (3.5) y (3.6) y simplificándolas se

obtiene por un lado:
   
−3(h2 + 8) h+8 1 n+1 h2 + 8 h+8
2
+ 2
+ um 1 + 2
− umn+1
3
2h (3α1 + 8) 4h Δt 2 6h (3hα1 + 8) 12h2 2

   
−3(h2 + 8) h+8 1 n h2 + 8 h+8
+ + − um 1 + − umn3 − (3.14)
2h2 (3hα1 + 8) 4h2 Δt 2 6h2 (3hα1 + 8) 12h2 2

h2 + 8 1  n+1 
− γ0 = F 1 + F 1n
h(3hα1 + 8) 2 2 2

39
 y por otro:
   
3 h+2 n+1 −1 h + 12 1 1
− um 1 + + + umn+1
3 − umn+1
5 +
2h(3hα1 + 8) 4h2 2 6h(3hα1 + 8) 12h2 Δt 2 2h2 2

   
3 h+2 n −1 h + 12 1
+ − um 1 + + − umn3 − (3.15)
2h(3hα1 + 8) 4h2 2 6h(3hα1 + 8) 12h2 Δt 2

1 n γ0 1  n+1 n

− 2 um 5 + = F3 + F3
2h 2 3hα1 + 8 2 2 2

Obsérvese que en ambas ecuaciones no aparecen los términos umn+1

0 y umn0 respectivos
al borde izquierdo.

Las ecuaciones correspondientes a los nodos centrales permanecen intactas ya que

en ellas no intervienen los factores um0 y umN , eso es:

 
1 1 1 1 1
− 2 umn+1
i− 21
+ 2
+ umn+1
i+ 1
− 2
umn+13 = umni− 1 + (3.16)
2h h Δt 2 2h i+ 2 2h2 2

 
1 1 n 1 n 1  n+1 n

+ − 2+ umi+ 1 + 2 umi+ 3 + F 1 + Fi+ 1
h Δt 2 2h 2 2 i+ 2 2

De manera similar, a partir de la ecuación:

 
1 4  n+1  3 1 3
α2 + umN + umnN − umn+1 N − 1 + umn+1
N − 3 − umnN − 1 + (3.17)
2 3h 2h 2 6h 2 2h 2

1 1  n+1 
+ umnN − 3 = FN + FNn
6h 2 2

sea
 
1 4  n+1  3 1 3
α2 + umN + umnN − umn+1
N − 1 + umn+1
N − 3 − umnN − 1 + (3.18)
2 3h 2h 2 6h 2 2h 2

40
 1
umnN − 3 = γ1
+
6h 2

 
Despejando el término umn+1
N + umnN de la ecuación (3.18) se obtiene:

 
  6h 3 1
umn+1
N + umnN = n+1 n+1
γ1 + umN − 1 − umN − 3 + (3.19)
3hα2 + 8 2h 2 6h 2

 
6h 3 1
+ umnN − 1 − umnN − 3
3hα2 + 8 2h 2 6h 2

Al sustituir la expresión (3.19) en la ecuaciones (3.8) y (3.9) se obtiene por un lado:

   
3 h+2 n+1 1 h + 12 1
− umN − 1 + − + + umn+1
N − 32
−
2h(3hα2 + 8) 4h2 2 6h(3hα2 + 8) 12h2 Δt

 
1 3 h+2
− 2 umn+1
N − 52
+ − umnN − 1 + (3.20)
2h 2h(3hα2 + 8) 4h2 2

 
−1 h + 12 1 1 γ1 1  n+1 
+ + 2
− umnN − 3 − umn
5+
N− 2
= F 3 + F n
3
N− 2
6h(3hα2 + 8) 12h Δt 2 2h2 3hα2 + 8 2 N− 2

y por otro:

   
−3(h2 + 8) h+8 1 h2 + 8 h+8
2
+ 2
+ umn+1
N − 12
+ 2
− umn+1
N− 3
+
2h (3hα2 + 8) 4h Δt 6h (3hα2 + 8) 12h2 2

   
−3(h2 + 8) h+8 1 h2 + 8 h+8
+ + − umnN − 1 + − umnN − 3 −
2h2 (3α2 + 8) 4h2 Δt 2 6h2 (3hα2 + 8) 12h2 2

(3.21)

41
 h2 + 8 1  n+1 
− γ1 = FN − 1 + FNn − 1
h(3hα2 + 8) 2 2 2

También en estas ecuaciones no aparecen los términos umn+1

N y umnN correspon-
dientes al borde derecho.
Las ecuaciones (3.14), (3.15), (3.20) y (3.21) se pueden reescribir como sigue:

     
1 4 n+1 1 2 1 n+1 1 2
− − 2 um0 + + + um 1 + − − umn+1
3 +
6 3h 4h h2 Δt 2 12h 3h2 2

     
1 4 n 1 2 1 n 1 2
+ − − 2 um0 − − − 2 + um 1 − + umn3 − +
6 3h 4h h Δt 2 12h 3h2 2

  
1  h2 + 8 1 4 3 1
− F 1n+1 + F 1n + α1 + umn+1
0
n+1 n+1
− um 1 + um 3 +
2 2 2 h(3hα1 + 8) 2 3h 2h 2 6h 2

  
h2 + 8 1 4 n 3 n 1 n
+ − α1 − um0 − um 1 + um 3 − γ0 = 0
h(3hα1 + 8) 2 3h 2h 2 6h 2

Esta expresión contiene la ecuación correspondiente al nodo x 1 y la aproximación

2

de la condición de borde.

   
1  n+1 n
 1 1 n+1 1 1 1
um0 + um0 + − − 2 um 1 + + + umn+1
3 −
6h 4h 2h 2 12h h2 Δt 2

   
1 1 1 1 1 1 1
− 2 umn+1
5 − + 2 n
um 1 − − − 2+ umn3 − umn5 −
2h 2 4h 2h 2 12h h Δt 2 2h2 2

42
   
1  1 1 4 3 1
− F 3n+1 + F 3n − α1 + umn+1
0
n+1 n+1
− um 1 + um 3 −
2 2 2 3hα1 + 8 2 3h 2h 2 6h 2

   
1 1 4 n 3 n 1 n
− − − α1 − um0 − f 1 + um 3 − γ0 = 0
3hα1 + 8 2 3h 2h 2 6h 2

De manera similar, esta ecuación contiene la expresión algebraica correspondiente

al nodo x 3 y la aproximación de la condición de borde tipo Robin.
2

   
1 n+1 1 1 n+1 1 n+1 1 1 1
umN + − − 2 umN − 1 − 2 umN − 5 + + 2+ umn+1
N − 32
+
6h 4h 2h 2 2h 2 12h h Δt

   
1 n 1 1 1 n 1 1
− 2 umN − 5 − − − 2− umN − 3 − + 2 umnN − 1 −
2h 2 12h h Δt 2 4h 2h 2

  
1  1 1 4 3 1
− FNn+1
− 3 + FN
n
− 3 − α2 + n+1 n+1 n+1
umN − umN − 1 + umN − 3 −
2 2 2 3hα2 + 8 2 3h 2h 2 6h 2

  
1 1 4 3 1
− α2 + umnN n n
− umN − 1 + umN − 3 − γ1 = 0
3hα2 + 8 2 3h 2h 2 6h 2

Esta ecuación está relacionada con el nodo xN − 3 . En ella se puede observar la

2

ecuación correspondiente a dicho nodo y la condición de borde.

     
1 4  n+1  1 2 n+1 1 2 1
− − 2 n
umN + uN + − − umN − 3 + + + umn+1
N − 21
−
6 3h 12h 3h2 2 4h h2 Δt

43
    
1 2 1 2 1 1  n+1 
− + 2 umn+1
N − 23
− − − 2+ umnN − 1 − FN − 1 + FNn − 1 +
12h 3h 4h h Δt 2 2 2 2

  
h2 + 8 1 4 3 n+1 1
+ α2 + umn+1
N
n+1
− fN − 1 + umN − 3 +
h(3hα2 + 8) 2 3h 2h 2 6h 2

  
h2 + 8 1 4 3 1
+ α2 + umnN n n
− umN − 1 + umN − 3 − γ1 = 0
h(3hα2 + 8) 2 3h 2h 2 6h 2

Esta ecuación es la que corresponde al nodo xN − 1 e igual que las anteriores contiene
2

la ecuación diferencial y la condición de borde en una misma ecuación.

Al calcular la expansión de Taylor se obtiene que:

para la ecuación centrada en el nodo x 1

2

∂u ∂ 2 u 1 ∂u3 1 ∂u3
− 2 − h 3 + Δt2 3 (3.22)
∂t ∂x 6 ∂x 24 ∂t

el error de truncación es de O(h) + O(Δt2 ).

Para la ecuación centrada en el nodo x 3

2

∂u ∂ 2 u 5 ∂u3 1 ∂u3
− 2 − h2 3 + Δt2 3 (3.23)
∂t ∂x 48 ∂x 24 ∂t

el error de truncación es de O(h2 ) + O(Δt2 ).

Para las ecuaciones centradas en los nodos xi+ 1

2

∂u ∂ 2 u 1 ∂u4 1 ∂u4
− 2 − h2 4 − Δt2 4 (3.24)
∂t ∂x 12 ∂x 8 ∂t

el error de truncación es de O(h2 ) + O(Δt2 ).

44
Los errores de truncación para la antepenúltima ecuación centrada en el nodo xN − 3 y
2

penúltima ecuación centrada en el nodo xN − 1 son iguales a los de las ecuaciones (3.23)
2

y (3.22) respectivamente.
En la siguiente tabla se resumen los órdenes de los errores de truncación de todos
los nodos que conforman la malla unidimensional:

Órdenes Nodos
O(h) + O(Δt2 ) x1 , xN − 1
2 2
2 2
O(h ) + O(Δt ) x3 , xi+ 1 , xN − 3
2 2 2

Tabla 3.1: Órdenes de los errores de truncación 1D

3.3.3. Estabilidad del esquema mimético tipo Crank-Nicolson

En esta sección se presenta el estudio de la estabilidad del esquema mimético tipo

Crank-Nicolson. La obtención de este resultado teórico se logrará a través del estudio
del radio espectral del sistema reducido el cual se obtiene a partir de las ecuaciones
(3.14), (3.15), (3.16), (3.20) y (3.21) detalladas en la sección anterior. Para ello se
demostrará que las matrices de ese sistema conmutan y por tanto son diagonalizables.
Finalmente se utilizará el teorema de Gershgorin [10, 15] para demostrar que el radio
espectral es menor igual ó igual que uno y concluir que el esquema es incondionalmente
estable.
A partir de las ecuaciones (3.14), (3.15), (3.16), (3.20) y (3.21) se obtiene el sistema
reducido:
Aun+1 = Bun (3.25)

en el cual la matriz A es:

45
⎛ 2
⎞
h2
− 3h
16 +
h
4 + 1
2 + 1
Δt 48 − h
12 − 1
2 0 ··· 0
⎜ ⎟
⎜ h
− 16 − 1 h
+1+ 1
− 12 0 ··· 0⎟
⎜ 2 16 Δt ⎟
⎜ ⎟
⎜ 0 − 12 1 + Δt 1
− 21 ··· 0 ⎟
⎜ ⎟
⎜ .. .. .. .. ⎟
⎜ . . . . ⎟
⎜ ⎟
⎜ 0 − 12 1+ 1
− 21 0 ⎟
⎜ Δt ⎟
⎜ ··· h
− 16 − 1 h 1
− 21 ⎟
⎝ 0 0 0 2 16 +1+ Δt ⎠
2 2
h h 1
0 ··· 0 0 0 48 − 12 − 2 − 3h
16 +
h
4 + 1
2 + 1
Δt

y la matriz B es:

⎛ ⎞
3h2 h 1 1 2

16 − 4 − 2 + Δt − h48 + h
12 + 1
2 0 ··· 0
⎜ ⎟
⎜ h
+ 1 h
− 16 −1+ 1 1
0 ··· 0 ⎟
⎜ 16 2 Δt 2 ⎟
⎜ 1 1 1 ⎟
⎜ 0 2 −1 + Δt 2 ··· 0 ⎟
⎜ ⎟
⎜ .. .. .. .. .. ⎟
⎜ . . . . . ⎟
⎜ ⎟
⎜ 0 1
−1 + 1 1
0 ⎟
⎜ 2 Δt 2 ⎟
⎜ ··· h 1 h
− 16 −1+ 1 1 ⎟
⎝ 0 0 0 16 + 2 Δt 2 ⎠
2
3h2
0 ··· 0 0 0 − h48 + h
12 + 1
2 16 − h
4 − 1
2 + 1
Δt

En vista de la complejidad de las matrices anteriormente descritas, es conveniente

descomponerlas repectivamente como sigue:
 
1 1
A= I + 2W (3.26)
Δt h

 
1 1
B= I − 2W (3.27)
Δt h
donde I es la matriz identidad y la matriz W es:

46
 ⎛ 2
⎞
h2
− 3h
16 +
h
4 + 1
2 48 − h
12 − 1
2 0 ··· 0
⎜ ⎟
⎜ h
− 16 − 1 h
+1 − 12 0 ··· 0⎟
⎜ 2 16 ⎟
⎜ ⎟
⎜ 0 − 21 1 − 12 ··· 0 ⎟
⎜ ⎟
⎜ .. .. .. .. ⎟
⎜ . . . . ⎟
⎜ ⎟
⎜ 0 − 12 1 − 21 0 ⎟
⎜ ⎟
⎜ ··· h
− 16 − 1 h
− 12 ⎟
⎝ 0 0 0 2 16 +1 ⎠
2 2
h h 1
0 ··· 0 0 0 48 − 12 − 2 − 3h
16 +
h
4 + 1
2

Al sustituir la descomposición de las matrices (3.26) y (3.27) en la ecuación (3.25) se

obtiene el nuevo sistema:

   
Δt n+1 Δt
I+ 2W u = I − 2 W un (3.28)
h h
Δt
Llamando r = h2
se tiene:

(I + rW ) un+1 = (I − rW ) un (3.29)

Se demostrará ahora que las matrices del sistema (3.29) conmutan y por tanto son
diagonalizables [6], según lo indı́ca el :

Teorema 3.2 Dos matrices son simultáneamente diagonalizables en una misma base
si y solo si conmutan.

En efecto, como las matrices conmutan:

(I + rW ) (I − rW ) = I−rW +rW −r 2W 2 = I+rW −rW −r 2W 2 = (I − rW ) (I + rW )
entonces son simultaneamente diagonalizables.
El resultado anterior implica que los autovalores de (I + rW )−1 (I − rW ) son de la
forma:
(1 − rλ)
μ= (3.30)
(1 + rλ)

47
donde λ representan los autovalores (reales o complejos) de la matriz W . Falta entonces
demostrar que esos autovalores tienen parte real mayor ó igual que cero para que
|μ| ≤ 1. Para lograr esto, se utilizará el:

Teorema 3.3 (Discos de Gershgorin)

Sea A = a(i, j), i, j = 1..n una matriz de orden n × n. Se define


n
r(i) = |a(i, j)|
j=1
i=j

para i = 1..n. Entonces cada autovalor (real o complejo) λ de A satisface la siguiente

desigualdad:

|λ − a(i, i)| ≤ r(i) (3.31)

A continuación se considerarán las filas de la matriz W para aplicarle el teorema

3.3 en los siguientes casos:

Caso 1: λ ∈ R
    2 
 3h2
h 1  h h 1 
λ − − + + ≤ − −  ⇒
 16 4 2   48 12 2

  
 3h2
h 1  2
 λ− − + +  ≤ −h + h + 1 ⇒ (3.32)
 16 4 2  48 12 2

−3h2 h 1 h2 h 1
+ + + − − ≤λ ⇒
16 4 2 48 12 2

h2 h
λ≥− + ≥0
16 6

48
siempre que 0 ≤ h ≤ 1.

Para la segunda fila, se tiene:

         
 h   1  h 1   h  1 h 1
λ − + 1  ≤ −  + − −  ⇒  λ − + 1  ≤ + +
 16   2   16 2  16 2 16 2

h h h h
+1− −1 ≤λ ≤ +1+ +1 ⇒ λ ≥ 0. (3.33)
16 16 16 16
Finalmente para las filas centrales aplicando nuevamente el teorema 3.3, se ob-
tiene:

   
 1  1
|λ − 1| ≤ −  + −  ⇒ |(λ − 1)| ≤ 1 ⇒ −1 + 1 ≤ λ ≤ 1 + 1 ⇒ λ ≥ 0.
2 2
(3.34)
Para la antepenúltima, penúltima y última filas el resultado es idéntico pues
existe simetrı́a en el problema.

Caso 2: λ ∈ C
    2 
 3h2
h 1  h h 1 
λ − − + + ≤ − −  ⇒
 16 4 2   48 12 2

   
 3h2
h 1  2
 x− − + + + iy  ≤ −h + h + 1 ⇒ (3.35)
 16 4 2  48 12 2

  2  2 2
3h2 h 1 2 h h 1
0≤ x− − + + +y ≤ − + + ⇒
16 4 2 48 12 2

  
3h2 h 1 3h2 h 1
x− − + + ≥0 ⇒x≥− + +
16 4 2 16 4 2
siempre que 0 ≤ h ≤ 1 se tiene que Re(λ) = x ≥ 0.

49
 Para la segunda fila, se tiene:
          
        1 h 1
λ − h + 1  ≤ − 1 +− h − 1   h
⇒ x− − +1 + iy  ≤ + +
 16   2   16 2  16 2 16 2

  2  2   
h 2 1 h 1 h
0≤ x− − +1 +y ≤ + + ⇒0≤ x− − +1
16 2 16 2 16
(3.36)
 
h
x≥ − +1 ⇒x≥0
16
siempre que 0 ≤ h ≤ 1.

Finalmente para las filas centrales se tiene:

   
 1  1
|λ − 1| ≤ −  + −  ⇒ |(x − 1) + iy| ≤ 1 ⇒ 0 ≤ (x − 1)2 + y 2 ≤ 1 (3.37)
2 2

⇒ (x − 1) ≥ 0 ⇒ x ≥ 1.

Nuevamente por la simetrı́a del problema para la antepenúltima, penúltima y

última filas el resultado es idéntico.

Luego, de manera general se tiene que:

    
 1 − rλ   (1 − rx) − riy 
 = (1 − rx) + r y .
2 2 2
 = (3.38)
 1 + rλ   (1 + rx) + riy  (1 + rx)2 + r 2 y 2
Ahora bien


(1 − rx)2 + r 2 y 2 (1 − rx)2 + r 2 y 2
≤1 ⇒ ≤1 ⇒ −2rx ≤ 2rx ⇒ x ≥ 0.
(1 + rx)2 + r 2 y 2 (1 + rx)2 + r 2 y 2
(3.39)
Finalmente, como |μ| ≤ 1 se puede concluir que analı́ticamente el esquema mimético
tipo Crank-Nicolson es incondicionalmente estable.

50
3.3.4. Estabilidad Numérica

El estudio de la estabilidad numérica del esquema mimético tipo Crank-Nicolson

que se está proponiendo se realizará a través de la observación del radio espectral de la
matriz asociada al sistema; en efecto después de haber realizados algunos experimentos
numéricos con sistemas de orden 5 × 5, 10 × 10 y 50 × 50 los resultados indican que
el radio espectral de las matrices asociadas a dichos sistemas siempre es menor o igual
que uno. La figura 3.9 muestra el módulo de los autovalores asociados a la matriz del
sistema reducido de orden 50 × 50 para un intervalo de tiempo Δt = 0, 01. Nótese que
todos los autovalores están concentrados en el intervalo (0.996,1).

Figura 3.9: Radio espectral para un sistema de órden 50 × 50

En la siguiente figura se muestran los resultados obtenidos en el experimento en

donde se ha reducido el tamaño de la malla a 5 × 5 y se ha aumentado el intervalo de
tiempo Δt. Como se puede observar el radio espectral de la matriz asociada al sistema

51
nuevamente es menor ó igual que uno. Sin embargo se puede notar que los módulos de
los autovalores están contenidos en el intervalo (0.2, 0.9) que es más grande respecto
al ejemplo anterior.

Figura 3.10: Radio espectral para un sistema de órden 5 × 5

Por tanto se puede concluir que independientemente del tamaño de la malla y del
paso del tiempo, el espectro de la matriz asociada al sistema permanece menor ó igual
que uno, ratificando una vez más que el esquema propuesto es incondicionalmente
estable.
Después de haber demostrado, en las secciones anteriores, que el esquema es consi-
stente e incondicionalmente estable, se utiliza el teorema 3.1 para afirmar que el nuevo
esquema mimético tipo Crank-Nicolson es convergente.

52
Capı́tulo 4

Resultados numéricos

4.1. Introducción
En este capı́tulo se dará una descripción del esquema de diferencias finitas tradi-
cionales (DF), luego se presentará un estudio comparativo entre éste y el método
mimético tipo Crank-Nicolson (MIMÉTICO); posteriormente se mostrarán los resul-
tados numéricos para la ecuación no estática de difusión en una y dos dimensiones.
Finalmente, se demostrará numéricamente la convergencia cuadrática del nuevo esque-
ma.

4.2. Método de diferencias finitas

Este método consiste en aproximar las derivadas que aparecen en las ecuaciones
diferenciales por las aproximaciones truncadas de sus desarrollos de Taylor en cada
uno de los nodos presentes en la malla de discretización [11]. En los nodos interiores las
discretizaciones que se realizan con este esquema son muy precisas pero en los bordes
esto no ocurre con tanta rigurosidad. En efecto, se procede a extender las mallas con

53
puntos fantasmas o imaginarios ubicados en el exterior del dominio de estudio. En
la siguiente figura se muestra la malla unidimensional de diferencias finitas, que se
diferencia de la malla mostrada en la figura 2.1 por tener puntos imaginarios (marcados
por cı́rculos) a la izquierda y derecha de los extremos del intervalo [0,1] respectivamente.

Figura 4.1: Malla referencial para el esquema de diferencias finitas 1D

El sistema asociado a la malla de la figura 4.1 está compuesto por una mayor can-
tidad de ecuaciones, en cambio en el sistema del método mimético existen N + 2 ecua-
ciones, N que provienen de los nodos internos y dos que pertenecen a las condiciones
de borde. Una situación similar ocurre para la malla bidimensional que se muestra a
continuación:

Figura 4.2: Malla referencial para el esquema de diferencias finitas 2D

54
La diferencia entre esta malla y la que se utiliza en el método mimético (figura 2.2)
radica en los bordes. En efecto la malla de diferencias finitas (figura 4.2) posee puntos
fantasmas en cada lado del cuadrado y nodos en las esquinas; este hecho influye de
manera importante en la cantidad de ecuaciones que conforman el sistema lo cual hace
que la resolución numérico-computacional se vuelva más pesada y menos eficiente.

4.3. Análisis 1D
El problema de contorno que se plantea resolver viene dado por la ecuación dife-
rencial:

∂u ∂ 2 u
− 2 = −e−t x (4.1)
∂t ∂x
definida en el intervalo (0, 1) y su solución debe satisfacer las condiciones de borde tipo
Robin:
∂u
u(0) − (0) = −e−t (4.2)
∂x
∂u
u(1) + (1) = 2e−t (4.3)
∂x

en los extremos del intervalo y condición inicial:

u(x, 0) = x (4.4)

Las ecuaciones (4.1)-(4.4) determinan un problema unidimensional bien definido cuya

solución analı́tica viene dada por:

u(x, t) = e−t x (4.5)

La ecuación (4.5) permitirá evaluar la calidad de las aproximaciones producidas por

los esquemas que se están estudiando. En la siguiente tabla se muestran los errores en
norma infinito cometidos por cada método.

55
 Métodos Tamaño de malla Errores
30 1.4022 ×10−5
MIMÉTICO 60 3.5048 ×10−6
100 1.2618 ×10−6
30 1.6500×10−2
DF 60 8.3000 ×10−3
100 5.0000×10−3

Tabla 4.1: Errores en norma infinito cometidos por ambos métodos

Se puede observar que a medida que se refina la malla el esquema mimético tipo
Crank-Nicolson produce errores con hasta tres órdenes de magnitud menos. Esto evi-
dencia la superioridad del esquema sobre el método tradicional de diferencias finitas.

Métodos Tamaño de malla Órdenes

30 2.0008
MIMÉTICO 60 2.0004
100 2.0002
30 1.1167
DF 60 1.0802
100 1.0607

Tabla 4.2: Órdenes de convergencia 1D

La tabla 4.2 resume los resultados de los órdenes de convergencia para distintos
tamaños de malla. En el caso del esquema mimético tipo Crank-Nicolson se obtiene
un órden de convergencia cuadrático, mientras que el esquema de diferencias finitas
produce un órden de convergencia lineal.

56
 La representación gráfica de los errores de ambos métodos se muestra en la figura
que sigue:

Figura 4.3: Comparación de los errores numéricos en norma infinito 1D

Los órdenes de convergencia vienen dados por las pendientes de las curvas y mientras
más abajo se encuentra una curva mejor será el esquema que la representa. La curva
de color verde representa los errores cometidos por el esquema de diferencias finitas
tradicional, mientras que la curva azul representa los errores del nuevo esquema siendo
éstos los de menor magnitud.

4.4. Análisis 2D
El problema de contorno que se plantea resolver viene dado por una extensión
natural del problema unidimensional estudiado en la sección anterior y es:
 2 
∂u ∂ u ∂2u
− + = e−t (x − y) (4.6)
∂t ∂x2 ∂y 2

57
en la región (0, 1) × (0, 1) y su solución debe satisfacer las condiciones de borde tipo
Robin:  
∂u
u− (x, 0, t) = −e−t x − e−t (4.7)
∂y
 
∂u
u+ (x, 1, t) = −e−t (x − 1) + e−t (4.8)
∂y
 
∂u
u+ (1, y, t) = −e−t (1 − y) − e−t (4.9)
∂x
 
∂u
u− (0, y, t) = e−t y + e−t (4.10)
∂x

para cada uno de los lados del cuadrado unitario y condición inicial:

u(x, y, 0) = −x + y (4.11)

Las ecuaciones (4.6)-(4.10) determinan un problema bidimensional bien definido cuya

solución analı́tica viene dada por:

u(x, y, t) = −e−t (x − y) (4.12)

En la siguiente tabla se muestran los errores en norma infinito cometidos por cada
método para distintos tamaños de malla.

58
 Métodos Tamaño de malla Errores
30 1.7302 ×10−6
MIMÉTICO 60 4.4296 ×10−7
100 1.5590×10−7
30 3.1559 ×10−4
DF 60 1.5686 ×10−4
100 9.3892 ×10−5

Tabla 4.3: Errores en norma infinito cometidos por ambos métodos

En ella se evidencia que a medida que se va refinando la malla ambos métodos

producen menores errores que el caso unidimensional, sin embargo los del esquema
mimético resultan ser los más pequeños. A continuación se muestra una tabla que
contiene los órdenes de convergencia de los métodos bajo estudio.

Métodos Tamaño de malla Órdenes

30 1.9613
MIMÉTICO 60 1.9922
100 1.9957
30 1.0518
DF 60 1.0294
100 1.0202

Tabla 4.4: Órdenes de convergencia 2D

Se puede observar que las tasas de convergencia para el nuevo método son casi
cuadráticas mientras que las del método de diferencias finitas tradicional son lineales.
Una interpretación gráfica de las tasas de convergencia se puede apreciar en la siguiente
figura:

59
 Figura 4.4: Comparación de los errores numéricos en norma infinito 2D

En ella se observa que al igual que el caso unidimensional la recta de color verde
está asociada al esquema de diferencias finitas tradicional y representa los errores
cometidos por ese esquema mientras que la recta de color azul representa los errores
cometidos por el nuevo esquema mimético los cuales nuevamente son de menor ma-
gnitud.

4.5. Conclusiones
En esta tesis se ha presentado un nuevo y original esquema mimético de diferencias
finitas tipo Crank-Nicolson para la resolución de la ecuación no estática de difusión
en una y dos dimensiones bajo las condiciones de borde tipo Robin. Se ha demostra-
do analı́ticamente que el nuevo esquema es consistente con la ecuación diferencial, es
incondicionalmente estable y por lo tanto utilizando el teorema de Lax, el esquema
propuesto es convergente. El método ha sido implementado en MATLAB en sus ver-

60
siones unidimensional y bidimensional para realizar un estudio numérico comparativo
con el esquema de diferencias finitas tradicional el cual ha sido esbozado brevemente.
Los resultados de esta comparación evidencian la superioridad del esquema mimético
tipo Crank-Nicolson tanto en la magnitud de los errores cometidos como en el órden
de convergencia. Esto permite confirmar una vez más que los métodos miméticos son
la mejor opción para resolver la ecuación no estática de difusión ya que permiten
tratar rigurosamente las condiciones de borde asociadas a las ecuaciones diferenciales
sin necesidad de utilizar puntos imaginarios o fantasmas.

61
Bibliografı́a

[1] Castillo J.E. y Grone S., A Matrix Approach to Higher Order Approximations for
Divergence and Gradient Satisfying a Global Conservation Law, SIAM Journal of
Matrix Analysis and Applications, 2003.

[2] Castillo J.E. y Yasuda M., A Comparison of two Matrix Operator Formulation for
Mimetic Divergence and Gradient Discretization, Journal of Mathematical Mo-
delling and Algorithm, 4, 2005.

[3] Duchateau P. y Zachmann D.W., Partial Differential Equations. McGraw Hill,

NY, 1988.

[4] Guevara-Jordan J.M., Rojas S., Freites-Villegas M. and Castillo J.E., Convergence
of a Mimetic Finite Difference Method for Static Diffusion Equations, Advances
in Difference Equations, 2007.

[5] Guevara-Jordan J.M. y Arispe J., A Second Order Mimetic Approach for Tracer
Flow in Oil Reservoirs. Society of Petroleum Engineers. Paper No. 107366. 2007.

[6] Hoffman K. y Kunze R., Álgebra lineal. Editorial Prentice Hall, 205-206, 1973.

[7] Hyman J.M. y Shashkov M.J., The Approximation of Boundary Conditions for
Mimetic Finite Difference Methods Operators, Computers and Mathematics with
Applications, 1998.

62
 [8] Hyman J.M. y Steinberg S., Mimetic Finite Difference for Diffusion Equations,
Computational Geoscience, 47,2002.

[9] Incropera, F. y DeWitt, D., Fundamentos de Transferencia de Calor, Prentice

Hall, 1999.

[10] Isaacson E. y Keller B.E., Analysis of Numerical Methods, John Wiley and Sons,
Inc., New York. 1966.

[11] Lapidus L. y Pinder G.F., Numerical Solution of Partial Differential Equations in

Science and Engineering, Wiley, NY, 1983.

[12] Mannarino I., Guevara-Jordan J.M. y Quintana Y., A numerical study of a

new mimetic scheme for unsteady heat equation. Sometido a evaluación en FA-
RAUTE.2007.

[13] Mannarino I., Guevara-Jordan J.M. y Quintana Y., A mimetic Crank-Nicolson

scheme for transient diffusion equation.Preprint 2007.

[14] Morton K.W. y Mayers D.F., Numerical Solution of Partial Differential Equation,
Oxford University Press, 1994.

[15] Saad Y., Iterative Methods for Sparse Linear Systems, 2000.

[16] Samarsky A.A. y Andreiev E.S., Método en Diferencias Finitas para Ecuaciones
Elı́pticas, Editorial Mir, 1979.

[17] Shashkov M., Conservative Finite Difference Methods on General Grids, CRC
Press, 1996.

[18] Striwerda J.C., Finite Difference Schemes and Partial Differential Equation,
Wadsworth-Brooks Cole, 1989.

63
[19] Universidad de Castilla la Mancha Resolución numércia de ecuaciones
en derivadas parciales. Disponible en http:// www.matematicas.uclm.es/ind-
cr/metnum. 2007.

[20] Villegas M., Guevara-Jordan J.M., Rojas O.R., Castillo J.E. y Rojas S.J., A
Mimetic Finite Difference Scheme for Solving the Steady State Diffusion Equation
with Singular Functions, Memorias del VII International Conference of Numerical
Methods in Engineering and Applied Science, 2004.

64