runge-kutta.h

Go to the documentation of this file.

/** @file runge-kutta.h
 */
/**
# Runge--Kutta time integrators
*/
 
static double update (scalar * ul, scalar * kl, double t, double dt,
          void (* Lu) (scalar * ul, double t, scalar * kl),
          scalar * dul, double w)
{
  scalar * u1l = list_clone (ul);
  for (int _i = 0; _i < _N; _i++) /* foreach */ {
    scalar u1, u, k;
    for (u1,u,k in u1l,ul,kl)
      u1[] = u[] + dt*k[];
  }
  Lu (u1l, t + dt, kl);
  for (int _i = 0; _i < _N; _i++) /* foreach */ {
    scalar du, k;
    for (du,k in dul,kl)
      du[] += w*k[];
  }
  delete (u1l), free (u1l);
  return w;
}
 
/**
The *runge_kutta()* function implements the classical first- (Euler),
second- and fourth-order Runge--Kutta time integrators for evolution
equations of the form
\f[
\frac{\partial\mathbf{u}}{\partial t} = L(\mathbf{u}, t)
\f]
with \f$\mathbf{u}\f$ a vector (i.e. list) of evolving fields and \f$L()\f$ a
generic, user-defined operator.
 
Given \f$\mathbf{u}\f$, the initial time *t*, a timestep *dt* and the
function \f$L()\f$ which should fill *kl* with the right-hand-side of the
evolution equation, the function below will return \f$\mathbf{u}\f$ at
time \f$t + dt\f$ using the Runge--Kutta scheme specified by *order*. */
 
void runge_kutta (scalar * ul, double t, double dt,
      void (* Lu) (scalar * ul, double t, scalar * kl),
      int order)
{
  scalar * dul = list_clone (ul);
  scalar * kl = list_clone (ul);
  Lu (ul, t, kl);
  for (int _i = 0; _i < _N; _i++) /* foreach */ {
    scalar du, k;
    for (du,k in dul,kl)
      du[] = k[];
  }
 
  double w = 1.;
  switch (order) {
  case 1: // Euler
    break;
  case 2:
    w += update (ul, kl, t, dt, Lu, dul, 1.);
    break;
  case 4:
    w += update (ul, kl, t, dt/2., Lu, dul, 2.);
    w += update (ul, kl, t, dt/2., Lu, dul, 2.);
    w += update (ul, kl, t, dt,    Lu, dul, 1.);
    break;
  default:
    assert (false); // not implemented
  }
  
  for (int _i = 0; _i < _N; _i++) /* foreach */ {
    scalar u, du;
    for (u,du in ul,dul)
      u[] += dt/w*du[];
  }
 
  delete (dul), free (dul);
  delete (kl), free (kl);
}