From: "Jochen Küpper" <jochen@fhi-berlin.mpg.de>
To: gsl-discuss@sourceware.org
Cc: Frank Reininghaus <frank78ac@googlemail.com>
Subject: Re: Bug report and fix for adaptive step size control in ODE solving
Date: Wed, 02 Jan 2008 14:18:00 -0000 [thread overview]
Message-ID: <0E55B479-624E-4F88-B511-CDFF93435D9E@fhi-berlin.mpg.de> (raw)
In-Reply-To: <4772D239.9070206@googlemail.com>
[-- Attachment #1: Type: text/plain, Size: 2107 bytes --]
Dear Frank, All,
On 26.12.2007, at 23:14, Frank Reininghaus wrote:
> I've found some problems with adaptive step size control which can
> lead to an infinite loop under special circumstances.
Nice write-up of this problem!
Although I have not looked through your cases in detail, I can
conform that we have had similar infinite-loop problems, for example,
on AMD Opteron systems running Linux and operated in 64bit mode
(gcc-4.x with -m64).
Because I did not have time to look into it in detail, and because
the q&d-hack below works for us, so far we use the following code
snippet to work around it:
while(t < tend) {
if(GSL_SUCCESS != gsl_odeiv_evolve_apply(evolve,
control, step,
&system, &t, tend, &h, y)) {
throw Error("Flight::operator(): integration
problem in GSL");
}
if(h < hmin) {
h = hinit;
}
}
with (for example)
hinit = tstep * 0.1
hmin = tstep * 0.01
We use this with gsl_odeiv_step_rkf45, gsl_odeiv_step_rkck, or
gsl_odeiv_step_rk8pd.
We would be quite interested in seeing this fixed in GSL, obviously.
I'll try to find some time to test your patch. However, why don't you
test what you really describe your comment? This would look like the
following (untested):
/* Step size h0 was decreased. Make sure it does not get too small:
It should be positive (i.e., h0 >= GSL_DBL_MIN) and large enough
to make t0 + h0 distinguishable from t0, i.e., at least of size
GSL_DBL_EPSILON*t0 */
h0 = GSL_MAX_DBL(GSL_DBL_MIN, (t0 != t0+h0) ? h0 :
GSL_DBL_EPSILON*t0);
This also avoids the multiplication in most cases.
Greetings,
Jochen
--
Einigkeit und Recht und Freiheit http://www.Jochen-
Kuepper.de
Liberté, Égalité, Fraternité GnuPG key: CC1B0B4D
Sex, drugs and rock-n-roll
[-- Attachment #2: This is a digitally signed message part --]
[-- Type: application/pgp-signature, Size: 186 bytes --]
next prev parent reply other threads:[~2008-01-02 14:18 UTC|newest]
Thread overview: 9+ messages / expand[flat|nested] mbox.gz Atom feed top
2007-12-26 22:14 Frank Reininghaus
2008-01-02 13:40 ` Kevin H. Hobbs
2008-01-02 14:18 ` Jochen Küpper [this message]
2008-01-02 18:51 ` Frank Reininghaus
2008-01-03 9:01 ` Unwanted compiler optimization (was: Bug report and fix for adaptive step size control in ODE solving) Jochen Küpper
2008-01-04 10:44 ` Brian Gough
2008-01-03 5:15 ` Bug report and fix for adaptive step size control in ODE solving Michael Stauffer
2008-01-03 15:54 ` Brian Gough
2008-01-04 9:46 ` Brian Gough
Reply instructions:
You may reply publicly to this message via plain-text email
using any one of the following methods:
* Save the following mbox file, import it into your mail client,
and reply-to-all from there: mbox
Avoid top-posting and favor interleaved quoting:
https://en.wikipedia.org/wiki/Posting_style#Interleaved_style
* Reply using the --to, --cc, and --in-reply-to
switches of git-send-email(1):
git send-email \
--in-reply-to=0E55B479-624E-4F88-B511-CDFF93435D9E@fhi-berlin.mpg.de \
--to=jochen@fhi-berlin.mpg.de \
--cc=frank78ac@googlemail.com \
--cc=gsl-discuss@sourceware.org \
/path/to/YOUR_REPLY
https://kernel.org/pub/software/scm/git/docs/git-send-email.html
* If your mail client supports setting the In-Reply-To header
via mailto: links, try the mailto: link
Be sure your reply has a Subject: header at the top and a blank line
before the message body.
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).