[Octave-bug-tracker] [bug #47596] odepkg v0.8.5 : constant mass matrix option is not implemented correctly

[Guillaume Bruand]

URL:
  <http://savannah.gnu.org/bugs/?47596>

                 Summary: odepkg v0.8.5 : constant mass matrix option is not
implemented correctly
                 Project: GNU Octave
            Submitted by: gbruand
            Submitted on: ven. 01 avril 2016 12:09:36 GMT
                Category: Octave Forge Package
                Severity: 3 - Normal
                Priority: 5 - Normal
              Item Group: Incorrect Result
                  Status: None
             Assigned to: None
         Originator Name: 
        Originator Email: 
             Open/Closed: Open
         Discussion Lock: Any
                 Release: 4.0.0
        Operating System: Any

    _______________________________________________________

Details:

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






    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?47596>

_______________________________________________
  Message posté via/par Savannah
  http://savannah.gnu.org/