runge-kutta.c

Convergence of the Runge–Kutta solvers

Convergence of the Runge–Kutta solvers

We solve numerically \displaystyle \frac{du}{dt} = tu with the initial condition u(0)=1. The solution is u(t)=e^{t^2/2}.

#include "runge-kutta.h"

This function returns the right-hand-side of the equation i.e. tu.

static void du (scalar * ul, double t, scalar * kl)
{
  scalar u = ul[0], k = kl[0];
  foreach()
    k[] = t*u[];
}

int main()
{
  init_grid (1);

  for (int order = 1; order <= 4; order *= 2)
    for (double dt = 1e-2; dt <= 8e-2; dt *= 2) {
      scalar u[];
      foreach()
	u[] = 1.; // the initial condition
      double emax = 0.;
      for (t = 0; t <= 2.; t += dt) {
	foreach() {
	  double e = fabs (u[] - exp(t*t/2.));
	  if (e > emax)
	    emax = e;
	  printf ("%g %g %g\n", t, u[], u[] - exp(t*t/2.));
	}
	runge_kutta ({u}, t, dt, du, order);
      }
      printf ("\n");
      fprintf (stderr, "%g %g %d\n", dt, emax, order);
    }
}

set xlabel 'dt'
set ylabel 'error'
set logscale
set key top left
set ytics format '%.0e'
plot "< grep '1$' log" pt 7 t '', 15.*x t '15 dt',  \
     "< grep '2$' log" pt 7 t '', 4.*x*x t '4 dt^2', \
     "< grep '4$' log" pt 7 t '', x**4/2. t 'dt^4/2'

Error convergence for different orders (script)