Retour
Dominique Lefebvre Mars 2006
La méthode de Runge-Kutta d'ordre 4
L'algorithme de Runge-Kutta d'ordre 4
L'algorithme de Runge-Kutta d'ordre 4
Les méthodes de Runge-Kutta (ou RK), d'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 autre chose 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 : yn+1 = yn + h*f(xn, yn), et xn+1 = xn + h
La méthode RK du deuxième ordre produit deux coefficients k1 et k2, qui permettent d'écrire (voir démo. dans le cours d'AN):
k1 = h*f(xn, yn)
k2 = h*f(xn + h/2, yn+ k1/2 )
yn+1 = yn + k2 + O(h3)
Cette méthode exige donc deux évaluations de f. L'erreur de consistance est en O(h3) et l'erreur globale de convergence est d'ordre O(h2)
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:
k1 = h*f(xn, yn)
k2 = h*f(xn + h/2, yn+ k1/2 )
k3 = h*f(xn + h/2, yn+ k2/2 )
k4 = h*f(xn + h, yn+ k3)
yn+1 = yn + k1/6 + k2/3 + k3/3 + k4/6 + O(h5)
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(h5) et l'erreur globale de convergence est d'ordre O(h4).
Voilà très rapidement pour la "théorie".
HOME
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 suffisament 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 lesC 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
C Typage de l'interface de ResolutionRK4
C
INTEGER NbEquation
REAL x, y(NbEquation), dx
C
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.
HOME
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 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
HOME