Another question on GSL ODE Solvers...

Another question on GSL ODE Solvers...

[Michael]

HI all,

If I want to vary the parameters "mu" with 10000 different values. Do
I really have to


for i = 1:1000

{

mu=Parameters(i);

gsl_odeiv_system sys = {func, jac, 2, &mu};

double t = 0.0, t1 = 100.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]);
         }

}

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

Is there a way to make things easier because I really only vary the
one of the parameters and nothing else changed at all, and I solve the
equation at the same time points t=5, 7, 9, 11.

Thank you!

MIchael.

Re: Another question on GSL ODE Solvers...

[Robert G. Brown]

On Wed, 6 Jun 2007, Michael wrote:

> HI all,
>
> If I want to vary the parameters "mu" with 10000 different values. Do
> I really have to
>
>
> for i = 1:1000
>
> {
>
> mu=Parameters(i);
>
> gsl_odeiv_system sys = {func, jac, 2, &mu};
>
> double t = 0.0, t1, tstep = 2.0;
>      double h = 1e-6;
>      double y[2] = { 1.0, 0.0 };
>

Just wrap up this segment in a loop.  Try the following:

      for(t1 = 5;t1 < 13;t1 += tstep){

>      while (t < t1) {
>
>          int status = gsl_odeiv_evolve_apply (e, c, s,
>                                               &sys,
>                                               &t, t1,
>                                               &h, y);
>
           if (status != GSL_SUCCESS){
             printf ("Oops.  Interval failed: %.5e %.5e %.5e %.5e\n", t, t1, y[0], y[1]);
             break;
           }
        }
        printf ("Results: y0(%.5e) = %.5e y1(%.5e) = %.5e\n", t, y[0], t, y[1]);
      }

What "should" happen is that you just get y0(5.00000e+00) = whatever,
y1(....) = .... for each of your requested time values.  The evolution
function may or may not be called a number of times to get from where
you start each step to where you finish, but you don't care.

Then try it for different tsteps and so on, try it for different
solvers as noted.

> Is there a way to make things easier because I really only vary the
> one of the parameters and nothing else changed at all, and I solve the
> equation at the same time points t=5, 7, 9, 11.

What's difficult?  You define the DE evaluator yourself (passing it the
one actual parameter any way you like, including in global memory or via
a macro if it is just a constant).  Then you call a simple loop to
output the result at exactly the points desired.  If you wanted to be
more general, you could create a vector t1[i] and fill it with desired
target times and then loop over i to cover the vector -- that's pretty
much what matlab does behind its pretty shell.  Then the "logic" would
be in filling t1[i].

There are also some amazing tricks one can play in C with pointers and
vectors for ODE systems with dimensionality and structure -- ways of
hiding whole Nth rank tensors of coupled ODEs behind the humble vector
and still writing the solver in terms of the tensor forms.  However I
sense that your problem is too simple for big guns like this.  Just
writing out a vector of parametrically varied solutions to a single ODE
is trivial -- write a vector of parameters (in mu or in global) and
initialize it, write your derivative routine to loop over that parameter
to fill in the vector of (pairwise?) INDEPENDENT derivatives and the
(diagonal?) jacobeans thereof, and use the very same loop.  The ODE
solver can handle a vector of derivatives -- you just have to tell it
the dimension and give it the vector of initial conditions.

    rgb

>
> Thank you!
>
> MIchael.
>

-- 
Robert G. Brown	                       http://www.phy.duke.edu/~rgb/
Duke University Dept. of Physics, Box 90305
Durham, N.C. 27708-0305
Phone: 1-919-660-2567  Fax: 919-660-2525     email:rgb@phy.duke.edu

Re: Another question on GSL ODE Solvers...

[Michael]

Hi RGB,

Thanks a lot for your answer.

Yes, it's not hard to write a loop. But it is hard to write such a
loop which is optimized for high performance.

Could you please show by example, which is the fastest way to organize
the loop to handle a vector of parameters(the coefficients of the ODE
systems)?

A rough estimate shows that I probably have to call the ODE45 solver
more than 10^9 times so the performance is critical.

Thanks!

Michael.

On 6/6/07, Robert G. Brown <rgb@phy.duke.edu> wrote:
> On Wed, 6 Jun 2007, Michael wrote:
>
> > HI all,
> >
> > If I want to vary the parameters "mu" with 10000 different values. Do
> > I really have to
> >
> >
> > for i = 1:1000
> >
> > {
> >
> > mu=Parameters(i);
> >
> > gsl_odeiv_system sys = {func, jac, 2, &mu};
> >
> > double t = 0.0, t1, tstep = 2.0;
> >      double h = 1e-6;
> >      double y[2] = { 1.0, 0.0 };
> >
>
> Just wrap up this segment in a loop.  Try the following:
>
>       for(t1 = 5;t1 < 13;t1 += tstep){
>
> >      while (t < t1) {
> >
> >          int status = gsl_odeiv_evolve_apply (e, c, s,
> >                                               &sys,
> >                                               &t, t1,
> >                                               &h, y);
> >
>            if (status != GSL_SUCCESS){
>              printf ("Oops.  Interval failed: %.5e %.5e %.5e %.5e\n", t, t1, y[0], y[1]);
>              break;
>            }
>         }
>         printf ("Results: y0(%.5e) = %.5e y1(%.5e) = %.5e\n", t, y[0], t, y[1]);
>       }
>
> What "should" happen is that you just get y0(5.00000e+00) = whatever,
> y1(....) = .... for each of your requested time values.  The evolution
> function may or may not be called a number of times to get from where
> you start each step to where you finish, but you don't care.
>
> Then try it for different tsteps and so on, try it for different
> solvers as noted.
>
> > Is there a way to make things easier because I really only vary the
> > one of the parameters and nothing else changed at all, and I solve the
> > equation at the same time points t=5, 7, 9, 11.
>
> What's difficult?  You define the DE evaluator yourself (passing it the
> one actual parameter any way you like, including in global memory or via
> a macro if it is just a constant).  Then you call a simple loop to
> output the result at exactly the points desired.  If you wanted to be
> more general, you could create a vector t1[i] and fill it with desired
> target times and then loop over i to cover the vector -- that's pretty
> much what matlab does behind its pretty shell.  Then the "logic" would
> be in filling t1[i].
>
> There are also some amazing tricks one can play in C with pointers and
> vectors for ODE systems with dimensionality and structure -- ways of
> hiding whole Nth rank tensors of coupled ODEs behind the humble vector
> and still writing the solver in terms of the tensor forms.  However I
> sense that your problem is too simple for big guns like this.  Just
> writing out a vector of parametrically varied solutions to a single ODE
> is trivial -- write a vector of parameters (in mu or in global) and
> initialize it, write your derivative routine to loop over that parameter
> to fill in the vector of (pairwise?) INDEPENDENT derivatives and the
> (diagonal?) jacobeans thereof, and use the very same loop.  The ODE
> solver can handle a vector of derivatives -- you just have to tell it
> the dimension and give it the vector of initial conditions.
>
>     rgb
>
> >
> > Thank you!
> >
> > MIchael.
> >
>
> --
> Robert G. Brown                        http://www.phy.duke.edu/~rgb/
> Duke University Dept. of Physics, Box 90305
> Durham, N.C. 27708-0305
> Phone: 1-919-660-2567  Fax: 919-660-2525     email:rgb@phy.duke.edu
>
>
>