Maxima scripts do not produce exakt quadrature coefficients

While working on the redesign of the quadraturerules I found out that the automatically generated 1d rules are not as exact as assumed, i.e. although the points and weights are printed with 100 digits the actual precision is only about 16-17 digits. This might be fine for double precision calculations, but when changing the ctype to long double or Float128 or even GMP it is not applicable.

I'm not an expert in Maxima, so I don't see the error in the script. It seems that the actual calculation is performed in double precision only and just the output is in high precision. @christi I think, you have uploaded the maxima scripts. Any ideas?

The error can be seen by simply summing up the weights. It should be 1 but the difference is about 1.e-16 for rules with order >= 8 at least.

To be mode precise: the quadrature weights are of low precision.

changed the description

Can you point out a particular rule (flavor, order) with this problem? Does the problem occur with the existing checked-in values, with your newly computed values or both?

The checked-in values are only exact up to 1.e-16. I just run the test-quadrature with long double or Float128 type, to check the accuracy of the 1d rules and it failed.

Examples: GaussLobatto, order 8/9: weight[2]

checked-in: 0.1892374781489234 2322840362420....
correct:    0.1892374781489234 9015830640410....

After checking the other rules again, it seems that only the gausslobatto rules are defect. The Jacobi1 and Jacobi2 are not tested in test-quadrature at all (its a bit more complicated due to the non-constant weight function)

OK, after reading a lot in the maxima documentations, I found the error in the gausslobatto.mac script:

replace

weight(n,xx) := ratsimp(2/((n*(n+1))*legendre_p(n,xx)^2))$

by

weight(n,xx) := at(ratsimp(2/((n*(n+1))*legendre_p(n,x)^2)),x=xx)$

Test for the Jacobi1 and Jacobi2 rules are nevertheless missing. There, the weights do sum up the volume time integral over the weight functions and the quadrature rule can not so easy be compared to the analytical solution, since a different integrand is evaluated. But, since the weight functions of Jacobi rules are also polynomials, a simple shift in the call to the analytical solution could be implemented.

Is this something worth to have, or is it not really interesting to test also the Gauss-Jacobi Quadrature rules?

Here is a quick and dirty perl script that might serve as a prototype for a test:

#!/usr/bin/perl

my @weights;

print "scale=100\n";

while(<>) {
    if(m{// order (.*)}) {
        print "\"$1: \"\n";
    }
    elsif(m{_weight\[\d+\] = (\d\..*);}) {
        push @weights, $1;
    }
    elsif(m{break;}) {
        if(scalar @weights) {
            print join('+', @weights)."-(4/3)\n";
            @weights = ();
        }
    }
}

Use as

./checkrules jacobi_2_0_imp.hh |BC_LINE_LENGTH=0 bc -l

Sample output:

-*- mode: compilation; default-directory: "~/Projekte/dune-master/dune-geometry/dune/geometry/quadraturerules/" -*-
Compilation started at Thu Jul 26 16:02:34

./checkrules jacobi_2_0_imp.hh |BC_LINE_LENGTH=0 bc -l
0,1: -.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000003
2,3: .0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001
4,5: .0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000002
6,7: -.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000206
8,9: .0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000030
10,11: .000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000023657
12,13: .000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000112392
14,15: .000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000120906
16,17: -.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000828679
18,19: .0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000035362095
20,21: -.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000022861897
22,23: -.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000003650446242
24,25: .0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000020507702936
26,27: -.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000016176064725
28,29: .00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001469200816525
30,31: .00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000010454639065128
32,33: -.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000024494910738105
34,35: -.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000320995800245057
36,37: .00000000000000000000000000000000000000000000000000000000000000000000000000000000000000002476268249601934
38,39: -.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000010423851549865927
40,41: .00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000387127417858493
42,43: .000000000000000000000000000000000000000000000000000000000000000000000000000000000000004132984708850974012
44,45: .000000000000000000000000000000000000000000000000000000000000000000000000000000000000008719233322855529539
46,47: .000000000000000000000000000000000000000000000000000000000000000000000000000000000000245534157558266465984
48,49: .000000000000000000000000000000000000000000000000000000000000000000000000000000000000539542683530248277922
50,51: .000000000000000000000000000000000000000000000000000000000000000000000000000000000001923603291442789745442
52,53: .000000000000000000000000000000000000000000000000000000000000000000000000000000000004903631338550060750439
54,55: .000000000000000000000000000000000000000000000000000000000000000000000000000000000021944223974032840723347
56,57: .000000000000000000000000000000000000000000000000000000000000000000000000000000000323690059411289760889261
58,59: .000000000000000000000000000000000000000000000000000000000000000000000000000000001278197698887824369351691
60,61: -.000000000000000000000000000000000000000000000000000000000000000000000000000000024234874565768724128284712
0,1: 2,3: 4,5: 6,7: 8,9: 10,11: 12,13: 14,15: 16,17: 18,19: 20,21: 22,23: 24,25: 26,27: 28,29: 30,31: 32,33: 34,35: 36,37: 38,39: 40,41: 42,43: 44,45: 46,47: 48,49: 50,51: 52,53: 54,55: 56,57: 58,59: 60,61: 
Compilation finished at Thu Jul 26 16:02:34

For rules other than jacobi 2, the (4/3) has to be replaced by 1.

The test would require the presence of bc, of course, plus perl, or whatever else is used to do the extraction.

Yes, these two rules are correct. I meant, should we integrate the test for the correct weights and correct points into the test-quadrature.cc implementation? The, the checkWeights() and checkQuadrature functions need to be modified slightly to incorporate the factors for the weight-functions (e.g. 4/3 in the case of jacobi2)

Testing the jacobi quadrature is definitely useful, and if the means the test functions need additional weight parameters, then so be it.

I'm just wondering how to reliably incorporate high-precision tests. These do not necessarily need to check everything, but at least they should check the weights of the 1-d rules. (Since they are wrong and we have no test that caught that, we need at least a sort of regression test).

• We could use long double or Float128, but those aren't really that precise.
• We could use GMPField, but then the test (or GMPField) needs adjusting, because GMPField has no support for std::numeric_limits, which is used to extract the precision. (I'm not 100% sure how much sense std::numeric_limits::epsilon() makes for GMPField.)
• Both Float128 and GMPField have the disadvantage that they are not necessarily available.
• We could use some script to extract the values (or get them some other way) and feed them to some high-precision shell utilities such as bc or dc for checking. But we probably cannot rely on them since they are both priority optional in Debian stretch.

Considering there are no high-precision shell utilities whose presence we can rely on, it is probably best to just do the test using GMPField. But in a separate test executable so we can skip it if gmpxx was not found.

So long as at least some of the CI system has GMP that should be fine

The problem with the high-precision tests is:

• the values for the 1d rules are provided with double literals only, meaning
• long double and Float128 get only double-precision coefficients and thus have also not the desired accuracy.
• Float128 is not considered for the construction from string case, since numeric_limits::is_specialized==true, although std::is_floating_point==false. I would suggest, to change the bool-test to is_floating_point, because it provides just double values and nothing else.
• Also, I would suggest, to use at least long double-literals in the 1d rules, by appending an L to the values.
• The Float128 can, on the other hand, also not be used in the construction-from-string implementation, since Float128 has no string constructor, only a operator>>. An idea would be, to use a generic string-to-number conversion function (as I have provided in one of the quadrature-rule MRs already) and use this instead of a simple ct(const char[]) constructor.
• The rules for simplices (and maybe also for prisms) are not in high precision (yet). In my high-precision tests, I have used a loop over types, that are different for tensor-product rules and simplex rules.

Another comment on the construction-from-string implementation:

GMPField<n> value1("0.12345"); // fine
Float128 value2("0.12345"); // ERROR
long double value3("0.12345"); // ERROR

but

GMPField<n> value1(0.12345..); // fine, but loss of precision
Float128 value2(0.12345..); // fine, but loss of precision
long double value3(0.12345..); // fine, but loss of precision

What about a generic function toNumber<ct>(const char*) or lexicalCast<ct>(const char*) or somthing similar?

GMPField<n> value1 = lexicalCast<GMPField<n>>("0.12345"); // fine
Float128 value2 = lexicalCast<Float128>("0.12345"); // fine
long double value3 = lexicalCast<long double>("0.12345"); // fine

And this function can be specialized for user-defined types and is provided for PODs.

mentioned in merge request !104 (merged)

I agree that we should probably make the numeric literals long double, and also use them for initializing long double quadrature rules (i.e. adjust the criterion that selects between the two implementations).

Float128 was introduced recently, and probably should just be extended to be constructible from a string (this holds whether or not that constructor is necessary to get quadrature rules working).

I may warm up to a generic conversion function, although it seems to be unnecessary for the types discussed here. However:

• many existing conversion functions are locale-dependent, which we want to avoid. So do not use the name lexicalCast(). Also, do not use from_string(): while no such facility exists, to_string() does exist and is locale-dependent.
• from_chars() is locale_independent, but the interface is too complicated to be directly useful in our case. It may be useful as a backend for a conversion function, but it is so new that it is unlikely extended precision types would provide an overload for it.
• The parameter tree in dune-common has conversion functions that should be locale-independent, based on the standard stringstreams and str.imbue(std::locale::classic()).

I think the easiest way to fix this without getting sidetracked is to add the string constructor to Float128, fix the numeric literals to be long double, and fix the criterion selecting between string and numeric literals.

mentioned in merge request dune-common!564 (closed)

mentioned in commit 5e6c08ec

closed via merge request !104 (merged)

closed via commit e228c90c

mentioned in commit e228c90c

reopened

Again the same error occurs. The newly generated quadrature coefficients are different from the ones stored in file. And the quadrature test fails with the new generated coefficients. Was there a change in maxima that causes this new error? These maxima scripts are really frustrating.

Try the following:

• (install maxima, e.g. sudo apt install maxima)
• run in dune/geometry/quadraturerules: maxima --batch gausslobatto.mac
• Compile the tests cmake --build build-cmake --target test-quadrature
• Run the test: build-cmake/dune/geometry/test/test-quadrature

I'm using the latest maxima 5.45 and get errors in the test.

For some reason, maxima produces wrong values for the roots of the legendre polynomials, when using bigfloat datatype.

This seems to be a bug in maxima. Current version is 5.45.1,.Correct values are computed in 5.38.1. Since 5.39.0 some legendre polynomial roots are wrong.

This issue needs to be documented in the maxima scripts. I will add this in a MR.

closed