# 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'