fixed classic rk4 step actually returns result for two half steps

[Krishna Myneni]

The function gsl_odeiv2_driver_apply_fixed_step() function in GSL 2.5,
using the rk4 step method actually returns the result for two half-
steps instead of the result for a single step. The attached code shows
the result for integrating the Lorenz system of ODEs with given
parameters and initial conditions, after 1 step at dt = 1e-4:

$ ./bench_gsl_rk4
1.59864268529822870e-03  9.99903676398808394e-01  7.99388902928489989e-
08

The above result for y[0] only has 12 significant digits in comparison
to the double precision result for a single step of the classical RK4
estimate for the step size used. Of course it is more accurate, but the
point is that it is not giving the expected result for 1 classical RK4
step. The problem appears to be with rk4_apply() which is using both
single step and two half steps for estimating the error. It leaves the
output state to be the result of the two half steps. 

Perhaps this is by design, but it will cause confusion when comparing
with other rk4 fixed step solutions.

Krishna Myneni


/*
 * Sample ODE solver using the GSL (GNU Scientific Library)
 * (Packages gsl and gsl-devel must be installed)
 *
 * This example is modified to solve the Lorenz equations with
 * a fixed step solver, from the example at following web page:
 * 
 * https://www.gnu.org/software/gsl/doc/html/ode-initval.html
 *
 * Lorenz Equations:
 *
 *   dx/dt = sig*(y - x)
 *   dy/dt = (r - z)*x - y
 *   dz/dt = x*y - bp*z
 *
 * Build:
 *
 *     gcc -o bench_gsl_rk4 -lgsl -lgslcblas -lm bench_gsl_rk4.c
 *
*/

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

int func (double t, const double y[], double f[],
      void *params)
{
  double* p = (double *) params;
  (void)(t); /* avoid unused parameter warning */
  double sig = p[0];
  double r =   p[1];
  double bp =  p[2];

  f[0] = sig*(y[1] - y[0]);
  f[1] = (r - y[2])*y[0] - y[1];
  f[2] = y[0]*y[1] - bp*y[2];
  return GSL_SUCCESS;
}

int jac (double t, const double y[], double *dfdy,
     double dfdt[], void *params)
{
  (void)(t); /* avoid unused parameter warning */
  double* p = (double*) params;
  double sig = p[0];
  double r   = p[1];
  double bp  = p[2];
  gsl_matrix_view dfdy_mat
    = gsl_matrix_view_array (dfdy, 3, 3);
  gsl_matrix * m = &dfdy_mat.matrix;
  gsl_matrix_set (m, 0, 0, -sig);
  gsl_matrix_set (m, 0, 1, sig);
  gsl_matrix_set (m, 0, 2, 0.0);
  gsl_matrix_set (m, 1, 0, r - y[2]);
  gsl_matrix_set (m, 1, 1, -1.0);
  gsl_matrix_set (m, 1, 2, -y[0]);
  gsl_matrix_set (m, 2, 0, y[1]);
  gsl_matrix_set (m, 2, 1, y[0]);
  gsl_matrix_set (m, 2, 2, -bp);
  dfdt[0] = sig;
  dfdt[1] = -1.0;
  dfdt[2] = 0.0;
  return GSL_SUCCESS;
}

int main (void)
{
  int nsteps = 1;
  double dt = 1e-4;
  double p[3];
  p[0] = 16.0 ; // sig
  p[1] = 45.92; // r
  p[2] = 4.0;   // bp

  gsl_odeiv2_system sys = { func, jac, 3, p };

  gsl_odeiv2_driver *d =
    gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk4,
                                   dt, 1e-11, 1e-11);

  double t = 0.0;
  double y[3] = { 0.0, 1.0, 0.0 };
  int i, j, s;

    for (i = 0; i < nsteps; i++)
    {
      s = gsl_odeiv2_driver_apply_fixed_step (d, &t, dt, 1, y);

      if (s != GSL_SUCCESS)
        {
          printf ("error: driver returned %d\n", s);
          break;
        }
    }

    /* print final values of solution */
    printf ("%.17e %.17e %.17e\n", y[0], y[1], y[2]);

  gsl_odeiv2_driver_free (d);
  return s;
}

signature.asc
Description: This is a digitally signed message part