This is the mail archive of the mailing list for the GSL project.

Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

Questions on bsimp

Hi, there,

I just started using GSL for my research.  I was wondering why for some
ODEs, rkf45 can give me very good results while the results from bsimp is

My code is as follows, the solution should be y1 = cos(x) and y2 = sin(x).
When I use bsimp as the stepper function, it seems it just jumps too each
step.  I need bsimp because the speed is an important issue in my project.
The following code is just a test.

Thanks a lot!



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

int func (double t, const double y[], double f[], void *params)
 f[0] = -y[1];
 f[1] =  y[0];
 return GSL_SUCCESS;

int jac (double t, const double y[], double *dfdy, double dfdt[],
         void *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, 1.0);
 gsl_matrix_set (m, 1, 1, 0.0);
 dfdt[0] = 0.0;
 dfdt[1] = 0.0;
 return GSL_SUCCESS;

int main (void)
 const gsl_odeiv_step_type * T
  = gsl_odeiv_step_rkf45;

 gsl_odeiv_step * s
  = gsl_odeiv_step_alloc (T, 2);
 gsl_odeiv_control * c
  = gsl_odeiv_control_standard_new(1e-30, 1e-10,
                                                 1.0, 1.0);
 gsl_odeiv_evolve * e
  = gsl_odeiv_evolve_alloc (2);

 gsl_odeiv_system sys = {func, jac, 2, NULL};

 double t = 0.0, t1 = 10.0;
 double h = 1e-6;
 double y[2] = {1.0, 0.0};

 while (t < t1)
  int status = gsl_odeiv_evolve_apply(e, c, s,
          &t, t1,
          &h, y);
  if (status != GSL_SUCCESS)
  printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);

 gsl_odeiv_evolve_free (e);
 gsl_odeiv_control_free (c);
 return 0;

Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]