Re: Another question on GSL ODE Solvers...

[Michael <comtech.usa@gmail.com>]

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