XXI Congreso de Ecuaciones Diferenciales y Aplicaciones

XI Congreso de Matemática Aplicada

Ciudad Real, 21-25 septiembre 2009
(pp. 1–8)

Resolución de la ecuación de advección-difusión en 2-D

utilizando diferencias finitas generalizadas. Consistencia y
Estabilidad.

F. Ureña1 , J.J. Benito2 , L. Gavete3

1
Dpto. de Matemáticas, Universidad de Castilla-La Mancha, Ciudad Real. E-mail:
[email protected].
2
Dpto. de Construcción y Fabricación, Universidad Nacional de Educación a Distancia, Madrid. E-mail:
[email protected].
3
Dpto. de Matemática Aplicada a los Recursos Naturales, Universidad Politécnica de Madrid. E-mail:
[email protected].

Palabras clave: diferencias finitas generalizadas, advección-difusión, método explı́cito, estrella.

Resumen

En esta comunicación se presenta la utilización del Método de Diferencias Finitas

Generalizadas para la resolución de la ecuación de advección-difusión, para 2-D. La
comunicación se inicia con la obtención de las expresiones explı́citas en diferencias
finitas generalizadas. A partir de estas expresiones se estudia el error de truncamien-
to, consistencia, estabilidad y convergencia. En la comunicación se incluyen algunos
resultados, de entre los numerosos casos analizados, como ejemplos representativos
de la resolución de la ecuación de advección-difusión que pretenden ilustrar el buen
comportamiento del método.

1. Introducción
El método de diferencias finitas generalizadas (GFDM), como generalización del méto-
do de diferencias finitas, permite la aplicación en dominios irregulares. Al desarrollo de
este método han contribuido los trabajos de Benito, Ureña y Gavete [1, 2].
Los artı́culos [3, 5] muestran la aplicación del método de diferencias finitas generalizadas
a la resolución de ecuaciones en derivadas parciales dependientes del tiempo.

1
 F. Ureña, J.J. Benito, L. Gavete

Figura 1: Estrella con nodo central (x0 , y0 )

2. Diferencias finitas generalizadas y método explı́cito

Se considera la resolución numérica de la ecuación de advección-difusión para la función
U (x, y, t)

∂U (x, y, t) ∂U (x, y, t) ∂U (x, y, t)

+ cx + cy =
∂t ∂x ∂y
∂ 2 U (x, y, t) ∂ 2 U (x, y, t)
α( + ) t > 0, (x, y) ∈ Ω ⊂ R2 (1)
∂x2 ∂y 2
con la condición inicial
U (x, y, 0) = f (x, y) (2)
y la condición de contorno
∂U (x0 , y0 , t)
aU (x0 , y0 , t) + b = g(t) in Γ (3)
∂n
siendo f (x, y) y g(t) dos funciones conocidas, a, b son constantes y Γ la frontera del dominio
Ω, donde α > 0 es el coeficiente de difusión y cx > 0, cy > 0 son velocidades constantes.
A continuación se van a obtener las fórmulas explı́citas en diferencias finitas generalizadas
para la aproximación de las derivadas parciales en los puntos del dominio. En primer lugar,
se discretiza el dominio Ω ∪ Γ. Se define el nodo central con un conjunto de nodos a su
alrededor, al conjunto de dichos nodos se le denomina estrella, estableciendo una relación
entre una estrella y su nodo central. De esta forma, a cada nodo en del dominio tiene
asociada una estrella, en la cual es el nodo central (ver figura 1).
Si U0 es el valor de la función en el nodo central de la estrella y Uj son los valores de las
funciones en el resto de los nodos, con j = 1, · · · , 8, entonces, de acuerdo con la serie de
expansión de Taylor

∂U0 ∂U0 h2j ∂U02 kj2 ∂U02 ∂U02

Uj = U0 + hj + kj + + + h j kj + ··· (4)
∂x ∂y 2 ∂x2 2 ∂y 2 ∂x∂y
donde (x0 , y0 ) son las coordenadas espaciales del nodo central, (xj , yj ) las coordenadas del
nodo j en la estrella, hj = xj − x0 , kj = yj − y0 .

2
 Resolución de la ecuación de advección-difusión en 2-D utilizando GFDM

Si en la ecuación 4 los términos de orden superior al segundo son eliminados, se obtiene

la aproximación de segundo orden para Uj , si se representa este valor por uj , entonces es
posible definir
8
X ∂u0 ∂u0 h2j ∂u20 kj2 ∂u20 ∂U02
B(u) = [(u0 − uj + hj + kj + + + h j kj )w(hj , kj )]2 (5)
∂x ∂y 2 ∂x2 2 ∂y 2 ∂x∂y
j=1

donde w(hj , kj ) es la función de ponderación.

Si la expresión 5 es minimizada con respecto a las derivadas parciales, se obtiene el siguiente
sistema de ecuaciones lineales
ADu = b (6)
Resolviendo el sistema 6 y teniendo en cuenta que hj = kj = h, se obtienen las siguientes
fórmulas finitas generalizadas para las derivadas parciales

∂un0 1
= [8un + un2 − un4 − 8un5 − un6 + un8 ] (7)
∂x 20h 1
∂un0 1 n
= [u + 8un3 + un4 − un6 − 8un7 − un8 ] (8)
∂y 20h 2
∂ 2 un0 ∂ 2 un0 1
2
+ 2
= 2 (−20un0 + 4un1 + un2 + 4un3 + un4 + 4un5 + un6 + 4un7 + un8 ) (9)
∂x ∂y 6h
Aproximando la derivada respecto del tiempo en el nodo central de la estrella por la
derivada de avance, y designando por un0 y un+1
0 los valores aproximados de la función
U (x, y, t) en el nodo central de coordenadas espaciales (x0 , y0 ) para los tiempos n4t y
(n + 1)4t respectivamente, se tiene

∂U (x0 , y0 , n4t) un+1 − un0

= 0 (10)
∂t 4t

Sustituyendo las ecuaciones 7, 8, 9 y 10 en la ecuación 1, se obtiene la ecuación lineal

un+1 − un0 1 1 n
0
+ cx [8un1 + un2 − un4 − 8un5 − un6 + un8 ] + cy [u + 8un3 + un4 − un6 − 8un7 − un8 ]
4t 20h 20h 2
1
+ α 2 (−20un0 + 4un1 + un2 + 4un3 + un4 + 4un5 + un6 + 4un7 + un8 ) (11)
6h
La expresión 11 relaciona el valor aproximado de la función en el nodo central de la es-
trella, en el paso de tiempo n + 1, con los valores aproximados de la función en los nodos
de la estrella, en el paso de tiempo n.

3. Convergencia
De acuerdo con el teorema de equivalencia de Lax [4], si la condición de consistencia
es satisfecha, la estabilidad es necesaria y suficiente para la condición de convergencia.

3
 F. Ureña, J.J. Benito, L. Gavete

3.1. Error de truncamiento

El error de truncamiento es la suma de los errores de truncamiento debido a las fórmulas

en diferencias finitas empleadas para las derivadas parciales temporales y espaciales. Si se
designan T Et y T E(x,y) dichos errores de truncamiento, respectivamente, se tiene

4t ∂ 2 U (x, y, t1 )
T Et = − + Θ((4t)2 ), n4t < t1 < (n + 1)4t (12)
2 ∂t2

Para obtener el error de truncamiento para las derivadas espaciales, en la serie de expansion
de Taylor se incluyen los términos hasta de cuarto orden. Si se designa por B ∗ (u) la
expresión 5 en la cual se han incluido los nuevos términos, y minimizando dicha expresión
respecto de las derivadas parciales de primer y segundo orden, se obtiene

2h4
 
0 0 0 0
 5 
 2h4 
 0 0 0 
5
 
5h2 h2
 
T E(x,y) = {−cx , −cy , α, α, 0}  ×
 − 0 
 3 3 
5h2
 
 
 SY M 0 
3
2h2
1 5 ∂ 3 U (x1 , y1 , t) 1 ∂ 3 U (x1 , y1 , t)
 
1
− 2( +3 ) − (0) − · · ·

 6h 2 ∂x3 2 ∂x∂y 2 24 

1 1 ∂ 3 U (x , y , t) 5 ∂ 3 U (x , y , t) 1
1 1 1 1
 
 − 2 (3 + ) − (0) − · · · 

 6h 2 ∂x2 ∂x 2 ∂y 3 24 

 1 4 4 4
− (0) − 1 ( 5 ∂ U (x1 , y1 , t) + 3 ∂ U (x1 , y1 , t) + 1 ∂ U (x1 , y1 , t) ) − · · · (13)

 12 48 2 ∂x 4 2
∂x ∂y 2 2 ∂y 4 
 4 4 4

 1 1 1 ∂ U (x 1 , y 1 , t) ∂ U (x 1 , y 1 , t) 5 ∂ U (x1 , y 1 , t) 
− (0) − ( + 3 + ) − · · ·
 12 48 2 ∂x4 ∂x2 ∂y 2 2 ∂y 4 
4 4
 
 1 1 1 ∂ U (x1 , y1 , t) 1 ∂ U (x1 , y1 , t) 
− (0) − (4 3
+4 3
) − ···
6 24 2 ∂x ∂y 2 ∂x∂y

2h4 1 5 ∂ 3 U (x1 , y1 , t) 1 ∂ 3 U (x1 , y1 , t)

T E(x,y) = −cx [− 2 ( + 3 ))]
5 6h 2 ∂x3 2 ∂x∂y 2
2h4 1 3 ∂ 3 U (x1 , y1 , t) 5 ∂ 3 U (x1 , y1 , t)
− cy [− 2 ( + ))]
5 6h 2 ∂x2 ∂y 2 ∂y 3
h2 ∂ 4 U (x1 , y1 , t) ∂ 4 U (x1 , y1 , t) ∂ 4 U (x1 , y1 , t)
−α ( + 2 + ) + Θ2 (h4 ) (14)
12 ∂x4 ∂x2 ∂y 2 ∂y 4

donde (x1 , y1 ) es un punto del interior del dominio definido por la estrella.
La ecuación 14 es el error de truncamiento para las derivadas espaciales.

4
 Resolución de la ecuación de advección-difusión en 2-D utilizando GFDM

Por tanto, el error de truncamiento total, T T E, viene dado por

4t ∂ 2 U (x, y, t1 ) 2 2h4 1 5 ∂ 3 U (x1 , y1 , t)

T T E = T Et + T E(x,y) = − + Θ 1 ((4t) ) − cx [− (
2 ∂t2 5 6h2 2 ∂x3
3
3 ∂ U (x1 , y1 , t) 2h 4 3 3
1 3 ∂ U (x1 , y1 , t) 5 ∂ U (x1 , y1 , t)
+ 2
))] − c2 [− 2 ( + ))]
2 ∂x∂y 5 6h 2 ∂x2 ∂y 2 ∂y 3
h2 ∂ 4 U (x1 , y1 , t) ∂ 4 U (x1 , y1 , t) ∂ 4 U (x1 , y1 , t)
−α ( 4
+2 + ) + Θ2 (h4 ) (15)
12 ∂x ∂x2 ∂y 2 ∂y 4

3.2. Consistencia
De acuerdo con la expresión 15 se tiene
lı́m TTE → 0 (16)
(4t,h)→(0,0)

Por tanto, el método es consistente.

A continuación se estudia la estabilidad utilizando el análisis de von Neumann.

3.3. Estabilidad
Si se define
Tx Tx
un0 = ξ n eiν 0
; unj = ξ n eiν j
(17)
donde ν = (νx , νy )T es el vector columna de los números de onda, x0 = (x0 , y0 ) es el vector
de las coordenadas del nodo central de la estrella y xj = (xj , yj ) son las coordenadas del
resto de los nodos de la estrella, con xj = x0 + hj .
además, ξ es denominado factor de amplificación. Si el módulo del factor de amplificación
es mayor que la unidad, (kξk > 1, el método es inestable.
Sustituyendo 17 into 11, se tiene
Tx Tx 1Tx
ξ n+1 eiν 0
= ξ n eiν 0
− ξ n eiν [cx (8eihνx + eih(νx +νy ) − eih(−νx +νy ) − 8e−ihνx −
0
4t
20h
e−ih(νx +νy ) + eih(νx −νy ) ) + cy (eih(νx +νy ) + 8eihνy + eih(−νx +νy ) − e−ih(νx +νy ) −
T α
8e−ihνy − eih(νx −νy ) )] + ξ n ei{νj } {x0 } 2 4t[−20 + 4eihνx + eih(νx +νy ) + 4eihνy + eih(−νx +νy )
6h
+ 4e−ihνx + e−ih(νx +νy ) + 4e−ihνy + eih(νx −νy ) ] (18)
T {x
Simplificando por ξ n ei{νj } 0} y operando se obtiene
4t
ξ =1−i [cx sen hνx (4 + cos hνy ) + cy (sen hνy (4 + cos hνy )]+
5h
α4t
[−20 + 8 cos hνx + 8 cos hνy + 4 cos hνx cos hνy )] (19)
6h2

α4t
− 1 ≤ Re(ξ) ≤ 1 ⇔ −1 ≤ 1 + [−20 + 8 cos hνx + 8 cos hνy + 4 cos hνx cos hνy )] ≤ 1
6h2
1
⇔ 0 ≤ 4t ≤ (20)
20
α 2
6h

5
 F. Ureña, J.J. Benito, L. Gavete

El módulo del factor de amplificación

α4t hνx hνy h(νx + νy ) h(νx − νy )
kξk = k1 − [16(sin2 ( ) + sin2 ( )) + 4(sin2 ( ) + sin2 ( ))]
6h2 2 2 2 2
4t
− i [cx sen hνx (4 + cos hνy ) + cy sen hνy (4 + cos hνx )]k (21)
5h
imponiendo la condición de estabilidad, se tiene
4t
( [cx sen hνx (4 + cos hνy ) + cy sen hνy (4 + cos hνx )])2
5h
α4t hνx hνy h(νx + νy ) h(νx − νy )
+ ( 2 [16(sin2 ( ) + sin2 ( )) + 4(sin2 ( ) + sin2 ( ))])2
6h 2 2 2 2
α4t hνx hνy h(νx + νy ) h(νx − νy )
≤ 2 2 [16(sin2 ( ) + sin2 ( )) + 4(sin2 ( ) + sin2 ( ))] (22)
6h 2 2 2 2
y operando, se tiene
4t[cx + cy ]2 10
≤ (23)
α 3
Por tanto el valor del paso de tiempo, 4t, para que el método sea estable viene dado por
1 10α
4t = min{ , } (24)
20 3[cx + cy ]2
α 2
6h

4. Resultados numéricos
En esta sección se muestran dos ejemplos de resolución numérica de´ecuaciones de
advección-difusión en 2-D. En ambos casos, la función de ponderación utilizada ha sido
1
w(hj , kj ) = q (25)
(h2j + kj2 )3

y el criterio de selección de los nodos el del cuadrante. El error global ha sido calculado
para cada paso de tiempo usando la siguiente norma
q PN T
2
j=1 (sol(j)−exac(j))
NT
Error global = × 100 (26)
|exacmax |
where sol(j) es el valor de la solución aproximada en el nodo j, exac(j) es la valor de
la solución exacta en el nodo j, exacmax es el máximo valor de la solución exacta en los
nodos interiores de de la malla considerada y N T es el número de nodos del interior.

4.1. Ejemplo 1

∂U (x, y, t) ∂U (x, y, t) ∂U (x, y, t)

+ + =
∂t ∂x ∂y
∂ 2 U (x, y, t) ∂ 2 U (x, y, t)
0,1( + ) t > 0, 0 < x, y < 1 (27)
∂x2 ∂y 2

6
 Resolución de la ecuación de advección-difusión en 2-D utilizando GFDM

o
Figura 2: Error global versus n nodos ; Error global versus 4t

con la condición inicial

U (x, y, 0) = ex+y (28)
y las condiciones de contorno Dirichlet, siendo la solución exacta

U (x, y, t) = e−1,8t+x+y (29)

En la figura 2 se muestra, manteniendo fijo el paso de tiempo (4t = 0,0001), la disminución

del error global al aumentar el número de nodos en la malla. También, en la figura 2 se
muestra la disminución del error global al disminuir el paso de tiempo para la malla de
441 nodos.

4.2. Ejemplo 2
La resolución de la ecuación 27 con la condición inicial

U (x, y, 0) = sen(x − y) (30)

y las condiciones de contorno Dirichlet, siendo la solución exacta

U (x, y, t) = e−0,2t sen(x − y) (31)

En la figura 3 se muestra, manteniendo fijo el paso de tiempo (4t = 0,0001), la disminución

del error global al aumentar el número de nodos en la malla. También, en la figura 3 se
muestra la disminución del error global al disminuir el paso de tiempo para la malla de
441 nodos.

5. Conclusiones
En esta comunicación se ha obtenido el error de truncamiento y, por tanto, la consis-
tencia ha sido demostrada. Igualmente, se ha obtenido el criterio de estabilidad utilizando
el análisis de von Neumann.
Los ejemplos resueltos, de los numerosos a los que se ha aplicado el GFDM, muestran su
buen comportamiento.

7
 F. Ureña, J.J. Benito, L. Gavete

Figura 3: Error global versus no nodos ; Error global versus 4t

Agradecimientos
Los autores agradecen la ayuda recibida del Ministerio de Ciencia e Innovación de
España en el proyecto TISMANCA, Ref.: CGL2008-01757/CLI.

Referencias
[1] J.J. Benito, F. Ureña, L. Gavete, Influence of several factors in the generalized finite difference
method. Applied Mathematical Modelling,2512,1039-1053(2001).
[2] J.J. Benito, F. Ureña, L. Gavete, R. Alvarez, An h-adaptive method in the generalized finite differ-
ences. Computer Methods in Applied Mechanics and Engineering, 192,735-759(2003).
[3] J.J. Benito, F. Ureña, L. Gavete, Solving parabolic and hyperbolic equations by Generalized Finite
Difference Method. Journal of Computational and Applied Mathematics, Vol 209, Issue 2, 15 De-
cember 2007, Pages 208-233.
[4] A.R. Mitchell, D.F. Griffiths, The Finite Difference Method in Partial Differential Equations. Inter-
national Journal for Numerical Methods in Engineering (1980).
[5] F. Ureña, J.J. Benito, L. Gavete, R. Alvarez, Resolución de ecuaciones diferenciales en derivadas
parciales dependientes del tiempo de segundo orden utilizando Diferencias Finitas Generalizadas.
Revista Internacional de Métodos Numéricos para cálculo y diseño en ingenierı́a. Vol. 19, 3, 331-340
(2003).