Aller au contenu

Schéma de Runge-Kutta d'ordre 4

Le schéma de Runge-Kutta d'ordre 4 est basée sur une approximation de la dérivée à un ordre supérieur (ordre 4). Le principe de la discrétisation est le même que pour la méthode d'Euler, mais on va faire quelques calculs supplémentaires pour approcher la dérivée:

  • on connait la fonction f, un point ti où on connait yi
  • on peut donc calculer f1=yi=f(yi,ti) (cf. méthode d'Euler=pente au point (ti,yi))
  • on calcule f2=f(yi+12Δtf1,ti+12Δt) (valeur estimée au milieu de l'intervalle, avec la pente prise en ti)
  • on calcule f3=f(yi+12Δtf2,ti+12Δt) (valeur estimée au milieu de l'intervalle ti+1/2, avec la pente prise en ti+1/2)
  • on calcule f4=f(yi+Δtf3,ti+Δt) (valeur estimée en ti+1, avec la pente prise en ti+1/2)
  • on a alors yi+1y1+16Δt(f1+2f2+2f3+f4)
  • on peut alors itérer (résoudre pas à pas) pour passer au point suivant. Le problème est initialisé en partant de t0 où on connait y0 (condition à la limite).

On voit clairement que la méthode est beaucoup plus précise. Même avec un pas de calcul beaucoup plus élevé, on approche correctement la solution, alors que la méthode d'Euler donne des résultats très éloignés de la solution exacte. On remarque aussi que, entre les points de discrétisation, on a approché la solution par des segments de droite (bien qu'on puisse avoir une interpolation plus évoluée si on a besoin de connaître des valeurs intermédiaires).

Le programme peut s'écrire de la façon suivante, avec le logiciel de calcul scientifique Scilab:

// Programme de resolution d'une equation differentielle
// Equation a resoudre: dy/dt=f(y,t)

// Definition de la fonction f
a=-0.1;
function z=f(y,t);z=a*y; endfunction

// Condition à la limite:
t0=0;
y0=1;

// Discretisation du temps
tmax=50;
dt=5;

// Nombre de pas de discretisation
N=(tmax-t0)/dt;
// Indices
ii=1:N+1;
t=(ii-1)*dt; // vecteur temps

// vecteur temps avec un pas plus fin, pour la solution exacte
t2=0:tmax;

// Solution par méthode d'Euler, notée ye
ye(1)=y0; // condition à la limite
for i=1:N
  ye(i+1)=ye(i)+dt*f(ye(i),t(i));
end

// Solution par méthode RK4, notée yrk4
yrk4(1)=y0; // condition à la limite
for i=1:N
  f1=f(yrk4(i),t(i));
  f2=f(yrk4(i)+f1*dt/2,t(i)+dt/2);
  f3=f(yrk4(i)+f2*dt/2,t(i)+dt/2);
  f4=f(yrk4(i)+f3*dt,t(i)+dt);
  yrk4(i+1)=yrk4(i)+dt*(f1+2*f2+2*f3+f4)/6;
end

// Tracé des solutions
scf(2)
plot(t2,exp(a*t2),'k-',t,ye,'kd-',t,yrk4,'ks-')
title('Resolution par méthode de Runge-Kutta d''ordre 4 - dy/dt=-0.1y')
xlabel('t')
ylabel('y')
legend('Solution exacte','Méthode Euler, dt=5','Méthode RK4, dt=5')