This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
Re: Questions on bsimp
- From: Di Xiao <dxiao at physics dot utexas dot edu>
- To: Slaven Peles <peles at cns dot physics dot gatech dot edu>
- Cc: <gsl-discuss at sources dot redhat dot com>
- Date: Fri, 2 Aug 2002 10:53:03 -0500 (CDT)
- Subject: Re: Questions on bsimp
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;
> > }
>