RK4.h

Runge-Kutta 4th order

Integration between t and t+dt of the state vector Y using the function f\_rk4 : \frac{\text{d} Y}{\text{d} t} = f\_rk4(t,Y).

// fixme: rename the structure my_vec into state_vec ?
typedef struct { double y0, y1;}  my_vec;   // vector (y,y') for RK4 integration
// why typedef struct and not just struct ? 
// It allows to wirte   my_vec Y;   instead of  struct my_vec Y;

my_vec RK4 (double t, double dt, my_vec Y, my_vec (* f_rk4) (double t, my_vec Y) ){

  my_vec k1 = f_rk4(t, Y);
  my_vec Y2;
  Y2.y0 = Y.y0+dt/2.*k1.y0;
  Y2.y1 = Y.y1+dt/2.*k1.y1;
  my_vec k2 = f_rk4(t+dt/2., Y2);
  Y2.y0 = Y.y0+dt/2.*k2.y0;
  Y2.y1 = Y.y1+dt/2.*k2.y1;
  my_vec k3 = f_rk4(t+dt/2., Y2);
  Y2.y0 = Y.y0+dt*k3.y0;
  Y2.y1 = Y.y1+dt*k3.y1;
  my_vec k4 = f_rk4(t+dt/2., Y2);

  Y2.y0 = Y.y0 + dt/6.*(k1.y0 + 2.*k2.y0 + 2.*k3.y0 + k4.y0);
  Y2.y1 = Y.y1 + dt/6.*(k1.y1 + 2.*k2.y1 + 2.*k3.y1 + k4.y1);
  
  return Y2;
}

Euler 1st order

my_vec Euler (double t, double dt, my_vec Y, my_vec (* f_rk4) (double t, my_vec Y) ){

  my_vec k1 = f_rk4(t, Y);
  my_vec Y2;
  Y2.y0 = Y.y0+dt*k1.y0;
  Y2.y1 = Y.y1+dt*k1.y1;
    
  return Y2;
}