You may want to crosspost this to the mpmath list. Aaron Meurer
On Thu, Jan 14, 2016 at 10:14 AM, Tobias Hartung <[email protected]> wrote: > Hi everyone, > > I am trying to integrate polynomials on an n-sphere and implemented (see > attached) the algorithm published in Alan Genz "Fully Symmetric > Interpolatory Rules for Multiple Integrals over Hyper-Spherical Surfaces" > Journal of Computational and Applied Mathematics 157: 187-195, 2003. On > double precision, everything is fine; however, I lose precision in some > integrals using more than double prec. > > I am using randomized samples on the sphere which (by construction) should > integrate the polynomials (that I tested) exactly. However, integrating all > polynomials up to degree 6 on a 5-dimensional sphere with points that should > integrate all polynomials up to degree 7 exactly, I obtain the following > test results : > > mp.prec = 1024 (aka mp.dps = 307) > > polynomial | degree of polynomial | log_10 ( difference of two > randomized integrations ) | log_10 ( relative error of weights ) | > analytic value of integral > > [0, 0, 0, 0, 0, 0] 0 -inf -306.962262 > 31.0062767 > [1, 0, 0, 0, 0, 0] 1 -308.111506 -306.962262 0.0 > [2, 0, 0, 0, 0, 0] 2 -14.519354 -306.962262 5.16771278 > [1, 1, 0, 0, 0, 0] 2 -14.8787358 -306.962262 0.0 > [3, 0, 0, 0, 0, 0] 3 -308.583548 -306.962262 0.0 > [2, 1, 0, 0, 0, 0] 3 -309.157806 -306.962262 0.0 > [1, 1, 1, 0, 0, 0] 3 -309.759866 -306.962262 0.0 > [4, 0, 0, 0, 0, 0] 4 -14.6442927 -306.962262 1.93789229 > [3, 1, 0, 0, 0, 0] 4 -15.3047045 -306.962262 0.0 > [2, 2, 0, 0, 0, 0] 4 -15.0673319 -306.962262 0.645964098 > [2, 1, 1, 0, 0, 0] 4 -16.0594155 -306.962262 0.0 > [1, 1, 1, 1, 0, 0] 4 -31.8146065 -306.962262 0.0 > [5, 0, 0, 0, 0, 0] 5 -309.37189 -306.962262 0.0 > [4, 1, 0, 0, 0, 0] 5 -309.481766 -306.962262 0.0 > [3, 2, 0, 0, 0, 0] 5 -309.739262 -306.962262 0.0 > [3, 1, 1, 0, 0, 0] 5 -310.095075 -306.962262 0.0 > [2, 2, 1, 0, 0, 0] 5 -310.474598 -306.962262 0.0 > [2, 1, 1, 1, 0, 0] 5 -310.714154 -306.962262 0.0 > [1, 1, 1, 1, 1, 0] 5 -311.46643 -306.962262 0.0 > [6, 0, 0, 0, 0, 0] 6 -14.7692314 -306.962262 0.968946146 > [5, 1, 0, 0, 0, 0] 6 -15.6057345 -306.962262 0.0 > [4, 2, 0, 0, 0, 0] 6 -15.4314091 -306.962262 0.193789229 > [3, 3, 0, 0, 0, 0] 6 -15.8275833 -306.962262 0.0 > [4, 1, 1, 0, 0, 0] 6 -16.5822942 -306.962262 0.0 > [3, 2, 1, 0, 0, 0] 6 -16.4753651 -306.962262 0.0 > [2, 2, 2, 0, 0, 0] 6 -16.0997045 -306.962262 > 0.0645964098 > [3, 1, 1, 1, 0, 0] 6 -32.3374853 -306.962262 0.0 > [2, 2, 1, 1, 0, 0] 6 -17.3552924 -306.962262 0.0 > [2, 1, 1, 1, 1, 0] 6 -33.4100353 -306.962262 0.0 > [1, 1, 1, 1, 1, 1] 6 -49.1751847 -306.962262 0.0 > > It should be noted that the polynomial are written in the following form: a > list p of length n represents the polynomial > > z_0^{p[0]} z_1^{p[1]} z_2^{p[2]} ... z_{n-1}^{p[n-1]}. > > If we look at the second to last column, we can see that the weights add are > exact with respect to mp.prec=1024 but do not add up to 1.0 exactly. Hence, > the loss in precision is unlikely to be caused by the weights being only on > double prec (and even if it is, then I don't understand why because they are > calculated with mp.prec=1024, as well). > > The really interesting column is the middle one. Here, I took two randomized > set of integration points (both should integrate all polynomials exactly up > to machine error) and printed the decadic logarithm of the absolute value of > their difference. In other words, we are expecting the middle column to be > populated with numbers in the vicinity of -307 (just like the fourth > column). However, we only get these values for odd degree polynomials (and > the constant 1 but that is just the sum of weights, i.e., the randomization > doesn't do anything). For even degree polynomials, we have significantly > reduced precision. > > Does anyone have an idea what might be the reason for such behavior? > > Thank you very much, > Tobias > > -- > You received this message because you are subscribed to the Google Groups > "sympy" group. > To unsubscribe from this group and stop receiving emails from it, send an > email to [email protected]. > To post to this group, send email to [email protected]. > Visit this group at https://groups.google.com/group/sympy. > To view this discussion on the web visit > https://groups.google.com/d/msgid/sympy/cf76f223-96a0-4789-8530-db224622f6bb%40googlegroups.com. > For more options, visit https://groups.google.com/d/optout. -- You received this message because you are subscribed to the Google Groups "sympy" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. To post to this group, send email to [email protected]. Visit this group at https://groups.google.com/group/sympy. To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAKgW%3D6LnZcT2SChLR3B_fXKt95eoH7xXQANkL8xwVwWk1csaBg%40mail.gmail.com. For more options, visit https://groups.google.com/d/optout.
