Escuela Técnica Superior de Ingeniería del Diseño

Universitat Politècnica de València

CÁLCULO DE TRAYECTORIAS DE
SATÉLITES

TRABAJO FIN DE GRADO

Autor: Diego Expósito Brioso

Tutor: Sergio Blanes Zamora

Grado en Ingeniería Aeroespacial

4º Curso, Junio 2017
 Resumen

En este proyecto se presenta el estudio necesario para calcular con ele-

vada precisión la trayectoria no ideal de los satélites, empleando métodos
numéricos específicos para este problema. Se comienza con una introducción
histórica al movimiento de los cuerpos celestes, a la cual sigue la definición del
problema abordado. Posteriormente se justifican las limitaciones que tienen
las herramientas de cálculo generales, lo que permite introducir los nuevos
métodos numéricos desarrollados recientemente, y la filosofía en la que se
basan. Tras ello, se aumenta la complejidad del estudio para reflejar sus
aplicaciones y desarrollos, y se expone un método muy preciso empleado ac-
tualmente para resolver nuestro problema. Por último, se generan archivos
de software que permiten la utilización de dichos métodos por parte de un
usuario general.

Palabras clave: Métodos numéricos, Simplecticidad, Hamiltoniano, Kepler,

Error, Paso, Perturbación.

3
Índice general

1. Introducción a la Mecánica Celeste 7

1.1. El comienzo . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2. La teoría de perturbaciones . . . . . . . . . . . . . . . . . . . 8
1.3. Los métodos numéricos . . . . . . . . . . . . . . . . . . . . . . 10
1.4. Finalidad del estudio . . . . . . . . . . . . . . . . . . . . . . . 11

2. Presentación del problema 13

2.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2. Achatamiento de la Tierra . . . . . . . . . . . . . . . . . . . . 16
2.3. Rozamiento con la atmósfera . . . . . . . . . . . . . . . . . . . 16
2.4. Oscilador armónico: el caso lineal . . . . . . . . . . . . . . . . 17
2.5. Planteamiento . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3. Métodos numéricos 21
3.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.2. Métodos numéricos generales . . . . . . . . . . . . . . . . . . . 22
3.2.1. Euler explícito . . . . . . . . . . . . . . . . . . . . . . . 22
3.2.2. Método Adams-Bashforth de tres pasos . . . . . . . . . 24
3.2.3. Método de Adams-Moulton de tres pasos . . . . . . . . 26
3.3. Comparación de los métodos . . . . . . . . . . . . . . . . . . . 28

4. Métodos simplécticos 31
4.1. Simplecticidad . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
4.2. Métodos simplécticos básicos . . . . . . . . . . . . . . . . . . . 34
4.3. Backward error analysis . . . . . . . . . . . . . . . . . . . . . 40

5. Aumento del orden 43

5.1. Estudio teórico . . . . . . . . . . . . . . . . . . . . . . . . . . 43
5.1.1. Composición . . . . . . . . . . . . . . . . . . . . . . . . 43
5.1.2. Splitting . . . . . . . . . . . . . . . . . . . . . . . . . . 44
5.1.3. Sistemas casi-integrables . . . . . . . . . . . . . . . . . 48

5
6 ÍNDICE GENERAL

5.2. Aplicaciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
5.2.1. Método NIA . . . . . . . . . . . . . . . . . . . . . . . . 53
5.2.2. Ode45 . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
5.3. Sistema caótico . . . . . . . . . . . . . . . . . . . . . . . . . . 56

6. Métodos actuales 61
6.1. Método ABAH844 . . . . . . . . . . . . . . . . . . . . . . . . 61
6.2. Evolución del sistema . . . . . . . . . . . . . . . . . . . . . . . 63
6.3. Conclusiones . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

A. Software 65
A.1. Kepler no perturbado . . . . . . . . . . . . . . . . . . . . . . . 65
A.2. Kepler perturbado . . . . . . . . . . . . . . . . . . . . . . . . 68

Bibliografía 75
Capítulo 1

Introducción a la Mecánica
Celeste

1.1. El comienzo
La astronomía es una de las mayores pasiones del ser humano, que aparece
en numerosas culturas a lo largo de toda la Historia. Siendo tan fundamental
la influencia del Sol en nuestra vida cotidiana, y a la vez, estando tan fuera
de nuestro alcance, el movimiento de los cuerpos celeste ha sido sujeto de
numerosas investigaciones, que incluyen la recogida de datos empíricos y la
formulación de teoremas matemáticos.
Uno de los mayores avances en este campo, una vez asumido el modelo
heliocéntrico defendido por Nicolás Copérnico (1473 − 1543), son lo que se
conocen como las leyes de Kepler. Tras acceder a los datos sobre órbitas
planetarias, recogidas durante años por Tycho Brahe (1546−1601), Johannes
Kepler (1571 − 1630) llegó a deducir lo que se conocen como las Tres Leyes
de Kepler:

1. Los planetas se mueven alrededor del Sol describiendo una elipse, es-
tando éste en uno de sus dos focos.

2. Las áreas barridas por los radios de los planetas son proporcionales al
tiempo empleado por éstos en recorrer el perímetro de dichas áreas.

3. El cuadrado de los períodos de la órbita de los planetas es proporcional

al cubo de la distancia al Sol.

Esta última ley es especialmente importante ya que interpreta, por pri-

mera vez, el movimiento entre dos astros.

7
8 Cálculo de trayectorias de satélites

En 1687, Isaac Newton (1643 − 1727) presenta en su Philosophiae Natu-

ralis Principia Mathematica la ley universal de la gravedad, que postula que
dos cuerpos de masas m1 , m2 situados a una distancia relativa r se atraen
con una fuerza:
m1 m2
F =G
r2
2
donde G es la constante gravitacional G = 6. 674 · 10−11 N kg
m
2.

Puesto que el mismo autor determina la fuerza con la ecuación diferencial:

d(mv)
F =
dt
se tiene que, por primera vez, el movimiento relativo entre dos astros viene
dado por una ecuación diferencial. Ésto supone el comienzo de la Mecánica
Celeste, donde el movimiento de los planetas se determina utilizando las
ecuaciones diferenciales, las cuales han demostrado ser herramientas muy
poderosas de cara a modelizar los fenómenos físicos de la realidad.
Agrupamos en la siguiente tabla las contribuciones más importantes a la
Mecánica Celeste desde Newton [5]:

Autor Contribución
Cálculo de órbitas de
Edmund Halley (1656 − 1742)
cometas
Teoría lunar para nave-
Leonhard Euler (1707 − 1783)
gación marítima
Introducción del méto-
Joseph Luis Lagrange (1736 − 1813) do de variación de pa-
rámetros
Teoría de perturbacio-
Pierre Simon de Laplace (1749 − 1827)
nes
Problema de tres cuer-
Henri Poincaré (1854 − 1912)
pos, teoría del caos

1.2. La teoría de perturbaciones

La ley gravitacional de Newton conecta el concepto de trayectoria con el
de fuerza en el movimiento orbital, al igual que ocurre, de forma más evidente,
en una bola empujada por una fuerza inicial. Los planetas giran alrededor del
Sol por esta ley gravitacional; sin embargo, puesto que depende únicamente
de las masas, los planetas son susceptibles de atraerse también entre ellos.
 Introducción a la Mecánica Celeste 9

La razón por la que las trayectorias se producen alrededor del Sol es

simplemente que éste tiene, con diferencia, la mayor masa de todo el Sistema
Solar (la masa del planeta más pesado, Júpiter, apenas es una milésima
parte que la del Sol). Las fuerzas que se producen entre los planetas son
mucho menores que la atracción al Sol, considerándose ”perturbaciones” de
la atracción ideal de dos cuerpos.
Sin considerar estas perturbaciones, el cálculo de las órbitas resulta ser
erróneo. Por ello, y debido a que la complejidad de las ecuaciones diferenciales
aumenta mucho al tener en cuenta las perturbaciones, se popularizó el uso de
expansiones en serie de funciones que permitían conocer una aproximación
analítica al problema. A falta de computadores, éste fue el método utilizado
para obtener la influencia de perturbaciones. Se produjo un gran desarrollo
del cálculo diferencial, expansiones en series trigonométricas,...
La forma más sencilla de obtener la representación de una función como
serie viene dada por la expansión en serie de Taylor, donde si se conocen los
valores de una función f y el de sus derivadas en un punto x0 , en un punto
siguiente x0 + h el valor de la función viene dada por:

∞
f (n) (x0 ) n
f (x0 + h) =
X
h
n=0 n!
Cuanto más términos se incluyen en la aproximación, tanto más acertada
ésta se vuelve; si se trunca el desarrollo hasta una serie de términos, cuanto
menor sea h, menor será el error cometido por la aproximación.
El cálculo de las derivadas de f es cada vez más costoso conforme se
va derivando f , lo cual limitaba la precisión que se podía obtener con este
método.
El siguiente suceso es una prueba del éxito que han llegado a alcanzar los
métodos numéricos y la teoría de perturbaciones [6]. A finales del siglo XVIII,
sólo se conocía la existencia de los planetas exteriores Júpiter, Saturno y
Urano. Alexis Bouvard (1767 − 1843), estudiante del famoso Pierre-Simon de
Laplace (1749 − 1827), fue otorgado por éste la tarea (poco grata) de calcular
numéricamente, haciendo uso de tablas, las órbitas de estos planetas. Las
ecuaciones de perturbaciones entre dos planetas ya fueron determinadas por
el propio Laplace, con lo que solo quedaba calcular las posiciones resultantes.
Bouvard consiguió calcular efectivamente las órbitas de Júpiter y Saturno;
la de Urano difería constantemente de sus cálculos, lo que indicaba que otro
cuerpo podía estar perturbando su órbita. El matemático John Couch Adams
(1819 − 1892) desarrolló por primera vez los métodos multipaso lineales, que
le permitieron describir, con mucha precisión, la órbita de un cuerpo exterior
a Urano. El descubrimiento de Neptuno se realizaría simultáneamente por el
10 Cálculo de trayectorias de satélites

francés Urbain Le Verrier (1811−1877), astrónomo del Observatorio de París,

adelantándose a Adams, simplemente porque este último perdió el interés por
el asunto.

Figura 1.1: Cálculo de órbitas de Neptuno por Adams y Le Verrier, [6]

1.3. Los métodos numéricos

Para tener en cuenta un número creciente de perturbaciones, y así obtener
una representación más fiel de la realidad, se requería de métodos numéricos
cada vez más precisos (además de máquinas capaces de realizar los cálculos,
lo cual pertenece al campo de la Informática y no será tratado aquí). El
método de Euler consiste en el truncamiento de la serie de Taylor al primer
término:

f (xn + h) ≈ f (xn ) + hf 0 (xn )

lo que representa una integración de f aproximada por la regla de cuadratura:
Z xn
f (xn ) dx ≈ (xn − xn−1 )f (xn−1 )
xn−1

La idea del matemático Carl David Tomé Runge (1856 − 1927) en su

publicación de 1895 fue desarrollar métodos numéricos basados en mejores
aproximaciones de integración que ya se conocían en la época, como el método
del trapecio y el método del punto medio. Karl Heun (1859−1929) aumentaría
 Introducción a la Mecánica Celeste 11

el orden de estos métodos en 1900. Más tarde, en 1901, Runge, conjuntamente

con Martin Wilhelm Kutta (1867 − 1944), publicarían lo que actualmente se
conocen como métodos Runge-Kutta, cuya utilidad para resolver ecuaciones
diferenciales los volvió muy populares [7].
De cara a calcular el movimiento de los planetas en intervalos de tiempos
extremadamente amplios (millones de años), se requería aumentar el orden
de los métodos numéricos, hasta tal punto, que estos debían estar diseñados
para resolver de forma específica el problema, para conservar sus propiedades.
Puesto que el movimiento celeste está basado en un problema Hamiltoniano,
donde se conservan una serie de propiedades geométricas, se empezaron a
desarrollar métodos numéricos específicos para problemas Hamiltonianos, lo
que marca el comienzo de lo que se conoce como ”Integración Geométrica”.
Un problema es Hamiltoniano cuando el flujo de sus soluciones es una
transformación simpléctica, la cual representa una propiedad geométrica del
sistema. Al utilizar un algoritmo numérico general, no se respeta el carácter
simpléctico de las soluciones y se pierde las propiedades geométricas, lo que
lleva a aumentar el error. Surge así la necesidad de construir métodos numé-
ricos tales que el flujo de las soluciones que aportan sigue siendo simpléctico.
Aunque este concepto fue introducido por De Vogelaere en 1956, no fue
hasta comienzos de los años ochenta cuando se empezaron a publicar inves-
tigaciones sobre este tema. El suizo F.M. Lasagni, el soviético Y.B. Suris y
el español J.M. Sanz-Serna fueron los pioneros en el estudio de estos méto-
dos numéricos. A día de hoy, se siguen desarrollando y expandiendo con la
finalidad de resolver con mayor eficacia una diversa gama de ecuaciones dife-
renciales [4]. Un ejemplo de aplicación práctica de esta teoría es [11], donde se
calcula la evolución temporal inversa del clima de la Tierra desde el periodo
actual hasta el Paleógeno, 66 millones de años atrás.

1.4. Finalidad del estudio

El presente trabajo constituye el Trabajo Final de Grado, así como el
proyecto de Beca de Colaboración 2017, del alumno Diego Expósito Brioso,
del Grado en Ingeniería Aeroespacial de la Universitat Politècnica de Va-
lència, bajo la tutela del Profesor Titular de la Universitat Politècnica de
València, Sergio Blanes Zamora, quien desarrolla proyectos de investigación
en integración geométrica en el Instituto de Matemática Multidisciplinar. Los
conocimientos presentados se extraen principalmente del libro [1].
El objetivo del trabajo es presentar la teoría de la integración geométrica
de una forma amena para el público científico general, y el cómo se aplica
dicha teoría para resolver un problema de ingeniería espacial, concretamente
12 Cálculo de trayectorias de satélites

el cálculo realista de órbitas de satélites alrededor de la Tierra. Se comien-

za con una presentación del problema a tratar, presentando el problema de
Kepler ideal para trayectorias no perturbadas, y presentando el comporta-
miento de las perturbaciones posteriormente. Acto seguido, se hace un repaso
de la teoría de los métodos numéricos generales, para más adelante, expli-
car el diseño de los métodos numéricos simplécticos y el cómo éstos pueden
superar a los métodos numéricos estándar. Se proporciona un ejemplo de
método simpléctico utilizado en la actualidad para obtener las trayectorias
en Mecánica Celeste. Finalmente, toda la teoría de la integración simpléctica
se condensa en un conjunto de paquetes de software para el cálculo numérico
de trayectorias disponibles para su uso como ”caja negra” por parte de un
usuario no versado en dicha teoría.
Capítulo 2

Presentación del problema

2.1. Introducción
El problema que vamos a estudiar es la trayectoria de un satélite alre-
dedor de un planeta considerando efectos reales que ejercen perturbaciones
sobre el caso ideal, conocido como problema no perturbado de Kepler. Es
especialmente interesante de cara a conocer la posición y velocidad de los
satélites en intervalos de tiempo largos, donde esperamos que nuestras solu-
ciones de posición y velocidad sean correctas por haber considerado todos los
efectos reales.
Recientemente, expertos del Department of Aerospace and Mechanical
Engineering de la Universidad de Lieja estudiaron el mismo problema [8],
pero con un enfoque distinto al nuestro. Considerando las perturbaciones del
rozamiento del satélite con la capa de la atmósfera más externa del planeta,
así como el achatamiento del propio planeta, decidieron resolver la ecuación
diferencial resultante mediante métodos analíticos:

r
r̈ + µ = fa + fr .
||r||3
Para ello, emplearon herramientas matemáticas avanzadas como el álge-
bra de Lie, transformaciones,... que les permitieron, asumiendo que la ex-
centricidad del planeta fuera pequeña, y utilizando un modelo de atmósfera
exponencial, obtener la solución analítica del problema.
Sin embargo, en dicho documento se detallan los inconvenientes de la
solución. Al haber obtenido un solución cerrada del problema, el efecto de
otras perturbaciones, como pueden ser la fuerza de la radiación del Sol sobre el
satélite, la no-homogeneidad de la capa alta de la atmósfera,... no se pueden
considerar; en apenas dos días, el error de posición estimada del satélite

13
14 Cálculo de trayectorias de satélites

alcanza los 10 kilómetros, lo que hace que simulaciones a largo plazo sean
inasumibles como representación correcta de la realidad.
Incorporar los nuevos efectos reales a la ecuación modifica completamente
la solución a obtener; si se quiere tener en cuenta cientos de efectos, el estudio
analítico no es práctico para obtener la solución del problema.

El planteamiento que nuestro estudio propone es el siguiente: suponga-

mos que queremos tener en cuenta n perturbaciones sobre el problema de
Kepler. Para un número n elevado, obtener la solución analítica es compli-
cado por la interferencia de los términos entre ellos; sin embargo, obtener la
solución analítica de solo una de ellas es sencillo y puede realizarse incluso
con software de cálculo simbólico. Entonces, si tenemos n soluciones exactas
de cada perturbación por separado, ofrecemos un método numérico que ob-
tenga la solución numérica del problema, combinando las n soluciones de las
perturbaciones. De esta forma, nuestro modelo puede cambiar con facilidad
con las perturbaciones que se quieran considerar. Nuestro desafío es, enton-
ces, el cómo se realiza dicha combinación para poder obtener el mínimo error
en la solución, con un paso de tiempo lo más alto posible, de cara a obtener
soluciones en intervalos de tiempo muy largos.

De ahora en adelante, implementaremos los métodos numéricos en Matlab,

y utilizaremos la nomenclatura de la mecánica Hamiltoniana. Ésta requiere
de dos coordenadas independientes (q, p), donde q representa la posición de
un cuerpo, y p su cantidad de movimiento. Si se estudia el problema en d
dimensiones, entonces dichas coordenadas son vectores q ∈ <d , p ∈ <d .
Si el problema no tiene en cuenta elementos perturbadores, como la in-
fluencia de otros planetas, el rozamiento con la atmósfera, efectos relativis-
tas,... se puede definir el Hamiltoniano como una función escalar que es la
suma de las energías cinética y potencial, que en nuestro caso tiene la forma:

1 1
H(q, p) = T (p) + V (q) = pT p − µ .
2 r
El Hamiltoniano es la base de todos los sistemas físicos que se denominan
Hamiltonianos, donde la energía total del sistema se conserva, como en los
problemas de dinámica molecular, mecánica cuántica, estadística clásica,...
La ley gravitacional que determina la trayectoria de dos masas puntales
que se atraen por el efecto de la gravedad es:
q
q̈ = −µ
,
r3
donde µ = GM , G es la constante gravitacional, M la suma de las masas
 Presentación del problema 15
q
y r = kqk= q T q es la distancia entre ambos cuerpos. La incógnita q es la
distancia relativa entre ambos cuerpos.
Considerando las condiciones iniciales (q(0), p(0)) = (q0 , p0 ), la solución
está proporcionada por [9]:

q(t) = f q0 + g p0 , p(t) = fq q0 + gp p0 ,

con

(cos x − 1)a sin x − x

f =1+ , g =t+ ,
r0 w

aw sin x cos x − 1
fp = − , gp = 1 + .
r0 (1 − σ cos x + ψ sin x) 1 − σ cos x + ψsinx

x se obtiene de

wt = x − σ sin x + ψ(1 − cos x),

la cual es una ecuación implícita, con lo que se puede obtener x con el método
de Newton-Raphson:
 
x[j−1] − σs[j−1] + ψ(1 − c[j−1] − wt) 
x[j] = x[j−1] −  ,
1 − σc[j−1] + ψs[j−1]

donde c[j−1] = cos(x[j−1] ), s[j−1] = sin(x[j−1] ).

El resto de parámetros son

u r0 µ µ 1 1
r
ψ= , σ = 1− , w= , a=− , E = pT0 p0 −µ ,
wa2 a a3 2E 2 r0

u = q0T p0 , r0 = kq0 k.

Esta solución teórica será utilizada para determinar la fiabilidad de los

distintos métodos numéricos, cuando los comparemos para resolver el proble-
ma sin perturbaciones. A la hora de calcular la trayectoria real, la solución
analítica no es válida debido a dichas perturbaciones. Analizamos las más
importantes y cómo afectan al sistema.
16 Cálculo de trayectorias de satélites

2.2. Achatamiento de la Tierra

El achatamiento de la Tierra provoca un efecto sobre la trayectoria del
satélite, de forma que el problema sigue siendo Hamiltoniano, pero con un
Hamiltoniano perturbado. Considerando un caso estrictamente bidimensio-
nal con q = (qx , qy )T , p = (px , py )T , la estimación de primer orden de la
perturbación viene dada por:

1 1 3q 2
!
εa
H = (p2x + p2y ) − µ − 3 1 − α 2x ,
2 r 2r r
donde εa es una constante cuyo valor para la Tierra es del orden de εa ∼ 10−3 ,
y α determina el plano en el cual se produce el movimiento. Si se produce
en el Ecuador se tiene que α = 1, mientras que si se produce en un plano
perpendicular al Ecuador, α = 0 [1].

2.3. Rozamiento con la atmósfera

Si la altura del satélite es lo suficientemente baja como para encontrar
un pequeño rozamiento con la atmósfera, como ocurre en las órbitas LEO
(Low Earth Orbit), el problema deja de ser conservativo, lo cual lleva a
replantear las ecuaciones Hamiltonianas.
Por aerodinámica, se sabe que la fuerza de arrastre sobre una esfera viene
dada por:

1
D = ρ(q)Sref CD ||q̇|| q̇
2
donde ρ(q) es la densidad del aire (que varía con la altura), Sref un área de
referencia que suele ser el área seccional de la esfera, y CD el coeficiente de
arrastre, que tiene en cuenta detalles como la geometría, fricción,...
La densidad ρ(q) variará con la altura dependiendo del modelo de atmós-
fera que se escoja. Tomando un modelo de atmósfera exponencial como en el
trabajo [8], ésta adopta la expresión:

1
!
Sref r − rref
ρ(q) = CD ρref exp −
2 m α
siendo m la masa del satélite, rref una altura de referencia (la distancia del
centro de la Tierra a la superficie), ρref la densidad a esa altura y α una
constante de escalado. Teniendo en cuenta que q̈ = −D/m − µ rq3 , se tiene
que la contribución del rozamiento al problema es:
 Presentación del problema 17

q 1 2 Sref
2
! !
r − rref q r−a
q̈ = −µ 3 − CD ρref exp − ||q̇|| q̇ = −µ 3 −εr exp − ||q̇|| q̇
r 4 m 2 α r b
S2
donde εr = 14 CD 2 ref
ρ
m2 ref
es una constante, a = rref y b = α determinan la
distribución exponencial.
No es estrictamente correcto emplear la nomenclatura Hamiltoniana para
este problema, ya que al tratarse de una fuerza no conservativa, este problema
ya no es Hamiltoniano, la energía no se conserva. Por ello, q y p representarán,
en esta situación, ecuaciones diferenciales q̇ = p, q̈ = ṗ. El conocido efecto
del rozamiento en la realidad es la disminución de la altura de los satélites,
que hay que tener en cuenta para corregir su trayectoria.

2.4. Oscilador armónico: el caso lineal

Puesto que en este trabajo se van a tratar temas avanzados sobre métodos
numéricos, aplicados a un problema Hamiltoniano no lineal, nos interesa, para
clarificar conceptos y poder seguir más fácilmente los razonamientos, observar
también cómo afectan los distintos análisis a un problema Hamiltoniano lineal
más sencillo.
El modelo más estudiado como banco de pruebas para métodos numéricos
es el del oscilador armónico simple, que constituye un sistema Hamiltoniano
lineal. En el caso de una dimensión, considerando una partícula de masa m
atada a un muelle elástico de constante k, la ecuación diferencial que lo rige
es:

d2 y
m 2 = mÿ = −ky.
dt
Pasando a coordenas Hamiltonianas, q = y, p = mv = mq̇, se tiene
p
q̇ = , ṗ = −kq,
m
con un Hamiltoniano
1 2 1 2
H(q, p) =p + kq .
2m 2
De forma matricial, las ecuaciones se expresan como:

0
! ! ! !
d 1
q q q
= m =A .
dt p −k 0 p p
con
18 Cálculo de trayectorias de satélites

0
!
1
A= m (2.1)
−k 0
La solución analítica de este problema es sencilla de obtener. Si conside-
ramos conocidos (q(0), p(0)) = (q0 , p0 ):

cos ωt sin ωt
! ! ! !
1
q(t) q0 q0
= ω = Mt ,
p(t) −ω sin ωt cos ωt p0 p0
q
con ω = k/m. Se puede observar que para el oscilador armónico, como la
dependencia entre variables es lineal, se puede emplear la notación matricial,
lo que facilita el estudio del sistema.
Con ello, estudiaremos primero cómo afectan los métodos numéricos a
este caso simplificado, y después nos enfrentaremos al problema de Kepler
con más seguridad.

2.5. Planteamiento
A continuación ofrecemos la base sobre la que estudiaremos el problema.
Decidimos trabajar en dos coordenadas (d = 2), lo que nos lleva a cuatro
variables canónicas, dos de posición y dos de velocidad: q1 , q2 , p1 , p2 , relacio-
nados por el siguiente sistema de ecuaciones para el problema de Kepler no
perturbado:
q
q̇ = p, ṗ = −µ ,
r3
q
donde r = q T q, q = (q1 , q2 )T , p = (p1 , p2 )T . Finalmente, nuestras variables
sobre las que aplicaremos el esquema son:
 
q1
q2
 
(q, p) =T
 
 ,

 p1 

p2
con lo que
   
q1 p1
d d  q2 p2 q T
   
(q, p)T = = = (p, −µ ) .
  
−µ rq13
  
dt dt 
 p1 




 r3
p2 −µ rq23
 Presentación del problema 19

Para realizar los estudios, consideraremos los siguientes valores de refe-

rencia de los parámetros del problema no perturbado, salvo si se especifican
otros datos:

1+e
s
q1 (0) = 1 − e, q2 (0) = 0, p1 (0) = 0, p2 (0) = ,
1−e
tf − t0
t0 = 0, tf = 100, N = 2000, h= , e = 0. 2, µ = 1.
N
(2.2)

Ésto indica que tratamos una órbita genérica con un satélite situado a una
distancia q1 (0) del cuerpo central, con una velocidad p2 (0), una excentricidad
e = 0. 2 y un balance de masas µ = 1. Este balance de masas nos permite
definir la escala de tiempo, de forma que el periodo del cuerpo móvil para
una revolución es T = 2π. Pasamos ahora a analizar los métodos numéricos
generales.
Capítulo 3

Métodos numéricos

3.1. Introducción
En esta sección se van a discutir algunos de los distintos métodos numé-
ricos generales, de amplio uso comercial, como son los métodos multipaso,
Runge-Kutta,... y cómo resuelven nuestro problema específico.
Para ello, primero introduciremos la filosofía sobre la que se basan los
métodos numéricos, basándonos en [2].
Un problema de valores iniciales consiste en, dada una ecuación diferen-
cial, obtener la función solución y(t) de la que se conoce su valor en el punto
inicial t0 . La derivada y 0 (t) es función de y y de t, con lo que el sistema queda
como:

y 0 = f (t, y),
y(t0 ) = y0 .
En caso de resolver el problema de forma analítica, los valores iniciales
permiten obtener las constantes de dicha función solución que la particulariza
para esas condiciones.
Si el problema se resuelve de forma numérica, dichas condiciones iniciales
son el punto de partida sobre el cual se van construyendo las soluciones apro-
ximadas en cada t. Así pues, conocemos las condiciones iniciales, y también
el tiempo tf al que queremos llegar. Subdividimos el intervalo [t0 , tf ] en N
t −t
subintervalos iguales, tomando tn = t0 + nh, n = 1, 2, ..., h = f N 0 .
Pasamos ahora a presentar las definiciones que clasifican los métodos
numéricos. Un método numérico puede ser:

Método de un paso: la aproximación a la solución exacta, yn+1 , se

calcula sólo a partir de la aproximación anterior, yn .

21
22 Cálculo de trayectorias de satélites

Método multipaso: la aproximación a la solución exacta, yn+1 , se obtie-

ne a partir de varias aproximaciones de pasos anteriores: yn+1 = ψh (t0 , y0 , ..., tn , yn ).

También se puede clasificar en:

Método explícito: la solución aproximada se calcula a partir de valores

conocidos con una ecuación explícita: yn+1 = ψh (tn , yn ).

Método implícito: la solución aproximada y(tn+1 ) se debe calcular re-

solviendo una ecuación implícita: y(tn+1 ) = ψh (tn , yn , tn+1 , yn+1 ).

De esta manera, ya tenemos la base sobre la que explicar los distintos

métodos numéricos estándar.

3.2. Métodos numéricos generales

Los métodos numéricos generales son aquellos que resuelven el sistema sin
considerar sus peculiaridades. Son los que se utilizan en software comercial
(como, por ejemplo, en Matlab) para resolver problemas de forma genérica.
Presentamos un método de alto error global, que se utiliza de forma más
académica que comercial, y dos métodos Adams multipaso, que representan
una mejor opción para la resolución de esos problemas.

3.2.1. Euler explícito

Es el método numérico más sencillo. Procede de considerar el desarrollo
de Taylor al primer orden:
∞
y (n) (t0 ) n
y(t0 + h) = h ≈ y(t0 ) + hf (t0 , y0 ).
X

n=0 n!
Así, se avanza para un n genérico:

yn+1 = yn + hf (tn , yn ), y0 = y(t0 ).

Si aplicamos el método para el problema del oscilador armónico, se tiene:

qn + h pmn 1
! ! ! ! !
h
qn+1 qn qn
= = m = (I + hA) ,
pn+1 pn − hkqn −hk 1 pn pn
donde I es la matriz identidad 2 × 2 y A es la matriz dada en (2.1). La matriz
(I + hA):
 Métodos numéricos 23

1
!
h
(I + hA) = m (3.1)
−hk 1
es el Jacobiano del método de Euler explícito respecto de q y de p.

Aplicando el método al problema de Kepler, se tiene:

! ! !
qn+1 qn pn
= +h .
pn+1 pn −µ qrn3
El error local de este método es de O(h2 ), con lo que al aumentar el
número de iteraciones (es decir, realizar cálculos prolongados en el tiempo),
el error numérico se dispara a un error global de O(h).
Para poder observarlo, realizamos la representación de los resultados de
la trayectoria obtenidos con este método y la solución real del problema no
perturbado de Kepler:

3
Euler explícito
Kepler
2

0
q2

−1

−2

−3

−4
−5 −4 −3 −2 −1 0 1 2
q1
Figura 3.1: Solución del método de Euler para el problema no perturbado de
Kepler, con los valores de (2.2). Se ha empleado un paso de tiempo h = 0. 05
para llegar a un instante tf = 100.
24 Cálculo de trayectorias de satélites

Se puede observar como, a partir de unas condiciones iniciales, este mé-

todo empieza a diverger de la solución real rápidamente. Para poder obtener
una aproximación realista, se requiere de un paso de tiempo h muy pequeño,
lo que aumenta mucho el coste computacional.

3.2.2. Método Adams-Bashforth de tres pasos

Los siguientes métodos son métodos multipaso que buscan evaluar la so-
lución aproximada utilizando varias soluciones calculadas anteriormente.
Si partimos del problema inicial:

y 0 = f (t, y), t0 ≤ t ≤ tn ,
integrando en el intervalo [tn , tn+1 ], se tiene
Z tn+1
y(tn+1 ) − y(tn ) = f (t, y(t)) dt.
tn

Ahora, supongamos que conocemos las aproximaciones al problema yj0 , yj1 , ..., yjk
de los tiempos tj0 , tj1 , ..., tjk , que se pueden calcular con otro método. Hace-
mos uso del polinomio de Lagrange Pk (t) que pasa por (tj0 , fj0 ), (tj1 , fj1 ), ..., (tjk , fjk ),
k
Pk (t) = Lk,l (t)fjl ,
X

l=0

siendo fn = f (tn , yn ). Lk,l (t) son los polinomios de Lagrange:

(t − t0 )(t − t1 ) · · · (t − tl−1 )(t − tl+1 ) · · · (t − tk )

Lk,l (t) = .
(tl − t0 )(tl − t1 ) · · · (tl − tl−1 )(tl − tl+1 ) · · · (tl − tk )

Se sustituye el problema por la aproximación:

Z tn+1 k Z tn+1
yn+1 = yn + Pk (t) dt = yn + Lk,l (t) dt.
X
f jl
tn l=0 tn

Llamando
Z tn+1
αk,l = Lk,l (t) dt,
tn

se obtienen los métodos multipaso o métodos Adams:

k
yn+1 = yn +
X
αk,l fjl .
l=0
 Métodos numéricos 25

Los métodos Adams explícitos se llaman Adams-Bashforth. Los méto-

dos implícitos incluyen el término (tn+1 , fn+1 ) para calcular el polinomio de
Lagrange, y se les conoce como métodos Adams-Moulton.
Considerando el método Adams-Bashforth de tres pasos se tiene:

h
yn+1 = yn + (23fn − 16fn−1 + 5fn−2 ),
12

con las condiciones iniciales

y0 = α0 , y1 = α1 , y2 = α2 .

Este método presenta un error local de O(h4 ).

Aplicando las ecuaciones para el caso del oscilador armónico, teniendo en
cuenta que:

!
qn
fn = A ,
pn

se tiene

! !  ! ! !
qn+1 qn h qn qn−1 qn−2
= + A 23 − 16 +5 .
pn+1 pn 12 pn pn−1 pn−2

Para el problema de Kepler, se tiene:

      
! !
qn+1 qn h   pn  pn−1 pn−2
= + 23 qn − 16  −µ qn−1  + 5  −µ qn−2 .

pn+1 pn 12 −µ r3
n
3
rn−1 3
rn−2

Con ello se consigue la siguiente figura para el problema de Kepler:

26 Cálculo de trayectorias de satélites

1
Adams Bashforth
0.8 Kepler

0.6

0.4

0.2
q2

−0.2

−0.4

−0.6

−0.8

−1
−1.5 −1 −0.5 0 0.5 1
q1
Figura 3.2: Solución del método de Adams-Bashforth para el problema no
perturbado de Kepler, con los valores de (2.2). Se ha empleado un paso de
tiempo h = 0. 05 para llegar a un instante tf = 100.

que refleja el aumento del orden del método, al disminuir la diferencia

entre el caso teórico y el computado.

3.2.3. Método de Adams-Moulton de tres pasos

Se trata también de un método multipaso de Runge-Kutta, aunque en
este caso el esquema se corresponde a una ecuación implícita y el error local
se reduce a O(h5 ). Con las condiciones iniciales igualmente definidas que
antes, se tiene:
h
(9fn+1 + 19fn − 5fn−1 + fn−2 ).
yn+1 = yn +
24
Con este método, la solución aproximada del oscilador armónico es:

! !  ! ! ! !
qn+1 qn h qn+1 qn qn−1 qn−2
= + A 9 + 19 −5 .
pn+1 pn 24 pn+1 pn pn−1 pn−2
 Métodos numéricos 27

Para el problema de Kepler, se tiene:

        
! !
qn+1 qn h pn+1 pn  pn−1 pn−2
= + 9  + 19  − 5  −µ qn−1  +  −µ qn−2 .
  
pn+1 pn 24 −µ rqn+1
3
qn
−µ r3 r3 r3
n+1 n n−1 n−2

Con este método obtenemos la mejor aproximación hasta ahora. Aumen-

tamos el paso de tiempo para obtener una deformación ”visible”:

1
Adams Moulton
0.8 Kepler

0.6

0.4

0.2
q2

−0.2

−0.4

−0.6

−0.8

−1
−1.5 −1 −0.5 0 0.5 1
q1
Figura 3.3: Solución del método de Adams-Moulton para el problema no
perturbado de Kepler, con los valores de (2.2). Se ha empleado un paso de
tiempo h = 0. 1 para llegar a un instante tf = 200.

Al ser un método implícito, el coste computacional de este método es

elevado ya que se tiene que resolver una ecuación implícita para cada paso
de tiempo.
 28 Cálculo de trayectorias de satélites

3.3. Comparación de los métodos

Es de esperar que los métodos de mayor orden obtengan una solución más
aproximada a la realidad. Definimos el error en posición como:

q
e(t) =(q1k − q1met )2 + (q2k − q2met )2 + (p1k − p1met )2 + (p2k − p2met )2 ,
(3.2)
donde el subíndice k hace referencia a la solución de Kepler, y met a la
solución computada por el método.

El error en el Hamiltoniano es:

eH (t) = |Hmet (t) − H0 | (3.3)

Comparamos los errores de los distintos métodos:

1 0

0 −1

−1 −2

−2 −3
log10(eH(t))
log10(e(t))

−3 −4

−4 −5

−5 −6

Euler explícito Euler explícito

−6 −7
Adams Moulton Adams Moulton
Adams Bashforth Adams Bashforth
Curva cuadrática Curva lineal
−7 −8
−1.5 −1 −0.5 0 0.5 1 1.5 2 −1.5 −1 −0.5 0 0.5 1 1.5 2
log10(t) log10(t)

Figura 3.4: Comparación logarítmica de los errores de los métodos, usando

la solución analítica del problema no perturbado de Kepler.
 Métodos numéricos 29

Se puede observar que el error en posición de estos métodos aumenta

alrededor de dos unidades logarítmicas en una unidad logarítmica de tiempo.
Ésto nos indica que el error en posición es proporcional al cuadrado de del
tiempo:

e(t) ∼ Ct2 .
Por ello, hemos dibujado dicha recta en la parte inferior izquierda, para
poder apreciar mejor dicha tendencia.
Se observa que el error en el Hamiltoniano crece de forma lineal con el
tiempo,

eH (t) ∼ Ct,
con lo que en este caso se ha dibujado la tendencia de una recta lineal.
Capítulo 4

Métodos simplécticos

4.1. Simplecticidad
Los métodos numéricos generales estudiados anteriormente no utiliza-
ban ninguna propiedad concreta de los sistemas Hamiltonianos que pudiese
mejorar su comportamiento. A base del estudio de lo que se conoce como
simplecticidad del flujo de los sistemas Hamiltonianos, se puede llegar a
construir métodos numéricos más eficientes.
El flujo ϕt (x0 ) de un sistema Hamiltoniano es una función <d → <d tal
que representa la solución de su sistema de ecuaciones diferenciales, para
un determinado t y unas condiciones iniciales x0 ∈ <d en t = 0. Definir la
solución de esta manera permite simplificar el estudio del sistema, ya que ϕt
viaja entre los mismos espacios vectoriales. Con ello, se realiza el siguiente
razonamiento [3]. Supongamos el espacio vectorial (q, p), q ∈ <d , p ∈ <d .
Definimos dos vectores ξ, η tales que:
! !
ξq ηq
ξ= , η=
ξp ηp
con ξ q , ξ p , η q , η p ∈ Rd . Estos vectores definen en el espacio (q, p) un paralelo-
gramo si consideramos:

P = {tξ + sη | 0 ≤ t ≤ 1, 0 ≤ s ≤ 1}.
Considerando el caso d = 1, se define el área orientada como:
!
ξ q ηq
[Link](P ) = det = ξ q ηp − ξ pηq .
ξ p ηp
Si d > 1, esta definición pasa a ser la suma de áreas orientadas de P
sobre los planos (qi , pi ):

31
32 Cálculo de trayectorias de satélites

d
ξiq ηiq d
!
ω(ξ, η) = det = (ξiq ηip − ξip ηiq ).
X X

i
ξip ηip i

La aplicación ω(ξ, η) es de vital importancia en sistemas Hamiltonianos.

Se trata de una aplicación bilineal sobre vectores en <2d . En notación matri-
cial, la aplicación se define como:

0 I
!
ω(ξ, η) = ξ Jη,
T
con J =
−I 0
El teorema de Poincaré (1899) demuestra que, considerando el flujo de
un sistema Hamiltoniano ϕt (x0 ), x0 = (q0 , p0 ), dicho flujo cumple que:

ϕ0t (x0 )T Jϕ0t (x0 ) = J.

Para poder apreciar mejor el significado de la expresión anterior, conside-
ramos una aplicación lineal A : <2d → <2d que también cumpla la ecuación
anterior:

AT JA = J
lo cual es equivalente a ω(Aξ, Aη) = ω(ξ, η) para cualquier ξ, η ∈ <2d . Ello
implica que A no modifica el área orientada definida anteriormente, lo
que se aprecia de forma más clara en la siguiente figura:

Figura 4.1: Conservación del área por una aplicación lineal, [3].

Por ello, para obtener una solución más realista, el método numérico
que emula la solución ϕt debe cumplir la misma propiedad que ésta. Una
aplicación diferenciable g : <2d → <2d se dice que es simpléctica si su
Jacobiano g 0 (q, p) cumple:

g 0 (q, p)T Jg 0 (q, p) = J

 Métodos simplécticos 33

en cuyo caso, conserva el área orientada. Dicha condición implica que |g 0 (q, p)| =
1.
Podemos observar que los métodos numéricos estándar definidos en el
capítulo anterior no cumplen esta propiedad; al no ser el determinante de su
Jacobiano igual a la unidad, deforman el área con cada paso de iteración. Para
el caso del oscilador armónico, el Jacobiano del método de Euler explícito se
obtuvo en (3.1):

1
!
h
(I + hA) = m .
−hk 1

Calculamos su determinante:

k
|I + hA| = 1 + h2 ,
m

que, como podemos observar, no cumple la condición de simplecticidad.

Para el problema de Kepler, con el método de Euler explícito se obtiene

un Jacobiano:

!
Id hId
N= .
− µh
r5
(Id r 2
− 3q q
n n
T
) Id

Donde Id es la matriz identidad de dimensión d × d. En un caso bidimen-

sional con q = (qx , qy ):

2h2 µ h4 µ2 3h2 µ 3h4 µ2

  !
 Id hId 
= 1 + + + − − ,

− µh (Id r2 − 3qn qnT ) Id  r3 r6 r3 r6
 
r5


con lo que tampoco se conserva el área en el problema de Kepler.

Para poder visualizarlo, decidimos calcular la solución del oscilador ar-
mónico con este método, con distintas condiciones iniciales, de forma que se
aprecie el área deformada [1]. Tomando k = m = 1, creamos una figura de
centro (q, p) = (3/2, 0) y r = 1/2 de radio, en los ángulos θ ∈ [−5π/6, 5π/6].
Aplicamos seis pasos h = π/6 y mostramos la solución en t = 0, t = π/2 y
t = π:
34 Cálculo de trayectorias de satélites

1
p

−1

−2

−3

−4
−4 −2 0 2 4
q
Figura 4.2: Área deformada por el método de Euler explícito para el oscila-
dor armónico ideal (h = π/6). Los círculos transparentes reflejan la solución
analítica, mientras que los oscuros se obtienen con el método de Euler explí-
cito.

Se aprecia que en este método en concreto, el área aumenta su valor.

4.2. Métodos simplécticos básicos

En base al estudio anterior, es interesante diseñar métodos numéricos que
sí sean simplécticos. Los primeros integradores simplécticos fueron propuestos
por Vogelaere (1956), Ruth (1983) y Feng Kang (1985) [3]. Presentamos ahora
los esquemas básicos de este tipo de métodos [1].
Supongamos que tenemos un Hamiltoniano H(q, p) separable en la suma
de energía cinética y energía potencial:

H(q, p) = T (p) + V (q).

Los sistemas Hamiltonianos cumplen las siguientes relaciones:
 Métodos simplécticos 35

ṗ = −∇q H(q, p) = −∇q V (q), q̇ = ∇p H(q, p) = ∇p T (p).

Como cada parte del Hamiltoniano depende de variables distintas, se
decide separar ambas variables, con lo que se tiene este sistema de ecuaciones:

q̇ = ∇p T (p) q̇ = 0
y
ṗ = 0 ṗ = −∇q V (q)

Al ser p y q constantes en el primer y segundo sistema de ecuaciones, res-

pectivamente, se pueden resolver fácilmente las ecuaciones de forma analítica,
con lo que se obtienen las siguientes soluciones:

q(t) = q0 + t∇p T (p0 ) q(t) = q0

ϕt :
[T ]
y
[V ]
ϕt :
p(t) = p0 p(t) = p0 − t∇q V (q0 )

[V ]
Ambas soluciones son simplécticas. Si se compone la aplicación ϕh desde
[T ]
las condiciones iniciales (qn , pn ) con ϕh , se obtiene el esquema:

χh ≡ ϕh ◦ ϕh : pn+1 = pn − h∇q V (qn ),

[T ] [V ]

qn+1 = qn + h∇p T (pn+1 ).

Éste es un método numérico de primer orden, que es simpléctico por ser

la composición de dos esquemas simplécticos [3]. Se le conoce como método
de Euler VT.
Si se intercambia el orden de la composición, se obtiene el método de
Euler TV:

χ∗h ≡ ϕh ◦ ϕh : qn+1 = qn + h∇p T (pn ),

[V ] [T ]

pn+1 = pn − h∇q V (qn+1 ).

El método Euler TV es el método adjunto del método Euler TV, que se

denota como:

χ∗h = χ−1
−h .

El adjunto de un método se puede obtener intercambiando h ↔ −h,

xn ↔ xn+1 , y permite aumentar el orden del método por composición.
36 Cálculo de trayectorias de satélites

Demostramos la simplecticidad del método de Euler VT para un problema

Hamiltoniano general. Para un problema unidimensional, derivamos respecto
de qn , pn para obtener el Jacobiano :

∂qn+1 ∂qn+1 ∂qn+1 ∂qn+1

= 1 + hHpq ; = h(Hpq + Hpp )
∂qn ∂qn ∂pn ∂pn
∂pn+1 ∂qn+1 ∂pn+1 ∂qn+1
= −hHpp ; = 1 − h(Hqq + Hqp ),
∂qn ∂qn ∂pn ∂pn
donde Hqq , Hqp , Hpp son las derivadas segundas parciales de los Hamiltonianos
evaluados en (qn+1 , pn ).
El Jacobiano del método queda como:

   
∂qn+1 ∂qn+1
1 + hHpq ∂q∂qn+1 h(Hpq ∂q∂pn+1 + Hpp )
N = ∂qn ∂pn = n n .
∂pn+1
∂qn
∂pn+1
∂pn
∂qn+1
−hHpp ∂qn 1 − h(Hqq ∂pn + Hqp )
∂qn+1

Finalmente, se observa que N cumple:

N T JN = J,
luego el método es simpléctico. Lo comprobamos con nuestros sistemas par-
ticulares.

Para el oscilador armónico simple, al ser un problema lineal, se puede

tratar qn , pn como si fueran escalares por no tener términos cruzados. Al
aplicar el método Euler VT se obtiene:

qn + h pn+1 1 − h2 m
! ! ! ! !
k h
qn+1 qn qn
= m = m = Mh .
pn+1 pn − hkqn −hk 1 pn pn

La matriz Mh representa el Jacobiano del sistema. Comprobamos que el

método sea simpléctico:

1 − h2 −hk 0 1 1 − h2 m
! ! !
k h
m =
m
h
1 −1 0 −hk 1
1 − h2 m 1 − h2 m
! !
h k k h
k m =
−1 h
m
−hk 1
0 1
!
,
−1 0
 Métodos simplécticos 37

con lo que el método es simpléctico. Efectivamente el determinante del Jaco-

biano es:

1 − h2 m
 
 k h 
|Mh | = m = 1.
 
−hk 1
 
 
Realizando el mismo procedimiento que en el apartado anterior, obser-
vamos que los métodos Euler VT y Euler TV deforman el área pero no
aumentan ni disminuyen su valor:

1
p

−1

−2

−3

−4
−4 −2 0 2 4
q
Figura 4.3: Deformación del área por los métodos Euler VT y TV para el
oscilador armónico (h = π/6). Los círculos transparentes reflejan la solución
analítica; los círculos ligeramente oscuros se obtienen con el método Euler
VT, y los oscuros con Euler TV.

Para el problema de Kepler, al ser no lineal, qn y pn se deben tratar como

vectores. Usando el método de Euler VT se obtiene:

qn
pn+1 = pn − hµ
rn3
qn+1 = qn + hpn+1 .
38 Cálculo de trayectorias de satélites

El Jacobiano de este método es:

 
2
Id − µh (Id r2 − 3qn qnT ) Id h 
N = r5 .
− r5 (Id r2 − 3qn qnT )
µh
Id

Se puede comprobar que cumple la simplecticidad:

 
µh2
!
 Id − r5
(Id r2 − 3qn qnT ) − µh
r5
(Id r2 − 3qn qnT )  Od Id
Id h Id −Id Od
   
2
µh2
Id − µh (Id r2 − 3qn qnT ) Id h   µh
(Id r2 − 3qn qnT ) Id − (Id r2 − 3qn qnT )
× r5 = r5 r5
− r5 (Id r2 − 3qn qnT )

µh
Id −Id Id h
 
2
(Id r2 − 3qn qnT )
!
Id − µh Id h  Od Id
× r5 = .
− µh
r5
(Id r2 − 3qn qnT ) Id −Id Od

y que el determinante de su Jacobiano es la unidad:

 
2
Id − µh (Id r2 − 3qn qnT )
 
In h 
= 1,

 r5

 − r5 (Id r2 − 3qn qnT )
µh
Id 

Si representamos el error como hicimos en el capítulo anterior, para el

problema no perturbado de Kepler, se puede apreciar su superioridad frente
al método de Euler explícito estándar:
 Métodos simplécticos 39

1 0
Euler explícito
Euler VT
Euler TV −1
0.5
Curva lineal

−2
0

−3
−0.5

log10(eH(t))
log10(e(t))

−4

−1

−5

−1.5
−6

−2
−7

−2.5
−8 Euler explícito
Euler VT
Euler TV
−3 −9
−1 −0.5 0 0.5 1 1.5 2 −1 −0.5 0 0.5 1 1.5 2
log10(t) log10(t)

Figura 4.4: Comparación logarítmica de los errores de los métodos, compa-

rados con la solución analítica del problema de Kepler. Valores de (2.2) con
un paso de tiempo h = 0. 05 y un instante tf = 100.

En el caso de los métodos simplécticos, el error en posiciones solo au-

menta de forma lineal con el tiempo, e(t) ∼ Ct, mientras que el error en
el Hamiltoniano no aumenta, sino que éste se encuentra acotado. Ésta es
una propiedad muy importante de los métodos simplécticos, que se conoce
como primera integral, que les otorga unas elevada fiabilidad para tiempos
de integración muy elevados.
La razón por la que, a pesar de ser métodos del mismo orden, se tiene
una diferencia en el error tan apreciable es lo que se estudia en el apartado
de backward error analysis.
Uno de los métodos simplécticos más empleados es el método de Störmer-
Verlet, el cual presenta un orden 2 al combinar el método de Euler VT con
el Euler TV de la siguiente manera:

[2] [V ] [T ] [V ]
Sh ≡ ϕh/2 ◦ ϕh ◦ ϕh/2 , (4.1)

lo que resulta en:

40 Cálculo de trayectorias de satélites

Figura 4.5: Backward error analysis, [3].

h
qn+1/2 = qn + pn ,
2
pn+1 = pn − h∇q V (qn+1/2 ),
h
qn+1 = qn+1/2 + pn+1 .
2
Este método se dice que es simétrico puesto que su adjunto es igual a
él mismo, χ∗h = χh , lo cual permite obtener métodos de ordenes muy eleva-
dos, como se verá posteriormente. También se puede componer el método de
Störmer-Verlet intercambiando T por V :
[2] [T ] [V ] [T ]
Sh ≡ ϕh/2 ◦ ϕh ◦ ϕh/2 .
La diferencia de comportamiento entre ambos métodos es mínima, y se
utilizará de ahora en adelante (4.1).

4.3. Backward error analysis

Una de las propiedades curiosas de los métodos numéricos es que, además
de resolver una ecuación diferencial con un cierto error, también represen-
tan la solución exacta de una ecuación diferencial ligeramente modificada
respecto de la que se intenta resolver. A este fenómeno se le conoce bajo
el nombre de backward error analysis. El documento [3] trata este asunto y
proporciona un método para obtener esa ecuación modificada.
Buscamos una ecuación diferencial del tipo:

ỹ˙ = f (ỹ) + hf2 (ỹ) + h2 f3 (ỹ) + ...,

 Métodos simplécticos 41

con la condición de que el método numérico ẏ = Φh (y), expandido como:

Φh (y) = y + hf (y) + h2 d2 (y) + h3 d3 (y) + ...,

sea la solución de dicha ecuación: ỹ(t + h) = Φh (y). Considerando una seme-
janza con la ecuación original y = ỹ(t) en un determinado punto t:

ỹ(t + h) =y + h(f (y) + hf2 (y) + h2 f3 (y) + ...)

h2 0
+ (f (y) + hf20 (y) + ...)(f (y) + hf2 (y) + ...) + ...,
2!
la igualdad nos proporciona las funciones f2 (y), f3 (y), ... necesarias para
obtener la ecuación diferencial modificada, en función del método numérico
empleado.
Los métodos numéricos simplécticos resuelven exactamente una ecuación
modificada que representa el Hamiltoniano real, más unas pequeñas contri-
buciones que aparecen en potencias de h mayores o iguales que el orden r del
método:

H̃ = H + hr Hr+1 + · · ·
Los métodos numéricos generales no son capaces de reproducir el Hamil-
toniano de esta manera. Ésta es la base de la superioridad de los métodos
simplécticos frente a métodos numéricos convencionales, y la razón por la
que se hace tanto hincapié en la propiedad simpléctica. Se puede demostrar
que, al hacer crecer exponencialmente el tiempo, el error que se observa en
el Hamiltoniano:

H(xn ) = H0 + O(hr ),
depende única y directamente del orden del método considerado, mientras
que en los métodos convencionales el error se acumula con cada paso de
tiempo.
En nuestro caso particular, ello supone que la trayectoria de satélites pue-
de ser calculada con precisión incluso en intervalos de tiempo muy amplios,
mientras que los métodos convencionales se quedarían obsoletos a partir de
un determinado instante donde el error acumulado se volvería intolerable.
Es una propiedad muy útil no solo de cara a la resolución numérica exacta
del problema, sino incluso para observar el comportamiento cualitativo de las
trayectorias a lo largo del tiempo: su estabilidad frente a diferentes posiciones
iniciales, la pérdida de energía por rozamiento,...
42 Cálculo de trayectorias de satélites

Se observa que la ecuación modificada cambiará según se cambie el paso h,

lo que hace que el método no pueda preservar las propiedades al ir cambiando
continuamente h en un esquema de paso variable. Ésta es la razón por la que
los métodos simplécticos deben emplearse con un paso h constante.
Capítulo 5

Aumento del orden

5.1. Estudio teórico

Existen dos técnicas popularmente utilizadas para aumentar el orden de
los métodos simplécticos descritos anteriormente: composición y splitting.
Además, para problemas Hamiltonianos ligeramente perturbados como
el nuestro, estas técnicas se pueden combinar de una manera especial, que
describimos posteriormente.

5.1.1. Composición
La técnica de composición consiste en, a partir de un método numérico
básico χh que cumple las propiedades simplécticas, pero no el orden requeri-
do, obtener un método ψh con la combinación del anterior de la forma:

ψh = χγs h ◦ χγs−1 h ◦ ... ◦ χγ1 h ,

donde γi son constantes elegidas para aumentar el orden del sistema.
La elección de las constantes no es arbitraria y su cálculo es tanto más
complejo cuanto mayor sea el orden que se requiera. Pero para aumentar el
orden del sistema siempre deben cumplir que, si s es el número de veces que
se combina el método básico, y r su orden,

γ1 + ... + γs = 1, γ1r+1 + ... + γsr+1 = 0.

De forma que el método global solo avanza un paso h.
Para calcular dichos coeficientes a partir de la exigencia de un determina-
do orden del método (order condition), se requieren ecuaciones adicionales.
Existen diversos procedimientos para calcularlas, pero no serán tratados en

43
44 Cálculo de trayectorias de satélites

este documento, ya que son en sí mismos trabajos de investigación avanzados.

Se puede encontrar más información al respecto en [1].

Composición de métodos simétricos

Los métodos simétricos de segundo orden pueden componerse, a su vez,
de forma simétrica, con lo que se puede obtener un método del orden par
que se desee de una forma sencilla. Considerando, por ejemplo, el método de
[2]
Störmer-Verlet Sh , se obtiene un método de orden 4 con:

[4] [2] [2] [2]1

Sh = Sαh ◦ Sβh ◦ Sαh , β = 1 − 2α.
con α =
2 − 21/3
Se puede generalizar el método a un orden general 2k + 2 realizando la
misma combinación:

[2k+2] [2k] [2k] [2k] 1

Sh = Sαh ◦ Sβh ◦ Sαh con α = , β = 1 − 2α.
2− 21/(2k+1)
A pesar de ello, esta forma de componer métodos es altamente ineficiente,
[4]
debido a los grandes pasos que se toman. Para Sh , se tienen α = 1. 35 . . .,
β = −1. 7 . . ., lo que lleva a pasos αh, βh muy grandes, que aumentan el
error local del método, ya que éste siempre es menor cuanto menor sea el
paso h.
Una manera más eficiente de construir métodos consiste en componer el
método básico χh con su adjunto χ∗h con diferentes pasos de tiempo:

ψh = χα2s h ◦ χ∗α2s−1 h ◦ · · · ◦ χα2 h ◦ χ∗α1 h (5.1)

eligiendo los valores reales apropiados (α1 , . . . , α2s ). Si se tiene que

α2s+1−i = αi , i = 1, 2, . . . , n
entonces el nuevo método ψh será simétrico. Los métodos simétricos son los
que más se están utilizando en la actualidad.
Compondremos más tarde un método con esta filosofía para observar su
ventaja frente a los métodos numéricos generales.

5.1.2. Splitting
En esta técnica, si recordamos nuestro campo de vectores:

ẋ = f (x), x(0) = x0 ∈ <d ,

 Aumento del orden 45

se subdivide f (x) en distintas funciones, que son las distintas contribuciones

al sistema: f (x) = f [1] (x) + f [2] (x) + ... + f [m] (x) de tal manera que cada
subproblema

ẋ = f [i] (x), x(0) = x0 ∈ <d i = 1, 2, ...m

[i]
sí pueda ser integrable con las soluciones exactas x(t) = ϕt . Se pueden
combinar las soluciones analíticas de cada subproblema, evaluadas en t = h,
de la siguiente manera usando composición:
[m] [2] [1]
χh = ϕh ◦ ... ◦ ϕh ◦ ϕh ,
para obtener una aproximación de muy elevada precisión. Se puede seguir
aumentando el orden utilizando composición con el adjunto. Teniendo en
[i] [i]
cuenta que (ϕh )−1 = ϕ−h , el adjunto de χh se puede calcular como:

[m] [2] [1] [1] [m−1] [m]

−h = (ϕ−h ◦ ... ◦ ϕ−h ◦ ϕ−h )
χ∗h = χ−1 −1
= ϕh ◦ ... ◦ ϕh ◦ ϕh

Este procedimiento es especialmente interesante en nuestro caso. Nues-

tro problema consiste en analizar la trayectoria del satélite, cuyo caso ideal
ya está resuelto, teniendo en cuenta distintas perturbaciones que hacen in-
calculable la solución de forma analítica. Con el splitting se pueden juntar
distintas soluciones numéricas, para proporcionar una primera aproximación
de por sí muy exacta, y después aumentar el orden y con ello la fiabilidad
del método.
Considerando el conjunto de perturbaciones, nuestro problema se formula
como:

∂Vεa (q)
!
q r−a
q̈ = −µ 3 − εa − εr exp − ||q̇|| q̇.
r ∂q b
Utilizando la nomenclatura Hamiltoniana se puede expresar el problema
de forma matricial:
   
0  0 
! !
q̇ p
= + εa   + εr 
−µ rq3 − ∂Vε∂qa (q) − exp − r−a

ṗ b
||p|| p

donde Vεa (q) es la perturbación del Hamiltoniano por el achatamiento, que

sólo depende de q.
Aún así, debe recordarse que cuanto más orden se requiera, mayor número
de veces se tendrán que combinar los métodos, lo que conlleva a reevaluar
las funciones más veces, aumentando el coste computacional.
46 Cálculo de trayectorias de satélites

Métodos Runge-Kutta-Nyström
Si consideramos un campo de vectores separable en 2 (m = 2):

f (x) = f [1] (x) + f [2] (x),

[2] [1]
podemos construir el esquema básico χh = ϕh ◦ ϕh con sus soluciones ana-
[i] [1] [2]
líticas ϕh , i = 1, 2. El método adjunto será χ∗h = ϕh ◦ ϕh , con lo que al
componerlos como en (5.1), se obtiene:

       
[2] [1] [1] [2] [2] [1] [1] [2]
ψh = ϕα2s h ◦ ϕα2s h ◦ ϕα2s−1 h ◦ ϕα2s−1 h ◦· · ·◦ ϕ α2 h ◦ ϕα2 h ◦ ϕα1 h ◦ ϕα1 h .

[i] [i]
Al ser ambas ϕh las soluciones exactas, cumplen la propiedad 1 ϕβh ◦ϕδh =
[i]
ϕ(β+δ)h . Ésto permite reducir el número de soluciones involucradas al juntar
[1] [2]
las soluciones ϕh , ϕh colindantes entre sí mismas, con lo que nos queda el
esquema:

[2] [1] [2] [2] [1] [2]

ψh = ϕbs+1 h ◦ ϕas h ◦ ϕbs h ◦ · · · ◦ ϕb2 h ◦ ϕa1 h ◦ ϕb1 h . (5.2)
Se puede aplicar el esquema anterior al problema no perturbado de Kepler
para comprobar la eficacia de estos métodos. Considerando nuestra ecuación
diferencial:

0
! ! !
q̇ p
ẋ = = + = f [1] (x) + f [2] (x)
ṗ 0 −µ rq3
[1]
las soluciones analíticas de las ecuaciones ẋ = f [1] (x) y ẋ = f [2] (x) son ϕh y
[2]
ϕh respectivamente. Al componerlas como se ha descrito anteriormente, se
obtienen los llamados métodos de Runge-Kutta-Nyström, N Bs[r] . Estos
métodos han demostrado ser especialmente eficaces para resolver sistemas
Hamiltonianos con energía cinética cuadrática, de la forma:

1
H(q, p) = pT M −1 p + V (q).
2
Para obtener un método de orden 4, se requiere s = 6, con los siguientes
coeficientes [12]:
1
La composición de una solución analítica consigo misma, en dos pasos de tiempo
distintos, equivale a avanzar la suma de esos pasos con la propia solución analítica. No es
válido con soluciones numéricas.
 Aumento del orden 47

[4]
Coeficientes para la composición N B6 dada por (5.2) con s = 6
b1 = 0. 0829844064174052 b2 = 0. 3963098014983681 b3 = −0. 039056304922348
b4 = 1 − 2(b1 + b2 + b3 ) b5 = b3 , b6 = b2 , b7 = b1 a1 = 0. 2452989571842710
a2 = 0. 6048726657110800 a3 = 1/2 − (a1 + a2 ) a4 = a3 ; a5 = a2 ; a6 = a1

Finalmente podemos comparar este método con los métodos Adams-

Moulton y Adams-Bashforth de orden similar que hemos descrito en el capí-
tulo correspondiente:

2 0
[4] [4]
N B6 N B6
Adams Moulton Adams Moulton
Adams Bashforth Adams Bashforth
0 −2

−2 −4
log10(eH(t))
log10(e(t))

−4 −6

−6 −8

−8 −10

−10 −12
−1.5 −1 −0.5 0 0.5 1 1.5 2 −1.5 −1 −0.5 0 0.5 1 1.5 2
log10(t) log10(t)

Figura 5.1: Comparación de métodos para el problema no perturbado de

Kepler con los valores (2.2), con h = 0. 05.

Se observa la gran mejoría que ofrece el método simpléctico. Hay que tener
en cuenta que el método N B6 realiza 6 evaluaciones del sistema, con lo que
es más costoso computacionalmente que el resto de métodos; sin embargo,
este coste se compensa con el bajo error que obtiene, lo que se demostrará
en apartados posteriores.
48 Cálculo de trayectorias de satélites

5.1.3. Sistemas casi-integrables

La situación de los sistemas casi-integrables es idéntica a la nuestra; al
considerar un sistema Hamiltoniano, se tiene una ligera perturbación del
mismo por efectos no ideales

H = H0 + εH1 ,
con 0 < ε  1, lo que transforma el sistema en

ẋ = f (x) = f [1] (x) + εf [2] (x),

y permite emplear el método de splitting para obtener una buena aproxi-
mación al sistema. Al hacerlo, en este caso especial, el error se ve bastante
reducido, ya que en él participa ε, que siempre es menor de 1. En esta clase
de métodos, el error viene dado por:

hE(h, ε) = O(εhr1 +1 + ε2 hr2 +1 + ... + εm hrm +1 ),

y se dice que el método tiene un orden generalizado de (r1 , r2 , ..., rm ), siendo
r1 el orden que más define al método, pues ε → 0. Este error es mucho más
pequeño que el producido por un método numérico general, puesto que éste
depende de εh, mucho menor que h.
Así pues, para cada uno de los casos, obtenemos la solución de la per-
turbación que aparece en nuestra ecuación diferencial, para componer un
método por splitting y compararlo con un método estándar de primer orden.

Achatamiento
Considerando sólo el primer término relevante, la ecuación del Hamilto-
niano perturbado por el achatamiento de un planeta es:

1 1 3q 2
!
ε
H = (p2x + p2y ) − µ − 3 1 − α 2x .
2 r 2r r
El problema sigue siendo Hamiltoniano, con lo que el Hamiltoniano per-
turbado deberá conservarse a lo largo del tiempo.
La ecuación diferencial perturbada se obtiene del Hamiltoniano pertur-
bado:
 
εa 3q 2
H(q, p)εa = − 3 1 − α ε2a ,x  = εa Va (q).
2rεa rεa
 Aumento del orden 49

Al ser el Hamiltoniano perturbado asimétrico para qεa ,x y qεa ,y , debemos

derivar respecto de cada variable para obtener las distintas ecuaciones:

d
q̇εa = H(q, p)εa = 0
dpεa
d 3qεa ,x (−qε2a ,y (1 + 2α) + qε2a ,x (−1 + 3α))
ṗεa ,x = − H(q, p)εa = εa
dqεa ,x 2rε7a
d 3qεa ,y (qε2a ,y + qε2a ,x (1 − 5α))
ṗεa ,y = − H(q, p)εa = −εa .
dqεa ,y 2rε7a

Puesto que en la perturbación las posiciones qεa se mantienen constantes

con el tiempo, ṗεa no depende del tiempo y se puede obtener la siguiente
solución analítica:

qεa (t) = qεa ,0 ,

3qεa ,x,0 (−qε2a ,y,0 (1 + 2α) + qε2a ,x,0 (−1 + 3α))
pεa ,x (t) = pεa ,x,0 + εa t,
2r7
3qεa ,y,0 (qε2a ,y,0 + qε2a ,x,0 (1 − 5α))
pεa ,y (t) = pεa ,y,0 − εa t.
2r7

[a] [k]
Denotando por ϕh la solución anterior en un paso h, y ϕh como la
solución analítica al problema de Kepler no perturbado con ese mismo paso,
tenemos el esquema:

[a] [k]
χh = ϕh ◦ ϕh ,

Tomamos el valor εa = 0. 001 y calculamos el error en el Hamiltoniano

como (3.3), comparándolo con el del método de Euler explícito:
50 Cálculo de trayectorias de satélites

6
Euler VT
Euler explícito
4

2
log10(eH(t))

−2

−4

−6

−8

−10
−1.5 −1 −0.5 0 0.5 1 1.5 2
log10(t)
Figura 5.2: Error en Hamiltonianos debido al achatamiento, con los valores
(2.2) y con εa = 0. 001, α = 1 y h = 0. 05.

Se puede comprobar que el método Euler VT mantiene acotado el error en

el Hamiltoniano perturbado, mientras que el explícito lo aumenta y alcanza
valores muy elevados.

Queremos observar el efecto del achatamiento sobre la trayectoria del sa-

télite. La calculamos tanto con el método de Euler VT, como con el ode45
de Matlab, con una tolerancia absoluta de 10−10 y una tolerancia relativa de
10−12 , para asegurarnos de que obtenemos un resultado correcto. La compa-
ración entre ambos métodos se realizará más adelante.
 Aumento del orden 51

Euler VT ode45

1 1

0.8 0.8

0.6 0.6

0.4 0.4

0.2 0.2
q2

q2
0 0

−0.2 −0.2

−0.4 −0.4

−0.6 −0.6

−0.8 −0.8

−1 −1
−1.5 −1 −0.5 0 0.5 1 −1.5 −1 −0.5 0 0.5 1
q1 q1

Figura 5.3: Efecto del achatamiento sobre el problema de Kepler, con εa =

0. 001, h = 0. 025, y tf = 500.

Para esta representación hemos decidido aumentar el tiempo de simula-

ción. Se observa que el efecto del achatamiento es girar la elipse a lo largo
del tiempo, aunque ésta se sigue manteniendo estable.

Rozamiento con la atmósfera

La ecuación diferencial de la parte perturbada se obtuvo matricialmente

como:
 
0 
!
q̇εr
= εr   ,
ṗεr − exp − rεrb−a ||pεr || pεr

lo que implica que la energía del sistema se irá perdiendo a lo largo del
tiempo. El valor de εr se puede asumir como bajo por los valores habituales
de CD, m2 ,... Con las condiciones iniciales de nuestro problema se tiene como
solución:
52 Cálculo de trayectorias de satélites

!
pεr ,0 rε ,0 − a
q(t)εr = qεr ,0 , p(t)εr = q , C = εr exp − r .
||pεr ,0 || (1/||pεr ,0 || + C t)2 b

[r] [k]
Denotando por ϕh la solución anterior en un paso h, y ϕh como la
solución analítica al problema de Kepler no perturbado con ese mismo paso,
tenemos el esquema de primer orden:
[r] [k]
χh = ϕh ◦ ϕh . (5.3)
Considerando una perturbación εr = 0. 001, decidimos calcular el error
en posición del método como (3.2) y compararlo con el método de Euler
explícito. Para ello, definimos como solución exacta a los resultados obtenidos
por (5.3) con un paso 10 veces menor que el empleado en dichos métodos,
con lo que finalmente se obtiene:

1
Euler VT
Euler Explícito
0

−1
log10(error(t))

−2

−3

−4

−5

−6
−1.5 −1 −0.5 0 0.5 1 1.5 2
log10(t)
Figura 5.4: Errores en posiciones para el problema de Kepler con rozamien-
to, con h = 0. 05, a = 0, b = 1 y con los valores (2.2).

Se observa la clara superioridad del método simpléctico frente al general.

 Aumento del orden 53

Observamos el efecto del rozamiento:

Euler VT ode45
1 1

0.8 0.8

0.6 0.6

0.4 0.4

0.2 0.2
q2

q2
0 0

−0.2 −0.2

−0.4 −0.4

−0.6 −0.6

−0.8 −0.8

−1 −1
−1.2 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8
q1 q1

Figura 5.5: Trayectoria para el problema de rozamiento, con εr = 0. 001,

a = 0, b = 1, h = 0. 05 y con los valores (2.2).

Efectivamente se comprueba que el rozamiento va disminuyendo la altura

del satélite.

5.2. Aplicaciones
Veremos ahora las aplicaciones de los conceptos presentados en los apar-
tados anteriores.

5.2.1. Método NIA

Nuestro sistema casi-integrable presenta la forma:

ẋ = f [k] (x) + εr f [r] (x) + εa f [a] (x) = f [k] (x) + f [ε] (x)
que es una estructura semejante a la estudiada en los métodos Runge-Kutta-
Nyström. Por ello, si ϕkh es la solución analítica de Kepler, y ϕεh un método
54 Cálculo de trayectorias de satélites

que resuelve las perturbaciones, se pueden componer con el mismo esquema

que dichos métodos. En este caso el error se ve reducido a εh y se obtiene los
métodos para sistemas casi-integrables, N IA:
[k] [ε] [k] [k] [ε] [ε]
N IA[s
s
1 ,s2 ,...]
= ϕas+1 h ◦ ϕbs h ◦ ϕas h ◦ · · · ◦ ϕa2 h ◦ ϕb1 h ◦ ϕa1 h .
Para obtener un orden generalizado (4, 2), se requiere s = 2, con lo que
se tiene:
[4,2] [k] [ε] [k] [ε] [k]
N IA2 = ϕa3 h ◦ ϕb2 h ◦ ϕa2 h ◦ ϕb1 h ◦ ϕa1 h
√
con a1 = a3 = (3 − 3)/6, a2 = 1 − 2a1 , b1 = b2 = 1/2. Este método es
especialmente adecuado ya que cumple que bi > 0, con lo que se garantiza
que el sistema vaya perdiendo energía con las perturbaciones. En métodos
con un orden más alto se puede llegar a tener bi < 0, con lo que se requiere
un tratamiento cuidadoso de los mismos.
[ε]
Como solución ϕh , decidimos componer las perturbaciones con el método
de Störmer-Verlet:
[ε] [a] [r] [a]
ϕ̃h = ϕh/2 ◦ ϕh ◦ ϕh/2
y si tuviésemos n perturbaciones a considerar, se podría considerar igual-
mente:

[ε] [ε ] [ε ] [ε ] [ε ] [ε ] [ε ] [ε ]
ϕ̃h = ϕh/2
1 2
◦ ϕh/2 n−1
◦ · · · ◦ ϕh/2 ◦ ϕh n ◦ ϕh/2
n−1 2
◦ · · · ◦ ϕh/2 1
◦ ϕh/2 . (5.4)

Componer las soluciones de las perturbaciones de esta manera no propor-

ciona la solución exacta del conjunto de perturbaciones, sino una aproxima-
ción tal que:
[ε] [ε]
ϕh = ϕ̃h + O(h2 εa εr )
que se puede considerar válida por el bajo valor de h2 εa εr .
Esta forma de tratar sistemas ligeramente perturbados genera métodos
muchos más eficaces que los ofrecidos en el software de cálculo general. En
[4,2]
el siguiente apartado se va a comparar el método N IA2 con un algoritmo
de orden similar, el ode45 de Matlab.

5.2.2. Ode45
El algoritmo de Matlab ode45 es uno de los métodos numéricos de bajo
orden más utilizados, al ser el primero que no requiere de condiciones especí-
ficas de la ecuación a resolver (matriz de masas constante, matriz de masas
 Aumento del orden 55

singular,...). Se basa en un método explícito de Runge-Kutta (4, 5) de paso

variable.

[4,2]
Decidimos comparar la eficiencia del método N IA2 con el algoritmo
ode45. Definimos la eficiencia de un método como su capacidad para alcanzar
un determinado orden de error con el menor coste computacional posible. El
coste computacional depende tanto del compilador como de la arquitectura
del ordenador, así como de la forma de implementación del programa. Por
ello, es habitual definir el coste de un método como el número de veces que
la parte más costosa del algoritmo debe ser ejecutada; en nuestro caso, dicha
parte es el campo vectorial f (t, y). Entonces, si el método debe evaluar n
veces la función f (t, y), se tiene que coste = n. Si calculamos el error en
posición como (3.2), entonces podemos comparar el algoritmo ode45 con el
[4,2]
método N IA2 . En este caso, solo tomaremos el error del valor final obtenido
(e(tf )), al que le corresponderá un determinado coste computacional; cuanto
menor sea el error exigido, mayor será el coste computacional.

Con ode45, obtenemos la solución final imponiendo tolerancias absolutas

de 10−i y tolerancias relativas de 10−i−2 , con i = 4, . . . , 13, y calculamos el
coste computacional asociado a cada una de ellas.

[4,2]
Para el método N IA2 , se define un número de iteraciones creciente
N = 10 · 2j , j = 1, . . . , 10, con lo que el paso de tiempo se irá reduciendo
t −t
con h = f N 0 , siendo t0 , tf y el resto de parámetros iguales a los valores
[4,2] [k]
de (2.2). El método N IA2 realiza 3 computaciones del término ϕh y 2
[ε]
computaciones del término ϕh por paso de tiempo, con lo que en el peor
caso su coste será coste (j) = 3 · N (j).

De esta manera, se obtiene el siguiente resultado:

56 Cálculo de trayectorias de satélites

−1

−2
log10(e(tf ))

−3

−4

−5

−6

−7 [4,2]
N IA2
ode45
−8
2 2.5 3 3.5 4
log10(coste)
[4,2]
Figura 5.6: Comparación del algoritmo ode45 con el método N IA2 .

[4,2]
La diferencia entre ambos métodos es muy notable; el método N IA2 es
capaz de alcanzar una solución más exacta con un menor coste computacio-
nal que el ode45. Esta diferencia se haría aún más grande si se integrase el
problema en tiempos finales mayores.
De esta manera, con el enfoque apropiado, se puede resolver mejor el
problema perturbado de Kepler con un método numérico ”casero” que con
un método comercial muy optimizado.

5.3. Sistema caótico

El estudio de soluciones numéricas precisas es especialmente importante
en mecánica celeste, debido a que la presencia de términos de perturbación
lleva a resultados muy distintos para condiciones iniciales ligeramente distin-
tas. Es lo que se conoce como caos.

Para ilustrarlo, vamos a realizar el siguiente experimento. Supongamos

las siguientes condiciones iniciales:
 Aumento del orden 57

1+e
s
q1 (0) = 1 − e, q2 (0) = 0, p1 (0) = 0, p2 (0) = ,
1−e (5.5)
t0 = 0, tf = 10, e = 0. 2, µ = 1.

Ahora supongamos que somos incapaces de saber si la posición exacta del

satélite se sitúa en el punto (1 − e, 0), o en algún punto situado en un radio
r = 1/100 alrededor de éste, con lo que tenemos el conjunto de posiciones
iniciales:

q1 (0) = (1 − e) + r cos(θ), q2 (0) = r sin(θ),

2π (5.6)
θi = i, j = 0, 1, . . . N, N = 50.
N

Decidimos, en primer lugar, integrar nuestro sistema con un ode45 muy

preciso, de tolerancia absoluta 10−10 y tolerancia relativa 10−12 . De esta forma
conseguimos la solución exacta del problema en tf = 10.

En segundo lugar, se integra usando ode45 con menor precisión (tolerancia

absoluta 10−2 y tolerancia relativa 10−4 ), de forma que se obtiene una solución
menos precisa, con un coste computacional menor que el anterior.

Se estima el número de iteraciones N que se requiere para calcular el

[4,2]
sistema con el método N IA2 y el mismo coste computacional que el ode45
poco preciso, N = costeode45 /3, y se calcula la solución con ese número de
iteraciones. El conjunto de soluciones que se obtiene es:
58 Cálculo de trayectorias de satélites

0.1

−0.1

−0.2

−0.3
q2

−0.4

−0.5

−0.6

−0.7

−0.8

−0.9
−1.5 −1 −0.5 0 0.5 1
q1
Figura 5.7: Conjunto de soluciones de posición para el problema perturbado
de Kepler considerando rozamiento y achatamiento, obtenidas con ode45
preciso (área poco oscura), ode45 poco preciso (área ligeramente oscura) y
[4,2] [4,2]
N IA2 (área oscura). Los resultados del ode45 preciso del método N IA2
se superponen.

Se extraen dos conclusiones de este resultado.

Por una parte, la falta de precisión del ode45 de elevada tolerancia penali-
za su solución de tal manera que ninguna coincide con las soluciones exactas
[4,2]
del problema, obtenidas con el método N IA2 y el ode45 preciso. Ello es
un indicativo de la necesidad de tener métodos de elevado orden de cara a
obtener resultados fiables.
Por otra parte, las soluciones exactas muestran un área fuertemente alar-
gada en comparación con el área circular de las condiciones iniciales. Se
produce por la elevada diferencia entre las soluciones exactas de cada punto
de posición inicial. Cuando se obtienen soluciones de este estilo, se dice que
el sistema sea caótico; pequeñas diferencias iniciales llevan a soluciones fi-
nales muy distintas. El origen de este caos son las fuerzas que actúan sobre
el sistema, que son muy sensibles a la posición y a la velocidad de satélite.
Se puede comprobar el mismo efecto para variaciones en la cantidad de
 Aumento del orden 59

movimiento. Considerando las condiciones iniciales (5.5), pero esta vez, con
diferentes cantidades de movimiento iniciales:

1+e
s
p1 (0) = r cos(θ), p2 (0) = + r sin(θ), θ ∈ [0, 2π],
1−e

con r = 1/100, y realizando el mismo procedimiento anterior, se obtienen las

siguientes soluciones de cantidad de movimiento:

1.5

0.5
p2

−0.5

−1
−0.2 0 0.2 0.4 0.6 0.8 1 1.2
p1
Figura 5.8: Conjunto de soluciones de cantidad de movimiento para el pro-
blema perturbado de Kepler considerando rozamiento y achatamiento, obte-
nidas con ode45 preciso (área blanca), ode45 poco preciso (área ligeramente
[4,2]
oscura) y N IA2 (área oscura) .

Se observan los mismos efectos que en el estudio de las soluciones de las

posiciones.
Capítulo 6

Métodos actuales

Tras haber visto todas las propiedades y ventajas de los métodos simpléc-
ticos, no es de extrañar que éstos sean los métodos que se utilizan actual-
mente para el cálculo de trayectorias en Mecánica Celeste. En este capítulo
vamos a describir esos métodos, que serán el producto final de todo el estudio
presentado sobre la integración geométrica.

6.1. Método ABAH844

En el documento [13] se presenta el desarrollo completo de los métodos
específicos de orden elevado para cálculos en Mecánica Celeste. De todos
estos métodos, estimamos que el método ABAH844 es el más adecuado de
presentar en este trabajo.
Sigue el mismo esquema (5.2), siendo específicamente de la forma:

[k] [ε] [k] [ε] [k] [ε] [k] [ε] [k] [ε] [k] [ε] [k]
ψh = ϕa1 h ◦ ϕ̃b1 h ◦ϕa2 h ◦ ϕ̃b2 h ◦ϕa3 h ◦ ϕ̃b3 h ◦ϕa4 h ◦ ϕ̃b3 h ◦ϕa3 h ◦ ϕ̃b2 h ◦ϕa2 h ◦ ϕ̃b1 h ◦ϕa1 h
[ε]
donde ϕ̃h es un método simétrico de segundo orden que resuelve las pertur-
baciones, como el definido en (5.4). Se ha diseñado específicamente con esa
condición para alcanzar un orden generalizado (8, 4) con estos coeficientes:

61
62 Cálculo de trayectorias de satélites

Coeficientes para el método ABAH844 dado por (5.2) con s = 6

a1 = 0. 2741402689434018761640565440378637101205
a2 = −0. 1075684384401642306251105297063236526845
a3 = −0. 04801850259060169269119541715084750653701
a4 = 0. 7628933441747280943044988056386148982021
b1 = 0. 6408857951625127177322491164716010349386
b2 = −0. 8585754489567828565881283246356000103664
b3 = 0. 7176896537942701388558792081639989754277
Con ésto, ya se pueden realizar estudios muy avanzados sobre trayectorias
reales de satélites.
Para comprobar su eficiencia, una vez más, comparamos su comporta-
miento al resolver el problema de Kepler perturbado con achatamiento y
[4,2]
rozamiento, con el ode45 y el N IA2 . En este caso, el coste del método
ABAH844 es costeABAH844 (j) = 7 · N (j). Reproduciendo las mismas condi-
ciones que en el capítulo anterior, obtenemos:

−1

−2

−3
log10(e(tf ))

−4

−5

−6

−7 ode45
[4,2]
N IA2
ABAH844
−8
2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 4
log10(coste)
Figura 6.1: Comparación del error en posición frente al coste computacional
[4,2]
de los métodos ode45, N IA2 y ABAH844 para el problema perturbado de
Kepler con achatamiento y rozamiento.
 Métodos actuales 63

[4,2]
La razón por la que se superponen los métodos N IA2 y ABAH844
en un punto es que el coste de ABAH844 es muy superior al del método
N IA. Así, para aproximaciones de bajo orden, se puede emplear el método
[4,2]
N IA2 , y para cálculos de elevada precisión, ABAH844.

6.2. Evolución del sistema

Ahora que se dispone de un método plenamente fiable para el cálculo
de la trayectoria del satélite, decidimos mostrar la evolución del sistema en
posiciones con las condiciones (5.5) modificados como (5.6), calculada con el
método ABAH844 para una vuelta tf = 2π. Mostramos los resultados para
t = i, i = 1, . . . , 5.
Procediendo de la misma manera que en apartados anteriores, acabamos
obteniendo:

0.5

0
q2

−0.5

−1

−1.5
−1.5 −1 −0.5 0 0.5 1
q1
Figura 6.2: Evolución de un satélite para el problema perturbado de Kepler
con achatamiento y rozamiento, calculada por el método ABAH844 para
distintas posiciones iniciales. Se han tomado instantes t = i, i = 1, . . . , 5,
con tf = 2π.
64 Cálculo de trayectorias de satélites

Se observa que conforme el satélite va avanzando, la diferencia entre las

condiciones iniciales se va ensanchando, como resultado del comportamiento
caótico del sistema que ya se explicó anteriormente.

6.3. Conclusiones
Una vez hemos llegado a este punto, nos queda evaluar los resultados que
hemos obtenido del trabajo, y cómo éstos han cumplido nuestro objetivo:

Presentamos el problema de Kepler, su forma de perturbarse, las di-

ficultades que presenta para su solución analítica, y particularizamos
nuestro estudio para dos perturbaciones, el achatamiento del planeta y
el rozamiento con la atmósfera, con tal de ofrecer un ejemplo concreto
de cómo tratar las perturbaciones.

Justificamos por qué los métodos numéricos generales no son aptos para
resolver nuestro problema, detallando sus inconvenientes.

Introducimos el concepto de simplecticidad, la propiedad de los siste-

mas Hamiltonianos que permite estimar el comportamiento de la solu-
ción real. Reflejamos que los métodos numéricos generales no cumplen
esta propiedad, lo que motivó la presentación de los métodos simpléc-
ticos. Estudiamos las ventajas de sus características.

Detallamos la manera de construir métodos simplécticos más eficaces

a través de las técnicas de composición y splitting tanto para distintos
casos como para nuestro problema concreto. Se seleccionó el método
[4,2]
N IA2 para comparar su eficiencia con el ode45 de Matlab, que cons-
tituye un buen ejemplo de herramienta para resolver numéricamente
ecuaciones diferenciales.

Finalmente, se aporta el método ABAH844 como ejemplo de méto-

do actual para resolver el problema perturbado de Kepler de forma
eficiente, con lo que cumplimos el objetivo de nuestro trabajo.

En la parte final de nuestro documento, construimos los métodos numé-

ricos como funciones de Matlab para tener un legado físico y práctico
del trabajo.
 Apéndice A

Software

Una vez presentadas las ventajas de los métodos simplécticos frente a

los métodos generales, recogemos aquí los códigos elaborados en Matlab con
dichos métodos, de forma que puedan ser utilizados por un usuario general.
Presentamos las funciones necesarias para resolver el problema no per-
turbado y el problema perturbado de Kepler, así como sus generalizaciones
para su uso en problemas diferentes.

A.1. Kepler no perturbado

La solución analítica de la parte V del problema no perturbado de Kepler,
[V ]
ϕh , viene dada por el script KeplerV.m:

1 f u n c t i o n [ q p]= KeplerV ( q0 , p0 , h )
2 % R e s u e l v e a n a l i t i c a m e n t e para un paso h
3 % l a p a r t e V d e l problema no p e r t u r b a d o de
4 % Ke pl er .
5 q=q0 ;
6 r 3 2 =(q0 ( 1 ) ^2+q0 ( 2 ) ^2) ^ ( 3 / 2 ) ;
7 p=p0−h∗q . / r 3 2 ;
8 end

De la misma manera, la parte T viene dada por KeplerT.m:

1 f u n c t i o n [ q p]= KeplerT ( q0 , p0 , h )
2 % R e s u e l v e a n a l i t i c a m e n t e para un paso h l a p a r t e T
3 % d e l problema no p e r t u r b a d o de Kepler .

65
 66 Cálculo de trayectorias de satélites

4 q=q0+h∗p0 ;
5 p=p0 ;
6 end

La función NB6.m produce la composición de ambas funciones para

obtener un método Runge-Kutta-Nyström de orden 4:

1 f u n c t i o n [ q p]=NB6( q0 , p0 , h )
2 % D e f i n i c i o n de c o e f i c i e n t e s .
3

4 b1 = 0 . 0 8 2 9 8 4 4 0 6 4 1 7 4 0 5 2 ; b2 = 0 . 3 9 6 3 0 9 8 0 1 4 9 8 3 6 8 1 ; b3
= −0.039056304922348;
5 b4=1−2∗(b1+b2+b3 ) ; b5=b3 ; b6=b2 ; b7=b1 ; a1
=0.2452989571842710;
6 a2 = 0 . 6 0 4 8 7 2 6 6 5 7 1 1 0 8 0 0 ; a3=1/2−(a1+a2 ) ; a4=a3 ; a5=a2 ; a6
=a1 ;
7

8 % Composicion .
9 [ q p]= KeplerV ( q0 , p0 , b1∗h ) ;
10 [ q p]= KeplerT ( q , p , a1 ∗h ) ;
11 [ q p]= KeplerV ( q , p , b2∗h ) ;
12 [ q p]= KeplerT ( q , p , a2 ∗h ) ;
13 [ q p]= KeplerV ( q , p , b3∗h ) ;
14 [ q p]= KeplerT ( q , p , a3 ∗h ) ;
15 [ q p]= KeplerV ( q , p , b4∗h ) ;
16 [ q p]= KeplerT ( q , p , a4 ∗h ) ;
17 [ q p]= KeplerV ( q , p , b5∗h ) ;
18 [ q p]= KeplerT ( q , p , a5 ∗h ) ;
19 [ q p]= KeplerV ( q , p , b6∗h ) ;
20 [ q p]= KeplerT ( q , p , a6 ∗h ) ;
21 [ q p]= KeplerV ( q , p , b7∗h ) ;
22

23 end

Para generalizar el método, se define NB6_General.m, que hace uso de

una lista de funciones indefinidas f , que el usuario deberá introducir como
solución de las partes que quiera componer:

1 f u n c t i o n [ q p]=NB6_General ( q0 , p0 , h , f )
 Software 67

2 % f e s un ’ ’ c e l l a r r a y o f f u n c t i o n handle ’ ’ que
3 % a l b e r g a dos f u n c i o n e s s o b r e l a s c u a l e s a p l i c a r
4 % e l a l g o r i t m o NB6 .
5

6 % D e f i n i c i o n de c o e f i c i e n t e s .
7

8 b1 = 0 . 0 8 2 9 8 4 4 0 6 4 1 7 4 0 5 2 ; b2 = 0 . 3 9 6 3 0 9 8 0 1 4 9 8 3 6 8 1 ; b3
= −0.039056304922348;
9 b4=1−2∗(b1+b2+b3 ) ; b5=b3 ; b6=b2 ; b7=b1 ; a1
=0.2452989571842710;
10 a2 = 0 . 6 0 4 8 7 2 6 6 5 7 1 1 0 8 0 0 ; a3=1/2−(a1+a2 ) ; a4=a3 ; a5=a2 ; a6
=a1 ;
11

12 % D e f i n i c i o n de f u n c i o n e s .
13 f u n c i o n 1=f { 1 } ;
14 f u n c i o n 2=f { 2 } ;
15

16 % Composicion .
17 [ q p]= f u n c i o n 2 ( q0 , p0 , b1∗h ) ;
18 [ q p]= f u n c i o n 1 ( q , p , a1 ∗h ) ;
19 [ q p]= f u n c i o n 2 ( q , p , b2∗h ) ;
20 [ q p]= f u n c i o n 1 ( q , p , a2 ∗h ) ;
21 [ q p]= f u n c i o n 2 ( q , p , b3∗h ) ;
22 [ q p]= f u n c i o n 1 ( q , p , a3 ∗h ) ;
23 [ q p]= f u n c i o n 2 ( q , p , b4∗h ) ;
24 [ q p]= f u n c i o n 1 ( q , p , a4 ∗h ) ;
25 [ q p]= f u n c i o n 2 ( q , p , b5∗h ) ;
26 [ q p]= f u n c i o n 1 ( q , p , a5 ∗h ) ;
27 [ q p]= f u n c i o n 2 ( q , p , b6∗h ) ;
28 [ q p]= f u n c i o n 1 ( q , p , a6 ∗h ) ;
29 [ q p]= f u n c i o n 2 ( q , p , b7∗h ) ;
30

31 end

Finalmente, en la página web [Link] se

puede encontrar una función que soluciona analíticamente el problema no
perturbado de Kepler, Kepler.m:

1 f u n c t i o n [ q p ] = Kepler ( q0 , p0 , h )
2 % Esta f u n c i o n s e puede o b t e n e r en
 68 Cálculo de trayectorias de satélites

3 % h ttp : / /www. g i c a s . u j i . e s /GNIBook . html

4 t o l =1. e −10;
5 r 0=norm ( q0 ) ; v02=dot ( p0 , p0 ) ; u=dot ( q0 , p0 ) ;
6 a=−1/(v02 −2/ r 0 ) ; En=1/ s q r t ( a∗a∗a ) ;
7 Ec=1−r 0 /a ; Es=u / (En∗a∗a ) ;
8

9 % Primer v a l o r para e l c a l c u l o de l a e c u a c i o n
implicita .
10 x=En∗h∗a/ r 0 ; x1=x ;
11

12 % C a l c u l o de l a e c u a c i o n i m p l i c i t a .
13 while 1
14 c=c o s ( x ) ; s=s i n ( x ) ;
15 x=x−(x−Ec∗ s+Es∗(1− c )−En∗h ) /(1−Ec∗ c+Es∗ s ) ;
16 i f abs ( x−x1 ) < t o l , break , end
17 x1=x ;
18 end
19 c=c o s ( x ) ; s=s i n ( x ) ;
20 FP=1 − Ec∗ c+Es∗ s ;
21 f f =(c −1)∗a/ r 0 +1. d0 ; g=h+(s−x ) /En ;
22 f f p=−a∗En∗ s / (FP∗ r 0 ) ; gp=(c −1)/FP+1;
23

24 f o r i =1:2
25 q ( i )= f f ∗ q0 ( i )+g∗p0 ( i ) ;
26 p ( i )=f f p ∗ q0 ( i )+gp∗p0 ( i ) ;
27 end

A.2. Kepler perturbado

La función Achatamiento.m proporciona la solución analítica a la per-
turbación producida por el achatamiento del planeta, con unos valores prefi-
jados de εa = 0. 001 y α = 1:

1 f u n c t i o n [ q p v a r a r g o u t ]= Achatamiento ( q0 , p0 , h , v a r a r g i n )
2 e p s i l o n =10^−3;
3 alpha =1;
4

5 % El achatamiento a f e c t a s o l o a p .
6 q=q0 ;
 Software 69

8 % La f o r m u l a d e l achatamiento no e s s i m e t r i c a y hay que

d e f i n i r l a por
9 % partes .
10 x=q0 ( 1 ) ;
11 y=q0 ( 2 ) ;
12

13 % En c a s o de que ya s e haya i n t r o d u c i d o l a r ,
aprovecharla .
14

15 switch nargin
16 case 4
17 r=v a r a r g i n ( 1 ) ;
18 r=c e l l 2 m a t ( r ) ;
19 otherwise
20 r =(x^2+y ^2) ^ ( 1 / 2 ) ;
21 v a r a r g o u t={r } ;
22 end
23 p ( 1 )=p0 ( 1 )+h∗3∗ x∗ e p s i l o n ∗(−y^2∗(1+2∗ alpha )+x^2∗(−1+3∗
a lp ha ) ) /(2 ∗ r ^ ( 7 ) ) ;
24 p ( 2 )=p0 ( 2 )−h∗3∗ y∗ e p s i l o n ∗ ( y^2+x^2∗(1 −5∗ alpha ) ) /(2 ∗ r ^ ( 7 )
);
25 end

En este caso hemos decidido introducir como entrada o salida opcional el

valor de r = ||q|| que se requiere para la función. Puesto que el achatamiento
y el rozamiento no alteran a q en sus ecuaciones, se puede utilizar el cálculo
de r en una de las llamadas de esta función para reintroducirlo
q en la siguiente
y reducir el coste computacional del método, puesto que q T q contiene un
elevado número de operaciones.
El rozamiento no requiere del cálculo de r, con lo que solo lo requeriría
como input opcional en la función Rozamiento.m:

1 f u n c t i o n [ q p]= Rozamiento ( q0 , p0 , h , v a r a r g i n )
2 e p s i l o n =10^−3;
3 a =0;
4 b=1;
5

6 % Solucion .
7
 70 Cálculo de trayectorias de satélites

8 switch nargin
9 case 4
10 r 0=v a r a r g i n ( 1 ) ;
11 r 0=c e l l 2 m a t ( r 0 ) ;
12 otherwise
13 r 0=norm ( q0 ) ;
14 end
15

16 q=q0 ;
17 c t e=e p s i l o n ∗ exp (−( r0−a ) /b ) ;
18 p=p0 / ( norm ( p0 ) ∗ s q r t ( ( 1 / norm ( p0 )+c t e ∗h ) ^2) ) ;
19 end

[ε]
Para calcular la aplicación ϕh , se compone ambas soluciones con el mé-
todo de Störmer-Verlet en Perturbaciones.m:

1 f u n c t i o n [ q p]= P e r t u r b a c i o n e s ( q0 , p0 , h )
2 % Composicion d i r e c t a de l a s p e r t u r b a c i o n e s
Achatamiento y Rozamiento ,
3 % con e l metodo de Stormer−V e r l e t .
4 [ q p r 0 ]= Achatamiento ( q0 , p0 , h / 2) ;
5 [ q p]= Rozamiento ( q , p , h , r 0 ) ;
6 [ q p]= Achatamiento ( q , p , h / 2 , r 0 ) ;
7 end

Ello permite componer la función NIA.m:

1 f u n c t i o n [ q p]=NIA( q0 , p0 , h )
2 % Composicion para e l s i s t e m a c a s i −i n t e g r a b l e .
3

4 % D e f i n i c i o n de l o s c o e f i c i e n t e s :
5 a1 =(3 −3^(1/2) ) / 6 ;
6 b1 =1/2;
7 a2=1−2∗a1 ;
8

9 % Composicion de l a s f u n c i o n e s .
10 [ q p]= Kepler ( q0 , p0 , a1 ∗h ) ;
11 [ q p]= P e r t u r b a c i o n e s ( q , p , b1∗h ) ;
12 [ q p]= Kepler ( q , p , a2 ∗h ) ;
 Software 71

13 [ q p]= P e r t u r b a c i o n e s ( q , p , b1∗h ) ;
14 [ q p]= Kepler ( q , p , a1 ∗h ) ;
15

16 q=q ’ ; p=p ’ ;
17 end

Si se quiere una aproximación más precisa, se puede utilizar ABAH844.m:

1 f u n c t i o n [ q p]=ABAH844( q0 , p0 , h )
2 % C o e f i c i e n t e s para e l metodo .
3

4 a1 = 0.2741402689434018761640565440378637101205;
5 a2 = −0.1075684384401642306251105297063236526845;
6 a3 = −0.04801850259060169269119541715084750653701;
7 a4 = 0.7628933441747280943044988056386148982021;
8 b1 = 0.6408857951625127177322491164716010349386;
9 b2 = −0.8585754489567828565881283246356000103664;
10 b3 = 0.7176896537942701388558792081639989754277;
11

12 % Composicion d e l metodo
13 [ q p]= Kepler ( q0 , p0 , a1 ∗h ) ;
14 [ q p]= P e r t u r b a c i o n e s ( q , p , b1∗h ) ;
15 [ q p]= Kepler ( q , p , a2 ∗h ) ;
16 [ q p]= P e r t u r b a c i o n e s ( q , p , b2∗h ) ;
17 [ q p]= Kepler ( q , p , a3 ∗h ) ;
18 [ q p]= P e r t u r b a c i o n e s ( q , p , b3∗h ) ;
19 [ q p]= Kepler ( q , p , a4 ∗h ) ;
20 [ q p]= P e r t u r b a c i o n e s ( q , p , b3∗h ) ;
21 [ q p]= Kepler ( q , p , a3 ∗h ) ;
22 [ q p]= P e r t u r b a c i o n e s ( q , p , b2∗h ) ;
23 [ q p]= Kepler ( q , p , a2 ∗h ) ;
24 [ q p]= P e r t u r b a c i o n e s ( q , p , b1∗h ) ;
25 [ q p]= Kepler ( q , p , a1 ∗h ) ;
26

27 q=q ’ ; p=p ’ ;
28 end

En caso de que se quieran abordar un mayor número de perturbaciones,

el usuario puede definir una lista de funciones indefinidas f e introducirlas
en las siguientes funciones generales:
 72 Cálculo de trayectorias de satélites

Perturbaciones_General.m

1 f u n c t i o n [ q p]= P e r t u r b a c i o n e s _ G e n e r a l ( q0 , p0 , h , f )
2 % f e s un ’ ’ c e l l a r r a y o f f u c n t i o n handles ’ ’ ,
3 % de forma que s e o b t i e n e e l metodo Stormer−V e r l e t para
4 % distintas perturbaciones a considerar .
5

6 % Primera p a r t e .
7 f o r i =1: l e n g t h ( f )−1
8 f u n c i o n=f { i } ;
9 [ q p]= f u n c i o n ( q0 , p0 , h / 2) ;
10 q0=q ; p0=p ;
11 end
12

13 % Funcion i n t e r m e d i a .
14 f u n c i o n=f { l e n g t h ( f ) } ;
15 [ q p]= f u n c i o n ( q0 , p0 , h ) ;
16 q0=q ; p0=p ;
17

18 % Segunda p a r t e .
19 f o r i =1: l e n g t h ( f )−1
20 f u n c i o n=f { l e n g t h ( f )−i } ; % Composicion de f u n c i o n e s
en s e n t i d o c o n t r a r i o .
21 [ q p]= f u n c i o n ( q0 , p0 , h / 2) ;
22 q0=q ; p0=p ;
23 end
24

25 end

NIA_General.m

1 f u n c t i o n [ q p]=NIA_General ( q0 , p0 , h , f )
2 % f e s un ’ ’ c e l l a r r a y o f f u n c t i o n handle ’ ’
3 % para i n s e r t a r l a s p e r t u r b a c i o n e s que s e d e s e e .
4

5 % D e f i n i c i o n de l o s c o e f i c i e n t e s :
6 a1 =(3 −3^(1/2) ) / 6 ;
7 b1 =1/2;
8 a2=1−2∗a1 ;
9
 Software 73

10 % Composicion de l a s f u n c i o n e s .
11 [ q p]= Kepler ( q0 , p0 , a1 ∗h ) ;
12 [ q p]= P e r t u r b a c i o n e s _ G e n e r a l ( q , p , b1∗h , f ) ;
13 [ q p]= Kepler ( q , p , a2 ∗h ) ;
14 [ q p]= P e r t u r b a c i o n e s _ G e n e r a l ( q , p , b1∗h , f ) ;
15 [ q p]= Kepler ( q , p , a1 ∗h ) ;
16

17 q=q ’ ; p=p ’ ;
18 end

ABAH844_General.m

1 f u n c t i o n [ q p]=ABAH844_General ( q0 , p0 , h )
2 % C o e f i c i e n t e s para e l metodo .
3

4 a1 = 0.2741402689434018761640565440378637101205;
5 a2 = −0.1075684384401642306251105297063236526845;
6 a3 = −0.04801850259060169269119541715084750653701;
7 a4 = 0.7628933441747280943044988056386148982021;
8 b1 = 0.6408857951625127177322491164716010349386;
9 b2 = −0.8585754489567828565881283246356000103664;
10 b3 = 0.7176896537942701388558792081639989754277;
11

12 % Composicion d e l metodo
13 [ q p]= Kepler ( q0 , p0 , a1 ∗h ) ;
14 [ q p]= P e r t u r b a c i o n e s _ G e n e r a l ( q , p , b1∗h ) ;
15 [ q p]= Kepler ( q , p , a2 ∗h ) ;
16 [ q p]= P e r t u r b a c i o n e s _ G e n e r a l ( q , p , b2∗h ) ;
17 [ q p]= Kepler ( q , p , a3 ∗h ) ;
18 [ q p]= P e r t u r b a c i o n e s _ G e n e r a l ( q , p , b3∗h ) ;
19 [ q p]= Kepler ( q , p , a4 ∗h ) ;
20 [ q p]= P e r t u r b a c i o n e s _ G e n e r a l ( q , p , b3∗h ) ;
21 [ q p]= Kepler ( q , p , a3 ∗h ) ;
22 [ q p]= P e r t u r b a c i o n e s _ G e n e r a l ( q , p , b2∗h ) ;
23 [ q p]= Kepler ( q , p , a2 ∗h ) ;
24 [ q p]= P e r t u r b a c i o n e s _ G e n e r a l ( q , p , b1∗h ) ;
25 [ q p]= Kepler ( q , p , a1 ∗h ) ;
26

27 q=q ’ ; p=p ’ ;
28 end
Bibliografía

[1] S. Blanes and [Link], A Concise Introduction to Geometric Numerical

Integration, Taylor & Francis Group, 2016.

[2] S. Blanes, D. Ginestar and M.D. Roselló, Introducción a los métodos nu-
méricos para ecuaciones diferenciales , Editorial Universitat Politècnica
de València, 2014.

[3] E. Hairer, C. Lubich and G. Wanner, Geometric Numerical Integration,

Springer, 2006.

[4] J.M. Sanz-Serna, Geometric Integration, The State of the Art in Nume-
rical Analysis (York, 1996) (New York), Inst. Math. Appl. Conf. Ser.
New Ser., vol. 63, Oxford Univ. Press, 1997, pp. 121 − 143.

[5] V. G. Szebehly and H. Mark, Adventures in Celestial Mechanics, 2nd

ed., Die Deutsche Bibliothek, 2004.

[6] J. Lequeux, Le Verrier - Magnificent and Detestable Astronomer, Sprin-

ger, 2013.

[7] J.C. Butcher, A history of Runge-Kutta methods, Elsevier, 1996.

[8] V. Martinusi, L. Dell Elce and G. Kerschen, First order analytic propa-
gation of satellites in the exponential atmosphere of and oblate planete,
Celest Mech Dyn Astr (2017) 127 : 451 − 476.

[9] J.M.A. Danby, Fundamentals of Celestial Mechanics, Willmann-Bell,

1988.

[10] F.R. Moulton, Introduction to Celestial Mechanics, The MacMillan

Company, 1960.

[11] J. Laskar, P. Robutel, F. Joutel, M. Gastineau, A.C.M. Correia and B.

Levrard, A long-term numerical solution for the insolation quantities of
the Earth, A&A 428, 261 − 285 (2004).

75
76 Cálculo de trayectorias de satélites

[12] S. Blanes and P.C. Moan, Practical symplectic partitioned Runge-Kutta

and Runge-Kutta-Nyström methods, J. Comput. Appl. Math. 142 (2002),
313 − 330.

[13] S. Blanes, F. Casas, A. Farrés, J. Laskar, J. Makazaga, and A. Murua,

New families of symplectic splitting methods for numerical integration
in dynamical astronomy. Appl. Numer. Math. 68 (2013), pp. 58 − 72.