Setting multiple ODE parameters

Setting multiple ODE parameters

[Elie Allouis]

Dear All,
 
I know this topic has already been covered to some extent here, but I unfortunately, it seems that the answers haven't been posted to the list and therefore, I could not find a clear answer to my queries:
 
I am currently implementing an ODE problem with around 6 equations and two structures. One with 8 constants parameters, and the other with some 13 variable parameters calculated elsewhere at each step.
 
I did implement an ODE, before, but I have found GSL much more elegant, but I am puzzled with the void *param in the function and jacobian definition. I would need some clues/advice/pointers into how to deal with these parameters into the ODE function. 
 
If you had any examples other than the Van der pol equation that would be great and much appreciated.
 
Thanking you in advance,
 
Kind regards,
 
Elie Allouis
 
 
------------------------
Elie Allouis
U.5
Surrey Space Centre
Surrey University
Guildford, UK
-------------------------
 

Re: Setting multiple ODE parameters

[Kevin H. Hobbs]

[-- Attachment #1: Type: text/plain, Size: 516 bytes --]

Versions of my projects are at http://crab-lab.zool.ohiou.edu/kevin/ I
pass all kinds of things to models this way.

I usually fill a structure with things that the model will need like
pointers to data, places where it can store results, and the parameters
of the run.  Then I assign a pointer to this structure to params. 
Inside the model I follow this pointer to get to the data.

If all you want to do is pass an array of parameters, you just use
params like you normally pass an array to a function.

[-- Attachment #2: This is a digitally signed message part --]
[-- Type: application/pgp-signature, Size: 189 bytes --]

Re: Setting multiple ODE parameters

[Robert G. Brown]

On Mon, 16 Aug 2004, Elie Allouis wrote:

> Dear All,
>
> I know this topic has already been covered to some extent here, but I
> unfortunately, it seems that the answers haven't been posted to the list
> and therefore, I could not find a clear answer to my queries:
>
> I am currently implementing an ODE problem with around 6 equations and
> two structures. One with 8 constants parameters, and the other with some
> 13 variable parameters calculated elsewhere at each step.
>
> I did implement an ODE, before, but I have found GSL much more
> elegant, but I am puzzled with the void *param in the function and
> jacobian definition. I would need some clues/advice/pointers into how to
> deal with these parameters into the ODE function.
>
> If you had any examples other than the Van der pol equation that would
> be great and much appreciated.

There are two general ways one can pass arguments to a subroutine in C.
One is via its call (inside the parantheses of its call instance in a
routine).  The other is via global or shared variables.

The contents of the data items inside function call parentheses are typically
pushed onto the stack, from which they can be retrieved inside the function
via a suitable displacement on the stack pointer (all magically linked up for
you by the linker -- one doesn't generally do this arithmetic by hand any
more although the *params method has been used plenty with hand-set
offsets in the old days).

Values passed in this way go away when the stack pointer is popped post call,
so you do not in general change the original values by changing them inside
the subroutine -- the passed variables are all "local".

An exception is when you pass to the subroutine the ADDRESS of a
variable, and load it into a POINTER to a variable inside the
subroutine.  In that case when you modify the contents of the variable
pointed to by the pointer, you modify the original variable (since there is
only one instance of the data and both functions share its address).

When writing a general purpose (e.g. ODE solver) subroutine, it is of course
impossible to know what or how many variables are used as arguments in the
function that evaluates the deriviative vector.  So instead one passes in a
generic argument that is a POINTER to the full argument list.  Ideally this
list would be packed into some sort of struct so that as an argument you just
pass the address of the struct, and inside the function you'd cast the void
pointer to the struct type so you could dereference it in a meaningful way.

Hopefully that is enough of an explanation so that the van der pol oscillator
example is clear -- double mu is the parameter, so its ADDRESS is passed to
the ODE solver in the setup line.  Inside func or jac, the contents (first *)
of the double address this points to (double * cast) are put into a local
double variable also named mu.  This KEEPS the value of mu in the func local.

A more interesting/useful (and not quite so obscure) example would be:

 typedef struct {

  double one;
  double two;

 } Mu;

 Mu mu;

 mu.one = 0.1;
 mu.two = 0.2;

Now one STILL passes the address of the mu struct into the ode solver setup
via:

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



Then in the func and jac routines you can do a cast such as (I think,
I'm not testing this as I go:-)

 Mu *mu = (Mu *)params;

which reads like "create local pointer to a Mu struct named mu, and put into 
it the contents of params (which should be the address of the mu variable
in the calling routine)".  You can then dereference the two componenents as

 mu->one;
 mu->two;

in the func and jac routines.  Obviously you can make the params struct as
complex as you like as long as its prototype is shared among the various
routines that use it, and by putting all the ode control variables into a
struct (which isn't necessarily the only way to proceed) you keep your code
tidy.  Note that if you set mu->one = 0.2 in func, this change DOES appear in
the calling routine's copy of mu, as again there is a single memory location
for the mu struct in the system, but both the calling routine and func now
have its address.

The alternative way to proceed requires no example -- just put your func and
jac arguments, variable or constant, into global (or at least shared)
variables.  This can easily be done with a suitable include file.  Set the
variables in the calling routine, use the variables in the func and jac
routines, nothing to it.  In this case the *params argument can be (should be,
I'd think) null, as you aren't passing arguments at all.  Note well that in
this case changes made to these variables in func or jac WILL affect the
values back in the calling routine, as there are only single copies of the
actual data shared across all the routines.

Which of these is "better" depends on your sense of style.  Many people eschew
global/shared variables to get data into or out of subroutines as "inelegant"
as they make the code less portable.  However, I personally think that they
are often both elegant and economical -- for example when portability isn't an
issue and you never plan to reuse your func and jac routines outside of the
specific application you are working on.  In this case, global variables are
simple and keep you from having to create several variables of the same name
(note that mu is defined in both the calling routine and in func and jac, all
with the same name). If your argument list is very long, putting the data in
global variables prevents the system from having to duplicate it all on every
call.  Although for small lists the difference is small, global memory is
invariably more efficient than passing arguments (which always have to be
duplicated and dereferenced off of the stack).  Efficiency is a valid
component of elegance.

Object oriented people would tend to go with passed arguments, procedural
programmers (especially old ones who have actually programmed in assembler and
are visualizing just what is going in when subroutines are called from the
description above) are a bit more flexible and would go either way as the
issue of portability and reuse and overall code structure vs efficiency
dictate.  I think, anyway -- this isn't an attempt to stimulate a religious
war but rather to describe a couple of solutions to your problem.

Hope this helps.

   rgb

>  
> Thanking you in advance,
>  
> Kind regards,
>  
> Elie Allouis
>  
>  
> ------------------------
> Elie Allouis
> U.5
> Surrey Space Centre
> Surrey University
> Guildford, UK
> -------------------------
>  
> 

-- 
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: Setting multiple ODE parameters

[Robert G. Brown]

On Tue, 17 Aug 2004, Elie Allouis wrote:

> I appreciate all the answers I have received. I want to thanks
> particularly Mr Roger Brown for his rather complete analysis of my
> query. I think it has clarified the topic.

Umm, that would be >>Robert<< Brown (or just plain rgb), and you're very
welcome...;-)

   rgb

-- 
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: Setting multiple ODE parameters

[Elie Allouis]

I appreciate all the answers I have received. I want to thanks particularly Mr Roger Brown for his rather complete analysis of my query. I think it has clarified the topic.
 
I am posting below one of the answer I received that has not been posted on the list, in case it can be useful to somebody in the future. This will provides an other practical example of Mr Brown analysis.
 
Thank you very much,
 
Rgds,
 
Elie Allouis
 
--------------------------------------------
Elie Allouis
U.5
Surrey Space Centre
Guildford, UK
--------------------------------------------

 
<<<- Edited--->>>
 
I define a structure as:
struct ODE_params
{
  double Q, drive_freq, r;
};

In the code I define a variable of type structure like 
  struct ODE_params rod_param={ 1.2577,    // Quality factor "Q"
1.17731,   // Drive frequency "a"
1.088 }; // Radii ratio "r"

I pass the pointer to rod_param to the ODE integrator as
  gsl_odeiv_system sys = {func, jac, ODE_DIM, &rod_param};

Functions 'func' and 'jac' may use some of the parameters. At the beginning of 
each function I add these several lines to recast void pointer into pointer 
to the structure and assign elements of the structure to local variables 
inside the function.

int
func (double t, const double y[], double f[],
      void *params)
{
  struct ODE_params *p 
    = (struct ODE_params *) params;

  double Q = p->Q;
  double a = p->drive_freq;
  double r = p->r;

(...)

Hope this helps. Also, have a look in gsl documentation under multiroot 
solvers. In the examples there they use the same trick. 
 
<<<----End--->>>>