Skip to content

Runge-Kutta diagram of order 4

The Runge-Kutta scheme of order 4 is based on an approximation of the derivative to a higher order (order 4). The principle of discretization is the same as for the Euler method, but we will make some additional calculations to approach the derivative:

  • we know the function \(f\), a point \(t_i\) where we know \(y_i\)
  • we can therefore calculate \(f_1=y'_i=f(y_i,t_i)\) (cf. Euler=pente method to the point (\(t_i,y_i\)))
  • we calculate \(f_2=f(y_i+\frac12 \Delta t f_1,t_i+\frac12 \Delta t)\) (estimated value in the middle of the interval, with the slope taken in \(t_i\))
  • we calculate \(f_3=f(y_i+\frac12 \Delta t f_2,t_i+\frac12 \Delta t)\) (estimated value in the middle of the interval \(t_{i+1/2}\), with the slope taken in \(t_{i+1/2}\)
  • we calculate \(f_4=f(y_i+ \Delta t f_3,t_i+ \Delta t)\) (estimated value in \(t_{i+1}\), with the slope taken in \(t_{i+1/2}\)
  • we then have \(y_{i+1}\simeq y_1 + \frac16 \Delta t(f_1+2f_2+2f_3+f_4)\)
  • you can then iterate (solve step by step) to move to the next point. The problem is initialized from \(t_0\) where we know \(y_0\) (boundary condition).

It is clear that the method is much more precise. Even with a much higher calculation step, the solution is approached correctly, whereas the Euler method gives results very far from the exact solution. We also note that, between the discretization points, the solution was approached by line segments (although we may have more advanced interpolation if we need to know intermediate values).

The program can be written as follows, using the Scilab scientific calculation software:

// Program for solving a differential equation
// Equation to be solved: dy/dt=f(y,t)

// Definition of the function f
a=-0.1;
function z=f(y,t);z=a*y; endfunction

// Boundary condition:
t0=0;
y0=1;

// Time discretization
tmax=50;
dt=5;

// Number of discretization steps
N=(tmax-t0)/dt;
// Indices
ii=1:N+1;
t=(ii-1)*dt; // time vector

// time vector with a finer step, for the exact solution
t2=0:tmax;

// Solution by Euler method, noted ye
ye(1)=y0; // boundary condition
for i=1:N
  ye(i+1)=ye(i)+dt*f(ye(i),t(i));
end

// Solution by RK4 method, noted yrk4
yrk4(1)=y0; // boundary condition
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

Solution mapping //
scf(2)
plot(t2,exp(a*t2),'k-',t,ye,'kd-',t,yrk4,'ks-')
title('Resolution by Runge-Kutta method of'order 4 - dy/dt=-0.1y')
xlabel('t')
ylabel('y')
legend('Exact solution','Euler method, dt=5','RK4 method, dt=5')