Loss of significance in a complicated recurrence in R, related to Gaussian quadrature

Question

I am not a programmer (this is my first post here), and don't have much experience with floating point arithmetic. I apologize if I missed something obvious.

I have been trying to find the parameters for Gaussian quadrature with custom weight function, using the general method described for example here. The method works, as checked for small number of points, when the parameters can be found by hand.

However, for large number of quadrature points it makes sense to compute the parameters numerically. The moments can be expressed through a hypergeometric function, which is given by the quickly converging series, which I am using here.

My algorithm for computing the necessary parameters an and bn involves finding the coefficients of the polynomials explicitly and using the formulas provided in the reference. In the end we have a complicated recurrence, which involves quite a few additions, subtractions, multiplications and divisions.

The problem is: I am pretty certain that in my case all an=0.5 exactly. But the algorithm I have made in R quickly loses the digits, giving 0.4999999981034791707302 instead on the 5th step. What can I change in the algorithm to avoid this problem?

Here's the code:

#Moments of sin(pi x) on [0,1] (hypergeometric function)
FIm <- function(n,N){   z <- -pi^2/4;
            f <- 1;
            k <- 0;
            a <- (n+2)/2;
            b <- 3/2;
            c <- (n+4)/2;
            while(k < N){f <- 1+f*z*(N-1-k+a)/(N-k)/(N-1-k+b)/(N-1-k+c);
                    k <- k+1}
            return(f*pi/(n+2))};
#Number of quadrature points
Nq <- 5;
n <- 0:(2*Nq+1);
#Moments
mu <- FIm(n,35);
#Recurrence parameters
an <- rep(0,Nq+1);
bn <- rep(0,Nq+1);
sn <- rep(0,Nq+1);
#Initial values
sn[1] <- mu[1];
an[1] <- mu[2]/sn[1];
#Coefficients of the orthogonal polynomials
Ank <- matrix(rep(0,(Nq+1)^2), nrow = Nq+1, ncol = Nq+1, byrow=TRUE);
#Initial values
Ank[1,1] <- 1;
Ank[2,1] <- - an[1];
Ank[2,2] <- 1;
#Starting recurrence
nn <- 2;
while(nn <= Nq){#Computing the coefficients of the squared polynomial
        Blj <- outer(Ank[nn,], Ank[nn,], FUN = "*");
        Cj <- rep(0,2*nn-1);
        j <- 1;
        while(j <= nn){l <- j;
                       while(l <= nn){if(j==l){Cj[j+l-1] <- Cj[j+l-1]+Blj[j,l]} else{Cj[j+l-1] <- Cj[j+l-1]+2*Blj[j,l]};
                l <- l+1};
        j <- j+1};
        #Computing the inner products and applying the recurrence relations
        sn[nn] <- sum(Cj*mu[1:(2*nn-1)]);
        an[nn] <- sum(Cj*mu[2:(2*nn)])/sn[nn];
        bn[nn] <- sn[nn]/sn[nn-1];
        k <- 1;
        while(k <= nn+1){if(k>1){Ank[nn+1,k] <- Ank[nn+1,k]+Ank[nn,k-1]};
                Ank[nn+1,k] <- Ank[nn+1,k]-an[nn]*Ank[nn,k]-bn[nn]*Ank[nn-1,k];
        k <- k+1};
nn <- nn+1};
#Computing the coefficients of the squared polynomial
Blj <- outer(Ank[nn,], Ank[nn,], FUN = "*");
Cj <- rep(0,2*nn-1);
j <- 1;
while(j <= nn){l <- j;
    while(l <= nn){if(j==l){Cj[j+l-1] <- Cj[j+l-1]+Blj[j,l]} else{Cj[j+l-1] <- Cj[j+l-1]+2*Blj[j,l]};
    l <- l+1};
j <- j+1};
#Computing the inner products and applying the recurrence relations
sn[nn] <- sum(Cj*mu[1:(2*nn-1)]);
an[nn] <- sum(Cj*mu[2:(2*nn)])/sn[nn];
bn[nn] <- sn[nn]/sn[nn-1];
an

The output I get for an is:

[1] 0.5000000000000000000000 0.5000000000000004440892 0.4999999999999593103261
[4] 0.4999999999963960495286 0.4999999998869631423482 0.4999999981034791707302

An obvious problem could be the computation of the moments, as it's done here, but increasing the number of terms N doesn't help, and more importantly, using the exact values for the moments doesn't change the output at all:

mu[1] <- 2/pi;
mu[2] <- 1/pi;
mu[3] <- 1/pi-4/pi^3;
mu[4] <- 1/pi-6/pi^3;
mu[5] <- (48 - 12 pi^2 + pi^4)/pi^5;
mu[6] <- (120 - 20 pi^2 + pi^4)/pi^5;
mu[7] <- (-1440 + 360 pi^2 - 30 pi^4 + pi^6)/pi^7;
mu[8] <- (-5040 + 840 pi^2 - 42 pi^4 + pi^6)/pi^7;
mu[9] <- (80640 - 20160 pi^2 + 1680 pi^4 - 56 pi^6 + pi^8)/pi^9;
mu[10] <- (362880 - 60480 pi^2 + 3024 pi^4 - 72 pi^6 + pi^8)/pi^9;
mu[11] <- (-7257600 + 1814400 pi^2 - 151200 pi^4 + 5040 pi^6 - 90 pi^8 + pi^10)/pi^11;
mu[12] <- (-39916800 + 6652800 pi^2 - 332640 pi^4 + 7920 pi^6 - 110 pi^8 + pi^10)/pi^11;

Using R for this task is my personal preference (as well as a learning opportunity), so if you think I need to use another language, I guess I will just do this in Mathematica, where the precision could be set arbitrarily high.

Answer 1

In the paper you cited:

However, the solution of the set of algebraic equations for the coefficients a_j and b_j in terms of the moments µ_k is extremely ill-conditioned: “Even in double precision it is not unusual to lose all accuracy by the time n = 12” [1].

The problem of solving a_j and b_j given µ_k is extremely ill-conditioned, and gets exponentially worse with increased number of points. In other words, a tiny change in µ_k (due to limited precision of floating point numbers) leads to a large change in corresponding a_j and b_j.

To get accurate results with this method, it is necessary to compute µ_k significantly more accurately. The paper you sited, for example, finds it necessary to compute µ_k to an accuracy of thousands of digits for n=64.