Hola a todos
Voy a dedicar este post a repasar brevemente varios conceptos relativos a lo que se denominan "ecuaciones hiperbólicas". Si recordáis, las ecuaciones hiperbólicas son aquellas que poseen una forma general como esta:
U_t + \nabla \cdot F(U) = 0
Para facilitar el ejemplo utilizaré el caso 1D con una función de flujo (a) constante:
U_t + aU_x = 0
con su correspondiente condición inicial:
U(x,0) = U_0
Las ecuaciones hiperbólicas generan como solución un "frente de onda", esto es, una solución que se propaga en el tiempo y en el espacio. De hecho, para este caso en concreto, la solución coincide con un valor desplazado en el tiempo y el espacio:
U(x,t) = U(x-at,0)
Aclararé algunos conceptos que son ligeramente confusos:
- Una ecuación hiperbólica puede generar numerosas familias de soluciones dependiendo de la condición inicial dada. Algunas tendrán significado físico y otras no.
- Existen dos tipos principales de soluciones: frentes de onda y ondas de rarefacción.
- Las familias de soluciones que poseen significado físico son por definición soluciones débiles: estas cumplen la forma de conservación requerida para que una solución sea conservativa:
\int_{x_0}^{x_1} U_t dx = - a \int_{x_0}^{x_1} U_x dx = -a(F(U(x1) - F(U(x0))
- Un frente de onda será solución débil si cumple la condición de Rankine - Hugonniot:
s=\frac{F(U_l) -F(U_r)}{U_l - U_r}
donde U_l y U_r son los valores de U a izquierda y derecha del frente de onda.
- Las ondas de rarefacción son por definición continuas y débiles.
- Cualquier método numérico para resolver un problema hiperbólico debe ser estable, esto es, debe asegurar que dicha solución no va a diverger a infinito (puesto que entonces la solución carece de significado físico).
- Para asegurar la estabilidad de una solución con un método numérico hay que comprobar su condición de entropía.
- La condición de entropía de una solución se puede comprobar utilizando la condición de Lax extendida a posteriori por Oleinik.
\frac{F(U) -F(U_l)}{U- U_l} \geq s \geq \frac{F(U) -F(U_r)}{U - U_r}
- Cualquier método que sea monótono es por definición entrópico. Una solución monótona es aquella que cumple que si dado un par de valores de la función en dos puntos x=x_1 y x=x_2 a tiempo t=0 (u_{1}^0 y u_{2}^0), con u_{1}^0 > u_{2}^0. Entonces se dice que la función es monótona si: u_{1}^t > u_{2}^t para cualesquiera puntos 1 y 2.
- El teorema de Godunov afirma que cualquier método monótono es a lo sumo de orden de exactitud 1. El método de Godunov es de orden 1.
- Existen otros métodos numéricos de mayor exactitud pero que NO son monótonos. En su lugar, preservan la Variación Total. Son los llamados métodos TVD (Total Variation Diminising). La variación total se puede definir como
TV = \int abs(Ux)dx
Como se puede observar, la variación total equivaldría a la cantidad de magnitud que entra en un espacio determinado. Si la TV crece o disminuye quiere decir que la magnitud se acumula o se reduce, lo cual carece de sentido si estamos intentando simular una ecuación que presume ser "conservativa", la cual no debería variar a lo largo de la integración en el tiempo.
- Un método TVD es aquel que cumple que para cualesquiera instantes de tiempo de integración, en cualquier iteración temporal del método se cumple que:
TV(u^{n+1}) \leq TV(u^{n})
Esta condición nos asegura que el método asegura la conservatividad de la solución simulada.
- Existen numerosos métodos de orden superior que utilizan lo que se denomina "limitador de flujo" para conseguir preservar la variación total.
- Cuando se escoge un método hay que elegir que método de reconstrucción de los puntos espaciales se van a utilizar para discretizar los distintos términos espaciales. Esto es, dado un método:
U_t = \frac{1}{\Delta x}\left(F(U_{i+1/2}) - F(U_{i-1/2}) \right)
Según como se discreticen los términos de F, estaremos escogiendo un método u otro. Dicha discretización puede ser constante a trozos (Problema de Riemman), lineal (Método de Godunov) o de orden superior (High resolution squeme).
Hasta aquí el post de hoy. Se profundizará más en este tema en futuros artículos.
jueves, 25 de octubre de 2012
martes, 7 de agosto de 2012
Modelos de autómata celular: Ejemplo
Hola
En este artículo (y aprovechando el periodo vacacional puesto que ahora tengo tiempo), voy a subir al blog una película donde se puede observar cual es el resultado cuando se lleva a cabo una modelización adecuada de un modelo de autómatas celulares.
En este caso se esta simulando el crecimiento bacteriano de un biofilm en presencia de un flujo que avanza de izqda a dcha. Las leyes de evolución implementadas son la reproducción, la erosión del biofilm y el "detachment".
Como punto a remarcar, observad que el segundo 00:07 se observa una forma tipo "champiñón", geometría que este tipo de organismos forma en el mundo real bajo cierto rango de concentración de nutrientes y de condiciones hidrodinámicas suaves.
Espero que sea ilustrativo.
Comentar que este vídeo ha sido elaborado por mi amigo y compañero de trabajo Baldvin Einarsson, excelente programador, que nos a ayudado a llevar adelante nuestras investigaciones sobre los biofilms.
En este artículo (y aprovechando el periodo vacacional puesto que ahora tengo tiempo), voy a subir al blog una película donde se puede observar cual es el resultado cuando se lleva a cabo una modelización adecuada de un modelo de autómatas celulares.
En este caso se esta simulando el crecimiento bacteriano de un biofilm en presencia de un flujo que avanza de izqda a dcha. Las leyes de evolución implementadas son la reproducción, la erosión del biofilm y el "detachment".
Como punto a remarcar, observad que el segundo 00:07 se observa una forma tipo "champiñón", geometría que este tipo de organismos forma en el mundo real bajo cierto rango de concentración de nutrientes y de condiciones hidrodinámicas suaves.
Espero que sea ilustrativo.
Comentar que este vídeo ha sido elaborado por mi amigo y compañero de trabajo Baldvin Einarsson, excelente programador, que nos a ayudado a llevar adelante nuestras investigaciones sobre los biofilms.
Modelos de autómatas celulares
Saludos a todos
Con este artículo no quiero decir que este tipo de modelos es mucho mejor que el resto. Sólo quiero darlos a conocer a todos aquellos que no conocen el campo y mostrar tanto sus puntos fuertes como sus debilidades.
Hoy quería introducir brevemente unas pequeñas nociones sobre una parte concreta del modelado matemático en el que he trabajado durante cierto tiempo y que considero una herramienta muy potente si se sabe utilizar con criterio. Estoy refiriéndome a los modelos de autómata celular.
Este tipo de modelos aplican un enfoque para abordar problemas bastante interesante. Si bien prácticamente todos los modelos físicos abordan los problemas utilizando las ecuaciones de conservación físicas ya sea en su forma diferencial o en su forma integrada (las cuales son deterministas), los modelos de autómata celular avanzan en otra dirección: estos modelos utilizan un enfoque plenamente estocástico (utilizando variables aleatorias). Se podría describir el modelado de este "approach" de la siguiente manera:
- Se discretiza en un mallado formado por celdas cuadradas (se podrían realizar otro tipo de celdas) el dominio que se quiere modelizar. Este dominio va a ser modificado por la fenomenología física a describir.
Ejemplos podrían ser la descripción de un fluido o un sólido sometido a fuerzas de diferente tipo que van a modificar la geometría del propio dominio.
- Cada celda puede adoptar distintos estados, normalmente asociados a os distintos subdominios de nuestro sistema. Por ejemplo, si nuestro sistema modeliza un fluido que se mueve por efecto del viento sobre un lecho de piedras, los estados que las celdas podrían adoptar podrían sera agua, aire y rocas.
- Dos dimensiones claves para modelizar van a ser el tamaño de la celda de la malla, la cual va a dar un significado físico a modelo que se va a llevar a cabo (y por ende determinará que procesos actúan a esa escala), y el paso temporal, el cual va a ser implícito e irá introducido a través de los procesos de iteración del modelo.
-Se observará el fenómeno físico determinado que se quiere modelar y se descompondrá la fenomenología en distintos "procesos" o "efectos" que van a actuar sobre el dominio. Al ser el modelo totalmente estocástico, cada proceso puede actuar o no sobre cada celda del dominio.
- Para determinar si un proceso ocurre o no, se propondrá una ley que permita calcular la probabilidad (valor entre 0 y 1) de ocurrencia de un fenómeno dado a partir de ciertos datos que se pueden calcular mediante otros procesos de cálculo ( campos de velocidad, de concentración , etc. determinados ya sea de forma convencional resolviendo PDEs que aplican en el dominio ). En este punto es conveniente recomendar que el cálculo de estos datos adicionales sea realizado utilizando hipótesis simplificatorias, ya que si no se corre el riesgo de generar un coste computacional demasiado elevado (lo cual es un inconveniente que se intenta evitar cuando se recurre a este tipo de modelos).
- Una vez obtenida la probabilidad, se realizarán sorteos aleatorios para cada celda para cada proceso. Esto provocará un efecto dentro del sistema, haciendo que cada celda cambie el valor de su estado en función del efecto producido por los fenómenos simulados. Por ejemplo, si el aire desplaza el agua, una casilla concreta que pudiera adoptar el estado de agua podría cambiar su estado a aire si la fuerza simulada ha desplazado el fluido de esa celda del dominio.
- En cuanto se ha actualizado todas las casillas, se incrementa el contador temporal del modelo y se vuelve a reiniciar el proceso. Este bucle simulará el paso del tiempo en el sistema.
Hay que destacar varios aspectos:
1 - El modelo es totalmente aleatorio: No hay 2 simulaciones idénticas, puesto que todos los eventos son estocásticos
2- La dinámica del sistema surge de forma espontánea a través de los distintos resultados que van modelando el dominio de celdas: "Nadie dice como tiene que evolucionar el sistema, sino que sólo se restringen los grados de libertad que este puede adoptar"
3- Este último punto nos hace deducir que sin una buena modelización de la física que se quiere simular, los resultados obtenidos van a resultar ilógicos y no realistas. Es imprescindible modelar de forma correcta las leyes de probabilidad que rigen los distintos eventos implementados en el modelo.
4- Muchos parámetros que se introduzcan en dichas leyes pueden ser inventados y sin embargo el modelo puede describir de forma realista un determinado fenómeno. El paso lógico será realizar un proceso de CALIBRACIÓN que permita obtener el rango de parámetros adecuado en el cual el modelo da resultados acordes con la realidad experimental.
5- Este tipo de modelos permite obtener resultados asombrosamente similares a los obtenidos en experimentación con un coste computacional muy por debajo de los métodos tradicionales de resolución de PDEs.
Con este artículo no quiero decir que este tipo de modelos es mucho mejor que el resto. Sólo quiero darlos a conocer a todos aquellos que no conocen el campo y mostrar tanto sus puntos fuertes como sus debilidades.
Conservatividad en las ecuaciones.
Un concepto que a mi juicio es importante en física y matemáticas que no se estudia con suficiente detalle en la universidad y que es importante para la comprensión de las ecuaciones matemáticas que rigen el comportamiento del mundo que nos rodean es el concepto de conservación.
u_t +F(u)_x = 0
Hablando toscamente, se puede afirmar que una ecuación se denomina conservativa cuando esta cumple un principio universal que, si bien a priori todo el mundo piensa que es "de cajón", tarda cierto tiempo en ser asimilado. El principio es el siguiente:
Variación en un dominio de una magnitud = entra - sale
Esto resulta trivial e incluso insultante. Pero detallemos en mayor profundidad las consecuencias de esto cuando trabajamos con ecuaciones diferenciales. Para este caso, y suponiendo que nuestro problema consta de dimensión temporal y espacial, podríamos decir que una ecuación es de tipo conservativo si se puede escribir de esta forma:
u_t + \nabla \cdot F(u) = 0
donde F(u) puede ser una función no lineal de u (no incluiría funciones donde intervengan términos en derivadas!)
A priori esta ecuación no esta muy clara que signifique lo anterior. Si trabajamos con 1 dimensión espacial la expresión se puede escribir como:
u_t +F(u)_x = 0
Pero si integramos la ecuación:
\int_{x1}^{x2}u_t dx = -[F(u_2) - F(u_1)]
Aquí se puede ver claramente lo que se ha comentado anteriormente: la ecuación realmente nos dice que la variación de la magnitud u con el tiempo en todo el dominio es igual a la cantidad de u que entra menos la que sale.
Este concepto, que quizá no parezca muy trascendental, es importante en el sentido de que las ecuaciones que miden magnitudes conservativas (energía, momento lineal, masa, etc.), todas obedecen esta estructura básica. Pueden llevar más términos o menos, pero esta estructura siempre esta incluida en las ecuaciones. Esto es coherente puesto que de otro modo se violaría la conservación de dichas magnitudes, caso que hasta el momento no se ha dado en ningún experimento.
La consevación tiene consecuencias más profundas. En concreto esta intimamente ligada a la resolución numérica de ecuaciones diferenciales. Los que os habéis informado un poco más en profundidad sobre métodos numéricos para resolver PDEs seguramente os habréis topado con el concepto de "formulación fuerte" y "formulación débil" de una ecuación diferencial. Estos conceptos están relacionados la conservatividad. Se dice que una ecuación satisface su formulación fuerte cuando la solución de la ecuación cumple estrictamente tanto con las condiciones de contorno como con la formulación diferencial en todo punto del dominio. Sin embargo, en muchos campos resolver la formulación fuerte de dichas ecuaciones resulta muy complicado (debido a que pueden aparecer discontinuidades o singularidades en la solución que hacen la ecuación no diferenciable aunque si sea continua), y en bastantes casos es hasta el momento imposible (ejemplo clásico son las ecuaciones de Navier - Stokes). Sin embargo, lo que se plantea es utilizar una formulación débil, la cual cumple con las condiciones de contorno pero no cumple la ecuación diferencial en si, sino su FORMULACIÓN INTEGRAL. La formulación integral garantiza el poder obtener soluciones incluso aunque estas no sean diferenciables, puesto que no hay necesidad de derivar. Puesto que al integrar en un dominio acotado, el resultado de la integral se evalúa entre 2 límites de integración y se restan, se puede derivar que cuando se aplica una formulación débil de una PDE donde tiempo y espacio son las dimensiones relevantes que describa un proceso físico, realmente estamos reafirmando que la magnitud SE CONSERVA.
viernes, 20 de julio de 2012
Tipos de PDEs
Hola a todos
Seguramente a muchos se nos viene a la cabeza preguntarnos que es una ecuación hiperbólica o parabólica. Voy a intentar detallar brevemente lo que son estos conceptos como los entiendo desde mi punto de vista.
Las ecuaciones en derivadas parciales (Partial differential equations en ingles, o PDE) se pueden clasificar atendiendo a varios criterios de clasificación, pero uno de los más utilizados divide las PDEs en 3 grandes grupos: parabólicas, elípticas e hiperbólicas. El criterio concreto para dividir los 3 grupos es la presencia o no en la ecuación de ciertos términos. Considerando que la ecuación se puede escribir de esta forma genérica:
A·u_{xx} + 2B·u_{xy} + C·u_{yy}+D·u_x+E·u_y+F = 0
Las ecuaciones se englobarán en una categoría u otra según los distintos coeficientes en letras mayúsculas en la ecuación sean cero o no. No obstante, mi criterio a la hora de recordar que tipo de ecuación tengo delante cuando la veo es guiarme por 3 criterios básicos y luego asociarlos.
Un método más intuitivo de aprender a indentificar cada tipo de ecuación a su forma matemática es tomando ejemplos clásicos de física, entendiendo lo que ocurre en la realidad y asociarle una forma matemática estándar que producirá dicho efecto numéricamente. Voy a poner un ejemplo para cada tipo:
· Ecuación parabólica: Un ejemplo clásico asociado a este tipo de ecuaciones es la ecuación del calor. Esta ecuación describe como se propaga el calor en un medio: si pensamos en un frigorífico cuyas paredes interiores están a una temperatura T_{cold} e introducimos un vaso de agua caliente a temperatura T_{hot}, la intuición nos dice que el calor fluirá del vaso a través del interior del frigorífico hacia las paredes (porque están más frías) de forma gradual y sostenida en el tiempo, de tal modo que al principio el agua perderá calor rápidamente pero conforme esté más templada tardará más tiempo en seguir enfriándose hasta que finalmente el agua alcanza la misma temperatura que la nevera, y el proceso se detiene. Vemos que en este proceso intervienen 2 variables dimensionales importantes: el tiempo y la distancia entre el vaso y las paredes (no es lo mismo que estén muy cerca que estén muy separadas del vaso). Juntando estas, la ecuación se puede escribir como:
\partial_t T = k \partial_{xx} T
con sus condiciones de contorno escritas de forma no ortodoxa (para que sea entendible):
T(x = pared )= T_{cold}
T(x = vaso) = T_{hot}
Lo que describe la ecuación es sencillamente como varía la temperatura en el interior de la nevera al pasar el tiempo en cada punto del interior de la nevera. Por lo tanto, si hablamos de ecuaciones parabólicas en procesos físicos donde las dimensiones espaciales y temporales están involucradas, entonces podemos afirmar que describen la evolución no estacionaria (hay un término con derivada temporal) de una magnitud (en este caso la temperatura) en un dominio dado (sería la zona del espacio donde quieres conocer la distribución de temperaturas). Esta ecuación da la temperatura para cada punto del espacio y para cada instante de tiempo.
· Ecuación elíptica: Un ejemplo para este tipo de ecuaciones es la ecuación de difusión. Consideremos una tubería rectangular llena de agua donde supongamos introducimos un tinte con una concentración C_{in} a través de uno de sus extremos. Supongamos también que por el otro lado de la tubería sale una concentración de tinte C_{out}. Si no hay corriente de agua alguna de líquido (el líquido esta en reposo), la pregunta que se nos plantea es como varíará la concentración de tinte a través de la tubería. La intuición nos vuelve a decir que el tinte "difundirá" progresivamente desde un lado de la tubería hasta el otro de tal manera que al observar la variación de la concentración a lo largo de la tubería, esta variará de forma progresiva.
La ecuación que describiría esta variación será:
D \partial_{xx} C= 0 T
con sus condiciones de contorno:
C(x = entrada )= C_{in}
C(x = salida) = C_{out}
Observad que en este caso, la ecuación NO DICE NADA SOBRE LA EVOLUCIÓN TEMPORAL DE LA CONCENTRACIÓN. Sólo da la DISTRIBUCIÓN ESPACIAL DE CONCENTRACIONES. Por lo tanto, las ecuaciones elípticas en problemas de física con dimensiones espaciales están asociadas a distribuciones estacionarias de magnitudes. De hecho, si tomamos la ecuación parabólica del calor (que es no estacionaria), y hacemos que el tiempo tienda a infinito, se podrá observar que el término temporal tiende a anularse, y la solución tiende a ser la misma que la de una ecuación elíptica.
· Ecuación hiperbólica: El ejemplo clásico que se suele poner en los libros de texto es el de la ecuación de ondas. Sin embargo, esto a priori no debe decir mucho a alguien no familiarizado con el mundo de las ondas.
Veamos un ejemplo inverso para entenderlo. Primero asumiremos que un ejemplo de una ecuación de ondas puede ser el siguiente:
H_t + c· H_x = 0
con condición inicial
H(x,0) = H_0(x)
Esta ecuación en el ejemplo que tratamos referirá a como varía en el tiempo la altura de una cuerda. Si nos fijamos en la ecuación, hay 2 magnitudes independientes: tiempo y espacio. Luego la solución debe ser función de estas. Esta ecuación expresa que la variación en el tiempo la altura de la cuerda es igual a como varía en dirección longitudinal.
Sabemos que H(x,0) es la condición inicial del sistema, por lo que satisface la ecuación. Es fácil comprobar que una solución idéntica espacialmente a esta pero trasladada en el tiempo exactamente c veces un intervalo de tiempo t satisface también la ecuación, esto es: H(x-ct, t) también satisfará la ecuación. Esto se podría entender mejor si efectivamente asumimos que la solución posee una forma de onda: si a tiempo 0 la forma de la onda es mi solución, si yo me desplazo a una velocidad EXACTAMENTE IGUAL a la que viaja la onda, esta no ha variado su forma desde mi punto de vista, por lo que en cualquier instante de tiempo también será solución. Efectivamente si me quedo estático podré observar como la cuerda se elevará y ddescenderá físicamente formando una geometría igual a la de la onda original. La solución es una onda que se propaga físicamente por la cuerda, pero también en el tiempo, ya que a tiempo suficiente da igual en que punto de la cuerda este, veré llegar la onda y la veré marcharse a una velocidad c con exactamente la misma geometría que a instante t = 0.
Espero que os sirva de ayuda.
Seguramente a muchos se nos viene a la cabeza preguntarnos que es una ecuación hiperbólica o parabólica. Voy a intentar detallar brevemente lo que son estos conceptos como los entiendo desde mi punto de vista.
Las ecuaciones en derivadas parciales (Partial differential equations en ingles, o PDE) se pueden clasificar atendiendo a varios criterios de clasificación, pero uno de los más utilizados divide las PDEs en 3 grandes grupos: parabólicas, elípticas e hiperbólicas. El criterio concreto para dividir los 3 grupos es la presencia o no en la ecuación de ciertos términos. Considerando que la ecuación se puede escribir de esta forma genérica:
A·u_{xx} + 2B·u_{xy} + C·u_{yy}+D·u_x+E·u_y+F = 0
Las ecuaciones se englobarán en una categoría u otra según los distintos coeficientes en letras mayúsculas en la ecuación sean cero o no. No obstante, mi criterio a la hora de recordar que tipo de ecuación tengo delante cuando la veo es guiarme por 3 criterios básicos y luego asociarlos.
Un método más intuitivo de aprender a indentificar cada tipo de ecuación a su forma matemática es tomando ejemplos clásicos de física, entendiendo lo que ocurre en la realidad y asociarle una forma matemática estándar que producirá dicho efecto numéricamente. Voy a poner un ejemplo para cada tipo:
· Ecuación parabólica: Un ejemplo clásico asociado a este tipo de ecuaciones es la ecuación del calor. Esta ecuación describe como se propaga el calor en un medio: si pensamos en un frigorífico cuyas paredes interiores están a una temperatura T_{cold} e introducimos un vaso de agua caliente a temperatura T_{hot}, la intuición nos dice que el calor fluirá del vaso a través del interior del frigorífico hacia las paredes (porque están más frías) de forma gradual y sostenida en el tiempo, de tal modo que al principio el agua perderá calor rápidamente pero conforme esté más templada tardará más tiempo en seguir enfriándose hasta que finalmente el agua alcanza la misma temperatura que la nevera, y el proceso se detiene. Vemos que en este proceso intervienen 2 variables dimensionales importantes: el tiempo y la distancia entre el vaso y las paredes (no es lo mismo que estén muy cerca que estén muy separadas del vaso). Juntando estas, la ecuación se puede escribir como:
\partial_t T = k \partial_{xx} T
con sus condiciones de contorno escritas de forma no ortodoxa (para que sea entendible):
T(x = pared )= T_{cold}
T(x = vaso) = T_{hot}
Lo que describe la ecuación es sencillamente como varía la temperatura en el interior de la nevera al pasar el tiempo en cada punto del interior de la nevera. Por lo tanto, si hablamos de ecuaciones parabólicas en procesos físicos donde las dimensiones espaciales y temporales están involucradas, entonces podemos afirmar que describen la evolución no estacionaria (hay un término con derivada temporal) de una magnitud (en este caso la temperatura) en un dominio dado (sería la zona del espacio donde quieres conocer la distribución de temperaturas). Esta ecuación da la temperatura para cada punto del espacio y para cada instante de tiempo.
· Ecuación elíptica: Un ejemplo para este tipo de ecuaciones es la ecuación de difusión. Consideremos una tubería rectangular llena de agua donde supongamos introducimos un tinte con una concentración C_{in} a través de uno de sus extremos. Supongamos también que por el otro lado de la tubería sale una concentración de tinte C_{out}. Si no hay corriente de agua alguna de líquido (el líquido esta en reposo), la pregunta que se nos plantea es como varíará la concentración de tinte a través de la tubería. La intuición nos vuelve a decir que el tinte "difundirá" progresivamente desde un lado de la tubería hasta el otro de tal manera que al observar la variación de la concentración a lo largo de la tubería, esta variará de forma progresiva.
La ecuación que describiría esta variación será:
D \partial_{xx} C= 0 T
con sus condiciones de contorno:
C(x = entrada )= C_{in}
C(x = salida) = C_{out}
Observad que en este caso, la ecuación NO DICE NADA SOBRE LA EVOLUCIÓN TEMPORAL DE LA CONCENTRACIÓN. Sólo da la DISTRIBUCIÓN ESPACIAL DE CONCENTRACIONES. Por lo tanto, las ecuaciones elípticas en problemas de física con dimensiones espaciales están asociadas a distribuciones estacionarias de magnitudes. De hecho, si tomamos la ecuación parabólica del calor (que es no estacionaria), y hacemos que el tiempo tienda a infinito, se podrá observar que el término temporal tiende a anularse, y la solución tiende a ser la misma que la de una ecuación elíptica.
· Ecuación hiperbólica: El ejemplo clásico que se suele poner en los libros de texto es el de la ecuación de ondas. Sin embargo, esto a priori no debe decir mucho a alguien no familiarizado con el mundo de las ondas.
Veamos un ejemplo inverso para entenderlo. Primero asumiremos que un ejemplo de una ecuación de ondas puede ser el siguiente:
H_t + c· H_x = 0
con condición inicial
H(x,0) = H_0(x)
Esta ecuación en el ejemplo que tratamos referirá a como varía en el tiempo la altura de una cuerda. Si nos fijamos en la ecuación, hay 2 magnitudes independientes: tiempo y espacio. Luego la solución debe ser función de estas. Esta ecuación expresa que la variación en el tiempo la altura de la cuerda es igual a como varía en dirección longitudinal.
Sabemos que H(x,0) es la condición inicial del sistema, por lo que satisface la ecuación. Es fácil comprobar que una solución idéntica espacialmente a esta pero trasladada en el tiempo exactamente c veces un intervalo de tiempo t satisface también la ecuación, esto es: H(x-ct, t) también satisfará la ecuación. Esto se podría entender mejor si efectivamente asumimos que la solución posee una forma de onda: si a tiempo 0 la forma de la onda es mi solución, si yo me desplazo a una velocidad EXACTAMENTE IGUAL a la que viaja la onda, esta no ha variado su forma desde mi punto de vista, por lo que en cualquier instante de tiempo también será solución. Efectivamente si me quedo estático podré observar como la cuerda se elevará y ddescenderá físicamente formando una geometría igual a la de la onda original. La solución es una onda que se propaga físicamente por la cuerda, pero también en el tiempo, ya que a tiempo suficiente da igual en que punto de la cuerda este, veré llegar la onda y la veré marcharse a una velocidad c con exactamente la misma geometría que a instante t = 0.
Espero que os sirva de ayuda.
lunes, 7 de mayo de 2012
Level - Set y método adjunto: comparativa
Estimados lectores
En este post voy a mostrar los resultados que se obtienen al utilizar dos métodos distintos para dibujar los resultados de un problema de optimización: en este caso se trata de la resolución de la "ecuación de radiación reducida", o también llamada ecuación de difusión (es el caso en el que el medio es ópticamente "grueso", esto es, que la absorción es despreciable frente al "scattering":
\nabla \cdot (D(r) \nabla I(r)) + \mu_a (r)I(r) = 0 \text{ in } \Omega
\mu_a (r)I(r) - \frac{2}{\pi} D(r) \nu(r_b) \nabla I(r)= f(r) \text{ in } \partial \Omega
He implementado dos soluciones: por un lado el método adjunto y por otro un método level-set. Las gráficas se anexan a continuación. Para el método level set se han utilizado dos pares de valores para el campo de valores de coeficiente de absorción del medio en la malla: por un lado un valor de 0.1 para el background y 0.03 para el objeto y para otro caso distinto un valor de 0.1 para el background y de 0.17 para el objeto. Esto se ha planteado así para analizar los resultados y ver las conclusiones.
En este post voy a mostrar los resultados que se obtienen al utilizar dos métodos distintos para dibujar los resultados de un problema de optimización: en este caso se trata de la resolución de la "ecuación de radiación reducida", o también llamada ecuación de difusión (es el caso en el que el medio es ópticamente "grueso", esto es, que la absorción es despreciable frente al "scattering":
\nabla \cdot (D(r) \nabla I(r)) + \mu_a (r)I(r) = 0 \text{ in } \Omega
\mu_a (r)I(r) - \frac{2}{\pi} D(r) \nu(r_b) \nabla I(r)= f(r) \text{ in } \partial \Omega
He implementado dos soluciones: por un lado el método adjunto y por otro un método level-set. Las gráficas se anexan a continuación. Para el método level set se han utilizado dos pares de valores para el campo de valores de coeficiente de absorción del medio en la malla: por un lado un valor de 0.1 para el background y 0.03 para el objeto y para otro caso distinto un valor de 0.1 para el background y de 0.17 para el objeto. Esto se ha planteado así para analizar los resultados y ver las conclusiones.
El resultado del método adjunto (resultado arriba ilustrado) muestra que en el dominio estudiado hay dos objetos
diagonalmente opuestos con unos coeficientes de difusión mayores y
menores que el valor de fondo del dominio, cuyo valor es de 0.1 . La
velocidad de convergencia fue muy lenta, del orden de horas.
Este resultado puede
confirmarse al desarrollar los métodos Level-Set:
- Al utilizar como
valores de referencia 0.1 para el fondo y 0.03 para el objeto, se
puede observar que el método level-set da como resultado un contorno
cercano al objeto con valor de coeficiente de absorción bajo en el
método adjunto, confirmando la ubicación predicha por el
primer método. El resultado se obtuvo en un tiempo de convergencia
mucho más bajo que en el caso del método adjunto. El resultado se ilustra a continuación:
- Al utilizar como
valores de referencia 0.1 para el fondo y 0.17 para el objeto se
visualiza que el método level-set da como resultado final el objeto
correspondiente a un coeficiente de absorción elevado obtenido en el
método adjunto. La velocidad de convergencia fue similar al
primer caso level-set.
Estos resultados permiten
concluir que:
- El método adjunto permite obtener en tu dominio todos los objetos que estén presentes
en el mismo, sea cual sea su coeficiente de absorción.
- El resultado de los
métodos level-set mostrarán un resultado u otro dependiendo del
valor escogido para el background y para el objeto. Si estos valores
no son cogidos de forma adecuada, puede que la solución obtenida sea
incompleta. Haría falta resolver el problema varias veces con
distintos valores para cerciorarse que el resultado obtenido es
correcto.
- El método adjunto es claramente más lento en términos de convergencia que los métodos
level-set, necesitando un número de iteraciones mucho mayor. Esto
concuerda con el hecho de que este método busca soluciones de coef.
de absorción en un campo de valores a priori infinito. Sin embargo,
el método level-set busca soluciones en un campo de coeficientes de
absorción limitado a dos valores: fondo u objeto. Por lo tanto se
espera que level- set sea mucho más rápido.
- El método adjunto da una solución difuminada en el espacio, esto es, la solución no
se sitúa en un punto sino en un área más o menos extensa.
- El método level-set
permite ubicar la posición del objeto de forma mucho más precisa,
ya que el campo de soluciones sólo permite dos valores para cada
punto de la malla.
- Como resumen de estas
conclusiones se puede comentar que el método level-set es un método
mucho más eficaz para si se conoce a priori el rango de valores que
puede adoptar el campo buscado. Si no se conoce información alguna,
el método adjunto puede resultar una opción acertada debido a
que permite averiguar todos los objetos con cualquier valor dentro
del dominio.
Si se deseara una buena
precisión, la opción preferida sería utilizar el método adjunto para acotar el valor que puede adoptar el campo estudiado y
posteriormente definir la posición de los objetos de forma mucho más
precisa con el método level-set.
Los resultados, como se pueden observar, no son exactos: como todo método de optimización, tu resultado final dependerá de lo preciso que sea el algoritmo de convergencia y de la estimación inicial que se realice. No obstante, este ejemplo permite visualizar la potencia de este tipo de técnicas.
Los resultados, como se pueden observar, no son exactos: como todo método de optimización, tu resultado final dependerá de lo preciso que sea el algoritmo de convergencia y de la estimación inicial que se realice. No obstante, este ejemplo permite visualizar la potencia de este tipo de técnicas.
Un saludo
lunes, 23 de abril de 2012
Métodos Level-Set
Estimados lectores
En esta ocasión voy a escribir sobre un tema sobre el cual estoy actualmente desarollando un trabajo y que me ha parecido ciertamente util: los métodos level-set.
Este tipo de métodos se engloban dentro de los métodos matemáticos de optimización, especialmente en aquellos destinados a la resolución de problemas inversos: ejemplos clásicos de problemas inversos pueden ser la reconstrucción de imágenes o la determinación de variables de campo implícitas mediante el uso de un número de datos obtenidos mediante la resolución del problema explícito (por ejemplo conocer el coeficiente de absorción de un dominio a partir de varios estímulos en la frontera con ciertas fuentes de luz de valor conocido).
Los métodos tipo level - set fueron planteados hace tiempo por Sethian y Onsager. Plantearon que resolviendo la siguiente ecuación, se podía conocer la evolución de la frontera de una función dada \phi conocido el campo de velocidades \vec{v} de cada uno de los puntos de la frontera:
\frac{\partial \phi}{\partial t} + \vec{v} \vec{\nabla} \phi = 0
con
\vec{\nabla} \phi = \vec{n} \mid \nabla \phi \mid
El método en si es más sencillo de lo que parece. El punto clave consiste en que dada una función F en una dimensión arbitraria, por ejemplo dimensión 2, el problema se puede generalizar utilizando una función de una dimensión adicional \phi (en el ejemplo dimensión 3) y buscar el valor de mi función cuando esta dimensión adicional adquiere un valor 0 "busco el nivel cero de la función \phi.
F(x,y,t) \longrightarrow \phi(x,y,z,t) \longrightarrow z = \phi(x,y,t)
\Gamma (x,y,t) \equiv \phi(x,y,0,t)
En este caso, resulta que los valores de la función que hagan cero mi función son los valores de la frontera que estoy buscando.
Una vez conocida mi función \phi de dimensión adicional, y conociendo el campo de velocidades \vec{v}, mediante la integración de la ecuación hiperbólica anteriormente descrita puedo hallar la evolución de mi frontera.
¿Como se traduce esto a efectos prácticos de resolver un problema concreto?
Es relativamente sencillo. Dada la función F de tu problema dentro del dominio estudiado, se propone una función \phi tal que \Gamma = \phi = 0
CUALQUIER FUNCIÓN PUEDE VALER. A priori vale cualquier forma siempre y cuando cumpla que en el instante inicial \Gamma = \phi = 0. Un ejemplo para el caso 3D puede ser un paraboloide de revolución clásico.
z=(\frac{x}{a})^2 + (\frac{y}{b})^2
Una vez se propone esta función, se comienza a integrar la PDE hiperbólica propuesta arriba utilizando el campo de velocidades del dominio. En cada paso "temporal", la función \phi se irá deformando, pero sólo su proyección en el plano dimensional del problema adicional es nuestra solución (\Gamma). Así que representando al final de cada paso iterativo \phi = 0, obtendremos la deformación de la frontera.
Este tipo de métodos pueden ser combinado conjuntamente con otros métodos de optimización, como el método del gradiente. Utilizando la función de optimización respectiva y ajustando el tamaño del incremento mediante un parámetro escalar variable en magnitud en cada paso, se pueden resolver problemas inversos como se ha comentado anteriormente (como una reconstrucción de imágenes).
Un saludo
En esta ocasión voy a escribir sobre un tema sobre el cual estoy actualmente desarollando un trabajo y que me ha parecido ciertamente util: los métodos level-set.
Este tipo de métodos se engloban dentro de los métodos matemáticos de optimización, especialmente en aquellos destinados a la resolución de problemas inversos: ejemplos clásicos de problemas inversos pueden ser la reconstrucción de imágenes o la determinación de variables de campo implícitas mediante el uso de un número de datos obtenidos mediante la resolución del problema explícito (por ejemplo conocer el coeficiente de absorción de un dominio a partir de varios estímulos en la frontera con ciertas fuentes de luz de valor conocido).
Los métodos tipo level - set fueron planteados hace tiempo por Sethian y Onsager. Plantearon que resolviendo la siguiente ecuación, se podía conocer la evolución de la frontera de una función dada \phi conocido el campo de velocidades \vec{v} de cada uno de los puntos de la frontera:
\frac{\partial \phi}{\partial t} + \vec{v} \vec{\nabla} \phi = 0
con
\vec{\nabla} \phi = \vec{n} \mid \nabla \phi \mid
El método en si es más sencillo de lo que parece. El punto clave consiste en que dada una función F en una dimensión arbitraria, por ejemplo dimensión 2, el problema se puede generalizar utilizando una función de una dimensión adicional \phi (en el ejemplo dimensión 3) y buscar el valor de mi función cuando esta dimensión adicional adquiere un valor 0 "busco el nivel cero de la función \phi.
F(x,y,t) \longrightarrow \phi(x,y,z,t) \longrightarrow z = \phi(x,y,t)
\Gamma (x,y,t) \equiv \phi(x,y,0,t)
En este caso, resulta que los valores de la función que hagan cero mi función son los valores de la frontera que estoy buscando.
Una vez conocida mi función \phi de dimensión adicional, y conociendo el campo de velocidades \vec{v}, mediante la integración de la ecuación hiperbólica anteriormente descrita puedo hallar la evolución de mi frontera.
¿Como se traduce esto a efectos prácticos de resolver un problema concreto?
Es relativamente sencillo. Dada la función F de tu problema dentro del dominio estudiado, se propone una función \phi tal que \Gamma = \phi = 0
CUALQUIER FUNCIÓN PUEDE VALER. A priori vale cualquier forma siempre y cuando cumpla que en el instante inicial \Gamma = \phi = 0. Un ejemplo para el caso 3D puede ser un paraboloide de revolución clásico.
z=(\frac{x}{a})^2 + (\frac{y}{b})^2
Una vez se propone esta función, se comienza a integrar la PDE hiperbólica propuesta arriba utilizando el campo de velocidades del dominio. En cada paso "temporal", la función \phi se irá deformando, pero sólo su proyección en el plano dimensional del problema adicional es nuestra solución (\Gamma). Así que representando al final de cada paso iterativo \phi = 0, obtendremos la deformación de la frontera.
Este tipo de métodos pueden ser combinado conjuntamente con otros métodos de optimización, como el método del gradiente. Utilizando la función de optimización respectiva y ajustando el tamaño del incremento mediante un parámetro escalar variable en magnitud en cada paso, se pueden resolver problemas inversos como se ha comentado anteriormente (como una reconstrucción de imágenes).
Un saludo
Disculpas a todos los lectores por la tardanza
Estimados lectores
Lamento profundamente no haber podido atender mi blog como es debido. Las obligaciones personales y con mi tesis doctoral no me han permitido actualizar los contenidos ni atender las preguntas a tiempo.
Espero en los próximos días ir actualizando con nuevo contenido mi blog.
Sin más, un saludo a todos
DAVID RODRIGUEZ
Lamento profundamente no haber podido atender mi blog como es debido. Las obligaciones personales y con mi tesis doctoral no me han permitido actualizar los contenidos ni atender las preguntas a tiempo.
Espero en los próximos días ir actualizando con nuevo contenido mi blog.
Sin más, un saludo a todos
DAVID RODRIGUEZ
lunes, 6 de febrero de 2012
Computación en paralelo I
Hola a todos
En esta ocasión voy a hablar muy escuetamente de lo que va a ser el futuro de la programación científica: la computación en paralelo.
Actualmente los procesadores tienden a tener un número de núcleos mayor de uno, tal y como era hace unos años: todos estamos familiarizados con los "Pentium core 2 duo" o los modernos "Pentium I3, I5, etc". Sin embargo, por lo que a mi respecta, la formación universitaria en ingeniería sólo cubre nociones básicas de programación en serie. Es obvio que, desde un punto de vista computacional, si se poseen varios núcleos de procesamiento, el realizar una programación capaz de utilizar simultáneamente la potencia de cálculo de esos cores va a ser mucho más ventajoso.
Sin embargo, hasta hace bien poco desconocía por completo un punto de partida que me pudiera orientar para empezar a caminar en este vasto campo. Pues bien, si queréis una piedra angular de referencia, vuestra palabra mágica es MPI.
El MPI se puede considerar como un anexo a los lenguajes de programación típicos en ciencias (Fortran, C) que permite procesar un código arrancando varios "threads", o hilos de computación. Realmente lo ideal es ejecutar tantos threads como núcleos tengas en tu cluster, pero no tiene por que ser así.
Para realizar un código en paralelo, hay que pasar de una mentalidad de cálculo en serie, y comenzar a pensar como se puede distribuir todas las tareas que tiene que hacer el código: habrá tareas, como una búsqueda de elementos en un fichero o matriz, que se puede "descomponer" en partes que pueden llevar a cabo cada núcleo. En cambio habrá otras que no se puedan descomponer de modo directo, como la escritura y lectura simultánea de un fichero.
Todas estas tareas deben ser coordinadas, esto es, normalmente se programa de forma que uno de los núcleos actúa como maestro y el resto como esclavos: el primero realizará tareas no paralelizables y se encargará de enviar y recibir información que los esclavos necesiten y/o proporcionen. Los cores esclavos se encargarán de realizar la computación respectiva, enviando o recibiendo información vital para el código.
MPI es una alternativa básica a la programación tradicional. Sin embargo, otras aplicaciones están desarrollándose rápidamente, con el afán de simplificar los procesos de comunicación entre cores y la distribución del trabajo: un ejemplo clásico es el OpenMP.
En próximos posts hablaré más en detalle de MPI y OpenMP.
Un saludo
En esta ocasión voy a hablar muy escuetamente de lo que va a ser el futuro de la programación científica: la computación en paralelo.
Actualmente los procesadores tienden a tener un número de núcleos mayor de uno, tal y como era hace unos años: todos estamos familiarizados con los "Pentium core 2 duo" o los modernos "Pentium I3, I5, etc". Sin embargo, por lo que a mi respecta, la formación universitaria en ingeniería sólo cubre nociones básicas de programación en serie. Es obvio que, desde un punto de vista computacional, si se poseen varios núcleos de procesamiento, el realizar una programación capaz de utilizar simultáneamente la potencia de cálculo de esos cores va a ser mucho más ventajoso.
Sin embargo, hasta hace bien poco desconocía por completo un punto de partida que me pudiera orientar para empezar a caminar en este vasto campo. Pues bien, si queréis una piedra angular de referencia, vuestra palabra mágica es MPI.
El MPI se puede considerar como un anexo a los lenguajes de programación típicos en ciencias (Fortran, C) que permite procesar un código arrancando varios "threads", o hilos de computación. Realmente lo ideal es ejecutar tantos threads como núcleos tengas en tu cluster, pero no tiene por que ser así.
Para realizar un código en paralelo, hay que pasar de una mentalidad de cálculo en serie, y comenzar a pensar como se puede distribuir todas las tareas que tiene que hacer el código: habrá tareas, como una búsqueda de elementos en un fichero o matriz, que se puede "descomponer" en partes que pueden llevar a cabo cada núcleo. En cambio habrá otras que no se puedan descomponer de modo directo, como la escritura y lectura simultánea de un fichero.
Todas estas tareas deben ser coordinadas, esto es, normalmente se programa de forma que uno de los núcleos actúa como maestro y el resto como esclavos: el primero realizará tareas no paralelizables y se encargará de enviar y recibir información que los esclavos necesiten y/o proporcionen. Los cores esclavos se encargarán de realizar la computación respectiva, enviando o recibiendo información vital para el código.
MPI es una alternativa básica a la programación tradicional. Sin embargo, otras aplicaciones están desarrollándose rápidamente, con el afán de simplificar los procesos de comunicación entre cores y la distribución del trabajo: un ejemplo clásico es el OpenMP.
En próximos posts hablaré más en detalle de MPI y OpenMP.
Un saludo
miércoles, 1 de febrero de 2012
Introducción sobre elementos finitos II
Hola a todos
Tras un largo parón navideño (de casi 2 meses), quería volver a retomar mi actividad hablando un poco más de los elementos finitos.
Como recordatorio breve, comentar que el método de los elementos finitos es una técnica matemática que permite obtener el cálculo aproximado de la solución de una ecuación determinada sin llegar a resolverla. Este método permite la implementación de la formulación débil de las ecuaciones diferenciales presentes en los modelos físicos: lo que quiero decir es que este método satisface forma integral de nuestro problema de contorno, y no el problema de contorno en si.
El método se basa, esencialmente, en aproximar nuestra solución por una combinación lineal de ciertas funciones que nosotros propondremos y que se denominan "funciones de base". Esto es, si nuestro problema es el siguiente:
L(u)=f
donde u es nuestra variable y L(u) es un operador aplicado a u, aproximamos la solución de u por \hat{u}:
\hat{u} = \Sigma a_i · N_i
donde N_i son las funciones de forma, que a priori son conocidas.
Existen una variedad muy grande de funciones de forma: normalmente, estas misteriosas funciones que aparecen en todos los libros como la panacea para resolver todos nuestros problemas por arte de magia no son más que el fruto del estudio de lo que se denomina teoría de espacios funcionales.
Este campo se dedica a hallar funciones, a priori desconocidas, que cumplan una serie de propiedades que puedan ser útiles para la resolución de problemas matemáticos: cumplir con condiciones de contorno, con condiciones de integrabilidad, topología de dimensión determinada.
En los libros suelen aparecer familias de funciones más o menos complejas: normamente poseen una forma de spline (funciones polinomiales con n puntos conocidos), si bien todas ellas suelen ser de cuadrado integrable: esta condición asegura que la integral de dicha función no va a ser infinito, lo cual es condición necesaria en la resolución de cualquier problema físico real. Ejemplos clásicos pueden ser los polinomios de Lagrange o los polinomios de Hermite. Sin embargo, otras funciones de base pueden ser utilizados en ciertos problemas, como las series de Fourier.
Por lo tanto, conocida la forma de dichas funciones, el objetivo del método, esencialmente, consistirá en buscar aquellos valores que cumplan que:
L(\hat{u}) - f \approx 0
Sabemos que, al no ser u solución exacta, no lo cumple. Como ya dijimos, entonces debe cumplir:
L(\hat{u}) - f = r
Con r, función residuo no nulo de la ecuación. La formulación integral del problema entra aquí. Nosotros vamos a imponer que el valor del residuo en promedio a lo largo de toda la función va a ser cero. Para ello, es necesario integrar el valor de r a lo largo del dominio.
\int w(u) r (u) d\Omega = 0
Esta condición se denomina en matemáticas "Método de los residuos ponderados".
La función w(u) se denomina función peso. La pregunta que surge es obvia: ¿Por que utilizamos esta fórmula con una función misteriosa w(u)? ¿Que motivo real existe para utilizar dicha expresión?. Bien, en el fondo la razón que subyace es la de que queremos que nuestra solución sea una buena solución. El método de elementos finitos busca una SOLUCIÓN APROXIMADA a la ecuación diferencial, cuyas soluciones existen en un subespacio concreto. Nosotros no sabemos cual es en concreto (sólo podemos conocer el ERROR que cometemos entre nuestra aproximación y el valor real), pero si podemos establecer un CRITERIO DE APROXIMACIÓN, esto es, establecer que aproximación puede ser válida para nosotros.
Esto se realiza matemáticamente definiendo la función residuo r(u) y la función peso w(u) de tal modo que cumpla nuestro criterio de aproximación: por ejemplo, si quiero seguir un criterio de mínima distancia geométrica entre mi solución real y la aproximada, definiré mi función residuo como el cuadrado de la diferencia (utilizaría el método de mínimos cuadrados), o si quiero que mi función valga exactamente cero en una serie de puntos, utilizaré unas funciones peso tipo delta de Dirac en esos puntos (método de colocación). Al establecer este criterio de aproximación estamos acotando matemáticamente un subespacio de soluciones válido dentro del subespacio que forman mis funciones de forma. El método se encargará de "ajustar" tus resultados a las condiciones de tu problema de la forma que hayas elegido.
En el caso del método de Galerkin (uno de los métodos más usados), las funciones peso que se escogen son las propias funciones de forma del subespacio que forma mi solución. Esta condición que se establece con el método de residuos ponderados en Galerkin es totalmente lógica: yo impongo que la mejor solución aproximada de la ecuación diferencial a resolver es su PROYECCIÓN ORTOGONAL sobre mi espacio de soluciones que me he definido con las funciones de forma que he escogido (esto matemáticamente se expresa como que su producto escalar debe ser nulo) . De este modo las funciones de forma se escogen para que sean ortogonales al propio error.
Todo este proceso, que en apariencia es relativamente sencillo, a la hora de realizar una implementación numérica (me refiero a utilizar un ordenador para realizar los cálculos) resulta hartamente más complejo. La razón es la propia naturaleza de los ordenadores: no están diseñados para trabajar con variables continuas. Por ello, lo que para nosotros en la teoría es una función continua, al ordenador hay que introducirlo como una serie discreta de puntos que representen dichas funciones. Esto obliga en la práctica a utilizar una notación matricial para resolver el problema.
En el próximo artículo voy a detallar un poco más del proceso de implementación de elementos finitos en un ordenador.
Tras un largo parón navideño (de casi 2 meses), quería volver a retomar mi actividad hablando un poco más de los elementos finitos.
Como recordatorio breve, comentar que el método de los elementos finitos es una técnica matemática que permite obtener el cálculo aproximado de la solución de una ecuación determinada sin llegar a resolverla. Este método permite la implementación de la formulación débil de las ecuaciones diferenciales presentes en los modelos físicos: lo que quiero decir es que este método satisface forma integral de nuestro problema de contorno, y no el problema de contorno en si.
El método se basa, esencialmente, en aproximar nuestra solución por una combinación lineal de ciertas funciones que nosotros propondremos y que se denominan "funciones de base". Esto es, si nuestro problema es el siguiente:
L(u)=f
donde u es nuestra variable y L(u) es un operador aplicado a u, aproximamos la solución de u por \hat{u}:
\hat{u} = \Sigma a_i · N_i
donde N_i son las funciones de forma, que a priori son conocidas.
Existen una variedad muy grande de funciones de forma: normalmente, estas misteriosas funciones que aparecen en todos los libros como la panacea para resolver todos nuestros problemas por arte de magia no son más que el fruto del estudio de lo que se denomina teoría de espacios funcionales.
Este campo se dedica a hallar funciones, a priori desconocidas, que cumplan una serie de propiedades que puedan ser útiles para la resolución de problemas matemáticos: cumplir con condiciones de contorno, con condiciones de integrabilidad, topología de dimensión determinada.
En los libros suelen aparecer familias de funciones más o menos complejas: normamente poseen una forma de spline (funciones polinomiales con n puntos conocidos), si bien todas ellas suelen ser de cuadrado integrable: esta condición asegura que la integral de dicha función no va a ser infinito, lo cual es condición necesaria en la resolución de cualquier problema físico real. Ejemplos clásicos pueden ser los polinomios de Lagrange o los polinomios de Hermite. Sin embargo, otras funciones de base pueden ser utilizados en ciertos problemas, como las series de Fourier.
Por lo tanto, conocida la forma de dichas funciones, el objetivo del método, esencialmente, consistirá en buscar aquellos valores que cumplan que:
L(\hat{u}) - f \approx 0
Sabemos que, al no ser u solución exacta, no lo cumple. Como ya dijimos, entonces debe cumplir:
L(\hat{u}) - f = r
Con r, función residuo no nulo de la ecuación. La formulación integral del problema entra aquí. Nosotros vamos a imponer que el valor del residuo en promedio a lo largo de toda la función va a ser cero. Para ello, es necesario integrar el valor de r a lo largo del dominio.
\int w(u) r (u) d\Omega = 0
Esta condición se denomina en matemáticas "Método de los residuos ponderados".
La función w(u) se denomina función peso. La pregunta que surge es obvia: ¿Por que utilizamos esta fórmula con una función misteriosa w(u)? ¿Que motivo real existe para utilizar dicha expresión?. Bien, en el fondo la razón que subyace es la de que queremos que nuestra solución sea una buena solución. El método de elementos finitos busca una SOLUCIÓN APROXIMADA a la ecuación diferencial, cuyas soluciones existen en un subespacio concreto. Nosotros no sabemos cual es en concreto (sólo podemos conocer el ERROR que cometemos entre nuestra aproximación y el valor real), pero si podemos establecer un CRITERIO DE APROXIMACIÓN, esto es, establecer que aproximación puede ser válida para nosotros.
Esto se realiza matemáticamente definiendo la función residuo r(u) y la función peso w(u) de tal modo que cumpla nuestro criterio de aproximación: por ejemplo, si quiero seguir un criterio de mínima distancia geométrica entre mi solución real y la aproximada, definiré mi función residuo como el cuadrado de la diferencia (utilizaría el método de mínimos cuadrados), o si quiero que mi función valga exactamente cero en una serie de puntos, utilizaré unas funciones peso tipo delta de Dirac en esos puntos (método de colocación). Al establecer este criterio de aproximación estamos acotando matemáticamente un subespacio de soluciones válido dentro del subespacio que forman mis funciones de forma. El método se encargará de "ajustar" tus resultados a las condiciones de tu problema de la forma que hayas elegido.
En el caso del método de Galerkin (uno de los métodos más usados), las funciones peso que se escogen son las propias funciones de forma del subespacio que forma mi solución. Esta condición que se establece con el método de residuos ponderados en Galerkin es totalmente lógica: yo impongo que la mejor solución aproximada de la ecuación diferencial a resolver es su PROYECCIÓN ORTOGONAL sobre mi espacio de soluciones que me he definido con las funciones de forma que he escogido (esto matemáticamente se expresa como que su producto escalar debe ser nulo) . De este modo las funciones de forma se escogen para que sean ortogonales al propio error.
Todo este proceso, que en apariencia es relativamente sencillo, a la hora de realizar una implementación numérica (me refiero a utilizar un ordenador para realizar los cálculos) resulta hartamente más complejo. La razón es la propia naturaleza de los ordenadores: no están diseñados para trabajar con variables continuas. Por ello, lo que para nosotros en la teoría es una función continua, al ordenador hay que introducirlo como una serie discreta de puntos que representen dichas funciones. Esto obliga en la práctica a utilizar una notación matricial para resolver el problema.
En el próximo artículo voy a detallar un poco más del proceso de implementación de elementos finitos en un ordenador.
Suscribirse a:
Entradas (Atom)