- 1 -

Capítulo 6 – Métodos Numéricos de Integração de Órbita

Os métodos de propagação de órbita podem ser classificados em: analítico,

semi-analítico, numérico, e de multi-revolução.

O método analítico é o método mais rápido de propagação de órbita; porém,

oferece poucas opções que possam torná-lo mais genérico. Neste método a derivação

da teoria envolve o uso de aproximações, expansões em séries e integração analítica.

Algumas características para o método analítico são: manipulação algébrica laboriosa

na dedução; modelos simples de forças perturbadoras; fácil visualização dos efeitos

perturbadores, o que permite a compreensão da tendência de evolução dos elementos

orbitais; pouco gasto de tempo computacional; e a órbita propagada num único passo

indiferente ao intervalo de propagação.

O método semi-analítico tenta fazer uso da rapidez computacional da teoria

analítica em conjunto com a precisão e flexibilidade de modelagem no método

numérico. Ele é baseado normalmente na formulação analítica em conjunto com a

assistência de métodos numéricos. Usualmente separa-se o movimento orbital de

longo período e secular pelo método das médias, para perturbações de origem

conservativa; e realiza-se média numérica para outras perturbações, normalmente

dissipativas, através de quadratura.

O método numérico, também chamado de método de perturbação especial, tem

a capacidade de lidar com qualquer tipo de situação com uma precisão arbitrária,

utilizando procedimentos de cálculo passo a passo. Este método é bastante utilizado

quando se requer precisão ou o número de revoluções da órbita a ser integrada é

limitado, porque sua maior desvantagem reside na lentidão do método devido à

sequência de cálculos passo a passo. Basicamente tais métodos tentam fazer

aproximações polinomiais da trajetória e, portanto, integram exatamente um

polinômio de certo grau a menos de erros de discretização.

Os métodos de multi-revolução essencialmente possuem a mesma forma dos

métodos numéricos clássicos; entretanto, devido à natureza periódica do movimento

- 2 -

orbital, consegue- se extrapolar muitas revoluções no processo de integração. Estes

métodos são baseados em aproximações polinomiais e, consequentemente, funções

polinomiais de variável independente discreta são integradas exatamente, a menos dos

erros de truncamento/arredondamento. Este método melhora a precisão para certos

problemas de órbita de satélites.

Existem duas metodologias consagradas entre o método numéricos: o método de

Cowell e o método de Encke que serão descritos brevemente a seguir.

6.1 - Método de Cowell

Uma das grandes vantagens deste método é a rapidez com que a técnica pode

ser usada. A sua vantagem é a simplicidade de formulação e implementação. Consiste

em escrever as equações do movimento do objeto a ser estudado, incluindo todas as

perturbações e então integrá-las passo-a-passo numericamente. As perturbações podem

ser incluídas todas ao mesmo tempo. Este método não requer espaço de

armazenamento nos computadores em relação aos outros métodos.

Entretanto, quando o movimento estiver próximo do corpo central, passos de

integração menores devem ser tomados o que afeta o tempo e erro acumulado do

processo. Quando as perturbações se tomam mais severas, exigem-se passos de

integração menores e mais cálculos, o que implica a provável ocorrência de maior erro

acumulado de discretização. Este método é várias vezes mais lento que o método de

Encke.

Para o problema de 2-corpos com perturbações, as equações seriam:

Pr

r =+ r3

µ&& (6.1)

onde r e v são os vetores posição e velocidade de um satélite com respeito ao corpo

central, µ a constante geo-gravitacional, e P o vetor de acelerações perturbadoras.

Para a integração numérica, a equação (1) pode ser reduzida a equações

diferenciais de 1a ordem:

- 3 -

rPv

vr

3r

µ−=

=

&

&

(6.2)

onde r e v são os vetores posição e velocidade de um satélite com respeito ao corpo

central. Para a equação de integração numérica, elas podem ser expressadas nas

componentes dos vetores:

−=

−=

−=

=

=

=

zr

Pv

yr

Pv

xr

Pv

vz

vy

vx

zz

yy

xx

z

y

x

3

3

3

,

µ

µ

µ

&

&

&

&

&

&

(6.3)

onde 222zyxr ++= .

6.2 - Método de Encke

Mais complexo do que o método de Cowell, o método de Encke realiza uma

integração de Equações Diferenciais Ordinárias (EDO) que é a diferença entre a

aceleração do problema dos 2-corpos, ou da referência adotada, e a aceleração de todas

as perturbações. O método tenta integrar as acelerações menores e, por isso, retem mais

algarismos significativos, enquanto cresce o tamanho do passo de integração.

Esquematicamente as EDOs ficariam na forma:

3

com ,

r

rrr

Pr

µ+≡∆

=∆

&&&&

&&

A vantagem do método de Encke reside no fato de que, como as perturbações

são em geral pequenas, r&&∆ e suas derivadas resultam igualmente pequenas e, em

consequência, é possível o uso de um intervalo (passo) de integração mais amplo do que

no método de Cowell. Nas vizinhanças do pericentro onde os intervalos de integração

têm de ser diminuídos quando se emprega o método de Cowell, o método de Encke

permanece invariável. Porém, cada etapa do método de Encke gasta mais tempo. Da

- 4 -

mesma forma, quando as perturbações se acumulam, r&&∆ contínua crescendo e as

equações podem tomar-se pouco manejáveis, caso em que se deve encontrar um novo

conjunto de elementos mais atualizados, com os quais se deve dar prosseguimento à

integração. Este processo é chamado de retificação. Logo, outro problema é a escolha

do(s) instante(s) de retificação, isto é, quando um novo r&&∆ deve ser inicializado para

recomeçar a integração.

6.3 - Métodos para resolução de Equações diferenciais ordinárias

O método de integração numérica é um dos mais poderosos métodos

conhecidos em mecânica celeste para calcular o movimento de qualquer corpo no

sistema solar em torno do corpo primário.

Os métodos de integração numérica podem ser divididos em métodos de passo

fixo, passos múltiplos (ou multi-passos), e extrapolação. Os métodos de passo fixo não

necessitam informações além das condições iniciais e das derivadas. Os métodos de

passos múltiplos exigem um procedimento de inicialização. Alguns dos métodos de

passo fixo conhecidos são Runge-Kutta, Gill, Euler-Cauchy e Bowie. Alguns dos

métodos de multi-passo são Milne, Adams-Moulton, Gauss-Jackson, Obrechkoff e

Adams-Bashforth.

As equações a serem integradas são Equações Diferenciais Ordinárias (EDOs)

de 1a ordem, com a condição inicial xo em to, ou seja,

),( xfx

tdt

d=

onde x e f são vetores, e f é uma função vetorial de t e x.

6.3.1 - Método de passo fixo “Runge-Kutta”

A técnica aproxima a função de derivadas em série de Taylor para cálculos das

primeiras derivadas em pontos dentro do intervalo de extrapolação. A ordem de um

- 5 -

membro particular da família é a ordem da potência mais alta do tamanho do passo, h,

na expansão da série de Taylor equivalente. A forma genérica do RK é dada por:

∑=

+ +=R

i

iinn w1

1 kxx

onde wi são pesos, R o número de avaliações, de f, e ki são funções dadas por:

++= ∑

−

=

1

1

,i

j

jijnii ahcth kxfk n

com Ric ,,2,1,01 K== . Logo:

( )( )( )

M

23213133

12122

1

,

,

,

kkxfk

kxfk

xfk

aahcth

ahcth

th

nn

nn

nn

+++=

++=

=

Para obter os valores para os parâmetros desconhecidos, expande-se 1+nx em potências

de h de forma a concordar com a solução das EDs até uma ordem específica da série de

Taylor. Assim, nota-se que o número de equações (não-lineares) de condição é menor

que o número de parâmetros desconhecidos. Esta é a razão pela qual existem diferentes

parâmetros para um Runge-Kutta de mesma ordem. Nota-se que os RKs fornecem a

EDO integrada em um passo, sem conhecimento prévio do tipo do problema, e com

avaliações sucessivas da função f dentro do passo de integração. O número de

avaliações depende da ordem do método; em geral, ordens maiores implicam maior

precisão às custas de maior tempo de processamento. Dentre os RKs utilizados em

problemas de mecânica orbital destacam-se os de Fehlberg (1968) e Shanks (1966).

Por exemplo, a formulação para o Runge-Kutta padrão de ordem 4 (RK4) é:

( )43211 226

kkkkxx ++++=+

hnn (6.4)

onde

- 6 -

( )

( )3

2

1

22

22

kxfk

kxfk

kxfk

xfk

4

3

2

1

h,ht

h,h

t

h,h

t

,t

nn

nn

nn

nn

++=

++=

++=

=

O método de Runge-Kutta é estável e não exige um procedimento de

inicialização. Ele é relativamente simples, fácil de implementar, produz um pequeno

erro de truncamento, e o tamanho do passo é fácil de ser mudado. Uma das

desvantagens é que não existe um caminho simples para determinar o erro de

truncamento, logo, é tarefa difícil determinar o tamanho do passo adequado. Uma

característica óbvia é que as informações usadas se referem somente ao passo de

integração atual. A idéia de usar informações de passos anteriores conduz ao método

de multi-passos ou preditor-corretor.

6.3.2 - Método de multi-passos

Os métodos de multi-passos usam informações de intervalos xk, xk-l, xk-2, ... para

calcular xk+1 em tk + h. A interpolação polinomial de Lagrange é utilizada para ajustar

os pontos aos valores correspondentes da função f, e o polinômio resultante é integrado

exatamente, a menos de erros de truncamento e arredondamento. Em geral, eles são

mais rápidos do que o método de passo fixo, embora à custa de maior complexidade e

uma necessidade de um procedimento de inicialização.

6.3.2.1 – Fórmulas de Adams-Bashforth

Alguns métodos calculam um valor predito para 1nx + e então substituem 1nx +

na equação diferencial para obter 1nx +& que é por sua vez usado para calcular uma

valor corrigido de 1nx + . Estes são naturalmente chamados de métodos do tipo

- 7 -

"preditor-corretor". As fórmulas de Adams são as mais utilizadas tradicionalmente em

métodos de multi-passos, principalmente quando a avaliação da função de derivadas f

é dispendiosa. Por exemplo, a fórmula de Adams-Bashforth, para integração de

EDOs, de ordem N em xn usando polinômios de grau k-1, interpolando via Lagrange a

função f calculado nos k pontos precedentes, é:

∑=

+ ∇+=N

k

n

k

knn h0

1 fxx α (6.5)

onde N é o número de termos desejados, h é o tamanho do passo, n é o número do passo

atual, α = 1, 1/2, 5/12, 3/8, 251/720, 95/288, 5,0 K=k , e

),(

0

111

nnn

nn

n

k

n

k

n

k

t xff

ff

fff

≡

≡∇

∇−∇≡∇ −−−

Os poucos primeiros passos são dados por:

nnn h fxx

∇+∇+∇+∇+∇++=+

54321 288

95

720

251

8

3

12

5

2

11

O método permite diminuir ou dobrar o tamanho do passo sem uma reinicialização

completa. Também o erro de truncamento pode ser avaliado.

6.3.2.2 – Fórmulas de Adams-Moulton

As fórmulas de Adams-Bashforth são um esquema de predição do método de

multipassos. Moulton adicionou uma fórmula de correção para o método. As fórmulas

de ordem k em xk usam o polinômio de grau k-1 , que interpola até a k-ésima derivada.

A fórmula em geral é dada por:

- 8 -

∑=

−++ ∇+=N

k

kn

k

knn h0

1*

1 fxx α (6.6)

Nota-se que é necessário uma estimativa de 1+nx para poder calcular o

polinômio interpolante, uma vez que 1+nf é também função de 1+nx . Poder-se-ia então

iterar 1+nx sucessivamente até se obter a convergência desta equação não-linear. A

maior complexidade computacional é tolerada porque estas fórmulas levam a

resultados mais precisos devido ao uso de um polinômio interpolante mais adequado.

6.3.2.3 – Preditor-Corretor

Uma fórmula que fornece a primeira aproximação para 1+nx é chamada de

preditor. Se o preditor é suficientemente acurado, pode ser desnecessário aplicar a

fórmula de correção de Moulton mais de uma vez. Por essa razão a fórmula de Adams-

Bashforth é recomendada para a composição de um preditor-corretor em conexão com

o método de Adams-Moulton. Neste caso, a fórmula para o preditor é

( )5321

)(

720

2519375955

24 dt

xdh

hnnnnn

ξ55p

1n ffffxx +−+−+= −−−+ (6.7)

e para o corretor é

( )5211

)(

720

195199

24 dt

xdh

hnnnnn

ξ55c

1n ffffxx ++−++= −−++ (6.8)

onde o último termo em cada uma das equações é o termo do erro truncamento. Na

equação (6.8), 1+nf é encontrado através do valor predito p

n 1+x :

( )p

nnn t 111 , +++ = xff

Para usar qualquer esquema de preditor-corretor, o valor predito é calculado

primeiro. Então, a derivada correspondente para o valor predito é encontrado e usado

para encontrar o valor corrigido usando a fórmula do corretor. É possível iterar a

- 9 -

fórmula do corretor até que não exista nenhuma mudança significativa em 1+nx nas

iterações sucessivas. O tamanho do passo pode ser mudado se o erro de truncamento

for maior do que o desejado. Este método, como outros métodos de multi-passos, exige

um método de passo fixo para iniciá-lo para obter K,, 1−nn ff etc. O método de Runge-

Kutta é normalmente sugerido como o inicializador. Desse modo, para métodos de

multi-passos de ordem N, requer-se uma tabela de N+1 avaliação das derivadas. Para

manter consistência, aplicam-se métodos de passo simples de ordem maior que ou igual

à ordem N, de modo que, se h é o passo, o erro local por passo ( )1+Ο Nh é o mesmo em

ambos os métodos.

A formulação Adams-Basforth-Moulton é geralmente o método de integração

multi-passos mais usado. Os métodos de multi-passos são frequentemente referidos

como métodos preditores-corretores.

6.3.2 - Método de extrapolação

Os métodos de extrapolação usam primeiramente um método de passo fixo com

vários comprimentos de passo para obter uma primeira aproximação e, em seguida,

através de uma extrapolação polinomial ou racional desses valores, calculam o estado

integrado. A essência do método é baseada na extrapolação de Richardson. A idéia é

simples; seja a expansão assintótica:

L+++= 221)( hAhAAhA o

Então as soluções para diversos passos h pode ser obtida por:

1/

)(11

111

−

−+=

=

+

−−+−

+

sms

m

s

m

sm

s

m

s

s

o

s

hh

xxxx

hAx

onde x é o valor extrapolado, o índice s se refere ao passo, e o índice superior m se

refere à extrapolação. Um exemplo simples seria:

)(2/)2/()(

)()(2

111

21

oooo

o

oooo

o

o

hhAAhAhAx

hhAAhAx

Ο++===

Ο++==

- 10 -

Nota-se que A( ho) e A(h1) contêm o erro de 1a ordem. Aplicando a equação de

extrapolação dada, fica:

[ ] [ ])(

)()(2/2

21/

)(

2

21

21

11

11

1

oo

oooooo

o

o

o

o

o

o

o

o

o

o

o

o

hA

hhAAhhAA

xxhh

xxxx

hAx

Ο+=

Ο++−Ο++=

−=−

−+=

=

e o termo de 1a ordem foi eliminado. Aplicando sucessivamente a fórmula de

extrapolação às várias sequências de passos hs, os termos principais de erros vão sendo

sucessivamente eliminados. Esta é a forma da extrapolação polinomial. Outra variante

se refere à chamada extrapolação racional, conforme descrita em Bulirsh e Stoer (1966).

Em geral, quando a avaliação das derivadas f não é muito dispendiosa, o melhor

método é o de extrapolação. Os custo de computação são menores para este método,

porém o método preditor-corretor requer, em geral, consideravelmente menor número

de avaliações das f.

6.4 - Formulações da equações diferenciais do movimento orbital

As equações que descrevem o movimento de um satélite sujeito ao campo

gravitacional terrestre são genericamente dadas por:

∂

∂−+=

rP

rr

3

U

r-µ&& (6.9)

onde P é o vetor que representa as perturbações orbitais dissipativas, e r∂∂ /U é o

gradiente dos potenciais perturbadores de origem conservativa.

Para um órbita puramente kepleriana, as equações do movimento seriam:

3

rr

r-µ=&& (6.10)

- 11 -

Quando a magnitude das perturbações é suficientemente forte para desviar

consideravelmente o movimento kepleriano ou as órbitas são muito excêntricas, as

equações diferenciais podem tornar-se difíceis de lidar devido à não-regularidade,

instabilidade e mesmo descontinuidades. Para contornar tais situações, pode-se utilizar

formulações das equações diferenciais que permitem adequar o tamanho do passo de

integração de forma analítica, por exemplo, no caso de se utilizar um sistema com o

tempo transformado, ou ainda, na utilização de sistemas estabilizados e regularizados ou

um sistema unificado, como se verá em seguida.

6.4.1 - Sistema convencional ou Newtoniano

As equações diferenciais para o sistema convencional são dadas pela Equação

6.9. As perturbações podem ser explicitadas via:

outrasalbedoprsdarrastomarésLuaSoltesseraiszonaisr

- PPPPPPPPr

r3

++++++++= −µ&&

como foi visto anteriormente, onde as 4 primeiras perturbações são de natureza

conservativa, e as 3 seguintes são dissipativas.

6.4.2 - Transformação do tempo

Em órbitas excêntricas, na região próxima ao perigeu a velocidade do satélite é

maior, e menor na região próxima ao apogeu; isto faz com que um comprimento fixo do

passo não divida a órbita em comprimentos iguais de arco quando a variável

independente é o tempo físico. Os diferentes comprimentos de arco tomados ao longo

de uma revolução elevam o nível dos erros obtidos ao término da propagação.

Entretanto, uma regulagem analítica do comprimento de arco pode ser feita, através de

uma transformação de tempo, ou seja, transforma-se a variável independente do sistema.

tempo fisico t, escrevendo-se um novo sistema equivalente ao primeiro com uma

- 12 -

variável independente s, conhecido como o tempo fictício. Esta transformação é

conhecida como transformação de Sundman; isto é:

dsgdt )(r= (6.11)

onde normalmente a função )(rg é da forma αrcg =)(r . A seguir é apresentada uma

tabela para várias transformações de tempo.

Tabela 6.1: Transformações do tempo

α c s corresponde a

0 1 tempo físico

1 µ/a anomalia excêntrica

3/2 ∫ −π

πµ0

2/1)cos1(/ dEEeanomalia intermediária

2 )1(/1 2ea −µ anomalia verdadeira

Assumiremos a transformação α = 1 e c = 1, ou seja, dsrdt = , onde s é

praticamente a anomalia excêntrica a menos de um fator de escala. Define-se então o

operador:

ds

d

rdt

d 1=

Logo, a transformação da formulação convencional nas equações diferenciais para uma

nova variável independente s é dada por:

=

=

ds

d

rdt

d

dt

d

ds

d

rdt

d

rr

rr

1

1

2

2

- 13 -

Aplicando o operador novamente, chega-se a:

−−=

=

2

2

22

2 '111

ds

d

ds

d

r

r

rds

d

rds

d

rdt

d rrrr

com ds

drr =' . Logo,

−−=

∂

∂−+−

2

2

23

'1

ds

d

ds

d

r

r

r

U

r

rr

rPr

µ

Reagrupando os termos, temos:

∂

∂−+

−==

rPr

rr

r Ur

ds

dr

rds

d 22

2

'1

" µ (6.12)

onde )",","(" zyx=r . O símbolo ' representa a derivada com respeito ao tempo fictício

s. Há entretanto a necessidade de se parar a integração em um determinado tempo

físico forçando assim o acréscimo de uma nova equação diferencial ao sistema:

rt =' (6.13)

As equações 6.12 e 6.13 mostram que, após a transformação de tempo, o sistema

passou a ser de sétima ordem, o que aumentará um pouco o tempo de processamento

na integração numérica. Porém, esse acréscimo no tempo de processamento é

compensado pelo expressivo ganho na precisão, devido à característica de regulagem

analítica de passo desse sistema de equações diferenciais.

6.4.3 - Estabilização

O sistema de equações diferenciais orbitais possui comportamento instável

segundo Liapunov. Tal instabilidade pode aumentar os erros globais em uma

propagação orbital, ou até tornar a solução divergente, quando se trabalha com órbitas

de alta excentricidade e com um tamanho de passo de integração considerado pequeno.

Para uma órbita circular de referência, tomando um raio, r2, levemente diferente

do raio de referência, r1, a diferença entre os ângulos polares formados pelos 2 raios ao

- 14 -

longo do tempo não se mantém limitada e tão pequena quanto se queria. De fato, pelas

leis de Kepler sobre a proporcionalidade entre períodos e áreas percorridos, uma vez

que os raios orbitais são diferentes é de esperar que em algum instante os ângulos

polares difiram de até 180°, ficando diametralmente opostos. Portanto não se mantém a

diferença entre ângulos tão pequena quando a desejada e em consequência, as equações

do movimento orbital são instáveis pelo conceito de Liapunov.

Em suma, a estabilidade segundo Liapunov é dada por:

Para algum 0>δ existe um 0>ε tal que:

εθθδ <−⇒<− )()( 2121 ttrr

para todo 0≥t .

Como a velocidade angular orbital depende do raio, o valor absoluto da

diferença dos ângulos tende a crescer no tempo, fato que caracteriza o sistema do

movimento kepleriano circular de Liapunov-instável.

Entretanto, a energia total é uma integral do movimento devido à natureza

conservativa do potencial gravitacional. Assim, imagina-se que um efeito estabilizante

poderia ser adicionado se a solução fosse vinculada para permanecer sobre a superfície

de energia constante. Isto pode ser simplesmente obtido ao introduzir um termo de

controle, em essência a energia, nas equações do movimento, as quais, em

consequência, tenderiam a manter a solução na vizinhança da integral do movimento.

Baumgarte (1972) retrata claramente em seu artigo os efeitos estabilizadores,

principalmente ao se trabalhar com órbitas altamente excêntricas (e > 0.8). Em seus

testes, Baumgarte usa a transformação de tempo dsrdt = em conjunto com o

procedimento de estabilização.

O princípio do método de estabilização caracteriza-se pela inclusão da equação

de energia total na equação da dinâmica do movimento perturbado tendo o tempo

fictício como variável independente.

- 15 -

A equação diferencial para transformação de tempo efetuada é dada pela Equação

6.12. A energia total do movimento é dada por:

02

1 2 =++− HUr

rµ

& (6.14)

sendo H o negativo da energia mecânica total. A equação da energia escrita em função

do tempo fictício s fica na forma:

0)'(

2

12

2

=++− HUrr

r µ(6.15)

Por ser nula, a Equação 6.15 multiplicada pelo vetor posição r e por uma constante w

pode ser introduzida na Equação 6.12 sem alterá-la:

( )

∂

∂−+

++−−−=

rPrrr'r

UrHU

rr

rwr

r

222

)'(2

1'

1"

µµ

Considerando 1=w (condição de estabilidade ), a Equação reduz-se a forma:

∂

∂−+

++−=

rPrr'r

UrHUr

rr

r

222

)'(2

1'

1" (6.16)

O sistema estabilizado possui o tempo fictício como variável independente e há

obviamente a necessidade de parar a integração em um tempo físico qualquer. Com

isso, toma-se necessária a integração da equação base da transformação do tempo:

rt ='

Observa-se também que a energia total do movimento passou a ser uma coordenada

generalizada do sistema, sendo variável no caso se existirem perturbações de natureza

não conservativa. Neste caso, há portanto a necessidade de se integrar a derivada da

equação da energia a fim de avaliá-la em cada passo de integração. A equação da

energia a ser integrada é dada por:

rP ⋅−='H (6.17)

Para se propagar o sistema estabilizado, é necessário que se integre as equações 6.16,

6.17 e 6.13. O sistema passa a ter 8 equações diferenciais de primeira ordem.

- 16 -

6.4.4 - Regularização

A regularização tem como objetivo tomar regular o sistema de equações

diferenciais do movimento orbital, eliminando a singularidade do mesmo na origem do

sistema de referência.

Basicamente tenta-se transformar o sistema em equações do tipo oscilador

harmônico cuja formulação é regular e se mantém válida mesmo para raios pequenos r

(colisão, ponto de singularidade). Desde que o conjunto de equações diferenciais do

oscilador harmônico é Liapunov-estável e apresenta, portanto, propriedades de

estabilidade e regularidade, algum tipo de transformação deve ser necessária para

alcançar formulação semelhante.

Um dos método de regularização é a chamada transformação KS

(Kustaanheimo-Stiefel). Este método propõe uma matriz de transformação em 4

dimensões, a matriz KS:

−−

−−

−−

=

1234

2143

3412

4321

)(

uuuu

uuuu

uuuu

uuuu

L u (6.18)

Com a transformação através da matriz KS, o espaço físico é mapeado para um espaço

de 4 dimensões, através da transformação:

=

4

3

2

1

3

2

1

)(

0 u

u

u

u

Lx

x

x

u (6.19)

Como consequência tem-se as equações:

- 17 -

)(2

)(2

42313

43212

24

23

22

211

uuuux

uuuux

uuuux

+=

−=

+−−=

(6.20)

As coordenadas dos espaços físico e paramétrico em 4 dimensões estão dispostas nos

respectivos vetores posições:

( )

( )T

T

uuuu

xxx

4321

321

=

=

u

r

A transformação possui então a forma compacta:

uur )(L=

É interessante observar que a matriz KS tem as seguintes propriedades:

1. Ortogonalidade: ( )[ ] [ ]11)()()()( rLLLL TT =⋅== uuuuuu

2. )()(' uu LL =

Após a regularização, o sistema de equações diferenciais orbitais assume a forma:

( )

( )rt

LH

)(LU

UH

T

T

=

⋅−=

+

∂

∂−++−=

'

')(2'

2

1

22

1"

2

uPu

Puu

uuu

(6.21)

onde H' é o negativo da energia total do movimento A energia total passou a ser uma

coordenada generalizada do sistema regularizado A variável independente do sistema

regularizado é o tempo fictício. A conhecida necessidade de se parar a integração em

um tempo físico final desejado obriga a integração da equação da derivada do tempo

físico em relação ao fictício, fazendo com que o tempo físico passe a ser uma

coordenada generalizada do sistema. Então, o sistema tem agora 10 equações

diferenciais de 1a ordem.

- 18 -

Considerando

+−

+−−

++−

++

=

322314

312413

342112

332211

)(

PuPuPu

PuPuPu

PuPuPu

PuPuPu

LT Pu

Os valores iniciais dos vetores u e u’ devem ser obtidos a partir dos vetores

posição e velocidade física r e r& no tempo físico inicial.

Das equações 6.20 e considerando que

24

23

22

21 uuuur +++=⋅== uur

pode-se escrever

( ) rxuu +=+ 124

212

uma vez que 1x e r são conhecidos, atribui-se um valor para 1u por exemplo,

encontrando-se consequentemente 4u .

As outras coordenadas do espaço paramétrico são obtidas de

rx

uxuxu

rx

uxuxu

+

−=

+

+=

1

42133

1

43122

Uma vez determinado o vetor posição inicial u, determina-se o vetor velocidade inicial

da seguinte forma:

')()(2')( uuuru LLL TT =

- 19 -

ou ainda,

')(Lr

')(L' TT ruruu

u2

1

2

12

==

Como 'r inicial é conhecido, 'u inicial também será conhecido .

6.4.5 – Sistema Unificado

O sistema unificado é um sistema regularizado mas não regulariza

analiticamente o tamanho do passo como os outros sistemas apresentados.

Os estudos de Altman (1972) sobe a teoria hodográfica (a hodógrafa de um vetor

é a trajetória produzida por sua derivada primeira; no espaço de velocidades, a

hodógrafa é um instrumento que facilita a visualização da aceleração) criaram um

método dinâmico que acoplou de forma genérica os dois tipos de movimento (dinâmica

de atitude e dinâmica de órbita). Este modelo unificado de estados, foi desenvolvido

para usar variáveis de estado e de coordenadas, cujas equações definem vínculos

dinâmicos que possuem propriedades atraentes para eficiência computacional. As

variáveis de estado são momentos (uma forma direta para a atitude e de forma

paramétrica para a órbita) e as variáveis de coordenadas são os parâmetros de Euler, que

representam as transformações de rotação dos eixos de referência. As variáveis de

estado e de coordenadas definem a dinâmica de trajetória e de atitude de forma simples,

comum e bem comportada, possibilitando o processamento separado de 2 conjuntos de

dados, um para a atitude e outro para a órbita, com o uso comum de várias funções

lógicas. No modelo unificado, os paramêtros utilizados são regularizados.

As equações da dinâmica da trajetória orbital são:

+

+−

−

=

3

2

1

2

1

0cos)1(sen

0sen)1(cos

00

e

e

e

f

f

a

a

a

p

p

p

R

R

C

dt

d

γγ

γγ

- 20 -

para os parâmetros hodográficos da órbita, onde a variável R foi decomposta em 2

componentes )( 21 ff RR no plano orbital instantâneo com 1fR ao longo do eixo x e

2fR normal a 1fR . Os parâmetros R e C representam, respectivamente, o raio e a

distância à origem. O significado físico dos parâmetros estão associados com o

momento angular e com a energia cinética. E para os quatérnions, tem-se

−−

−

−=

04

03

02

01

31

31

13

13

04

03

02

01

00

00

00

00

2

1

e

e

e

e

ww

ww

ww

ww

e

e

e

e

dt

d

Juntamente com as seguintes relações auxiliares:

2eV

Cp =

−+

=

2

1

2

1

cossen

sencos0

f

f

e

e

R

R

CV

V

γγ

γγ

−+=

203

204

0403

204

203

21

cos

sen

ee

ee

eeγ

γ

ωχγ +=

2

31

e

e

V

aw =

µ2

3eCV

w =

- 21 -

[ ] ∑=

= i

e

e

e

a

a

a

a

3

2

1

a

onde o vetor a representa o vetor aceleração e ∑ é o somatório das acelerações

perturbadoras. O sub-índices "0" signíficam orbital, os "ei" indicam as direções 1 radial,

2 transversal, e 3 normal. Os 1eV e 2eV são as velocidades radial e tangencial.

Problemas

1 – Considere a transformação dsrdt = . Prove que 2/1)/(2 µπ as = , quando o tempocorresponde ao período orbital no movimento kepleriano; onde a é o semi-eixo maior, eµ é a constante geo-gravitacional.

2 – Considere os seguintes elementos orbitais:

a = 34869261me = 0.8i = 15°Ω = 45°ω = 30°M = 0°

e as constantes 2314 /sm109860064.3 ×=µ , raio da Terra R = 6378139m.

a) Achar a energia do movimento kepleriano,b) Achar as variáveis u e u’ da transformação KS correspondentes à posição e

velocidade,c) Achar a EDO da energia em termos das variáveis u,d) Achar o valor do passo s∆ correspondente a T/20, onde T é o período

orbital. Assuma movimento kepleriano puro.

3 - Escrever as EDOs de primeira ordem do movimento orbital para os sistemasconvencional, transformado, estabilizado e regularizado KS. Assuma transformação

dsrdt = . Em seguida, simplificá-las para o movimento kepleriano puro.

- 22 -

Dica: Para aqueles que não se lembram mais como realizar a transformação deelementos keplerianos para posição e velocidade, as coordenadas de posição evelocidade correspondentes aos elementos keplerianos acima são:

x = 1888980.04103698my = 6652209.67475597mz = 902482.883545056mvx = -9585.79511076297m/svy = 2413.57051166562m/svx = 2273.50409709003m/s

4 - Para os mesmos elementos orbitais anteriores, integrar numericamente a órbita, emprecisão dupla, para um período orbital T, considerando o movimento kepleriano. Usarum integrador Runge-Kutta de 4a ordem com passo fixo T/20, ou o equivalente notempo fictício, para os seguintes casos:

e) Sistema convencional,f) Sistema transformado,g) Sistema estabilizado,h) Sistema regularizado KS.

Montar uma tabela de erros comparativos de posição e velocidade em relação à soluçãoexata (analítica). Considere ),,( aaaa zyx=r e ),,( aaaa zyx &&&=v a solução analítica; e

),,( nnnn zyx=r e ),,( nnnn zyx &&&=v . Calcule os erros da seguinte forma:

Erro em posição = na rr −

Erro em velocidade = na vv −

Referências

Altman, S.P., “A unified state model of orbital trajectory and attitude dynamics,”

Celestial Mechanics, No. 6, 1972, pp. 425-446.

Fehlberg, E. Classical fifth, sixth, seventh, and eigth-order Runge-Kutta formulas with

stepsize control. Washington, DC, NASA 1968 (Nasa Technical Report R-287).

- 23 -

Shanks, E.B. Solutions of differential equations by evaluation of functions.

Mathematics of Computation, 20(93): 1-38, 1966.

Liu, 1974

Liu, 1983

Kolvalevsky, 1967

Kaula, 1966

Pines 1973

Cunningham, 1969

Heiskanen, W. A.; Moritz, H. Physical Geodesy, 1967

Silva, W. C. C.; Ferreira, L. D. D. Satélite artificial – movimento orbital. São José dos

Campos, INPE, XXXX (INPE-3163-RPE/458)

Kondapalli, R. R. Um estudo dos métodos de perturbação na determinação de órbita de

satélites artificiais de baixa altitude. São José dos Campos, INPE, XXXX (INPE-3781-

RPI/151).

Wylie, C. R. Advanced engineering mathematics. 1975.

Kuga, H. K.; Medeiros,V. M.; Carrara, V. Cálculo recursivo da aceleração do

geopotencial. São José dos Campos, INPE, maio de 1983 (INPE-2735-RPE/433).