[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Bug-gsl] Numerical integration - Missing warning in doc, chap. 17.1 "In

From: Matthias Nagel
Subject: [Bug-gsl] Numerical integration - Missing warning in doc, chap. 17.1 "Introduction"
Date: Fri, 20 Jan 2012 19:41:11 +0100
User-agent: KMail/1.13.7 (Linux/2.6.39-gentoo-r3; KDE/4.6.5; x86_64; ; )


I followed the advice in the documentation, chapter 17.1 "Numerical Integration 
Introduction" that reads

"To compute to a specified relative error, set epsabs to zero."

Thereupon my code sometimes crashed with the error

"gsl: qag.c:148: ERROR: cannot reach tolerance because of roundoff error on 
first attempt".

Tracking down the error I       found that this always happens, if the true 
integral equals zero. This behavior is explainable, because the documentation 

| RESULT - I | <= ABSERR <= max( epsabs, epsrel |I| )

If "epsabs" is set to zero and the true integral "I" itself is zero, too, the 
admissible approximation error equals zero. In most cases this cannot be 
obtained and the function "gsl_integration_qag" aborts with the above error. If 
the true integral might become zero and one does not know that beforehand, the 
choice epsabs = 0 is not a good advice. I have two suggestions for improvement:

1) Put a warning into the documentation, that the advice "To compute to a 
specified relative error, set epsabs to zero" is dangerous if the true integral 
might be zero. At least this saves people from long lasting troubleshooting.

2) Do not abort the program execution, but return with the best possible 
approximation and GSL_EROUND. Actually that is, what I expected after reading 
the documentation. Does the documentation mention that "gsl_integration_qag" 
might abort the program execution, if the approximation bound cannot be 
reached? If yes, this should be more visible and and someone could add an 
explanation of the difference between GSL_EROUND and a program abortion.


reply via email to

[Prev in Thread] Current Thread [Next in Thread]