bugGNU Octave - Bugs: bug #38499, odepkg: solvers sometimes fail...

 
 

bug #38499: odepkg: solvers sometimes fail when error term gets too small

Submitter:  None
Submitted:  Mon 11 Mar 2013 05:28:01 AM UTC
   
 
Category:  Octave Package Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  cdf
Originator Name:  Matthew Chapman Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * 3.6.4
Operating System:  * Any Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Fri 23 Oct 2015 03:10:52 PM UTC, comment #3: 

The ode45 function has been moved into core Octave in the development branch (it will be in Octave 4.2) and has undergone a lot of modifications.

With the current version I see:


>> [t,u] = ode45(@(t,u) [0 1], [0 500], [0 0], odeset ());
warning: Option 'RelTol' not set, new value 0.001000 is used
warning: called from
    ode45 at line 181 column 5
warning: Option 'AbsTol' not set, new value 0.000001 is used
warning: called from
    ode45 at line 193 column 5
warning: option 'InitialStep' not set, estimated value 0.000100 is used
warning: called from
    ode45 at line 270 column 5
warning: Option 'MaxStep' not set, new value 50.000000 is used
warning: called from
    ode45 at line 281 column 5



so it seems this bug is not there anymore.
i am closing this ticket, please post a new bug if you see other problems with ode45.


Carlo de Falco <cdf>
Group Member
Fri 02 Oct 2015 04:14:06 PM UTC, comment #2: 

Hi,
Can someone comment about whether this issue still relevant with respect to the current version of ode45?
c.

Carlo de Falco <cdf>
Group Member
Mon 18 Mar 2013 09:52:16 PM UTC, comment #1: 

Your patch looks OK, but I haven't had time to test it yet.
Would it be possible for you to add also some embedded tests
to check that the solver behavior is correct after the patch
application?

Thanks,
c.

Carlo de Falco <cdf>
Group Member
Mon 11 Mar 2013 05:28:01 AM UTC, original submission:  

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

Anonymous

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #27596:  odepkg-min-vtau.patch added by None (3KiB - text/x-patch)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by cdf (Posted a comment)
  • -email is unavailable- added by cdf
  • -email is unavailable- added by None (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 6 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2015-10-23 cdf StatusNone Fixed
        Open/ClosedOpen Closed
    2013-03-18 cdf Assigned toNone cdf
        Releasedev 3.6.4
        Carbon-Copy- Added sebastian schöps <schoeps@gsc.tu-darmstadt.de>
    2013-03-11 None Attached File- Added odepkg-min-vtau.patch, #27596

    Back to the top

    Powered by Savane 3.13-d3ae.
    Corresponding source code