Hola a todos
En este artículo voy a hablar sobre programación avanzada en FORTRAN. Cuando me refiero a programación "avanzada" en FORTRAN quiero remarcar implícitamente que las estructuras sintácticas utilizadas para realizar nuestro programa no van a ser aquellas que permitan ejecutar correctamente el programa, sino aquellas que pueden ofrecer grandes ventajas en cuanto a TIEMPO DE COMPUTACIÓN Y ALMACENAMIENTO EN MEMORIA.
Una estructura de computación avanzada es extremadamente útil cuando se necesita elaborar un programa que por sus características requiere de un esfuerzo computacional muy elevado. Ejemplos clásicos de este tipo de programas pueden ser los programas CFD (donde los solvers utilizados requieren una malla de puntos muy extensa para llevar a cabo la computación de la ecuación a resolver en el dominio estudiado) o programas que manipulen listas o matrices de datos MUY grandes.
Dejando a un lado el caso de que seas un experto programador (en cuyo caso lo que voy a decir en este artículo seguramente ya lo sabes puesto que de facto yo no soy tal), a la hora de computar un programa no siempre conviene diseñar EL MEJOR PROGRAMA, sino aquel programa que te de lo que necesitas en el menor tiempo posible DESDE QUE COMIENZAS A DISEÑARLO HASTA QUE OBTIENES LA RESPUESTA: de nada sirve elaborar un programa extremadamente eficiente que ahorra un 12% tu tiempo de computación si has empleado un esfuerzo en depurarlo y modificarlo al diseñarlo que sobrepasa con creces dicha cantidad de tiempo.
A continuación anotaré una serie de consejos en cuanto al manejo de FORTRAN buscando una estructura avanzada de programación:
1 - Conseguir elaborar tu programa intentando optimizar estos requisitos es una tarea hartamente complicada si no conoces casi a fondo las estructuras sintácticas que permite FORTRAN y las funciones intrínsecas que tiene implementadas. Échale un vistazo a las funciones intrínsecas que tiene implementado el programa y siempre que puedas utilízalas en vez de hacerte unas propias.
2 - Evita el uso de bucles DO siempre que puedas, sustituyéndolas por funciones intrínsecas. Por ejemplo en vez de recorrer una matriz elemento a elemento y realizar una operación específica, utiliza notación matricial que afecta simultáneamente a todos los elementos de la misma y es mucho más rápida.
Ejemplo: A = A +5, donde A es una matriz de tamaño M x N.
3 - Si quieres filtrar ciertos elementos de la matriz que te interesan y que son fácilmente distinguibles del resto (por ejemplo, elementos menores que 0 o que cumplan algún requisito lógico), la declaración WHERE te permite seleccionarlos instantáneamente mucho más rápidamente.
4 - FORTRAN odia los mapeados no lineales de elementos, prefiriendo seleccionar datos de forma lineal. Ejemplo: Si b1 = [1 3 8] y b2 = [2 4 6] son dos vectores con las posiciones de ciertos elementos dentro de otro vector más grande, entonces FORTRAN tarda sensiblemente más en acceder a A(b1) que a A(b2), puesto que el acceso a b2 se puede introducir en notación intrínseca de FORTRAN como A(2:6:2).
5 - Utiliza los punteros. Para los que no saben que es esto (yo lo desconocía hasta hace poco), un puntero es un objeto que no tiene una forma predefinida, esto es, es un objeto que puede tomar una dimensión variable dentro de su clase. Lo que quiero decir es que un puntero es un "alias" del objeto al que apunta, de tal forma que cuando hago referencia al puntero estoy accediendo en memoria física al objeto que se le ha asignado. Ejemplo: si se apunta con un puntero llamado Punt a un array a1=[1 0 25 0 0]
Punt => a1
Entonces si se referencia Punt(3), entonces FORTRAN toma el valor 25. Pero a continuación se puede reasignar Punt a otro vector, b1=[2 3 4 5 6 7 8 9 0]. Ahora, Punt(3) = 4. Como se puede ver, los punteros son extremadamente eficientes cuando se realizan actualizaciones de matrices grandes a partir de datos del paso anterior, ya que estos me permiten en cada paso calcular el nuevo valor de mi matriz a partir de los datos anteriores y asignárselo al puntero (lo que cambio es la dirección en memoria) en vez de borrar la matriz del paso anterior y volver a copiarla de la ya calculada en el paso actual.
6 - Graba en tu fichero de datos sólo los datos que necesites. Parece un consejo estúpido, pero cuando se hacen las cosas por inercia se tienden a cometer fallos de bulto. Por ejemplo, si estoy manipulando una matriz 8000 x 8000, pero mi información realmente está en unos cuantos puntos de la misma, en vez de grabar la matriz entera, copia sólo los elementos que se necesiten y su posición en la matriz (si es relevante).
Eso te evitará generar ficheros extremadamente pesados cuando hablamos de matrices de datos grandes.
7 - Evita la salida de datos por pantalla. Si el ordenador tiene que escribir por pantalla la iteración en la que esta a cada instante y datos de la simulación, es tiempo que esta empleando en ofrecer información que quizá sea irrelevante una vez el programa ya funciona. Reduce la salida por pantalla hasta dejar sólo lo imprescindible.
8 - El uso de variables dinámicas (junto con ALLOCATE y DEALLOCATE) permite gestionar mucho más eficientemente la memoria disponible por el programa, la cual es finita. Esto se hace evidente cuando tienes que gestionar un volumen de datos muy grande. Utiliza este tipo de variables para aquellas variables tipo dummy, las cuales las quieres sólo para recoger una información que luego se guarda en otro sitio, quedando dicha variable inútil.
9 - Procura compactar tu código al máximo. Un código de 500 líneas es infinitamente más sencillo de repasar, corregir y mejorar que uno de 2000 líneas.
10 - Existen un sin fin de tutoriales y webs donde se describen todas las funciones intrínsecas y consejos sobre programación. Es fácil caer en el error de que uno sabe hacer bien las cosas, pero por norma general siempre se pueden hacer mejor. Sólo es cuestión de informarse adecuadamente.
Espero que estos consejos os ayuden en vuestros proyectos de programación.
Un saludo
domingo, 7 de abril de 2013
jueves, 25 de octubre de 2012
Ecuaciones hiperbólicas
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.
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.
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
Suscribirse a:
Comentarios (Atom)