src/test/runge-kutta.c

Convergence of the Runge–Kutta solvers

We solve numerically dudt=tu with the initial condition u(0)=1. The solution is u(t)=et2/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);
    }
}
Error convergence for different orders

Error convergence for different orders