Simular una órbita con el método de Euler

La segunda ley de Newton como ecuación diferencial resuelta a pasos: el método de Euler, por qué hace crecer la órbita (la energía aumenta como h²) y cómo lo arregla Euler–Cromer.

⏱ 165 min derivada como razón de cambioecuaciones diferencialesvectores en el planoaproximación linealsucesiones definidas por recurrenciaorden de un método numéricoconservación de la energía revision

Las leyes de Kepler describen cómo se mueven los planetas. Newton explicó por qué: una fuerza produce una aceleración. Pero pasar de «conozco la aceleración» a «conozco la órbita» exige resolver una ecuación diferencial, y para la gravedad no hay una fórmula elemental. En esta lección vamos a hacerlo como lo hacen los computadores: a pasitos. Y vamos a descubrir que la forma más ingenua de hacerlo tiene una trampa que un físico detectó, y arregló, en 1981.

1. Repaso: la ley de Newton es una ecuación diferencial

Ilustración: Repaso: la ley de Newton es una ecuación diferencial

En la lección anterior vimos que la velocidad es la derivada de la posición, $\vec v = \vec r\,'$ , y la aceleración la derivada de la velocidad, $\vec a = \vec v\,'$ . En el plano usamos vectores: $\vec r = (x, y)$ , $\vec v = (v_x, v_y)$ , y cada componente se deriva por separado.

Un planeta de masa $m$ alrededor de una estrella de masa $M$ siente la fuerza de gravedad, dirigida hacia la estrella y de tamaño $GMm/r^2$ , donde $r = |\vec r|$ es la distancia. Con $\vec F = m\vec a$ , la masa del planeta se cancela:

\vec a = -\frac{GM}{r^3}\,\vec r .

(El factor $\vec r/r$ es un vector de tamaño 1 que apunta de la estrella al planeta; dividir por $r^2$ más y poner el signo menos lo convierte en la aceleración que apunta a la estrella.) Como $\vec a$ es la segunda derivada de $\vec r$ , esto es una ecuación diferencial: relaciona la función desconocida $\vec r(t)$ con sus derivadas.

Para no cargar con números enormes usaremos unidades donde $GM = 1$ . Entonces una órbita circular de radio $r = 1$ tiene velocidad $v = 1$ (porque $v^2/r = GM/r^2$ ) y periodo $T = 2\pi \approx 6.28$ .

2. El método de Euler: seguir la tangente

Ilustración: El método de Euler: seguir la tangente

Si $h$ es pequeño, la definición de derivada nos dice que

x(t + h) \approx x(t) + h\,x'(t).

Es seguir la recta tangente durante un tiempito $h$ . Aplicado a la posición y a la velocidad, obtenemos el método de Euler: a partir del estado en el paso $n$ , calculamos el paso $n+1$ con

\vec r_{n+1} = \vec r_n + h\,\vec v_n, \qquad \vec v_{n+1} = \vec v_n + h\,\vec a_n, \qquad\text{con}\quad \vec a_n = -\frac{\vec r_n}{r_n^3}.

Es una sucesión definida por recurrencia: cada término se calcula con el anterior.

Ejemplo resuelto: tres pasos a mano

Empezamos en la órbita circular: $\vec r_0 = (1, 0)$ , $\vec v_0 = (0, 1)$ , y usamos $h = 0.1$ .

Paso 1. $r_0 = 1$ , así que $\vec a_0 = -(1, 0)/1^3 = (-1, 0)$ .

\vec r_1 = (1, 0) + 0.1\,(0, 1) = (1,\ 0.1), \qquad \vec v_1 = (0, 1) + 0.1\,(-1, 0) = (-0.1,\ 1).

Paso 2. $r_1 = \sqrt{1^2 + 0.1^2} = \sqrt{1.01} \approx 1.0050$ , de modo que $\vec a_1 = -(1, 0.1)/1.0050^3 \approx (-0.9852,\ -0.0985)$ .

\vec r_2 = (1, 0.1) + 0.1\,(-0.1, 1) = (0.99,\ 0.2), \qquad \vec v_2 \approx (-0.1, 1) + 0.1\,(-0.9852, -0.0985) = (-0.1985,\ 0.9901).

Paso 3. Repite tú el procedimiento. Resumen (redondeado a 4 decimales):

$n$	$x_n$	$y_n$	$v_{x,n}$	$v_{y,n}$	$r_n$	$E_n$
0	1.0000	0.0000	0.0000	1.0000	1.0000	−0.5000
1	1.0000	0.1000	−0.1000	1.0000	1.0050	−0.4900
2	0.9900	0.2000	−0.1985	0.9901	1.0100	−0.4802
3	0.9701	0.2990	−0.2946	0.9707	1.0152	−0.4705
4	0.9407	0.3961	−0.3873	0.9422	1.0207	−0.4609

La última columna es la energía por unidad de masa,

E = \frac{v^2}{2} - \frac{GM}{r},

que en la órbita real se conserva (vale $-1/2$ todo el tiempo). Pero en la tabla el radio crece y la energía sube en cada paso. Tras una vuelta completa (63 pasos con $h = 0.1$ ), el planeta simulado ya está a $r \approx 2.1$ : ¡la órbita se abrió en espiral!

3. ¿Por qué Euler hace crecer la órbita?

Ilustración: ¿Por qué Euler hace crecer la órbita?

Mira el primer paso. En una órbita circular la velocidad es perpendicular al radio. Euler mueve el planeta en línea recta por la tangente, así que el nuevo radio es la hipotenusa de un triángulo rectángulo de catetos $1$ y $h$ :

r_1^2 = 1 + h^2 > 1 .

Y la aceleración es perpendicular a la velocidad, así que también por Pitágoras $v_1^2 = 1 + h^2$ . Entonces

E_1 = \frac{1 + h^2}{2} - \frac{1}{\sqrt{1 + h^2}} \approx -\frac{1}{2} + h^2 ,

usando $1/\sqrt{1 + u} \approx 1 - u/2$ para $u$ pequeño. Cada paso añade energía (del orden de $h^2$ ). Una tangente siempre queda por fuera de un círculo, y Euler siempre sigue tangentes.

¿Y si hacemos $h$ más pequeño?

Cada paso añade $\sim h^2$ , pero una vuelta necesita $2\pi/h$ pasos: el error por vuelta es $\sim h^2 \times 2\pi/h = 2\pi h$ , proporcional a $h$ . Se dice que Euler es un método de orden 1: dividir $h$ entre 10 solo divide el error entre 10.

$h$	Pasos por vuelta	$r$ tras una vuelta (Euler)	Error
0.1	63	2.09	109 %
0.01	628	1.121	12 %
0.001	6283	1.0125	1.25 %

Para bajar el error a una parte en un millón harían falta mil millones de pasos por vuelta. Hace falta una idea mejor, no un computador más rápido.

4. Una corrección pequeña: Euler–Cromer

Ilustración: Una corrección pequeña: Euler–Cromer

Alan Cromer (1981) notó que basta con cambiar el orden: primero actualiza la velocidad y luego usa la velocidad nueva para mover el planeta:

\vec v_{n+1} = \vec v_n + h\,\vec a_n, \qquad \vec r_{n+1} = \vec r_n + h\,\vec v_{n+1}.

El primer paso da $\vec v_1 = (-0.1,\ 1)$ y $\vec r_1 = (0.99,\ 0.1)$ : el planeta ya «gira hacia adentro» un poco.

$n$	$x_n$	$y_n$	$v_{x,n}$	$v_{y,n}$	$r_n$	$E_n$
0	1.0000	0.0000	0.0000	1.0000	1.0000	−0.5000
1	0.9900	0.1000	−0.1000	1.0000	0.9950	−0.5000
2	0.9700	0.1990	−0.2005	0.9898	0.9902	−0.4999
3	0.9399	0.2959	−0.3004	0.9694	0.9854	−0.4999
4	0.9000	0.3898	−0.3986	0.9384	0.9808	−0.4998

La energía oscila un poquito pero no se desboca: tras una vuelta con $h = 0.1$ , el radio es $1.0006$ (error del 0.06 %, frente al 109 % de Euler con el mismo costo); tras diez vueltas, $1.006$ ; y con $h = 0.01$ el error tras una vuelta es de 2 partes en 100 000. Euler–Cromer pertenece a una familia especial, los integradores simplécticos, que respetan la geometría del movimiento: su error de energía queda acotado en lugar de acumularse. En la universidad se estudia por qué (el módulo metodos-numericos lo retoma con el método leapfrog, que es Euler–Cromer «simetrizado»).

5. Compara con la solución exacta

Ilustración: Compara con la solución exacta

El simulador siguiente no usa pasos: resuelve la ecuación de Kepler, que es la solución exacta de la ecuación diferencial para dos cuerpos. Ponlo con $e = 0$ (órbita circular) y compara: el planeta vuelve siempre al mismo punto y $r$ no cambia. Ese es el comportamiento que un buen método numérico debe reproducir. Con $e = 0.9$ , observa lo rápido que cambia la velocidad cerca del perihelio: allí cualquier método con paso fijo sufre.

Semieje mayor a = 1.00 UAExcentricidad e = 0.500Velocidad = 0.30 años/s

Periodo T = 1.000 años · t = 0.000 años · r = 0.000 UA · v = 0.00 km/s

Para curiosos: un simulador con dos integradores

Este segundo simulador (del nivel universitario) integra numéricamente, como tú en la tabla, y grafica en escala logarítmica el error de energía. Elige el escenario «2 cuerpos, e = 0.8», el integrador «Ambos superpuestos» y mira la gráfica: el error del método simpléctico (leapfrog) oscila sin crecer, mientras que el de RK4, un método más preciso pero no simpléctico, crece poco a poco. Después sube el paso $h$ y observa cómo se agrandan los errores.

EscenarioIntegradorPaso h = 1.0e-2Velocidad = 1.5 unidades de t / s Perturbar x₁ en δ₀ = 10⁻⁶ (copia vecina)

tres masas iguales, T ≈ 6.3259 · t = 0.00

Actividad: la órbita en una hoja de cálculo

Materiales: computador o celular con hoja de cálculo (LibreOffice, Excel, Google Sheets).

Crea las columnas t, x, y, vx, vy, r, E. En la fila 2 pon el estado inicial: 0, 1, 0, 0, 1, con r = =RAIZ(B2^2+C2^2) y E = =(D2^2+E2^2)/2-1/F2.
En la fila 3 escribe las fórmulas de Euler con $h = 0.1$ : x = =B2+0.1*D2, y = =C2+0.1*E2, vx = =D2+0.1*(-B2/F2^3), vy = =E2+0.1*(-C2/F2^3), y copia r y E. Arrastra hacia abajo 200 filas.
Haz un gráfico de dispersión de y contra x: verás la espiral. Grafica también E contra t.
Copia la hoja y cambia solo las fórmulas de posición para que usen la velocidad de la misma fila (x = =B2+0.1*D3, y = =C2+0.1*E3): es Euler–Cromer. Compara las dos gráficas.
Cambia $h$ a 0.05 y a 0.01 (necesitarás más filas) y anota $r$ tras una vuelta ( $t \approx 6.28$ ) en cada caso.

Pregunta final: ¿cómo cambia el error de Euler al dividir $h$ entre 10? ¿Y el de Euler–Cromer? ¿Cuál de los dos usarías para simular el Sistema Solar durante mil años?

Ejercicios

Básico

Con $GM = 1$ , ¿cuál es la aceleración de un planeta en $\vec r = (2, 0)$ ? ¿Y en $\vec r = (0, -0.5)$ ?

Ver solución

$\vec a = -(2, 0)/8 = (-0.25, 0)$ ; $\vec a = -(0, -0.5)/0.125 = (0, 4)$ . Siempre apunta a la estrella y es 4 veces mayor si la distancia se reduce a la mitad.

Verifica el paso 1 de la tabla de Euler y calcula $E_1$ .

Ver solución

$\vec r_1 = (1, 0.1)$ , $\vec v_1 = (-0.1, 1)$ ; $E_1 = 1.01/2 - 1/\sqrt{1.01} = 0.505 - 0.99504 \approx -0.4900$ .

Haz el primer paso de Euler–Cromer con $h = 0.1$ y comprueba que $E_1 \approx -0.5$ .

Ver solución

$\vec v_1 = (-0.1, 1)$ , $\vec r_1 = (0.99, 0.1)$ , $r_1 = \sqrt{0.9901} \approx 0.99504$ ; $E_1 = 0.505 - 1/0.99504 = 0.505 - 1.00498 \approx -0.49998$ .

Medio

Completa a mano el paso 3 de Euler (partiendo de la fila 2) y compáralo con la tabla.

Ver solución

$r_2 \approx 1.0100$ , $\vec a_2 \approx -(0.99, 0.2)/1.0100^3 \approx (-0.9609, -0.1941)$ ; $\vec r_3 = (0.99 - 0.01985,\ 0.2 + 0.09901) \approx (0.9701, 0.2990)$ ; $\vec v_3 \approx (-0.1985 - 0.0961,\ 0.9901 - 0.0194) = (-0.2946, 0.9707)$ . ✓

Haz el primer paso de Euler con $h = 0.05$ . ¿Cuánto vale $r_1$ ? ¿Cuánto aumentó la energía? Compara con la predicción $\Delta E \approx h^2$ .

Ver solución

$r_1 = \sqrt{1.0025} \approx 1.00125$ ; $E_1 = 0.50125 - 0.99875 \approx -0.4975$ : aumentó $\approx 0.0025 = 0.05^2$ . ✓

Usa la tabla de la sección 3 para estimar el error de Euler tras una vuelta con $h = 0.0001$ . ¿Cuántos pasos serían? ¿Y cuántos necesitaría Euler–Cromer con $h = 0.1$ para el mismo error, si su error oscila sin crecer?

Ver solución

Orden 1: al dividir $h$ entre 10 el error se divide entre 10: $\approx 0.125\,\%$ , con 62 832 pasos. Euler–Cromer con $h = 0.1$ ya tiene un error del 0.06 % con 63 pasos, mil veces menos trabajo.

Reto

¿Por qué no basta con hacer $h$ diez veces más pequeño para que Euler funcione bien durante 1000 órbitas? Piensa en cuántos pasos hay y cuánta energía añade cada uno.

Ver solución

Con $h$ diez veces menor, cada paso añade unas 100 veces menos energía, pero hay 10 veces más pasos: el error por vuelta baja solo unas 10 veces (orden 1). Para 1000 órbitas el error acumulado sigue creciendo sin límite; hace falta un método mejor, no solo un paso menor.

Un estudiante simula la Tierra con Euler durante 100 años y concluye que «la Tierra se aleja del Sol». ¿Qué le responderías? ¿Cómo podría comprobar si el resultado es físico o numérico?

Ver solución

Que puede ser un error numérico: Euler siempre añade energía. Puede repetir con $h$ más pequeño (si el resultado cambia, es numérico), vigilar la energía (debe conservarse) o usar Euler–Cromer. Un resultado que depende del paso no es física.

En el método de Euler, ¿qué papel cumple la aproximación $x(t + h) \approx x(t) + h\,x'(t)$ ? Relaciónala con la recta tangente y explica en qué casos el error es grande.

Ver solución

Es seguir la recta tangente durante un tiempo $h$ . El error es grande si $h$ es grande comparado con el tiempo en que cambia la velocidad: por ejemplo, cerca del Sol en una órbita muy excéntrica, donde la curva se dobla rápido (míralo en el simulador con $e = 0.9$ ).

Recuerda

Newton: $\vec a = -GM\,\vec r/r^3$ es una ecuación diferencial; con $GM = 1$ , la órbita circular tiene $r = v = 1$ y $T = 2\pi$ .
Euler: $\vec r_{n+1} = \vec r_n + h\vec v_n$ , $\vec v_{n+1} = \vec v_n + h\vec a_n$ . Sigue tangentes, así que la órbita crece: $\Delta E \approx h^2$ por paso.
Es de orden 1: el error por vuelta es proporcional a $h$ .
Euler–Cromer: primero $\vec v$ , luego $\vec r$ con la $\vec v$ nueva. Mismo costo, energía acotada.
La conservación de la energía es el termómetro de una simulación: si $E$ deriva, el resultado es numérico, no físico.

🍎 Guía docente

Tiempo: 3 sesiones de 55 min
Estándar MEN: Matemáticas 10.º a 11.º, pensamiento variacional: «Utilizo las técnicas de aproximación en procesos infinitos numéricos» e «Interpreto la noción de derivada como razón de cambio y como valor de la pendiente de la tangente a una curva y desarrollo métodos para hallar las derivadas de algunas funciones básicas en contextos matemáticos y no matemáticos»; pensamiento métrico: «Justifico resultados obtenidos mediante procesos de aproximación sucesiva, rangos de variación y límites en situaciones de medición». Ciencias Naturales 10.º a 11.º: «Modelo matemáticamente el movimiento de objetos cotidianos a partir de las fuerzas que actúan sobre ellos», «Establezco relaciones entre las diferentes fuerzas que actúan sobre los cuerpos en reposo o en movimiento rectilíneo uniforme y establezco condiciones para conservar la energía mecánica», «Propongo modelos para predecir los resultados de mis experimentos y simulaciones» y «Relaciono la información recopilada con los datos de mis experimentos y simulaciones».
Evidencia de aprendizaje: Ejecuta a mano tres pasos de Euler y uno de Euler–Cromer con sus energías, explica con Pitágoras por qué Euler añade ≈ h² de energía por paso, programa ambos métodos en una hoja de cálculo y mide cómo cambia el error tras una vuelta al dividir h entre 10.

Materiales

Calculadora científica
Computador por pareja con hoja de cálculo
Fotocopia de las dos tablas de pasos con casillas en blanco para completar
Computador o tableta con los simuladores

Secuencia sugerida

Sesión 1. La ecuación y los pasos (55 min). Repaso de 10 min (lección anterior). Deducción de $\vec a = -\vec r/r^3$ y unidades $GM = 1$ (10 min). Los tres pasos a mano repartidos entre grupos (cada grupo hace un paso y lo pasa al siguiente): se hace visible que cada paso depende del anterior (25 min). Cierre: ejercicios 1 y 2.

Sesión 2. La trampa y el arreglo (55 min). Pitágoras y $\Delta E \approx h^2$ (15 min); tabla de orden 1 (10 min); Euler–Cromer con el primer paso a mano y ejercicio 3 (15 min); comparación con la solución exacta y, si hay tiempo, el simulador de dos integradores (15 min).

Sesión 3. Hoja de cálculo (55 min). Actividad por parejas (40 min): espiral de Euler, Euler–Cromer y medición del error con tres valores de $h$ . Mini-prueba (10 min). Retos 7–9 como tarea. Conecta con el proyecto del grado: programar la órbita de la Tierra con Euler.

Errores frecuentes

Usar r en lugar de r³ en el denominador (olvidar que r/r ya es unitario). Verificar con el ejercicio 1.
En Euler–Cromer, usar la velocidad vieja al actualizar la posición. En la hoja de cálculo, señalar la referencia D3 frente a D2.
Interpretar la espiral como un efecto físico («los planetas se alejan»). Ejercicio 8: la energía debe conservarse.
Creer que un paso 10 veces menor arregla todo. Tabla de la sección 3: orden 1.
Confundir la solución exacta del simulador de Kepler con «otra simulación a pasos».

1. Repaso: la ley de Newton es una ecuación diferencial