Solución numérica de las ecuaciones del movimiento

prev.gif (1231 bytes)home.gif (1232 bytes)next.gif (1211 bytes)

Dinámica celeste

Leyes de Kepler
El descubrimiento de
la ley de la gravitación
Fuerza central y
conservativa
Ecuación de la trayectoria
marca.gif (847 bytes)Solución numérica de
 las ecuaciones 
Trayectorias hiperbólicas
Órbita de transferencia
Encuentros espaciales
Trayectoria espiral
Encuentro de una sonda
espacial con Júpiter
Orbitas de la misma
energía
Trayectoria de un 
proyectil (I)
Trayectoria de un 
proyectil (II)
Movimiento relativo
Caída de un satélite en
órbita hacia la Tierra.
Los anillos de un planeta
Movimiento bajo una
fuerza central y una
perturbación
El problema de Euler
Viaje a la Luna
Solución numérica de las ecuaciones del movimiento

Fuerza central y conservativa

Ejemplos

java.gif (886 bytes) Actividades

Código fuente

En esta página, estudiamos el movimiento de un cuerpo de masa m cuando se lanza desde un punto situado en el eje X a una distancia r1 del centro fijo de fuerzas, con velocidades crecientes v1 perpendiculares al radio vector.

Las leyes de Kepler describen el movimiento de los planetas en torno del Sol, sin indagar en las causas que producen tal movimiento.

1.-Los planetas describen órbitas elípticas estando el Sol en uno de sus focos.

2.-El vector posición de cualquier planeta con respecto del Sol, barre áreas iguales de la elipse en tiempos iguales.

3.-Los cuadrados de los periodos de revolución son proporcionales a los cubos de los semiejes de la elipse.

Las leyes de Newton, no solamente explican las leyes de Kepler sino que predicen otras trayectorias para los cuerpos celestes: las parábolas y las hipérbolas. En general, un cuerpo bajo la acción de la fuerza de atracción gravitatoria describirá una trayectoria plana que es una cónica.

Como se ha comentado, las propiedades central y conservativa de la fuerza de atracción entre un cuerpo celeste y el Sol, determinan un sistema de dos ecuaciones diferenciales de primer orden, que cuando se expresan en coordenadas polares, conducen a la ecuación de la trayectoria, una cónica.

El programa interactivo procede de otro modo, calcula las componentes de la aceleración a lo largo del eje X, y a lo largo del eje Y, dando lugar a un sistema de dos ecuaciones diferenciales de segundo orden que se resuelven por procedimientos numéricos

 

Solución numérica de las ecuaciones del movimiento

Supongamos que una partícula de masa m (un planeta) es atraída por un cuerpo masivo de masa M (el Sol) Supondremos que la la influencia de la partícula sobre el cuerpo es despreciable permaneciendo en reposo en el origen.

La partícula está sometida a una fuerza atractiva F cuya dirección es radial y apuntando hacia el centro del Sol.  El módulo de la fuerza viene dado por la ley de la Gravitación Universal

Kepler2.gif (2069 bytes) Siendo r la distancia entre la partícula y el centro de fuerzas, y x e y su posición respecto del sistema de referencia cuyo origen está situado en el Sol.

Las componentes de la fuerza son

Aplicando la segunda ley de Newton, y expresando la aceleración como derivada segunda de la posición, tenemos un sistema de dos ecuaciones diferenciales de segundo orden.

Dadas las condiciones iniciales (posición y velocidad inicial), el sistema de dos ecuaciones diferenciales se puede resolver aplicando el procedimiento numérico de Runge-Kutta.

Escalas

Antes de resolver el sistema de ecuaciones diferenciales por procedimientos numéricos, es conveniente prepararlas para que el ordenador no maneje números excesivamente grandes o pequeños.

Establecemos un sistema de unidades en el que la longitud se mide en unidades astronómicas, la distancia media entre el Sol y la Tierra. L=una UA=1.496·1011 m y el tiempo en unidades de año, P=un año= 365.26 días=3.156·107 s.

En el nuevo sistema de unidades x=Lx’, t=P·t’, la primera ecuación diferencial se escribe

Como L es el semieje mayor de la órbita de la Tierra alrededor del Sol, P es el periodo o tiempo que tarda en dar una vuelta completa, y M es la masa del Sol. Por la tercera ley de Kepler, el término

Volviendo a la notación x e y para la posición y t para el tiempo en el nuevo sistema de unidades. El sistema de ecuaciones diferenciales se escribe

Principio de conservación de la energía

La energía total de la partícula es una constante del movimiento. La energía de la partícula de masa m en el instante inicial t=0 es

Cuando E0<0 la partícula permanece confinada en el espacio que rodea a los dos cuerpos. Cuando E0≥0 la partícula escapa al infinito

La energía de la partícula en el instante t es igual a

En el nuevo sistema de unidades establecido

v=v’·L/P, x=x’·L, y=y’·L, d=d’·L

Volviendo a la notación previa. Definimos una nueva energía e por unidad de masa en este sistema de unidades

El programa interactivo evalúa en cada instante el cociente

que denominaremos tanto por ciento de error. Cuando la energía e difiere de e0 de modo que el cociente es mayor que la unidad el programa interactivo se detiene, la trayectoria calculada puede que se desvíe significativamente de la real, y esto suele ocurrir cuando la partícula está muy cerca del centro de fuerzas

 

Fuerza central y conservativa

En la figura, se muestra un cuerpo que describe una trayectoria elíptica alrededor del centro de fuerzas situado en uno de sus focos. La distancia de máximo acercamiento es r1 y la de máximo alejamiento r2. Las velocidades que lleva el cuerpo en estas dos posiciones extremas son v1 y v2 respectivamente. La constancia del momento angular y de la energía  nos permite relacionar estas cuatro magnitudes.

Se pueden plantear dos problemas :

  • Conocido r1 y r2 calcular v1 y v2.

Se despeja v1 en la primera ecuación y se sustituye en la segunda. Obtenemos

Algunos estudiantes erróneamente, calculan v1 ó v2 a partir de la dinámica del movimiento circular uniforme, pensando que el centro de fuerzas coincide el centro de curvatura de la elipse. Si R es el radio de curvatura de la elipse en la posición más cercana al centro de fuerzas, la ecuación de la dinámica del movimiento circular uniforme se escribe

Como vemos R coincide con el parámetro d, que interviene en la ecuación de la elipse en coordenadas polares.

  • Conocido r1 y v1 calcular r2 y v2.

Conocida la posición r1 y la velocidad v1 en el momento del lanzamiento, aplicando la constancia de la energía y del momento angular calculamos la velocidad v2 y distancia r2 al centro de fuerzas. Después de algunas operaciones, obtenemos la velocidad v2 en función de r1 y v1

       (1)

Velocidad de escape

Se denomina velocidad de escape ve de una partícula que está a una distancia r1 del centro de fuerzas, a la velocidad que hemos de proporcionarle para que llegue al infinito con velocidad nula

Por ejemplo, la velocidad de escape de la superficie de la Tierra es r1=6.37·106 m, M=5.98·1024 kg es ve=11190.7 m/s

Podemos expresar la velocidad v2 en términos de la velocidad de escape ve

Obtenemos la distancia r2 al centro de fuerzas a partir de la constancia del momento angular

Movimiento circular

Si el satélite es lanzado con una velocidad vc tal que

describirá una órbita circular de radio r.

Podemos comprobar en la fórmula (1), que si el satélite sigue una órbita circular entonces v2=v1

En el caso de los satélites artificiales que circundan la Tierra

  • Si la velocidad de lanzamiento v1 es inferior a vc el punto de lanzamiento es el apogeo.

  • Si la velocidad de lanzamiento v1 es mayor que vc el punto de lanzamiento es el perigeo.

En la figura, la trayectoria de color rojo es circular.

 

Ejemplos:

Ejemplo 1. La Tierra

La Tierra describe una órbita aproximadamente circular de radio r=1.496·1011 m=1 UA. Aplicando la dinámica del movimiento circular uniforme calculamos la velocidad de al Tierra. Sabiendo que la masa del Sol es M=1.98·1030 kg, la velocidad de la Tierra es v=29711.8 m/s. En el sistema de unidades establecido en el primer apartado

  • La longitud se mide en L=1.496·1011 m
  • El tiempo se mide en P=un año=365.26 días=3.156·107 s.

la velocidad v' vale v=v'·L/P

v'=6.268 UA/año

El tiempo que tarda en dar una vuelta completa es P=1 año

Cuando introducimos en el programa interactivo x=1.0 y vy=6.27, obtenemos una trayectoria circular alrededor del Sol.

Ejemplo 2. Marte

Los datos de Marte son
  • Semieje mayor a=1.524 UA,
  • Excentricidad ε=0.093.
  • La distancia más cercana al Sol (perihelio) es r1=a-c=a-εa=1.382 UA=2.068·1011 m

  • La distancia más alejada del Sol (afelio) es  r2=a+c=a+εa=1.666 UA=2.492·1011 m

Conocidas las distancias máxima y mínima r1 y r2 calculamos la velocidad en el perihelio v1

v1=26420.7 m/s=5.573 UA/año

Cuando introducimos en el programa interactivo x=1.382 y vy=5.573, obtenemos la trayectoria elíptica de Marte alrededor del Sol, y el tiempo que tarda en dar una vuelta completa P=1.86 años.

 

Actividades

En el applet que contiene esta página, se trazarán las trayectorias que describen los cuerpos celestes. Se comprobará la constancia  de la energía, se verificará que el momento angular es constante en las posiciones de máxima proximidad o de máximo alejamiento, y por último, se comprobará la tercera ley de Kepler, midiendo el periodo y el semieje mayor de la elipse.

Se introducen la posición y la velocidad inicial del cuerpo celeste::

  • La posición inicial x, la ordenada y=0
  • La componente Y de la velocidad inicial vy, la componente X de la velocidad vx=0.

Se pulsa el botón titulado Empieza,

Se traza la trayectoria del móvil, al mismo tiempo se muestra en la parte izquierda del applet, como van cambiando los valores de la posición y velocidad a medida que transcurre el tiempo. Observaremos que la energía y el momento angular permanecen constantes.

En la parte inferior izquierda, se muestra en color azul el tanto por ciento de error. Cuando es mayor que la unidad el programa interactivo se detiene. Como podemos comprobar, los mayores porcentajes de error se obtienen cuando la partícula pasa muy cerca del centro fijo de fuerzas.

Se pulsa el botón titulado Pausa, para parar el movimiento, por ejemplo, cuando el planeta pasa por la posición más próxima o más alejada, para medir el semieje mayor, la velocidad en dicha posición, y el semiperiodo (la mitad del tiempo que tarda el cuerpo celeste en dar una vuelta completa).

Se pulsa el mismo botón que ahora se titula Continua, para reanudar el movimiento..

Se pulsa varias veces en el botón Paso, para mover el cuerpo paso a paso, se usa para acercarnos a la posición deseada.

Cuando se ha completado una órbita se introduce la velocidad inicial en UA/año de un nuevo cuerpo sin cambiar la posición, y se pulsa el botón Empieza. Su trayectoria se traza en un color distinto.

Finalmente, se cambia la posición inicial x y se vuelven a introducir varios valores de la velocidad vy.

Cuando se han acumulado varias trayectorias, se pulsa en el botón Borrar para limpiar el área de trabajo del applet

Ejercicio:

kepler3.gif (2071 bytes)

Introduciendo la posición inicial del móvil Rp y la velocidad inicial Vp, se traza las sucesivas posiciones del planeta a intervalos fijos de tiempo.

Pulsando en los botones Pausa y Paso, se tomarán los siguientes datos y se completará una tabla como la siguiente.

Rp

Vp

Rp Vp

Ra

Va

Ra Va

a

P

P2/a3

                 
                 
                 

1.-Anotar en la primera y segunda columna las condiciones iniciales introducidas en los controles de edición: la distancia al perihelio Rp y la velocidad Vp

2.-Anotar en la cuarta columna de la tabla la distancia al afelio Ra en UA

3.-Anotar en la quinta columna de la tabla la velocidad en el afelio Va en UA/año

4.-Comprobar que (tercera y sexta columna de la tabla)

5.-Obtener el semieje mayor de la elipse a  a partir de las medidas de Rp y Ra, en la figura  2a=Rp+Ra

6.-Anotar el periodo P, el tiempo en años que tarda en dar una vuelta completa. Comprobar la tercera ley de Kepler (el cociente P2/a3 debe ser aproximadamente constante) en  la última columna de la tabla.

KeplerApplet1 aparecerá en un explorador compatible con JDK 1.1.

 

Código fuente

public abstract class RungeKutta {
	double h;
RungeKutta(double h){
	this.h=h;
}
void setPaso(double dt){
	this.h=dt;
}
public void resolver(Estado e){
	//variables auxiliares
	double k1, k2, k3, k4;
	double l1, l2, l3, l4;
	double q1, q2, q3, q4;
	double m1, m2, m3, m4;
//estado inicial
	double x=e.x;
	double y=e.y;
	double vx=e.vx;
	double vy=e.vy;
	double t=e.t;

	k1=h*vx;
	l1=h*f(x, y, vx, vy, t);
	q1=h*vy;
	m1=h*g(x, y, vx, vy, t);

	k2=h*(vx+l1/2);
	l2=h*f(x+k1/2, y+q1/2, vx+l1/2, vy+m1/2, t+h/2);
	q2=h*(vy+m1/2);
	m2=h*g(x+k1/2, y+q1/2, vx+l1/2, vy+m1/2, t+h/2);

	k3=h*(vx+l2/2);
	l3=h*f(x+k2/2, y+q2/2, vx+l2/2, vy+m2/2, t+h/2);
	q3=h*(vy+m2/2);
	m3=h*g(x+k2/2, y+q2/2, vx+l2/2, vy+m2/2, t+h/2);

	k4=h*(vx+l3);
	l4=h*f(x+k3, y+q3, vx+l3, vy+m3, t+h);
	q4=h*(vy+m3);
	m4=h*g(x+k3, y+q3, vx+l3, vy+m3, t+h);

	x+=(k1+2*k2+2*k3+k4)/6;
	vx+=(l1+2*l2+2*l3+l4)/6;
	y+=(q1+2*q2+2*q3+q4)/6;
	vy+=(m1+2*m2+2*m3+m4)/6;
	t+=h;

//estado final
	e.x=x;
	e.y=y;
	e.vx=vx;
	e.vy=vy;
	e.t=t;
}
abstract public double f(double x, double y, double vx, double vy, double t);
abstract public double g(double x, double y, double vx, double vy, double t);

}
public class Planeta extends RungeKutta{
Planeta(double h){
	super(h);
}
public double f(double x, double y, double vx, double vy, double t){
	double r=Math.sqrt(x*x+y*y);
	double z=-4*Math.PI*Math.PI*x/(r*r*r);
	return z;
}
public double g(double x, double y, double vx, double vy, double t){
	double r=Math.sqrt(x*x+y*y);
	double z=-4*Math.PI*Math.PI*y/(r*r*r);
	return z;
}
public double energia(double x, double y, double vx, double vy){
	double r=Math.sqrt(x*x+y*y);
	double z=(vx*vx+vy*vy)/2-4*Math.PI*Math.PI/r;
	return z;
}
}
//Objetos de la clase Planeta
	Estado estado=new Estado(0.0, x, 0.0, 0.0, vy); //estado inicial
	Planeta planeta=new Planeta(0.005);
 	eInicial=planeta.energia(estado.x, estado.y, estado.vx, estado.vy);
//rutina que calcula la trayectoria paso a paso
	planeta.resolver(valor);
	double energia=planeta.energia(estado.x, estado.y, estado.vx, estado.vy);
	double error=Math.abs((energia-eInicial)*100/eInicial);
	if(error>1.0){
		//detiene el movimiento
	}