El problema de Kepler

Supongamos un planeta, como la Tierra que se mueve en una órbita elíptica de semieje mayor a y de excentricidad ε=c/a con el Sol situado en el foco F. Sea el periodo, o tiempo que tarda en dar una vuelta completa, P. Al cabo de un cierto tiempo t, el planeta se encuentra en una posición (r, θ). El ángulo θ se denomina anomalía verdadera en el instante t.

Como hemos deducido en la página, la ecuación de la trayectoria en coordenadas polares es

r=( L 2 GM m 2 ) 1 1+εcosθ ε= 1+ 2 L 2 E m 3 G 2 M 2

La distancia más cercana al Sol r1 se obtiene en la ecuación de la trayectoria con θ=0 , y la más alejada r2 con θ=π. Como vemos en la figura r1+r2=2a

r 1 =( L 2 GM m 2 ) 1 1+ε r 2 =( L 2 GM m 2 ) 1 1ε 2a= r 1 + r 2 a=( L 2 GM m 2 ) 1 1 ε 2

La ecuación de la trayectoria en coordenadas polares se escribe

r= a( 1 ε 2 ) 1+εcosθ

Por ser la fuerza de atracción una fuerza central, el momento angular L es constante, y por ser conservativa la energía E es constante. El momento angular en coordenadas polares se escribe

L=m r 2 dθ dt

Combinando ambas expresiones, obtenemos la expresión de la posición angular θ en función del tiempo t.

0 θ dθ ( 1+εcosθ ) 2 = G 2 M 2 m 3 L 3 t

Tomamos el tiempo t=0, cuando el planeta pasa por el perihelio θ=0.

La integral de la izquierda se puede resolver con el cambio de variable

cosθ= 1 x 2 1+ x 2 sinθ= 2x 1+ x 2 dθ= 2 1+ x 2 dxx=tan( θ 2 ) dθ ( 1+εcosθ ) 2 = 2 ( 1+ε ) 2 ( 1+ x 2 ) ( 1+a x 2 ) 2 dxa= 1ε 1+ε 2 ( 1+ε ) 2 { 1 a dx ( 1+a x 2 ) + a1 a dx ( 1+a x 2 ) 2 }= 2 ( 1+ε ) 2 { 1 a dx ( 1+a x 2 ) + a1 a [ dx ( 1+a x 2 ) a x 2 ( 1+a x 2 ) 2 dx ] }= 2 ( 1+ε ) 2 { 1 a 3 2 tan 1 ( a x )+ a1 a ( 1 2 a tan 1 ( a x )+ x 2( 1+a x 2 ) ) }= 1 ( 1+ε ) 2 { a+1 a 3 2 tan 1 ( a x )+ a1 a x ( 1+a x 2 ) }

El primer paso es hacer el cambio de variable de cosθ y a x y dx. Después hay que descomponer la fracción integrando en dos fracciones, e integrar por partes.

Comprobamos con MATLAB Symbolic Toolbox que obtenemos la expresión entre corchetes sin multiplicar por 2/(1+ε)2

>> syms a x e fi;
>> y=int('(1+x^2)/(1+a*x^2)^2',x)
>> y =(atan(a^(1/2)*x)*(a + 1))/(2*a^(3/2)) + 
(x*(a - 1))/(2*a*(a*x^2 + 1))       

El resultado después de deshacer los cambios de a a ε y de x a θ es

1 ( 1 ε 2 ) 3 2 ( 2 tan 1 ( 1ε 1+ε tan θ 2 ) ε 1 ε 2 sinθ 1+εcosθ )= G 2 M 2 m 3 L 3 t

El resultado que nos proporciona MATLAB es apenas legible, incluso utilizando la función pretty

>> yy=(2/(1+e)^2)*subs(y,{a,x},{(1-e)/(1+e),tan(fi/2)});
>> simplify(yy)
ans =-(2*(atan(tan(fi/2)*(-(e - 1)/(e + 1))^(1/2))/
((-(e - 1)/(e + 1))^(1/2)*(e - 1)) - (e*sin(fi)*(e + 1))
/(2*(e - 1)*(e*cos(fi) + 1))))/(e + 1)^2

que se puede escribir como

2 tan 1 ( 1ε 1+ε tan θ 2 ) ε 1 ε 2 sinθ 1+εcosθ = ( 1 ε 2 ) 3 2 G 2 M 2 m 3 L 3 t

Que es una ecuación implícita θ=θ(t)

El periodo P es el tiempo que tarda en dar una vuelta completa un planeta, que describe una órbita elíptica

P= 2π a 3 2 GM P=2π L 3 G 2 M 2 m 3 1 ( 1 ε 2 ) 3 2

Escribimos la ecuación implícita θ=θ(t) en términos del periodo P

2 tan 1 ( 1ε 1+ε tan θ 2 ) ε 1 ε 2 sinθ 1+εcosθ =2π t P

El miembro izquierdo se denota como Me (anomalía media)

M e =2π t P

Representamos Me en función de la posición angular θ, para los siguientes valores de la excentricidad ε=0, 0.2, 0.5 0.8 y 0.9

x=(0:360)*pi/180;
hold on
for e=[0,0.2,0.5,0.8,0.9]
    y=2*atan(sqrt((1-e)/(1+e))*tan(x/2))-e*sqrt(1-e^2)*sin(x).
/(1+e*cos(x));
    yy=y.*(y>0)+(2*pi+y).*(y<0);
    plot(x,yy, 'displayName',num2str(e))
end
legend('-DynamicLegend','Location','northwest')
set(gca,'XTick',0:pi/2:2*pi)
set(gca,'XTickLabel',{'0','\pi/2','\pi','3\pi/2','2\pi'})
set(gca,'YTick',0:pi/2:2*pi)
set(gca,'YTickLabel',{'0','\pi/2','\pi','3\pi/2','2\pi'})
hold off
grid on
xlabel('\theta')
ylabel('M_e')
title('La ecuación del tiempo')

Trazamos la trayectoria elíptica de semieje mayor a y excentricidad ε en color rojo. c=ε·a es la semidistancia focal entre O y el foco F donde se encuentra el Sol. Trazamos una circunferencia de radio a en color azul y dibujamos una recta perpendicular al eje mayor que corta a la circunferencia en el punto Q. El ángulo de la recta OQ se denomina E, anomalía excéntrica del planeta en el instante t. Obtenemos las siguientes relaciones geométricas: acosE=c+rcosθ

acosE=aε+ a( 1 ε 2 ) 1+εcosθ cosθ cosE= ε+cosθ 1+εcosθ sin 2 E+ cos 2 E=1sinE= 1 ε 2 sinθ 1+εcosθ

A partir de la primera relación entre cosE y cosθ obtenemos la ecuación de la trayectoria en función de E

r= a( 1 ε 2 ) 1+εcosθ r= a( 1 ε 2 ) 1+ε( εcosE εcosE1 ) =a( 1εcosE )

Utilizando la relación trigonométrica de la tangente del ángulo mitad, llegamos a una relación entre los ángulos E y la posición angular θ

tan 2 ( E 2 )= 1cosE 1+cosE tan 2 ( E 2 )= 1ε 1+ε 1cosθ 1+cosθ = 1ε 1+ε tan 2 ( θ 2 ) tan( E 2 )= 1ε 1+ε tan( θ 2 ) E=2 tan 1 ( 1ε 1+ε tan( θ 2 ) )

Escribimos la ecuación implícita θ=θ(t)

M e =2 tan 1 ( 1ε 1+ε tan θ 2 ) ε 1 ε 2 sinθ 1+εcosθ M e =EεsinE

Esta ecuación mucho más simple se denomina ecuación de Kepler. Representamos Me en función del ángulo E.

x=(0:360)*pi/180;
hold on
for e=[0,0.2,0.4,0.6,0.8,1.0]
    y=x-e*sin(x);
    plot(x,y, 'displayName',num2str(e))
end
legend('-DynamicLegend','Location','northwest')
set(gca,'XTick',0:pi/2:2*pi)
set(gca,'XTickLabel',{'0','\pi/2','\pi','3\pi/2','2\pi'})
set(gca,'YTick',0:pi/2:2*pi)
set(gca,'YTickLabel',{'0','\pi/2','\pi','3\pi/2','2\pi'})
grid on
hold off
xlabel('E')
ylabel('M_e')
title('La ecuación del tiempo')

Dada la posición angular θ del satélite determinar el tiempo t

Dada la posición angular θ, determinamos el ángulo E mediante la relación

E=2 tan 1 ( 1ε 1+ε tan( θ 2 ) )

Sustituyendo en la ecuación de Kepler, obtenemos Me

M e =EεsinE

Conocido Me obtenemos el tiempo

t= M e 2π P

Ejemplo

Un satélite artificial describe una órbita elíptica con cuyas distancias máxima y mínima al centro de la Tierra son, respectivamente r1=9.6·106 m y r2=21·106. Determinar el instante t en el que la posición del satélite es θ=2π/3 (120°)

2a= r 2 + r 1 a=15.3 10 6 m c=a r 1 ε= c a = r 2 r 1 r 2 + r 1 =0.37 E=2 tan 1 ( 1ε 1+ε tan( θ 2 ) ) E=2 tan 1 ( 10.37 1+0.37 tan( π 3 ) )=1.728 M e =EεsinE M e =1.7280.37sin1.728=1.36 P 2 = 4 π 2 GM a 3 = 4 π 2 6.67 10 11 5.98 10 24 ( 15.3 10 6 ) 3 P=18828s t= M e 2π P t= 1.36 2π 18828=4076s

El instante en el que se alcanza la posición angular θ=2π/3 es t=4076 s

Dado el tiempo t determinar la posición angular θ del satélite

Dado el tiempo t, calculamos Me

M e =2π t P

Sustituyendo Me en la ecuación de Kepler, obtenemos una ecuación trascendente cuya raíz tenemos que determinar por procedimientos numéricos

EεsinE= M e

Un procedimiento apropiado para resolver esta ecuación trascendente es el método de Newton. Sea f(x)=0. El proceso iterativo para calcular la raíz es

x i+1 = x i f( x i ) f'( x i ) E i+1 = E i E i εsin E i M e 1εcos E i

Se repite el proceso iterativo hasta que se encuentra la raíz E con el nivel deseado de precisión. Una vez determinado E se calcula la posición angular θ mediante la relación

tan( θ 2 )= 1+ε 1ε tan( E 2 )

Ejemplo

Encontrar posición angular θ del satélite del ejercicio anterior en el instante t=10800 s (3 horas)

M e =2π t P M e =2π 10800 18828 =3.60

Resolvemos la ecuación trascendente E-ε·sinE=Me por el procedimiento de Newton tomando como valor inicial de partida E0=Me/2=3.418

x0=3.418;
e=0.3725;
Me=3.604;
while(1)
    x=x0-(x0-e*sin(x0)-Me)/(1-e*cos(x0));
    if abs((x-x0)/x0)<1e-6
        break;
    end
    x0=x;
end
fprintf('La raíz buscada es %1.3f\n',x0)
La raíz buscada es 3.480

Si E=3.480, a partir de la relación

tan( θ 2 )= 1+0.37 10.37 tan( 3.48 2 )θ=3.372

obtenemos la posición angular θ=3.372 del satélite en el instante t=10800 s (3 horas)

Desarrollo en serie

Dado que no existe una solución analítica de la ecuación de Kepler, se han propuesto desarrollos en serie de infinitos términos que en la práctica hay que limitar a un número finito N.

El propuesto por Lagrange converge para excentricidades que no excedan el valor límite 0.6627. Para la Tierra que tiene una excentricidad muy pequeña es suficiente con los primeros tres términos N=3 del desarrollo en serie

E= M e + n=1 a n ε n E M e +εsin M e + ε 2 2 sin( 2 M e )+ ε 3 8 ( 3sin( 3 M e )sin M e )

Me=0:pi/180:2*pi;
e=0.6;
E=Me+e*sin(Me)+e^2*sin(2*Me)/2+e^3*(3*sin(3*Me)-sin(Me))/8;
x=0:pi/180:2*pi;
y=x-e*sin(x);

plot(E,Me,'r',x,y,'b')
xlabel('E')
ylabel('M_e');
legend('aproximado','exacto','Location', 'northwest')
title('Ecuación del tiempo')

Otra alternativa que no tiene las limitaciones del método de Lagrange, es el desarrollo en serie cuyos términos son funciones Jn(x) de Bessel

E= M e + n=1 N 2 n J n ( nε ) sin( n M e )

Me=0:pi/180:2*pi;
e=0.6;

N=3;
E=Me;
for n=1:N
    E=E+2*besselj(n,n*e)*sin(n*Me)/n;
end
x=0:pi/180:2*pi;
y=x-e*sin(x);

plot(E,Me,'r',x,y,'b')
xlabel('E')
ylabel('M_e');
legend('aproximado','exacto','Location', 'northwest')
title('Ecuación del tiempo')

Referencias

Orbital Mechanics for Engineering Students. Howard D. Curtis. Elsevier, Capítulo 3