Trabajo final de Análisis Numérico

Sobre el método de los elementos finitos

en una dimensión
Juan Camilo Garzón Upegui
Profesor: Mauricio A. Osorio

Universidad Nacional de Colombia sede Medellı́n

Escuela de Matemáticas
Mayo 28 de 2016

1. Introducción
La ecuación de Stokes (también llamada Stokes estacionaria) es un sistema de ecua-
ciones diferenciales parciales lineales derivado de las ecuaciones de Navier-Stokes, las
cuales pretenden determinar el flujo de un fluido viscoso (denominado fluido Newtoniano)
cuando el número de Reynolds es muy bajo, aplicando los principios de la mecánica y la
termodinámica a un volumen fluido.
(El número de Reynolds (Re) es un número adimensional utilizado en mecánica de flui-
dos, diseño de reactores y fenómenos de transporte, para caracterizar el movimiento de
un fluido, este relaciona la densidad, viscosidad, velocidad y dimensión tı́pica de un flujo
en una expresión adimensional, que interviene en numerosos problemas de dinámica de
fluidos).
Esta ecuación se presenta en la atmósfera, corrientes de agua, hidráulica, etc; su relación
con las ecuaciones de Navier-Stokes se basa en que la ecuación de Stokes es una lineali-
zación estacionaria de estas.
Para esta ecuación no se dispone de una solución analı́tica general, por lo que recurrimos
al “Método de Elementos Finitos (MEF)” para encontrar una aproximación de la solu-
ción.
Para el desarrollo de MEF, se presenta la ecuación de Stokes en su formulación fuerte
dada por:

Encontrar ~u = (u1 , u2 ) y p tal que :



 −µ∆~u + ∇p = f~ para x ∈ Ω

div~ u = 0 para x ∈ Ω



R
Ω
p = 0 condicion de solubilidad
~u(x) = h~D (x) para x ∈ ΓD ⊂ ∂Ω





∇~u(x) · ~η = h~ (x) para x ∈ Γ ⊂ ∂Ω

N N

1
 Donde x = (x1 , x2 ), Ω ⊂ R2 , ~u = (u1 , u2 ) es la velocidad, p es la presión, µ
es el coeficiente de viscosidad del fluido, y f es una fuerza externa que afecta al fluido.
La condición de frontera la separamos en dos partes, la condición de frontera de Dirichlet
denotada ΓN , y la condición de frontera de Neumann denotada ΓD , donde ∂Ω = ΓD ∪ΓN .

Este es un ejemplo clásico en la actualidad de la necesidad de un metodo numerico como

el MEF para resolver problemas cotidianos en ciencias e ingenierı́a.
Otro ejemplo es la ecuación de advección-difusión.

La ecuación de advección-difusión, es una EDP de tipo parabólico, y describe el pro-

ceso de la advección y la difusión de una sustancia en un fluido, donde la advección es
una especie de transporte de la sustancia, es decir, la sustancia es transportada por el efec-
to de un campo vectorial, generado por una energı́a, por ejemplo, el transporte de una
sustancia en un rı́o, y la difusión es que tanto se propaga la sustancia en el medio que la
rodea.
Esta ecuación aparece en el estudio de la propagación de sustancias contaminantes en di-
ferentes medios, su solución representa la concentración de la sustancia que se dispersa
y es una información muy importante en el estudio de diversos problemas de diferente
naturaleza (fı́sica, quı́mica, económica, ambiental, etc.).
Debido a la complejidad de estos problemas no es posible en general obtener una solución
analı́tica para la ecuación de advección-difusión, por lo que usamos el MEF para hallar
una solución aproximada.
La ecuación de advección-difusión esta dada por:

Encontrar u : Ω ⊂ R2 → R tal que :

 ∂u ∂u ∂u

 ∂t
− µ∆u + β1 ∂x 1
+ β2 ∂x 2
= f para x ∈ Ω

u(x, t) = g (x) para x ∈ Γ ⊂ ∂Ω, t > 0
D D


 ∇u(x, t)η(x) = gN (x) para x ∈ ΓN ⊂ ∂Ω, t > 0
u(0, x) = u0 (x) para x ∈ Ω


∂u
donde x = (x1 ; x2 ), ∂t es el diferencial a través del tiempo, ∆u es el termino co-
rrespondiente a la difusión, con coeficiente de difusión µ > 0 y las primeras derivadas
parciales corresponden a la advección en la dirección β = (β1 , β2 ).
Para la aplicación del MEF a la ecuación de advección-difusión, se denomina a la defini-
ción anterior como la formulación fuerte de ésta, se halla su respectiva formulación debil,
y ası́ se sigue con el método para llegar a la solución de la ecuación. Podemos ver que
esta ecuación depende de el tiempo y del espacio.
Estos ejemplos nos muestran la importancia y necesidad de poder recurrir a los méto-
dos numéricos, tanto dentro de la matemática, como en la fı́sica, quı́mica, ingenierı́a y la
ciencia en general. Veamos ahora en más detalle el método de elementos finitos MEF.

2
2. Método de los elementos finitos
Al estudiar ecuaciones diferenciales parciales, existen distintos métodos para resol-
verlas numéricamente, pero en este documento nos enfocaremos en el “método de los
elementos finitos (MEF)”, el cual, en general, requiere de tres pasos:

1) Debemos construir una “formulación débil (o variacional)” de la ecuación diferencial

en estudio, esto implica que, la ecuación debe ser planteada en un espacio de funciones
adecuado, que frecuentemente es un espacio de Hilbert, y luego verificar la existencia de
la solución de la formulación débil para este problema.

2) Luego de escribir el problema en un espacio de funciones adecuado, debemos introdu-

cir espacios de elementos finitos, es decir, debemos aproximar el problema obtenido en
el primer paso a un espacio de dimension finita, de aquı́, surge un sistema linear, ya que
introducimos unas bases para hallar una aproximación de elementos finitos.

3) Por último, debemos resolver el sistema lineal obtenido en el paso 2, el cual se so-
luciona por un método directo o iterativo.

2.1. Conceptos preliminares

Para obtener una forma precisa de una formulación débil de alguna ecuación diferen-
cial parcial (EDP), o para estudiar el comportamiento del MEF, es necesario introducir
algunos espacios de funciones. Para esto, recordemos algunas definiciones y resultados.

Definición. Dado un espacio vectorial V, una función k·k : V → R+ ∪ {0} define una
norma, si satisface las siguientes propiedades:
i) Dado v ∈ V, kvk = 0 si y solamente si v = 0.
ii) Para todo α ∈ R, y v ∈ V se tiene que kαvk = |α| kvk.
iii) Dados u, v ∈ V se cumple, ku + vk ≤ kuk + kvk.
Y ası́ decimos que el par (V, k·k) forma un espacio vectorial normado.

Definición. Decimos que un espacio vectorial normado (V, k·k) es completo, si toda su-
cesión de Cauchy en V es convergente.

Teorema. Sea (V, k·k) espacio vectorial normado, entonces existe un único espacio vec-
torial completo (H, k·k) tal que:
i) V ⊂ H.
ii) Dado cualquier elemento v ∈ H, existe una sucesión {vn }∞
n=1 ⊂ V , tal que
lı́mn→∞ vn = v.
En este caso decimos que V es denso en H; o en otras palabras que H es el cierre o com-
pleción de V .

3
 Usando estas herramientas, podemos construir los siguientes espacios, denominados
Espacios de funciones tipo Sobolev.

Definición. Definimos los espacios C(Ω), C ∞ (Ω) y C0∞ (Ω) como:

i) C(Ω) el conjunto de las funciones reales continuas con dominio Ω.
ii) C ∞ (Ω) el conjunto de las funciones continuamente diferenciables, es decir, si existen
todas sus derivadas en Ω.
iii) C0∞ (Ω) = {v ∈ C ∞ (Ω) : v = 0 en Γ ⊂ ∂Ω} .

Ahora, con los resultados anteriores, podemos definir el espacio L2 ası́:

La norma L2 de una función v : Ω → R esta dada por

Z 1
2
2
kvk0 = ( v(x) dx) .
Ω
El espacio L2 (Ω) esta dado por el compleción del espacio C(Ω) con relación a la
norma L2 , el espacio L2 (Ω) es el conjunto de las funciones cuadrado integrable, es decir
si v ∈ C ∞ (Ω), entonces
Z
2
v ∈ L (Ω) si y solo si v 2 < ∞.
Ω

Es decir, v es cuadrado integrable si la integral de la función al cuadrado es finita.

Los espacios H 1 y H01 los definimos asi:

La norma H 1 estará dada por
Z 1
2
2 2
kvk1 = ( v(x) + |∇v(x)| dx) .
Ω
1
Los espacios H (Ω) y H01 (Ω),son la compleción de los espacios C ∞ (Ω) y C0∞ (Ω) con
relación a la norma kvk1 . El espacio H 1 (Ω) es un espacio de Hilbert, apropiado para la
mayorı́a de problemas considerados en este documento. Note que,
H 1 (Ω) := v ∈ L2 (Ω) : v ′ ∈ L2 (Ω) .


También usaremos el espacio de funciones H01 (Ω), definido por,

H01 (Ω) := v ∈ H 1 (Ω) : v = 0 en Γ ⊂ ∂Ω ,


tal que, H01 (Ω) ⊂ H 1 (Ω).

Finalmente definimos la norma del espacio H 2 asi:

Z 1
2
2 2 2
kvk2 = ( v(x) + |∇v(x)| + |∆v(x)| dx) .
Ω
1
Donde el espacio H (Ω) es definido como la compleción del espacio C ∞ (Ω) con respecto
a la norma k·k2 . Decimos que H 2 etsa definido como,
H 2 (Ω) := v ∈ L2 (Ω) : v ′ , v ′′ ∈ L2 (Ω) .


4
 Podemos ver que,
H 2 (Ω) ⊂ H 1 (Ω) ⊂ L2 (Ω).
Ahora introduciremos V ⊂ Rd con d = 1; 2, el espacio de funciones lineales por partes,
el cual es uno de los ejemplos mas simples de espacios de elementos finitos. Primero,
debemos conocer la noción de Espacio de funciones de dimensión finita.

Definición. (Espacio de dimensión finita). Un espacio de funciones V tiene dimensión

finita n, si existe un conjunto de funciones φi : Ω → R, i ∈ {1, ..., n}, tal que cualquier
función ϕ ∈ V puede ser escrita como una única
P∞ combinación lineal de las funciones φi .
Esto es, existen n constantes αi tal que, ϕ = i=1 αi φi .

Consideremos el caso unidimensional, (d = 1), es decir cuando Ω ⊂ R.

Como Ω ⊂ R podemos suponer, sin perdida de generalidad que Ω = (a, b). Entenderemos
como una partición de Ω, una división de este dominio en subintervalos de la siguiente
forma S
[a, b] = N i=1 [zi−1 , zi ], donde z0 = a, zN = b y zi < zi+1 . Los valores zi son llamados
vértices de la partición.
Vamos a definir el espacio de las funciones lineales por partes sobre la partición introdu-
cida como, 
V = f ∈ C(Ω) : f |(zi ,zi+1 ) es lineal ,
Lo que quiere decir que, f |(zi ,zi+1 ) = ai x + bi , siendo ai y bi constantes apropiadas.
Podemos ver que el conjunto de funciones φi ∈ V, i ∈ {0, 1, ..., N } donde,
(
1 si i = j
φi (zj ) =
0 si i 6= j

forman una base para el espacio V . En otras palabras, para cualquier función v ∈ V puede
ser escrita como,
N
X
v(x) = c0 φ0 (x) + c1 φ1 (x) + ... + cN φN (x) = ci φi (x),
i=0

donde ci = v(xi ) con xi = zi vértices de la partición.

5
 Existen distintas condiciones de frontera, que se usan como su nombre lo indica, para
restringir el valor de una EDP en su frontera, de acuerdo a la situación del problema. Es-
tas condiciones se pueden introducir independientes una de la otra, o se pueden combinar,
pero en este documento solo trabajaremos con la condición de Dirichlet, y la condición
de Neumann, las cuales presentaremos a continuación.

Condición de Dirichlet
Esta condición indica que la solución, restringida a la frontera, es requerida con un valor
determinado; la enunciaremos para el caso de una dimension.
Sabemos que Ω = (a, b) ⊂ R para el caso de una dimension, y que su frontera son los
puntos a y b.
Encontrar u : (a, b) → R tal que
(
EDP en estudio, para x ∈ (a, b)
.
u(x) = g(x) para x = a, x = b

Condición de Neumann
En una dimensión tenemos que Ω = (a, b) ⊂ R, y que su frontera son los puntos a y b.

Encontrar u : (a, b) → R tal que



 EDP en estudio para x ∈ (a, b)

κ(x)u′ (x) = g(x) para x = a, x = b
Rb .

 g(a) − g(b) = a
f (condicion de computabilidad)
R b

a
u = 0 (condicion de solubilidad)

Algunas EDP’s que se usarán para introducir el MEF serán enunciadas a continuación.

EDP básica
Consideremos la siguiente EDP, en esta sección esta enunciada con la condición de Diri-
chlet, pero esta condición cambia de acuerdo al problema. La enunciaremos para el caso
de una dimension. En este caso tenemos que Ω = (a, b) ⊂ R, y consideramos la EDP
básica como,
Encontrar u : (a, b) → R tal que
(
−(κ(x)u′ (x))′ = f (x) para x ∈ (a, b)
(F ormulacion F uerte)
u(x) = g(x) para x = a, x = b
donde 0 < κmin ≤ κ(x) ≤ κmax para todo x ∈ (a, b).

EDP de Laplace
Esta EDP la usaremos, aludiendo a su fácil estudio, y comodidad para trabajar, acá defi-
niremos su denominada formulación fuerte, para el caso de una dimension, con condición
de Dirichlet, la cual puede variar de acuerdo al problema.

6
 En una dimensión tenemos que Ω = (a, b) ⊂ R, y consideramos la EDP de Laplace,

Encontrar u : (a, b) → R tal que

(
−u′′ (x) = f (x) para x ∈ (a, b)
(F ormulacion F uerte)
u(x) = g(x) para x = a, x = b

2.2. MEF en una dimensión

Veamos el MEF en una ecuación diferencial en una dimensión (tomaremos, sin perdi-
da de generalidad, el intervalo (a,b)), encontraremos los espacios adecuados para resolver
el problema, sus formulaciones, y miraremos algunos ejemplos.

Formulación Débil
Para deducir la formulación débil de una ecuación diferencial en el intervalo (a,b) tene-
mos que,
1) Suponer que existe u solución de la ecuación diferencial, multiplicar los dos lados de la
ecuación por una función de prueba v ∈ C0∞ (a, b) e integrar a ambos lados de la igualdad.

2) Usar la fórmula de integración por partes y las condiciones iniciales (de contorno)
para obtener expresiones que necesiten solamente de derivadas del menor orden posible.

3) Cambiar v ∈ C0∞ (a, b) por v ∈ V donde el espacio de funciones de prueba V sea

el más grande posible. También escoger U (espacio de funciones de forma), tal que la
solución u esté naturalmente en U .

Los espacios de funciones U y V más adecuados para estas condiciones son los “Es-
pacios de funciones tipo Sobolev” (Espacios de Hilbert) definidos en (a, b).
Los espacios U y V adecuados deben cumplir:

• La elección de estos espacios U y V es independiente una de otra, aunque por comodi-

dad, se puede tomar U = V , y su ventaja es que se obtienen sistemas lineales simétricos,
y de aquı́ podemos llegar, a calcular más fácil su solución.

• Las integrales que se obtuvieron deben estar definidas en estos espacios.

• Para toda u ∈ U se tiene que poder imponer la condición de contorno en estudio.

Para la elección de los espacios de funciones, podemos hacer Ω = (a, b) ⊂ R y tomamos

los espacios mencionados anteriormente; Por ejemplo, los espacios de Hilbert, H 1 (a, b),
H01 (a, b), H 2 (a, b). Al obtener dicha formulación débil, podemos garantizar la existencia
de la solución de esta formulación, que se conoce como “solución débil” de la ecuación
original. La solución del problema original se denomina “solución fuerte”, y en este caso
se cumplen las igualdades en el problema para todo x ∈ (a, b), y por lo tanto se pueden
calcular sus derivadas.

7
Obs: Toda solución fuerte de una ecuación es también solución débil. Si los coeficien-
tes de la ecuación son regulares y una solución débil es regular, entonces esta es solución
fuerte. Regular, quiere decir que la solución tiene tantas derivadas continuas como el or-
den de la ecuación diferencial.
Para continuar con la aplicacion del MEF, es necesario utilizar la formulación de
Galerkin, la cual omitiremos por extensa, pero se puede consultar en detalle, en [1] de las
referencias.

2.2.1. Ejemplo aplicativo: Ecuación de advección-difusión

Anteriormente, al inicio del documento, nos referimos al ejemplo de la ecuación de
advección-difusión, junto con el de la ecuación de Stokes.
Para aplicar el MEF al problema de advección-difusión, se comienza por hallar la formu-
lacion debil de la ecuación, definida en la introducción de éste documento.
comenzaremos suponiendo que u es una solución a la ecuación de advección-difusión, y
tomamos v ∈ H01 (Ω) como función de prueba.
Omitiendo la prueba, que no es muy larga, obtenemos la formulación débil de la ecuación
de advección-difusión que esta dada ası́,
R ∂u
v + µ Ω ∇u∇v + Ω β∇uv = Ω f v para toda v ∈ H01 (Ω)
R R R

 Ω ∂t

u(x, t) = g (x) para x ∈ Γ ⊂ ∂Ω, t > 0
D D
∇u(x, t)η(x) = gN (x) para x ∈ ΓN ⊂ ∂Ω, t > 0


u(0, x) = u0 (x) para x ∈ Ω


Ahora, con la ecuación anterior, debemos elegir un espacio de dimensión finita adecuado,
para construir la formulación de Galerkin, la cual es el paso a seguir en el MEF.
La formulación de Galerkin es un método para discretizar el espacio Ω en el que estamos
trabajando, pero el tiempo sigue siendo continuo. Luego de realizar esta discretización
del espacio, se puede discretizar el tiempo, usando la serie de Taylor de u.

Escribimos un código en Matlab para resolver la ecuación de advección-difusión en una

dimensión, en donde usamos la formulación matricial de esta, armando las matrices lo-
calmente, y luego ensamblando las matrices para obtener la aproximación. En Matlab,
creamos estas matrices con variables tipo estructura, las cuales asignan a cada elemento
de la partición Ki toda la información necesaria para construir la matriz.
Debemos saber que para este código trabajamos en el intervalo [0, 1], con el espacio de
funciones lineales por partes P10 (T h ), en vez de P20 (T h ), y que vamos a usar la condición
inicial u0 = 0. Ası́ que lo primero que hacemos es definir una serie de elementos que
vamos a usar a lo largo del código, entonces creamos las siguientes variables:

8
Donde N es el número de elementos de la partición, dt es el intervalo de tiempo, u0 es
la condición inicial, y los términos terminados en aux, son usados para la integración
numerica, por medio de la cuadratura de Gauss, la cual consiste en dar una aproximación
de una integral definida, de una función dada. Una cuadratura n de Gauss, selecciona unos
puntos de evaluación xi llamados puntos de cuadratura y unos coeficientes o pesos ωi para
i = 1, ..., n. El dominio de tal cuadratura es [−1, 1] y esta dada por,
Z 1 n
X
f (x)dx ≈ ωi f (xi ),
−1 i=1

pero como nuestras integrales deben ser calculadas en un dominio [a, b] = [0, 1] diferente
a [−1, 1] entonces debemos hacer el cambio,

b−a 1
Z b  
b−a a+b
Z
f (x)dx = f x+ dx
a 2 −1 2 2

y al aplicar cuadratura obtenemos,

b n  
b−aX b−a a+b
Z
f (x)dx = ωi f xi + ,
a 2 i=1 2 2

y esta es la aproximación que usamos en el código para calcular las integrales.

A continuación, se calcula la información correspondiente a los elementos finitos y las
funciones base de elementos finitos, ası́ como los datos necesarios para calcular las inte-
grales locales que aparecen en la formulación de Galerkin. Tenemos el código, a conti-
nuación, creamos las matrices locales, que en Matlab se le asigna elemento por elemento
la información necesaria para calcular las integrales, estas variables se les llama estructu-
ras, las cuales asignan a cada elemento Ki una parte de la información de la matriz. Con
el siguiente código.

9
 Ahora, creamos las matrices locales de M, A, V y el vector b, que en Matlab se le
asigna elemento por elemento la información necesaria. Para esto usamos el siguiente
código,

10
donde escribimos la integral para calcularla numéricamente localmente por el método de
la cuadratura de Gauss. Para este ejemplo tomamos los valores de parámetros µ = 0,05,
β = 0,5 y f = 0. Ahora para crear las matrices globales, creamos matrices dispersas
nulas, y por cada elemento vamos añadiendo a cada matriz su contribución, que se ve en
las siguientes lineas,

11
Con esto obtenemos las matrices globales que definen el sistema global y cuya solu-
ción aproxima la solución de la ecuación de advección-difusión. Aun falta discretizar el
tiempo, observe que la derivada temporal aparece en el argumento de la matriz M , de
la formulación matricial implı́cita, por lo que creamos una matriz nueva con el siguiente
código, y usando la condición de Dirichlet u(0) = 1 y u(1) = 0 la reescribimos:

Donde la última linea impone la condición de Dirichlet. Ahora, como tenemos condición
de Dirichlet solo necesitamos resolver para los grados de libertad interiores, y por tanto
solo necesitamos de la submatriz de grados de libertad interiores. Esto se extrae con,

Con esto, solo creamos una iteración, que reemplaza u0 por el último resultado para lograr
notar su desarrollo a través del tiempo, imprimiendo su resultado en una figura con las
siguientes lineas,

Con esto logramos la de la ecuación de advección-difusión por elementos finitos en una

dimension, las siguientes figuras F.1, F.2, F.3 nos muestran la evolución de la aproxima-
ción de la solución en los tiempos t = 0, 1, t = 1, 5 y t = 5 (tiempo final) respectivamente.

12
 F.1

F.1 es el plot de la aproximación de la ecuación de advección-difusión en el instante t =

∆t, con incremento de tiempo ∆t = 0,1. Para esto, usamos la formulación de Galerkin
cuando µ = 0,05, β = 0,5, f = 0.

F.2

F.2 es la aproximación de la ecuación de advección-difusión en el instante t = 1,5, con

incremento de tiempo ∆t = 0,1. Para esto, usamos la formulación de Galerkin cuando
µ = 0,05, β = 0,5 y f = 0.

13
 F.3

Con esto, tenemos una aproximación de la ecuación de advección-difusión, con t = 5,

coeficiente de difusión µ = 0,05 y velocidad de transporte β = 0,5, con intervalos de
tiempo ∆t = 0,1.
Concluimos entonces éste ejemplo ilustrativo de una aplicación directa del método de los
elementos finitos en una dimension, para ver un panorama mas general y detallado, ver
[1] de las referencias.

≪ Tomado de la tesis de grado;

Galeano V. Jonathan D.
Solución numérica de ecuaciones
diferenciales parciales parabólicas,
Universidad Nacional de Colombia, Bogotá. ≫

14
3. Referencias

[1] J. David Galeano V. (2013). Solución numérica de ecuaciones diferenciales par-

ciales parabólicas; método de elementos finitos aplicado a las ecuaciones de Stokes y de
advección-difusión . (Tesis de grado en matemáticas).
Universidad Nacional de Colombia, Bogotá.

[2] Claes Johnson. (1987). Numerical solution of partial diferential equations by the
finite element method;
Cambridge University Press, Cambridge, England.

15