Introducción a los elementos finitos

Hola a todos

Hoy voy a realizar una brevísima introducción sobre uno de los métodos numéricos de mayor uso dentro del mundo de la ciencia e ingeniería: el método de elementos finitos.

Para los que no estais familiarizados, el método de elementos finitos consiste en realizar una aproximación discreta de una función continua que sea solución de una ecuación o sistema de ecuaciones que detallan el comportamiento de un sistema físico o matemático.

Voy a ilustrar un ejemplo unidimensional muy intuitivo para ver cual es la mecánica que se sigue:

Supongamos que se desea resolver una ecuación como esta:

$ \partial_x (p \partial_x y)+ (q+\lambda \sigma) \partial_x y= 0$

Esta ecuación es una expresión general para el problema de Sturm-Liouville, pero en este artículo esto es intrascendental. Se puede ver como una ecuación diferencial de segundo orden con coeficientes variables.

El método de elementos finitos propone la división de la función "y" en n elementos geométricos. Los extremos de cada elemento donde se juntan 2 o más elementos se denominan nodos. Se va a considerar que la función "y" es conocida EXACTAMENTE en todos los nodos. El valor de la misma en los elementos debe ser interpolado. Para realizar esto, se propone para cada elemento una "función nodal" determinada: esta función puede tener diversas formas, pero siempre debe cumplir que sea 1 en el nodo al que corresponde dicha función y 0 en el resto. Puesto que se ha discretizado la función en n+1 nodos donde el valor de la función es conocida, se puede definir la función y como la suma de cada uno de las funciones nodales multiplicadas por el valor de la función en cada una de ellas:

$y=\Sigma_{n=0}^{n+1} \phi_n y_n$

Puesto que esta solución (aunque aproximada) es solución de la ecuación diferencial, podemos sustituirla en la ecuación dada:

$ \partial_x (p \partial_x (\Sigma_{n=0}^{n+1} \phi_n y_n))+ (q+\lambda \sigma) \partial_x (\Sigma_{n=0}^{n+1} \phi_n y_n)= E $

Notad que ahora la ecuación NO esta igualada a 0, lo cual es lógico puesto que hemos "aproximado" la solución real por una solución aproximada y discreta. Este error de aproximación viene dado por la función error E.

Se puede observar que, al hacer esta aproximación se ha podido definir EXACTAMENTE cual es la expresión de la función error de nuestra propuesta numérica. Ahora, sólo hay que variar la función error con un criterio adecuado para hacer tender la función error a 0. Si se consigue, la solución aproximada se acercará a la real. Para hacer esto se recurrirá a lo que se denomina "método de los pesos ponderados".

Este método consiste en hayar una serie de funciones peso tal que la integral de su producto escalar con respecto a la función error E sea cero.

$\int w(x) E dx = 0$

Dicho de otro modo, se quiere encontrar una serie de funciones ortogonales a la función error. La propiedad de ortogonalidad y las funciones peso son importantes: si se elige como función peso las mismas funciones de interpolación nodales que se han introducido en la solución aproximada tal que sean ortogonales a la función peso, entonces es inmediato observar que las funciones de interpolación hayadas harán que la función error E tenga un valor muy próximo a cero y por lo tanto nuestra aproximación sea válida como solución final.

El concepto es muy simple en su fundamento, pero su implementación computacional se muestra mucho más tediosa. Sabemos que la solución aproximada será tanto más exacta cuanto mayor número de nodos se consideren en la solución. Sin embargo, se observa que habrá que resolver un sistema de ecuaciones algebráico con un número de ecuaciones igual al número de nodos que se consideren (¡¡ habrá que cumplir la condición de ortogonalidad para cada función peso !!), aumentando por lo tanto el esfuerzo de computación.

En posts siguientes explicaré más en profundidad algún ejemplo y otras observaciones más detalladas acerca de este método.