J'étais simplement curieux de savoir pourquoi les méthodes Runge – Kutta d'ordre élevé (c'est-à-dire supérieures à 4) ne sont presque jamais discutées / utilisées (du moins à ma connaissance). Je comprends que cela nécessite un temps de calcul plus long par étape (par exemple RK14 avec une étape intégrée de 12e ordre ), mais y a-t-il d'autres inconvénients à utiliser des méthodes Runge – Kutta d'ordre supérieur (par exemple, des problèmes de stabilité)? Lorsqu'elles sont appliquées à des équations avec des solutions très oscillantes sur des échelles de temps extrêmes, de telles méthodes d'ordre supérieur ne seraient-elles pas généralement préférées?

Il existe des milliers d'articles et des centaines de codes utilisant des méthodes Runge-Kutta de cinquième ordre ou plus. Notez que l'intégrateur explicite le plus couramment utilisé dans MATLAB est ODE45, qui fait avancer la solution en utilisant une méthode Runge-Kutta du 5e ordre.

Exemples de méthodes Runge-Kutta d'ordre élevé largement utilisées

Le papier de Dormand & Prince donnant une méthode du 5ème ordre a plus de 1700 citations selon Google Scholar . La plupart d'entre eux sont des articles utilisant leur méthode pour résoudre un problème. Le document sur la méthode Cash-Karp compte plus de 400 citations . La méthode d'ordre supérieur à 5 la plus utilisée est peut-être la méthode du 8ème ordre de Prince-Dormand qui compte plus de 400 citations sur Google Scholar . Je pourrais donner beaucoup d'autres exemples; et gardez à l'esprit que beaucoup (sinon la plupart) des personnes utilisant ces méthodes ne citent jamais les papiers.

Notez également que l'extrapolation d'ordre élevé et les méthodes de correction différée sont des méthodes de Runge-Kutta .

Méthodes d'ordre élevé et erreur d'arrondi

Si votre précision est limitée par des erreurs d'arrondi, vous devez utiliser une méthode d'ordre supérieur . En effet, les méthodes de niveau supérieur nécessitent moins d'étapes (et moins d'évaluations de fonctions, même s'il y a plus d'évaluations par étape), de sorte qu'elles commettent moins d'erreurs d'arrondi. Vous pouvez facilement le vérifier vous-même avec des expériences simples; c'est un bon problème de devoirs pour un premier cours d'analyse numérique.

Les méthodes du dixième ordre sont extrêmement utiles en arithmétique double précision. Au contraire, si tout ce que nous avions était la méthode d'Euler, l'erreur d'arrondi serait un problème majeur et nous aurions besoin de nombres à virgule flottante de très haute précision pour de nombreux problèmes où les solveurs d'ordre supérieur fonctionnent très bien.

Les méthodes de haut niveau peuvent être tout aussi stables

@RichardZhang a référencé la deuxième barrière Dahlquist, mais cela ne s'applique qu'aux méthodes en plusieurs étapes. La question publiée ici concerne les méthodes Runge-Kutta, et il existe des méthodes Runge-Kutta de chaque ordre qui sont non seulement stables, mais aussi stables$A$$B$ (une propriété de stabilité utile pour certains problèmes non linéaires). Pour en savoir plus sur ces méthodes, voir par exemple le texte de Hairer & Wanner.

Méthodes d'ordre élevé en mécanique céleste

Tu demandes

Lorsqu'elles sont appliquées à des équations avec des solutions très oscillantes sur des échelles de temps extrêmes, de telles méthodes d'ordre supérieur ne seraient-elles pas généralement préférées?

Tu as parfaitement raison! Un bon exemple de ceci est la mécanique céleste. Je ne suis pas un expert dans ce domaine. Mais cet article , par exemple, compare des méthodes de mécanique céleste et ne considère même pas un ordre inférieur à 5. Il conclut que les méthodes d'ordre 11 ou 12 sont souvent les plus efficaces (avec la méthode Prince-Dormand d'ordre 8 également souvent très efficace).

Tant que vous utilisez l'arithmétique à virgule flottante double précision standard, aucune méthode d'ordre très élevé n'est nécessaire pour obtenir une solution avec une grande précision en un nombre raisonnable d'étapes. Dans la pratique, je trouve que la précision de la solution est normalement limitée à une erreur relative de 1,0e-16 par la représentation en virgule flottante à double précision plutôt que par le nombre / la longueur des étapes qui sont effectuées avec RKF45.

Si vous passez à un schéma arithmétique à virgule flottante d'une précision supérieure à la double précision, l'utilisation d'une méthode du 10e ordre pourrait valoir la peine.

Pour ajouter à l'excellente réponse de Brian Borcher, de nombreuses applications réelles admettent des ODE ou DAE très rigides. Intuitivement, ces problèmes subissent des changements brusques et brusques au fil du temps, ils sont donc mieux modélisés à l'aide de polynômes de faible ordre répartis finement sur de courtes tailles de pas, par opposition aux polynômes de haut niveau étirés sur de longues tailles de pas. De plus, la stabilité nécessite souvent l'utilisation de méthodes implicites , pour lesquelles la pénalité de calcul des méthodes d'ordre supérieur est beaucoup plus prononcée.

Plus rigoureusement, les méthodes d'ordre supérieur sont moins stables que les méthodes d'ordre inférieur pour les problèmes rigides. Nous avons par exemple les barrières Dahlquist pour les méthodes linéaires à plusieurs étapes.

$r\le2$

Des déclarations similaires (mais beaucoup plus compliquées) peuvent être faites pour la stabilité L dans les formules RK. Dans tous les cas, l'augmentation de l'ordre ne conduit pas toujours à des solutions plus précises. Ce qui suit est un extrait de l'article fondateur de Prothero et Robinson de 1974:

En utilisant des méthodes en une étape stables pour résoudre de grands systèmes d'équations différentielles non linéaires rigides, nous avons constaté que
(a) certaines méthodes stables en A donnent des solutions très instables, et
(b) la précision des solutions obtenues lorsque les équations sont raide semble souvent sans rapport avec l'ordre de la méthode utilisée.

Pour des traitements encore plus rigoureux de ce sujet, voir le texte classique de Hairer & Wanner, "Résolution des équations différentielles ordinaires II: Rigides et différentielles - Problèmes algébriques", 1991.

En pratique, les équations rigides sont presque toujours résolues en utilisant la règle trapézoïdale ou la formule TR-BDF2 (fonctions ode23t et ode23tb dans MATLAB). Ces deux méthodes sont implicites de second ordre. Bien sûr, lorsque la stabilité n'est pas un problème (c'est-à-dire dans les équations non rigides), nous sommes libres de choisir parmi un certain nombre d'options; Le RK45 est le choix le plus courant.

La configuration de référence

Dans le logiciel Julia DifferentialEquations.jl, nous avons implémenté de nombreuses méthodes d'ordre supérieur, y compris les méthodes Feagin. Vous pouvez le voir dans notre liste de méthodes , puis il y a des tonnes d'autres que vous pouvez utiliser comme tableaux fournis . Étant donné que toutes ces méthodes sont regroupées, vous pouvez facilement les comparer. Vous pouvez voir les références que j'ai en ligne ici , et voir qu'il est très simple de comparer de nombreux algorithmes différents. Donc, si vous voulez prendre quelques minutes pour exécuter les benchmarks, allez-y. Voici un résumé de ce qui sort.

Tout d'abord, il est important de noter que, si vous regardez chacun des benchmarks, vous verrez que nos DP5(Ordre Dormand-Prince 5) et nos DP8méthodes sont plus rapides que les codes Hairer Fortran ( dopri5et dop853), et donc ces implémentations sont très bien optimisées . Ceux-ci montrent que, comme indiqué dans un autre fil, la surutilisation des méthodes Dormand-Prince est due au fait que les méthodes sont déjà écrites, et non pas parce qu'elles sont toujours les meilleures. La vraie comparaison entre les implémentations les plus optimisées est donc entre les méthodes Tsitorous, les méthodes Verner et les méthodes Feagin de DifferentialEquations.jl.

Les resultats

En général, les méthodes d'un ordre supérieur à 7 ont un coût de calcul supplémentaire qui n'est généralement pas compensé par l'ordre, compte tenu des tolérances choisies. Une des raisons à cela est que les choix de coefficients pour les méthodes d'ordre inférieur sont plus optimisés (ils ont de petits "coefficients d'erreur de troncature de principe", qui importent plus lorsque vous n'êtes pas asymétriquement petits). Vous pouvez voir que dans de nombreux problèmes comme ici, les méthodes Verner Efficient 6 et 7 fonctionnent très bien, mais des méthodes telles que Verner Efficient 8 peuvent avoir une pente plus faible. En effet, les «gains» d'ordre supérieur s'accumulent à des tolérances plus faibles, il y a donc toujours une tolérance où les méthodes d'ordre supérieur seront plus efficaces.

Cependant, la question est alors, à quel point? Dans une implémentation bien optimisée, cela devient assez bas pour deux raisons. La première raison est que les méthodes d'ordre inférieur implémentent quelque chose appelé FSAL (premier comme dernier). Cette propriété signifie que les méthodes d'ordre inférieur réutilisent une évaluation de fonction de l'étape précédente à l'étape suivante, et ont donc effectivement une évaluation de fonction en moins. Si cela est utilisé correctement, alors quelque chose comme une méthode du 5ème ordre (Tsitorous ou Dormand-Prince) prend en fait 5 évaluations de fonctions au lieu des 6 que les tableaux suggèrent. Cela est également vrai pour la méthode Verner 6.

L'autre raison est due aux interpolations. Une raison d'utiliser une méthode d'ordre très élevé est de prendre moins de pas et d'interpoler simplement des valeurs intermédiaires. Cependant, afin d'obtenir les valeurs intermédiaires, la fonction d'interpolation peut nécessiter plus d'évaluations de fonction que celles utilisées pour effectuer l'étape. Si vous regardez les méthodes Verner, il faut 8 évaluations de fonctions supplémentaires pour la méthode Order 8 pour obtenir un interpolant Order 8. Plusieurs fois, les méthodes de faible ordre fournissent un interpolant "libre", par exemple la plupart des méthodes de 5e ordre ont une interpolation de 4e ordre gratuite (pas d'évaluation de fonction supplémentaire). Cela signifie donc que si vous avez besoin de valeurs intermédiaires (dont vous aurez besoin pour un bon tracé si vous utilisez une méthode d'ordre élevé), il y a des coûts cachés supplémentaires. Tenez compte du fait que ces valeurs interpolées sont vraiment importantes pour la gestion des événements et la résolution des équations différentielles de retard et vous voyez pourquoi les coûts d'interpolation supplémentaires sont pris en compte.

Qu'en est-il des méthodes Feagin?

Vous verrez donc que les méthodes Feagin manquent étrangement aux benchmarks. Ils sont bons, les tests de convergence fonctionnent sur des nombres de précision arbitraires, etc., mais pour les faire bien fonctionner, vous devez demander des tolérances assez absurdes. Par exemple, j'ai trouvé dans des benchmarks non publiés que les Feagin14surperformances Vern9(la méthode efficace Verner du 9ème ordre) à des tolérances comme 1e-30. Pour les applications avec une dynamique chaotique (comme dans les Pléides ou les problèmes d'astrophysique à 3 corps), vous pouvez souhaiter cette quantité de précision en raison de la dépendance sensible (les erreurs dans les systèmes chaotiques se composent rapidement). Cependant, la plupart des gens calculent probablement avec des nombres à virgule flottante double précision, et je n'ai pas trouvé de référence où ils surpassent dans ce domaine de tolérance.

De plus, il n'y a pas d'interpolant pour accompagner les méthodes de Feagin. Donc, ce que je fais, c'est simplement de mettre une interpolation Hermite du troisième ordre sur eux afin que cette façon existe (et cela fonctionne étonnamment bien). Cependant, s'il n'y a pas de fonction d'interpolation standard, vous pouvez utiliser la méthode Hermite récursive (utilisez cette interpolation pour obtenir le point médian, puis effectuez une interpolation du 5e ordre, etc.) pour obtenir une interpolation d'ordre élevé, mais cela est très coûteux et le résultat l'interpolation n'a pas nécessairement un terme d'erreur de troncature de principe bas (donc c'est seulement bon quand dtc'est vraiment petit, ce qui est exactement l'opposé du cas que nous voulons!). Donc, si jamais vous avez besoin d'une très bonne interpolation pour correspondre à votre précision, vous devez au moins revenir à quelque chose comme Vern9.

Remarque sur l'extrapolation

Notez que les méthodes d'extrapolation sont simplement des algorithmes pour générer des méthodes Runge-Kutta d'ordre arbitraire. Cependant, pour leur ordre, ils prennent plus d'étapes que nécessaire et ont des coefficients d'erreur de troncature de principe élevé, et ils ne sont donc pas aussi efficaces qu'une méthode RK bien optimisée à un ordre donné. Mais compte tenu de l'analyse précédente, cela signifie qu'il existe un domaine de tolérance extrêmement faible où ces méthodes feront mieux que les méthodes RK "connues". Mais dans chaque benchmark que j'ai couru, il semble que je ne sois pas descendu aussi bas.

Remarque sur la stabilité

Le choix n'a vraiment rien à voir avec les problèmes de stabilité. En fait, si vous parcourez les tableaux DifferentialEquations.jl (vous pouvez juste plot(tab)pour les régions de stabilité), vous verrez que la plupart des méthodes ont des régions de stabilité étrangement similaires. C'est en fait un choix. Habituellement, lors de la dérivation des méthodes, l'auteur effectue généralement les opérations suivantes:

1. Trouver les coefficients d'erreur de troncature de principe les plus bas (c'est-à-dire les coefficients pour les termes d'ordre suivants)
2. Sous réserve des contraintes de commande
3. Et rendez la région de stabilité proche de celle de la méthode Dormand-Prince Order 5.

Pourquoi la dernière condition? Eh bien, parce que cette méthode a tendance à être toujours stable avec la façon dont les choix de taille adaptatifs contrôlés par PI sont effectués, c'est donc une bonne barre pour les régions de stabilité "assez bonnes". Ce n'est donc pas un hasard si les régions de stabilité ont toutes tendance à être similaires.

Conclusion

Il y a des compromis dans chaque choix de méthode. Les méthodes RK d'ordre le plus élevé ne sont tout simplement pas aussi efficaces à des tolérances inférieures à la fois parce qu'il est plus difficile d'optimiser le choix des coefficients, et parce que le nombre d'évaluations de fonctions se compose (et croît encore plus rapidement lorsque des interpolations sont impliquées). Cependant, si la tolérance devient suffisamment faible, ils l'emportent, mais les tolérances requises peuvent être bien inférieures aux applications "standard" (c'est-à-dire qu'elles ne s'appliquent vraiment qu'aux systèmes chaotiques).