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 : 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).

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

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:

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:

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:

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:

PenduleRK4.prj

RoutineRK4.for

SysPendule.for

penduleRK4.for

PlotPendule.p

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


Contenu et design par Dominique Lefebvre - www.tangenteX.com mars 2013   Licence Creative Commons   Contact : PhysiqueX ou

Cette œuvre est mise à disposition selon les termes de la Licence Creative Commons Attribution - Pas d’Utilisation Commerciale - Pas de Modification 3.0 France.