This is the mail archive of the gsl-discuss@sourceware.org mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

Re: Another question on GSL ODE Solvers...


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



Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]