[Octave-bug-tracker] [bug #38499] odepkg: solvers sometimes fail when error term gets too small

[anonymous]

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

                 Summary: odepkg: solvers sometimes fail when error term gets
too small
                 Project: GNU Octave
            Submitted by: None
            Submitted on: Mon 11 Mar 2013 05:28:01 AM UTC
                Category: Octave Forge Package
                Severity: 3 - Normal
                Priority: 5 - Normal
              Item Group: Incorrect Result
                  Status: None
             Assigned to: None
         Originator Name: Matthew Chapman
        Originator Email: address@hidden
             Open/Closed: Open
         Discussion Lock: Any
                 Release: dev
        Operating System: Any

    _______________________________________________________

Details:

It appears that when using a variable step size, if the estimated error
'vdelta' gets too small (e.g. if your solution is converging to a constant),
the step size sometimes starts getting smaller rather than larger, until the
solution fails.  Here is a trivial example that demonstrates this behaviour in
the ode45 function...

Example
-------


octave:1> [t,u] = ode45(@(t,u) [0 1], [0 500], [0 0]);
warning: Option "RelTol" not set, new value 0.000001 is used
warning: Option "AbsTol" not set, new value 0.000001 is used
warning: Option "InitialStep" not set, new value 50.000000 is used
warning: Option "MaxStep" not set, new value 50.000000 is used
error: Solving has not been successful. The iterative integration loop exited
at time t = 267.144937 before endpoint at tend = 500.000000 was reached. This
may happen if the stepsize grows smaller than defined in vminstepsize. Try to
reduce the value of "InitialStep" and/or "MaxStep" with the command "odeset".



Reason
------

It seems the code that causes this is around line 440 of ode45.m:


      %# It could happen that max (vdelta) == 0 (ie. that the original
      %# vdelta was 0), in that case we double the previous vstepsize
-->   vdelta(find (vdelta == 0)) = max (vtau) .* (0.4 ^ (1 / vpow));

      if (vdirection == 1)
        vstepsize = min (vodeoptions.MaxStep, ...
           min (0.8 * vstepsize * (vtau ./ vdelta) .^ vpow));


Before this code vdelta = [0;0] and vtau = [1.0000e-06;2.6714e-04] (since
u=[0;267.14]).

After replacing the zero values with max(vtau).*(0.4^(1/vpow)), vdelta =
[2.7356e-06;2.7356e-06].

Now the value of 0.8*(vtau./vdelta).^vpow) is [0.65415;1.99999], since we take
the min() of that we multiply vstepsize by a factor of 0.65415 - this is
different from.the comment which says we should double vstepsize in this case
(which makes much more sense).


Suggested solution
------------------
I would suggest the simplest solution is to amend line 440 to replace the zero
values with min(vtau).*(0.4^(1/vpow)) rather than max(vtau).*(0.4^(1/vpow)). 
Then, since subsequent code also uses min() [or max() of the negative value],
the code will always have the expected effect of doubling the vstepsize.  I am
attaching a patch to this effect for all the affected functions.


Thanks,
Matt




    _______________________________________________________

File Attachments:


-------------------------------------------------------
Date: Mon 11 Mar 2013 05:28:01 AM UTC  Name: odepkg-min-vtau.patch  Size: 3kB 
 By: None

<http://savannah.gnu.org/bugs/download.php?file_id=27596>

    _______________________________________________________

Reply to this item at:

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

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/