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

[Hg200]

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

@Markus: Thanks for pointing out the sparse case.

For the record: The IDA interface of ___ode15___() expects a cell with the two
matrices "{df/dy, df/dy'}" for the Jacobi matrix argument (see [1] p.76). I
was asking myself really a long time why the mass matrix is written to
"m_dfdyp" in the function

++
IDA& set_jacobian (SparseMatrix *dy, SparseMatrix *dyp,
                   DAEJacCellSparse j)
--

Given now is a differential equation system with Mass matrix "M", where
"M(t,y) * y' = f(t,y)". If one differentiates left and right sides, one gets
"df/dy' = M". I think this is the reason why for example in ode15s.m in

++
Function [jac, jact] = wrapjacfun (t, y, yp, Jac, Mass,
                                   havetimedep, havejacfun)
--

the Jacobi matrix and the Mass matrix are put together into the "options"
argument. Concerning our bug in ode15s.m line 274: If no Mass matrix is
specified in the ode15s() call, then we can argue that "I * y' = f(t, y)"
where "I" denotes the identity, and therefore we must put "{df/dy, I}" into
the "options" argument. Therefore it is a bug in the ode15s.m file and not in
___ode15___.c.

And in order to catch the sparse case, it is probably best to evaluate
"options.havejacsparse" and then use "speye (n)" instead of "eye (n)". 

If everyone agrees, then I will prepare the patch.

[1]: https://github.com/LLNL/sundials/tree/master/doc/ida

    _______________________________________________________

Reply to this item at:

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

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