Re: Questions on bsimp

[Slaven Peles <peles@cns.physics.gatech.edu>]

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

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