Por que os métodos Runge – Kutta de ordem superior não são usados ​​com mais frequência?

17

Eu estava curioso para saber por que os métodos Runge-Kutta de alta ordem (ou seja, maiores que 4) quase nunca são discutidos / empregados (pelo menos que eu saiba). Entendo que isso requer maior tempo computacional por etapa (por exemplo, RK14 com etapa incorporada de 12ª ordem ), mas existem outras desvantagens do uso de métodos Runge – Kutta de ordem superior (por exemplo, problemas de estabilidade)? Quando aplicados a equações com soluções altamente oscilantes em escalas de tempo extremas, esses métodos de ordem superior não seriam normalmente os preferidos?

Mathews24
fonte
2
Eu acho que essa é uma pergunta muito subjetiva. A maior desvantagem, como você já observou, é o custo da computação. Geralmente tentamos equilibrar entre precisão e tempo computacional. Nos PDE, quando as pessoas falam sobre ordem superior, geralmente pensam em terceira ou quarta ordem. E o intervalo de tempo também é mantido na mesma ordem.
21716 Vikram
3
No PDE, um esquema de precisão de alta ordem para dependência temporal não faz sentido se a precisão espacial for pior. De fato, a precisão da dependência espacial é principalmente de 2ª ou 3ª ordem, especialmente quando se trabalha em malhas não estruturadas. As pessoas precisam controlar o truncamento global de erros com o menor custo, portanto, considera o Runge-Kutta com uma ordem de precisão alta o suficiente em casos específicos.
tqviet
@tqviet Se usando aproximações de diferença central ou para trás até a ordem 8 para as derivadas espaciais, RK8 seria adequado, não? Em geral, há problemas de precisão ou estabilidade com o uso de aproximações de diferenças finitas de alta ordem das derivadas espaciais?
Mathews24
11
@ Mathews24: Eu não mencionei a estabilidade, que depende fortemente da equação. Quando um esquema altamente preciso é aplicada a dependência espacial, adotamos RK a dependência temporal com pelo menos a mesma ordem de precisão, mas a condição de estabilidade pode exigir um valor menor de . Δt
tqviet

Respostas:

17

Existem milhares de documentos e centenas de códigos por aí, usando os métodos Runge-Kutta de quinta ordem ou superior. Observe que o integrador explícito mais usado no MATLAB é o ODE45, que avança a solução usando um método Runge-Kutta de 5ª ordem.

Exemplos de métodos Runge-Kutta de alta ordem amplamente utilizados

O artigo de Dormand & Prince, dando um método de 5ª ordem, tem mais de 1700 citações, de acordo com o Google Scholar . A maioria deles são documentos usando seu método para resolver algum problema. O artigo do método Cash-Karp tem mais de 400 citações . Talvez o método de ordem mais amplamente utilizado, superior a 5, seja o método de 8ª ordem de Prince-Dormand, que possui mais de 400 citações no Google Scholar . Eu poderia dar muitos outros exemplos; e lembre-se de que muitas (se não a maioria) das pessoas que usam esses métodos nunca citam os documentos.

Observe também que os métodos de extrapolação de alta ordem e correção diferida são métodos de Runge-Kutta .

Métodos de alta ordem e erro de arredondamento

Se sua precisão é limitada por erros de arredondamento, você deve usar um método de ordem superior . Isso ocorre porque os métodos de ordem superior exigem menos etapas (e menos avaliações de função, mesmo que haja mais avaliações por etapa), portanto eles cometem menos erros de arredondamento. Você pode verificar isso facilmente com experimentos simples; é um bom problema de lição de casa para um primeiro curso de análise numérica.

Os métodos de décima ordem são extremamente úteis na aritmética de precisão dupla. Pelo contrário, se tudo o que tínhamos fosse o método de Euler, o erro de arredondamento seria um problema importante e precisaríamos de números de ponto flutuante de alta precisão para muitos problemas em que os solucionadores de alta ordem funcionam bem.

Métodos de alta ordem podem ser igualmente estáveis

UMAB (uma propriedade de estabilidade útil para alguns problemas não lineares). Para aprender sobre esses métodos, consulte, por exemplo, o texto da Hairer & Wanner.

Métodos de alta ordem em mecânica celeste

Você pergunta

Quando aplicados a equações com soluções altamente oscilantes em escalas de tempo extremas, esses métodos de ordem superior não seriam normalmente os preferidos?

Você está exatamente certo! Um excelente exemplo disso é a mecânica celeste. Eu não sou um especialista nessa área. Mas este artigo , por exemplo, compara métodos para a mecânica celeste e nem considera a ordem inferior a 5. Conclui que os métodos da ordem 11 ou 12 são frequentemente os mais eficientes (com o método da ordem 8 de Prince-Dormand também muitas vezes eficiente).

David Ketcheson
fonte
Ketchson: você poderia fornecer alguma evidência ou explicação sobre esta afirmação: "métodos de extrapolação de alta ordem e correção diferida são métodos de Runge-Kutta"? Especialmente os "métodos de correção diferidos". Obrigado.
tqviet
@ David Ketcheson Você pode discutir como sua resposta mudaria se usasse técnicas de computação validadas (verificadas), como intervalo arredondado para fora ou aritmética radial? Que tal se fosse utilizado um intervalo arredondado para fora mais alto que a precisão dupla ou uma aritmética radial? O que acontecerá com a quebra e a dependência à medida que a ordem Runge-Kutta for aumentada e, apenas por diversão, digamos que a ODE é muito rígida.?
Mark L. Stone
@ MarkL.Stone Esse é um conjunto de perguntas completamente diferente. Se você quiser perguntar a eles, envie-os como perguntas separadas. No entanto, não sou especialista nessas coisas e não poderei responder.
David Ketcheson
11
@tqviet Veja este artigo para obter uma explicação.
David Ketcheson
12

Desde que você esteja usando a aritmética padrão de ponto flutuante de precisão dupla, métodos de ordem muito alta não são necessários para obter uma solução com alta precisão em um número razoável de etapas. Na prática, acho que a precisão da solução é normalmente limitada a um erro relativo de 1.0e-16 pela representação de ponto flutuante de precisão dupla, em vez do número / comprimento das etapas executadas com o RKF45.

Se você mudar para um esquema aritmético de ponto flutuante de precisão mais alta que dupla, o uso de um método de 10ª ordem pode valer a pena.

Brian Borchers
fonte
5
Eu acho que essa resposta é enganosa. Os métodos de alta ordem levam a muito menos erro de arredondamento, enquanto os métodos de baixa ordem sofrem com o erro de arredondamento ser dominante quando a precisão necessária é grande ou o intervalo de tempo é longo; veja minha resposta abaixo.
David Ketcheson 18/11/2016
2
O ponto é que, no ponto flutuante de precisão dupla, você não pode sequer representar uma solução com precisão relativa superior a 1,0e-16. Em muitas situações práticas, o bom e velho RKF45 o levará a esse nível de precisão durante o período em que você estiver interessado, sem a necessidade de pequenos passos. Pode não ser uma boa escolha para sistemas rígidos ou situações em que é necessário um integrador simplético, mas um método Runge Kutta de ordem superior também não é uma ótima solução para essas situações. Concordo que, por períodos muito longos, os métodos Runge Kutta de ordem superior podem fazer algum sentido.
precisa
10

Apenas para acrescentar à excelente resposta de Brian Borcher, muitos aplicativos da vida real admitem ODEs ou DAEs altamente rígidos. Intuitivamente, esses problemas experimentam mudanças bruscas e não suaves ao longo do tempo; portanto, são melhor modelados usando polinômios de ordem inferior, distribuídos finamente por tamanhos curtos, em oposição a polinômios de alta ordem estendidos por tamanhos longos. Além disso, a estabilidade frequentemente requer o uso de métodos implícitos , para os quais a penalidade computacional de métodos de ordem superior é muito mais acentuada.

Mais rigorosamente, os métodos de ordem superior são menos estáveis ​​que os métodos de ordem inferior para problemas rígidos. Temos, por exemplo, as barreiras de Dahlquist para métodos lineares de várias etapas.

r2

Declarações semelhantes (mas muito mais complicadas) podem ser feitas para a estabilidade de L nas fórmulas RK. Em todos os casos, o aumento na ordem nem sempre leva a soluções mais precisas. A seguir, um trecho do artigo seminal de 1974 de Prothero e Robinson:

Ao usar métodos de uma etapa estáveis ​​A para resolver grandes sistemas de equações diferenciais não lineares rígidas, descobrimos que
(a) alguns métodos estáveis ​​A fornecem soluções altamente instáveis ​​e
(b) a precisão das soluções obtidas quando as equações são rígido normalmente parece não ter relação com a ordem do método usado.

Para tratamentos ainda mais rigorosos deste tópico, consulte o texto clássico de Hairer & Wanner, "Resolvendo equações diferenciais ordinárias II: Problemas rígidos e diferenciais - algébricos", 1991.

Na prática, equações rígidas são quase sempre resolvidas usando a regra trapezoidal ou a fórmula TR-BDF2 (funções ode23t e ode23tb no MATLAB). Ambos são métodos implícitos de segunda ordem. Obviamente, onde a estabilidade não é um problema (ou seja, em equações não rígidas), somos livres para escolher entre várias opções; RK45 é a escolha mais comum.

Richard Zhang
fonte
Muito interessante. Existe alguma explicação (intuitiva) para o motivo pelo qual o pedido deve ser menor ou igual a 2 para que ele seja um método de etapas múltiplas estável em A? E apenas para esclarecer, quando você diz que declarações semelhantes podem ser feitas para fórmulas RK, é da ordem 2 mais uma vez?
Mathews24
Mas para os métodos Runge-Kutta, existem métodos A-estáveis ​​de ordem arbitrária.
David Ketcheson
@DavidKetcheson Sim, mas eles não são fortemente estáveis ​​em A (ou seja, estáveis ​​em L). Eles têm muitos problemas quando usados ​​para resolver DAEs, por exemplo, simular circuitos simples de transistor. De fato, o TR é famoso por causar toque artificial no SPICE, que foi o que motivou o desenvolvimento do TR-BDF2.
Richard Zhang
@DavidKetcheson Para referência, consulte doi.org/10.1090/S0025-5718-1974-0331793-2 . A noção de estabilidade A não é suficientemente forte para DAEs, e os métodos estáveis ​​A de alta ordem geralmente produzem resultados estranhos quando usados ​​para resolver DAEs.
Richard Zhang
Claro, mas a questão não é sobre DAEs ou sobre métodos de várias etapas.
David Ketcheson
9

A configuração de benchmark

No software Julia DifferentialEquations.jl , implementamos muitos métodos de ordem superior, incluindo os métodos Feagin. Você pode vê-lo em nossa lista de métodos e, em seguida, existem muitos outros que podem ser usados como tabelas fornecidas . Como todos esses métodos são reunidos, é possível comparar facilmente entre eles. Você pode ver os benchmarks que eu tenho online aqui e ver que é muito simples comparar muitos algoritmos diferentes. Portanto, se você quiser levar alguns minutos para executar os benchmarks, faça isso. Aqui está um resumo do que sai.

Primeiramente, é importante observar que, se você olhar para cada um dos benchmarks, verá que nossos métodos DP5e ordem (Ordem 5 de Dormand-Prince) DP8são mais rápidos que os códigos Hairer Fortran ( dopri5e dop853) e, portanto, essas implementações são muito bem otimizadas . Isso mostra que, conforme observado em outro segmento, o uso excessivo dos métodos Dormand-Prince é porque os métodos já foram escritos, não porque ainda são os melhores. Portanto, a comparação real entre as implementações mais otimizadas é entre os métodos Tsitorous, Verner e Feagin do DifferentialEquations.jl.

Os resultados

Em geral, os métodos de uma ordem superior a 7 têm um custo computacional adicional que geralmente não é compensado pela ordem, dadas as tolerâncias escolhidas. Uma razão para isso é que as opções de coeficiente para métodos de ordem inferior são mais otimizadas (elas têm pequenos "coeficientes de erro de truncamento de princípio", que importam mais quando você não é assimtopicamente pequeno). Você pode ver que em muitos problemas, como aqui, os métodos Verner Efficient 6 e 7 se saem extremamente bem, mas métodos como o Verner Efficient 8 podem ter uma inclinação mais baixa. Isso ocorre porque os "ganhos" de ordem superior são compostos com tolerâncias mais baixas, portanto sempre há uma tolerância em que os métodos de ordem superior serão mais eficientes.

No entanto, a questão é então, quão baixo? Em uma implementação bem otimizada, isso fica muito baixo por dois motivos. A primeira razão é porque os métodos de ordem inferior implementam algo chamado FSAL (primeiro igual ao anterior). Essa propriedade significa que os métodos de ordem inferior reutilizam uma avaliação de função da etapa anterior na próxima etapa e, portanto, têm efetivamente uma avaliação de função a menos. Se isso for usado corretamente, algo como um método de 5ª ordem (Tsitorous ou Dormand-Prince) está na verdade fazendo 5 avaliações de função em vez das 6 que os tablóides sugerem. Isso também se aplica ao método Verner 6.

A outra razão é devido a interpolações. Uma razão para usar um método de ordem muito alta é executar menos etapas e simplesmente interpolar valores intermediários. No entanto, para obter os valores intermediários, a função de interpolação pode precisar de mais avaliações de função do que as usadas para dar o passo. Se você olhar para os métodos Verner, são necessárias oito avaliações de funções extras para o método Pedido 8 para obter um interpolante Pedido 8. Muitas vezes, os métodos de baixa ordem fornecem um interpolante "gratuito", por exemplo, a maioria dos métodos de 5ª ordem tem uma interpolação de 4ª ordem gratuita (sem avaliações de funções extras). Portanto, isso significa que, se você precisar de valores intermediários (que serão necessários para um bom gráfico, se você estiver usando um método de alta ordem), haverá alguns custos ocultos extras. Considere o fato de que esses valores interpolados são realmente importantes para o tratamento de eventos e a solução de equações diferenciais de atraso e você verá por que o custo extra da interpolação é fator importante.

Então, e os métodos Feagin?

Portanto, você verá que os métodos Feagin estão faltando suspeitosamente nos benchmarks. Eles são bons, os testes de convergência funcionam em números de precisão arbitrários etc., mas para realmente fazê-los funcionar, você precisa pedir tolerâncias absurdamente baixas. Por exemplo, descobri em benchmarks não publicados que o Feagin14desempenho supera Vern9(o Método Eficiente Verner de 9ª ordem) em tolerâncias semelhantes 1e-30. Para aplicações com dinâmica caótica (como nos problemas de Pleides ou astrofísica de três corpos), você pode querer essa quantidade de precisão devido à dependência sensível (erros nos sistemas caóticos são compostos rapidamente). No entanto, a maioria das pessoas provavelmente está computando com números de ponto flutuante de precisão dupla, e eu não encontrei uma referência em que eles tenham desempenho superior nesse domínio de tolerância.

Além disso, não há interpolante para ir junto com os métodos Feagin. Então, o que eu faço é simplesmente colocar uma interpolação Hermite de terceira ordem sobre eles, para que assim exista (e funciona surpreendentemente bem). No entanto, se não houver uma função de interpolação padrão, você poderá executar o método Hermite recursivo (use essa interpolação para obter o ponto médio e, em seguida, faça uma interpolação de 5ª ordem, etc.) para obter uma interpolação de alta ordem, mas isso é muito caro e o resultado interpolação não tem necessariamente um termo de erro de truncamento de princípio baixo (portanto, só é bom quando dté realmente pequeno, que é exatamente o oposto do caso que queremos!). Portanto, se você precisar de uma interpolação realmente boa para corresponder à sua precisão, precisará pelo menos voltar a algo assim Vern9.

Nota sobre extrapolação

Observe que os métodos de extrapolação são simplesmente algoritmos para gerar métodos Runge-Kutta de ordem arbitrária. No entanto, por sua ordem, eles tomam mais medidas do que o necessário e têm altos coeficientes de erro de truncamento de princípio, e, portanto, não são tão eficientes quanto um método RK bem otimizado em uma determinada ordem. Mas, dada a análise anterior, isso significa que existe um domínio de tolerância extremamente baixa em que esses métodos se saem melhor do que os métodos RK "conhecidos". Mas em todos os benchmarks que corri, parece que não cheguei tão baixo.

Nota sobre estabilidade

A escolha realmente não tem nada a ver com problemas de estabilidade. De fato, se você passar pelos quadros DifferentialEquations.jl (você pode apenas plot(tab)para as regiões de estabilidade), verá que a maioria dos métodos possui regiões de estabilidade suspeitamente semelhantes. Esta é realmente uma escolha. Geralmente, ao derivar os métodos, o autor geralmente faz o seguinte:

  1. Encontre os menores coeficientes de erro de truncamento do princípio (ou seja, os coeficientes para os próximos termos do pedido)
  2. Sujeito às restrições do pedido
  3. E torne a região de estabilidade próxima à do método Dormand-Prince Order 5.

Por que a última condição? Bem, como esse método tende a ser sempre estável com a maneira como são feitas as escolhas adaptativas de etapas controladas por PI, é uma boa opção para regiões de estabilidade "suficientemente boas". Portanto, não é por acaso que todas as regiões de estabilidade tendem a ser semelhantes.

Conclusão

Existem trocas em todas as opções de método. Os métodos RK de ordem mais alta simplesmente não são tão eficientes em tolerâncias mais baixas, porque é mais difícil otimizar a escolha de coeficientes e porque o número de avaliações de funções se compõe (e cresce ainda mais rápido quando há interpolação). No entanto, se a tolerância ficar baixa o suficiente, elas vencem, mas as tolerâncias necessárias podem estar muito abaixo das aplicações "padrão" (isto é, realmente aplicáveis ​​apenas a sistemas caóticos).

Chris Rackauckas
fonte