Introducción. Métodos de Volúmenes Finitos y Elementos Finitos
Métodos Numéricos Avanzados. 2023-2024
David Gómez-Castro 
Modelización. Problemas en forma de divergencia
Ecuación de continuidad
\[ \newcommand{\bm}[1]{\mathbf{#1}} \newcommand{\diff}{\,\mathrm{d}} \DeclareMathOperator{\diver}{div} \newcommand{\flux}{{F}} \]
Algunos problemas clásicos en física se escriben en forma de divergencia.
Sea densidad \(\rho\) dea una cantidad que se transporta por un flujo \(\bm \flux\).
Si \(\omega \subset \mathbb R^d\) es un conjunto genérico (regular) \[ \frac{\diff}{\diff t} \int_\omega \rho(t,x) \diff x = - \int_{\partial \omega} \bm \flux (t,x) \cdot \bm n (x) \diff S(x) \]
Aplicando el teorema de la divergencia, tenemos \[ \frac{\diff}{\diff t} \int_\omega \rho(t,x) \diff x = - \int_{\omega} \diver \bm \flux(t,x) \diff x. \]
\[ \int_\omega \left( \frac{\partial \rho}{\partial t} (t,x) + \diver \bm \flux(t,x) \right) \diff x = 0 \]
Finalmente obtenemos la ecuación de continuidad
\[ \tag{C} \label{eq:ec continuidad} \frac{\partial \rho}{\partial t} (t,x) + \diver \bm \flux(t,x) = 0. \]
Ejemplos
Hay múltiples ejemplos de esta ecuación:
La ecuación de calor: en el modelo más sencillo el calor se transporta en dirección opuesta al gradiente \(\bm \flux = -\nabla \rho\) \[ \tag{H} \label{eq-calor} \frac{\partial \rho}{\partial t} - \Delta \rho = 0. \]
Si el espacio no es homogéneo, ha de modelizarse que el coeficiente de difusión depende del punto espacial \(\bm \flux = - A(x) \nabla \rho\) y se llega a la ecuación de difusión \[\begin{equation} \frac{\partial \rho}{\partial t} - \diver(A(x) \nabla \rho) = 0. \end{equation}\]
La ecuación de transporte (en un campo de velocidades): tomando \(\bm \flux = \rho(t,x) \bm v(t,x)\) donde \(\bm v\) es la velocidad del transporte. Tenemos \[\begin{equation} \frac{\partial \rho}{\partial t} + \diver( \rho \bm v) = 0. \end{equation}\]
Observación. En ocasiones la difusión del de forma no lineal de \(\rho\), y \(\bm \flux = - \nabla \varphi (\rho)\). Por ejemplo con \(\varphi(\rho) = \rho^m\), \(m > 0\). Esto da lugar a la famosa Ecuación de Medios Porosos.
Ley de Gauss para el campo eléctrico
Si \(\omega \subset \mathbb R^d\) es un conjunto genérico (regular), \(\bm E\) es el campo eléctrico, \(Q_\omega\) es la carga eléctrica en \(\omega\) y \(\varepsilon_0\) es la permitividad del vacío, entonces la carga eléctrica es proporcional al campo en el borde: \[ \int_{\partial \omega} \bm E(t,x) \cdot \bm n(x) \diff S(x)= \frac{Q_\omega (t)}{\varepsilon_0} \]
Frecuentemente, la cantidad de carga en \(\omega\) se recupera integrando la densidad de carga \(\rho\) (habitualmente conocida), de modo que \[ \int_{\partial \omega} \bm E(t,x) \cdot \bm n(x) \diff S(x) = \frac{1}{\varepsilon_0} \int_\omega \rho(t,x) \diff x. \]
Reproduciendo el argumento anterior \[ \diver\bm E (t,x) = \frac{\rho(t,x)}{\varepsilon_0}. \]
En electroestática, el campo \(\bm E\) no depende del tiempo, y se puede calcular un potencial \(V\) tal que \(\bm E = \nabla V\). Deducimos la ecuación de Laplace \[\begin{equation} -\Delta V (x) = f(x) \end{equation}\] con \(f\) conocida.
Ecuación de Laplace
Volvemos sobre el viejo ejemplo dado por la Ecuación de Laplace \[-u_{xx} (x) = f (x)\]
Método de Diferencias Finitas
Tomando el desarrollo de Taylor \[ \begin{aligned} u(x_{i+1}) &= u(x_i) + u_x(x_i) (x_{i+1}-x_i) + \frac 12 u_{xx}(x_i) (x_{i+1} - x_i)^2 + O(|x_{i+1}-x_i|^3)\\ u(x_i) &= u(x_i) \\ u(x_{i-1}) &= u(x_i) + u_x(x_i) (x_{i-1}-x_i) + \frac 12 u_{xx}(x_i) (x_{i-1} - x_i)^2 + O(|x_{i-1}-x_i|^3) \end{aligned} \]
Malla uniforme. De modo que si \(x_i = a + i\Delta x\) tenemos \[ (\Delta x)^2 u_{xx}(x_i) = \Big( {u(x_{i+1} ) + u(x_{i-1}) - 2u(x_i)} \Big)+ O((\Delta x)^3). \]
Escribimos \[ - \frac{u(x_{i+1} ) + u(x_{i-1}) - 2u(x_i)}{\Delta x^2} = f(x_i) + O(\Delta x). \]
Definimos el problema numérico \[ - \frac{u_{i+1} + u_{i-1} - 2u_i}{(\Delta x)^2} = f(x_i). \] Con condiciones Dirichlet homogéneas \(u_0 = u_{N+1} = 0\) tenemos la matriz \[\begin{equation} \label{eq:Laplaciano diferencias finitas} A = \frac{1}{(\Delta x)^2} \begin{pmatrix} 2 & - 1 \\ -1 & 2 & -1 \\ &&\ddots \\ &&-1 & 2 \end{pmatrix} \end{equation}\]
Extensiones:
- Difusión no homogénea: si el problema es \(-( a(x) u_x ) _x = f\), no es tan fácil escribir la ecuación
Una nueva idea
Si integramos la ecuación en un intervalo \([x_{i-\frac 1 2}, x_{i+\frac 1 2}]\) tenemos \[ u_x (x_{i+\frac 1 2}) - u_x (x_{i-\frac 1 2}) = \int_{x_{i-\frac 1 2}}^{x_{i+\frac 1 2}} u_{xx}(x) \diff x = -\int_{x_{i-\frac 1 2}}^{x_{i+\frac 1 2}} f(x) \diff x \]
Volviendo sobre Taylor, podemos elegir la aproximación \[ u_x (x_{i+\frac 1 2}) = \frac{u(x_{i+1}) - u(x_i)}{x_{i+1}-x_i} + O(|x_{i+1}-x_i|). \]
De modo que \[ \frac{u(x_{i+1}) - u(x_i)}{x_{i+1}-x_i} - \frac{u(x_{i}) - u(x_{i-1})}{x_{i}-x_{i-1}} = -\int_{x_{i-\frac 1 2}}^{x_{i+\frac 1 2}} f(x) \diff x + O(|x_{i+1}-x_i|) + O(|x_{i}-x_{i-1}|) \]
Podemos escribir el método numérico en malla genérica \[ \frac{u_{i+1} - u_i}{x_{i+1}-x_i} - \frac{u_i - u_{i-1}}{x_{i}-x_{i-1}} = -\int_{x_{i-\frac 1 2}}^{x_{i+\frac 1 2}} f(x) \diff x \]
Malla uniforme. Cuando \(x_i = a + ih\) tenemos \[ - \frac{u_{i+1} + u_{i-1} - 2u_i}{h^2} = \int_{x_{i-\frac 1 2}}^{x_{i+\frac 1 2}} f(x) \diff x. \]
El método de Volúmenes Finitos
Volúmenes Finitos
Tomemos la ecuación de calor \(\eqref{eq-calor}\) como ejemplo, y \(d = 1\). Si integramos sobre intervalos de la forma \(K_i = [ x_{i-\frac 1 2} , x_{i+\frac 1 2} ]\) \[ \frac{\diff}{\diff t} \int_{x_{i-\frac 12}}^{x_{i+\frac 1 2}} \rho(t,x) \diff x = - \Big( \flux(t,x_{i+\frac 1 2}) - \flux(t, x_{i-\frac 1 2}) \Big). \]
Esto significa que puedo buscar ecuaciones para \[ \rho_i (t) \approx \frac{1}{x_{i+\frac 1 2}- x_{i-\frac 1 2}} \int_{x_{i-\frac 12}}^{x_{i+\frac 1 2}} \rho(t,x) \diff x. \]
El motivo de esta elección de índices es que \(\rho(t,x_i) \sim \rho_i (t)\). Construimos el método a base de elegir flujos en los bordes del intervalo \[ \frac{\diff}{\diff t} \overline \rho_i (t) = - \frac{ \flux_{i+\frac 1 2} (t) - \flux_{i-\frac 1 2} (t) }{x_{i+\frac 1 2} - x_{i-\frac 1 2}} \]
Cuando depende de \(\rho\) (es decir \(\flux(t,x) = \flux(\rho(t,x))\)), debemos buscar una aproximación en función de las medias a ambos lados, que típicamente será de la forma \[ \flux_{i + \frac 1 2} (t) = G_{i+\frac 1 2}( x, \rho_i (t) , \rho_{i+1} (t) ). \]
Ejemplo. Para la ecuación de calor, \(\bm \flux = -\nabla \rho\).
Usando diferencias centradas, podemos tomar \[ \flux _{i+\frac 1 2} = -\frac{\rho_{i+1} - \rho_i}{x_{i+1} - x_i} \]
Cuando la malla es uniforme recuperamos la semi-discretización espacial usual \[ \frac{\diff}{\diff t} \rho_i = \frac{\rho_{i+1} + \rho_{i-1} - 2 \rho_i}{h^2}. \]
La discretización temporal se puede elegir de forma independiente.
Problemas en dominios acotados
Para poder computar soluciones numéricas, debemos trabajar con una cantidad finita de nodos y buscar \(\rho_1, \cdots, \rho_N\).
Esto es equivalente a que nuestra ecuación esté planteada en un dominio acotado \(\Omega\).
Debemos introducir condiciones de contorno en nuestro problema.
Condiciones de contorno
Condiciones de flujo nulo Dado que \[ \frac {\diff}{\diff t} \int_\Omega \rho(t) \diff x = - \int_{\partial \Omega} \bm F(t) \cdot \bm n \diff x, \] si imponemos que \[ \bm \flux (x) \cdot \bm n (x) = 0 , \qquad \forall x \in \partial \Omega. \] entonces se conservará la masa. Este tipo de condiciones son fáciles de implementar, puesto que simplemente tomaremos \(\flux_{i+\frac 1 2} = 0\) donde nos sea conveniente. Pero hay que tener cuidado con la física del problema.
Condiciones periódicas, donde se supone que el dato inicial se ha repetido periódicamente por \(\mathbb R^d\) (esto debe ser compatible con la geometría de \(\Omega\)), y se espera que la solución mantenga esta periodicidad.
Este tipo de métodos de volúmenes finitos funcionan muy bien para ecuaciones hiperbólicas, entre ellas las que modelizan fluidos. Por eso, son los más utilizados en Ingeniería Aeronáutica y Espacial.
Método de Elementos Finitos
Métodos de tipo Galerkin
Estos métodos se basan en elegir funciones espaciales \(\varphi_i\) y buscar soluciones aproximadas \[ u(t,x) = \sum_{i=1}^N u^h_i (t) \varphi^h_i (x). \]
Son útiles cuando entendemos una llamada “formulación débil” del problema.
Estos método nacen para tratar problemas de tipo calor / elasticidad, y fueron introducidos en Ingeniería de Caminos, Canales, y Puertos.
Para problemas elípticos \(\frac{\partial u}{\partial t} + \mathcal L u = f\) escribimos ecuaciones para \(i = 1, \cdots, N\) \[ \left\langle \frac{\partial u}{\partial t} , \varphi^h_i \right\rangle + \langle \mathcal L u , \varphi^h_i \rangle = \langle f (t) , \varphi^h_i \rangle \]
Si \(\mathcal L\) es lineal, esto se reduce a \[ \sum_{j=1}^N \frac{\diff u_j^h}{\diff t} \left\langle \varphi_j^h , \varphi^h_i \right\rangle + \sum_{j=1}^N u_j^h \langle \mathcal L \varphi_j^h , \varphi^h_i \rangle = \langle f (t) , \varphi^h_i \rangle \]
En forma matricial \[\begin{equation} \label{eq:Galerkin sistema EDOs} M \frac{\diff }{\diff t} \begin{pmatrix} u^h_1 \\ \vdots \\ u^h_N \end{pmatrix} + A \begin{pmatrix} u^h_1 \\ \vdots \\ u^h_N \end{pmatrix} = \begin{pmatrix} \langle f(t) , \varphi^h_1 \rangle \\ \vdots \\ \langle f(t) , \varphi^h_N \rangle \end{pmatrix} \end{equation}\]
donde \(M\) es la llamada matriz de masa y \(A\) la llamada matriz de rigidez \[ \begin{aligned} M &= \begin{pmatrix} \langle \varphi^h_1 , \varphi^h_1 \rangle & \cdots & \langle \varphi^h_N , \varphi^h_1 \rangle \\ \vdots & \ddots & \vdots \\ \langle \varphi^h_1 , \varphi^h_N \rangle & \cdots & \langle \varphi^h_N , \varphi^h_N \rangle \end{pmatrix} \\ A&=\begin{pmatrix} \langle \mathcal L \varphi^h_1 , \varphi^h_1 \rangle & \cdots & \langle \mathcal L \varphi^h_N , \varphi^h_1 \rangle \\ \vdots & \ddots & \vdots \\ \langle \mathcal L \varphi^h_1 , \varphi^h_N \rangle & \cdots & \langle \mathcal L \varphi^h_N , \varphi^h_N \rangle \end{pmatrix}. \end{aligned} \] Como es natural, los problemas estaciones corresponden a eliminar la derivada temporal. Si \(\mathcal L\) no es lineal, la situación es más difícil, como veremos en un ejemplo.
Método de Elementos Finitos
En el método de elementos finitos, esta base de funciones \(\varphi_i\) se eligen a partir de una malla del dominio. A cada uno de los dominios de esta malla se los llama elemento. Un ejemplo clásico es el siguiente.
Ejemplo. Dado el intervalo \([-1,1]\), tomamos puntos una sucesión creciente de puntos \(x_i\). Tomamos \[ \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}, \]
Así \(\varphi_{i} (x_j) = \delta_{ij}\). Por tanto \(u(t,x_i) = a_i(t)\). Si tomamos \(\mathcal L = -\Delta\) y condiciones de contorno Dirichlet, entonces integrando por partes \[ \langle \mathcal L \varphi_i, \varphi_j\rangle = \int_0^1 \frac{\partial \varphi_i} {\partial x}\frac{\partial \varphi_j}{\partial x} \diff x . \]
Como las funciones son constantes a trozos, no es difícil calcular estas integrales analíticamente. Si \(x_i - x_{i-1} = \Delta x\), recuperamos la matriz de diferencias finitas.
Estos métodos tienen algunas ventajas.
Cuando \(\mathcal L\) es lineal, definir \(a(u,v) = \langle \mathcal L u, v \rangle = (Au) \cdot v\) nos da un operador bilineal.
Las propiedades de la matriz de la forma bilineal \((a (\varphi_i, \varphi_j))\), hereda las buenas propiedades del problema continuo: simetría si \(\mathcal L\) es autoadjunto, coercividad si la da \(\mathcal L\),…
Si \(\mathcal L\) no es lineal pero tiene propiedades (caso del \(p\)-Laplaciano, por ejemplo) el problema discreto también las hereda: monotonía, etc…
Julia
El curso puede seguirse en cualquier lenguaje de programación.
Durante el curso, los algoritmos presentados serán mostrados en Julia.