Re: Another question on GSL ODE Solvers...

["Robert G. Brown" <rgb@phy.duke.edu>]

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