Quero apenas saber se o método numérico de Dormand-Prince ou o método numérico de Cash-Karp é mais preciso.

Como acabei de otimizar muitos deles em um software, DifferentialEquations.jl , decidi fazer uma comparação dos principais métodos do Order 4/5. O método Fehlberg foi deixado de fora porque é geralmente conhecido por ser menos eficiente que o método DP5.

Backstories

Dormand-Prince 4/5

O método Dormand-Prince foi desenvolvido para ser preciso como um par 4/5 com uso de extrapolação local (ou seja, passo com o par da ordem 5. Isso ocorre porque ele foi projetado para estar próximo do coeficiente de erro de truncamento do princípio ideal (ou seja, mínimo) (sob a restrição de também ter o número mínimo de etapas para atingir a ordem 5. Ele possui uma interpolação de ordem 4 que é gratuita, mas precisa de etapas extras para uma interpolação de ordem 5.

Cash-Karp 4/5

O método Cash-Karp foi desenvolvido para satisfazer diferentes restrições, ou seja, para lidar melhor com problemas não suaves. Eles escolheram o , a percentagem do instante temporal no th passo (ou seja, é o tempo que o th passo é calculada a) para ser tão uniforme quanto possível, mas ainda alcançar a ordem 5. Em seguida, ele também foi derivado de ter incorporado métodos de 1ª, 2ª, 3ª e 4ª ordem com essa uniformidade do$c_i$$i$$t+c_i \Delta t$$i$$c_i$. Eles são espaçados de maneira que você possa descobrir onde começa uma parte rígida pela qual a diferença é grande. Além disso, observe que quanto mais rígida a equação, pior é o método de ordem superior (porque ele precisa limitar as derivadas mais altas). Então eles desenvolvem uma estratégia que usa os 5 métodos incorporados para "sair mais cedo": ou seja, se você detectar rigidez, pare no estágio$i<6$para diminuir o número de chamadas de função e economizar tempo. Portanto, no final, esse "par" foi desenvolvido com muitas outras restrições em mente e, portanto, não há razão para esperar que seja "mais preciso", pelo menos como um par de 4/5. Se você adicionar todo esse outro equipamento, em problemas (semi-) rígidos, ele será mais preciso (mas, nesse caso, convém usar um método diferente, como o método W-Rosenbrock). Essa é uma das razões pelas quais esse par não se tornou padrão em relação ao DP5, mas ainda pode ser útil (talvez seja bom para um método híbrido que alterna para um solucionador rígido quando a rigidez é encontrada?).

Bogacki & Shampine 4/5

Para completar a resposta, vamos discutir o par Bogacki & Shampine mencionado no comentário. O método BS5 elimina a restrição de "usar o mínimo de chamadas de função" (usa 8 em vez de 6) para fazer duas coisas:

• Obtenha coeficientes de erro de truncamento de princípio realmente baixos.

• Produza uma interpolação da ordem 5 com coeficientes de erro mais baixos.

Esses coeficientes são tão baixos que, para muitos problemas com tolerâncias que os usuários provavelmente usam, ele mede como se fosse de sexta ordem. O artigo deles mostra que, para chamadas de funções baratas, isso pode ser mais eficiente que o DP5, aproximadamente o mesmo valor no DP5 era superior ao RKF5 (o método Fehlberg).

Você pode juntar dois e dois e ver: espere um segundo, Shampine é a mesma pessoa que desenvolveu a suíte MATLAB ODE; isso foi depois que o artigo do par BS5 foi publicado, por que o MATLAB não ode45usa o par BS5? Um dos motivos é que isso foi feito antes do relés do par BS5. A outra razão é porque a ode45função foi desenvolvida para minimizar o tempo. Enquanto o par BS5 é mais eficiente (ou seja, obtém menor precisão), o objetivo deode45é ter um erro bom o suficiente para fazer uma plotagem boa o suficiente. Isso significa que, para lidar com os grandes passos, ele também produz duas soluções interpoladas extras entre cada passo. Para o método DP5, há uma interpolação de ordem "livre" 4 e, portanto, isso é muito mais rápido do que usar o BS5. Como também é "preciso o suficiente" em tolerâncias moderadas, esse método é definido como padrão, pois fornece uma experiência de usuário padrão melhor do que o BS5 ao fazer computação interativa (portanto, essa opção foi específica ao contexto).

Tsitorous 4/5

Aqui está um a menos que as pessoas conhecem. É derivado neste artigo . É derivado usando menos suposições do que o método DP5 e tenta obter um par com coeficientes de erro de truncamento de princípio mais baixos. Em seus testes, afirma que consegue isso. Ele também possui uma interpolação de ordem gratuita 4 como o método DP5.

Testes numéricos

Eu escrevi o pacote numérico DifferentialEquations.jl para ser um conjunto bastante abrangente de solucionadores para Julia. Ao longo do caminho de guerra, implementei mais de 100 métodos Runge-Kutta e otimizei bastante a mão. Três dos integradores otimizados para mão são os métodos DP5, BS5 e Tsit5 (eu não fiz o CK5 porque, como observado na história anterior, o principal caso é de problemas que são meio rígidos. Acho que a melhor maneira de lidar com eles é usar o DP5 / BS5 e alternar para solucionadores rígidos, conforme necessário, de uma maneira como LSODE, mas essa é uma história para um tempo diferente) (uma maneira de ver que eles estão perto do ideal é que esses métodos são mais rápidos que as dopri5implementações da Hairer , então eles são pelo menos implementações decentes). Testes entre muitos métodos de Runge-Kutta em equações não rígidaspode ser encontrado na pasta benchmarks . Estou adicionando mais conforme vou avançando, mas você pode ver nos diagramas de precisão do ODE linear e do problema dos três corpos, medimos os métodos DP5 e Tsit5 para ter uma eficiência quase idêntica, superando o método BS5 no ODE linear, enquanto DP5 e BS5 são quase idênticos no problema dos três corpos com o Tsit5 atrás. A partir dessas informações, pelo menos por enquanto, decidi pelo método DP5 como padrão, correspondendo às recomendações anteriores. Isso pode mudar com testes futuros (ou você pode adicionar pontos de referência! Sinta-se à vontade para contribuir ou marque a repo como estrela para dar mais suporte a esse esforço).

Conclusão

Em conclusão, os pares da Ordem 5 são assim:

• O par Dormand-Prince 4/5 é um bom par útil, pois é bem otimizado em termos de coeficiente de erro de truncamento de princípio e possui uma interpolação barata de ordem 4, o que torna mais rápido a produção de plotagens decentes.

• O par Cash-Karp possui mais restrições para lidar melhor com equações rígidas. No entanto, para obter o benefício completo, você precisará usar o algoritmo completo com os 5 métodos incorporados.

• O método Bogacki e Shampine Order 5 pode ser o mais eficiente em termos de produção de erro por chamada de função (ele possui um estimador de erro duplo; portanto, em problemas mais difíceis, provavelmente se sai melhor), mas isso permite que ele demore mais tempo. No entanto, se você deseja apenas produzir uma plotagem suave, é necessário reagir com este método: use uma tolerância mais baixa (para que demore mais que DP5, mas com menos erro) ou use etapas mais interpoladas. No final, isso significava que talvez não fosse melhor para aplicativos interativos, embora fosse melhor para alguns aplicativos de computação científica.

• O Tsitorous 4/5. Foi desenvolvido recentemente (2011) para superar o DP5 em uma comparação frente a frente. Meus testes não me dão uma razão para acreditar que é muito melhor que o DP5 que agora deve ser considerado como o novo método padrão, mas testes futuros podem começar a ficar a seu favor.

Editar

Eu melhorei a implementação do Tsit5. Agora ele se sai melhor que o DP5 na maioria dos testes, nas implementações DifferentialEquations.jl e Hairer dopri (embora se possa surpreender que as implementações DifferentialEquations.jl sejam realmente mais rápidas, o que naturalmente ajuda a implementação do Tsit5). Agora eu o recomendo como o método padrão da ordem 4/5.

Se você estiver interessado em comparar dois integradores, resolva

com valores iniciais e para ambos. Em seguida, plote os erros$(q,p) = (1,0)$$dt = 0.1$

para um intervalo razoável de (0 a 100, mil etapas no total) e escolha o integrador com o menor erro. Este teste do oscilador harmônico mostrará os erros de fase e amplitude com muito pouco esforço.$t$

• PJ Channell e C Scovel. Integração simplética de sistemas hamiltonianos . Journal of Computational Physics 73, 468 (1973).