/** # Convergence of the Runge--Kutta solvers We solve numerically $$\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); } } /** ~~~gnuplot Error convergence for different orders 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' ~~~ */