test.c

math library from gnu

    }

  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]");
}

/* Test cases from Frank Reininghaus <[e-mail]> */

int rhs_stepfn (double t, const double * y, double * dydt, void * params) {
  if (t >= 1.0)
    dydt [0] = 1;
  else
    dydt [0] = 0;

  return GSL_SUCCESS;
}

void test_stepfn (void) {
  /* infinite loop for epsabs = 1e-18, but not for 1e-17 */
  double epsabs = 1e-18;
  double epsrel = 1e-6;

  const gsl_odeiv_step_type * T = gsl_odeiv_step_rk2;
  gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 1);
  gsl_odeiv_control * c = gsl_odeiv_control_y_new (epsabs, epsrel);
  gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (1);
  gsl_odeiv_system sys = {rhs_stepfn, 0, 1, 0};
     
  double t = 0.0;
  double h = 1e-6;
  double y = 0.0;
  int i = 0;
  int status;

  while (t < 2.0 && i < 1000000) {
    status = gsl_odeiv_evolve_apply (e, c, s, &sys, &t, 2, &h, &y);
#ifdef DEBUG
    printf("i=%d status=%d t=%g h=%g y=%g\n", i, status, t, h, y);
#endif
    if (status != GSL_SUCCESS)
      break;
    
    i++;
  }

  gsl_test_abs(t, 2.0, 1e-16, "evolve step function, t (stepfn/rk2)");
  gsl_test_rel(y, 1.0, epsrel, "evolve step function, y (stepfn/rk2)");
       
  gsl_odeiv_evolve_free (e);
  gsl_odeiv_control_free (c);
  gsl_odeiv_step_free (s);
}

int rhs_stepfn2 (double t, const double * y, double * dydt, void * params) {
  if (t >= 0.0)
    dydt [0] = 1e300;
  else
    dydt [0] = 0;

  return GSL_SUCCESS;
}

void test_stepfn2 (void) {
  /* infinite loop for epsabs = 1e-25, but not for 1e-24 */
  double epsabs = 1e-25;
  double epsrel = 1e-6;

  const gsl_odeiv_step_type * T = gsl_odeiv_step_rk2;
  gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 1);
  gsl_odeiv_control * c = gsl_odeiv_control_y_new (epsabs, epsrel);
  gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (1);
  gsl_odeiv_system sys = {rhs_stepfn2, 0, 1, 0};
     
  double t = -1.0;
  double h = 1e-6;
  double y = 0.0;

  int i = 0;
  int status;

  while (t < 1.0 && i < 10000) {
    status = gsl_odeiv_evolve_apply (e, c, s, &sys, &t, 1.0, &h, &y);
#ifdef DEBUG
    printf("i=%d status=%d t=%g h=%g y=%g\n", i, status, t, h, y);
#endif
    if (status != GSL_SUCCESS)
      break;

    i++;
  }

  gsl_test_abs(t, 1.0, 1e-16, "evolve big step function, t (stepfn2/rk2)");
  gsl_test_rel(y, 1e300, epsrel, "evolve big step function, y (stepfn2/rk2)");
     
  gsl_odeiv_evolve_free (e);
  gsl_odeiv_control_free (c);
  gsl_odeiv_step_free (s);
}

int rhs_stepfn3 (double t, const double * y, double * dydt, void * params) {

  static int calls = 0;

  if (t >= 0.0)
    dydt [0] = 1e300;
  else
    dydt [0] = 0;

  calls++;

  return (calls < 100000) ? GSL_SUCCESS : -999;
}


void test_stepfn3 (void) {
  /* infinite loop for epsabs = 1e-26, but not for 1e-25 */
  double epsabs = 1e-26;
  double epsrel = 1e-6;

  const gsl_odeiv_step_type * T = gsl_odeiv_step_rkf45;
  gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 1);
  gsl_odeiv_control * c = gsl_odeiv_control_y_new (epsabs, epsrel);
  gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (1);
  gsl_odeiv_system sys = {rhs_stepfn3, 0, 1, 0};
     
  double t = -1.0;
  double h = 1e-6;
  double y = 0.0;

  int i = 0;
  int status;

  while (t < 1.0 && i < 10000) {
    status = gsl_odeiv_evolve_apply (e, c, s, &sys, &t, 1.0, &h, &y);
#ifdef DEBUG
    printf("i=%d status=%d t=%g h=%g y=%g\n", i, status, t, h, y);
#endif
    if (status != GSL_SUCCESS)
      break;

    i++;
  }

  gsl_test_abs(t, 1.0, 1e-16, "evolve big step function, t (stepfn3/rkf45)");
  gsl_test_rel(y, 1e300, epsrel, "evolve big step function, y (stepfn3/rkf45)");
     
  gsl_odeiv_evolve_free (e);
  gsl_odeiv_control_free (c);
  gsl_odeiv_step_free (s);
}

int rhs_cos (double t, const double * y, double * dydt, void * params) {
  dydt [0] = cos (t);
  return GSL_SUCCESS;
}

int jac_cos (double t, const double y[], double *dfdy, double dfdt[],
             void *params)
{
  dfdy[0] = 0.0;
  dfdt[0] = -sin(t);

  return GSL_SUCCESS;
}

/* Test evolution in negative direction */

void
test_evolve_negative_h (const gsl_odeiv_step_type * T, double h, double err)
{ 
  /* Tolerance factor in testing errors */
  const double factor = 10;

  gsl_odeiv_step * step = gsl_odeiv_step_alloc (T, 1);
  gsl_odeiv_control * c = gsl_odeiv_control_standard_new (err, err, 1.0, 0.0);
  gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (1);
  gsl_odeiv_system sys = {rhs_cos, jac_cos, 1, 0};

  double t = 0;
  double t1 = -4.0;

  double y = 0.0;
  double yfin = sin (t1);

  /* Make initial h negative */
  h = -fabs(h);

  while (t > t1) {
    int status = gsl_odeiv_evolve_apply (e, c, step, &sys, &t, t1, &h, &y);
    
    if (status != GSL_SUCCESS) 
      {
	gsl_test(status, "%s evolve_apply returned %d for negative h",
		 gsl_odeiv_step_name (step), status);
	break;
      }
  }

  gsl_test_abs (y, yfin, factor * e->count * err,
		"evolution with negative h (using %s)", 
                gsl_odeiv_step_name (step));

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

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_evolve_negative_h (p[i].type, p[i].h, 1e-7);
    }

  test_compare_vanderpol();
  test_compare_oregonator();
  
  test_stepfn();
  test_stepfn2();
  test_stepfn3();
 
  exit (gsl_test_summary ());
}