Las paradojas del C.S.I.

En esta nueva entrada de la serie de las Paradojas (ver las otras entradas aquí ,
aquí y aquí) vamos a tratar de explicar paradojas que nacen de la mala interpretación de la probabilidad.

Todo el mundo conoce la serie C.S.I. y sus secuelas. En la serie un grupo de investigadores de la policía se dedican a buscar pruebas para resolver casos que sin el uso de avanzadas técnicas científicas y la participación de estos peritos forenses sería imposible solucionar. Los capítulos de la serie cuentan el asesinato y cómo los miembros del C.S.I. han recabado las pistas y acaban con la detención del culpable. Sin embargo hay una última parte que nunca tratan. Me refiero al juicio. Y es que, aunque tal cual está en la serie parece obvia la culpabilidad dadas las pruebas que se han encontrado, la valoración de las evidencias es un tema peliagudo en extremo.

El problema si lo escribimos de manera matemática involucra probabilidades condicionadas. La definición de este objeto (sacada de la wikipedia) es, dado un espacio de probabilidad (\Omega, \mathcal F, P) y dos eventos o sucesos A, B\in \mathcal F con P(B)>0, la probabilidad condicional de ‘A’ dado ‘B’ está definida como:

P(A|B) = \frac{P(A \cap B)}{P(B)}.

Hemos de pensar en B como en que ha ocurrido algo que modifica la información que tenemos disponible.

Así, si dado un caso de los de la serie de televisión, consideramos los sucesos

1) “se encuentran evidencias”. Este suceso se escribe E

2) “el sospechoso es culpable”. Este suceso se escribe C.

3) “el sospechoso es inocente”. Este suceso se escribe I.

4) “el sospechoso es inocente dado que hemos encontrado la prueba E”. Este suceso lo escribimos como P(I|E).

y

5) “el sospechoso es culpable dado que hemos encontrado la prueba E”. Este suceso lo denotamos por P(C|E).

Así la probabilidad que le interesa al juez con sobrecarga de trabajo que le toque es P(I|E) o bien P(C|E).

Vamos a motivar nuestra nueva entrada sobre paradojas con un caso famosísimo en su época: el caso Dreyfus (ver esta página para más detalles). Esta novela por entregas tiene de todo, un bueno inocentón, el propio Dreyfus, un malo malísimo con multitud de problemas, Esterházy, y hasta cameos de personajes famosos como por ejemplo Zola. Un brevísimo resumen: tras la guerra franco-prusiana del siglo XIX se detectó una filtración en el ejército francés. Se detuvo al Capitán Alfred Dreyfus. Se tenía como principal evidencia una carta. Lo que ocurrió es que un “perito” estimó la probabilidad de coincidencia entre la letra de la carta y la letra del Capitán Dreyfus, suponiendo que Dreyfus fuese inocente, como 0.2. Como este mismo perito había observado 4 coincidencias y éstas se asumían independientes (es decir, que la probabilidad de las cuatro es el producto de las cuatro probabilidades) el perito conluyó que Dreyfus era culpable porque, según él, “la probabilidad de que fuese inocente dada la evidencia (las coincidencias entre la carta y la letra de Dreyfus) es muy pequeña”. El lector avispado se habrá dado cuenta del truco. El perito “estima” (que esa es otra cosa, ¿por qué la probabilidad es 0.2?) P(E|I)=(0.2)^4 y entonces concluye que P(I|E)=P(E|I)=(0.2)^4.

Veamos qué implicaciones tiene la asunción P(I|E)=P(E|I). Se tiene que
P(I|E)=\frac{P(I \cap E)}{P(E)}=\frac{P(E \cap I)}{P(I)}\frac{P(I)}{P(E)}=P(E|I)\frac{P(I)}{P(E)},
y, similarmente,
P(C|E)=\frac{P(C \cap E)}{P(E)}=\frac{P(E \cap C)}{P(C)}\frac{P(C)}{P(E)}=P(E|C)\frac{P(C)}{P(E)}.
Ahora consideramos la cantidad
\frac{P(I|E)}{P(C|E)}=\frac{P(E|I)\frac{P(I)}{P(E)}}{P(E|C)\frac{P(C)}{P(E)}}.
Si ahora usamos la hipotesis de que P(I|E)=P(E|I) y P(C|E)=P(E|C) se concluye que
P(I)=P(C), lo que contradice la presunción de inocencia en casi todos los casos.

¿Cuál es la moraleja? Pues que la probabilidad no es intuitiva nunca (o casi nunca) y que todos los argumentos probabilísticos deben ser cuidadosamente repasados. Si no lo hacemos así nos arriesgamos a repetir tristes historias (ver esta entrada de la  wikipedia o esta otra).

–Referencias: 1) Estadística y Evaluación de la Evidencia para Expertos Forenses. Segunda Edición. C. Aitken y F. Taroni. (Traducido por: J. Lucena, L. Gil y R. Granero), Ed. Dykinson.

— Con esta entrada participamos en la edición 2.7 del Carnaval de Matemáticas, que está organizado por La Aventura de la Ciencia.

Las paradojas de la probabilidad

En esta nueva entrada de la serie de las paradojas (las primeras estradas son http://scientiapotentiaest.ambages.es/?p=244 y http://scientiapotentiaest.ambages.es/?p=266) nos vamos a centrar en las paradojas que vienen de la probabilidad. Y es que la probabilidad pese a ser algo bastante mencionado en la vida diaria no es entendida por mucha gente.

Voy a comenzar tratando el problema de Monty Hall. Este problema es muy divertido, porque todo el mundo al que se lo contemos dirá que es muy intuitivo. ¡El problema es que para cada uno será intuitiva una respuesta distinta!

Supongamos que estamos en un programa de estos de la tele. En este programa hay tres puertas. Una de las tres tiene un premio y las otras dos no tienen nada. La puerta que tiene el premio se elige al azar de manera equiprobable. El juego consiste en elegir una puerta. Tras nuestra elección, el presentador, que sabe dónde está el premio, abre una de las dos puertas que no hemos elegido y que no tiene nada. Tras esto nos pregunta si queremos cambiar de puerta. La pregunta viene ahora, ¿qué nos sale mejor como jugadores? ¿cambiar nuestra primera elección o dejarla estar?

Observamos que el premio no cambia de puerta mientras nosotros jugamos.

La respuesta correcta es que siempre nos sale mejor cambiar de puerta. Sí, amigos lectores, es mejor cambiar. ¿Qué no me creéis? Vamos a ver si os convenzo. Como el premio está en una de las puertas de manera equiprobable la probabilidad de que el premio se encuentre en la puerta que hemos elegido al comenzar a jugar es 1/3. Ahora el presentador abre una de las otras dos puertas sin descubrir el premio. ¿Cuál es la probabilidad de que nuestra puerta esconda el premio? LA MISMA, 1/3. Por lo tanto si cambiamos de puerta nuestra nueva elección tendrá probabilidad 2/3. Vamos a ilustrar esto un poco más mirando los posibles casos. Vamos a suponer que elegimos siempre la puerta 1, como todo es simétrico todos los casos estarán reflejados igualmente.

1) El coche está en la puerta 1.
1.a) Cambiamos de puerta y perdemos.
1.b) No cambiamos y ganamos.
2) El coche está en la puerta 2.
2.a) Cambiamos de puerta y ganamos.
2.b) No cambiamos y perdemos.
3) El coche está en la puerta 3.
3.a) Cambiamos de puerta y ganamos.
3.b) No cambiamos y perdemos.

Hagamos el recuento: Si cambiamos ganamos 2 veces frente a una que perdemos. Si no cambiamos la situación es recíproca. Como hay 3 casos posibles la probabilidad de ganar cambiando es 2/3.

Otra manera de convencerse es aumentar el número de puertas. Supongamos que tenemos 100 puertas. Entonces la probabilidad de que el premio esté en nuestra puerta es 1/100. Tras elegir, el presentador abre todas las puertas menos dos, la que tu elegiste y otra, de manera que ninguna de las abiertas mostrase el premio. Ahora parece claro que lo mejor es cambiar.

Otra paradoja que invariablemente se estudia en los cursos iniciales de probabilidad es la paradoja de los hijos (a falta de un nombre mejor en castellano). Supongamos que la probabilidad de que nazca un chico o una chica es la misma, 50%, y que el sexo de cada hijo es independiente del sexo de los demás hermanos. Vamos con el problema. La señora López tiene dos vástagos. Supongamos, a modo de calentamiento, que el hijo mayor de la señora López es un niño, ¿cuál es la probabilidad de que el menor sea niña? La respuesta, claro, es 50%. Ahora se nos plantea un problema ligeramente distinto: Supongamos que uno de los hijos de la señora López es niño, ¿cuál es la probabilidad de que el otro hijo también sea niño? ¿50%? Veamos los casos tal cuál está escrito (observando que la familia está fija para los más puntillosos) y ordenando los niños por edad:

1) Chico-chico
2) Chica-chico
3) Chico-chica

Hay 3 casos posibles y sólo 1 es el que nos preguntan. Como son equiprobables la probabilidad de que ambos hijos de la señora López sean niños es 1/3. El truco aquí está en que al decir que uno de los hijos es chico perdemos información sobre su edad. En la primera pregunta tenemos más información y podemos descartar más casos.

La probabilidad se va descubriendo un poco farragosa y “extraña” en ocasiones. Sin embargo todavía no hemos tratado la principal fuente de problemas al tratar con la probabilidad. Me refiero a la ambigüedad. Las palabras “al azar” no tienen un significado preciso y se usan con más asiduidad de la recomendable. Esto quedó claramente expuesto en la obra de J. Bertrand “Calcul des probabilités”. En este manual propuso un enunciado y dio varias respuestas al problema, todas bien lógicas y correctas. El problema era que en cada uno la definición de “al azar” era distinta. Es ahora cuando nos ponemos técnicos. En la entrada anterior ya mencionamos cosas como “medidas” de conjuntos (ver la entrada anterior). Eso nos volverá a ser útil, pues la manera matemática de dar sentido a las palabras “al azar” utiliza esas ideas. En matemáticas un “espacio de probabilidad” es un espacio de sucesos posibles en nuestro experimento, una manera de agruparlos en conjuntos y una manera de “medir la probabilidad” cada uno de estos conjuntos de sucesos. Claro está que la definición precisa es mucho más técnica. Así la idea que quiero dejaros es que esas tres cosas abstractas es lo que da sentido a la palabra “azar”.

El ejemplo que más me gusta de esto es la paradoja de Bertrand. Consideremos un círculo con un triángulo equilatero inscrito. La pregunta ahora es ¿cuál es la probabilidad de que una cuerda trazada “al azar” sea más larga que los lados del triángulo? Veamos las tres maneras clásicas de calcular dicha probabilidad:

1) Supongamos, sin pérdida de generalidad, que uno de los extremos de la cuerda coincide con uno de los vértices del triángulo. En ese caso al quedar la circunferencia dividida en 3 trozos iguales y coincidir uno con el conjunto donde la cuerda es más larga que el lado del triángulo concluimos que la probabilidad pedida es 1/3. Aquí el espacio de sucesos es el conjunto de puntos de la circunferencia y la medida de probabilidad asociada es la longitud del arco considerado dividido la longitud total de la circunferencia.

2) Consideremos ahora un radio perpendicular a uno de los lados del triángulo. Esto de nuevo nos facilita la vida pero no perdemos generalidad. Ahora trazamos la cuerda de manera perpendicular a dicho radio por un punto aleatorio del mismo. La probabilidad de que la cuerda sea mas larga que el lado es justamente 1/2 si la trazamos con estas reglas. Aquí el espacio de sucesos es el conjunto de puntos del radio y la medida de probabilidad es la longitud del segmento considerado dividido por la longitud total.

3) Ahora nuestro experimento consiste en elegir el punto medio de la cuerda. Si trazamos la cuerda con esta regla la probabilidad de que la longitud de la cuerda sea mayor que la del lado del triángulo es la misma que la probabilidad de que el punto medio de la cuerda esté en un círculo concentrico inscrito en el triángulo de radio la mitad (y área 1/4 del área del círculo original). Por lo tanto la probabilidad será 1/4. Aquí el espacio de sucesos es el conjunto de puntos del círculo y la medida de probabilidad es el área del conjunto de puntos considerado dividido entre el área del círculo original.

La paradoja aquí estriba en que para una pregunta tenemos tres respuestas. Además, eligiendo la definición de “azar” todas son correctas.

Creo que ha quedado claro que cuando decimos cosas como “aleatorio” o “al azar” no estamos en realidad diciendo nada y que debemos entrar, aunque no queramos, en tecnicismos para evitar este tipo de paradojas.

— Nota: Las imágenes han sido obtenidas de la Wikipedia en inglés

— Nota: Con esta entrada participamos en el Carnaval de Matemáticas en su edición 2.6 (albergado por “La vaca esférica” , http://lavacaesferica.com/).

Cuántica e incertidumbre (Parte I)

Hoy vamos a hablar, procurando no ser muy rigurosos (aparecen algunas fórmulas pero se pueden saltar), de por qué la mecánica cuántica es necesaria y del principio de incertidumbre de Heisenberg. La entrada la hemos escrito de manera conjunta Rafa y yo mismo.
Werner Heisenberg

Werner Heisenberg, uno de los padres de la Mecánica Cuántica

Principios de incertidumbre hay muchos, básicamente no es mas que una desigualdad de cierto tipo, donde un término controla a otro. En el caso concreto que nos atañe afirma que es imposible conocer con precisión la posición y la velocidad de una partícula cuántica (para fijar ideas un electrón). Los físicos (experimentales al menos) razonan diciendo que eso es porque para detectar un electrón hay que golpearlo con un fotón y entonces cambias su velocidad.
Esto nos deja con la duda: ¿habrá algún método que permita conocer las dos cosas?
La respuesta es que no. Y el motivo, enraizado en las matemáticas, tiene que ver con la transformada de Fourier. El ‘teorema’ aquí es que si tienes una función “grande” su transformada de Fourier es “pequeña” y viceversa (se puede intuir usando la desigualdad de Haussdorf-Young por ejemplo.). ¿Cómo afecta esto a nuestro electrón?. Los impacientes estarán pensando que se me está yendo la olla.
Pero comencemos aclarando por qué todas estas cosas tan raras de la cuántica son necesarias: supongamos que tenemos un átomo de hidrógeno (un electrón y un protón), y que tanto nuestro electrón como el núcleo del átomo se comportan como partículas clásicas. Entonces, al tener el núcleo carga eléctrica positiva y nuestro electrón carga eléctrica negativa deberían atraerse por la fuerza coulombiana entre ellos. Cuando un cuerpo con carga se acelera emite radiación electromagnética (véase la fórmula de Larmor [aquí inglés]) perdiendo así energía.
Trayectoria en espacio fásico

En una dimensión, la posición del electrón, x, oscilaría y caería hacia el origen de coordenadas. Su velocidad, p, también oscilaría y lo llevaría hacia el centro, radiando energía en el proceso.

Según este modelo el electrón se vería atraído irremisiblemente hacia el núcleo y acabaría chocando con éste. Pero esto implicaría ¡que no habría electrones!. En jerga científica diríamos que la materia no sería estable. Puesto que nosotros somos materia y estamos aquí las hipótesis de nuestro modelo no pueden ser correctas.
En otra entrada explicaremos la descripción de los estados en mecánica clásica y cuántica, pero permitidnos ahora que hagamos algunos supuestos. En la mecánica clásica la evolución de un sistema se describe como la trayectoria de un punto en el estado de fases. Sin embargo, en la mecánica cuántica esto no es posible, dado que el estado de un sistema solo puede venir dado como una “probabilidad”, de modo que no existe trayectoria de un solo punto que se adecue a la evolución temporal de un estado. Aquí, pues, la incógnita será una función (con argumentos en el espacio-tiempo y con valores en los complejos) \Psi (x,y,z,t) cuyo módulo al cuadrado (que es un número real) nos dará la densidad de probabilidad de que un sistema se encuentre en un cierto estado.
Esta función, que llamaremos “función de onda”,  evoluciona de acuerdo con la ecuación de Schrödinger. Esta es una ecuación en derivadas parciales (EDP) “cualitativamente” hiperbólica (ver la entrada anterior ). Para mayor simplicidad vamos a suponer aquí que no depende del tiempo, con lo que ahora nuestra ecuación de Schrödinger es elíptica y se puede escribir como la función que minimiza un funcional de acción (ya hablamos antes de este concepto)  J(\Psi)= \int |\nabla \Psi(x)|^2 dx +\int V(s)|\Psi(x)|^2 dx
Si utilizamos, como representación de la función de onda,  la base de POSICIONES, la densidad de probabilidad asociada a nuestra función de onda, \int_A |\Psi(x)|^2dx,  nos dará la probabilidad de que el sistema (en este caso, nuestro electrón) esté en una zona A del espacio.
Si nos interesa el MOMENTO (i.e. la velocidad), su función de onda es la transformada de Fourier de la función de onda de posiciones, \hat{\Psi}(p),  y esta es la madre del cordero.
En la mecánica clásica la posición y el momento eran independientes. Por eso, podríamos realizar medidas sobre una y otra cantidad sin afectar a la otra.
Ahora, no obstante, tenemos una ligadura: hay una relación que conecta la posición y el momento, y ha de cumplirse. Podemos, pues, escribir nuestro funcional en términos de las dos funciones de onda J(\Psi,\hat{\Psi})= \int |\hat{ \Psi}(p)|^2|p^2| dp +\int V(s)|\Psi(x)|^2 dx.
La formulación del principio de incertidumbre de Heisenberg es, pues, \int x^2|\Psi(x)|^2 dx \int p^2|\hat{\Psi}(p)|^2 dp=\int |\nabla_x \Psi(x)|^2 dx \int |\nabla_p \hat{\Psi}(p)|^2 dp\geq C donde C es cierta constante que podemos escribir explícitamente pero que hacerlo nos da una pereza superlativa por lo que lo dejamos ‘para el lector interesado’ ;). La constante es una ensalada de constante de Planck \hbar, \pi, algún dos… Puede escribirse también involucrando sólo una función de onda de la misma manera que se hace con el funcional.
Hemos dicho que \Psi(x,t) nos da una probabilidad, con función de densidad |\Psi(x,t)|^2dx. Así observamos que en realidad el término \int x^2|\Psi(x)|^2 dx es la varianza de nuestra variable aleatoria |\Psi(x,t)|^2 (que asumimos tiene media x=0, es decir, que nuestro electrón está más o menos rondando el origen de coordenadas espaciales).
Entonces si entendemos esas ‘normas’ de los gradientes como los “errores” (utilizaremos esta palabra aquí a modo explicativo, pero hay que ser cautos con la terminología), veremos que un error pequeño en alguna de nuestras mediciones, para fijar ideas, la posición, implica que el error del momento (es decir, de la velocidad) tiene que ser alto para que el producto sea mayor que una constante. Es más: cuanto mejor midamos una (menor error) mayor es el error de la otra. ¡Qué vida más dura la nuestra!
No he hablado nada claro en este último párrafo, y los exigentes se habrán quedado pensando que me explico muy mal (quizá tengan razón). Hemos calculado un par de ejemplos ilustrativos usando python. Hemos supuesto que es unidimensional y que nuestras funciones de onda toman valores en los reales y no en los complejos. No lo hacen, pero sed generosos con nosotros que esto es sólo una afición.
Supongamos que nuestra onda es como la de esta figura
Función de onda con momento definido: seno

Si la funcion de onda es como un seno, la longitud de onda (o el momento) está bien definida.

Entonces por la fórmula de De Broglie, que conecta la visión de partículas con la de ondas, tenemos una velocidad definida (con poco o ningún error), que depende de la longitud de onda (¡otra entrada por hacer!). Esto de nuevo enlaza con nuestra transformada de Fourier, dado que este tipo de transformadas nos llevan desde el espacio de posiciones al espacio de momentos. Sin embargo, no conocemos nada de la posición, pues la probabilidad no se decanta por ninguna zona en particular y salvo puntos donde nunca estarán los electrones por lo demás no podemos decir nada. (Para el lector interesado, decir que esta función de onda no está bien definida en el espacio de posiciones: no es normalizable)
Otro caso interesante es algo así como:
Función de onda de un paquete gaussiano

Un paquete gaussiano es una elección "popular" de función de onda porque tiene incertidumbre mínima

En esta segunda imagen tenemos una función tal que al hacer el cuadrado obtenemos una zona que acapara casi toda la probabilidad, un entorno del origen (como la vida misma, unos pocos acaparan casi toda la proba…, digo, billetes). Es decir, casi sin error podemos saber su posición; sin embargo, no podemos usar la fórmula de De Broglie para calcular su velocidad, pues ¿quién me dice su longitud de onda?…
Para ir concluyendo ¿cómo enlaza esto con los gradientes (es decir, las derivadas)?, observamos que una función como la de la figura 1 tiene un gradiente ‘pequeño’, mientras que una función cómo la de la figura 2 tiene un gradiente ‘grande’. Sólo hay que ver que ese pico tiene una derivada bien grande. Así tenemos que entender que los gradientes me dan una idea del error, pero cambiando la variable. Es decir, una gradiente alto en las x me dice un error grande en las p y un error pequeño en las x. Aunque no hemos hablado de ello, esto tiene que ver con la estructura geométrica de las ecuaciones del movimiento Hamiltonianas (ya llegará, ya llegará…).
Es sorprendente (pero usual) cómo la física a veces acaba teniendo que ver con ideas abstractas de las matemáticas (al revés también ocurre).
Hay más que contar, pero se está insinuando una entrada muy larga, así que en la segunda parte de esta entrada escribiremos más acerca de los desarrollos matemáticos de la Transformada de Fourier y de si se puede obtener alguna propiedad física desde las matemáticas.
=======================
Esta entrada es participante en la XVIII Edición del Carnaval de la Física, alojado por “La Aventura , en la mecánica cuántica esto no es posible, de la Ciencia.

¿Podemos ir hacia atrás (matemáticamente) en el tiempo?

Me comentó David que sería interesante que explicase un poquito de lo que significa la reversibilidad temporal desde el punto de vista de las ecuaciones en derivadas parciales (EDPs de ahora en adelante). Es decir, sin mencionar nada de entropías. Este es un tema muy interesante e importante aunque a primera vista parezca una tontería. Lo hemos puesto como una serie de dos entradas. La primera (esta que estás leyendo) presentará las EDPs más fáciles y estudiará sus propiedades de cara a la reversibilidad temporal. En la segunda veremos la derivación de unos modelos más complicados y trataremos de entender dónde aparece la irreversibilidad temporal. Lo cierto es que en cierto modo es llamativo, pues las leyes de Newton son reversibles en tiempo y muchas (casi todas) de las ecuaciones de las que hablaremos surgen de ellas.

Viendo qué ocurre:
Para entender qué quiere decir y qué implica la reversibilidad (o irreversibilidad) temporal hemos de comprender primero los ejemplos más básicos de EDPs. Como dato para los técnicos supondré que todos los problemas están puestos para (t,x)\in [0,T]\times\mathbb{R}. Inmerso en el texto hay imágenes .jpg, imágenes animadas .gif (pinchad en ellas para que empiecen a moverse) y código Matlab (¡usadlo si podéis!).

La ecuación del transporte:

Comenzaremos con la ecuación del transporte unidimensional con coeficientes constantes. Esta es la EDP más sencilla que podemos poner.

\partial_t u + c \partial_x u=0, u(0,x)=f(x).

Supondremos que f es una función derivable con derivada continua una vez. La solución de esta ecuación es u(t,x)=f(x-ct) (¡comprobadlo!). Esta ecuación se llama ‘del transporte’ porque lo que hace es eso ‘mueve’ nuestra distribución inicial f. Si queremos cambiar el sentido del tiempo hemos de hacer el cambio t por -t. Entonces la nueva ecuación es

-\partial_t u + c \partial_x u=0, u(0,x)=f(x).

Observamos que el cambiar el tiempo de sentido es equivalente a cambiar el signo de c. Si utilizamos el código que propongo con varios valores de c observamos que este parámetro es una velocidad. Por lo tanto, parece natural que cambiar el sentido del tiempo cambie el sentido del movimiento. Es decir, que si para c>0 íbamos a la derecha, para tiempos negativos (o equivalentemente -c) tenemos que ir a la izquierda. Concluímos así que la ecuación del transporte es reversible en tiempo y que la reversibilidad es muy natural si partimos del proceso físico que se modela con esta ecuación.


Otra consecuencia, ésta mucho más sutil, de la reversibilidad temporal es que nuestra solución NUNCA va a ser mejor que nuestro dato inicial f. Esto es obvio en este caso porque tenemos una solución explícita, pero es cierto en general. Si u tuviese más derivadas que f entonces dando la vuelta al tiempo tendríamos una contradicción.

function [u,x,t]=transporte(dx,dt,f,c)
%%
%Funcion que me aproxima la solucion exacta (conocida) de la
%ecuacion del transporte u_t+cu_x=0 con dato inicial
% u(0,x)=f(x), paso espacial dx y paso temporal dt.
%f sera una funcion
% Rafael Granero Belinchon
%%

%Definicion de parametros:
T=10; %El tiempo final
t=0:dt:T; %el vector de tiempos
x=-pi:dx:pi; %el vector de espacio donde queremos
%nuestra aproximacion. No necesitamos condiciones de borde ¿por que?
u=zeros(length(x),length(t));
F=feval(f,x);

%Calculo de la solucion:
%La solucion del problema anterior es u(x,t)=f(x-ct)
u(:,1)=F;
for j=2:length(t)
u(:,j)=feval(f,x-c*t(j));
plot(x,u(:,j-1));%Representacion de los resultados
drawnow
end

function f=prueba(x)
f=sin(x);


La ecuación de ondas:

La siguiente ecuación es el paradigma de ecuación hiperbólica. Me refiero a la ecuación de ondas. Viendo el nombre está claro qué proceso físico quiere describir ¿no?.

Esta ecuación se escribe

\partial_t ^2 u= c^2\partial_x^2 u, u(0,x)=f(x), \partial _t(0,x)=g(x).

Visualmente parece mucho más complicada que la ecuación del transporte… sin embargo en realidad es igual (al menos en un cierto sentido). Vamos a escribirla como un sistema. Para ello definimos el sistema

\partial_t u=c\partial_x v, \partial_t v=c\partial_x u.

Si ahora derivamos en tiempo la ecuación para \partial_t u y utilizamos la segunda ecuación obtenemos

\partial_t^2 u=c\partial_x \partial_tv=c^2\partial_x ^2 u.

Es decir, que la ecuación de ondas no es más que dos transportes acoplados. Sin embargo todavía podemos hacerlo mejor. Podemos darnos cuenta de que el operador diferencial se puede escribir como

\partial_t^2-c^2\partial_x^2=(\partial_t + c\partial_x)(\partial_t-c\partial x),

y por lo tanto si tenemos u=u_1+u_2 con

\partial_t u_1 + c\partial_x u_1=0; \partial_t u_2 - c\partial_x u_2=0

tenemos una solución de la ecuación original. Concluímos que, como la ecuación del transporte era reversible, la ecuación de ondas, que se puede escribir como un par de ecuaciones del transporte debe ser reversible también.

La ecuación del calor:

Esta ecuación es parabólica. Se escribe

\partial_t u= \partial_x^2 u, u(0,x)=f(x).

Visualmente parece estar a medio camino entre la ecuación del transporte y la ecuación de ondas, sin embargo su comportamiento en radicalmente distinto. Para convencernos de ello podemos ‘jugar’ un poco con el código Matlab que adjunto. Los datos iniciales por defecto son los mismos, pero os animo a cambiarlos.

function [u,x,t,mx]=heatff(N,dt,K)
%%
%Funcion que me aproxima la solucion de la
%ecuacion del calor con dato inicial seno
%con condiciones periodicas (para usar FFT)
%N es el numero de nodos espaciales que se quieren
%usar. dt es el paso temporal que se quiere.
%K es la constante de difusion.
%Devuelve el espacio, el tiempo, la aproximacion de
%la solucion y el maximo de dicha solucion en cada tiempo
%Rafael Granero Belinchon
%%
T=5;%Tiempo final
dx=2*pi/(N-1);
x=-pi:dx:pi;%Espacio
t=0:dt:T;%Tiempo
uo=sin(x);%dato inicial
for k=1:N/2 %Operador laplaciano en espacio de fourier
L(k)=(k-1)*(k-1);
L(k+N/2)=(N/2-k+1)*(N/2-k+1);
end
L=K*L;
u(:,1)=uo’;
mx(1)=max(uo);
for l=1:length(t)
u(:,l+1)=ifft(exp(-L*dt*l).*fft(uo))’;%solucion
mx(l+1)=max(u(:,l+1));%Evolucion del maximo
plot(x,u(:,l));axis([-pi,pi,-1,1]);%Representacion de los resultados
drawnow
end
end

Para estudiar esta ecuación vamos a utilizar la transformada de Fourier. Para la transformada de Fourier de u(x) usaremos la notación \hat{u}(k). Así si transformamos la ecuación en espacio obtenemos las ecuaciones diferenciales ordinarias (EDOs a partir de ahora) indexadas en la longitud de onda k siguientes

\frac{d}{dt}\hat{u}(t,k)=-k^2\hat{u}(t,k), \hat{u}(0,k)=\hat{f}(k).

Esta ecuación la podemos resolver explícitamente

\hat{u}(t,k)=e^{-k^2t}\hat{f}(k).

Observemos ahora qué quiere decir el cambio del sentido del tiempo. De nuevo hagamos el cambio t por -t. La ecuación nos queda

\partial_t u= -\partial_x^2 u, u(0,x)=f(x).

No se ve nada claro, sin embargo, si buscamos los efectos del cambio en la solución explícita tenemos

\hat{u}(t,k)=e^{k^2t}\hat{f}(k),

de manera que cuando invirtamos la transformada de Fourier estamos calculando una convolución con una función que no está acotada, ni tiene ninguna potencia integrable… Vamos, que nuestra solución (que existe explícitamente) no está en ningún espacio razonable ni con propiedades físicas razonables. Por ejemplo, si u es la temperatura, entonces su integral (que es el calor) debe ser finita. Pues si damos la vuelta al tiempo obtenemos calor infinito para cualquier tiempo. Concluímos que la ecuación del calor NO es reversible en tiempo.

Una propiedad que a veces se da en las ecuaciones irreversibles y que es bien interesante es el ‘efecto regularizante’. Es decir, tu dato inicial f es continuo (por ejemplo), pero tu solución u es infinitamente derivable para todo tiempo (positivo). Como ya hemos mencionado antes, este comportamiento difiere del de las ecuaciones hiperbólicas usuales. La prueba de esto se puede hacer sin más que multiplicar por u e integrar por partes en espacio (¡comprobadlo!). Después basta observar que la ecuación es invariante por derivación tanto en tiempo como en espacio (¡Concluid el argumento!).

Vistos estos 3 ejemplos parece que hay una relación entre la ‘simetría’ del problema y su reversibilidad temporal. Quiero decir que, al menos de momento, las ecuaciones que tienen el mismo número de derivadas temporales que espaciales han resultado ser reversibles, mientras que las que no las tienen son irreversibles.

Otra cosa que se nos puede ocurrir es que las ecuaciones reversibles sean las que ‘no tiendan a nada’. Así vemos que la ecuación del calor tiende a ser idénticamente cero (necesita tiempo infinto para llegar a serlo) mientras que la ecuación del transporte sólo se movía por el espacio.

Veamos otro ejemplo:

La ecuación de Schrödinger:

Esta ecuación, clave en mecánica cuántica, se escribe

\partial _t u= i\partial_x^2 u, u(0,x)=f(x).

Si repetimos el análisis que hicimos para la ecuación del calor obtenemos que la solución es

\hat{u}(t,k)=e^{-ik^2t}\hat{f}(k),

que tiene un comportamiento oscilatorio. Por lo tanto, pese a tener una derivada en tiempo y dos en espacio se parece más a una ecuación de ondas que a una ecuación del calor. Observamos que el hecho de que aparezca la unidad imaginaria hace que u no sea real, sino compleja. Por lo tanto tiene una función conjugada. Si ahora cambiamos el sentido del tiempo observamos que para la función conjugada \bar{u} la ecuación es la misma. Por lo tanto si u es nuestra solución con el tiempo hacia delante, \bar{u} es una solución con el tiempo hacia atrás. Por lo tanto la ecuación de Schrödinger es reversible. Este ejemplo desmonta la hipótesis de que als reversibles debían tener el mismo número de derivadas en espacio y en tiempo.

function [u,x,t,L2,mx]=schrodinger(N,dt)
%%
%Funcion que me aproxima la solucion de la
%ecuacion de schrodinger con dato inicial seno
%con condiciones periodicas (para usar FFT)
%N es el numero de nodos espaciales que se quieren
%usar. dt es el paso temporal que se quiere.
%Devuelve el espacio, el tiempo, la aproximacion de
%la solucion, la norma L^2 de dicha solucion en cada tiempo
%y el maximo de la solucion en todo tiempo.
%Rafael Granero Belinchon
%%
T=5;%Tiempo final
dx=2*pi/(N-1);
x=-pi:dx:pi;%Espacio
t=0:dt:T;%Tiempo
uo=sin(x);%dato inicial
for k=1:N/2 %Operador laplaciano en espacio de fourier
L(k)=(k-1)*(k-1);
L(k+N/2)=(N/2-k+1)*(N/2-k+1);
end
L=i*L;
u(:,1)=uo’;
L2(1)=norm(uo)*dx;
mx(1)=max(abs (uo));
for l=1:length(t)
u(:,l+1)=ifft(exp(-L*dt*l).*fft(uo))’;%solucion
L2(l+1)=norm(u(:,l+1))*dx;%Evolucion de la norma L^2
mx(l+1)=max(abs(u(:,l+1)));%Evolución del máximo
plot(x,real(u(:,l)));axis([-pi,pi,-1,1]);%Representacion de los resultados
drawnow
end
end

Como dato anécdotico de esta ecuación hacemos notar que no puede reflejar efectos relativistas (¿por qué?).

En la próxima entrada, enlazando con esta, trataré la derivación desde la mecánica hamiltoniana de los modelos que se utilizan en mecánica de fluidos y trataré de explicar dónde aparece la irreversibilidad en el proceso. Sin embargo, esto es algo que no está ‘completamente entendido’ todavía.

La teoría del caos y el efecto mariposa

Para empezar sólo la voy a llamar teoría del caos en el título. Era para captar al atención, y este nombre lo consigue.  Es mejor decirle ‘de los sistemas dinámicos’ o algo así. Por desgracia, la fama de todo esto ha hecho que se encuentre de todo en internet. Por ejemplo hoy he visto una página donde se afirmaba que el efecto mariposa era los efectos de… ¡viajar en el tiempo! Impresionante. Y es que aficionarse a la ciencia está bien, pero sin perder el norte.

Por el nombre (caos) hay gente que piensa que un comportamiento caótico es aleatorio e impredecible. Vamos, un cisco bueno. Lo cierto es que lo primero es falso y lo segundo, como casi siempre, es ‘depende’.

El caos es ‘determinista‘, que quiere decir que dado un estado inicial, el comportamiento a largo plazo está determinado sin error posible y es único. Esto es, que está ‘determinado’ por el estado inicial. Consideremos un sistema discreto, es decir, una ley de recurrencia, por ejemplo la también famosa ley de Fibonacci, pero sin estados iniciales. Entonces la dinámica (la ley que sigue el sistema) es f_n=f_{n-1}+f_{n-2} Dados dos estados iniciales, por ejemplo 1,1 conocemos todos los valores de f. Además, si realizamos el experimento dos veces con los mismos valores iniciales los resultados serán idénticos. Eso quiere decir determinista.

Diferentes son los sistemas probabilistas. En estos sistemas hay un componente azaroso que impide conocer el largo plazo. Pero lo que de verdad los caracteriza es que para el mismo dato inicial podemos obtener resultados de los experimentos completamente distintos. Por ejemplo (muy poco válido como veremos ahora), una moneda. Consideramos el experimento tirar una vez la moneda. A mismas condiciones unas veces saldrá cara y otras cruz. Digo que es un mal ejemplo, porque este modelo es probabilista sólo por nuestro desconocimiento, pues si conociésemos la dirección y la fuerza exactas del lanzamiento sabríamos si saldrá cara o cruz. La gravedad es determinista. Esto nos podría llevar a pensar hasta que punto existe el azar, o si puede ser la probabilidad sólo una herramienta útil dada nuestra ignoracia de la realidad completa. Pensar en un mundo completamente determinista ya lo hizo Laplace. Y tiene una frase famosa por ello.

Ya hemos entendido la palabra determinista (si no es así tenéis que releerlo). Veamos el ‘depende’.

Es aquí donde entra el ‘efecto mariposa‘, que es el nombre que le puso Edward Lorenz a la sensibilidad a los datos iniciales. ¿A que es exótico?. Creo que su idea era atraer la atención hacia su conferencia. El efecto mariposa viene a decir que cualquier cambio minúsculo acaba teniendo repercusiones enormes, y por lo tanto nuestra aproximación (predicción) será una chapuza completa. O exóticamente

Si una mariposa batiese sus alas en Pekín provocaría un tornado en Texas un mes siguiente.

O algo parecido. Bueno, no matéis a todas las mariposas para evitar los tornados. No hay que cogerlo tan literal. En realidad hay que hacer una interpretación de casi todo. Las mariposas no provocan tornados. Los tornados surgen de un conjunto de factores que los hacen posibles, esto es, todas las mariposas del mundo, nuestros aviones, nosotros corriendo, la humedad en mi pueblo… Ahora, si pudiesemos tener dos planetas Tierra, con exactamente las mismas condiciones salvo una mariposa, entonces los climas serían distintos. ¿En qué sentido (ahora viene el depende)?. Bueno, no distintos en el sentido de que en Valencia helase por las noches. Llamemos a este tipo de cambios bruscos cambios de tipo 1. No, serían cambios en el orden y en el tiempo (tiempo-temporal, no tiempo-clima). Por ejemplo, un tornado que apareciese en un planeta el día 1 de Julio en el otro no aparecería y aparecería uno el 18 de Agosto. O una tormenta en mi pueblo no caería, caería en Albaladejo. Estos son los cambios de tipo 2.

Los cambios de tipo 1 son cambios bruscos que no quiere nadie. Estos cambios van asociados a cambios muy profundos en el sistema. Ahora tengo que ponerme técnico, lo siento. En un sistema dinámico, hay asociado un espacio de fases, que es un sitio donde viven las características del sistema. En este lugar de posiciones y velocidades (si es físico el sistema) o en el caso del clima de humedades y temperaturas, existen ‘cosas’ que atraen. Además son cosas raras normalmente en los casos caóticos. De hecho son fractales. Como no tenían muchas ganas de buscar un nombre exótico los llamaron ‘atractores extraños’. Posee la virtud de la simpleza. Los comportamientos del sistema cambian bruscamente entre unos atractores y otros (el mismo sistema puede tener varios al variar los parámetros). Es decir, para seguir con el ejemplo del clima, tenemos nuestro porcentaje de CO2 en la atmósfera en un nivel x. Nuestro sistema entonces tiene x como un parámetro. Si aumenta el porcentaje, digamos a 2x entonces nuestro sistema cambia de parámetro, pudiéndose producir un cambio de atractor, con el consiguiente nevazo en la Malvarrosa. Podemos ver esto del plano de fases y los atractores como un par de platos hondos y una aceituna, tenemos los platos juntos de manera que la aceituna reposa en el borde de ellos. La aceituna caerá rodando a uno de ellos. Los platos son los atractores, el conjunto de los dos es el espacio de fases y nuestra aceituna es el estado del sistema.

Dicho esto, está claro que hay que evitar los cambios de tipo 1. Los de tipo 2 son mejores en general.

Resumiendo, si tenemos un estado inicial, este nos viene dado por un conjunto de mediciones que hemos hecho. Estas mediciones no tienen (ni pueden) tener una precisión infinita (en cuyo caso el estado del sistema estaría determinado siempre) por lo que aparecen pequeños errores entre nuestro estado inicial para realizar los cálculos y la predicción y el verdadero estado inicial. Es decir, en nuestro espacio de fases hay dos puntos distintos, el de las medidas y el real. Al ser el comportamiento caótico, al avanzar el tiempo las curvas que tracen estos puntos se separarán. Y consecuentemente nuestra predicción a largo tiempo fallará.

Concluyendo, el ‘depende’ significa que podemos predecir el corto (quizá muy corto) plazo con poco error y podemos predecir si habrá cambio de atractor o no. Esto es comparable a decir que en Cuenca en Julio hará calor pero no poder decir habrá 36º C a las 15 de la tarde del día 2 de Julio. Decir que hará calor es decir el atractor, decir la temperatura exacta a la hora cabal es una predicción a largo plazo.

El comportamiento caótico existe, además es muy común. Aparece en todas las ramas del saber, física, biología medicina… Necesita algunas cosas para que se pueda dar. El sistema de ecuaciones diferenciales ha de ser no lineal y tener una dimensión mayor que dos. Pero esto no es nada raro en la cruda realidad fuera de las ‘oscilaciones pequeñas’ y cosas por el estilo.

He hecho un programita en Matlab para ver la sensibilidad a los datos iniciales:

function [x1,y1,z1,Y1,x2,y2,z2,Y2]=caoslogistica
%Este codigo estudia varios casos de 2 sistemas dinamicos discretos
%que conducen a un comportamiento caotico para ciertos valores
%de un parametro. Asi mismo dibuja unos diagramas de bifurcacion.
%Para el primer caso se puede poner k=4 para ver el comportamiento
%caotico. Para el segundo se puede poner k=1.6.
%Rafael Granero Belinchon.

disp(‘Comenzamos con el sistema dinamico Xn+1=k1Xn(1-Xn)’)
disp(‘—Primer Experimento:—‘)
k1=input(‘Dame una constante entre 0 y 4:’);
x0=input(‘Dame un valor inicial entre 0 y 1:’);
x1(1)=x0;
for i=1:100
x1(i+1)=k1*x1(i)*(1-x1(i));
end
disp(‘—Segundo Experimento:—‘)
k2=input(‘Dame una constante entre 0 y 4:’);
y0=input(‘Dame un valor inicial entre 0 y 1:’);
y1(1)=y0;
for i=1:100
y1(i+1)=k2*y1(i)*(1-y1(i));
end
disp(‘—Tercer Experimento:—‘)
k3=input(‘Dame una constante entre 0 y 4:’);
z0=input(‘Dame un valor inicial entre 0 y 1:’);
z1(1)=z0;
for i=1:100
z1(i+1)=k3*z1(i)*(1-z1(i));
end
subplot(4,1,1)
plot(x1);title(‘Primer experimento’);
subplot(4,1,2)
plot(y1);title(‘Segundo experimento’);
subplot(4,1,3)
plot(z1);title(‘Tercer experimento’);
subplot(4,1,4)
plot(z1,’r’);
hold on
plot(x1);
hold on
plot(y1,’k’);title(‘Todos juntos’);
hold off
a=input(‘Presiona cualquier tecla para continuar:’);
clear a;
disp(‘Comenzamos con el sistema dinamico Xn+1=k1Xn^2-1’)
disp(‘—Primer Experimento:—‘)
k1=input(‘Dame una constante entre 0 y 4:’);
x0=input(‘Dame un valor inicial entre 0 y 1:’);
x2(1)=x0;
for i=1:100
x2(i+1)=k1*x2(i)^2-1;
end
disp(‘—Segundo Experimento:—‘)
k2=input(‘Dame una constante entre 0 y 4:’);
y0=input(‘Dame un valor inicial entre 0 y 1:’);
y2(1)=y0;
for i=1:100
y2(i+1)=k2*y2(i)^2-1;
end
disp(‘—Tercer Experimento:—‘)
k3=input(‘Dame una constante entre 0 y 4:’);
z0=input(‘Dame un valor inicial entre 0 y 1:’);
z2(1)=z0;
for i=1:100
z2(i+1)=k3*z2(i)^2-1;
end
figure
subplot(4,1,1)
plot(x2);title(‘Primer experimento’);
subplot(4,1,2)
plot(y2);title(‘Segundo experimento’);
subplot(4,1,3)
plot(z2);title(‘Tercer experimento’);
subplot(4,1,4)
plot(z2,’r’);
hold on
plot(x2);
hold on
plot(y2,’k’);title(‘Todos juntos’);
hold off
a=input(‘Presiona cualquier tecla para continuar:’);
clear a;
disp(‘Vamos a dibujar ahora un diagrama de bifurcacion para el primer sistema’)
K=0:0.01:4;
X=zeros(length(K),5000);;
X(:,1)=0.3;
for j=1:length(K);
for i=1:5000
X(j,i+1)=K(j)*X(j,i)*(1-X(j,i));
end
end
Y1=X(:,4000:end);
figure
plot(Y1);title(‘Diagrama de bifurcacion para el primer sistema’)
a=input(‘Presiona cualquier tecla para continuar:’);
clear a;
disp(‘Vamos a dibujar ahora un diagrama de bifurcacion para el segundo sistema’)
K=0:0.01:4;
X=zeros(length(K),5000);
X(:,1)=0.3;
for j=1:length(K);
for i=1:5000
X(j,i+1)=K(j)*X(j,i)^2-1;
end
end
Y2=X(:,4000:end);
figure
plot(Y2);title(‘Diagrama de bifurcacion para el segundo sistema’)

Primer experimento (k=4)

La diferencia entre los datos iniciales es de 0.01 entre el de arriba y el segundo y de 0.001 entre el primero y el tercero. Como se puede ver en la gráfica de abajo, las trayectorias se separan más o menos por llegado un tiempo, pero al principio iban bien juntitas. Además podemos ver que todas estan en el mismo atractor. Eso se ve por los ‘patrones’ característicos que tienden a producirse. Me refiero a las oscilaciones grandes seguidas de varias muy pequeñas.  Esto se entiende fácil si se imagina uno nuestro estado del sistema como una mosca cansina. Estará dando vueltas a tu alrededor, quizá de forma complicada, pero antes o despues va a volver a pasar por tu oreja, quizá no como la vez anterior, pero muy cerca. Eso es lo que produce estos patrones. Que nuestro estado pasa cerca de una cierta parte del atractor.

Segundo experimento (k=1.6)

Arriba se muestran los resultados de ejecutar el código poniendo k=4 para el primer sistema y k=1.6 para el segundo. Veamos cómo cambia el carácter de la solución al mover k:

Otra cosa que resultaba ‘rara’ del comportamiento caótico era que reglas aparentemente muy sencillas (como puede ser elevar al cuadrado y restar 1) daban comportamientos muy complicados. También ocurre al revés, reglas aparentemente complicadas dan comportamiento simple.

Quizá otro día hable del sistema de Lorenz y de su atractor famoso.