[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

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

From: Jeroen Meidam
Subject: Re: [Help-gsl] ODE solver: oscillation around expected value
Date: Mon, 01 Aug 2011 13:50:35 +0100
User-agent: Mozilla/5.0 (Windows; U; Windows NT 6.1; en-GB; rv: Gecko/20110616 Lightning/1.0b2 Thunderbird/3.1.11

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);
    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:

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


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?

reply via email to

[Prev in Thread] Current Thread [Next in Thread]