Re: Help with ODE solver

[Marco Antoniotti]

Thanks a lot for the answer.

Let me give you some background on this little project of mine.

I have a set of ODE in textual form which I parse to generate
Octave/Matlab code.  Generating the simple XDOT function to pass lsode
is easy.  It is how to set up experiments with "discontinuity" in
parameters that is unclear to me.

> From: "John W. Eaton" <address@hidden>
> Date: Thu, 12 Sep 2002 19:18:38 -0500
> Cc: address@hidden, address@hidden
> 
> On 12-Sep-2002, Marco Antoniotti <address@hidden> wrote:
> 
> | Pardon my obtuseness.  Do you mean that I will have to have a big 'if'
> | in the RHS function like
> | 
> |    if (t < t_crit_1)
> |       ...
> |    elseif (t_crit_1 <= t && t < t_crit_2)
> |       ...
> |    elseif (t > t_crit_k)
> |       ...
> |    end
> |    
> | and call
> | 
> |     result = lsdot('RHS', intial_conditions, timesteps, [t_crit_1; t_crit_2 
> ... t_crit_k]);
> | 
> | ?
> 
> 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.

> Perhaps someone else on the list can help you with the programming
> details.

I would appreciate this.  Thanks in advance.

> My reason for responding in this thread was to point out that using
> tcrit with lsode is important if you have a discontinuity in the
> derivative due to changing parameter values.  Simply breaking up the
> integration into intervals may not be sufficient because if you don't
> specify tcrit, lsode may integrate past your final time point and
> interpolate backward to find the final state values.  If you have a
> discontinuity there, then the integrator may still have trouble,
> respond by reducing the step size to zero, and eventually fail.

Ok, I think I at least understand that.

Cheers

-- 
Marco Antoniotti ========================================================
NYU Courant Bioinformatics Group        tel. +1 - 212 - 998 3488
715 Broadway 10th Floor                 fax  +1 - 212 - 995 4122
New York, NY 10003, USA                 http://bioinformatics.cat.nyu.edu
                    "Hello New York! We'll do what we can!"
                           Bill Murray in `Ghostbusters'.



-------------------------------------------------------------
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
-------------------------------------------------------------