Discretización de ecuaciones diferenciales parciales

2 de Septiembre de 2020
Universidad Nacional de Colombia, sede Medellín
Juan Esteban González Sierra, jegonzals@[Link]

1. Introducción
En el siguiente trabajo se pretende dar una introducción al método de Discretización de las
Ecuaciones Diferenciales Parciales(EDP) para aproximar las soluciones a dichas ecuaciones. Para
ello se hará una breve descripción del proceso de Discretización, y a continuación se presentaran
3 métodos utilizados para aproximar las soluciones, que emplean Discretización. Los métodos son
Diferencias Finitas, Elementos Finitos y de Volúmenes Finitos. El trabajo estará apoyado en el
software matemático MATLAB, y anexaran las rutinas que acompañaran a cada método. Por ultimo
en las conclusiones se darán algunas pautas que ayudaran al lector a tomar la decisión sobre que
método emplear, según las características del problema a resolver.

2. Discretización
La forma mas común de afrontar un problema que se describe mediante EDP, es discretizarlas
dichas ecuaciones, es decir, aproximar las ecuaciones mediante un número nito de incógnitas. La
discretización es un proceso matemático mediante el cual vamos a obtener resultados aproximados
de la ecuación diferencial del problema, mediante la creación de divisiones en el dominio, para
generar sub-intervalos, en los cuales se cambiara las EDP por funciones que ayudaran aproximar las
soluciones. Esta división es mejor conocida como malla, cada método que emplea la Discretización,
tiene una forma particular de malla, que ayuda optimizar su aproximación a la solución.

3. Método Diferencias Finititas(MDF)

Este método fue implementado por el matemático Brook Taylor, en su obra Methodus Incre-
mentorum Directa et Inversa (Londres, 1715), y consiste en que se sustituyen las derivadas por
diferencias. Es un método nodal, es decir depende de una maya.
Supongamos que tenemos la EDP:
∂u ∂u
= (1)
∂t ∂x
El primer paso para realizar la Discretización de las EDP por MDF, es realizar las mallas,
Consiste en seccionar el dominio en n, m puntos en los ejes x y el eje y respectivamente, como se
muestra en la gura 1.
Sea Vij la aproximación de u en el punto de la malla (xi ; tj ), y denimos la variable discreta
vectorial V j = [V1j , V2j , ..., Vij ] para nombrar el vector de aproximaciones en el tiempo j . Para MDF
existen varios métodos, que a continuación veremos:
MDF explícito o regresivo

Este método como su nombre lo indica, consiste en que para aproximar Vij+1 se utiliza los términos
anteriores de la siguiente manera:
 (a) n=0,m=0 (b) n=3,m=5 (c) n=9,m=5

Figura 1: Mayas Rectangulares

j
Vij+1 −Vij Vi+1 −2Vij +Vi−1
j

k = h2

donde h = n+1
b−a
con a y b condiciones iniciales y n es el número de sub-intervalos. Luego de operar
llegamos a la siguiente formula:
Vij+1 = Vij + r(Vi+1
j
− 2Vij + Vi−1
j
) (2)
Por lo que para hallar V j+1 se debe resolver un sistema de n ecuaciones, o expresado de manera
matricial, con A = trid(r; 1 − 2r; r) la matriz tridiagonal de orden n con elementos diagonales 1 − 2r
y elementos super y sub diagonales iguales a r. Tenemos que:
V j+1 = AV j (3)
lo que indica que
V j = Aj V 0 (4)
Los iterados consecutivos se consiguen por aplicación de potencias de A a la condición inicial.
MDF Implícito o progresivo

Este método como su nombre lo indica, consiste en que para aproximar Vij+1 se utiliza los términos
anteriores y el mismo Vij+1 de la siguiente manera:

j+1
Vij+1 −Vij Vi+1 −2Vij +Vi−1
j+1

k = h2

Luego de operar llegamos a la siguiente formula:

Vij = (1 − 2r)Vij+1 − r(Vi+1
j+1 j+1
+ Vi−1 ) (5)
Por lo que para hallar V j+1 se debe resolver un sistema de n ecuaciones, o expresado de manera
matricial, con A = trid(−r; 1 + 2r; −r) la matriz tridiagonal de orden n con elementos diagonales
1 + 2r y elementos super y sub diagonales iguales a −r. Tenemos que:

AV j+1 = V j (6)
En cada j , se debe resolver un sistema n de ecuaciones lineales con matriz simétrica, tridiagonal y
diagonalmente dominante (por tanto no singular).
MDF Combinado o Mixto

2
 (a) n=0 (b) n=1 (c) n=39

Figura 2: Solución Anlitica vs Solución MDF

Este método como su nombre lo indica, consiste en combinar los dos métodos anteriores mediante
un parametro α ∈ [0, 1] de la siguiente manera:

j+1
Vij+1 −Vij Vi+1 −2Vij +Vi−1
j+1 j
Vi+1 −2Vij +Vi−1
j

k =α h2 + (1 − α) h2

Cuando el α = 1
2 e se llama de Crank-Nicolson y esta dado por la igualdad
j+1
− rVi−1 + 2(1 + r)Vij+1 − rVi+1
j+1 j
= rVi−1 + 2(1 − r)Vij + rVi+1
j
(7)
Por lo que para hallar V j+1 se debe resolver un sistema de n ecuaciones, o expresado de manera
matricial, con A = trid(−r; 2(1 + r); −r) y B = trid(r; 2(1 − r); r) matrices tridiagonal de orden n
con elementos diagonales. Tenemos que:
AV j+1 = BV j (8)
Para hallar V j+1 este es el sistema que debe resolverse para cada iteración temporal.

Veamos un ejemplo en 1-D y otro en 2-D.

Sea Ut (x, t) = cUxx , en la R = {(x, t) : 0 ≤ x ≤ 1, 0 ≤ t ≤ 0,5}, f (x, t) = sin(πx)(1 + 2cos(πx))

La Solucción analitíca de la EDP esta dada por sin(πx)e−tπ + sin(2πx)e−4tπ , Para resolver la
2 2

EDP se utilizo la función Adi_CP_1D, en la gura 1 se observa la el gráco con la solución analítica
vs la aproximación de la solución mediante MDF, para tres n = 1, 2, 39. Se observa que para n = 0
las gracas son muy cercanas, y a medida que se va iterando, la aproximación MDF, se aleja de la
solución analítica, para n = 39 se observa una diferencia de 10−3 , que equivale al error.

Ahora veamos un ejemplo 2-D. Sea la EDP:

∂2u ∂2u
U= ∂x2
+ ∂y 2

Condiciones de frontera Dirichlet homogeneas y condiciones iniciales de 4C para t = 0, y tomaremos

26 nodos en la maya, y n = 25.

Para resolver la EDP se utilizo la función aditemp, en la gura 2 se observa la el graco la

aproximación de la solución mediante MDF, para tres n = 0, 4, 25.

3
 (a) n=0 (b) n=4 (c) n=25

Figura 3: Solución Solución MDF 2-D

4. Método Elementos Finitos(MEF)

El método de elementos nitos se empezó a desarrollar en la década de 1940, con los trabajos
de Hrenniko en 1941 y McHenry en 1943, en las soluciones de tensiones en solidos continuos, mas
adelante se introdujo la interpolación por partes sobre las funciones triangulares en subregiones
(Elementos) nitas que componen el conjunto de la región como un método para obtener soluciones
aproximadas.

Para resolver una EDP por medio de MEF se siguen los siguientes pasos:

1. Reformulación del problema en forma variacional.

2. Dividir el dominio de las variables independientes en subregiones, llamados elementos, asocia-
dos a un espacio vectorial de dimensión nita, estos elementos pueden ser triángulos, rectán-
gulos, pentágonos, etc.
3. Obtener la proyección del problema original sobre el espacio de elementos nitos, lo que lo
convierte en un sistema de ecuaciones lineales de grandes dimensiones que se resuelve.
4. Se resuelve el sistema de ecualciones.

Veamos que como se desarrolla estos pasos de manera general, para un mayor detalle en la de-
ducción se puede consultar 1 y 2 .

Supongamos que tenemos la EDP parabólica del calor:

ut (x, t) − cuxx (x, t) = f (x) 0 < x < a, 0<t≤b (9)
Condición de frontera de Dirichlet y condición inicial u(x, 0) = g(x) , 0 < x < a , t = 0.

El primer paso es reformulación del problema en forma variacional, para ello tomamos:
um −um−1
ut =
R a4t
F(v) = cA(u, v) + 0 um −u
4t
m−1
v(x)dx
Ra
M(um−i , v) = 0 um−i v(x)dx con i = 0, 1, 2...

La forma variacional queda:

1
Yousef Saad, Iterative Methods for Sparse Linear Systems, capitulo 2 sección 3
2
M. G. Larson and F. Bengzon, The nite element method:theory, capitulo 2

4
 (a) Solución Analitíca (b) Solución Aproximada (c) n=0

Figura 4: Solución Anlitica vs Solución MEF

1 1
, ∀v ∈ V, um− , um−1 ∈ U


 4t M(um , v) + cA(um , v) = F(v) + 4t M(um−1 , v)



u(x, 0) = g(x) , para0 < x < a




0<t≤b

u(t, 0) ,
Para hallar V j+1 este es el sistema que debe resolverse para cada iteración temporal.

Veamos un ejemplo en 1-D

Sea ut (x, t) − uxx = f (x), en la R = {(x, t) : 0 ≤ x ≤ 1, 0 ≤ t ≤ 1},

f (x, t) = sin(πx)(1 + 2cos(πx))
La Solucción analitíca de la EDP esta dada por sin(πx)e−tπ + sin(2πx)e−4tπ , Para resolver
2 2

la EDP se utilizo la función Mef_calor_1D, en la gura 2 se observa la el gráco con la solución

analítica vs la aproximación de la solución mediante MEF. La solución analítica tiene una maya de
[19x19] en cambio la solución mediante MEF solido un malla de [19x6].

5. Método Volúmenes Finitos(MEF)

El método de volúmenes nitos (MVF) fue desarrollado en la década de 1970 por McDonald,
MacCormack, y Paullay, su principal uso es para resolver problemas de mecánica de uidos y de las
leyes de conservación.

Para resolver una EDP por medio de MVF se siguen los siguientes pasos:
1. Dividir el dominio de las variables independientes en subregiones con volúmenes de control, a
diferencia de los métodos anteriores, los vértices de las subregiones no son los nodos, los nodos
de están dentro de cada subregion, como se observa en la gura 5, a la izquierda se observa la
malla para el MDF y a la derecha la malla para el MVF.
2. Formular las ecuaciones integrales de conservación para cada volumen de control.
3. Hacer la aproximación de las integrales, mediante algún método numérico.
4. Aproximar los valores de las variables en las caras y las derivadas con la información de las
variables nodales.

5
 Figura 5: Malla MDF vs Malla MVF

5. resolver el sistema algebraico obtenido.

A continuación veamos un ejemplo de la implementación de este método tomado de 3 .

5.1. Burgers Equation

Es una ecuación escalar parabólica:

∂w ∂w ∂2w
+w −v 2 =0 (10)
∂t ∂x ∂x
donde v es el coeciente de viscosidad. Se centró principalmente en la ecuación obtenida como límite
cuando te tiende a cero, que, en forma conservadora, es la siguiente:
w2
 
∂w ∂
+ =0 (11)
∂t ∂x 2

En la gura 6 se puede observar la particularidad de las soluciones de la ecuación sistema hiper-

bólico lineal en una variable espacial, en un tiempo k la solucion rompe contra un eje.
Las siguientes ecuaciones son las soluciones aproximadas asociada con la ecuación de Burgers:

n+1 n 4t n n n ) , si w n ≥ 0
 wj = wj −
 4x wj (wj − wj−1 j

 wn+1 = wn − 4t n n − wjn ) , si wjn ≤ 0


j j 4x wj (wj+1

Existen cuatro Esquemas conservadores para aproximar las soluciones:

Godunov
Lax Friedrichs
Q-scheme
Non conservative

Estos esquemas se pueden consultar en 4 , y la comparación de estos se realiza a travez de la

función BEconservation.
3
Maria Elena Vazquez Cendon, Hyperbolic Conservation Laws. Basic Concepts and Examples, capitulo 2 pagina
20
4
M. G. Larson and F. Bengzon, The nite element method:theory, capitulo 6

6
 (a) n=0 (b) n=3 (c) n=k

Figura 6: Burgers Equation

6. Conclusiones
El MDF es fácil de comprender y de programar, su uso se restringe a problemas con regiones
geométricas simples. Su estabilidad depende de los valores de los parámetros, sobre todo en el
MDF explícito.
El MEF En la actualidad es el mas recomendado para los diversos problemas, ya que su
versatilidad a la hora de implementarse, su en mallado es muy exible y se puede adaptar a
problemas con regiones de mas complejas.
Los MEF son ampliamente utilizados por la comunidad de análisis numérico para estudiar
métodos numéricos para el ujo de uidos. Existe una gran cantidad de trabajos descritos en
publicaciones cientícas sobre CFD y MEF, así como trabajos más recientes para los métodos
nodales discontinuos de Galerkin (DG), que son métodos de elementos nitos con funciones de
base discontinua.
MVF puede ser más exacto que el MEF, pero al tener mallas jas, en algunos problemas sobre
tensión el material pasa atravez de estas mallas, aumentando el error.
No hay un unico criterio para escoger MEF ó MVF depende del problema y la preferencia del
usuario.

Referencias
[1] Larson, M. and F. Bengzon, The nite element method: theory, implementation and applications,
Springer, 2013, capitulo 2.
[2] Randall J. LeVeque, Numerical Methods for Conservation Laws, capitulos 1 al 4.
[3] Maria Elena Vazquez Cendon, Hyperbolic Conservation Laws. Basic Concepts and Examples.
[4] Carlos Enrique Mejía Salazar, Paper Invitación al Análisis Numérico, capitulo 7.