NumPy para estados, trayectorias y mallas

Forma de los arreglos y convención del modelo

NumPy es una biblioteca pública para cálculo con arreglos. La instrucción import numpy as np la carga bajo el alias convencional np. Primero se construye un vector pequeño y se inspecciona:

import numpy as np

state = np.array([1.0, 0.0], dtype=np.float64)
print(state)
print(state.shape)   # (2,)
print(state.dtype)   # float64

np.array(datos, dtype=tipo) recibe los datos como primer argumento y el tipo como argumento nombrado. shape y dtype son atributos del arreglo, por eso se consultan sin paréntesis. En cambio, state.copy() es un método y sí se llama con paréntesis.

NumPy no adivina qué representa cada eje. Una trayectoria con $N$ tiempos y $d$ variables puede almacenarse como $(N,d)$ o como $(d,N)$; ambas convenciones son válidas, pero mezclarlas produce resultados numéricamente plausibles y conceptualmente falsos. En los ejemplos que siguen, cada fila será un tiempo y cada columna una variable:

times = np.linspace(0.0, 40.0, 8001, dtype=np.float64)
states = np.empty((times.size, 3), dtype=np.float64)
states[0] = (1.0, 1.0, 1.0)

assert states.shape == (8001, 3)
assert times.shape == (8001,)

La guía inicial oficial de NumPy explica dimensiones, forma, tamaño, tipo y operaciones por ejes. Conviene predecir la forma de cada expresión antes de ejecutarla. El código debe poder responder: ¿qué eje enumera tiempos?, ¿cuál variables?, ¿qué unidad tiene cada columna?

np.linspace(inicio, fin, cantidad) crea muestras equiespaciadas e incluye ambos extremos. np.empty(forma) reserva memoria sin fijar sus valores; por eso debe llenarse antes de leerse. Estas funciones proceden de NumPy y tienen un contrato documentado. La función validate_trajectory del apartado siguiente será código propio y añadirá las comprobaciones científicas de este curso.

Seleccionar states[:, 0] produce una serie de longitud $N$. Seleccionar states[2000:] descarta un transitorio por índice, pero solo representa un tiempo físico si se conoce el muestreo. Es más seguro calcular una máscara:

retained = states[times >= 10.0]

Así, cambiar el paso o el número de muestras no cambia inadvertidamente el transitorio físico.

Práctica interactiva. En una celda de Jupyter o del navegador, crea los arreglos [1, 2, 3] y [[1, 2, 3], [4, 5, 6]]. Predice ndim, shape, size y el resultado de sumar 10 antes de ejecutar. Transpón la matriz y explica por qué cambian los ejes, pero no el número de elementos. El ejercicio usa solo NumPy y no requiere archivos ni hardware externo.

Tipos, finitud y escala

El tipo determina precisión y memoria. Una trayectoria de diez millones de estados tridimensionales ocupa aproximadamente $240$ MB en precisión doble y $120$ MB en precisión simple. Reducir memoria no es gratis: puede cambiar el horizonte antes de que el redondeo domine una dinámica sensible.

Se valida antes de analizar:

def validate_trajectory(times, states, dimension):
    times = np.asarray(times, dtype=np.float64)
    states = np.asarray(states, dtype=np.float64)
    if times.ndim != 1:
        raise ValueError("times debe ser un vector")
    if states.shape != (times.size, dimension):
        raise ValueError("states debe tener forma (N, dimension)")
    if not np.isfinite(times).all() or not np.isfinite(states).all():
        raise FloatingPointError("Hay NaN o infinito")
    if not np.all(np.diff(times) > 0):
        raise ValueError("El tiempo no es estrictamente creciente")
    return times, states

La conversión no repara el experimento; normaliza el contrato. Un infinito puede señalar paso inadecuado, parámetros erróneos o escape real. Se registra la causa antes de eliminar la fila.

Vistas, copias y cambios silenciosos

Muchas rebanadas son vistas de la misma memoria:

x_view = states[:, 0]
x_copy = states[:, 0].copy()
x_view[0] = 99.0

La primera asignación también cambia states[0, 0]; la copia permanece independiente. Esta propiedad ahorra memoria, pero puede corromper datos primarios si una función “limpia” una serie en sitio. Una función de análisis debería documentar si modifica la entrada. Por defecto, esta ruta conserva los datos y devuelve un resultado nuevo.

El indexado booleano suele crear copia:

mask = times >= 10.0
tail = states[mask]

No se memoriza una tabla aislada de reglas. Se usa np.shares_memory cuando la independencia sea importante y se añade una prueba que asegure que el análisis no altera la trayectoria.

Broadcasting como álgebra sobre estados

El broadcasting oficial alinea dimensiones desde la derecha. Si states tiene forma $(N,3)$ y equilibrium forma $(3,)$, la resta

equilibrium = np.array([0.0, 0.0, 0.0])
displacements = states - equilibrium
distances = np.linalg.norm(displacements, axis=1)

trata el equilibrio como una fila que se aplica a todos los tiempos. El argumento axis=1 calcula una norma por estado; usar axis=0 produciría una norma por variable. El resultado puede tener números razonables en ambos casos, por lo que la prueba debe revisar también la forma esperada $(N,)$.

Para comparar $M$ condiciones iniciales con $K$ equilibrios:

delta = initial[:, None, :] - equilibria[None, :, :]
distance = np.linalg.norm(delta, axis=-1)
assert distance.shape == (M, K)
nearest = np.argmin(distance, axis=1)

Los ejes de longitud uno no almacenan copias completas; permiten que NumPy combine formas compatibles. Sin embargo, el arreglo resultado $(M,K,3)$ sí puede ser enorme. Si $M$ y $K$ crecen, se procesa por bloques.

Representación columnar de un oscilador

Considérese una trayectoria del oscilador amortiguado con columnas $(x,v)$ y masa unitaria. La energía mecánica instantánea es $E(t)=v(t)^2/2+kx(t)^2/2$. Si states tiene forma $(N,2)$, el cálculo mantiene una salida por tiempo:

x = states[:, 0]
v = states[:, 1]
k = 4.0
energy = 0.5 * v**2 + 0.5 * k * x**2
assert energy.shape == times.shape

La aserción no es decorativa. Si se hubiera almacenado la trayectoria como $(2,N)$, las mismas rebanadas seleccionarían dos instantes completos y el resultado tendría longitud dos. La fórmula seguiría ejecutándose porque NumPy opera con cualquier forma compatible; solo el significado físico revelaría el error. Por eso una interfaz de datos debe declarar tanto el nombre como la posición de cada eje.

Para un sistema sin amortiguamiento, la variación relativa $|E(t)-E(0)|/E(0)$ ayuda a revisar el integrador. Para el sistema amortiguado no se exige energía constante: se verifica que su tendencia concuerde con $\dot E=-c v^2$. Una pequeña subida aislada puede proceder del error de discretización; una subida sostenida sugiere signo equivocado, paso excesivo o columnas intercambiadas. El arreglo no decide cuál explicación es correcta, pero permite formular indicadores que separen esas posibilidades.

El tiempo merece el mismo cuidado que el estado. La expresión np.diff(times) produce los anchos de $N-1$ intervalos. Si todos coinciden dentro de tolerancia, una FFT puede asociarse con una frecuencia de muestreo única; si el integrador entregó una malla adaptativa, aplicar una FFT directamente sería conceptualmente incorrecto aunque el programa aceptara el arreglo. En ese caso se evalúa la solución densa sobre una malla uniforme o se elige un análisis adecuado para muestreo irregular.

También hay que distinguir normalización numérica de cambio de unidades. Restar la media y dividir por la desviación estándar puede ser útil para una comparación estadística, pero ya no deja coordenadas en metros y metros por segundo. Se conserva el arreglo original y se devuelve, junto con el normalizado, el desplazamiento y la escala usados. Así puede invertirse la transformación y el gráfico puede etiquetar honestamente sus ejes.

Supongamos que se tienen $M$ trayectorias del mismo oscilador con forma $(M,N,2)$. La energía se calcula sin copiar coeficientes sobre el tiempo:

weights = np.array([k, 1.0])
ensemble_energy = 0.5 * np.sum(weights * trajectories**2, axis=-1)
assert ensemble_energy.shape == (M, N)

Aquí el broadcasting aplica dos pesos sobre el último eje y la suma elimina precisamente ese eje. Escribir axis=1 mezclaría tiempos; el tamaño resultante podría parecer correcto para algún caso cuadrado. Una prueba con $M\ne N\ne2$ evita que la coincidencia accidental de dimensiones oculte el defecto.

Mallas y conservación de ejes

Una sección bidimensional de condiciones iniciales se construye así:

x0 = np.linspace(-2.0, 2.0, 161)
v0 = np.linspace(-2.0, 2.0, 121)
X0, V0 = np.meshgrid(x0, v0, indexing="xy")
initial_grid = np.stack((X0, V0), axis=-1)

assert X0.shape == (121, 161)
assert initial_grid.shape == (121, 161, 2)

El primer eje de la malla corresponde a $v_0$ y el segundo a $x_0$. Para un integrador que recibe una lista de estados se aplana solo la parte espacial:

flat_initial = initial_grid.reshape(-1, 2)
labels_flat = classify_many(flat_initial)
labels = labels_flat.reshape(V0.shape)

Se conserva la forma original para reconstruir la imagen. Invertir ejes, transponer etiquetas o usar una extensión gráfica incongruente puede reflejar la cuenca respecto de la diagonal sin producir error de Python.

El experimento de cuencas se inspira en el tratamiento geométrico de Tél y Gruiz Tél, 2006 Tamás Tél y Márton Gruiz (2006) Chaotic Dynamics: An Introduction Based on Classical Mechanics Cambridge University Press Ubicación consultada: §§1.2.2 y 6.8 . Esa obra aporta la pregunta dinámica; NumPy aporta la representación y las operaciones con las que se ejecuta.

Álgebra lineal para lotes de jacobianos

La referencia oficial de álgebra lineal incluye solución de sistemas, autovalores, normas y descomposiciones. Para clasificar un equilibrio de un sistema bidimensional:

jacobian = np.array([[0.0, 1.0], [1.0, -0.3]])
eigenvalues = np.linalg.eigvals(jacobian)

Los autovalores no se ordenan de manera científicamente significativa. Se clasifican por partes reales, no por posición en el arreglo. Para una familia de $M$ jacobianos con forma $(M,d,d)$, varias rutinas aceptan los dos últimos ejes como matrices:

eigenvalues = np.linalg.eigvals(jacobians)
assert eigenvalues.shape == (M, d)
max_real = eigenvalues.real.max(axis=1)

Un valor cercano a cero requiere tolerancia y control del parámetro. No se redondea primero para decidir estabilidad. En matrices no normales, los autovalores tampoco describen todo crecimiento transitorio; la interpretación se limita al criterio calculado.

Para resolver $A\mathbf{x}=\mathbf{b}$ se usa np.linalg.solve(A, b), no se forma la inversa. Se comprueba el residuo $\|A\mathbf{x}-\mathbf{b}\|$ y el condicionamiento. Un residuo pequeño no garantiza error pequeño cuando $A$ está mal condicionada.

Vectorización y legibilidad del algoritmo

Vectorizar elimina bucles de Python cuando todos los elementos siguen la misma operación. No convierte automáticamente un método en correcto ni más legible. Una primera versión explícita puede servir como referencia:

reference = np.array([
    np.linalg.norm(state - equilibrium)
    for state in states
])
vectorized = np.linalg.norm(states - equilibrium, axis=1)
np.testing.assert_allclose(vectorized, reference)

Después de validar, se mide tiempo y memoria. Una expresión compacta que crea varios arreglos temporales gigantes puede ser peor que un bucle por bloques. La elección se documenta con el tamaño real del experimento.

Todavía no se ha ejecutado el programa.

Ensayos de forma, memoria y ejes

1. Para una trayectoria $(5000,4)$, calcula media, mínimo y máximo por variable.
2. Demuestra con np.shares_memory la diferencia entre vista y copia.
3. Predice la forma de tres operaciones de broadcasting y confírmala con aserciones.
4. Resuelve un sistema lineal y compara residuo con condicionamiento.
5. Procesa distancias por bloques y comprueba que coincidan con la versión vectorizada completa.

Cada arreglo del informe debe poder describirse con una frase científica y una tupla de formas. Si ambas descripciones no coinciden, todavía no es seguro usarlo en una figura o en una conclusión cuantitativa.