[Octave-bug-tracker] [bug #60240] Jacobian for ode15s as matrix fails

[Hg200]

Follow-up Comment #5, bug #60240 (project octave):

___ode15___.cc is an interface to IDA that provides more use cases than are
actually used in ode15s. ode15s is a subset and I think a general default of
df/dy' to the identity in ___ode15___.cc is wrong. It holds for "0 = y' - f(t,
y)", but not for the IDA scheme "0 = f(t, y, y')". For example, in ode15i(),
the Jacobi matrix has the form J(t, y, y'). Here you may have the case where
df/dy' in {df/dy, df/dy'} becomes unequal to the identity. If df/dy' is
specified empty (because someone forgets), this should be treated as an error,
because it could be a required non-identity matrix input.

>> Since this would mean a reduction of potential code duplication, ...

It looks to me that the two files ode15i.m and ode15s.m are written this way.
___ode15___.cc is treated as internal file (a 1:1 interface to IDA that does
not know the actual use case) and the ode15i.m and ode15s.m are wrapper
scripts that have the task to deliver correct and plausible input data.
Examples:

ode15s.m: If no mass matrix is specified and the Jacobi matrix is fun, the
identity is generated with speye () and passed in the else branch:
++
function [jac, jact] = wrapjacfun (t, y, yp, Jac, Mass,
                                   havetimedep, havejacfun)

  if (! isempty (Mass) && havetimedep)
    jact = Mass (t);
  elseif (! isempty (Mass))
    jact = Mass;
  else
    jact = speye (numel (y));
  endif

endfunction
--

In ode15i.m an error is thrown if the number of jacobians is not two:
++
options.havejac = true;
if (iscell (options.Jacobian))
  if (numel (options.Jacobian) == 2)
    J1 = options.Jacobian{1};
    J2 = options.Jacobian{2};
    if (   ! issquare (J1) || ! issquare (J2)
       || rows (J1) != n || rows (J2) != n
       || ! isnumeric (J1) || ! isnumeric (J2)
       || ! isreal (J1) || ! isreal (J2))
      error ("Octave:invalid-input-arg",
            'ode15i: "Jacobian" matrices must be real square matrices');
    endif
    if (issparse (J1) && issparse (J2))
       options.havejacsparse = true;  # Jac is sparse cell
    endif
  else
    error ("Octave:invalid-input-arg",
           'ode15i: invalid value assigned to field "Jacobian"');
endif
--

I have added a patch that should solve the problem in the .m domain so we can
see the changes for better understanding. The patch also includes a new BIST.

Maybe we can have more opinions on that topic.






(file #51096, file #51097)
    _______________________________________________________

Additional Item Attachment:

File name: patch_ode15s.diff              Size:1 KB
    <https://file.savannah.gnu.org/file/patch_ode15s.diff?file_id=51096>

File name: testcases.m                    Size:1 KB
    <https://file.savannah.gnu.org/file/testcases.m?file_id=51097>



    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?60240>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/