Hi all, I'm trying to solve the Van der Pol problem that is in the gsl documentation using the step type (aka solver) rk1imp - Euler backward - with fixed step size and low level functions (evolve, step, control). So far all I get is that the status from the evolve_apply_fixed_step is 3. Looking that the gsl code I found that this means that there is pointer problem. Now, if I change the solver to, say, rk2 the code runs. The strange thing is that if I use the driver with fixed step and rk1imp, the code runs. Any ideas of what may be wrong? Any help would be appreciated.
-- Juan PS. Here is my code using the low level functions: #include <stdio.h> #include <gsl/gsl_errno.h> #include <gsl/gsl_matrix.h> #include <gsl/gsl_odeiv2.h> int func (double t, const double y[], double f[], void *params) { double mu = *(double *)params; f[0] = y[1]; f[1] = -y[0] - mu*y[1]*(y[0]*y[0] - 1); return GSL_SUCCESS; } int jacobian (double t, const double y[], double *dfdy, double dfdt[], void *params) { double mu = *(double *)params; gsl_matrix_view dfdy_mat = gsl_matrix_view_array (dfdy, 2, 2); gsl_matrix * m = &dfdy_mat.matrix; gsl_matrix_set (m, 0, 0, 0.0); gsl_matrix_set (m, 0, 1, 1.0); gsl_matrix_set (m, 1, 0, -2.0*mu*y[0]*y[1] - 1.0); gsl_matrix_set (m, 1, 1, -mu*(y[0]*y[0] - 1.0)); dfdt[0] = 0.0; dfdt[1] = 0.0; return GSL_SUCCESS; } int main (void) { const double epsAbs = 1e-4, epsRel = 1e-4; double mu = 10.; double y[2] = { 1.0, 0.0 }; double t = 0.0, tf = 10.0, h = 1e-4, hstart = h; int status = 0; gsl_odeiv2_system system = { func, jacobian, 2, &mu }; const gsl_odeiv2_step_type * solvingAlgorithm = gsl_odeiv2_step_rk1imp; gsl_odeiv2_step * step = gsl_odeiv2_step_alloc(solvingAlgorithm, 2); gsl_odeiv2_control * control = gsl_odeiv2_control_y_new(epsAbs, epsRel); gsl_odeiv2_evolve * evolve = gsl_odeiv2_evolve_alloc(2); while (t < tf) { status = gsl_odeiv2_evolve_apply_fixed_step(evolve, control, step, &system, &t,hstart, y); if (status != GSL_SUCCESS) { printf ("error: evolve returned %d\n", status); break; } printf ("%.5e %.5e %.5e\n", t, y[0], y[1]); } gsl_odeiv2_evolve_free(evolve); gsl_odeiv2_control_free(control); gsl_odeiv2_step_free(step); return 0; } PSS And here my code using the driver: #include <stdio.h> #include <gsl/gsl_errno.h> #include <gsl/gsl_matrix.h> #include <gsl/gsl_odeiv2.h> int func (double t, const double y[], double f[], void *params) { double mu = *(double *)params; f[0] = y[1]; f[1] = -y[0] - mu*y[1]*(y[0]*y[0] - 1); return GSL_SUCCESS; } int jac (double t, const double y[], double *dfdy, double dfdt[], void *params) { double mu = *(double *)params; gsl_matrix_view dfdy_mat = gsl_matrix_view_array (dfdy, 2, 2); gsl_matrix * m = &dfdy_mat.matrix; gsl_matrix_set (m, 0, 0, 0.0); gsl_matrix_set (m, 0, 1, 1.0); gsl_matrix_set (m, 1, 0, -2.0*mu*y[0]*y[1] - 1.0); gsl_matrix_set (m, 1, 1, -mu*(y[0]*y[0] - 1.0)); dfdt[0] = 0.0; dfdt[1] = 0.0; return GSL_SUCCESS; } int main (void) { const double hstart = 1e-5, epsAbs = 1e-4, epsRel = 1e-4; double mu = 10.; const unsigned long numberOfSteps = 10000; const int problemDimension = 2; gsl_odeiv2_system sys = { func, jac, problemDimension, &mu }; gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk1imp, hstart, epsAbs, epsRel); double t = 0.0; double y[2] = { 1.0, 0.0 }; int i, s; for (i = 0; i < numberOfSteps; i++) { s = gsl_odeiv2_driver_apply_fixed_step (d, &t, hstart, numberOfSteps, y); if (s != GSL_SUCCESS) { printf ("error: driver returned %d\n", s); break; } printf ("%.5e %.5e %.5e\n", t, y[0], y[1]); } gsl_odeiv2_driver_free (d); return s; } _______________________________________________ Help-gsl mailing list Help-gsl@gnu.org https://lists.gnu.org/mailman/listinfo/help-gsl