bugGNU Octave - Bugs: bug #60240, Jacobian for ode15s as matrix fails

 
 

bug #60240: Jacobian for ode15s as matrix fails

Submitter:  Andreas Stahel <andreasstahel>
Submitted:  Tue 16 Mar 2021 10:04:01 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Unexpected Error or Warning
Status:  Fixed Assigned to:  None
Originator Name:  Andreas Stahel Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Wed 31 Mar 2021 03:30:32 PM UTC, comment #7: 

The test cases provided in this report seem to be working and there were no complaints about any regressions in the last 10 days.

Closing report.

Markus Mützel <mmuetzel>
Group administrator
Sun 21 Mar 2021 02:02:42 PM UTC, comment #6: 

Thanks for clearing that up for me.
You've convinced me that you've been right from the start. :-)

I pushed your patch here:
https://hg.savannah.gnu.org/hgweb/octave/rev/ebc3f80673f0

Marking as ready for test.

Markus Mützel <mmuetzel>
Group administrator
Sun 21 Mar 2021 01:01:01 PM UTC, comment #5: 

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

Hg200 <hg200>
Sun 21 Mar 2021 11:02:01 AM UTC, comment #4: 

I think I can follow most of your reasoning. Apart from the last point where you conclude that `ode15s.m` needs to be fixed.

Imho, it is an arbitrary requirement that the mass matrix must be set if the Jacobian is a matrix. If the mass matrix is omitted, the user probably expects that the system of differential equations `y' = f(t,y)` is solved (instead of `M(t,y) * y' = f(t,y)`).
Like you wrote, that means that the mass matrix could be set to the identity matrix to solve that system.
In principle, it doesn't matter if we default to the identity matrix in `ode15s.m` or in `__ode15__.cc`. However, if we decided to do the change in `ode15s.m` and left `__ode15__.cc` as it is, that would mean that every function that calls `__ode15__` would need to be checked if it sets the correct default value for the mass matrix.
If instead we decided to set the mass matrix to the identity matrix (if it is empty) in `__ode15__.cc`, we (and anyone else using `__ode15__`) wouldn't need to worry about implementing that check. Since this would mean a reduction of potential code duplication, I'd argue to change `__ode15__.cc`.

Markus Mützel <mmuetzel>
Group administrator
Sat 20 Mar 2021 10:44:46 PM UTC, comment #3: 

@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

Hg200 <hg200>
Sat 20 Mar 2021 01:22:06 PM UTC, comment #2: 

@hg200: Both solutions you propose sound fine to me.

Because `__ode15__` might be called by other functions as well (and IIRC it is used by at least one Octave Forge package), I'm leaning towards defining a reasonable default value for the mass matrix in that function (maybe using `octave::identity_matrix`).

Maybe the `set_jacobian` function would be a good place to do that?

Could the same occur for sparse matrices?

Could you please prepare a changeset?

Markus Mützel <mmuetzel>
Group administrator
Thu 18 Mar 2021 07:37:39 PM UTC, comment #1: 

The mentioned code in the OP runs with Matlab, but fails in Octave devel and also fails with Octave stable. If we stop in ode15s.m, we see that for the code in the OP, the mass matrix in the option field is empty. Therefore an empty mass matrix is copied to m_dfdyp in __ode15__.cc in


IDA&
set_jacobian (matrix *dy, matrix *dyp, DAEJacCellDense j)


In __ode15__.cc:688 the following statement then crashes:


if (IDASolve (m_mem, tend, &tsol, yy, yyp, IDA_ONE_STEP) != 0)


To test this, the attached patch replaces the empty mass matrix with the identity. With this change, ode15s() returns the correct result if i run the code from the OP.

The question now would be how to initialize the IDA interface when no mass matrix is specified. But maybe passing the identity in ode15s.m is also a good fix? At least sound reasonable.


(file #51087)

Hg200 <hg200>
Tue 16 Mar 2021 10:04:01 AM UTC, original submission:  

When calling ode15s() with a Jacobian passed as function handle it works as expected.
When the Jacobian is passed as matrix, it fails.


%% works well if Jacobian passed as function handle
opt = odeset ("RelTol", 1e-8, "AbsTol", 1e-6,'Jacobian',@(t,y)[98, 198;-99, -199],'Stats','on');
%% does not return results of Jacobian passed as matrix
%opt = odeset ("RelTol", 1e-8, "AbsTol", 1e-6,'Jacobian',[98, 198;-99, -199],'Stats','on');
f = @(t,y)[98, 198;-99, -199]*(y-[1;0]);
y_exact = @(t)2*exp(-t)-exp(-100*t)+1;

disp('results with Jacobian')
[t15sJ, y15sJ] = ode15s (f ,[0, 5],[2;0], opt);  y15sJ = y15sJ(:,1);

Andreas Stahel <andreasstahel>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #51096:  patch_ode15s.diff added by hg200 (2KiB - text/x-patch)
file #51097:  testcases.m added by hg200 (1KiB - text/x-objcsrc)
file #51087:  identity_for_mass.patch added by hg200 (548B - text/x-patch)
file #51072:  ODE15s_Question.m added by andreasstahel (981B - application/vnd.wolfram.mathematica.package)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by hg200 (Updated the item)
  • -email is unavailable- added by andreasstahel (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 9 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2021-03-31 mmuetzel StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2021-03-21 mmuetzel StatusConfirmed Ready For Test
    2021-03-21 hg200 Attached File- Added patch_ode15s.diff, #51096
        Attached File- Added testcases.m, #51097
    2021-03-20 mmuetzel StatusNone Confirmed
        Operating SystemGNU/Linux Any
    2021-03-18 hg200 Attached File- Added identity_for_mass.patch, #51087
    2021-03-16 andreasstahel Attached File- Added ODE15s_Question.m, #51072

    Back to the top

    Powered by Savane 3.13-d3ae.
    Corresponding source code