help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Help-gsl] Re: Another question on GSL ODE Solvers...


From: Michael
Subject: [Help-gsl] Re: Another question on GSL ODE Solvers...
Date: Wed, 6 Jun 2007 12:31:03 -0700

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 <address@hidden> 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:address@hidden







reply via email to

[Prev in Thread] Current Thread [Next in Thread]