Fundamentos matemáticos de la mecánica de los medios continuos
5º matemáticas - Curso 2007/08
Movimiento de partículas sometidas a la acción de un potencial

ALUMNO:
	Apellidos y nombre: 
	D.N.I.:

Esta práctica es voluntaria y se considerará como evaluación complementaria. 
Se ha elegido el programa Mathematica para su implementación, aunque cada alumno puede optar por usar este u otro programa que considere oportuno.

La introducción que sigue es complementaria de la proporcionada durante las clases; lea detenidamente su contenido antes de pasar a la resolución numérica del problema.

El plazo (voluntario) de entrega finaliza el 8 de noviembre y se hará en formato electrónico.

INTRODUCCIÓN: Curvas características y ecuación de Liouville (Vlasov)

	En esta práctica pretendemos resolver la ecuación de Liouville, que modela el movimiento de un sistema de partículas sometidas a la acción de un campo de fuerzas. Recordemos brevemente algunas leyes físicas elementales que nos ayudarán, no sólo a deducir la ecuación, sino a hallar un método de resolución muy efectivo. 
	Consideramos una partícula (de masa 1 para no arrastrar constantes) que se mueve en el plano (ℝ^2) (también podemos trabajar en (ℝ^3) pero el modelo en (ℝ^2) basta para comprender el modelo). Llamamos t al tiempo transcurrido desde el momento inicial (que normalmente se toma como 0). Entonces, si en el instante inicial t = 0 nuestra partícula ocupaba una posición inicial (Subscript[p, 0]) = ((Subscript[x, 0]), (Subscript[y, 0])) y tenía una velocodad inicial (Subscript[v, 0]) = ((Subscript[v1, 0]), (Subscript[v2, 0])) conocidas, queremos hallar su trayectoria (al hacer actuar sobre ella un campo de fuerzas), es decir, su posición en cada instante de tiempo posterior, que notaremos como p(t) = (X(t), Y(t)).
	Antes de continuar con la partícula veamos qué es un campo de fuerzas; no es más que una aplicación vectorial:
	        \(TraditionalForm\`F : ℝ\^2 -> \ ℝ\^2\),  F(x, y) = (f1(x, y), f2(x, y)).

F(x,y) se interpreta como la fuerza (por unidad de masa) que actúa sobre una partícula situada en el punto (x,y). Por lo tanto, cuando nuestra partícula esté en un punto p(t), sobre ella actuará una fuerza F(p(t)). Bien, veamos cómo usar este hecho para determinar p(t); usaremos la Segunda Ley de Newton. 
	Ya que la velocidad es la derivada de la posición con respecto al tiempo y la aceleración a su vez es la derivada de la velocidad con respecto al tiempo, podemos escribir:
                               p'(t) = v(t)     y     v'(t) = a(t).
Como decíamos, la Segunda Ley de Newton (fuerza = masa x aceleración), nos permite cerrar estas dos ecuaciones como sigue 
                               p'(t) = v(t)     y     v'(t) = a(t) = F(p(t)).
Si escribimos estas EDO's (ecuaciones diferenciales ordinarias) en cada componente nos queda:
    	X'(t) = v1(t),
    	Y'(t) = v2(t),
    	v1'(t) = f1(X(t),Y(t)),              (ECUACIONES CARACTERÍSTICAS)
    	v2'(t) = f2(X(t),Y(t)),
con condiciones iniciales en t = 0:
    	X(0) = (Subscript[x, 0]),
    	Y(0) = (Subscript[y, 0]),
    	v1(0) = (Subscript[v1, 0]),
    	v2(0) = (Subscript[v2, 0]).
Por lo tanto si resolvemos este sistema de 4 EDO's de primer orden con estas 4 condiciones iniciales, hallaremos la trayectoria seguida por nuestra partícula, que llamaremos CURVA CARACTERÍSTICA asociada al campo F que pasa por (Subscript[p, 0]) = ((Subscript[x, 0]), (Subscript[y, 0])) con velocidad (Subscript[v, 0]) = ((Subscript[v1, 0]), (Subscript[v2, 0])) en el instante t = 0. Esto será lo que hagamos en el ejercicio que hay más abajo.

Rápida deducción de la ecuación de Liouville
	En la práctica se tiene una cantidad muy elevada de partículas (del orden del número de Avogadro que es aproximadamente 10^23) por lo que el número de ecuaciones es enorme y además las posiciones y velocidades iniciales son muy difíciles de obtener, por no decir imposible. Es entonces cuando la mecánica estadística entra en juego. En lugar de imaginarnos las posiciones de cada una de las partículas, nos imaginamos una función de probabilidad f(t,p,v) que representa la probabilidad de encontrar una partícula en el instante t en la posición p y con velocidad v. Por coherencia con la representación de las curvas caracteríaticas se supone que "la probabilidad de que en t = 0 encontremos una partícula con posición (Subscript[p, 0]) y velocidad (Subscript[v, 0]) es la misma que la de encontrar una partícula en en instante t con posición p(t) y velocidad v(t) siendo éstas últimas las curvas características asociadas".  Matemáticamente escribimos:
	
	f(0, (Subscript[x, 0]), (Subscript[v, 0])) = f(t, p(t), v(t)).
	
De esta menera es equivalente conocer f a conocer las curvas características. Finalmente, si derivamos en la igualdad anterior con respecto a la variable t obtendremos:

0 = ∂ (f/∂) t + (Subscript[∇, p]) f  . 
              p' (t) + (Subscript[∇, v]) f  . v' (t).

Si usamos que p' = v   y que v' = F(p) obtenemos finalmente la ecuación de Liouville:

     ∂ (f/∂) t + 
          v . (Subscript[∇, p]) f  + F(p) . (Subscript[∇, v]) f = 
        0,      con  f(0, x, v)  conocida. 
     
que es, quizás, la ecuación de evolución (se conose lo que pasa en un instante y se quiere ver cómo evoluciona en el tiempo) más sencilla de la mecánica clásica.

En el primer ejercicio que nos proponemos consideramos sólamente 20 partículas sometidas a la acción del siguiente campo de fuerzas:
F(x,y) = (f1(x,y), f2(x,y))  donde:

   		f1(x, 
          y) = -((x (-54 + 20 (x^2) + 
                    (x^4)))/(((10 + (x^2)))^2)), 
f2(
          x, y) = ((1 - 2 y )/2).
NOTA: Este campo es conservativo, es decir, es el gradiente de alguna función V llamada potencial. En este caso 
	\(TraditionalForm\`F = \(\(-∇V\)\ \ \ donde\ \(V(x, 
              y)\) = \(\((x + 2)\) \((x + 1)\) \((x - 1)\) \((x - 2)\)\)\/\(20 + 2  x\^2\) + \(\(\(y\ \((y - 1)\)\)\/2\)\(.\)\)\)\)

Las posiciones iniciales de cada una de las 20 partículas quedan descritas más abajo, en la lista puntos y consideraremos que inicialmente están en reposo, es decir, los (Subscript[v1, 0]) y (Subscript[v2, 0]) de todas ellas son cero.
Pretendemos resolver numéricamente (de manera aproximada) las ecuaciones características:
    	X'(t) = v1(t),
    	Y'(t) = v2(t),
    	v1'(t) = f1(X(t),Y(t)), 
    	v2'(t) = f2(X(t),Y(t)),
para cada una de las partículas, es decir, poniendo las 20 condiciones iniciales asociadas a cada punto:
    	{(Subscript[x, 0]), (Subscript[y, 0])} = puntos[[i]],  para cada i=1,2,...,20,
    	{(Subscript[v1, 0]),(Subscript[v2, 0])} = {0, 0}   (inicialmente suponemos que están en reposo)
El objetivo principal de este ejercicio es dibujar las sucesivas posiciones de las 20 partículas en los instantes t = 0,(1/2),1,(3/2),2,...,30 para animarlas luego a modo de película. Por ello fijaremos una ventana con x entre -2.5 y 2.5 e y entre -0.2 y 1.5 con la opción PlotRange.
Veamos cómo resolver aproximadamente estas 4 ecuaciones.

A. Discretización en tiempo
Primero hacemos una partición del intervalo de tiempo [0,30] en 3000 trozos iguales, es decir, tomamos h = (1/100)como paso temporal y llamamos (Subscript[t, i]) = i h,  i = 0, 1,  ... , 3000 a los sucesivos instantes de tiempo con distancia h entre ellos. Así los instantes de tiempo en los que queremos dibujar 

	t = 0,(1/2),1,(3/2),2,...,30    corresponden con i = 0,50,100,150,...,3000.
	
B. Paso temporal general ((Subscript[t, i - 1]) ->  (Subscript[t, i])). (Método de Splitting)
Lo que vamos a hacer es:
-conocidos p(0) y v(0), aproximamos p((Subscript[t, 1])) y v((Subscript[t, 1])),
-conocidos p((Subscript[t, 1])) y v((Subscript[t, 1])), aproximamos p((Subscript[t, 2])) y v((Subscript[t, 2])),
...
-conocidos p((Subscript[t, i - 1])) y v((Subscript[t, i - 1])), aproximamos p((Subscript[t, i])) y v((Subscript[t, i])),
...
-conocidos p((Subscript[t, 2999])) y v((Subscript[t, 2999])), aproximamos p((Subscript[t, 3000])) y v((Subscript[t, 3000])).
es decir, dar los 3000 pasos uno a uno. Veamos cómo dar un paso genérico ((Subscript[t, i - 1]) ->  (Subscript[t, i])):
     Datos conocidos en t = (Subscript[t, i - 1]):       ->   	Datos buscados  en t = (Subscript[t, i - 1]) + h= (Subscript[t, i])
     posición:  	(Subscript[x, i - 1]) e (Subscript[y, i - 1]),			posición:  	(Subscript[x, i]) e (Subscript[y, i]),
     velocidad:  	(Subscript[v1, i - 1]) y (Subscript[v2, i - 1]),     		velocidad:  	(Subscript[v1, i]) y (Subscript[v2, i]).
     
B1. Medio paso tomando campo constante
Resolvemos 	v1'(t) = f1(X(t),Y(t)),      v1((Subscript[t, i - 1])) = (Subscript[v1, i - 1])
		v2'(t) = f2(X(t),Y(t)),      v2((Subscript[t, i - 1])) = (Subscript[v2, i - 1]) 
como si x e y no variasen en un tiempo h/2 (la mitad del paso), es decir, son constantemente (Subscript[x, i - 1]) e (Subscript[y, i - 1]), lo que hace que el campo, f1 y f2, no varíe. Por lo tanto la solución en (Subscript[t, i - 1]) + (h/2) es:
 	(Subscript[v1, i]) = v1((Subscript[t, i - 1]) + (h/2)) = 
          (Subscript[v1, i - 1]) + (h/2) f1((Subscript[x, i - 1]), (Subscript[y, i - 1]))
	(Subscript[v2, i]) = v2((Subscript[t, i - 1]) + (h/2)) = 
          (Subscript[v2, i - 1]) + (h/2) f2((Subscript[x, i - 1]), (Subscript[y, i - 1]))
 
B2. Medio paso tomando velocidad constante 
Ahora resolvemos 	X'(t) = v1(t),      X((Subscript[t, i - 1])) = (Subscript[x, i - 1])
			Y'(t) = v2(t),      Y((Subscript[t, i - 1])) = (Subscript[y, i - 1]) 
como si v1 y v2 no variasen en un tiempo h/2 (la otra mitad del paso) y se mantuviesen constantemente los obtenidos en B1: \(TraditionalForm\`\(v1\^*\)\) y \(TraditionalForm\`\(v2\^*\)\), lo que hace que la solución en (Subscript[t, i - 1]) + (h/2) + (h/2) = (Subscript[t, i]) sea:
 	(Subscript[x, i]) = X((Subscript[t, i])) = 
          (Subscript[x, i - 1]) + (h/2) (Subscript[v1, i]),
	 (Subscript[y, i]) = Y((Subscript[t, i])) = 
          (Subscript[y, i - 1]) + (h/2) (Subscript[v2, i]).

C. EJERCICIO.
C1. Programa el paso temporal general descrito en B como una función llamada paso[{x_,y_,v1_,v2_}]:=..., cuyas entradas {x,y,v1,v2} sean las correspondientes a posición y velocidad en un instante y cuya salida sea también una lista con 4 elementos correspondientes a posiciones y velocidades en el instante aumentado en h.
C2. Repite los 3000 pasos, pero dibujando (en 2D) la distribución de puntos cada 50 pasos (es decir cada 1/2 unidad de tiempo) en el mismo rango [-2.5,2.5] para x, y [-0.2,1.5] para y, y así poder hacer una película con la secuencia de las sucesivas posiciones.
C3. Repite el paso anterior pero dibujando en una gráfica 3D simultáneamente la gráfica del potencial V(x,y) y los puntos sobre la misma (x,y,V(x,y)). 
C4. Repite el paso anterior pero elige tú varios potenciales V que desees y la distribución inicial de patículas (su cantidad, sus posiciones y sus velocidades).

(*lista inicial de puntos*)
puntos = 
          Table[{(-1)^i, (2  i - 1 - (-1)^i)/40}, {i, 1, 
              20}];

      ListPlot[puntos, PlotStyle -> PointSize[0.02], 
          PlotRange -> {{-2.5, 2.5}, {-0.2, 
                1.5}}];

Solución:
Desarrolla las soluciones a partir de aquí

(*construcción de F a partir de V*)
V[x_, 
            y_] := (x + 2) (x + 1) (x - 
                  1) (x - 2)/(20 + 2  x^2) + 
            y (y - 1)/2;

      Clear[x, y]

      f1[x_, y_] = -D[V[x, y], x] // Simplify;

      f2[x_, y_] = -D[V[x, y], y] // Simplify;

      Print["f1(x,y)=", f1[x, y]]

      Print["f2(x,y)=", f2[x, y]]

(*Gráficas del potencial*)
\(\(Plot[V[x, 0], {x, \(-2.3\), 2.3}, 
          PlotLabel -> "perfil del potencial para y=0"];\)

      \(Plot[V[0, y], {y, 0, 1}, 
          PlotLabel -> "perfil del potencial para x=0"];\)

      \(Plot3D[V[x, y], {x, \(-2.3\), 2.3}, {y, 0, 1}, Mesh -> False, 
          PlotLabel -> "potencial en 3D"];\)\)

(*Ejercicio C1*)
h = 1/100. ; 

      paso[{x_, y_, v1_, 
            v2_}] := {x + (h/2)*(v1 + (h/2)*f1[x, y]), 
          y + (h/2)*(v2 + (h/2)*f2[x, y]), 
          v1 + (h/2)*f1[x, y], v2 + (h/2)*f2[x, y]}

(*Ejercicio C2*)
n = 3000;

      borde = Table[{i, puntos}, {i, n}];

      For[p = 1, 
        p ≤ 
          20, p++, 
{v1, v2} = {0, 
            0}; 
{x, y} = 
          borde[[1, 2, p]]; 
For[i = 1, 
          i ≤ n - 1, i++, 
{x, y, v1, v2} = 
            paso[{x, y, v1, v2}]; 
borde[[i, 2, 
                p]] = {x, 
              y};
]
]


      For[i = 1, i ≤ 3000, 
        i = i + 50, 
ListPlot[borde[[i, 2]], 
          PlotStyle -> PointSize[0.02], 
          PlotRange -> {{-2, 2}, {-0.5, 1.5}}]]