Hace tiempo que la colección “Física con Excel” no aparece en esta web. Ahora es el momento de hacerlo, ya que no solamente ha comenzado el curso académico sino que tenemos un caso de candente actualidad donde nos va a venir muy bien.
El inductor a ello ha sido Michael Atiyah. Tranquilo, yo tampoco había oído hablar de él hasta ahora. Por lo que veo en la Wikipedia, la palabra genio se le queda pequeña. Ha ganado la medalla Fields, el premio Abel, y a los noventa años ha decidido saltar la banca y demostrar uno de los problemas del milenio: la conjetura de Riemann. Esto ya es en sí algo gordo, pero resulta que uno de los puntos en los que se apoya su demostración es aún más extraordinario, a mi parecer, y es la determinación de la constante de estructura fina. Esa constante relaciona diversas constantes físicas como la permitividad del vacío (ε0,) la carga del electrón (e), la velocidad de la luz en el vacío (c) y la constante de Planck (h/2π):
La constante de estructura fina α nos indica la fuerza de la interacción electromagnética entre partículas. Tiene dos particularidades interesantes. La primera es que se trata de una constante adimensional, es decir, es un número sin unidades; la segunda es que su valor no es muy grande ni muy pequeño, que es lo que suele pasar con las constantes físicas. Al principio se creía que alfa (la constante de estructura fina) era exactamente igual a 1/137. Ahora hemos medido mejor y tenemos un valor tal que así: 1/α=137,035 999 139…
No se sabe por qué tiene ese valor, pero la cosa es que la constante de estructura fina se ha intentado relacionar con dimensiones del Universo, en una tarea que tiene más de numerología que de ciencia verdadera. Ahora Michael Atiyah afirma que hay una ecuación que permite obtener la constante de estructura fina. De ser cierto, resultaría extraordinario, más aún (en mi opinión) que la propia conjetura de Riemann. Hay que recordar que α no se obtiene a partir de ecuaciones o relaciones sencillas, y para calcularla hace falta meter nada menos que electrodinámica cuántica.
Bien, para ser más correctos, Atiyah proporciona dos modos de calcular alfa. Uno de ellos, el difícil, involucra funciones y manipulaciones tan complicadas que en estos momentos los expertos están intentando averiguar si es algo serio o tan sólo una patochada. Una segunda forma, mucho más sencilla, pasa por calcular una sencilla integral llamada ч. Esa integral la multiplicamos por el número pi, la dividimos por la constante de Euler-Mascheroni (conocida como γ), y ya está. ¡Sencillo!
Bueno, vamos a poner todo esto en forma de ecuación:
Esa integral es la que nos va a dar guerra, pero vamos a domarla pronto. Vamos a reescribir esta ecuación de la siguiente forma:donde Aj es el paréntesis de la ecuación anterior, esto es, (1 – integral). Vamos a simplificar la tarea mediante una relación de recurrencia. Se puede comprobar fácilmente que A1=1, ya que la integral es cero (sus dos límites de integración coinciden). Para un valor de j distinto se puede demostrar que se cumple Aj=Aj-1-I1j-I2j, dondeEsas son las integrales que vamos a calcular. Para ello usaremos la regla de Simpson compuesta, que suele ir bastante bien y da buenos resultados. La idea es tomar el intervalo de integración, dividirlo en n trocitos y hacer esto:
Bien, pues vamos allá. Nuestra hoja de cálculo (como siempre, enlace al final del post) consta de tres hojas. La hoja 2 calcula I1 mediante la regla de Simpson. La columna A nos da el valor de j, y las columnas B-D tienen los valores de a, b y h, respectivamente (dichos valores dependen de j, claro). La fila 5 tiene coeficientes de la regla de Simpson (1, 4, 2, 4, 2…4, 1). Sumando todos los elementos de la columna 8, fila E-DA, obtenemos la integral I1 para j=1; la fila siguiente lo mismo para j=2, y así hasta que nos aburramos (yo paré en j=50). De modo similar, la hoja 3 calcula la integral I2.
Ahora volvamos a la hoja 1. Lo primero que quiero indicar es que hemos fijado el valor de N, el número de puntos de integración, en la celda E1. He permitido un máximo de cien; si quieres aumentar esa cifra, tendrás que incluir más columnas en las hojas 2 y 3, pero no hay problema. De modo similar, si quieres calcular para valores de j mayores de 50 no hay más que incluir más filas en la hoja 1.
Los valores de I1, I2 para cada valor de j están en las columnas E y F respectivamente. A partir de ahí, obtenemos la integral Aj (columna D), multiplicamos por el coeficiente 1/2j (que está en la columna B), y finalmente tenemos todos los términos de la sumatoria hecha (columna C). Hacemos la suma y dividimos por dos en la celda H5 (ojo: si hacéis j mayor que 50, tenéis que hacer la suma de todos esos términos, no sólo de esos 50 como he hecho yo).
Ya solo tenemos que multiplicar por pi (celda H2), dividir por gamma (celda H1) y ya tenemos el valor de la constante de estructura fina α. Por comodidad haremos 1/α (celda H4), a ver si nos hemos acercado mucho al valor ese de 1/137 y pico. ¡Ni de coña! Tenemos en nuestra hoja un valor de α = 6,239867885. Nos hemos columpiado tres órdenes de magnitud. Tela marinera.
Tenemos un problema, eso es indudable. Quizá sea la integración, que no esté bien hecha, pero creo que sí lo está. Lo sé porque, al integrar con N=98 y luego con N=100 obtengo casi el mismo resultado con una diferencia de una parte en diez mil millones. ¿Acaso me he equivocado y he tomado el logaritmo neperiano? No, porque me he cuidado mucho de hacerlo bien. La función LN(A) te devuelve el logaritmo neperiano de la celda A, pero si quieres usar otra base para los logaritmos puedes hacerlo mediante la función LOG. Yo he hecho LOG(A,2) para obtener el logaritmo en base 2.
Resulta, por otro lado, que tengo confirmación de haberlo hecho bien. Un hilo en Reddit, aparecido precisamente sobre esta cuestión, muestra valores de α muy similares a los míos: 6,239867897. Han usado phyton3 (ojo: phyton2 da valores incorrectos, ya que al parecer divide mal los números enteros). Daniel Duque, desde el frente tuitero, obtiene por su parte un valor de 6,239867897632185. Que a todos nos salga prácticamente el mismo número (salvo por pequeños cambios en las últimas cifras decimales) no creo que sea casualidad.
Eso es lo que hay por el momento. En cuanto a la segunda forma, más complicada, me declaro incompetente. De hecho, parece que es de ahí de donde Atiyah sacó su conclusión sobre la conjetura de Riemann, y por lo que he visto los matemáticos la están haciendo picadillo (ver p. ej. este artículo de @twalmar), así que no pienso ni intentarlo. Me conformo con ver que una ecuación que hubiera partido la pana de ser cierta puede comprobarse fácilmente (en este caso, refutarse fácilmente) con una sencilla hoja de cálculo).
ACTUALIZACIÓN 26/9/18: Me dicen por el pinganillo (gracias, Daniel Duque) que la integral que aparece en la ecuación de marras es integrable. Esto quiere decir que no hay que meterse en líos de integración numérica como he hecho yo. Al integrar analíticamente, resulta un valor para la constante de estructura fina de 6,23986789761156 para una sumatoria hasta j=50. Sumando algunos términos más tenemos el valor final de 6,23986789763231. Sigue estando muy lejos del valor real.
[Hoja de cálculo disponible: Constante de estructura fina de Atiyah]
No es que python2 divida mal los numeros enteros, es que la operacion q realiza es distinta segun sean los numeros en coma flotante o enteros. Si los numeros son de tipo coma flotante, hace la division convencional, y sale un numero con decimales. Pero si ambos numeros son de tipo entero, hace la operacion division entera, devuelve un entero. Es decir,
4.0/3.0 = 1.3333..
4/3 = 1
Este comportamiento puede ser adecuado en algunos casos, pero ha originado grandes desastres inesperados, asi que en python3 han decidido ser coherentes. El operador / hace siempre la division convencional, intependientemente del tipo de operandos. Con los operadores // y % se hace siempre la division entera y el resto.
Python fue creado por un matematico, y tiene muchos tipos numericos, como numeros decimales, racionales o complejos.
Código en Matlab/Octave en doble precisión (la integral es evaluable analíticamente)
gamma = 0.57721566490153286060651209008240243104215933593992;
g = @(j) 1-((log(j) + 1)./j +j.*(log(j) – 1))/log(2); %1- int(log2(x),x,1/j,j)
j = 1:1e6; j = j(end:-1:1);
a = g(j);
alphaInv = (pi/(2*gamma)* sum(a.*2.^(-j)))
alpha = 1./alphaInv
Me devuelve alpha = 6.23986789763232
python2 divide perfectamente números enteros, lo q pasa es el operador “/” es interpretado de forma diferente en python2 y python3: en python3, “/” es siempre división verdadera (7/2==3.5; 7.0/2.0==3.5), pero en python2 “/” cuando los opera con números enteros, es división de parte entera (7/2==3; 7.0/2.0==3.5).
Es más, si a un script de python2 (a partir de python2.5, me suena) se añade la siguiente línea al inicio del script, tendrá el mismo comportamiento q python3:
from __future__ import división