Códigos para generar curvas

  • Categoría de la entrada:Artículos

La semana pasada di un curso en Tenerife sobre cómo explicar conceptos matemáticos de forma visual y sencilla. Uno de los temas que tratamos fue la teoría local de curvas planas, y aunque el nombre suena chungo, es la base para crear animaciones como esta:

Quiero compartir los código de Python con los que generé la animación anterior. La idea es crear cada fotograma por separado y juntarlos luego con algún programa de edición de vídeo. Suelo hacerlo así para poder editar los colores o el fondo, pero se puede exportar directamente en formato vídeo modificando alguna línea del código que os adjunto.

Lo primero es lo primero: importar las bibliotecas.

import matplotlib.pyplot as plt
import numpy as np

Una curva cualquiera

Una curva plana se puede entender como la trayectoria que siguen un punto que ese mueve por el plano. Ese punto tendrá coordenadas x(t) y y(t), que cambiarán según avanza el tiempo t\in\mathbb{R}^{+}. En general, escribimos una curva en el plano como:

c(t) = \left(x(t), y(t)\right)

En la animación que abre este artículo, el símbolo de la lemniscata, podemos describir la trayectoria como aquella que tiene de coordenadas

x(t) = \frac{\cos t}{1+\sin^2t}, ~~~~~~~~~~y(t) = \frac{\sin t\cos t}{1+\sin^2t}.

A la hora de escribir el código que la genera, podemos tomar t\in[0,2\pi]; el resto es una repetición periódica.

En código que adjunto, podéis elegir las expresiones que siguen las dos coordenadas de una curva cualquiera. También tenéis que seleccionar el intervalo [a,b] para la t, y el número de fragmentos en los que dividimos este intervalo. Cuantos más fragmentos, más suave se verá la curva.

# Dividir el intervalo [a, b] en n trozos de igual longitud
n = 300
a = 0
b = 2*np.pi

# Añadir dos puntos adicionales en los extremos para poder calcular las segundas derivadas
a -= 2*(b-a)/n
b += 2*(b-a)/n
n = n+4
delta = (b-a)/n

# Funciones de las coordenadas
def x(t):
    return np.cos(t)/(1+np.sin(t)**2)
    
def y(t):
    return np.sin(t)*np.cos(t)/(1+np.sin(t)**2)

# Curva total
def c():
    t = np.linspace(a, b, n)
    xl, yl = x(t), y(t)
    return xl, yl

# Derivada de la curva
def dc(x,y):
    return np.gradient(x)/delta, np.gradient(y)/delta

Parámetros de la curva

A partir de los códigos anteriores podemos calcular los parámetros típicos de la curva.

# Curva
x, y = c()

# Velocidad
dx,dy = dc(x,y)

# Vector tangente
velocity = np.array([[dx[i], dy[i]] for i in range(dx.size)])
ds = np.sqrt(dx*dx + dy*dy)
tan = np.array([1/ds] * 2).transpose() * velocity
tan_x = tan[:, 0]
tan_y = tan[:, 1]

#Vector normal
deriv_tan_x = np.gradient(tan_x)/delta
deriv_tan_y = np.gradient(tan_y)/delta
dT_dt = np.array([ [deriv_tan_x[i], deriv_tan_y[i]] for i in range(deriv_tan_x.size)])
length_dT_dt = np.sqrt(deriv_tan_x * deriv_tan_x + deriv_tan_y * deriv_tan_y)
normal = np.array([1/length_dT_dt] * 2).transpose() * dT_dt

# Adicionales
d2s_dt2 = np.gradient(ds)/delta
d2x_dt2 = np.gradient(dx)/delta
d2y_dt2 = np.gradient(dy)/delta

curvature = np.abs(d2x_dt2 * dy - dx * d2y_dt2) / ((dx * dx + dy * dy)**1.5)
t_comp = np.array([d2s_dt2] * 2).transpose()*tan
n_comp = np.array([curvature * ds * ds] * 2).transpose()*normal

En primer lugar, calculamos los puntos de la curva (representados en amarillo en la imagen). Justo después, computamos las derivadas de cada coordenada, que utilizaremos para hacer cálculos más adelante.

En segundo lugar, calculamos el vector tangente a la curva en cada punto. Este vector es el que aparece en rojo en la imagen. Lo normal es calcularlo como un vector unitario:

\vec T(t) = \frac{c'(t)}{\|c'(t)\|}

En la animación, yo utilizo una fórmula ligeramente distinta, con un factor que indica qué velocidad lleva el punto en cada momento:

\vec T(t) = \frac{c'(t)}{\|c'(t)\|}~v(t)

Así podemos hacernos una idea de lo rápido que se desplaza el punto mirando la longitud del vector. Creo que es más útil y más visual.

En tercer lugar, calculamos el vector normal a la curva en cada punto. Este es el vector que representamos en azul en la imagen. La dirección del vector normal es perpendicular a la curva y (de momento) supondremos que es unitario:

\vec N(t) = \frac{\vec T'(t)}{\|\vec T'(t)\|}

Finalmente, calculamos algunos parámetros de la curva, como la curvatura o las aceleraciones normal y tangencial. Por destacar una de ellas, la curvatura se calcula como:

k(t) = \frac{x'y'' - y'x''}{\left[ (x')^2+(y')^2  \right]^{3/2}}

Crear el gráfico de la función

Ya tenemos todos los ingredientes para crear la animación. Como ya he dicho, este código genera y guarda cada fotograma por separado. Luego hay que juntarlos con algún programa externo. Hay tantos fotogramas como fragmentos en los que hemos dividido el intervalo [a,b].

Habrás visto que mis gráficos no se parecen a los típicos de Python. Eso es porque borro todo el formato por defecto y genero los elementos visuales uno a uno, desde el fondo a los ejes coordenados. En mi caso, uso una biblioteca personal que creé hace años para hacer los gráficos bonitos. En vuestro caso, he hecho un apaño para que podáis tener un formato similar al mío:

# Límite del gráfico
lim_eje_x = [-2.05, 2.05]
lim_eje_y = [-2.05, 2.05]

for f in range(2, n-2):
    #Borrar el formato de matplotlib
    plt.gca().set_aspect('equal')
    plt.axis('off')
    
    # Plotear los ejes coordeandos
    plt.plot([-2,2],[0,0],linewidth=1,color='#C0C0C0',solid_capstyle='round')
    plt.plot([0,0],[-2,2],linewidth=1,color='#C0C0C0',solid_capstyle='round')
    for i in range(20):
        plt.plot([0.25*i-2,0.25*i-2],[-0.01,0.01],linewidth=1,color='#C0C0C0',solid_capstyle='round')
    for i in range(20):
        plt.plot([-0.01,0.01],[0.25*i-2,0.25*i-2],linewidth=1,color='#C0C0C0',solid_capstyle='round')
    plt.xlim(lim_eje_x)
    plt.ylim(lim_eje_y)


    # Plotear la curva
    plt.plot(x,y,linewidth=2,color='#FAC205',solid_capstyle='round')
    
    #Plotear punto que se desplaza
    plt.plot(x[f],y[f],'.',color='#EF4026',zorder=3.5*2)

    # Vector tangente
    factor = 2 # Para ajustar el tamaño en el dibujo
    if(True):
        plt.arrow(x[f], y[f], dx[f]/factor, dy[f]/factor, width=0.01, color='#EF4026',length_includes_head=True,zorder=4)
   
    # Vector normal
    if(True):
        cu = curvature[f]
        plt.arrow(x[f], y[f], normal[f,0]/cu, normal[f,1]/cu, width=0.01, color='#04B8C2',length_includes_head=True,zorder=3)
        plt.plot(x[f],y[f],'.',color='#04B8C2',zorder=2)
        
    # Círculo tangente
    if(True):
        circle=plt.Circle((x[f]+normal[f,0]/cu, y[f]+normal[f,1]/cu), 1/cu, color='#04B8C2',fill=False,zorder=3)
        plt.gca().add_patch(circle)
    

    plt.savefig("Curva"+str(f-2)+".png",transparent=True,dpi=500)
    #plt.show()
    plt.close()

Hay tres elementos opcionales para añadir al gráfico: el vector tangente, el vector normal y la circunferencia osculatriz. Aparecerán en el dibujo los que tengan un True en su sección.

También hay una variable que he llamado factor. Sirve para escalar los vectores y que tengan un tamaño visualmente adecuado. A veces salen muy grandes, a veces muy pequeños. Juega con él hasta encontrar el tamaño que más bonito quede.

En este caso, he decidido que el vector normal tenga la longitud del radio de la circunferencia osculatriz, de forma que el vector normal siempre apunta al centro de la circunferencia. La relación que hay entre el radio de la circunferencia y la curvatura de la curva en un punto es sencilla:

R(t) = \frac{1}{|k(t)|}

Curva con significado físico

Si imaginamos la curva como la trayectoria que sigue una hormiga sobre un plano, muchos de los elementos que han parecido en las secciones anteriores se pueden interpretar de forma física. Por ejemplo, el vector tangente nos dice la velocidad de la hormiga: tanto la dirección como el módulo.

Si atendemos a la cinemática del movimiento, hay dos parámetros muy útiles para describir la trayectoria de un cuerpo: la aceleración normal y la tangente. Aunque ambas actúan a la vez, podemos decir que la aceleración tangente es la que nos indica cómo cambia la velocidad del cuerpo según avanza por la curva. A veces va más rápido, a veces frena un poco.

La aceleración normal es la que curva el movimiento. Hace el papel de una fuerza que tuerce el cuerpo hacia la derecha o la izquierda, lo desvía de seguir un camino recto. Cuanto más grande sea esta aceleración, más se curvará la trayectoria del cuerpo. Podemos dibujar estas dos aceleraciones con el siguiente código.

# Aceleraciones
lim_eje_x = [-2.05, 2.05]
lim_eje_y = [-2.05, 2.05]

for f in range(2, n-2):
    #Borrar el formato de matplotlib y ajustes
    plt.gca().set_aspect('equal')
    plt.axis('off')
    
    # Plotear los ejes coordeandos
    plt.plot([-2,2],[0,0],linewidth=1,color='#C0C0C0',solid_capstyle='round')
    plt.plot([0,0],[-2,2],linewidth=1,color='#C0C0C0',solid_capstyle='round')
    for i in range(20):
        plt.plot([0.25*i-2,0.25*i-2],[-0.01,0.01],linewidth=1,color='#C0C0C0',solid_capstyle='round')
    for i in range(20):
        plt.plot([-0.01,0.01],[0.25*i-2,0.25*i-2],linewidth=1,color='#C0C0C0',solid_capstyle='round')
    plt.xlim(lim_eje_x)
    plt.ylim(lim_eje_y)
    

    plt.plot(x,y,linewidth=2,color='#FAC205',solid_capstyle='round')

    # Vector tangente
    factor = 1
    plt.arrow(x[f], y[f], t_comp[f,0]/factor, t_comp[f,1]/factor, width=0.01, color='#EF4026',length_includes_head=True,zorder=4)
    plt.plot(x[f],y[f],'.',color='#EF4026',zorder=3.5)
    
    # Vector normal
    plt.arrow(x[f], y[f], n_comp[f,0]/factor,  n_comp[f,1]/factor, width=0.01, color='#04B8C2',length_includes_head=True,zorder=3)
    plt.plot(x[f],y[f],'.',color='#04B8C2',zorder=2)


    plt.savefig("Aceleraciones"+str(f)+".png",transparent=True,dpi=500)
    #plt.show()
    plt.close()