Procurando a 8ª ordem de Runge-Kutta em C / C ++

Eu gostaria de usar o método de 8ª ordem Runge-Kutta (89) em um aplicativo de mecânica / astrodinâmica celeste, escrito em C ++, usando uma máquina Windows. Portanto, pergunto-me se alguém conhece uma boa biblioteca / implementação documentada e livre para usar? Tudo bem se estiver escrito em C, desde que não haja nenhum problema de compilação esperado.

Até agora encontrei esta biblioteca (mymathlib) . O código parece bom, mas não encontrei nenhuma informação sobre licenciamento.

Você pode me ajudar, revelando algumas das alternativas que você talvez conheça e que atendam ao meu problema?

EDIT:
Vejo que realmente não existem muitos códigos-fonte C / C ++ disponíveis como eu esperava. Portanto, uma versão do Matlab / Octave também ficaria bem (ainda precisa ser livre para usar).

A Biblioteca Científica GNU (GSL) (C) e o Boost Odeint (C ++) apresentam métodos Runge-Kutta de 8ª ordem.

Ambos são de código aberto e, no linux e mac, devem estar diretamente disponíveis no gerenciador de pacotes. No Windows, provavelmente será mais fácil usar o Boost do que o GSL.

A GSL é publicada sob a licença GPL e o Boost Odeint sob a licença Boost.

Edit: Ok, o Boost Odeint NÃO possui o método Runge-Kutta 89, apenas o 78, mas fornece uma receita para fazer steppers arbitrários do Runge-Kutta.

Os métodos de 8ª ordem são bastante altos e provavelmente um exagero para o seu problema.

Prince-Dormand refere - se a um tipo específico de Runge-Kutta e não está diretamente relacionado à ordem, mas o mais comum é 45. Matlabs ode45, que é o algoritmo ODE recomendado, é uma implementação do Prince-Dormand 45. Esse é o mesmo algoritmo implementado no Boost Odeint Runge_Kutta_Dopri5 .

Se você estiver fazendo mecânica celeste em longas escalas de tempo, o uso de um integrador Runge-Kutta clássico não economizará energia. Nesse caso, o uso de um integrador simplético provavelmente seria melhor. O Boost.odeint também implementa um esquema Runge-Kutta simplético de quarta ordem que funcionaria melhor por longos intervalos de tempo. A GSL não implementa nenhum método simplético, até onde eu sei.

resumindo alguns pontos:

1. Se for uma integração de longo prazo de um modelo não dissapativo, é o que você procura por um integrador simplético.
2. Caso contrário, como é uma equação de movimento, os métodos Runge-Kutta Nystrom serão mais eficientes do que uma transformação em um sistema de primeira ordem. Existem métodos RKN de alta ordem devido ao DP. Existem algumas implementações, como aqui em Julia, elas estão documentadas e aqui está uma MATLAB .
3. Os métodos Runge-Kutta de alta ordem são necessários apenas se você deseja uma solução de alta precisão. Se houver tolerâncias mais baixas, um RK de 5ª ordem provavelmente será mais rápido (para o mesmo erro). A melhor coisa a fazer se você precisar resolver isso com frequência é testar vários métodos diferentes. Neste conjunto de parâmetros de referência para problemas de três corpos , vemos que (para o mesmo erro) os métodos RK de alta ordem são apenas uma melhoria marginal na velocidade, embora como erro -> 0, você pode ver que a melhoria já vai para> 5x contra Dormand -Príncipe 45 ( DP5) quando você olha com 4 dígitos de precisão (as tolerâncias são muito mais baixas para isso. As tolerâncias são apenas uma estimativa em qualquer problema). À medida que você reduz as tolerâncias, diminui ainda mais a melhoria de um método RK de ordem superior, mas pode ser necessário começar a usar números de precisão mais altos.
4. O algoritmo de ordem 7/8 de Dormand-Prince possui um quadro de oitava ordem diferente do método DP853 de Hairer dop853e DifferentialEquations.jl DP8(que são iguais). O último método 853 não pode ser implementado na versão padrão do método Runge-Kutta, pois o estimador de erros não é padrão. Mas esse método é muito mais eficiente e eu não recomendaria nem o uso dos métodos Fehlberg 7/8 ou DP 7/8 mais antigos.
5. Para métodos RK de alta ordem, os métodos Verner "Efficient" são o padrão-ouro. Isso aparece nos benchmarks que eu vinculei. Você pode codificá-las no Boost você mesmo ou usar um dos 2 pacotes que as implementam, se você quiser os mais fáceis (Mathematica ou DifferentialEquations.jl).

Gostaria de acrescentar que, embora o que Geoff Oxberry sugira para integração a longo prazo (usando integradores simpléticos) seja verdadeiro, em alguns casos não funcionará. Mais especificamente, se você possui forças dissipativas, seu sistema não preserva mais a energia e, portanto, você não pode recorrer a integradores simpléticos nesse caso. A pessoa que fez a pergunta estava falando sobre órbitas baixas da Terra, e essas órbitas exibem uma grande quantidade de arrasto atmosférico, que é uma força dissipativa que impede o uso de tais integradores simpléticos.

Nesse caso específico (e nos casos em que você não pode usar / não tem acesso / não deseja usar integradores simpléticos), eu recomendaria o uso do integrador Bulirsch-Stoer se precisar de precisão e eficiência em prazos longos. Funciona bem por experiência e também é recomendado pelas Receitas Numéricas (Press et al., 2007).