Re: Help with ODE solver

[John W. Eaton]

On 13-Sep-2002, Marco Antoniotti <address@hidden> wrote:

| > Yes, that would be one way, but then you would need to make all those
| > t_crit_N values global.
| 
| I thing I cannot avoid it.  As it is, and because I need to be
| selective about what I plot, I alread make the initial_contitions
| global.
| 
| > It might make more sense to use a vector of
| > values and a loop in your RHS function.
| 
| Sorry, but I do not understand this.  Why would you want to do that?
| 
| > Or, you could vectorize it depending on how your time-dependent
| > parameter values are stored. 
| 
| Again I do not understand this.

OK, suppose you have a vector of N parameters and another vector of N
points that mark the beginnings of the intervals for which those
parameters are valid:

       p1     p2     ..     pN
    +------+------+------+------>
    t1     t2     ..     tN    


  p = [p1, p2, ..., pN];
  tcrit = [t1, t2, ..., tN];

then you can write

  tmp = find (t >= tcrit);
  i = tmp (length (tmp));

then p(i) is the parameter you want.

For example, use p = 10*(1:9), tcrit = 1:9, and t >= 0 and evaluate
the expressions

  t >= tcrit
  tmp = find (t >= tcrit)
  length (tmp)
  i = tmp (length (tmp))
  p(i)

to see how this will work.  Note that you could also use:

  i = sum (t >= tcrit)

instead of find and length, or if you don't actually need i, you can
write the entire expression as

  p (sum (t >= tcrit))

This is going to be much faster than a set of if/elseif conditions,
especially as the number of parameters becomes large.

The tcrit vector should be the same set of values that you pass to
lsode as the tcrit argument.  I would make it global and define it in
one place.  I would probably also make p global and define it
somewhere near the definition of tcrit so that their relationship is
obvious.

jwe



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------