bugGNU Octave - Bugs: bug #47596, [octave forge] (odepkg) constant...

 
 

bug #47596: [octave forge] (odepkg) constant mass matrix option is not implemented correctly

Submitter:  Guillaume Bruand <gbruand>
Submitted:  Fri 01 Apr 2016 12:09:36 PM UTC
   
 
Category:  Octave Package Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  In Progress Assigned to:  cdf
Originator Name:  Open/Closed:  * Open
Release:  * 4.0.0 Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Fri 25 Nov 2016 08:46:22 PM UTC, comment #5: 

Unfortunately this bug is also present in the ode45 function released with Octave 4.2.0.

I have attached a somewhat different patch against the 4.2.0 version of ode45.m.

(file #39074)

Bill Greene <billgreene>
Thu 07 Jul 2016 06:43:56 AM UTC, comment #4: 

Jacopo,

Was this change ever actually applied in odepkg/core?


c.

Carlo de Falco <cdf>
Group Member
Tue 05 Apr 2016 09:02:55 AM UTC, comment #3: 

Guillaume, hanks for checking!

@Jacopo could you please apply the suggested changes
to ode45 and ode23?

rather than computing the inverse, especially for a sparse
Jacobian, I would recommend computing and storing the factorization.


Carlo de Falco <cdf>
Group Member
Mon 04 Apr 2016 07:58:40 AM UTC, comment #2: 

The problem is also present in the development version; we must care to not forget the mass matrix when computing the solution. In ode45.m, line 262 :


 if (! isempty (odeopts.Mass) && isnumeric (odeopts.Mass))
    havemasshandle = false;
    mass = odeopts.Mass;  # constant mass
  elseif (isa (odeopts.Mass, "function_handle"))
    havemasshandle = true;    # mass defined by a function handle
  else  # no mass matrix - creating a diag-matrix of ones for mass
    havemasshandle = false;   # mass = diag (ones (length (init), 1), 0);
  endif


the mass matrix is defined in the constant case but not used. I suggest to use the parameter "fun" the same way a dynamic mass matrix is handled, e.g. line 301, add the bottom lines :


if (havemasshandle)   # Handle only the dynamic mass matrix,
    if (massdependence) # constant mass matrices have already
      mass = @(t,x) odeopts.Mass (t, x, odeopts.funarguments{:});
      fun = @(t,x) mass (t, x, odeopts.funarguments{:}) ...
             \ fun (t, x, odeopts.funarguments{:});
    else                 # if (massdependence == false)
      mass = @(t) odeopts.Mass (t, odeopts.funarguments{:});
      fun = @(t,x) mass (t, odeopts.funarguments{:}) ...
             \ fun (t, x, odeopts.funarguments{:});
    endif

# here the suggested fix

else
     invmass = inv(mass);
     fun = @(t,x) invmass * fun (t, x, odeopts.funarguments{:});
endif


Do not forget to uncomment the mass matrix line 268 (it will be equal to identity)

Guillaume Bruand <gbruand>
Sat 02 Apr 2016 05:47:03 AM UTC, comment #1: 

thanks for reporting this.

the development version of odepkg has undergone
a significant rewrite and many solvers have been
refactored to avoid code duplication.

furthermore some solvers are being moved to core
octave, and ode45 that you mention in your report
is one of them.

can you please check whether the bug also applies
to the development version of core octave?

you can find the relevant souces here:
http://hg.savannah.gnu.org/hgweb/octave/file/fa58fcb7c51d/scripts/ode/


Carlo de Falco <cdf>
Group Member
Fri 01 Apr 2016 12:09:36 PM UTC, original submission:  

Package : odepkg v0.8.5
Functions : any solver (ode45, ode54, etc) used in the package

Description :
Any use of the 'Mass' option with one of the solver (though the settings function 'odeset') will fail if the mass matrix is constant; the function runs without error/warning but the result is incorrect.

This problem does not occur when a function handle @mass is passed as the mass option.

Fix :
The change must be done in every solver. For example in ode54.m, line 232, mass definition must be uncommented :

vhavemasshandle = false; %# vmass = diag (ones (length (vinit), 1), 0);

to become :

vhavemasshandle = false;
vmass = diag (ones (length (vinit), 1), 0);


and then used line  367 :

vk(:,j) = feval ...
          (vfun, vthetime, vtheinput, vfunarguments{:});

must be changed into

vk(:,j) = vmass \ feval ...
          (vfun, vthetime, vtheinput, vfunarguments{:});


Important remark : in the case of a constant matrix we dont need to perform a matrix division at each solver step, because the computation time will greatly increase. It is better to compute the inverse of the mass matrix before the main loop, and then do a multiplication (rather than a matrix division). This gives :

l.227 :

if (~isempty (vodeoptions.Mass) && isnumeric (vodeoptions.Mass))
    vhavemasshandle = false;
    invmass = inv(vodeoptions.Mass); %# constant mass
  elseif (isa (vodeoptions.Mass, 'function_handle'))
    vhavemasshandle = true; %# mass defined by a function handle
  else %# no mass matrix - creating a diag-matrix of ones for mass
    vhavemasshandle = false;
    invmass = diag (ones (length (vinit), 1), 0);
  end


and line 356 :

if (vhavemasshandle)
...
...
...
        vk(:,j) = vmass \ feval ...
          (vfun, vthetime, vtheinput, vfunarguments{:});
      else
        vk(:,j) = invmass * feval ...
          (vfun, vthetime, vtheinput, vfunarguments{:});
      end


Guillaume Bruand <gbruand>

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #39074:  ode45.patch added by billgreene (1KiB - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by cdf
  • -email is unavailable- added by jwe (Updated the item)
  • -email is unavailable- added by billgreene (Updated the item)
  • -email is unavailable- added by cdf (Posted a comment)
  • -email is unavailable- added by gbruand (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 7 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2017-08-14 cdf Assigned toNone cdf
    2017-08-14 cdf Carbon-Copy- Added -email is unavailable-
    2017-08-13 jwe Summaryodepkg v0.8.5 : constant mass matrix option is not implemented correctly [octave forge] (odepkg) constant mass matrix option is not implemented correctly
    2016-11-25 billgreene Attached File- Added ode45.patch, #39074
    2016-04-05 cdf Carbon-Copy- Added -email is unavailable-
    2016-04-05 cdf StatusNeed Info In Progress
    2016-04-02 cdf StatusNone Need Info

    Back to the top

    Powered by Savane 3.13-3230.
    Corresponding source code