"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, ¤t_time_instant, final_time_instant, ¤t_time_step, ¤t_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" _______________________________________________ Help-gsl mailing list Help-gsl@gnu.org http://lists.gnu.org/mailman/listinfo/help-gsl