Elementos Finitos en \(d=1\)

Métodos Numéricos Avanzados. 2023-2024

Author
Affiliation

David Gómez-Castro

Universidad Complutense de Madrid

Formulación débil

Integración por partes

Recordamos que \[ \frac{\partial}{\partial x} \left( v \frac{\partial u}{\partial x} \right) = \frac{\partial v}{\partial x} \frac{\partial u}{\partial x} + v \frac{\partial^2 u}{\partial x^2} \]

Integrando \[ v(1) \frac{\partial u}{\partial x}(t,1) - v(-1) \frac{\partial u}{\partial x}(t,-1) = \int_{-1}^1 \frac{\partial v}{\partial x} (x) \frac{\partial u}{\partial x}(t,x) \diff x + \int_{-1}^1 v(x) \frac{\partial^2 u}{\partial x^2} (t,x) \diff x \]

Un momento de teoría

Si \(v\) es una función Lipschitz, entonces es derivable en casi todo punto. Y denotamos \(\frac{\diff v}{\diff x}\).

Esta derivada (a veces llamada “débil”) sigue teniendo la propiedad de integración por partes: si \(u\) es de clase \(C^2\) entonces \[ v(1) \frac{\partial u}{\partial x}(t,1) - v(-1) \frac{\partial u}{\partial x}(t,-1) = \int_{-1}^1 \frac{\partial v}{\partial x} (x) \frac{\partial u}{\partial x}(t,x) \diff x + \int_{-1}^1 v(x) \frac{\partial^2 u}{\partial x^2} (t,x) \diff x \]

El espacio de funciones Lipschitz es un ejemplo de los llamados espacios de Sobolev. Denotaremos \[ W^{1,\infty}_0 (-1,1) = \{ v \in C([-1,1]) \text{ Lipschitz tal que } v(-1) = v(1) = 0 \}. \]

Ejemplo. Ecuación de calor con condiciones Dirichlet

Volvemos sobre la ecuación de calor \[ \tag{C} \label{eq:calor} \frac{\partial u}{\partial t} = \Delta u + f(t,x) \]

Tomemos \(\eqref{eq:calor}\) en \((-1,1)\) con condiciones de Dirichlet homogéneas \(u(t,-1) = u(t,1) = 0\). La solución exacta de la ecuación satisface que para toda \(v\) Lipscthitz se tiene

\[ \int_{-1}^1 \frac{\partial u}{\partial t} v = \int_{-1}^1 (\Delta u) v + \int_{-1}^1 f v \]

Si \(v(\pm 1) = 0\) podemos integrar por partes llegamos a la llamada llamada formulación débil \[ \int_{-1}^1 \frac{\partial u}{\partial t} v = - \int_{-1}^1 \frac{\partial u}{\partial x} \frac{\partial v}{\partial x} + \int_{-1}^1 f v , \qquad \forall v \in W_0^{1,\infty} (\Omega). \]

Ejemplo. Ecuación de calor con condiciones Neumann

Tomemos \(\eqref{eq:calor}\) en \((-1,1)\) con condiciones de Neumann homogéneas \(\frac{\partial u}{\partial x}(t,-1) = \frac{\partial u}{\partial x}(t,1) = 0\). La solución exacta de la ecuación satisface que para toda \(v\) Lipscthitz se tiene

\[ \int_{-1}^1 \frac{\partial u}{\partial t} v = \int_{-1}^1 (\Delta u) v + \int_{-1}^1 f v \]

Podemos integrar por partes incluso si \(v(\pm 1) \ne 0\) a la llamada llamada formulación débil \[ \int_{-1}^1 \frac{\partial u}{\partial t} v = - \int_{-1}^1 \frac{\partial u}{\partial x} \frac{\partial v}{\partial x} + \int_{-1}^1 f v , \qquad \forall v \in W^{1,\infty} (\Omega). \]

Método de Galerkin

Un problema en dimensión finita

Elegimos funciones de ``base’’ \[ \varphi^h_1, \cdots, \varphi^h_N \]

y queremos buscar una aproximación \(u(t,x) \sim u^{h} (t,x)\) donde \[ \tag{U} \label{eq:uh} u^{h} (t,x)= \sum_{j=1}^N u_j^h (t) \varphi^h_j (x). \]

Estas funciones \(\varphi^h_i\) generan un sub-espacio vectorial en el espacio de funciones, que llamamos \(V^h\). Se busca \[ u^h (t) \in V^h = \mathrm{span} \{ \varphi^h_1, \cdots, \varphi^h_N \}. \]

Sin pérdida de generalidad consideramos que son linealmente independiente, eliminando en las que sean necesario para construir un base, y cambiándolo los índices.

Ejemplo. Ecuación de calor con condiciones Dirichlet

Volvemos sobre \(\eqref{eq:calor}\) en \(t > 0\) y \(x \in (-1,1)\) con condiciones de Dirichlet homogéneas \(u(t,-1) = u(t,1) = 0\).

La solución exacta es la única curva \(t \mapsto u(t) \in C^2([-1,1] ) \cap C_0([-1,1])\) tal que \[ \int_{-1}^1 \frac{\partial u}{\partial t} v = - \int_{-1}^1 \frac{\partial u}{\partial x} \frac{\partial v}{\partial x} + \int_{-1}^1 f v , \qquad \forall v \in W_0^{1,\infty}([-1,1]) \]

Como aproximación numérica buscamos un espacio \(V^h \subset W_0^{1,\infty} (-1,1)\) y la única curva \(t \mapsto u^h(t) \in V^h\) tal que \[ \int_{-1}^1 \frac{\partial u^h}{\partial t} v = - \int_{-1}^1 \frac{\partial u^h}{\partial x} \frac{\partial v}{\partial x} + \int_{-1}^1 f v \diff x , \qquad \forall v \in V^h. \]

Como \(V^h\) es un espacio de dimensión finita, esto es equivalente a que

\[ \int_{-1}^1 \frac{\partial u^h}{\partial t} \varphi_i^h = - \int_{-1}^1 \frac{\partial u^h}{\partial x} \frac{\partial \varphi_i^h}{\partial x} + \int_{-1}^1 f \varphi_i^h \diff x , \qquad \forall i = 1, \cdots, N. \]

Volviendo sobre \(\eqref{eq:uh}\) tenemos \[ \sum_{j=1}^N \frac{\diff u_j^h}{\diff t} (t) \int_{-1}^1 \varphi_j^h \varphi_i^h = - \sum_{j=1}^N u_j^h(t) \int_{-1}^1 \frac{\partial \varphi_j^h}{\partial x} \frac{\partial \varphi_i^h}{\partial x} + \int_{-1}^1 f(t,x) \varphi_i^h (x) \diff x. \]

Llegamos a un sistema de la forma \[ M^h \frac{\diff }{\diff t} \begin{pmatrix} u_1^h \\ \vdots \\ u^N_h \end{pmatrix} = - A^h \begin{pmatrix} u_1^h \\ \vdots \\ u^N_h \end{pmatrix} + b^h (t) \]

donde \[ M^h_{ij} = \int_{-1}^1 \varphi_j^h \varphi_i^h , \qquad A_{ij}^h = \int_{-1}^1 \frac{\diff \varphi_j^h}{\diff x} \frac{\diff \varphi_i^h}{\diff x}, \qquad b_i^h (t) = \int_{-1}^1 f(t,x) \varphi_i^h (x) \diff x . \]

Dado que los elementos de la base son linealmente independiente, la matriz \(M^h\) es no singular. Por tanto, podemos escribir \[ \frac{\diff }{\diff t} \begin{pmatrix} u_1^h \\ \vdots \\ u^N_h \end{pmatrix} = - (M^h)^{-1} A^h \begin{pmatrix} u_1^h \\ \vdots \\ u^N_h \end{pmatrix} + (M^h)^{-1} b^h(t) \] Este problema evolutivo tiene siempre solución.

Concluímos así la semi-discretización espacial. Vamos a estudiar las buenas propiedades que tiene este método, empezando por el problema estacionario asociado.

Ejemplo: ecuación de Laplace con condiciones Dirichlet

Consiste en encontrar soluciones estacionarias (e.d. independientes del tiempo) de la ecuación de calor, es decir buscamos un único vector \(u^h \in \mathbb R^N\) tal que \[ A^h \begin{pmatrix} u_1^h \\ \vdots \\ u^N_h \end{pmatrix} = b^h \]

Marco abstracto

Sea \(A: H \to H\) un operador lineal, donde \(H\) tiene un producto escalar \((\cdot, \cdot)_H\).

Para estudiar el problema evolutivo

\[ \frac{\diff}{\diff t} u = - A u \]

Siempre podemos escoger una base y escribir \(\eqref{eq:uh}\).

Repitiendo el razonamiento anterior, donde \(Au = - \Delta u - f\), tenemos \[ M^h \frac{\diff }{\diff t} \begin{pmatrix} u_1^h \\ \vdots \\ u^N_h \end{pmatrix} = - A^h \begin{pmatrix} u_1^h \\ \vdots \\ u^N_h \end{pmatrix} \]

La matriz \(M^h\) se suele llamar matriz de masa, y \(A^h\) se llama matriz de rigidez. Esto tiene que ver con los problemas de cuarto orden que aparecen en elasticidad.

donde \[ M^h_{ij} = (\varphi_j^h, \varphi_i^h)_H , \qquad A_{ij}^h = - (A \varphi_j, \varphi_i)_H. \]

Observación Se puede extender esta formulación a un problema no-lineal.

Método de elementos finitos

El método de elementos finitos consiste en:

  • Dividir el dominio en poliedros (es decir, crear una malla)

  • Tomar como funciones de base funciones que son polinomios de un grado \(g\) dado en cada elemento

Para esta sección recomendamos la bibliografía: Ramos [2012], Süli [2020]

En dimensión 1

En dimensión 1, esto significa tomar puntos \[ a = x_0^h < x_1^h < \dots < x_{M+1}^h = b. \]

Sobre esta malla se toma el espacio de las funciones de polinomios de grado \(g\) a trozos.

Si introducimos condiciones de tipo Dirichlet, entonces añadimos esta condiciones al espacio: \[ V^h = \{ v^h \in C_0^p ([-1,1]) : v^h|_{[x_i, x_{i+1}]} \text{ es polinomio de grado } g \} \]

En las proximas secciones tomaremos \(p = 0\) (es decir funciones continuas, pero no necesariamente diferenciables), \(g = 1\) (es decir lineal en cada intervalo).

La base de funciones de funciones lineales a trozos

Es muy habitual tomar \(g = 1\), es decir funciones continuas y lineales a trozos.

Como las funciones son lineales, para determinar una de ellas nos basta con determinar los valores \(v^h(x_i)\).

Elegimos la base de funciones tales que \(\varphi_i(x_j) = \delta_{ij}\).

En tal caso, se cogen las funciones de base lineales a trozos \[ \varphi_i (x) = \begin{dcases} \left( \frac{x-x_{i-1}}{x_i - x_{i-1}} \right)_+ & x \le x_i , \\ \left( \frac{x_{i+1}-x}{x_{i+1} - x_i} \right)_+ & x > x_i , \\ \end{dcases}, \]

Matriz de rigidez para el Laplaciano con condiciones Dirichlet

Como las condiciones son Dirichlet podemos tomar el espacio \[ V^h = \mathrm{span} \{ \varphi_1, \cdots, \varphi_M \}. \]

Si la malla es uniforme \(x_{i+1} - x_i = \Delta x\) y tomamos condiciones Dichlet \[ \int_{-1}^1 \frac{\partial \varphi_i^h}{\partial x} \frac{\partial \varphi_j^h}{\partial x} = \frac 1 {\Delta x} \begin{cases} 2 & i = j \\ {-1} & |i-j| = 1, \\ 0 & |i-j| > 1. \end{cases} \]

Acabamos con la famosa matriz para el problema de Dichilet \[ A^h = \frac{1}{\Delta x} \begin{pmatrix} 2 & - 1 & \\ -1 & 2 & -1 \\ && \ddots \\ && -1 & 2 \end{pmatrix} \]

No calcularemos analíticamente la matriz \(M^h\) ahora.

Observación. En el método de diferencias finitas la matriz tiene factor \((\Delta x)^2\). Pero el vector \[ b_i^h(t) = \int_{-1}^1 f(t,x) \varphi_i^h(x) \diff x \sim \Delta x \]

Matriz de rigidez para el Laplaciano con condiciones Neumann

Como las condiciones son Neumann podemos tomar el espacio \[ V^h = \mathrm{span} \{ \varphi_0, \varphi_1, \cdots, \varphi_M, \varphi_{M+1} \}. \]

Para los puntos interiores los la integrales son iguales. Para los puntos de borde \[ \int_{-1}^1 \frac{\partial \varphi_0}{\partial x}\frac{\partial \varphi_0}{\partial x} = \frac 1 {\Delta x} \qquad \int_{-1}^1 \frac{\partial \varphi_0}{\partial x}\frac{\partial \varphi_1}{\partial x} = \frac {-1} {\Delta x} \]

Se tiene la matriz \[ A^h = \frac{1}{\Delta x} \begin{pmatrix} 1 & - 1 & \\ -1 & 2 & -1 \\ && \ddots \\ && -1 & 1 \end{pmatrix} \]

Cálculo simbólico

El paquete SymPy.jl

Utilizaremos una librería de Python sympy. Debemos tener instalado python en el sistema

En MacOS / Linux ejecutar en terminal

$ pip3 install sympy

Y luego en la interfaz en Julia

using Pkg; Pkg.add("SymPy")

En Windows si Python está instalado mediante Anaconda, el cálculo simbólico explicado abajo sólo funcionará en libreta de Jupyter. En una libreta de Jupyter escribir

using Pkg; Pkg.add("SymPy")

Función lineal a trozos para puntos equiespaciados

Son un desplazamiento de

using SymPy, Plots, LaTeXStrings

x  = symbols("x", real=true)

Φ = sympy.Piecewise( 
        ( 1 + x, And(Ge(x,-1),Le(x, 0))) , 
        ( 1 - x, And(Le(x,1),Gt(x, 0))), 
        (0, True))

display(Φ)

plot(Φ,xlim=(-2.0, 2.0), label=L"\Phi", thickness_scaling=1.5, linewidth=2)

\(\begin{cases} x + 1 & \text{for}\: x \geq -1 \wedge x \leq 0 \\1 - x & \text{for}\: x \leq 1 \wedge x > 0 \\0 & \text{otherwise} \end{cases}\)

Problema con condiciones Dirichlet

Funciones de base lineales a trozos con condiciones de Dirichlet

Tomamos la malla \[ x_0 = -1 < x_1 = -1+\Delta x < \cdots < x_{M+1} = 1 \] Las condiciones de Dirichlet homogéneas nos dicen que los valores \(u^h(t,\pm 1)=0\), luego las incógnitas son \(u_1^h, \cdots, u_N^h\).

Construímos las funciones de base lineales a trozos en los puntos que son incógnitas

Δx = 2.0^(-2)
xx = (-1.0 + Δx):Δx:(1.0-Δx)
NN = length(xx);

φ = subs.(Φ , x ,  x/Δx .- xx/Δx )

display(φ[1]);

\(\begin{cases} 4.0 x + 4.0 & \text{for}\: 4.0 x + 3.0 \geq -1 \wedge 4.0 x + 3.0 \leq 0 \\- 4.0 x - 2.0 & \text{for}\: 4.0 x + 3.0 \leq 1 \wedge 4.0 x + 3.0 > 0 \\0 & \text{otherwise} \end{cases}\)

p = plot(xlim=(-1,1),thickness_scaling=1.5,legend=:outerright,
    xlabel=L"x",ylabel=L"u")
for i=1:length(xx)
    p=plot!(φ[i], label=latexstring("\\varphi_{$(i)}"),linewidth=2)
end
plot(p)