Re: Questions on bsimp

[Di Xiao <dxiao at physics dot utexas dot edu>]

Thanks a lot!  It works.

However, I still have a question.  We use adaptive stepzie control because 
we want to deduce the computational effort.  The algorithm should be smart 
to know at which region there is the smooth uninteresting curve, while in 
some region many small steps are needed.  My question is why bsimp chose 
some unreasonably big steps in my code.  It can jump from one period to 
another period without showing the details of each one.  If you  
run the code, you will find from the results givin by bsimp, it is hard 
to tell what function it likes. I tried different gsl_odeiv_control, but none 
of them gave me good results.  I guess the reason is when there are 
oscillations, bsimp is not a good one.  Can anybody explain it?

Thanks.

Di





On Thu, 1 Aug 2002, Slaven Peles wrote:

> On Thursday 01 August 2002 12:19 pm, you wrote:
> > 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.
> 
> When using adaptive step size different integrators will produce different 
> number of intermediate values during integration from t to t1. bsimp produces 
> fewer intermmediate values so that is why it seems to you like your results 
> are messy. Have look at the second example in the documentation at 
> http://sources.redhat.com/gsl/ref/gsl-ref_25.html#SEC381
> It's explained there how to "fix" this problem ;-).
> 
> Slaven
> 
> > --------------------------------------------------------------
> >
> > #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,
> >           &sys,
> >           &t, t1,
> >           &h, y);
> >   if (status != GSL_SUCCESS)
> >    break;
> >   printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
> >  }
> >
> >  gsl_odeiv_evolve_free (e);
> >  gsl_odeiv_control_free (c);
> >  gsl_odeiv_step_free(s);
> >  return 0;
> > }
>