Questions on bsimp

Questions on bsimp

[Di Xiao]

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
messy.

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!

Di

--------------------------------------------------------------

#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;
}

Re: Questions on bsimp

[Brian Gough]

Di Xiao writes:
 >  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?

Bsimp can take very large steps and still be highly accurate.  See the
reference to the paper in the manual for details.

Re: Questions on bsimp

[Di Xiao]

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;
> > }
> 

Re: Questions on bsimp

[Slaven Peles]

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;
> }

-------------------------------------------------------