[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Octave-bug-tracker] [bug #38499] odepkg: solvers sometimes fail when er
From: |
anonymous |
Subject: |
[Octave-bug-tracker] [bug #38499] odepkg: solvers sometimes fail when error term gets too small |
Date: |
Mon, 11 Mar 2013 05:28:03 +0000 |
User-agent: |
Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.22 (KHTML, like Gecko) Chrome/25.0.1364.160 Safari/537.22 |
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/
- [Octave-bug-tracker] [bug #38499] odepkg: solvers sometimes fail when error term gets too small,
anonymous <=