#174 Error in simplex quadrature rules of order 6 and 7

Metadata

Property	Value
Reported by	unknown (unknown)
Reported at	Sep 6, 2006 12:16
Type	Bug Report
Version	Git (pre2.4) [autotools]
Operating System	Unspecified / All
Closed by	Oliver Sander (oliver.sander@tu-dresden.de)
Closed at	Oct 4, 2006 12:41
Closed in version	Unknown
Resolution	Fixed
Comment

Description

This is NOT a duplicate of #162 (closed)! The triangle quadrature rules of order 6 and 7 produce wrong results. This is easily checked by integrating y^p:

#include

#include "dune/config.h" #include "dune/grid/common/quadraturerules.hh"

int main() { using namespace Dune;

for (int p=0; p<12; ++p) { QuadratureRule<double,2> const& qr = QuadratureRules<double,2>::rule(GeometryType(GeometryType::simplex,2),p); double integ = 0; for (size_t g=0; g<qr.size(); ++g) { FieldVector<double,2> const& x = qr[g].position(); double weight = qr[g].weight();

  // integrating x[0]^p instead yields correct result!
  // Probably the y-coordinate of an integration point is wrong.
  integ += weight*std::pow(x[1],p);
}
double exact = 1.0/((p+2)*(p+1));
std::cout.precision(16);
std::cout << "p=" << p << " integ=" << integ << " (should be " << exact << ") error = " << integ-exact << "\n";

} return 0; }

This program produces the following output:

p=0 integ=0.5 (should be 0.5) error = 0 p=1 integ=0.16666666665 (should be 0.1666666666666667) error = -1.666666804567285e-11 p=2 integ=0.08333333333333331 (should be 0.08333333333333333) error = -1.387778780781446e-17 p=3 integ=0.05 (should be 0.05) error = 0 p=4 integ=0.03333333333333333 (should be 0.03333333333333333) error = 0 p=5 integ=0.02380952380952381 (should be 0.02380952380952381) error = 3.469446951953614e-18 p=6 integ=0.01788018253696759 (should be 0.01785714285714286) error = 2.303967982473024e-05 p=7 integ=0.01389590970262738 (should be 0.01388888888888889) error = 7.020813738494933e-06 p=8 integ=0.01111111111111111 (should be 0.01111111111111111) error = 1.734723475976807e-18 p=9 integ=0.009090909090909089 (should be 0.00909090909090909) error = -1.734723475976807e-18 p=10 integ=0.00757575757575758 (should be 0.007575757575757576) error = 4.336808689942018e-18 p=11 integ=0.006410256410256407 (should be 0.00641025641025641) error = -2.602085213965211e-18

In particular, order 6 and 7 are quite wrong, but also order 1 seems to be less accurate than one should like.

Martin Weiser