test.c

This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY without ev

  double err_target;

  /* linear */
  h = 1e-1;
  err_target = 1e-10;
  y[0] = 0.58;
  yfin[0] = y[0] + 2 * h;
  test_odeiv_stepper (T, &rhs_func_lin, h, 0.0, "linear",
		      y, yfin, err_target);
  
  /* exponential */
  h = 1e-4;
  err_target = 1e-8;
  y[0] = exp(2.7);
  yfin[0] = exp(2.7 + h);
  test_odeiv_stepper (T, &rhs_func_exp, h, 2.7, "exponential",
		      y, yfin, err_target);
  /* cosine-sine */
  h = 1e-3;
  err_target = 1e-6;
  y[0] = cos(1.2);
  y[1] = sin(1.2);
  yfin[0] = cos(1.2 + h);
  yfin[1] = sin(1.2 + h);
  test_odeiv_stepper (T, &rhs_func_sin, h, 1.2, "cosine-sine",
		      y, yfin, err_target);

  /* classic stiff */
  h = 1e-7;
  err_target = 1e-4;
  y[0] = 1.0;
  y[1] = 0.0;

  {
    const double e1 = exp (-h);
    const double e2 = exp (-1000.0 * h);
    yfin[0] = 2.0 * e1 - e2;
    yfin[1] = -e1 + e2;
  }

  test_odeiv_stepper (T, &rhs_func_stiff, h, 0.0, "classic_stiff",
		      y, yfin, err_target);  
}

void
test_evolve_system (const gsl_odeiv_step_type * T,
                    const gsl_odeiv_system * sys,
                    double t0, double t1, double hstart,
                    double y[], double yfin[],
                    double err_target, const char *desc)
{
  /* Tests system sys with stepper T. Step length is controlled by
     error estimation from the stepper.
  */     
  
  int steps = 0;
  size_t i;

  double t = t0;
  double h = hstart;

  /* Tolerance factor in testing errors */
  const double factor = 10;

  gsl_odeiv_step * step = gsl_odeiv_step_alloc (T, sys->dimension);

  gsl_odeiv_control *c =
    gsl_odeiv_control_standard_new (err_target, err_target, 1.0, 0.0);

  gsl_odeiv_evolve *e = gsl_odeiv_evolve_alloc (sys->dimension);

  while (t < t1)
    {
      int s = gsl_odeiv_evolve_apply (e, c, step, sys, &t, t1, &h, y);

      if (s != GSL_SUCCESS && sys != &rhs_func_xsin) 
	{
	  gsl_test(s, "%s evolve_apply returned %d",
		   gsl_odeiv_step_name (step), s);
	  break;
	}

      if (steps > 100000)
	{
	  gsl_test(GSL_EFAILED, 
		   "%s evolve_apply reached maxiter",
		   gsl_odeiv_step_name (step));
	  break;
	}

      steps++;
    }

  /* err_target is target error of one step. Test if stepper has made
     larger error than (tolerance factor times) the number of steps
     times the err_target */

  for (i = 0; i < sys->dimension; i++)
    {
      gsl_test_abs (y[i], yfin[i], factor * e->count * err_target,
		    "%s %s evolve(%d)",
		    gsl_odeiv_step_name (step), desc, i);
    }

  gsl_odeiv_evolve_free (e);
  gsl_odeiv_control_free (c);
  gsl_odeiv_step_free (step);
}

int
sys_driver (const gsl_odeiv_step_type * T,
	    const gsl_odeiv_system * sys,
	    double t0, double t1, double hstart,
	    double y[], double epsabs, double epsrel,
	    const char desc[])
{
  /* This function evolves a system sys with stepper T from t0 to t1.
     Step length is varied via error control with possibly different
     absolute and relative error tolerances.
  */
  
  int s = 0;
  int steps = 0;

  double t = t0;
  double h = hstart;

  gsl_odeiv_step * step = gsl_odeiv_step_alloc (T, sys->dimension);

  gsl_odeiv_control *c =
    gsl_odeiv_control_standard_new (epsabs, epsrel, 1.0, 0.0);
  gsl_odeiv_evolve *e = gsl_odeiv_evolve_alloc (sys->dimension);

  while (t < t1)
    {
      s = gsl_odeiv_evolve_apply (e, c, step, sys, &t, t1, &h, y);

      if (s != GSL_SUCCESS) 
	{
	  gsl_test(s, "sys_driver: %s evolve_apply returned %d",
		   gsl_odeiv_step_name (step), s);
	  break;
	}

      if (steps > 1e7)
	{
	  gsl_test(GSL_EMAXITER, 
		   "sys_driver: %s evolve_apply reached maxiter at t=%g",
		   gsl_odeiv_step_name (step), t);
	  s = GSL_EMAXITER;
	  break;
	}

      steps++;
    }

  gsl_test(s, "%s %s [%g,%g], %d steps completed", 
	   gsl_odeiv_step_name (step), desc, t0, t1, steps);

  gsl_odeiv_evolve_free (e);
  gsl_odeiv_control_free (c);
  gsl_odeiv_step_free (step);

  return s;
}

void
test_compare_vanderpol (void)
{
  /* Compares output of van Der Pol oscillator with several steppers */
  
  /* system dimension */
  const size_t sd = 2;
  
  const gsl_odeiv_step_type *steppers[20];
  const gsl_odeiv_step_type **T;

  /* Required error tolerance for each stepper. */
  double err_target[20];

  /* number of ODE solvers */
  const size_t ns = 11;

  /* initial values for each ode-solver */
  double y[11][2];
  double *yp = &y[0][0];

  size_t i, j, k;
  int status = 0;

  /* Parameters for the problem and stepper  */
  const double start = 0.0;
  const double end = 100.0;
  const double epsabs = 1e-8;
  const double epsrel = 1e-8;
  const double initstepsize = 1e-5;

  /* Initialize */

  steppers[0] = gsl_odeiv_step_rk2;
  err_target[0] = 1e-6;
  steppers[1] = gsl_odeiv_step_rk4;
  err_target[1] = 1e-6;
  steppers[2] = gsl_odeiv_step_rkf45;
  err_target[2] = 1e-6;
  steppers[3] = gsl_odeiv_step_rkck;
  err_target[3] = 1e-6;
  steppers[4] = gsl_odeiv_step_rk8pd;
  err_target[4] = 1e-6;
  steppers[5] = gsl_odeiv_step_rk2imp;
  err_target[5] = 1e-5;
  steppers[6] = gsl_odeiv_step_rk2simp;
  err_target[6] = 1e-5;
  steppers[7] = gsl_odeiv_step_rk4imp;
  err_target[7] = 1e-6;
  steppers[8] = gsl_odeiv_step_bsimp;
  err_target[8] = 1e-7;
  steppers[9] = gsl_odeiv_step_gear1;
  err_target[9] = 1e-2;
  steppers[10] = gsl_odeiv_step_gear2;
  err_target[10] = 1e-6;
  steppers[11] = 0;

  T = steppers;

  for (i = 0; i < ns; i++) 
    {
      y[i][0] = 1.0;
      y[i][1] = 0.0;
    }
  
  /* Call each solver for the problem */

  i = 0;
  while (*T != 0)
    {
      {
	int s = sys_driver (*T, &rhs_func_vanderpol,
			    start, end, initstepsize, &yp[i], 
			    epsabs, epsrel, "vanderpol");
	if (s != GSL_SUCCESS)
	  {
	    status++;
	  }
      }
      
      T++;
      i += sd;
    }

  if (status != GSL_SUCCESS)
    {
      return;
    }

  /* Compare results */
      
  T = steppers;

  for (i = 0; i < ns; i++)
    for (j = i+1; j < ns; j++)
      for (k = 0; k < sd; k++)
	{
	  const double val1 = yp[sd * i + k];
	  const double val2 = yp[sd * j + k];
	  gsl_test_abs (val1, val2, 
			( GSL_MAX(err_target[i], err_target[j]) ),
			"%s/%s vanderpol",
			T[i]->name, T[j]->name);
	}

}

void
test_compare_oregonator (void)
{
  /* Compares output of the Oregonator with several steppers */
  
  /* system dimension */
  const size_t sd = 3;
  
  const gsl_odeiv_step_type *steppers[20];
  const gsl_odeiv_step_type **T;

  /* Required error tolerance for each stepper. */
  double err_target[20];

  /* number of ODE solvers */
  const size_t ns = 2;

  /* initial values for each ode-solver */
  double y[2][3];
  double *yp = &y[0][0];

  size_t i, j, k;
  int status = 0;
  
  /* Parameters for the problem and stepper  */
  const double start = 0.0;
  const double end = 360.0;
  const double epsabs = 1e-8;
  const double epsrel = 1e-8;
  const double initstepsize = 1e-5;

  /* Initialize */

  steppers[0] = gsl_odeiv_step_rk2simp;
  err_target[0] = 1e-6;
  steppers[1] = gsl_odeiv_step_bsimp;
  err_target[1] = 1e-6;
  steppers[2] = 0;

  T = steppers;

  for (i = 0; i < ns; i++) 
    {
      y[i][0] = 1.0;
      y[i][1] = 2.0;
      y[i][2] = 3.0;
    }
  
  /* Call each solver for the problem */

  i = 0;
  while (*T != 0)
    {
      {
	int s = sys_driver (*T, &rhs_func_oregonator,
			    start, end, initstepsize, &yp[i], 
			    epsabs, epsrel, "oregonator");

	if (s != GSL_SUCCESS)
	  {
	    status++;
	  }
      }

      T++;
      i += sd;
    }

  if (status != GSL_SUCCESS)
    {
      return;
    }
      
  /* Compare results */
  
  T = steppers;

  for (i = 0; i < ns; i++)
    for (j = i+1; j < ns; j++)
      for (k = 0; k < sd; k++)
	{
	  const double val1 = yp[sd * i + k];
	  const double val2 = yp[sd * j + k];
	  gsl_test_rel (val1, val2, 
			( GSL_MAX(err_target[i], err_target[j]) ),
			"%s/%s oregonator",
			T[i]->name, T[j]->name);
	}

}

void
test_evolve_linear (const gsl_odeiv_step_type * T, double h, double err)
{
  double y[1];
  double yfin[1];

  y[0] = 1.0;
  yfin[0] = 9.0;
  test_evolve_system (T, &rhs_func_lin, 0.0, 4.0, h, y, yfin, err,
		      "linear[0,4]");
}

void
test_evolve_exp (const gsl_odeiv_step_type * T, double h, double err)
{
  double y[1];
  double yfin[1];

  y[0] = 1.0;
  yfin[0] = exp (2.0);
  test_evolve_system (T, &rhs_func_exp, 0.0, 2.0, h, y, yfin, err,
		      "exp[0,2]");
}

void
test_evolve_sin (const gsl_odeiv_step_type * T, double h, double err)
{
  double y[2];
  double yfin[2];

  y[0] = 1.0;
  y[1] = 0.0;
  yfin[0] = cos (2.0);
  yfin[1] = sin (2.0);
  test_evolve_system (T, &rhs_func_sin, 0.0, 2.0, h, y, yfin, err,
		      "sine[0,2]");
}

void
test_evolve_xsin (const gsl_odeiv_step_type * T, double h, double err)
{
  double y[2];
  double yfin[2];

  y[0] = 1.0;
  y[1] = 0.0;
  yfin[0] = cos (2.0);
  yfin[1] = sin (2.0);
  test_evolve_system (T, &rhs_func_xsin, 0.0, 2.0, h, y, yfin, err,
		      "sine[0,2] w/errors");
}


void
test_evolve_stiff1 (const gsl_odeiv_step_type * T, double h, double err)
{
  double y[2];
  double yfin[2];

  y[0] = 1.0;
  y[1] = 0.0;
  {
    double arg = 1.0;
    double e1 = exp (-arg);
    double e2 = exp (-1000.0 * arg);
    yfin[0] = 2.0 * e1 - e2;
    yfin[1] = -e1 + e2;
  }
  test_evolve_system (T, &rhs_func_stiff, 0.0, 1.0, h, y, yfin, err,
                      "stiff[0,1]");
}

void
test_evolve_stiff5 (const gsl_odeiv_step_type * T, double h, double err)
{
  double y[2];
  double yfin[2];

  y[0] = 1.0;
  y[1] = 0.0;
  {
    double arg = 5.0;
    double e1 = exp (-arg);
    double e2 = exp (-1000.0 * arg);
    yfin[0] = 2.0 * e1 - e2;
    yfin[1] = -e1 + e2;
  }
  test_evolve_system (T, &rhs_func_stiff, 0.0, 5.0, h, y, yfin, err,
                      "stiff[0,5]");
}


int
main (void)
{
  int i;

  struct ptype
  {
    const gsl_odeiv_step_type *type;
    double h;
  }
  p[20];

  p[0].type = gsl_odeiv_step_rk2;
  p[0].h = 1.0e-3;
  p[1].type = gsl_odeiv_step_rk4;
  p[1].h = 1.0e-3;
  p[2].type = gsl_odeiv_step_rkf45;
  p[2].h = 1.0e-3;
  p[3].type = gsl_odeiv_step_rkck;
  p[3].h = 1.0e-3;
  p[4].type = gsl_odeiv_step_rk8pd;
  p[4].h = 1.0e-3;
  p[5].type = gsl_odeiv_step_rk2imp;
  p[5].h = 1.0e-3;
  p[6].type = gsl_odeiv_step_rk2simp;
  p[6].h = 1.0e-3;
  p[7].type = gsl_odeiv_step_rk4imp;
  p[7].h = 1.0e-3;
  p[8].type = gsl_odeiv_step_bsimp;
  p[8].h = 1.0e-3;
  p[9].type = gsl_odeiv_step_gear1;
  p[9].h = 1.0e-3;
  p[10].type = gsl_odeiv_step_gear2;
  p[10].h = 1.0e-3;
  p[11].type = 0;

  gsl_ieee_env_setup ();

  for (i = 0; p[i].type != 0; i++)
    {
      test_stepper(p[i].type);
    }

  for (i = 0; p[i].type != 0; i++)
    {
      test_evolve_linear (p[i].type, p[i].h, 1e-10);
      test_evolve_exp (p[i].type, p[i].h, 1e-6);
      test_evolve_sin (p[i].type, p[i].h, 1e-8);
      test_evolve_xsin (p[i].type, p[i].h, 1e-8);
      test_evolve_stiff1 (p[i].type, p[i].h, 1e-7);
      test_evolve_stiff5 (p[i].type, p[i].h, 1e-7);
    }

  test_compare_vanderpol();
  test_compare_oregonator();

  exit (gsl_test_summary ());
}