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

[Guillaume Bruand]

Follow-up Comment #2, bug #47596 (project octave):

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)

    _______________________________________________________

Reply to this item at:

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

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