impetu 0.2.0

Motor de física clásica sin dependencias: N cuerpos simpléctico, Dormand-Prince adaptativo y cuerpo rígido con cuaterniones. Energía conservada a nivel de redondeo.

ímpetu

Motor de física clásica zero-dependency. Solo la stdlib de Python (math), sin NumPy ni nada más.

No es otro wrapper de Euler. ímpetu trae integradores simplécticos que conservan la energía de sistemas hamiltonianos durante millones de pasos —donde un RK4 ingenuo hace que las órbitas espiraleen— más un Dormand–Prince RK5(4) adaptativo con control de error para lo no conservativo, y una mecánica de cuerpo rígido con cuaterniones que integra las ecuaciones de Euler y reproduce el efecto Dzhanibekov (teorema de la raqueta de tenis) conservando energía y momento angular a nivel de redondeo.

pip install impetu

Por qué

Un integrador de Euler o incluso RK4 sobre un problema de dos cuerpos acumula error de energía sin cota: la órbita se abre o colapsa. Los integradores simplécticos (Verlet, Yoshida-4) conservan un hamiltoniano sombra, así que la energía queda acotada para siempre, no crece.

La evidencia, tras 200 órbitas keplerianas (dt=0.01):

Integrador	orden	|ΔE/E₀| máx
euler	1	7.6e-01 (diverge)
rk4	4	3.5e-09
verlet	2	2.5e-09
yoshida4	4	8.6e-14 (precisión de máquina)

Yoshida-4 es de 4º orden y simpléctico: la precisión de RK4 con energía conservada a nivel de redondeo.

Uso

import impetu as im

# Sistema Sol–Tierra en unidades naturales (G=1)
sol    = im.Body(1.0,  im.Vec3(0, 0, 0))
tierra = im.Body(3e-6, im.Vec3(1, 0, 0), im.Vec3(0, 1, 0))

sis  = im.NBody([sol, tierra], G=1.0)
traj = sis.run(dt=0.01, steps=628, method="yoshida4")

print(traj.max_energy_drift())   # ~1e-14

Álgebra vectorial

Vec3 es un dataclass inmutable con slots (hashable, ligero):

v = im.Vec3(3, 4, 0)
v.norm                      # 5.0
v.normalized()              # Vec3(0.6, 0.8, 0)
im.X_AXIS.cross(im.Y_AXIS)  # Vec3(0, 0, 1)
v.rotate(im.Z_AXIS, 3.14159)  # rotación de Rodrigues

Integradores de propósito general

Para cualquier EDO dy/dt = f(t, y):

from impetu import DormandPrince

dp = DormandPrince(rtol=1e-9, atol=1e-12)
ts, ys = dp.integrate(f, y0=(1.0, 0.0), t0=0.0, t1=10.0,
                      t_eval=[2.5, 5.0, 7.5])   # aterriza exacto en esos t

Pasos de un solo tick disponibles sueltos: euler_step, rk4_step, velocity_verlet_step, yoshida4_step.

Constantes físicas (SI, CODATA 2018 / IAU)

from impetu import constants as c

c.c        # velocidad de la luz, 2.998e8
c.G        # gravitación, 6.674e-11
c.hbar     # Planck reducida
c.M_sun    # masa solar
c.au       # unidad astronómica

Fórmulas cerradas

Orbitales (orbital_period, circular_velocity, escape_velocity, vis_viva, specific_orbital_energy, semimajor_axis) y cinemática de proyectiles (projectile_range, projectile_time_of_flight, projectile_max_height).

im.escape_velocity(mu=c.G * c.M_earth, r=c.R_earth)   # ~11.2 km/s

Fuerzas y colisiones

gravitational_force, spring_force, drag_force, elastic_collision_1d, resolve_collision (impulso con coeficiente de restitución).

Cuerpo rígido y cuaterniones

Quat es un cuaternión unitario (escalar primero, w,x,y,z) inmutable con slots. Convención activa cuerpo→mundo: v_mundo = q.rotate(v_cuerpo).

q = im.Quat.from_axis_angle(im.Z_AXIS, math.pi / 2)
q.rotate(im.X_AXIS)          # Vec3(0, 1, 0)
q.to_matrix()                # matriz 3×3 ortonormal
a.slerp(b, 0.5)              # interpolación esférica sobre el arco corto
im.Quat.between(u, v)        # rotación mínima que lleva u a v

RigidBody integra la traslación (Euler semi-implícito, simpléctico) y la rotación resolviendo las ecuaciones de Euler I·ω̇ = τ − ω×(Iω) con RK2, avanzando la orientación por el mapa exponencial del cuaternión. La inercia son los tres momentos principales (marco cuerpo); hay helpers para sólidos comunes.

El caso estrella es la inestabilidad del eje intermedio (efecto Dzhanibekov): un cuerpo con I₁<I₂<I₃ que gira en torno al eje intermedio tumba caóticamente, aunque energía y momento angular se conserven exactos.

rb = im.RigidBody(
    mass=1.0,
    inertia=im.Vec3(1.0, 2.0, 3.0),          # gira sobre el eje intermedio (y)
    angular_velocity=im.Vec3(0.2, 5.0, 0.1),
)
for _ in range(200_000):
    rb.step(1e-4)
# ΔE/E₀ ≈ 2e-9 y ΔL/L₀ ≈ 1e-9, pero el eje de giro tumba: la energía
# "fuga" del eje intermedio a los otros dos y vuelve, una y otra vez.

Helpers de inercia: inertia_solid_sphere, inertia_hollow_sphere, inertia_solid_box, inertia_solid_cylinder.

Métodos de integración de NBody

method	orden	simpléctico	uso
yoshida4	4	sí	por defecto; largo plazo
verlet	2	sí	rápido, estable
rk4	4	no	precisión a corto plazo
euler	1	no	didáctico / comparación

La gravitación usa suavizado de Plummer y fuerzas por pares simétricas, así que el momento lineal se conserva exacto salvo redondeo con cualquier integrador.

Diseño

• Cero dependencias. Portable, auditable, sin cadena de suministro que romperse.
• Tipado completo (py.typed), from __future__ import annotations.
• Precisión numérica con math.fsum donde importa.
• Python ≥ 3.11.

Verificación

python verify.py corre una batería de oráculos con solo la stdlib: valores de referencia, órdenes de convergencia empíricos, conservación de energía a 200 órbitas, cierre metamórfico de órbitas, invariantes de N cuerpos y colisiones, y el trompo asimétrico libre de par (conservación de E y L bajo el régimen inestable de Dzhanibekov). La suite completa (pytest) añade propiedades con Hypothesis y las anclas de regresión de las auditorías.

Benchmarks y precisión (medidos)

Campañas con preregistro, oráculos y scripts reproducibles en benchmarks/ (Ryzen 5 5600G, Python 3.13, scipy 1.15, numpy 2.3).

Contra la competencia (RESULTADOS.md):

Benchmark	Resultado
Mismo tableau DP 5(4) vs scipy.solve_ivp(RK45)	precisión idéntica (nfev 1.00×), 1,8× más rápido (sistemas pequeños)
Conservación de energía a igual tiempo de CPU (10³ órbitas)	583× mejor que RK45, 15,8× mejor que DOP853
Throughput N cuerpos vs NumPy vectorizado (mismo algoritmo)	gana hasta N≈4; pierde ~10× para N≥64

Precisión (PRECISION.md, oráculos exactos con Fraction/Decimal):

• Álgebra en el redondeo correcto: norm a 0 ulps, normalized ≤2 ulps, to_matrix ortonormal a 9e-16.
• Órdenes de convergencia medidos 1.00 / 2.00 / 3.99 / 4.00; suelos de redondeo en ~1e-15.
• DormandPrince: error global ≈ 100 × rtol, estable de 1e-8 a 1e-13.
• Largo plazo: energía acotada (~1e-10 tras 10⁶ pasos); la posición solo deriva en fase, lineal y predecible (5,8e-8/órbita a dt=0.01) — 99,94 % de exactitud tras 10.000 órbitas con la energía intacta a 9 cifras.

Licencia

MIT © 2026 Raul Cruz Acosta

Autores: Esraderey · omri