Ir al contenido principal

Nicolás' blog

Reto de métodos numéricos: Día 2 – Regula falsi

Nicolás Guarín-Zapata

Comentarios

Durante octubre (2017) estaré escribiendo un programa por día para algunos métodos numéricos famosos en Python y Julia. Esto está pensado como un ejercicio, no esperen que el código sea lo suficientemente bueno para usarse en la "vida real". Además, también debo mencionar que casi que no tengo experiencia con Julia, así que probablemente no escriba un Julia idiomático y se parezca más a Python.

Regula falsi

El segundo método a considerar es el método de posición falsa, o regula falsi. Este método se utiliza para resolver la ecuación \(f(x) = 0\) para \(x\) real, y \(f\) continua. Se comienza con un intervalo \([a, b]\), donde \(f (a)\) y \(f (b)\) deben tener signos opuestos. El método es similar al método de bisección pero en lugar de dividir a la mitad el original intervalo toma como nuevo punto del intervalo la intersección de la línea que conecta la función evaluada en los dos puntos extremos. Entonces, el nuevo el punto se calcula a partir de

\begin{equation*} c = \frac{a f(b) - b f(a)}{f(b) - f(a)} \end{equation*}

Como en el método de bisección, mantenemos el intervalo donde la solución aparece, es decir, donde cambia el signo de la función.

Usaremos la función \(f(x) = \cos(x) - x\) para probar los códigos, y como intervalo inicial \([0.5, \pi/4]\).

A continuación se presentan los códigos.

Python

from __future__ import division, print_function
from numpy import abs, cos, pi

def regula_falsi(fun, a, b, niter=50, ftol=1e-12, verbose=False):
    if fun(a) * fun(b) > 0:
        c = None
        msg = "The function should have a sign change in the interval."
    else:
        for cont in range(niter):
            qa = fun(a)
            qb = fun(b)
            c = (a*qb - b*qa)/(qb - qa)
            qc = fun(c)
            if verbose:
                print("n: {}, c: {}".format(cont, c))
            msg = "Maximum number of iterations reached."
            if abs(qc) < ftol:
                msg = "Root found with desired accuracy."
                break
            elif qa * qc < 0:
                b = c
            elif qb * qc < 0:
                a = c
    return c, msg

def fun(x):
    return cos(x) - x

print(regula_falsi(fun, 0.5, 0.25*pi))

Julia

function regula_falsi(fun, a, b, niter=50, ftol=1e-12, verbose=false)
    if fun(a) * fun(b) > 0
        c = nothing
        msg = "The function should have a sign change in the interval."
    else
        for cont = 1:niter
            qa = fun(a)
            qb = fun(b)
            c = (a*qb - b*qa)/(qb - qa)
            qc = fun(c)
            if verbose
                println("n: $(cont), c: $(c)")
            end
            if abs(fun(c)) < ftol
                msg = "Root found with desired accuracy."
                break
            elseif qa * qc < 0
                b = c
            elseif qb * qc < 0
                a = c
            end
            msg = "Maximum number of iterations reached."
        end
    end
    return c, msg
end

function fun(x)
    return cos(x) - x
end

println(regula_falsi(fun, 0.5, 0.25*pi))

Comparación

Respecto al número de líneas de código tenemos: 29 en Python y 32 en Julia. La comparación en tiempo de ejecución se realiza con el comando mágico de IPython %timeit y con @benchmark en Julia.

Para Python:

%timeit regula_falsi(fun, 0.5, 0.25*pi)

con resultado

10000 loops, best of 3: 35.1 µs per loop

Para Julia:

@benchmark regula_falsi(fun, 0.5, 0.25*pi)

con resultado

BenchmarkTools.Trial:
  memory estimate:  48 bytes
  allocs estimate:  2
  --------------
  minimum time:     449.495 ns (0.00% GC)
  median time:      464.371 ns (0.00% GC)
  mean time:        493.785 ns (0.52% GC)
  maximum time:     9.723 μs (92.54% GC)
  --------------
  samples:          10000
  evals/sample:     198

Reto de métodos numéricos: Día 1 – Método de bisección

Nicolás Guarín-Zapata

Comentarios

Por el próximo mes estaré escribiendo un programa por día para algunos métodos numéricos conocidos en Python y Julia. Está destinado a ser un ejercicio, entonces no espere que el código sea lo suficientemente bueno para un uso real. Además, debo mencionar que casi no tengo experiencia con Julia, por lo que probablemente no será Julia idiomático sino Julia más parecido a Python.

Día 1: Método de bisección

El primer método a considerar es el método de bisección. Este método se usa para resolver la ecuación \(f(x) = 0\) para \(x\) real, y \(f\) continua. Se empieza con un intervalo \([a,b]\), donde \(f(a)\) y \(f(b)\) deben tener signos opuestos. El método procede partiendo a la mitad el intervalo y seleccionando el subintervalo en donde aparece la solución, es decir, el signo de la función cambia.

Usaremos la función \(f(x) = \cos(x) - x^2\) para probar los códigos, y el intervalo inicial es [0, 1].

A continuación están los códigos:

Python

from __future__ import division, print_function
from numpy import log2, ceil, abs, cos

def bisection(fun, a, b, xtol=1e-6, ftol=1e-12):
    if fun(a) * fun(b) > 0:
        c = None
        msg = "The function should have a sign change in the interval."
    else:
        nmax = int(ceil(log2((b - a)/xtol)))
        for cont in range(nmax):
            c = 0.5*(a + b)
            if abs(fun(c)) < ftol:
                msg = "Root found with desired accuracy."
                break
            elif fun(a) * fun(c) < 0:
                b = c
            elif fun(b) * fun(c) < 0:
                a = c
            msg = "Maximum number of iterations reached."
    return c, msg

def fun(x):
    return cos(x) - x**2

print(bisection(fun, 0.0, 1.0))

Con resultado

(0.824131965637207, 'Maximum number of iterations reached.')

Julia

function bisection(fun, a, b, xtol=1e-6, ftol=1e-12)
    if fun(a) * fun(b) > 0
        c = nothing
        msg = "The function should have a sign change in the interval."
    else
        nmax = ceil(log2((b - a)/xtol))
        for cont = 1:nmax
            c = 0.5*(a + b)
            if abs(fun(c)) < ftol
                msg = "Root found with desired accuracy."
                break
            elseif fun(a) * fun(c) < 0
                b = c
            elseif fun(b) * fun(c) < 0
                a = c
            end
            msg = "Maximum number of iterations reached."
        end
    end
    return c, msg
end

function fun(x)
    return cos(x) - x^2
end

println(bisection(fun, 0.0, 1.0))

Con resultado

(0.824131965637207, "Maximum number of iterations reached.")

En este caso, ambos códigos están bastante cerca. El código de Python tiene 25 líneas, mientras que el de Julia 27. Como se esperaba, los resultados son iguales.

Edición (2017-10-02)

Como sugirió Edward Villegas, decidí comparar los tiempos de ejecución. Usé %timeit para IPython y @benchmark (de BenchmarkTools) para Julia.

En IPython, tenemos

%timeit bisection(fun, 0.0, 2.0)

con resultado

10000 loops, best of 3: 107 µs per loop

Y en Julia, tenemos

@benchmark bisection(fun, 0.0, 2.0)

con resultado

BenchmarkTools.Trial:
  memory estimate:  48 bytes
  allocs estimate:  2
  --------------
  minimum time:     1.505 μs (0.00% GC)
  median time:      1.548 μs (0.00% GC)
  mean time:        1.604 μs (0.00% GC)
  maximum time:     38.425 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     10

Parece que la versión de Julia es alrededor de 100 veces más rápida que su equivalente en Python.

Solución de la ecuación de Schrödinger usando el método de Ritz

Nicolás Guarín-Zapata

Comentarios

En esta publicación describimos la solcuión de la ecuación de Schrödinger usando el método de Ritz y base de funciones de Esta base parece ser una buena elección para la ecuación de Schrödinger ya que es una base ortogonal sobre \((-\infty, \infty)\).

Transformando la ecuación en una algebraica

Queremos resolver la ecuación de Schrödinger independiente del tiempo

\begin{equation*} \left[-\frac{1}{2}\nabla^2 + V(x)\right] \psi = E\psi\, , \end{equation*}

donde estamos usando unidades naturales ya que son una buena elección para cálculos numéricos.

Resolver esta ecuación es equivalente a resolver la siguiente ecuación variacional

\begin{equation*} \delta \Pi[\psi] = 0\, , \end{equation*}

con

\begin{equation*} \Pi[\psi] = \frac{1}{2} \langle \nabla \psi, \nabla\psi\rangle + \langle \psi, V(x) \psi\rangle - E\langle \psi, \psi\rangle \, , \end{equation*}

con \(\psi\) la función de onda, \(V(x)\) el potencial y \(E\) la energía. Esta formulación variacional es equivalente a la ecuación de Schrödinger independiente del tiempo, y \(E\) funciona como un multiplicador de Lagrange que garantiza que la probabilidad sobre todo el dominio sea 1.

Podemos expandir la función de onda en la base ortogonal

\begin{equation*} \psi = \sum_{n=0}^N c_n u_n(x)\, , \end{equation*}

donde \(u_n(x) \equiv \mu_n H_n(x) e^{-x^2/2}\) es una función de Hermite normalizada, \(\mu_n\) es el inverso de la magnitud del \(n\)-ésimo polinomio de Hermite

\begin{equation*} \mu_n = \frac{1}{\sqrt{\pi^{1/2} n! 2^n}}\, , \end{equation*}

y \(c_n\) son los coeficientes de la expansión. Esta representación es exacta en el límite \(N \rightarrow \infty\).

Si remplazamos la expansión en el funcional, obtenemos

\begin{equation*} \Pi_N = \sum_{m=0}^N\sum_{n=0}^N c_m c_n\left[ \frac{1}{2} \langle \nabla u_m, \nabla u_n\rangle + \langle u_m, V(x) u_n\rangle - E^N \delta_{mn}\right]\, . \end{equation*}

El integrando que involucra las dos derivadas sería

\begin{align*} u_m' u_n' =& \left[2m \frac{\mu_{m-1}}{\mu_m}u_{m-1} - x u_m\right] \left[2n \frac{\mu_{n-1}}{\mu_n}u_{n-1} - x u_n\right]\\ =& 4mn\frac{\mu_{m-1} \mu_{n-1}}{\mu_m \mu_n} u_{n-1} u_{m-1} - 2m\frac{\mu_{m-1}}{\mu_{m}}x u_{m-1} u_n\\ &- 2n\frac{\mu_{n-1}}{\mu_{n}}x u_{n-1} u_m + x^2 u_m u_n \end{align*}

Entonces, el término de la energía cinética sería

\begin{align*} \langle \nabla u_m, \nabla u_n \rangle =& 4mn\frac{\mu_{m-1} \mu_{n-1}}{\mu_m \mu_n} \langle u_{n-1}, u_{m-1}\rangle - 2m\frac{\mu_{m-1}}{\mu_{m}} \langle u_{m-1}, x u_n\rangle\\ &- 2n\frac{\mu_{n-1}}{\mu_{n}} \langle u_{m}, x u_{n - 1}\rangle + \langle u_m, x^2 u_n\rangle\\ =& 4mn \frac{\mu_{m-1}^2}{\mu_m^2}\delta_{mn} - 2m \frac{\mu_{m-1}}{\mu_m} \alpha_{m-1, n} - 2n \frac{\mu_{n-1}}{\mu_n} \alpha_{m, n-1} + \beta_{mn} \, , \end{align*}

con

\begin{equation*} \alpha_{m,n} \equiv \langle u_{m}, x u_n\rangle = \begin{cases} \sqrt{\frac{n + 1}{2}} & m=n +1\\ \sqrt{\frac{n}{2}} & m=n - 1\\ 0 & \text{otherwise}\end{cases}\, , \end{equation*}

y

\begin{equation*} \beta_{m,n} \equiv \langle u_{m}, x^2 u_n\rangle = \begin{cases} \frac{\sqrt{n(n-1)}}{2} & m = n - 2\\ \frac{2n + 1}{2} & m = n \\ \frac{\sqrt{(n+1)(n+1)}}{2} & m = n + 2 \\ 0 & \text{otherwise}\end{cases}\, . \end{equation*}

El funcional se reescribe como

\begin{align*} \Pi_N =& \sum_{m=0}^N \sum_{n=0}^N c_m c_n \left[2mn \frac{\mu^2_{m-1}}{\mu^2_m}\delta_{mn} - m\frac{\mu_{m-1}}{\mu_m}\alpha_{m - 1, n} - n\frac{\mu_{n-1}}{\mu_n}\delta_{m, n-1} \right.\nonumber \\ &\left. - \frac{1}{2}\beta_{mn} + \langle u_m, V(x) u_n\rangle - E^N \delta_{mn}\right] \, . \end{align*}

Tomando la variación

\begin{equation*} \delta \Pi_N = 0\, , \end{equation*}

que en este caso es equivalente a

\begin{equation*} \frac{\partial \Pi_N}{\partial c_i} = 0\quad \forall i=0, 1, \cdots N\, , \end{equation*}

lleva a

\begin{equation*} [H]\lbrace c\rbrace = E^N\lbrace c\rbrace \, , \end{equation*}

donde las componentes de la matriz \([H]\) están dadas por

\begin{equation*} H_{mn} = 2mn \frac{\mu^2_{m-1}}{\mu^2_m}\delta_{mn} - m\frac{\mu_{m-1}}{\mu_m}\alpha_{m - 1, n} - n\frac{\mu_{n-1}}{\mu_n}\delta_{m, n-1} - \frac{1}{2}\beta_{mn} + \langle u_m, V(x) u_n\rangle\, . \end{equation*}

La última integral se puede calcular usando la cuadratura de Gauss-Hermite. Y necesitaremos más puntos de Gauss si queremos integrar polinomios de orden alto. Este método funciona bien para funciones que pueden ser aproximadas por polinomios.

Ejemplos

Una implementación en Python de este método se puede encontrar en este repo.

Para todos los ejemplos usamos las siguientes importaciones

from __future__ import division, print_function
import numpy as np
from hermite_QM import *

Oscilador armónico cuántico

En este caso el potencial (normalizado) está dado por

\begin{equation*} V(x) = \frac{1}{2} x^2 \end{equation*}

y los valores propios exactos normalizados son

\begin{equation*} E_n = n + \frac{1}{2} \end{equation*}

El siguiente bloque de código calcula los primeros 10 valores propios y grafica los estados propios correspondientes

potential = lambda x: 0.5*x**2
vals, vecs = eigenstates(potential, nterms=11, ngpts=12)
print(np.round(vals[:10], 6))
fig, ax = plt.subplots(1, 1)
plot_eigenstates(vals, vecs, potential);

con resultados

[ 0.5  1.5  2.5  3.5  4.5  5.5  6.5  7.5  8.5  9.5]
/images/hermite_ritz_harmonic.svg

Potencial de valor absoluto

potential = lambda x: np.abs(x)
vals, vecs = eigenstates(potential)
print(np.round(vals[:10], 6))
fig, ax = plt.subplots(1, 1)
plot_eigenstates(vals, vecs, potential, lims=(-8, 8));

con resultados

[ 0.811203  1.855725  2.57894   3.244576  3.826353  4.38189   4.895365
  5.396614  5.911591  6.421015]
/images/hermite_ritz_abs.svg

Potencial cúbico

potential = lambda x: np.abs(x/3)**3
vals, vecs = eigenstates(potential)
print(np.round(vals[:10], 6))
fig, ax = plt.subplots(1, 1)
plot_eigenstates(vals, vecs, potential, lims=(-6, 6));

con resultados

[ 0.180588  0.609153  1.124594  1.681002  2.272087  2.889805  3.530901
  4.191962  4.871133  5.566413]
/images/hermite_ritz_cubic.svg

Oscilador armónico con perturbación cuártica

potential = lambda x: 0.5*x**2 + 0.1*x**4
vals, vecs = eigenstates(potential, nterms=20, ngpts=22)
print(np.round(vals[:10], 6))
fig, ax = plt.subplots(1, 1)
plot_eigenstates(vals, vecs, potential, lims=(-5, 5))

con resultados

[  0.559146   1.769503   3.138624   4.628884   6.220303   7.899789
   9.658703  11.489094  13.38638   15.361055]
/images/hermite_ritz_pert_harm.svg

Un notebook de Jupyter con los ejemplos se puede encontrar acá.

Usando estilos predefinidos en matplotlib

Nicolás Guarín-Zapata

Comentarios

Podemos dar formato a los gráficos de forma simple usando el paquete de estilo en matplotlib. La idea principal es crear un archivo con algunos parámetros predefinidos (esto también puede hacerse a través de rcParams).

Esta publicación no es un tutorial en cómo usar estos archivos, para estos puedes mirar la página style sheet reference. Acá, quiero jugar un poco con algunos de los parámetros para crear tres estilos diferentes. Los primeros dos tienen el estilo de un software (infame para algunos), que es usado por la mayoría de las personas como su plataforma de visualización, el tercero es un estilo limpio. Todos los archivos usados se pueden descargar aquí.

Para todos los ejemplos se realizan las siguientes importaciones:

from __future__ import division, print_function
import numpy as np
import matplotlib.pyplot as plt

Primer ejemplo: MS 2003

En nuestro primer ejemplo queremos reproducir el estilo que solíamos ver como la opción por defecto hace más de una década.

Este es el contenido del archivo MS2003.mplstyle

font.family : sans-serif

axes.facecolor : c0c0c0
axes.edgecolor : black
axes.prop_cycle : cycler('color',['000080', 'FF00FF', 'FFFF00', '00FFFF','800080', '800000', '008080', '0000FF'])
axes.grid : True

axes.spines.left   : True
axes.spines.bottom : True
axes.spines.top    : True
axes.spines.right  : True

grid.color : black
grid.linestyle : -

lines.linewidth : 1

figure.figsize : 5, 3

legend.fancybox : False
legend.frameon : True
legend.facecolor : white
legend.edgecolor : black
legend.loc : center left

El siguiente bloque de código usa este estilo

style = "MS2003.mplstyle"
with plt.style.context(style):
    x = np.linspace(0, 4, 100)
    y = np.sin(np.pi*x + 1e-6)/(np.pi*x + 1e-6)
    fig = plt.figure()
    ax = plt.subplot(111)
    for cont in range(5):
        plt.plot(x, y/(cont + 1), label=cont)

    plt.gca().xaxis.grid(False)
    box = ax.get_position()
    ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
    plt.legend(bbox_to_anchor=(1, 0.5))

y este es el resultado

/images/MS2003_style.svg

Segundo ejemplo: MS 2007

En el segundo ejemplo queremos reproducir la prole del primer estilo en este ejemplo. Este estilo es una mejora respecto al anterior, pero no es perfecto.

El siguiente es el contenido del archivo MS2007.mplstyle

font.family : sans-serif

axes.facecolor : white
axes.edgecolor : 4d4d4d
axes.prop_cycle : cycler('color',['4573a7', 'aa4644', '89a54e', '71588f','4298af', 'db843d', '93a9d0', 'd09392'])
axes.grid : True
axes.linewidth : 0.5

axes.spines.left   : True
axes.spines.bottom : True
axes.spines.top    : False
axes.spines.right  : False

lines.linewidth : 2

grid.color : 4d4d4d
grid.linestyle : -
grid.linewidth : 0.5

figure.figsize : 5, 3

legend.fancybox : False
legend.frameon : False
legend.facecolor : white
legend.edgecolor : 4d4d4d
legend.loc : center left

El siguiente código usa este estilo

style = "MS2007.mplstyle"
with plt.style.context(style):
    x = np.linspace(0, 4, 100)
    y = np.sin(np.pi*x + 1e-6)/(np.pi*x + 1e-6)
    fig = plt.figure()
    ax = plt.subplot(111)
    for cont in range(5):
        plt.plot(x, y/(cont + 1), label=cont)

    plt.gca().xaxis.grid(False)
    box = ax.get_position()
    ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
    plt.legend(bbox_to_anchor=(1, 0.5))

y este es el resultado

/images/MS2007_style.svg

Tercer ejemplo: un estilo limpio

El último ejemplo es un estilo limpio que usa una paleta de colores tomada de ColorBrewer.

Este es el contenido del archivo clean_style.mplstyle

font.family : sans-serif

axes.facecolor : white
axes.prop_cycle : cycler('color',['e41a1c', '377eb8', '4daf4a', '984ea3', 'ff7f00', 'ffff33', 'a65628', 'f781bf'])
axes.linewidth : 0.0
axes.grid : True

lines.linewidth : 1.5

xtick.direction : in
ytick.direction : in

grid.color : c7dedf
grid.linestyle : -
grid.linewidth : 0.3

figure.figsize : 6, 4

legend.fancybox : False
legend.frameon : False
legend.loc : best

El siguiente código usa este estilo

style = "clean.mplstyle"
with plt.style.context(style):
    x = np.linspace(0, 4, 100)
    y = np.sin(np.pi*x + 1e-6)/(np.pi*x + 1e-6)
    fig = plt.figure()
    ax = plt.subplot(111)
    for cont in range(5):
        plt.plot(x, y/(cont + 1), label=cont)

    plt.legend()

y este es el resultado

/images/clean_style.svg

También podemos usar archivos que están almacenado remotamente. Por ejemplo, podríamos usar la siguiente URL:

style = "https://raw.githubusercontent.com/nicoguaro/matplotlib_styles/master/styles/clean.mplstyle"

Recursos

Como mencioné anteriormente, el objetivo de esta publicación era crear algunos archivos de estilo simples para matplotlib y verlos en acción.

Esta publicación no permite ser una guía para buenos gráficos/visualizaciones. Para este propósito sugiero mirar la siguiente referencia:

Además, encuentro muy útiles las siguientes herramientas:

  • ColorBrewer2 permite elegir mapas de colores con diferentes criterios (cuantitativo/cualitativo, apto para impresión, apto para daltónicos).

  • ColRD tiene muchas paletas de colores. También permite generar paletas a partir de imágenes.

  • Colorgorical es una herramienta para crear paletas de colores categóricas (cualitativas) para visualización de información.

Interpolación de Hermite a trozos

Nicolás Guarín-Zapata

Comentarios

En esta entrada encontramos las funciones de interpolación de Hermite para el dominio [-1, 1]. Luego, las usamos para una interpolación a trozos. Observe que esta interpolación tiene continuidad \(C^1\) y no \(C^0\) como en la interpolación de Lagrange.

Para calcular los polinomios explícitamente podemos usar sympy.

Queremos encontrar una base con funciones que sastisfagan

\begin{align*} N_1(x_1) &= 1\\ N_1(x_2) &= 0\\ N_2(x_1) &= 1\\ N_2(x_2) &= 0\\ N'_3(x_1) &= 1\\ N'_3(x_2) &= 0\\ N'_4(x_1) &= 1\\ N'_4(x_2) &= 0 \end{align*}

Esto puede escribirse como

\begin{equation*} \begin{bmatrix} 1 &x_1 &x_1^2 &x_1^3\\ 1 &x_2 &x_2^2 &x_2^3\\ 0 &1 &2x_1 &3x_1^2\\ 0 &1 &2x_2 &3x_2^2 \end{bmatrix} \begin{bmatrix} a_{11} &a_{12} &a_{13} &a_{14}\\ a_{21} &a_{22} &a_{23} &a_{24}\\ a_{31} &a_{32} &a_{33} &a_{34}\\ a_{41} &a_{42} &a_{43} &a_{44} \end{bmatrix} = \begin{bmatrix} 1 &0 &0 &0\\ 0 &1 &0 &0\\ 0 &0 &1 &0\\ 0 &0 &0 &1 \end{bmatrix} \end{equation*}

La fórmula para la interpolación sería

\begin{equation*} f(x) \approx N_1(x) u_1 + N_2(x) u_2 + |J|(N_3(x) u'_3 + N_4(x) u'_4)\quad \forall x\in [a, b] \end{equation*}

con \(|J| = (b - a)/2\) el determinante jacobiano de la transformación.

\begin{equation*} \left[\begin{matrix}1 & x_{1} & x_{1}^{2} & x_{1}^{3}\\ 1 & x_{2} & x_{2}^{2} & x_{2}^{3}\\ 0 & 1 & 2 x_{1} & 3 x_{1}^{2}\\ 0 & 1 & 2 x_{2} & 3 x_{2}^{2}\end{matrix}\right] \end{equation*}

Podemos ver que la matriz anterios e una matriz de Vandermonde confluente [1]. Es imilar a una matriz de Vandermonde para los primeros \(n\) nodos y las derivadas de cada columna para los siguientes.

Los coeficientes de nuestros polinomios están dados por la inversa de esta matriz.

\begin{equation*} \left[\begin{matrix}\frac{x_{2}^{2} \left(3 x_{1} - x_{2}\right)}{x_{1}^{3} - 3 x_{1}^{2} x_{2} + 3 x_{1} x_{2}^{2} - x_{2}^{3}} & \frac{x_{1}^{2} \left(x_{1} - 3 x_{2}\right)}{x_{1}^{3} - 3 x_{1}^{2} x_{2} + 3 x_{1} x_{2}^{2} - x_{2}^{3}} & - \frac{x_{1} x_{2}^{2}}{x_{1}^{2} - 2 x_{1} x_{2} + x_{2}^{2}} & - \frac{x_{1}^{2} x_{2}}{x_{1}^{2} - 2 x_{1} x_{2} + x_{2}^{2}}\\- \frac{6 x_{1} x_{2}}{x_{1}^{3} - 3 x_{1}^{2} x_{2} + 3 x_{1} x_{2}^{2} - x_{2}^{3}} & \frac{6 x_{1} x_{2}}{x_{1}^{3} - 3 x_{1}^{2} x_{2} + 3 x_{1} x_{2}^{2} - x_{2}^{3}} & \frac{x_{2} \left(2 x_{1} + x_{2}\right)}{x_{1}^{2} - 2 x_{1} x_{2} + x_{2}^{2}} & \frac{x_{1} \left(x_{1} + 2 x_{2}\right)}{x_{1}^{2} - 2 x_{1} x_{2} + x_{2}^{2}}\\\frac{3 x_{1} + 3 x_{2}}{x_{1}^{3} - 3 x_{1}^{2} x_{2} + 3 x_{1} x_{2}^{2} - x_{2}^{3}} & - \frac{3 x_{1} + 3 x_{2}}{x_{1}^{3} - 3 x_{1}^{2} x_{2} + 3 x_{1} x_{2}^{2} - x_{2}^{3}} & - \frac{x_{1} + 2 x_{2}}{x_{1}^{2} - 2 x_{1} x_{2} + x_{2}^{2}} & - \frac{2 x_{1} + x_{2}}{x_{1}^{2} - 2 x_{1} x_{2} + x_{2}^{2}}\\- \frac{2}{x_{1}^{3} - 3 x_{1}^{2} x_{2} + 3 x_{1} x_{2}^{2} - x_{2}^{3}} & \frac{2}{x_{1}^{3} - 3 x_{1}^{2} x_{2} + 3 x_{1} x_{2}^{2} - x_{2}^{3}} & \frac{1}{x_{1}^{2} - 2 x_{1} x_{2} + x_{2}^{2}} & \frac{1}{x_{1}^{2} - 2 x_{1} x_{2} + x_{2}^{2}}\end{matrix}\right] \end{equation*}

Y remplazamos los valores para los nodos extremos, -1 y 1.

\begin{equation*} \left[\begin{matrix}\frac{1}{2} & \frac{1}{2} & \frac{1}{4} & - \frac{1}{4}\\ - \frac{3}{4} & \frac{3}{4} & - \frac{1}{4} & - \frac{1}{4}\\ 0 & 0 & - \frac{1}{4} & \frac{1}{4}\\ \frac{1}{4} & - \frac{1}{4} & \frac{1}{4} & \frac{1}{4}\end{matrix}\right] \end{equation*}

Los polinomios son el producto de los coeficientes y los monomios

\begin{equation*} \left[\begin{matrix}\frac{1}{4} \left(x - 1\right)^{2} \left(x + 2\right)\\- \frac{1}{4} \left(x - 2\right) \left(x + 1\right)^{2}\\\frac{1}{4} \left(x - 1\right)^{2} \left(x + 1\right)\\\frac{1}{4} \left(x - 1\right) \left(x + 1\right)^{2}\end{matrix}\right] \end{equation*}

La base interpolante sería

\begin{align*} N_1 (x) &= \frac{1}{4} (x - 1)^2 (2 + x)\\ N_2 (x) &= \frac{1}{4} (x + 1)^2 (2 - x)\\ N_3 (x) &= \frac{1}{4} (x - 1)^2 (x + 1)\\ N_4 (x) &= \frac{1}{4} (x + 1)^2 (x - 1)\, , \end{align*}

y la siguiente función calcula la interpolación para una función y derivada dadas

En este caso, interpolamos la función sinc

El siguiente bloque de código calcula la interpolación y la grafica.

/images/sinc_interp.svg

Referencias

Herramientas (auxiliares) para Investigación

Nicolás Guarín-Zapata

Comentarios

En esta entrada quiero hablar de algunas herramientas (...o algo así) que son útiles en el día a día de la investigación, pero que no suelen ser tan populares por ser en cierta forma tangenciales al área específica en la que cada uno trabaja.

Scripting

Un script  es un programa usualmente simple que realiza una serie de tareas. ¡Lo siento! trato de ceñirme al castellano pero no me gustan las traducciones para este término (archivo de órdenes, archivo de procesamiento por lotes o guion). Si existe una tarea que debe ser hecha más de... cinco veces (el número varía de acuerdo a la paciencia), entonces creo que es algo que podemos pedirle al computador que haga esta tarea por nosotros. En otras palabras: podemos automatizar ese trabajo. Algunas tareas que pueden ser una buena alternativa para automatizar son: renombrar 100 archivos, convertir un archivo de un formato a otro (p. ej., STL a OBJ), leer 387 archivos con información sobre el clima y graficar su evolución temporal (temperatura mínima, máxima y promedio). Estas tareas pueden ser fáciles de hacer a mano, pero por la cantidad de trabajo que involucran son tediosas.

Lo primero es hacernos con un lenguaje de scripting. Algunas opciones son Python, Bash, Julia, Matlab/Octave, Scilab. Dejándome llevar por el sesgo, recomiendo usar Python.

Gráficos y esquemas

Una imagen vale más que mil palabras, o eso dice el dicho. Personalmente, me parece absolutamente cierto y trato de hacer dibujitos para poder entender mejor algo o lograrlo explicar mejor. Lo primero que me gustaría mencionar es la diferencia entre imágenes de mapas de bits (o ráster) e imágenes vectoriales.

  • Imagen de mapa de bits: es una imagen que está representada por un arreglo (o rejilla rectangular) de pixeles. Dicho de otro modo, se almacena la información de color que hay en cada punto de la imagen. Los formatos más populares almacenan la información comprimida. Para gráficos con alto contraste (como esquemas o diagramas) el mejor formato es PNG. Si se tiene una animación, GIF sería preferible. Y para el caso de fotografías es mejor utilizar JPG.

  • Imagen vectorial: es una imagen que está formada por entidades geométricas. En esta no se almacena la información punto a punto sino la construcción de las formas que la constituyen. Por esta razón, estas imágenes no se pixelan pues la información que se tiene es de cómo construirla. Este tipo de imágenes es la mejor opciones para esquemas y diagramas, pues sólo se almacena la información de los trazos y texto que se añadan en ellas (ver Figura 1). El estándar de facto para este tipo de imágenes es PDF —es el que suelo para incluir en mis documentos \(\LaTeX\), aunque existe forma de embeber SVG en \(\LaTeX\) pero es algo que aún no exploro—. A pesar de que PDF sea el estándar, el formato a preferir es SVG (Scalable Vector Graphics) que es un estándar a través de la internet y la mayoría de navegadores modernos permiten visualizar.

Recapitulando, deberíamos usar imágenes JPG para fotografías y SVG para esquemas/diagramas. Otro atributo que puede ser de utilidad es el manejo de capas, SVG permite esto.... y de formatos ráster tenemos la opción de usar TIFF.

En cuanto a software para generar/editar este tipo de imágenes debo decir que existen gran cantidad de programas que permiten exportar a estos formatos: Python/Matplotlib, Matlab, Inkscape, Adobe Illustrator, GIMP, Photoshop, LibreOffice. Si el gráfico es generado a partir de un cálculo o de una serie de datos yo uso Matplotlib. Si en cambio, queremos hacer un esquema o —como suelo decir— un muñeco mi herramienta favorita es Inkscape. Este programa pretende ser una alternativa libre a programas como Adobe Illustrator —y sí que lo consigue. Obviamente, se podría usar Illustrator o Corel Draw para esta tarea. Si el único uso sería hacer esquemas técnico, creo que sería un desperdicio.

/images/sensor_hall.png

Figura 1. Esquema de sensado de campo magnético de imanes permanentes por efecto Hall. Sacado de mi trabajo de grado de Ingeniería Física.

Tomar notas

Supongo que para algunos parecería un poco trivial hablar de "tomar notas" y mucho más viniendo de alguien que no tenía cuadernos en el bachillerato, pero como soy un poco terco creo que igual escribiré un poco sobre esto. Lo primero que me gustaría mencionar es que recuerdo que me hablaran de este tema en el colegio, pero nunca hubo nada formal respecto a desarrollar estas habilidades. Buscando en la web, hay gran cantidad de información. Incluso el artículo de Wikipedia en inglés es interesante. No hay nada mejor para escribir que tener una pluma y un papel de buen gramaje, por eso es que aún uso un cuaderno en donde llevo cuenta de lo que hago en mi investigación y tomo notas. Sin embargo, este esquema es bastante lineal y deja por fuera oportunidades más contemporáneas. Es decir ¿por qué conformarse con un documento en este tiempo de hiper-documentos? Las ventajas de tomar notas digitalmente saltan a la vista, en un hiper-documento se pueden tener enlaces, embeber imágenes, video y sonido.

En cuanto a herramientas, acá incluyo una pequeña lista:

  • Evernote es probablemente la herramienta más popular para tomar notas. Es multiplataforma, Freemium (funcionalidad básica gratuita y avanzada paga), y cuenta con muchas opciones. Yo la uso, pero no mucho en mi investigación.

  • Zim es una wiki offline. Tiene gran cantidad de opciones como calendario, ecuaciones con código \(\LaTeX\), imágenes... en fin. El pero que le encuentro es que no he logrado configurar las ecuaciones en Windows (y en mi oficina debo usar Windows :-/).

  • Docear esta es una herramienta pensada, principalmente, para manejar bibliografía. Sin embargo, permite tomar notas y, en general, manejar la información de la investigación. La característica más (o menos) atractiva es que funciona alrededor de mapas mentales.

  • Zotero también es una herramienta para manejar bibliografía, aunque permite manejar algo de toma de notas (al menos alrededor de la bibliografía).

  • Mendeley es muy parecido al anterior, aunque con más funcionalidades. El pero más grande que le encuentro es que en 2013 fue comprado por Elsevier.

En cuanto a manejo de bibliografía también me gustaría mencionar EndNote que es el programa con mayor trayectoria y  JabRef que es el que he usado por más tiempo. Algunas referencias interesantes comparando manejadores de bibliografía están acá: [A] [B] [C].

Reconstrucción de gráficos

Es común encontrarse con información presentada en forma de gráficos. También es común que queramos tener los datos numéricos de estos gráficos para poder compararlos con los nuestros. Para saber si nuestras medidas/simulaciones/métodos dan resultados similares a otros presentados en la literatura. Para esta tarea se pueden usar potentes software de procesamiento de imágenes, u otros más modestos diseñados específicamente para este fin.

/images/pointplot.jpg

Figura 2. Gráfico original.
/images/digitized_pts.png

Figura 3. Gráfico procesado en Engauge Digitizer. Algunos puntos fueron seleccionados (automáticamente) para obtener sus coordenadas.

  • Digitizer of XY chart este es un plugin para Libreoffice/OpenOffice y exporta el resultado a la hoja de cálculo actual, es simple y fácil de usar.

  • Engauge Digitizer, es el que normalmente uso cuando necesito hacer esta tarea (ver Figuras XX). Es de código abierto (y libre) y tiene una buena cantidad de opciones para facilitar la tarea.

  • Plot Digitizer no tengo mucha información sobre este (pues nunca lo he usado), salvo que está escrito en Java.

  • ImageJ este es un (completo) programa de procesamiento de imágenes que está escrito en Java. No lo he usado para esta tarea de manera regular, pero podría usarse para ello.

Visualización científica

La visualización científica que se encarga de generar gráficos que permitan visualizar "datos científicos" para facilitar el entendimiento que hay detrás de los datos. Para esta labor muchos hemos usado lenguajes de scripting como Matlab/Octave, Scilab o Python (con Matplotlib o Mayavi). Sin embargo, como la visualización se trata de algo visual —¿como si no?—, es bueno contar con una herramienta que permita generar y cambiar los gráficos de manera interactiva, aunque siempre debemos automatizar la mayor cantidad de trabajo posible (la pereza siempre ha sido uno de los móviles más grandes de la humanidad, hay que aceptarlo).

  • MayaVi, este es un programa escrito en Python que usa VTK. Es una herramienta muy versátil y la gran ventaja que tiene es que puede usarse dentro de scripts de Python.

  • Paraview, este programa también está basado en VTK y permite paralelizar las labores (para los computadores con múltiple núcleo y los clústers). Abajo incluyo un video generado en Paraview para mostrar sus capacidades.

  • Visit, este programa también está basado en VTK, nunca lo he usado pero quise incluirlo porque la gente dice que puede ser más intuitivo que Paraview.

  • Tecplot, este programa es muy popular en Purdue. Creo que inicialmente fue pensado para CFD, pero se ha expandido mucho. En cuanto a gráficos 3D no me parece mejor que ParaView, sin embargo, las capacidades de gráficos 2D (gráficos XY, y demás) lo hacen atractivo.

  • Scavis, este está escrito en Java. No lo conocía hasta que inicie la escritura de esta entrada pero me llamó la atención y quise incluirlo en la lista. Algo que me llamó la atención es que permite scripting en varios lenguajes: Java, Python, Ruby, BeanShell y Matlab/Octave.

  • Origin, nunca lo he usado pero no lo quise dejar por fuera porque siempre he oído hablar maravillas de él (probablemente compárandolo con Excel... pero no puedo opinar).

Simulación de fuego cruzado en Vimeo.

Manejo de versiones

El manejo de versiones es la administración de cambios en documentos, código fuente y otro tipo de información. Esto puede hacerse de forma manual, pero es fácil cometer errores o remplazar la versión de un código fácilmente, y por esto es recomendable usar un software que facilite el trabajo. La idea es tener un lugar (repositorio) en donde se almacenan las versiones y los cambios, y llevar un registro de estos. De esta forma se puede volver a una versión anterior de los documentos y varias personas pueden trabajar conjuntamente.  Existen dos paradigmas (o arquitecturas) para el manejo de versiones: centralizada y distribuida. En la primera existe un repositorio centralizado en donde se encuentra toda la información. En la arquitectura distribuida cada usuario tiene una copia del respositorio. Personalmente sólo he usado Git, que entra en la categoría distribuida y es uno de los software de manejo de versiones más populares actualmente; lo usan compañías como Google, Facebook y Netflix.

Un ejemplo puede verse en este repositorio, en donde está el documento de trabajo de grado de Santiago Echeverri, el cual tuve la oportunidad de asesorar. Este documento lo editamos conjuntamente mientras él estaba en Medellín  y yo me encontraba en Estados Unidos. El documento se hizo en el lenguaje de marcadores \(\LaTeX\).

Además de tener un control sobre las versiones y poder acceder a versiones anteriores, es útil poder almacenar la información en un lugar accesible desde cualquier lugar del mundo con una conexión a internet. Esto puede lograrse con un servidor propio, obviamente, o también a través de un proveedor externo. Dos proyectos que  son muy populares para alojar repositorios y manejar sus versiones son (comparación entre Github y BitBucket):

  • Github  es el más popular en este momento. Permite tener proyectos con un número ilimitado de colaboradores. Para tener un repositorio privado es necesario pagar.

  • BitBucket tiene como principal ventaja que permite tener repositorios privados sin la necesidad de pagar. Sólo es gratuito para proyectos con 5 colaboradores o menos (o para proyectos académicos).

Enlaces sugeridos

  1. Software Carpentry. http://software-carpentry.org/

  2. Python Scientific Lecture Notes. https://scipy-lectures.github.io/

Seguro dejé mucho temas por fuera así como herramientas dentro de algún tópico. Si ese es el caso, agradecería que me lo digan en los comentarios.

La zorra y las uvas

Nicolás Guarín-Zapata

Comentarios

Mirando en internet no encontré una versión en español de la fábula de Esopo "La zorra y las uvas" que me gustara (o que al menos rimara :P). De hecho, el artículo de Wikipedia en español es bastante más malo que su equivalente en inglés.  Por esto, y por petición de una amiga, hice mi versión con rima asonante. Y reza así

La Zorra y las uvas

La zorra que con ansias uvas buscaba
no pudo alcanzarlas por lo alto que estaban.
Giró y exclamó, con risa lastimera:
¡Están verdes! No valen mi espera.

Tomando notas de clase con mi tableta

Nicolás Guarín-Zapata

Comentarios

El año pasado el niño dios me regaló una Samsung Galaxy Note 10.1 (modelo 2012). Esta tableta viene diseñada para tomar notas, pues trae un S-Pen (fabricado por Wacom, que creo que fabrican la tecnología para la mayoría de tabletas comerciales) y varias aplicaciones pensadas para ello. Como: S Note, PhotoShop Touch, Crayon Physics (tienen que jugarlo es muy divertido).

Decidí entonces hacer un experimento tomando notas de clase con mi tableta completamente en las dos materias que cursé en el "semestre" de primavera de 2014. En Optimización usé el lapiz que viene con la tableta y para Física Molecular usé un Bamboo Stylus que, se supone, es diseñado específicamente para esta tableta. Las notas las tomé con el software por defecto (S Note) y convertí a PDF (las notas fueron hechas enteramente con el software):

Aunque la sensación al tacto del Bamboo Stylus es mucho más agradable, el desempeño del S Pen es mucho mejor. En ambos casos, la sensación es muy diferente a la de usar papel y lápiz (lapicero/pluma), en gran medida porque la fricción es muy poca en comparación (algo que es mejor con los lapices conductores). La experiencia de usuario del lápiz varía mucho entre una aplicación y otra, lo que genera buenas expectativas. En cuanto a tomar notas: es muy recomendado para tomar notas cortas (como reuniones); no tanto para tomar notas largas (como notas de clase), a menos que se quiera tener una copia digital fácilmente.

Apps

La tecnología detrás de la tableta es la parte más importante. De manera simple y rápida, el lapiz tiene una antena y la tableta una rejilla de antenas que están sintonizadas en la misma frecuencia (531 kHz, en este enlace explican su funcionamiento un poco mejor). Sin embargo, también es importante tener en cuenta el software que se usa para tomar nota:

S Note

  • Es la aplicación nativa para tomar notas.

  • Tiene reconocimiento de formas y ecuaciones. Las ecuaciones deben ser simples, no esperen tomar notas de Mecánica Cuántica o Mecánica del Medio Continuo usado esta herramienta.

  • Permite insertar imágenes y grabar sonidos dentro de la aplicación.

  • Tiene "palm rejection".

  • El tamaño de papel no es configurable.

Papyrus

  • Tiene "palm rejection" y permite configurar el dedo como borrador (que puede ser útil).

  • Permite tener varios tamaños de papel y diferentes fondos.

  • Permite hacer rectángulos y elipses.

  • Los pinceles son muy básicos.

Note Anytime

  • Permite tener varios tamaños de papel.

  • Almacena los trazos vectorialmente y se pueden editar (cambiar color, tamaño o estilo) posteriormente.

  • Tiene varios pinceles y además son configurables.

  • No puede convertir las notas en PDF.

  • No tiene "palm rejection".

Quill

  • Es open source, aunque si lo adquieren a través de la "Google Store" tiene costo.

  • El algoritmo de reconocimiento de trazos es mejor que el de los otros software (y es ajustable).

  • Aún está en una etapa preliminar.

Lástimosamente no hay una versión de Paper (de FiftyThree) para Android, y, aunque existe una versión de Bamboo Paper (de Wacom) para Android no es compatible con esta tableta.

En conclusión

Remplazar la toma de notas "tradicionales" por notas  "digitales" puede ser un poco apresurado, pero es realizable. La aplicación que viene por defecto (S Note) permite realizar el trabajo "out of the box", pero podría ser mejor. Las aplicaciones que me parecen más promisorias son: Note Anytime y Quill.

Palabras terminadas en "ar", "er" o "ir"

Nicolás Guarín-Zapata

Comentarios

Cuando estaba en el colegio, en ocasiones mencionaban las palabras terminadas en "ar", "er" o "ir" —iguales al final de nuestros verbos en infinitivo— que no eran verbos. Eran pocos ejemplos los que mencionaban, así que algunos años después me dio el arrebato de hacer una lista… y hela aquí:

AR

AR

AR

ER

IR

agar

dispar

pajar

alfiler

cachemir

ajuar

ejemplar

paladar

alquiler

casimir

alar

escolar

palatar

amater

elixir

albañar

espectacular

palomar

anteayer

emir

albar

estelar

papilar

antier

faquir

alfar

extracurricular

par

atelier

kefir

alizar

familiar

paramilitar

ayer

menhir

almiar

fular

particular

bachiller

nadir

almimbar

globular

peculiar

bereber

porvenir

alminar

hangar

pendular

brigadier

sir

altar

hogar

peninsular

canciller

souvenir

alveolar

ijar

perpendicular

cualquier

visir

antealtar

impar

pilar

doquier

zafir

antisolar

insular

pinar

dossier

ar

intercelular

planar

enser

auricular

intraocular

polar

ester

avatar

jaguar

policelular

menester

axilar

juglar

pulgar

mercader

azahar

lagar

pulmonar

micer

azar

lanar

quasar

neceser

bacillar

lar

radar

palier

bar

limonar

radicular

placer

bazar

llar

rectangular

plumier

bicelular

lobular

reticular

premier

bienestar

lodazar

samovar

primer

bifilar

lugar

secular

rosicler

biliar

lunar

seglar

salmer

billar

malar

semicircular

somier

binocular

maleolar

semilunar

sumiller

bipolar

malvar

similar

taller

broncopulmonar

mandibular

solar

tercer

bulevar

manglar

subcelular

veguer

calamar

manillar

supraclavicular

calcañar

manjar

supralunar

capilar

mar

telar

capsular

maxilar

testicular

caviar

melocotonar

tisular

celular

melonar

trifilar

centenar

membrillar

trufar

ciliar

miliar

tubular

clavicular

militar

tular

collar

millar

ultramar

consular

molar

uvular

coplanar

mollar

vascular

crepuscular

muladar

vehicular

corpuscular

monocelular

vermicular

cuadrangular

mular

vesicular

cular

muscular

vulgar

curricular

nuclear

yugular

dactilar

ocular

zar

dinar

olivar

Si tienes alguna sugerencia de palabra que haga falta la puedes dejar en los comentarios.

Última actualización: Octubre 21, 2022.

Descargar la lista en formato CSV o en formato XLS.

Respecto a la teletransportación en Star Trek

Nicolás Guarín-Zapata

Comentarios

Alguna vez una amiga me preguntó sobre cuál era más fácil de dos métodos para teletransportar personas. Uno consistía en "desconfigurar" y "reconfigurar" a la persona (como información, similar a Start Trek y la otra en "crear un atajo" entre dos lugares del espacio-tiempo para pasar a la persona. Obviamente, en ambos casos suponiendo que tales cosas se pueden hacer.

Respecto a la creación del atajo en el espacio-tiempo, creo que no he estudiado las ecuaciones de Relatividad General para poder hacer las cuentas, aunque tal vez en Física Pasión nos regalen estos cálculos.

La primera de las alternativas consiste, básicamente, en la representación de una persona como información y su procesamiento de "desconfiguración", transmisión y procesamiento de "reconfiguración". Ya que los humanos estamos hechos, mayormente, de agua asumiremos una persona de agua para los cálculos. Además tomaremos como masa de referencia \(70\ \mbox{kg}\).

Algunos datos que usaremos son:

  • La masa molar del agua es \(M_{H_2O}=18,01528\ \mbox{g/mol}\);

  • El número de Avogadro es \(N_A=6,022 \times 10^{23}/\mbox{mol}\).

Entonces tendríamos que el número de moléculas en una persona es

\begin{equation*} m_\text{persona} = \frac{70\ \mbox{kg}\ \mbox{moles}} {18,015\times 10^{-3} \text{kg}} =3885,7\ \mbox{moles}\, , \end{equation*}

y el número de moléculas sería

\begin{equation*} N_\text{moleculas} = N_A m_\text{persona} = 2,34\times 10^{27}\ \mbox{moleculas}\, . \end{equation*}

El agua tiene 12 modos vibracionales y 6 traslacionales, y además 40 números cuánticos electrónicos. Lo que suma un total de 58 grados de libertad por cada molécula. Así que el número total de grados de libertad es

\begin{equation*} N_{GDL} = 1,36\times 10^{29}\, . \end{equation*}

En 2011 hubo un record de transmisión de 186 Gb/s, considerando esta tasa el tiempo que tardaría en transmitirse la totalidad de los datos (asumiendo datos representados como reales de 32 bits) tardaría \(2,3\times 10^{16}\ \mbox{s}\) o, equivalentemente, \(1,70\times 10^{10}\) años (la edad estimada del universo es de \(1,37\times 10^{10}\) años).

Si quisiéramos que esto pasara realmente rápido, digamos en 1 milisegundo, la tasa de transmisión debería ser \(2,3\times 10^{19}\) más rápida que el récord que alcanzaron en el CERN… algo que es inconcebible para nosotros actualmente.