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 $t_i$ où on connait $y_i$
on peut donc calculer $f_1=y'_i=f(y_i,t_i)$ (cf. méthode d'Euler=pente au point ( $t_i,y_i$ ))
on calcule $f_2=f(y_i+\frac12 \Delta t f_1,t_i+\frac12 \Delta t)$ (valeur estimée au milieu de l'intervalle, avec la pente prise en $t_i$ )
on calcule $f_3=f(y_i+\frac12 \Delta t f_2,t_i+\frac12 \Delta t)$ (valeur estimée au milieu de l'intervalle $t_{i+1/2}$ , avec la pente prise en $t_{i+1/2}$ )
on calcule $f_4=f(y_i+ \Delta t f_3,t_i+ \Delta t)$ (valeur estimée en $t_{i+1}$ , avec la pente prise en $t_{i+1/2}$ )
on a alors $y_{i+1}\simeq y_1 + \frac16 \Delta t(f_1+2f_2+2f_3+f_4)$
on peut alors itérer (résoudre pas à pas) pour passer au point suivant. Le problème est initialisé en partant de $t_0$ où on connait $y_0$ (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')