You are currently browsing the monthly archive for julio 2017.

Siguiendo con las pruebas numéricas, suponemos el siguiente sistema de ODEs:

\begin{cases} \frac{d}{dm}u = u^2 v  \\  \frac{d}{dm}v = u v^2 \end{cases}

que es un sistema sencillo  con fuentes no lineales acopladas entre si del que conocemos su solución analítica:

u(m) = e^m, \,\,\,\, v(m) = e^{-m}

Tratamos de resolverlo tal y como hicimos en el post  Derivando sistemas de ODEs: derivando en ambos lados del sistema:

\begin{cases} \frac{d^2}{dm^2}u = \frac{d}{dm}(u^2 v)  \\  \frac{d^2}{dm^2}v = \frac{d}{dm}(u v^2) \end{cases}

De esta manera, utilizamos una discretización en diferencias finitas de las derivadas segundas de los lados izquierdos de las ecuaciones del sistema; mientras que derivamos numéricamente los términos fuente del lado derecho del sistema. Además, por ser no lineal, tenemos un bucle externo adicional en el que vamos actualizando los términos fuente cada vez que encontramos una solución para las fuentes que teníamos en el paso anterior. Las fronteras son las mismas que teniamos antes de derivar:

u(0)=1,\,\, u(1)=e,\,\, v(0)=1,\,\, v(1)=1/e

La siguiente secuencia de imágenes muestra como nos acercamos correctamente a la solución tras muy pocas iteraciones (tres) de la parte debida a la no linealidad:

Anuncios

Tal y como acabamos de ver en el post anterior, podemos derivar en ambos lados de una ODE para incrementar el orden de las derivadas. Además, hemos realizado la derivación del término fuente de manera numérica.

¿Qué pasa si en lugar de una única ODE tenemos un sistema de ODEs? Pues que todo sigue funcionando como se puede apreciar en la siguiente figura (el sistema de ecuaciones que se está resolviendo es trivialmente inferible):

system

Finalmente, si el sistema no es lineal, añadimos una iteración exterior que, dado que estamos asumiendo suavidad, convergerá en un número finito de iteraciones.

Otra opción que tenemos, si la fuente de la ecuación es lo suficientemente suave, es aplicar lo que comentamos en el post Derivando derivadas.

Para el mismo problema resuelto en Probando esquemas tipo Euler, lo que nos queda es \frac{\partial^2}{\partial m^2}f = \frac{\partial}{\partial m} \cos(m) con las mismas condiciones de frontera f(0)=0 y f(\pi) = 0.

Ahora tenemos un problema donde es aplicable un método iterativo, puesto que tenemos derivadas segundas. Realizamos la derivada del termino fuente de manera numérica. Las soluciones que obtenemos son (para 7 y 32 puntos respectivamente):

Volvemos a recuperar la solución f(m) = \sin (m) y con muchos menos puntos pues tenemos orden cuadrático.

Tal y como comenté hace pocos posts, las discretizaciones de ecuaciones con primeras derivadas llevan a esquemas tipo Euler (que, por tanto, podrían sustituirse por otros de mas alto orden, tipo Runge-Kutta por ejemplo).

Si resolvemos numéricamente el problema \frac{d}{dm}f = \cos(m) en el intervalo [0,\pi] con f(0) = 0 y f(\pi)=0, obtenemos las siguientes soluciones en función de si utilizamos esquemas hacia adelante, hacia atrás o mixto empezando por ambos lados (solución analítica en azul y numérica en amarillo):

El último es el único que asegura los valores conocidos de la solución en la frontera, dejando el error para el centro del dominio. Además, éste se reduce con el número de puntos de la malla utilizada (10, 50 y 1000 puntos):

Así pues, estábamos en lo cierto y, efectivamente, todo funciona como comentamos en el post Esquemas tipo Euler para ODEs en problemas de frontera

Partiendo del mismo problema tratado en otro post

\frac{d}{dm}f = g

con f(0) = a y f(M)=b, suponiendo que 0 \leq m \leq M. Suponiendo que las funciones sean lo suficientemente suaves, entonces podemos derivar a ambos lados de la expresion

\frac{d^2}{dm^2} f = \frac{d}{dm}g

Podriamos derivar analíticamente \frac{d}{dm}g. Sin embargo, si discretizamos, y utilizamos una derivación numérica del lado derecho de la expresion del mismo orden que la discretización, nos queda

\frac{f_{i-1}-2f_i+f_{i+1}}{h^2} = \frac{g_{i+1}-g_{i-1}}{2h}

Los valores de la frontera son condiciones sobre la solución, por lo que podemos utilizar las mismas. Solo hemos cambiado la manera de obtener ésta. ¿Es todo esto correcto?

 

Tratando de discretizar el problema

\frac{d}{dm}f = k_1

sobre una malla de cuatro puntos con f(0) = k1 y f(3) = 3k_1+k_2

nos encontramos con el esquema

f_i = h k_1 + f_{i-1}

si utilizamos la derivada hacia atrás, y con el esquema

f_i = f_{i+1} - h k1

si utilizamos la derivada hacia adelante.

Observamos que la información fluye de una frontera hacia la otra o viceversa. Si empezamos con el valor de la frontera f(0) nos llega f(3)=4k_1 \neq 3k_1+k_2 y si empezamos con la frontera f(3) = 3k_1+k_2 con el otro esquema, nos llega f(0)=k2 \neq k1.

En realidad, estamos obteniendo un esquema de Euler para ODEs. Son problemas de valor inicial. Allí solo sabemos el valor en t_0, en una de las fronteras, y el mayor error se produce en el tiempo final t_f, la otra frontera. Como sabemos el orden del método, sabemos el orden de este error producido. Ahora, como es un problema de frontera, sabemos el valor en ambos lados. Podríamos empezar dos esquemas de este tipo, uno desde cada frontera, y parar a la mitad, donde habría una discontinuidad, pero del orden del error del método de Euler. Sin embargo, ¿como se realizan estos recorridos en mas de una dimensión? (¿O mes mas sencillo propagar todo y realizar la media?)

Supongamos que tenemos la ecuación

\frac{d}{dx} = f(u,x)

y que queremos hacer un estudio de su estabilidad de von Neumann. Para empezar, necesitamos la discretización con factor de relajación y, además, asumimos que el termino fuente no afecta al mismo. Lo que nos queda es

u_{i}^{k+1} = (1-\omega) u_i^k +  \omega u_{i-1}^k

de donde tenemos que

\varepsilon_{i}^{k+1} = \varepsilon_i^k - \omega ( \varepsilon_i^k - \varepsilon_{i-1}^k)

Como lo que nos interesa es el factor de amplificación G = \frac{\varepsilon_i^{k+1}}{\varepsilon_i^k}, lo que nos queda es

G = 1 - \omega (1-e^{-i k h})

Siendo k los números de onda. Así como en las discretizaciones para el Laplaciano derivaba en expresiones donde podíamos combinar las exponenciales para formar funciones trigonométricas, aquí no pasa eso.

Si añadimos discretizaciones de mas puntos, no centradas para obtener el punto actual, seguimos teniendo el mismo problema. Por ejemplo, para el caso de no centrada de orden dos, lo que nos queda es

G = 1- \omega (-2 + 2e^{-i k h} + \frac{1}{2}- \frac{1}{2} e^{-i k 2h})

 

En el post anterior hemos discretizado un sistema de ecuaciones donde únicamente aparecen primeras derivadas.

Para que quede claro lo que queremos decir, vamos a trabajar con la ecuación genérica

\frac{d}{dx}u = f(u,x)

Ya en el punto de la discretización necesitamos comentar dos cosas. Por una parte, dado que queremos escribir el método iterativamente, necesitamos que en la discretización aparezca un término con el punto actual. En caso contrario, no podemos despejar la variable en el punto actual. De esta manera, para la discretización de la primera derivada a orden uno hacia atrás, tenemos:

\frac{u_{i} - u_{i-1}}{h} = f(u,x) \Leftrightarrow u_i = h f(u,x) + u_{i-1}

Y en diferencias hacia adelante, nos encontramos con:

\frac{u_{i+1} - u_i}{h} = f(u,x) \Leftrightarrow u_i = u_{i+1} - h f(u,x)

Sin embargo, si utilizamos la aproximación a orden dos centrada, tenemos

\frac{u_{i+1} - u_{i-1}}{2h} = f(u,x)

nos encontramos con que no podemos despejar u_i para implementar la iteración ya que este termino no aparece en la ecuación.

Pero no solo eso. Además, la discretización que utilicemos debería dar lugar a una matriz con todos los autovalores con el mismo signo, para que una iteración tipo Jacobi pueda aplicarse. La primera discretización y la segunda cumplen esta propiedad (dan lugar a una matriz con 1 en la diagonal y -1 arriba o abajo de la diagonal en función de la discretización escogida). La última, sin embargo, da lugar a una matriz con todos los autovalores complejos (0 en la diagonal, 1 y -1 arriba y abajo de esta).

 

Asumiendo una estrella con simetría esférica y en complete equilibrium (hydrostatic equilibrium + thermally adjusted), las ecuaciones diferenciales básicas que nos darán su estructura son:

\frac{\partial}{\partial m}r = \frac{1}{4 \pi r^2 \rho}

\frac{\partial}{\partial m}P = - \frac{G m}{4 \pi r^4}

\frac{\partial}{\partial m}l = \varepsilon_n + \varepsilon_\nu

\frac{\partial}{\partial m}T = - \frac{G m T}{4 \pi r^4 P} \nabla

donde \nabla = \frac{d \ln T}{d \ln P}. Queremos encontrar funciones solución para la distancia radial r(m), para la presión P(m), para la energía l(m) y para la temperatura T(m) donde 0 \leq m \leq M, siendo M la masa total de la estrella. Necesitamos conocer los valores en la frontera de las funciones incógnita, esto es, r(0), P(0), l(0), T(0),  r(M), P(M), l(M), T(M). A través de la relación de operadores

\frac{\partial}{\partial m} = \frac{1}{4 \pi r^2 \rho} \frac{\partial}{\partial r}

deberíamos de recuperar las ecuaciones de la hidrodinámica en steady-state expresadas en su forma habitual.

Podemos discretizarlas a primer orden:

\frac{r_{i} - r_{i-1}}{\Delta m} = \frac{1}{4 \pi r_{i/2}^2 \rho_i}

\frac{P_{i} - P_{i-1}}{\Delta m} = - \frac{G m_i}{4 \pi r_{i/2}^4}

\frac{l_{i} - l_{i-1}}{\Delta m} = (\varepsilon_n)_i + (\varepsilon_\nu)_i

\frac{T_{i} - T_{i-1}}{\Delta m} = - \frac{G m_{i} T_{i/2}}{4 \pi r_{i/2}^4 P_{i/2}} \nabla_{i/2}

Si despejamos las variables en el punto i de la discretización, obtenemos:

r_{i}^{*} = \Delta m \frac{1}{4 \pi r_{i/2}^2 \rho_i} + r_{i-1}

P_{i}^{*} = - \Delta m \frac{ G m_i}{4 \pi r_{i/2}^4} + P_{i-1}

l_{i}^{*} = \Delta m [ (\varepsilon_n)_i + (\varepsilon_\nu)_i] + l_{i-1}

T_{i}^{*} = - \Delta m \frac{G m_{i} T_{i/2}}{4 \pi r_{i/2}^4 P_{i/2}} \nabla_{i/2} + T_{i-1}

y con esto, podemos obtener esquemas con factor de relajación:

r_{i}^{k} = (1 - \omega) r_{i}^{k-1} + \omega r_{i}^{*}

P_{i}^{k} = (1 - \omega) P_{i}^{k-1} + \omega P_{i}^{*}

l_{i}^{k} = (1 - \omega) l_{i}^{k-1} + \omega l_{i}^{*}

T_{i}^{k} = (1 - \omega) T_{i}^{k-1} + \omega T_{i}^{*}

 

Este es el título de mi tesis para optar al grado de Doctor en Física. La defensa será el próximo día 19 de Julio a las 12h en el Salón de Grados “Manuel Valdivia” de la Facultat de Ciències Matemàtiques de la Universitat de València. Estáis todos invitados 😉

julio 2017
L M X J V S D
« Feb   Ago »
 12
3456789
10111213141516
17181920212223
24252627282930
31