test.c

该文件为c++的数学函数库!是一个非常有用的编程工具.它含有各种数学函数,为科学计算、工程应用等程序编写提供方便！

  gsl_test (s, "%s, sine [0,pi], max relative error = %g",
            gsl_odeiv_step_name (stepper), delmax);

  delmax = 0.0;
  for (; t < 100.5 * M_PI; t += h)
    {
      gsl_odeiv_step_apply (stepper, t, h, y, yerr, 0, 0, &rhs_func_sin);
      del = fabs (y[1] - sin (t));
      delmax = GSL_MAX_DBL (del, delmax);
      count++;
    }

  if (del > count * 2.0 * base_prec)
    {
      s++;
      printf ("  SIN(%22.18g)  %22.18g  %22.18g  %10.6g\n", t + h, y[1],
              sin (t), del);
    }

  gsl_test (s, "%s, sine [pi,100.5*pi], max absolute error = %g",
            gsl_odeiv_step_name (stepper), delmax);

  gsl_odeiv_step_free (stepper);
}


void
test_stepper_exp (const gsl_odeiv_step_type * T, double h, double base_prec)
{
  int s = 0;
  double y[2];
  double yerr[2];
  double t;
  double del, delmax = 0.0;
  int count = 0;

  gsl_odeiv_step *stepper = gsl_odeiv_step_alloc (T, 2);

  y[0] = 1.0;
  y[1] = 1.0;

  for (t = 0.0; t < 20.0; t += h)
    {
      double ex = exp (t + h);
      gsl_odeiv_step_apply (stepper, t, h, y, yerr, 0, 0, &rhs_func_exp);
      del = fabs ((y[1] - ex) / y[1]);
      delmax = GSL_MAX_DBL (del, delmax);
      if (del > (count + 1.0) * 2.0 * base_prec)
        {
          printf ("  EXP(%20.17g)  %20.17g  %20.17g  %8.4g\n", t + h, y[1],
                  ex, del);
          s++;
        }
      count++;
    }

  gsl_test (s, "%s, exponential [0,20], max relative error = %g",
            gsl_odeiv_step_name (stepper), delmax);

  gsl_odeiv_step_free (stepper);
}

void
test_stepper_stiff (const gsl_odeiv_step_type * T, double h, double base_prec)
{
  int s = 0;
  double y[2];
  double yerr[2];
  double t;
  double del, delmax = 0.0;
  int count = 0;

  gsl_odeiv_step *stepper = gsl_odeiv_step_alloc (T, 2);

  y[0] = 1.0;
  y[1] = 0.0;

  for (t = 0.0; t < 20.0; t += h)
    {
      gsl_odeiv_step_apply (stepper, t, h, y, yerr, NULL, NULL,
                            &rhs_func_stiff);

      if (t > 0.04)
        {
          double arg = t + h;
          double e1 = exp (-arg);
          double e2 = exp (-1000.0 * arg);
          double u = 2.0 * e1 - e2;
          /* double v = -e1 + e2; */
          del = fabs ((y[0] - u) / y[0]);
          delmax = GSL_MAX_DBL (del, delmax);

          if (del > (count + 1.0) * 100.0 * base_prec)
            {
              printf ("  STIFF(%20.17g)  %20.17g  %20.17g  %8.4g\n", arg,
                      y[0], u, del);
              s++;
            }
        }
      count++;
    }

  gsl_test (s, "%s, stiff [0,20], max relative error = %g",
            gsl_odeiv_step_name (stepper), delmax);

  gsl_odeiv_step_free (stepper);
}

void
test_evolve_system_flat (gsl_odeiv_step * step,
                         const gsl_odeiv_system * sys,
                         double t0, double t1, double hstart,
                         double y[], double yfin[],
                         double err_target, const char *desc)
{
  int s = 0;
  double frac;

  double t = t0;
  double h = hstart;

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

  while (t < t1)
    {
      gsl_odeiv_evolve_apply (e, NULL, step, sys, &t, t1, &h, y);
    }

  frac = fabs ((y[1] - yfin[1]) / yfin[1]) + fabs ((y[0] - yfin[0]) / yfin[0]);

  if (frac > 2.0 * e->count * err_target)
    {
      printf ("FLAT t = %.5e  y0 = %g y1= %g\n", t, y[0], y[1]);
      s++;
    }

  gsl_test (s, "%s, %s, evolve, no control, max relative error = %g",
            gsl_odeiv_step_name (step), desc, frac);

  gsl_odeiv_evolve_free (e);
}


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)
{
  int s = 0;
  double frac;

  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 (0.0, err_target, 1.0, 1.0);
  gsl_odeiv_evolve *e = gsl_odeiv_evolve_alloc (sys->dimension);

  while (t < t1)
    {
      gsl_odeiv_evolve_apply (e, c, step, sys, &t, t1, &h, y);
      /* printf ("SYS  t = %.18e h = %g y0 = %g y1= %g\n", t, h, y[0], y[1]); */
    }

  frac = fabs ((y[1] - yfin[1]) / yfin[1]) + fabs ((y[0] - yfin[0]) / yfin[0]);

  if (frac > 2.0 * e->count * err_target)
    {
      printf ("SYS  t = %.5e h = %g y0 = %g y1= %g\n", t, h, y[0], y[1]);
      s++;
    }

  gsl_test (s, "%s, %s, evolve, standard control, relative error = %g",
            gsl_odeiv_step_name (step), desc, frac);

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


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

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

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_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-03;
  p[1].type = gsl_odeiv_step_rk2imp;
  p[1].h = 1.0e-03;
  p[2].type = gsl_odeiv_step_rk4;
  p[2].h = 1.0e-03;
  p[3].type = gsl_odeiv_step_rk4imp;
  p[3].h = 1.0e-03;
  p[4].type = gsl_odeiv_step_rkf45;
  p[4].h = 1.0e-03;
  p[5].type = gsl_odeiv_step_rk8pd;
  p[5].h = 1.0e-03;
  p[6].type = gsl_odeiv_step_rkck;
  p[6].h = 1.0e-03;
  p[7].type = gsl_odeiv_step_bsimp;
  p[7].h = 0.1;
  p[8].type = gsl_odeiv_step_gear1;
  p[8].h = 1.0e-03;
  p[9].type = gsl_odeiv_step_gear2;
  p[9].h = 1.0e-03;
  p[10].type = 0;

  gsl_ieee_env_setup ();

  for (i = 0; p[i].type != 0; i++)
    {
      test_stepper_linear (p[i].type, p[i].h, GSL_SQRT_DBL_EPSILON);
      test_stepper_sin (p[i].type, p[i].h / 10.0, GSL_SQRT_DBL_EPSILON);
      test_stepper_exp (p[i].type, p[i].h / 10.0, GSL_SQRT_DBL_EPSILON);
      test_stepper_stiff (p[i].type, p[i].h / 10.0, GSL_SQRT_DBL_EPSILON);
    }

  for (i = 0; p[i].type != 0; i++)
    {
      test_evolve_exp (p[i].type, p[i].h, GSL_SQRT_DBL_EPSILON);
      test_evolve_sin (p[i].type, p[i].h, GSL_SQRT_DBL_EPSILON);
      test_evolve_stiff1 (p[i].type, p[i].h, 1e-5);
      test_evolve_stiff5 (p[i].type, p[i].h, 1e-5);
    }

  exit (gsl_test_summary ());
}