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

[John W. Eaton]

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