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.
Mostrando entradas con la etiqueta Matemáticas. Mostrar todas las entradas
Mostrando entradas con la etiqueta Matemáticas. Mostrar todas las entradas
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
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.
domingo, 20 de noviembre de 2011
Condiciones de contorno
Hola
Voy a escribir un post acerca de una parte fundamental en la resolución de ecuaciones diferenciales. Se trata de las condiciones de contorno.
Como ya sabréis, las condiciones de contorno permiten definir una ecuación diferencial a un caso concreto, esto es, la ecuación pasa de tener infinitas familias de soluciones a tener una sola. Las condiciones de contorno son definidas de esta forma porque típicamente describen el comportamiento de un sistema en la zona donde se acaba el sistema que se quiere analizar: por ejemplo, si quiero describir la transferencia de calor a través de una placa metálica, el volumen en el cual me interesa analizar como fluye el calor utilizando mi ecuación diferencial será una zona determinada que yo me trace, digamos por simplicidad, un rectángulo o cuadrado que abarque todo el espesor de la placa. Se introduce un dibujo para ilustrar el ejemplo:
En la imagen se muestra en linea sólida el contorno físico de nuestra placa metálica, y en linea discontinua el límite de nuestro "volumen de control" que vamos a considerar. Típicamente en la literarura se denomina al dominio que se va a analizar con una ecuación Ω, y a cada uno de las superficies que delimitan dicho volumen se las denomina con la letra Γ . En este caso, se ha denotado como Γ1, Γ2, Γ3 y Γ4 a las superficies que delimitan el contorno.
En este punto uno se plantea la pregunta de cuantas condiciones de contorno debe imponer y que forma tienen que tener. Esta pregunta no tiene una respuesta inmediata: depende de que fenómeno estemos describiendo.
Las ecuaciones diferenciales de la forma $u_t = K \nabla^2 u$, con $\nabla^2 = \frac{\partial}{\partial x_i} $ son denominadas ecuaciones parabólicas. Este tipo de ecuaciones se caracteriza porque debemos definir una condición de contorno por cada derivada que aparezca para cada una de las variables que este derivada.
Por ejemplo, si tenemos $u_t = D(u_{xx} + u_{zz})$, entonces tendremos que definir 4 condiciones de contorno, puesto que tenemos 2 variables espaciales independientes ("x" y "z") derivadas 2 veces. Estas condiciones deben estar dadas de forma que cubran TODAS las fronteras que posee nuestro volumen. Además, la variable tiempo también es una variable independiente del problema y esta derivada una vez. Por lo tanto necesitamos una condición más para esta variable. Típicamente se denomina a las condiciones de contorno para la variable tiempo "condiciones iniciales", puesto que marcan el estado de un sistema en un punto del tiempo, normalmente en el estado inicial de la simulación (t=0). Ejemplos clásicos de este tipo de ecuaciones son la ley de Fourier del calor o la ley de Fick de transferencia de masa.
Las ecuaciones del tipo $u_{tt}=u_{xx}$ se denominan ecuaciones hiperbólicas. Notad que la diferencia con respecto al caso parabólico es que la variable tiempo esta derivada dos veces. Para la resolución de estas ecuaciones sólo son necesarias dos condiciones iniciales.
Las ecuaciones tipo $\nabla^2 u = f(u_i) $ son denominadas ecuaciones elípticas. Estas ecuaciones necesitan definir i condiciones de contorno, esto es, una condición de contorno para cada variable dimensional independiente. Notad que en este tipo de problemas no hay variable tiempo.
¿Que forma deben tener las condiciones de contorno?. En la literatura estan clasificadas en tres tipos:
- Condición de contorno tipo "Dirichlet": Son aquellas que dan el valor de una variable en una de las fronteras. Ej: $T(y=0)=T_0$
- Condición de contorno tipo "Newman": Son aquellas que dan el valor de LA DERIVADA de una variable en una de las fronteras. Ej: $\frac{dT}{dx}(x=0)=0$
- Condición de contorno tipo "Robin" o "mixta": Son aquellas que dan el el valor de LA SUMA del valor de una función y de su derivada. Ej: $T(x=0)+\frac{dT}{dx}(x=0)=0$
Aplicaremos todo lo explicado al ejemplo que se trataba anteriormente. En primer lugar, el problema es un problema bidimensional, puesto que el dibujo es 2D y no se especifica nada de 3 dimensiones. Al ser una ecuación de Fourier, la ecuación es parabólica, por lo que exige definir condiciones de contorno en Γ1, Γ2, Γ3 y Γ4 . Estableceremos el origen de coordenadas en el extremo inferior izquierdo para determinar las condiciones apropiadamente.
$$T_t= K\left( T_{xx}+T{zz} \right)$$
Recordad que nuestro objetivo es hallar T(x,z,t)
Conocemos las temperaturas de la placa superior e inferior. Esto nos da dos condiciones de contorno:
$$T(x,0,t)=T_{cold}$$
$$T(x,H,t)=T_{hot}$$
siendo H el espesor de la placa.
Falta definir 2 condiciones de contorno para los bordes Γ3 y Γ4. No conocemos el valor de la temperatura en dichos contornos, pero si podemos suponer que si nuestra placa se extiende uniformemente en la dirección x, es de esperar que el calor fluya para cualquier x en dirección z. Esto se traduce en que para una misma altura z, el calor va a fluir de la misma forma sea cual sea el punto en el eje x. Eso implica que no existirá variación de la temperatura al movernos en el eje x:
$$\frac{\partial T}{\partial x}(0,z,t)=0$$
$$\frac{\partial T}{\partial x}(L,z,t)=0$$
siendo L la longitud del trozo de placa que hemos tomado.
Puesto que la ley de Fourier depende del tiempo, debemos establecer una condición inicial. Podemos suponer que la temperatura en todo el dominio $\Omega$ sea una temperatura $T_0$ que tiene que ser dada, la cual puede ser una constante o una función.
$$T(x,z,0)=T_0$$
Con esto el problema quedaría cerrado y podría ser resuelto.
Un detalle más. Se trata de la diferencia entre condición de contorno homogenea y condición de contorno inhomogenea. Una condición de contorno es homogenea cuando esta igualada a cero (ej: T(0,y)=0). En caso contrario se denomina inhomogenea (Ej: T(0,y)=f(x,y)).
Un punto imporante que debe ser entendido es que las condiciones de contorno inhomogeneas introducen UNA FUENTE DE EXCITACIÓN EN EL SISTEMA. Esto es, cuando se dice que la temperatura en una frontera toma un valor $$T_0$$, estoy asumiendo que estoy introduciendo calor al sistema en ese punto. Esto de un modo u otro va a suponer que el sistema RESPONDERÁ ante dicha entrada con un comportamiento dado, que en este ejemplo será que el perfil de temperaturas en el volumen de control será uno específico. Una condición de contorno homogenea NO introduce excitación en el problema.
Voy a escribir un post acerca de una parte fundamental en la resolución de ecuaciones diferenciales. Se trata de las condiciones de contorno.
Como ya sabréis, las condiciones de contorno permiten definir una ecuación diferencial a un caso concreto, esto es, la ecuación pasa de tener infinitas familias de soluciones a tener una sola. Las condiciones de contorno son definidas de esta forma porque típicamente describen el comportamiento de un sistema en la zona donde se acaba el sistema que se quiere analizar: por ejemplo, si quiero describir la transferencia de calor a través de una placa metálica, el volumen en el cual me interesa analizar como fluye el calor utilizando mi ecuación diferencial será una zona determinada que yo me trace, digamos por simplicidad, un rectángulo o cuadrado que abarque todo el espesor de la placa. Se introduce un dibujo para ilustrar el ejemplo:
En este punto uno se plantea la pregunta de cuantas condiciones de contorno debe imponer y que forma tienen que tener. Esta pregunta no tiene una respuesta inmediata: depende de que fenómeno estemos describiendo.
Las ecuaciones diferenciales de la forma $u_t = K \nabla^2 u$, con $\nabla^2 = \frac{\partial}{\partial x_i} $ son denominadas ecuaciones parabólicas. Este tipo de ecuaciones se caracteriza porque debemos definir una condición de contorno por cada derivada que aparezca para cada una de las variables que este derivada.
Por ejemplo, si tenemos $u_t = D(u_{xx} + u_{zz})$, entonces tendremos que definir 4 condiciones de contorno, puesto que tenemos 2 variables espaciales independientes ("x" y "z") derivadas 2 veces. Estas condiciones deben estar dadas de forma que cubran TODAS las fronteras que posee nuestro volumen. Además, la variable tiempo también es una variable independiente del problema y esta derivada una vez. Por lo tanto necesitamos una condición más para esta variable. Típicamente se denomina a las condiciones de contorno para la variable tiempo "condiciones iniciales", puesto que marcan el estado de un sistema en un punto del tiempo, normalmente en el estado inicial de la simulación (t=0). Ejemplos clásicos de este tipo de ecuaciones son la ley de Fourier del calor o la ley de Fick de transferencia de masa.
Las ecuaciones del tipo $u_{tt}=u_{xx}$ se denominan ecuaciones hiperbólicas. Notad que la diferencia con respecto al caso parabólico es que la variable tiempo esta derivada dos veces. Para la resolución de estas ecuaciones sólo son necesarias dos condiciones iniciales.
Las ecuaciones tipo $\nabla^2 u = f(u_i) $ son denominadas ecuaciones elípticas. Estas ecuaciones necesitan definir i condiciones de contorno, esto es, una condición de contorno para cada variable dimensional independiente. Notad que en este tipo de problemas no hay variable tiempo.
¿Que forma deben tener las condiciones de contorno?. En la literatura estan clasificadas en tres tipos:
- Condición de contorno tipo "Dirichlet": Son aquellas que dan el valor de una variable en una de las fronteras. Ej: $T(y=0)=T_0$
- Condición de contorno tipo "Newman": Son aquellas que dan el valor de LA DERIVADA de una variable en una de las fronteras. Ej: $\frac{dT}{dx}(x=0)=0$
- Condición de contorno tipo "Robin" o "mixta": Son aquellas que dan el el valor de LA SUMA del valor de una función y de su derivada. Ej: $T(x=0)+\frac{dT}{dx}(x=0)=0$
Aplicaremos todo lo explicado al ejemplo que se trataba anteriormente. En primer lugar, el problema es un problema bidimensional, puesto que el dibujo es 2D y no se especifica nada de 3 dimensiones. Al ser una ecuación de Fourier, la ecuación es parabólica, por lo que exige definir condiciones de contorno en Γ1, Γ2, Γ3 y Γ4 . Estableceremos el origen de coordenadas en el extremo inferior izquierdo para determinar las condiciones apropiadamente.
$$T_t= K\left( T_{xx}+T{zz} \right)$$
Recordad que nuestro objetivo es hallar T(x,z,t)
Conocemos las temperaturas de la placa superior e inferior. Esto nos da dos condiciones de contorno:
$$T(x,0,t)=T_{cold}$$
$$T(x,H,t)=T_{hot}$$
siendo H el espesor de la placa.
Falta definir 2 condiciones de contorno para los bordes Γ3 y Γ4. No conocemos el valor de la temperatura en dichos contornos, pero si podemos suponer que si nuestra placa se extiende uniformemente en la dirección x, es de esperar que el calor fluya para cualquier x en dirección z. Esto se traduce en que para una misma altura z, el calor va a fluir de la misma forma sea cual sea el punto en el eje x. Eso implica que no existirá variación de la temperatura al movernos en el eje x:
$$\frac{\partial T}{\partial x}(0,z,t)=0$$
$$\frac{\partial T}{\partial x}(L,z,t)=0$$
siendo L la longitud del trozo de placa que hemos tomado.
Puesto que la ley de Fourier depende del tiempo, debemos establecer una condición inicial. Podemos suponer que la temperatura en todo el dominio $\Omega$ sea una temperatura $T_0$ que tiene que ser dada, la cual puede ser una constante o una función.
$$T(x,z,0)=T_0$$
Con esto el problema quedaría cerrado y podría ser resuelto.
Un detalle más. Se trata de la diferencia entre condición de contorno homogenea y condición de contorno inhomogenea. Una condición de contorno es homogenea cuando esta igualada a cero (ej: T(0,y)=0). En caso contrario se denomina inhomogenea (Ej: T(0,y)=f(x,y)).
Un punto imporante que debe ser entendido es que las condiciones de contorno inhomogeneas introducen UNA FUENTE DE EXCITACIÓN EN EL SISTEMA. Esto es, cuando se dice que la temperatura en una frontera toma un valor $$T_0$$, estoy asumiendo que estoy introduciendo calor al sistema en ese punto. Esto de un modo u otro va a suponer que el sistema RESPONDERÁ ante dicha entrada con un comportamiento dado, que en este ejemplo será que el perfil de temperaturas en el volumen de control será uno específico. Una condición de contorno homogenea NO introduce excitación en el problema.
sábado, 19 de noviembre de 2011
Sobre la divergencia y el rotacional (II)
Hola a todos
Voy a realizar alguna aclaración más acerca de estos dos operadores, de gran importancia en el mundo de la ciencia y la ingeniería como la divergencia y el rotacional. La idea me ha surgido al realizar unos ejercicios de comprensión sobre el significado de estos, y pienso que será util para aquellos que tenga que tratar con ellos o bien durante su formación o bien en su trabajo.
Se comentó en el anterior post que la divergencia tiene que ver con un balance de flujo a través de una superficie. Quizá lo más riguroso sería acudir a su definición matemática estricta:
$$\nabla \cdot F = \frac{1}{\bigtriangleup V}\int_S \vec{F} \cdot \vec{n} \cdot dS $$
Como se puede observar, la divergencia se define como el flujo de campo neto que atraviesa una superficie asociada a un volumen determinado. Vamos a poner un ejemplo de cual sería la divergencia para un campo concreto. En concreto, utilizaremos un campo en coordenadas esféricas:
$\vec{F}=r^2 \vec{r}$
Aplicando la definición de divergencia para un sistema de referencia de coordenadas esféricas, se puede obtener el valor exacto:
$\nabla \cdot \vec{F}=\frac{1}{r^2}\frac{\partial}{\partial r} \left( r^2 \vec{F} \right) = 4r $
El resultado nos dice que la divergencia es siempre positiva en todo el campo, salvo en el centro que debe ser nula. Efectivamente, cuanto más nos acercamos las lineas de campo son cada vez menores.
Se puede analizar las zonas del campo en las cuales la divergencia es positiva, negativa o nula. Recordad que la divergencia tiene que ver con un balance en un volumen determinado. Como se puede observar, la divergencia es siempre positiva en todo el campo. Conforme nos alejamos del origen de coordenadas, el flujo a través de las distintas superficies que englobarían volúmenes de tipo esférico (escogemos estos por ser un sistema de coordenadas esféricas) dentro del campo (en el plano 2D serían circunferencias concéntricas) cada vez es mayor.
Al representar gráficamente un detalle de este campo y pintar las lineas de campo a una distancia r y otra r +dr, se observa que efectivamente las lineas de campo son mayores al alejarnos del centro y por lo tanto mayores en la superficie S+dS que en la superficie S. La divergencia aplicada a este diferencial de volumen que abarca ambas superficies debe ser mayor que cero.
Pero este balance visual puede llevar a error algunas veces, por lo que hay que ser precavido a la hora de evaluar. Lo ilustraremos ahora realizando el mismo análisis con otro campo distinto:
$\vec{F}=\frac{1}{r^2} \vec{r}$
La representación del campo es:
Como se aprecia, las lineas de campo son decrecientes conforme nos alejamos del centro. Si hacemos un análisis visual basándonos en el ejemplo anterior, a priori podría pensarse que la divergencia de este campo es siempre menor que cero para cualquier r, haciéndose 0 en el infinito. Si ampliamos en imagen un detalle en una zona del campo:
Parece que las lines de campo son menores en S+dS que en S en todo el plano, por lo que el balance sería menor que cero. Sin embargo, esto no es así. Calculando la divergencia analíticamente:
$\nabla \cdot \vec{F}=\frac{1}{r^2}\frac{\partial}{\partial r} \left( r^2 \vec{F} \right) = 0 $
El resultado del cálculo analítico arroja un resultado 0 para cualquier punto del campo. ¿Que es lo que ha fallado en nuestra deducción?
La respuesta a esta pregunta no es trivial. Si observamos de nuevo la definición de divergencia, se puede apreciar que es la integral del flujo superficial a través de una superficie. Al escoger un volumen como el del detalle anterior, se observa que el valor del campo es mayor en un punto de la superficie S que en una posición de la superficie S+dS (porque al tener mayor radio, el valor del campo es menor). En concreto la variación del campo es inversamente proporcional a $s^2$. Pero hay que destacar que LA SUPERFICIE A TRAVÉS DE LA QUE CIRCULA EL FLUJO DEL CAMPO TAMBIÉN ES MAYOR EN S+dS QUE EN S. Concretamente, es mayor proporcionalmente con el cuadrado de la distancia.
Se tiene en conjunto que, LA DISMINUCIÓN DEL VALOR DEL CAMPO AL AUMENTAR EL RADIO SE VE COMPENSADA CON EL AUMENTO DE SUPERFICIE AL AUMENTAR DICHO RADIO. Digamos que, aunque F es más intenso en la superficie S que en S+dS, la superficie que atravesa es menor que la de S+dS, compensando así los valores en ambos puntos. En este caso, esta compensación es exactamente igual, y por tanto el flujo a través de ambas superficies es exactamente el mismo. El balance neto entonces es exactamente cero.
Por lo tanto, para resumir:
- La divergencia de un campo vectorial es un escalar.
- La divergencia implica un BALANCE NETO de flujo en un volumen determinado a través de su superficie.
- Al calcular la divergencia en un volumen dado (por ejemplo un rectángulo, o un sector circular) hay que tener en cuenta tanto el valor de campo en cada superficie como el tamaño de las superficies que atraviesa. Recordar que es un flujo superficial.
El próximo post introduciré un ejemplo similar para el operador rotacional
Voy a realizar alguna aclaración más acerca de estos dos operadores, de gran importancia en el mundo de la ciencia y la ingeniería como la divergencia y el rotacional. La idea me ha surgido al realizar unos ejercicios de comprensión sobre el significado de estos, y pienso que será util para aquellos que tenga que tratar con ellos o bien durante su formación o bien en su trabajo.
Se comentó en el anterior post que la divergencia tiene que ver con un balance de flujo a través de una superficie. Quizá lo más riguroso sería acudir a su definición matemática estricta:
$$\nabla \cdot F = \frac{1}{\bigtriangleup V}\int_S \vec{F} \cdot \vec{n} \cdot dS $$
Como se puede observar, la divergencia se define como el flujo de campo neto que atraviesa una superficie asociada a un volumen determinado. Vamos a poner un ejemplo de cual sería la divergencia para un campo concreto. En concreto, utilizaremos un campo en coordenadas esféricas:
$\vec{F}=r^2 \vec{r}$
Aplicando la definición de divergencia para un sistema de referencia de coordenadas esféricas, se puede obtener el valor exacto:
$\nabla \cdot \vec{F}=\frac{1}{r^2}\frac{\partial}{\partial r} \left( r^2 \vec{F} \right) = 4r $
El resultado nos dice que la divergencia es siempre positiva en todo el campo, salvo en el centro que debe ser nula. Efectivamente, cuanto más nos acercamos las lineas de campo son cada vez menores.
Se puede analizar las zonas del campo en las cuales la divergencia es positiva, negativa o nula. Recordad que la divergencia tiene que ver con un balance en un volumen determinado. Como se puede observar, la divergencia es siempre positiva en todo el campo. Conforme nos alejamos del origen de coordenadas, el flujo a través de las distintas superficies que englobarían volúmenes de tipo esférico (escogemos estos por ser un sistema de coordenadas esféricas) dentro del campo (en el plano 2D serían circunferencias concéntricas) cada vez es mayor.
Al representar gráficamente un detalle de este campo y pintar las lineas de campo a una distancia r y otra r +dr, se observa que efectivamente las lineas de campo son mayores al alejarnos del centro y por lo tanto mayores en la superficie S+dS que en la superficie S. La divergencia aplicada a este diferencial de volumen que abarca ambas superficies debe ser mayor que cero.
$\vec{F}=\frac{1}{r^2} \vec{r}$
La representación del campo es:
$\nabla \cdot \vec{F}=\frac{1}{r^2}\frac{\partial}{\partial r} \left( r^2 \vec{F} \right) = 0 $
El resultado del cálculo analítico arroja un resultado 0 para cualquier punto del campo. ¿Que es lo que ha fallado en nuestra deducción?
La respuesta a esta pregunta no es trivial. Si observamos de nuevo la definición de divergencia, se puede apreciar que es la integral del flujo superficial a través de una superficie. Al escoger un volumen como el del detalle anterior, se observa que el valor del campo es mayor en un punto de la superficie S que en una posición de la superficie S+dS (porque al tener mayor radio, el valor del campo es menor). En concreto la variación del campo es inversamente proporcional a $s^2$. Pero hay que destacar que LA SUPERFICIE A TRAVÉS DE LA QUE CIRCULA EL FLUJO DEL CAMPO TAMBIÉN ES MAYOR EN S+dS QUE EN S. Concretamente, es mayor proporcionalmente con el cuadrado de la distancia.
Se tiene en conjunto que, LA DISMINUCIÓN DEL VALOR DEL CAMPO AL AUMENTAR EL RADIO SE VE COMPENSADA CON EL AUMENTO DE SUPERFICIE AL AUMENTAR DICHO RADIO. Digamos que, aunque F es más intenso en la superficie S que en S+dS, la superficie que atravesa es menor que la de S+dS, compensando así los valores en ambos puntos. En este caso, esta compensación es exactamente igual, y por tanto el flujo a través de ambas superficies es exactamente el mismo. El balance neto entonces es exactamente cero.
Por lo tanto, para resumir:
- La divergencia de un campo vectorial es un escalar.
- La divergencia implica un BALANCE NETO de flujo en un volumen determinado a través de su superficie.
- Al calcular la divergencia en un volumen dado (por ejemplo un rectángulo, o un sector circular) hay que tener en cuenta tanto el valor de campo en cada superficie como el tamaño de las superficies que atraviesa. Recordar que es un flujo superficial.
El próximo post introduciré un ejemplo similar para el operador rotacional
domingo, 6 de noviembre de 2011
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.
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.
domingo, 30 de octubre de 2011
Metodos perturbativos: Brevísima introducción
Hola a todos
Hoy voy a hablar de un tema que me ha interesado mucho y que creo es importante conocer por aquellos que deseen aprender e introducirse en el mundo de la ciencia de vanguardia: Me refiero a los metodos perturbativos.
Estas tecnicas matematicas fueron desarrolladas aproximadamente durante el siglo pasado por cientificos que, o bien se dedicaban a la matematica de los sistemas dinamicos, o bien trabajaban en areas de la fisica o la ingenieria donde debian tratar a diario con ecuaciones o modelos matematicos no resolubles analiticamente (como por ejemplo la mecanica de fluidos).
Los metodos perturbativos, tambien llamados asintoticos, tratan de obtener informacion sobre el sistema de ecuaciones a resolver sin proceder a su resolucion. Para ello estas tecnicas introducen en las ecuaciones "soluciones perturbadas": se presupone una solucion para la ecuacion que depende de una variable pequeña epsilon y otros parametros deseados. La forma de estas soluciones varia segun el tipo de perturbacion que se desee introducir en el sistema: pued ser una serie de potencias de epsilon, una funcion exponencial u otro tipo de funcion analitica dependiente de epsilon.
Puesto que la solucion propuesta es efectivamente, una solucion de la ecuacion, satisfara el sistema de ecuaciones. La sustitucion de la misma en las ecuaciones dara lugar a una serie de ecuaciones en funcion de los parametros incluidos en nuetra solucion, las cuales seran estudiadas matematicamente utilizando teoria de dinamica no lineal, esto es, se procedera al estudio de su estabilidad (convergencia de la solucion a tiempo infinito, soluciones oscilatorias, etc.)
Como se puede observar, los metodos perturbativos son de extremada utilidad, si bien su aplicacion no resulta tan inmediata como parece. En los casos de estudio de ecuaciones diferenciales con condiciones de contorno, este proceso es cuanto menos pesado, si no cuasi ininteligible para aquellos no muy acostumbrados a trabajar de forma fluida con este tipo de metodos. Es por ello por lo que don pocos los expertos sobre esta materia.
En proximos post hablare un poco mas en detalle sobre este tipo de metodos, para que veais como se aplican en la practica.
Hoy voy a hablar de un tema que me ha interesado mucho y que creo es importante conocer por aquellos que deseen aprender e introducirse en el mundo de la ciencia de vanguardia: Me refiero a los metodos perturbativos.
Estas tecnicas matematicas fueron desarrolladas aproximadamente durante el siglo pasado por cientificos que, o bien se dedicaban a la matematica de los sistemas dinamicos, o bien trabajaban en areas de la fisica o la ingenieria donde debian tratar a diario con ecuaciones o modelos matematicos no resolubles analiticamente (como por ejemplo la mecanica de fluidos).
Los metodos perturbativos, tambien llamados asintoticos, tratan de obtener informacion sobre el sistema de ecuaciones a resolver sin proceder a su resolucion. Para ello estas tecnicas introducen en las ecuaciones "soluciones perturbadas": se presupone una solucion para la ecuacion que depende de una variable pequeña epsilon y otros parametros deseados. La forma de estas soluciones varia segun el tipo de perturbacion que se desee introducir en el sistema: pued ser una serie de potencias de epsilon, una funcion exponencial u otro tipo de funcion analitica dependiente de epsilon.
Puesto que la solucion propuesta es efectivamente, una solucion de la ecuacion, satisfara el sistema de ecuaciones. La sustitucion de la misma en las ecuaciones dara lugar a una serie de ecuaciones en funcion de los parametros incluidos en nuetra solucion, las cuales seran estudiadas matematicamente utilizando teoria de dinamica no lineal, esto es, se procedera al estudio de su estabilidad (convergencia de la solucion a tiempo infinito, soluciones oscilatorias, etc.)
Como se puede observar, los metodos perturbativos son de extremada utilidad, si bien su aplicacion no resulta tan inmediata como parece. En los casos de estudio de ecuaciones diferenciales con condiciones de contorno, este proceso es cuanto menos pesado, si no cuasi ininteligible para aquellos no muy acostumbrados a trabajar de forma fluida con este tipo de metodos. Es por ello por lo que don pocos los expertos sobre esta materia.
En proximos post hablare un poco mas en detalle sobre este tipo de metodos, para que veais como se aplican en la practica.
viernes, 28 de octubre de 2011
Magnitudes escalares y vectoriales
Hola a todos
En este post quiero realizar una serie de aclaraciones que personalmente me trajeron quebraderos de cabeza cuando empecé a introducirme seriamente en el mundo de las matemáticas. Hablo de los conceptos de magnitud escalar y magnitud vectorial.
Lo normal y típico cuanto alguien termina sus estudios en ingeniería (yo lo pensaba) es que sepa que la diferencia entre ambas es la presencia de dirección y sentido en las magnitudes vectoriales. Eso es correcto, pero a la hora de trabajar con ecuaciones no tiene mucha utilidad. Voy a explicar un poco más en detalle sus diferencias, e incluiré un ejemplo de cada una relacionado con el campo de mecánica de fluidos.
Magnitud escalar: Las magnitudes escalares son aquellas que se pueden definir matemáticamente como una función de las variables que intervienen en la descripción de un sistema. Son funciones que, evaluadas en un punto del dominio en el que estan definidas (el "volumen de control" que se estudia típicamente en ingeniería) dan como resultado un número (2, 0, 3E-6, etc.). Ejemplos típicos de magnitudes escalares son la temperatura o la presión.
Ejemplo: Presión en un punto de un volumen de control
$$ P(x,y,z)=4x^2+3y^2-8.4xy-2z+1 $$
Como se observa, la presión es una función de las coordenadas espaciales z, x e y. Sólo depende de la posición que se considere.
Magnitud vectorial: Son aquellas que para su definición deben ser expresadas en lo que conocemos como base vectorial, esto es, estan referenciadas a un sistema de coordenadas direccional. Estas magnitudes poseen distintos valores para un mismo punto dependiendo de la dirección que se considere (es lo que se denomina componentes del vector). Cada componente será a su vez una expresión que dependerá de las variables del sistema. Notad que, al necesitar una base vectorial que defina dichas magnitudes, aparece de forma natural el concepto de "dirección" y "sentido" de la magnitud.
Ejemplo: Velocidad de una partícula fluida en un punto del volumen de control
$$u=[u_x,u_y,u_z]=\left[4x^2-9y+z,3z^4-2,2x+4y^3\right]$$
Como se puede observar, el vector u posee tres componentes, una por cada dirección espacial. Además cada componente es a su vez una función de z,x e y. La dirección final de la velocidad de la partícula fluida en cualquier punto será la suma de ambas componentes en forma vectorial, o la raíz cuadrada del cuadrado de cada componente como escalar, más su vector director, que será aquel vector unitario que indique la dirección y sentido del vector u.
Quizá esta diferencia es obvia, pero es la que marca la diferencia en el cálculo diferencial o tensorial de los distintos operadores: no es lo mismo el laplaciano de una magnitud escalar que el laplaciano de una magnitud vectorial.
En el próximo post escribiré una pequeña tablita con los resultados para distintas magnitudes vectoriales y escalares.
En este post quiero realizar una serie de aclaraciones que personalmente me trajeron quebraderos de cabeza cuando empecé a introducirme seriamente en el mundo de las matemáticas. Hablo de los conceptos de magnitud escalar y magnitud vectorial.
Lo normal y típico cuanto alguien termina sus estudios en ingeniería (yo lo pensaba) es que sepa que la diferencia entre ambas es la presencia de dirección y sentido en las magnitudes vectoriales. Eso es correcto, pero a la hora de trabajar con ecuaciones no tiene mucha utilidad. Voy a explicar un poco más en detalle sus diferencias, e incluiré un ejemplo de cada una relacionado con el campo de mecánica de fluidos.
Magnitud escalar: Las magnitudes escalares son aquellas que se pueden definir matemáticamente como una función de las variables que intervienen en la descripción de un sistema. Son funciones que, evaluadas en un punto del dominio en el que estan definidas (el "volumen de control" que se estudia típicamente en ingeniería) dan como resultado un número (2, 0, 3E-6, etc.). Ejemplos típicos de magnitudes escalares son la temperatura o la presión.
Ejemplo: Presión en un punto de un volumen de control
$$ P(x,y,z)=4x^2+3y^2-8.4xy-2z+1 $$
Como se observa, la presión es una función de las coordenadas espaciales z, x e y. Sólo depende de la posición que se considere.
Magnitud vectorial: Son aquellas que para su definición deben ser expresadas en lo que conocemos como base vectorial, esto es, estan referenciadas a un sistema de coordenadas direccional. Estas magnitudes poseen distintos valores para un mismo punto dependiendo de la dirección que se considere (es lo que se denomina componentes del vector). Cada componente será a su vez una expresión que dependerá de las variables del sistema. Notad que, al necesitar una base vectorial que defina dichas magnitudes, aparece de forma natural el concepto de "dirección" y "sentido" de la magnitud.
Ejemplo: Velocidad de una partícula fluida en un punto del volumen de control
$$u=[u_x,u_y,u_z]=\left[4x^2-9y+z,3z^4-2,2x+4y^3\right]$$
Como se puede observar, el vector u posee tres componentes, una por cada dirección espacial. Además cada componente es a su vez una función de z,x e y. La dirección final de la velocidad de la partícula fluida en cualquier punto será la suma de ambas componentes en forma vectorial, o la raíz cuadrada del cuadrado de cada componente como escalar, más su vector director, que será aquel vector unitario que indique la dirección y sentido del vector u.
Quizá esta diferencia es obvia, pero es la que marca la diferencia en el cálculo diferencial o tensorial de los distintos operadores: no es lo mismo el laplaciano de una magnitud escalar que el laplaciano de una magnitud vectorial.
En el próximo post escribiré una pequeña tablita con los resultados para distintas magnitudes vectoriales y escalares.
miércoles, 26 de octubre de 2011
Escribe Latex en tu blog
Hola a todos
Para todos aquellos que utiliceis Latex en vuestros escritos científicos, pero no sabeis como implementar dicha escritura en vuestro blog, os remito a esta web que he encontrado que te explica como hacerlo (al menos en blogger/blogspot).
http://watchmath.com/vlog/?p=438
El proceso es muy sencillo: solo es copiar y pegar el script que deja el autor del blog tal y como se indica en su web, si bien hay que realizarun cambio en la dirección web que esta incluida en el script:
cambiar
por
Y ya está. A disfrutar de vuestra creatividad.
Para los que no sabeis que es Latex, os muestro un ejemplo para que veais de que va el tema. Y si os interesa, solo teneis que buscar información sobre el mismo en la web. Hay cientos de manuales y foros que hablan sobre Latex. Hace falta un poco de esfuerzo para adaptarse a la escritura, pero una vez controlas los comandos la escritura es muy rápida e intuitiva. Yo he dejado de usar word !!
$$ F(x,y)=\int_{a}^{b} g(x,y) dx $$
$$\nabla \cdot \gamma = \frac{\partial \gamma}{\partial x}+ \frac{\partial \gamma}{\partial y} $$
ATENCIÓN!! HA HABIDO CAMBIO DE SERVIDOR. ESTA DIRECCIÓN YA NO FUNCIONA. MIRAR POST POSTERIORES PARA ENCONTRAR ALTERNATIVA PARA ESCRITURA CON LATEX
Para todos aquellos que utiliceis Latex en vuestros escritos científicos, pero no sabeis como implementar dicha escritura en vuestro blog, os remito a esta web que he encontrado que te explica como hacerlo (al menos en blogger/blogspot).
http://watchmath.com/vlog/?p=438
El proceso es muy sencillo: solo es copiar y pegar el script que deja el autor del blog tal y como se indica en su web, si bien hay que realizarun cambio en la dirección web que esta incluida en el script:
cambiar
http://www.watchmath.com/cgi-bin/mathtex3.js
por
http://www.watchmath.com/main/cgi-bin/mathtex3.js
Y ya está. A disfrutar de vuestra creatividad.
Para los que no sabeis que es Latex, os muestro un ejemplo para que veais de que va el tema. Y si os interesa, solo teneis que buscar información sobre el mismo en la web. Hay cientos de manuales y foros que hablan sobre Latex. Hace falta un poco de esfuerzo para adaptarse a la escritura, pero una vez controlas los comandos la escritura es muy rápida e intuitiva. Yo he dejado de usar word !!
$$ F(x,y)=\int_{a}^{b} g(x,y) dx $$
$$\nabla \cdot \gamma = \frac{\partial \gamma}{\partial x}+ \frac{\partial \gamma}{\partial y} $$
ATENCIÓN!! HA HABIDO CAMBIO DE SERVIDOR. ESTA DIRECCIÓN YA NO FUNCIONA. MIRAR POST POSTERIORES PARA ENCONTRAR ALTERNATIVA PARA ESCRITURA CON LATEX
domingo, 23 de octubre de 2011
Sobre la divergencia y el rotacional (I)
Para aquellos que estudien matemática tensorial, seguramente se hayan topado con algunos de estos dos operadores matemáticos. Sin embargo cuando se habla con un sentido físico se observa que estos operadores poseen cierto significado físico:
- La divergencia esta siempre relacionado con la componente normal de un flujo (esto es, variación de magnitud en el tiempo por unidad de area). Si se hace un análisis de lineas de corriente de flujo se puede observar que si la divergencia es cero, las lineas de corriente que cruzan un volumen de control se conservan: esto es: entran el mismo número de lineas al volumen que las que salen. Algún ejemplo:
· En electromagnetismo, un imán genera un campo magnético que se puede representar con un campo vectorial. Si se toma un volumen de control alrededor de dicho imán se observa que, puesto que el propio imán es la fuente del campo, y que según las leyes del electromagnetismo, la divergencia del campo magnético es cero, se deduce que las lineas de corriente deben ser cerradas (esto es, nacen en un punto del imán pero deben morir en otro punto del mismo), de forma que la ecuación de la divergencia se cumpla.
· En mecánica de fluidos, la divergencia nula sobre un volumen de control implica la conservación de la masa en el volumen implicado. En fluidos cuasi incompresibles como el agua, implica también conservación de caudal.
Nota: La divergencia de un campo vectorial es un ESCALAR.
- El rotacional esta relacionado con la componente tangencial de un flujo. Este operador da como resultado un vector tangente a los implicados en dicho operador.
Ejemplo:
· En mecánica de fluidos, el rotacional del campo de velocidades del flujo esta relacionado con la circulación y la vorticidad del fluido.
Nota: El resultado del rotacional de dos campos vectoriales es otro campo vectorial perpendicular a los anteriores.
- La divergencia esta siempre relacionado con la componente normal de un flujo (esto es, variación de magnitud en el tiempo por unidad de area). Si se hace un análisis de lineas de corriente de flujo se puede observar que si la divergencia es cero, las lineas de corriente que cruzan un volumen de control se conservan: esto es: entran el mismo número de lineas al volumen que las que salen. Algún ejemplo:
· En electromagnetismo, un imán genera un campo magnético que se puede representar con un campo vectorial. Si se toma un volumen de control alrededor de dicho imán se observa que, puesto que el propio imán es la fuente del campo, y que según las leyes del electromagnetismo, la divergencia del campo magnético es cero, se deduce que las lineas de corriente deben ser cerradas (esto es, nacen en un punto del imán pero deben morir en otro punto del mismo), de forma que la ecuación de la divergencia se cumpla.
· En mecánica de fluidos, la divergencia nula sobre un volumen de control implica la conservación de la masa en el volumen implicado. En fluidos cuasi incompresibles como el agua, implica también conservación de caudal.
Nota: La divergencia de un campo vectorial es un ESCALAR.
- El rotacional esta relacionado con la componente tangencial de un flujo. Este operador da como resultado un vector tangente a los implicados en dicho operador.
Ejemplo:
· En mecánica de fluidos, el rotacional del campo de velocidades del flujo esta relacionado con la circulación y la vorticidad del fluido.
Nota: El resultado del rotacional de dos campos vectoriales es otro campo vectorial perpendicular a los anteriores.
Suscribirse a:
Entradas (Atom)