Re: ode45 gives solution in the "Dark Side," but not in Octave

[John B. Thoo]

On Sep 2, 2009, at 5:30 PM, John W. Eaton wrote:

> On  2-Sep-2009, John B. Thoo wrote:
>
> | On Sep 1, 2009, at 3:00 PM, John W. Eaton wrote:
> |
> | > On  1-Sep-2009, John B. Thoo wrote:
> | >
> | > | Hi.  I have a question about using ode45 in Octave 3.2.2 (and  
> OdePkg
> | > | on Mac OS X 10.4.11) vs. in the "Dark Side" 7.4.0 (R2007a) (on a
> | > | server ... Linux?).  My code can be found here at the moment:
> | > |
> | > | <http://ms.yccd.edu/~jb2/PickUp/twoway.zip>
> | > |
> | > | "run_twoway.m" executes the code.
> | > |
> | > | Both Octave and the Dark Side solve the problem (i.e., do not  
> return
> | > | an error) to produce a solution  u,  and both give
> | > |
> | > |  > size (u)
> | > | ans =
> | > |
> | > |     8192    102
> | > |
> | > | However, while I can plot the solution in the Dark Side, in  
> Octave I
> | > | can plot only the first four frames.  In fact, in Octave I get
> | > NaN +
> | > | 0.00000i  from  u(:,5)  on.  Why is that?
> | > |
> | > | I would appreciate any hints.  Thanks.
> | >
> | > Have you tried using lsode in Octave instead of ode45 (which is  
> not
> | > part of Octave, BTW).
> |
> | I haven't tried lsode.  The reason for my using ode45 is my advisor
> | uses the Dark Side only, and we need to run the same code.  Thanks.
>
> Here's a simple wrapper for lsode that I think will provide a simple
> ode45 interface.  It doesn't cover all the options, but it will
> probably work in a lot of cases, provided that you just need a
> solution, and not specifically a Runge-Kutta method.  It shoudl be
> fairly easy to fix the interface if ti isn't compatible enough for
> you.
>
> jwe
>
> function [t_out, x] = ode45 (f, t, x0, options)
>
>   if (nargin == 3 || nargin == 4)
>
>     global __ode45s_rhs_fcn__;
>     __ode45s_rhs_fcn__ = f;
>
>     if (nargin == 4)
>       if (isfield (options, "RelTol") && ! isempty (options.RelTol))
>         lsode_options ("relative tolerance", options.RelTol)
>       endif
>       if (isfield (options, "AbsTol") && ! isempty (options.AbsTol))
>         lsode_options ("absolute tolerance", options.AbsTol)
>       endif
>     endif
>
>     x = lsode (@__lsode_rhs__, x0, t);
>     t_out = t;
>
>   else
>     print_usage ();
>   endif
>
> endfunction
>
> %  
> ----------------------------------------------------------------------
>
> function xdot = __lsode_rhs__ (x, t)
>
>   global __ode45s_rhs_fcn__;
>
>   xdot = __ode45s_rhs_fcn__ (t, x);
>
> endfunction

Hi, jwe.  Thanks for taking time on this.

I'm sorry, but I'm not familiar with using a "wrapper," so I haven't  
tried it, yet.  Do I save the above to a .m file?  Then how do I call  
on it?

In any case, I tried lsode directly in my code, in which, making a  
guess for how to proceed, I replaced the lines in "run_twoway.m"

options = odeset ();
[t, w] = ode45 (@F, tspan, w0, options, k, lx, ltrunc, epsilon, mu);

with the line

[w, istate, msg] = lsode (@F, tspan, w0, k, lx, ltrunc, epsilon, mu);

But when I executed the code, I got this:

octave-3.2.2:5> run_twoway
error: Invalid call to lsode.  Correct usage is:

 -- Loadable Function: [X, ISTATE, MSG] = lsode (FCN, X_0, T, T_CRIT)


Additional help for built-in functions and operators is
available in the on-line version of the manual.  Use the command
`doc <topic>' to search the manual index.

Help and information about Octave is also available on the WWW
at http://www.octave.org and via the address@hidden
mailing list.
octave-3.2.2:5>

Thanks again.

---John.