Re: [Help-gsl] ODE solver: oscillation around expected value

[Jeroen Meidam]

In python I used the scipy.odeint solver.
I can give an example of what I am doing in the case of an adaptive stepper.
In the main I create the vector with initial conditions, calculated 
before then I calculate an integration time, I give it a sampling rate 
and from that I calculate the number of calculations that the solver needs.
  double y[4] = { r0, phi0, pr0, pphi0 } ;
  calculate_all(y, var, con) ;  /*calculate all variables at t=0.*/
  con->tfs =  25.0 * pow( ( con->M/(2.8*con->MO) ), -5.0/3.0) * 
pow(con->f0_GW/40.0), -8.0/3.0) ; /*Integration time in seconds*/
  con->tf = con->tfs/con->M ;                  /*Integration time*/
  con->SR = 8192.0 ;                             /*sample rate [s^-1]*/
  con->N = (int)(con->SR * con->tfs) ; /*number of steps*/
  con->delta = con->tf/con->N ;                /*stepsize [t/M]*/

Next I call this function which contains the solver:
void ode_solver_adaptivestep(double* y, struct variables* var, struct 
constants* con, FILE *fp) {

  void *params = var ;
  gsl_odeiv2_driver *d ;

  double err = con->accuracy ;
  double hstart = con->delta ; /* The initial stepsize */


  gsl_odeiv2_system sys = {func, jac, 4, params};

  switch (con->simulation_type) { /* To easily access different stepper 
functions*/
    case 5:
      printf ("Starting simulation:  Adaptive timestep, rk2...\n") ;
      d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk2, hstart,
                                                            err, err);
      break ;
    case 6:
      printf ("Starting simulation:  Adaptive timestep, rk4...\n") ;
      d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk4, hstart,

            /* And there are more options, but I want to save space */
  }

  int i,ti;
  double t = 0.0 ;
  double amplitude_h, omg ;
  int evaluations ;

  evaluations = con->N;

  for (i = 1; i <= evaluations; i++) {
    ti = i * con->tf / evaluations;
    calculate_all(y, var, con) ; /* This function calculates all 
variables needed for each step*/
    amplitude_h = pow(var->omega,1.0/3.0)*cos(2.0*y[1]) ;
    omg = var->omega ;
    int status = gsl_odeiv2_driver_apply (d, &t, ti, y);

    if (status != GSL_SUCCESS)
    {
      printf ("error, return value=%d\n", status);
      break;
    }
    if (y[0] < con->rLR) {
      printf ("Simulation has reached light ring at time %.3f and 
iteration %d.\n", t, i) ;
      break ;
    }
    fprintf(fp, "%.5f;%.5f;%.5f;%.5f;%.5f;%.5f;%.5f\n", t, amplitude_h, 
y[0], y[1], y[2], y[3], omg);
  }

  gsl_odeiv2_driver_free (d);

}

Does that give any insight into what might be a problem?


On 1-8-2011 13:32, Benjamin wrote:
> On 08/01/2011 02:08 PM, Jeroen Meidam wrote:
> > Hi,
> >
> > I'm using the gsl ode solver to solve a system of 4 coupled 
> > differential equations. I've done the exact same thing in python, but 
> > needed to rewrite it in c. There is one component of the plotted 
> > solutions that shows an unexpected oscillatory behaviour around the 
> > value that python gives as aresult.
> > I have no idea what could cause this. I have tried all the available 
> > stepping functions and none appear to get rid of this oscillation.
> >
> > Also, I would expect more differences between the different solving 
> > methods. If I plot for example the rk2 , rk4, rkf45 and rk8pd 
> > solutions in one graph, I can't see the slightest bit of difference. 
> > Is that to be expected?
> >
> >
> >
> > _______________________________________________
> > Help-gsl mailing list
> > address@hidden
> > https://lists.gnu.org/mailman/listinfo/help-gsl
> Hi,
>
> I am not an expert on this, but if you use the adaptive stepping 
> algorithm from odeiv the different solving methods could give very 
> similar results because of the adaptive control of errors.
>
> I would suspect that if you would plot their respective step sizes you 
> should see a difference.
>
> As for the oscillation I have no idea but I think it would in general 
> a good idea if you'd also post your source code/a working example so 
> that it will be easier to suggest something.
>
> Maybe you can also mention how you implemented the solution in python, 
> did you write your own solver routine or did you for example use the 
> scipy.odeint solver?