TangenteX.com

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, qui est une méthode RK1.

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.
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 :
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 en FORTRAN

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 en FORTRAN

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

Implémentation de l'algorithme RK4 en C

L'implémentation en C de la méthode RK4 est assez semblable à celle en FORTRAN. Voici un exemple de code que j'utilise :

void RK4(double yin[], double yout[], double h)
{
   int i;
   double dydx[N], dydxt[N], yt[N], k1[N], k2[N], k3[N], k4[N];

   Derivee(yin, dydx);
   for (i = 0; i < N; i++)
    {
      k1[i] = h*dydx[i];
      yt[i] = yin[i] + 0.5*k1[i];
    }

   Derivee(yt, dydxt);
   for (i = 0; i < N; i++)
    {
      k2[i] = h*dydxt[i];
      yt[i] = yin[i] + 0.5*k2[i];
    }

   Derivee(yt, dydxt);
   for (i = 0; i < N; i++)
    {
      k3[i] = h*dydxt[i];
      yt[i] = yin[i] + k3[i];
    }

   Derivee(yt, dydxt);
   for (i = 0; i < N; i++)
    {
      k4[i] = h*dydxt[i];
      yout[i] = yin[i] + k1[i]/6. + k2[i]/3. + k3[i]/3. + k4[i]/6.;
    }

  return;
}

La routine Derivee() définit l'équation ou le système d'équations différentielle(s), par exemple, pour l'EDO du pendule simple non linéaire :

void Derivee(double Y[], double dY[])
{
   dY[0] = Y[1];
   dY[1] = -omega2*sin(Y[0]);
   return;
}

Vous trouverez différents exemples d'utilisation de RK4 en C sur TangenteX, par exemple ici.

Implémentation de l'algorithme RK4 en Python

L'implémentation de la méthode RK4 en Python est très similaire à celle du FORTRAN ou du C. Voici un exemple d'implémentation de RK4, que j'utilise pour intégrer un système de deux EDO de premier ordre (le système qui définit le mouvement d'un pendule simple non linéaire dans mon exemple) :

# définition de la fonction RK4
def RK4Sys2(fonction, dom, xini, yini, h):
""" 27 juillet 2017
    @author: dominique lefebvre - tangentex.com
    Cette fonction intégre le système d'EDO d'ordre 2 avec les paramètres suivants:
    - fonction : désigne la fonction de définition du système
    - dom : désigne le domaine d'intégration du système
    - xini, yini : conditions initiales du système
    - h : pas d'intégration sur le domaine La fonction retourne deux listes contenant les (xi,yi)
          de la courbe intégrale
"""
    xi = [];yi = []
    r = array([xini,yini])
    for t in dom:
        xi.append(r[0])
        yi.append(r[1])
        k1 = h*fonction(r,t)
        k2 = h*fonction(r + k1/2.0, t + h/2.0)
        k3 = h*fonction(r + k2/2.0, t + h/2.0)
        k4 = h*fonction(r + k3, t + h)
        r += (k1 + 2.0*k2 + 2.0*k3 + k4)/6.0
    return (xi,yi)

Exemple d'application du RK4 en Python

Reprenons l'exemple du pendule simple non linéaire que nous allons cette fois résoudre avec un script Python. Comme en FORTRAN ou en C, je construis la fonction qui définit l'équation différentielle à intégrer. Elle est très similaire, mais Python est moins verbeux que FORTRAN ou C :

# Définition du système Pendule simple
def Pendule(dxdy,t):
    x,y = dxdy   
    return array([y,-omega0**2*sin(x)])

Puis je définie les paramètres physiques du système (ici la pulsation propre du pendule) et les conditions initiales du système (angle initale et vitesse angulaire initiale) :

# paramètres du système
omega0 = 1.0      # pulsation propre du pendule
# rangement des conditions initiales dans un tableau
theta0 = pi/6.0   # angle initial
dtheta0 = 0.00    # vitesse angulaire initiale

Je définie le domaine d'intégration :

# définition du domaine d'intégration
xinf = 0.0
xsup = 4.0*pi
dt = 0.01
nbpas = int((xsup - xinf)/dt)
domaine = arange(xinf,xsup,dt)

Puis j'intègre l'équation différentielle avec la fonction RK4 programmée ci-dessus :

# intégration du système par RK4
theta,dtheta = RK4Sys2(Pendule,domaine,theta0,dtheta0,dt)

Python me permet d'afficher simplement le résultat, sans passer par GnuPlot, avec les intructions suivantes :

# affichage de la courbe des résultats
plt.figure()
plt.grid(True)
plt.xlim(xinf,xsup)
plt.title("Pendule simple - RK4")
plt.xlabel('temps')
plt.ylabel('angle')
plt.plot(domaine,theta,'b')
plt.show()

Voici le résultat :

Pendule simple RK4 Python

Vous constatez que l'usage de Python est un peu plus simple que le code en FORTRAN... Le résultat obtenu est heureusement le même, mais les moyens à mettre en oeuvre sont beaucoup plus légers. C'est sans doute pour cela que pour les calculs pas trop lourds, Python a détroné tous les langages compilés, C ou FORTRAN. d'autant qu'il est possible d'appeler des routines C ou FORTRAN dans un script Python...

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