与高斯求积相关的 R 中复杂递归的意义丧失

Question

我不是程序员（这是我第一次 post 来这里），对浮点运算没有太多经验。如果我遗漏了一些明显的东西，我深表歉意。

我一直在尝试使用 here 中描述的一般方法来寻找具有自定义权重函数的高斯求积参数。当可以手动找到参数时，该方法有效，因为检查了少量点。

但是，对于大量的正交点，以数字方式计算参数是有意义的。这些矩可以通过超几何函数来表示，它由我在这里使用的快速收敛级数给出。

我计算必要参数 an 和 bn 的算法涉及显式查找多项式的系数并使用参考资料中提供的公式。最后我们有一个复杂的递归，其中涉及相当多的加减乘除。

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?

代码如下：

#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

我为 an 得到的输出是：

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

一个明显的问题可能是矩的计算，就像这里所做的那样，但是增加项的数量 N 没有帮助，更重要的是，使用精确的矩值没有帮助根本不改变输出：

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;

使用 R 来完成这个任务是我个人的偏好（也是一个学习机会），所以如果你认为我需要使用另一种语言，我想我会在 Mathematica 中这样做，在那里精度可以设置任意高。

Answer 1

在您引用的论文中：

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].

给定μ_k求解a_j和b_j的问题极其ill-conditioned，并且随着点数的增加呈指数级恶化。换句话说，µ_k 的微小变化（由于浮点数的精度有限）会导致相应的 a_j 和b_j.

要使用此方法获得准确的结果，有必要更准确地计算 µ_k。例如，您找到的论文发现有必要将 µ_k 计算到 n=64 的数千位精度。

与高斯求积相关的 R 中复杂递归的意义丧失

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

r

floating-accuracy

significant-digits