La méthode de Runge-Kutta d'ordre 4

L'algorithme de Runge-Kutta d'ordre 4

Les méthodes de Runge-Kutta (ou RK), l'ordre 2 ou 4, sont très couramment utilisées pour la résolution d'équations différentielles ordinaires (EDO). Ce sont des méthodes à pas unique, directement dérivées de la méthode d'Euler. Elles ont l'avantage d'être simples à programmer et d'être assez stables pour les fonctions courantes de la physique. Sur le plan de l'analyse numérique, elles ont surtout l'immense avantage de ne pas nécessiter autrechose que la connaissance des valeurs initiales. Elles démarrent toutes seules.
Elles ont quand même un inconvénient, surtout la méthode RK d'ordre 4, de son petit nom RK4: elles sont assez consommatrices en temps de calcul. On peut donc les employer lorsque le temps de calcul n'est pas trop grand. Dans le cas contraire, il vaut mieux se tourner vers une méthode à prédicteur/correcteur (Adams par exemple). Si la précision requise est très importante, vous devrez choisir la méthode RK4 adaptative ou mieux encore vers la méthode de Bulirsch-Stoer.
Voici une description rapide de la méthode RK4. Je ne présente que les résultats. Si vous souhaitez une démonstration de la chose, vous pouvez vous reporter à votre cours d'analyse numérique. A défaut, un excellent bouquin d'introduction : "Introduction à l'analyse numérique - Applications sous Matlab" de Jérôme Bastien et Jean-Noël Martin chez DUNOD (pub gratuite...).
Donc, l'algorithme RK4:
On part de la formule d'Euler, qui donne : y_n+1 = y_n + h*f(x_n, y_n), et x_n+1 = x_n + h
La méthode RK du deuxième ordre produit deux coefficients k₁ et k₂, qui permettent d'écrire (voir démo. dans le cours d'AN):
k₁ = h*f(x_n, y_n)
k₂ = h*f(x_n + h/2, y_n + k₁/2 )
y_n+1 = y_n + k₂ + O(h³)

Cette méthode exige donc deux évaluations de f. L'erreur de consistance est en O(h³) et l'erreur globale de convergence est d'ordre O(h²)
Pour obtenir plus de précision, mais en doublant le temps de calcul puisqu'on procéde à 4 évaluations de f, voici la méthode RK4:

k₁ = h*f(x_n, y_n)

k₂ = h*f(x_n + h/2, y_n + k₁/2 )

k₃ = h*f(x_n + h/2, y_n + k₂/2 )

k₄ = h*f(x_n + h, y_n + k₃)

y_n+1 = y_n + k₁/6 + k₂/3 + k₃/3 + k₄/6 + O (h⁵)

La méthode RK4 exige donc 4 évaluations de f, ce qui peut être gênant si f est compliquée. L'erreur de consistance est en O(h⁵) et l'erreur globale de convergence est d'ordre O(h⁴).

Implémentation de l'algorithme RK4

L'implémentation en FORTRAN de l'algorithme RK4 décrit ci-dessus est pratiquement immédiate. Il s'agit d'écrire une routine qui soit suffisamment générale pour être réutilisable dans tous nos programme de physique numérique. Cela signifie qu'elle ne doit pas contenir de données propres aux programmes et que ses paramètres doivent être suffisamment complets pour supporter tous les échanges nécessaires entre le programme et la routine. Ce qui donne:

SUBROUTINE ResolutionRK4(x,y,dx,NbEquation,Derivee)

C x : abscisse

C y() : A l'appel ce tableau contient les valeurs y initiales. Au retour, il contient les

C valeurs calculees

C dx : pas de calcul

C NbEquation : nombre d'equations dans le systeme dif.

C Derivee : nom de la fonction decrivant le systeme dif.

IMPLICIT NONE

EXTERNAL Derivee ! fonction de calcul de la derivee

C Typage de l'interface de ResolutionRK4

INTEGER NbEquation

REAL x, y(NbEquation), dx

C Declaration des variables locales

REAL pred1(NbEquation),pred2(NbEquation),pred3(NbEquation),1 pred4(NbEquation), ytemp(NbEquation), halfdx

INTEGER i

halfdx = dx/2

C Premiere prediction

CALL Derivee(x,y,pred1,NbEquation)

DO i=1, NbEquation

ytemp(i) = y(i) + pred1(i)*halfdx

ENDDO

C Seconde prediction

CALL Derivee(x+halfdx,ytemp,pred2,NbEquation)

DO i=1, NbEquation

ytemp(i) = y(i) + pred2(i)*halfdx

ENDDO

C Troisieme prediction

CALL Derivee(x+halfdx,ytemp,pred3,NbEquation)

DO i=1, NbEquation

ytemp(i) = y(i) + pred3(i)*dx

ENDDO

C Quatrieme prediction

CALL Derivee(x+dx,ytemp,pred4,NbEquation)

C Estimation de y pour le pas suivant

DO i=1, NbEquation

y(i)=y(i)+(dx/6.0)*(pred1(i)+2.0*pred2(i)+2.0*pred3(i) + pred4(i))

ENDDO

END

Quelques commentaires:

Le vecteur y de l'interface de la routine est utilisé à la fois pour transmettre les valeurs initiales de y et pour retourner les valeurs calculées de y. On rencontre souvent des implémentations où le programmeur sépare les y d'entrée et les y de sortie. Cela est affaire d'habitude. Je considère que l'écriture ci-dessus est plus concise.
La subroutine Derivee, qui est un paramètre de la routine RK4, procède à l'évaluation de f. Lorsque vous écrirez un programme de physique numérique, c'est cette routine qu'il faudra adapter pour résoudre votre système.
Le programme colle à l'algorithme. Il n'y a pas de piège caché...

Cette routine pourra être employée sans précaution particulière d'ordre informatique dans tous les programmes.

Exemple d'application du RK4

Pour illustrer l'emploi de la méthode RK4, reprenons l'exemple de résolution du système du pendule simple, que nous avons résolu en utilisant la méthode d'Euler.
Le programme PenduleRK4 sera associé à deux subroutines Derivee et ResolutionRK4. Cette dernière est décrite ci-dessus. Le code source de la routine Derivee, qui décrit le système différentiel:

SUBROUTINE Derivee(x,y,dy, NbEquation)

IMPLICIT None

REAL l, g

COMMON /Pendule/ l

DATA g/9.81/

INTEGER nbEquation

REAL x, y(NbEquation), dy(NbEquation)

C Description du systeme differentiel de l'oscillateur

dy(1) = y(2)

dy(2) = -(g/l)*SIN(y(1))

RETURN

END

Quelques commentaires:

Notez le COMMON qui permet de partager des données, ici la longueur du pendule avec le programme principal, sans modifier l'interface de la routine.
La variable NbEquation vaut 2, car comme vous le constatez, le système différentiel comprend deux équations.

Le programme PenduleRK4 maintenant:

PROGRAM PenduleRK4

IMPLICIT None

C Declaration des routines externes

C ResolutionRK4 : routine d'implementation du schema RK4

C Derivee : description du systeme differentiel de l'oscillateur

EXTERNAL ResolutionRK4

EXTERNAL Derivee

C Declaration des parametres RK4

C Le systeme differentiel comporte deux equations

INTEGER NbEquation

PARAMETER (NbEquation = 2)

C Declaration des constantes

REAL g, PI

DATA g/9.81/, PI/3.141592654/

C Declaration des variables

INTEGER i, NbPas

REAL y(NbEquation), t, l, PasTemps, a0

C Declaration des variables communes

COMMON /Pendule/l

C Saisie des parametres du mouvement

WRITE(*,'(a,$)') 'Longueur du pendule (en m) : '

READ(*,*) l

WRITE(*,'(a,$)') 'angle initial du pendule (en degres) : '

READ(*,*) a0

WRITE(*,'(a,$)') 'Pas temporel de calcul (en s) : '

READ(*,*) PasTemps

WRITE(*,'(a,$)') 'Nombre de pas : '

READ(*,*) NbPas

C Declaration des conditions initiales pour l'integration

y(1) = a0*PI/180.

! angle initial du pendule

y(2) = 0.0

! vitesse angulaire initiale

C Ouverture du fichier de stockage des resultats de calcul

C Ce fichier sera trace par GnuPlot

OPEN(10, FILE='PenduleRK4.dat')

C Ecriture des conditions initiales

t = 0.0

WRITE(10,* ) t, y(1)

C Boucle de calcul - le coeur du programme

DO i=1, NbPas

t = i*PasTemps

CALL ResolutionRK4(t, y, PasTemps, NbEquation, Derivee)

WRITE(10,*) t, y(1)*180./PI

! notez la conversion en degrés

ENDDO

C Fermeture du fichier de donnees

CLOSE(10)

STOP

END

Quelques commentaires:

L'instruction EXTERNAL permet de dire au programme qu'il doit chercher ailleurs que dans son fichier source le code source des routines Derivee et Resolution RK4. Il ira donc les chercher dans d'autres fichiers sources du projet.
Notez le DATA qui permet d'initialiser des données, au lieu de l'usage du PARAMETER.
Le plus gros du programme consiste à déclarer les variables et à saisir les valeurs des données. Le coeur du programme ne fait que 5 lignes!

Pour faire fonctionner ce programme, il faut créer un projet PenduleRK4 avec le programme VFORT (voir Introduction au FORTRAN), puis créer les 3 fichiers sources PenduleRK4, ResolutionRK4 et Derivee.
Après exécution, on peut visualiser les résultats en utilisant le fichier de commandes GnuPlot suivant:

# Fichier de tracé du mouvement d'un pendule

#selon la methode RK4

set data style line

set title "Variation de l'angle en fonction du temps"

set xlabel "Temps"

set ylabel "Angle"

plot 'PenduleRK4.dat'

Pour télécharger le projet, les fichiers sources et le fichier de commandes GnuPlot:

Nota : le fichier source SysPendule.for contient la subroutine Derivee