#include <base/quadrature_lib.h>
#include <cmath>
+#include <limits>
// please note: for a given dimension, we need the quadrature formulae
Quadrature<1> (n)
{
const unsigned int m = (n+1)/2;
- long double z;
- long double pp;
- long double p1, p2, p3;
- for (unsigned int i=1;i<=m;++i)
+ // tolerance for the Newton
+ // iteration below. we need to make
+ // it adaptive since on some
+ // machines (for example PowerPC)
+ // long double is the same as
+ // double -- in that case we can
+ // only get to a certain multiple
+ // of the accuracy of double there,
+ // while on other machines we'd
+ // like to go further down
+ const long double tolerance
+ = std::max (static_cast<long double>(std::numeric_limits<double>::epsilon() / 100),
+ static_cast<long double>(std::numeric_limits<long double>::epsilon() * 5));
+
+ for (unsigned int i=1; i<=m; ++i)
{
- z = std::cos(deal_II_numbers::PI * (i-.25)/(n+.5));
+ long double z = std::cos(deal_II_numbers::PI * (i-.25)/(n+.5));
- // Newton-iteration
+ long double pp;
+ long double p1, p2, p3;
+
+ // Newton iteration
do
{
// compute L_n (z)
pp = n*(z*p1-p2)/(z*z-1);
z = z-p1/pp;
}
- while (abs(p1/pp) > 1.e-19);
+ while (abs(p1/pp) > tolerance);
double x = .5*z;
this->quadrature_points[i-1] = Point<1>(.5-x);
<h3>base</h3>
<ol>
+ <li> <p> Fixed: The constructor of the <code class="class">QGauss</code> class
+ computed positions and weights of quadrature points in long double accuracy.
+ However, on machines where long double is the same as double, it therefore
+ never reached the requested accuracy, in effect leading to an infinite loop.
+ This is now fixed.
+ <br>
+ (WB 2002/09/19)
+ </p>
+
<li> <p> New: The <code class="class">Function</code> class now
exports the value of its template argument through the static
member variable <code class="member">dimension</code>.