Re: [Help-gsl] [ODE] GSL and Octave comparison

[Marco Maggi]

"Brian Gough" wrote:
>Marco Maggi wrote:
>> and then recomputing with GSL step function rk2
>> control function 'y' eps_abs = eps_rel = 0.001,
>
> Those are large tolerances -- they are applied
> to each step, so the global error is
> O(Nsteps * tolerance).  Also rk2 is the least
> accurate method, try one of the others,
> e.g. rk4 at least.

As the following test demonstrates this is not
the problem. For such a simple evolution even
rk2 with error control parameters set to 0.1
yields recognisable-as-correct results.

There is some error in my program that I will
have to discover.

/* odetest.c */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv.h>


const double    E   = 100.0;
const double    tau = 0.2;
const double    y0  = 2.0;

/*
------------------------------------------------------------ */

double
function (double time, double y)
{
  return (-y/tau + E/tau);
}
double
jacobian (double time, double y)
{
  return (-1.0/tau);
}
double
time_derivative (double time, double y)
{
  return 0.0;
}
double
solution (double time)
{
  return ((y0-E) * exp(-time/tau) + E);
}

/*
------------------------------------------------------------ */

int
ode_function (double t, const double y[], double dy_dt[],
void * params)
{
  dy_dt[0] = function(t, y[0]);
  return GSL_SUCCESS;
}
int
ode_derivatives (double t, const double y[],
                 double * df_dy, double df_dt[], void * params)
{
  df_dy[0] = jacobian(t, y[0]);
  df_dt[0] = time_derivative(t, y[0]);
  return GSL_SUCCESS;
}

/*
------------------------------------------------------------ */

int
main (int argc, const char *const argv[])
{
  gsl_odeiv_step *      step_function;
  gsl_odeiv_control *   control_function;
  gsl_odeiv_evolve *    evolve;
  gsl_odeiv_system      system;
  size_t                dimension = 1;
  double                eps_abs, eps_rel, a_y, a_dydt;
  double                current_time_instant, current_time_step = 0.01;
  double                final_time_instant, current_value = y0;
  int                   e;


  eps_abs = eps_rel = a_y = a_dydt = 0.1;

  step_function =
    gsl_odeiv_step_alloc(gsl_odeiv_step_rk2, dimension);
  control_function =
    gsl_odeiv_control_standard_new(eps_abs, eps_rel, a_y,
a_dydt);
  evolve = gsl_odeiv_evolve_alloc(dimension);

  system.function       = ode_function;
  system.jacobian       = ode_derivatives;
  system.dimension      = dimension;
  system.params         = NULL;

  current_time_instant = 0.0;
  fprintf(stdout, "y(%f) = %f\n", current_time_instant,
current_value);
  for (final_time_instant=0.2;
       final_time_instant <= 1.0;
       final_time_instant += 0.2)
    {
      while (current_time_instant < final_time_instant)
        {
          e = gsl_odeiv_evolve_apply(evolve, control_function,
                                     step_function, &system,
                                     &current_time_instant, final_time_instant,
                                     &current_time_step, &current_value);
          if (e)
            {
              fprintf(stderr, "%s: %s\n", argv[0], gsl_strerror(e));
              exit(EXIT_FAILURE);
            }
        }
      fprintf(stdout, "y(%f) = %f\t\tsol(%f) = %f\n",
              current_time_instant, current_value,
              current_time_instant, solution(current_time_instant));
    }

  gsl_odeiv_evolve_free (evolve);
  gsl_odeiv_control_free(control_function);
  gsl_odeiv_step_free   (step_function);
  exit(EXIT_SUCCESS);
}

/* end of file */


--
Marco Maggi

"They say jump!, you say how high?"
Rage Against the Machine - "Bullet in the Head"