Hola a todos
En esta ocasión voy a hablar muy escuetamente de lo que va a ser el futuro de la programación científica: la computación en paralelo.
Actualmente los procesadores tienden a tener un número de núcleos mayor de uno, tal y como era hace unos años: todos estamos familiarizados con los "Pentium core 2 duo" o los modernos "Pentium I3, I5, etc". Sin embargo, por lo que a mi respecta, la formación universitaria en ingeniería sólo cubre nociones básicas de programación en serie. Es obvio que, desde un punto de vista computacional, si se poseen varios núcleos de procesamiento, el realizar una programación capaz de utilizar simultáneamente la potencia de cálculo de esos cores va a ser mucho más ventajoso.
Sin embargo, hasta hace bien poco desconocía por completo un punto de partida que me pudiera orientar para empezar a caminar en este vasto campo. Pues bien, si queréis una piedra angular de referencia, vuestra palabra mágica es MPI.
El MPI se puede considerar como un anexo a los lenguajes de programación típicos en ciencias (Fortran, C) que permite procesar un código arrancando varios "threads", o hilos de computación. Realmente lo ideal es ejecutar tantos threads como núcleos tengas en tu cluster, pero no tiene por que ser así.
Para realizar un código en paralelo, hay que pasar de una mentalidad de cálculo en serie, y comenzar a pensar como se puede distribuir todas las tareas que tiene que hacer el código: habrá tareas, como una búsqueda de elementos en un fichero o matriz, que se puede "descomponer" en partes que pueden llevar a cabo cada núcleo. En cambio habrá otras que no se puedan descomponer de modo directo, como la escritura y lectura simultánea de un fichero.
Todas estas tareas deben ser coordinadas, esto es, normalmente se programa de forma que uno de los núcleos actúa como maestro y el resto como esclavos: el primero realizará tareas no paralelizables y se encargará de enviar y recibir información que los esclavos necesiten y/o proporcionen. Los cores esclavos se encargarán de realizar la computación respectiva, enviando o recibiendo información vital para el código.
MPI es una alternativa básica a la programación tradicional. Sin embargo, otras aplicaciones están desarrollándose rápidamente, con el afán de simplificar los procesos de comunicación entre cores y la distribución del trabajo: un ejemplo clásico es el OpenMP.
En próximos posts hablaré más en detalle de MPI y OpenMP.
Un saludo
lunes, 6 de febrero de 2012
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, 4 de diciembre de 2011
Cambio en server de LATEX
Hola
Debido a que ha habido cambios en el server de LATEX del cual suministré su enlace, he buscado otro server para poder mostrar LATEX en tu propia web sin cambios. El nuevo script que tenéis que introducir (sustituyendo al que ya tenías en caso de que siguierais las instrucciones tal y como expliqué en el anterior post) es el siguiente:
<script type="text/javascript" src="http://cdn.mathjax.org/mathjax/latest/MathJax.js">
MathJax.Hub.Config({
extensions: ["tex2jax.js","TeX/AMSmath.js","TeX/AMSsymbols.js"],
jax: ["input/TeX", "output/HTML-CSS"],
tex2jax: {
inlineMath: [ ['$','$'], ["\\(","\\)"] ],
displayMath: [ ['$$','$$'], ["\\[","\\]"] ],
},
"HTML-CSS": { availableFonts: ["TeX"] }
});
</script>
Con esto ya podréis escribir con LATEX.
Nota:
El script esta diseñado para que si se pone un dolar a cada lado de la escritura en latex ($ texto $), se mostrará el texto en línea, esto es, sin espacios. Si se escriben dobles dólares ($$ texto $$), se hace un punto y aparte y se separa del resto del cuerpo. Ejemplos:
Esto es una prueba de este texto $f=\sqrt{x^2+y^2}$ dentro de linea.
Esto es una prueba
$$ g=\frac{\partial}{\partial x}\left( z^2+6x \right) $$
de texto separado de cuerpo
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.
Suscribirse a:
Entradas (Atom)