我不是程序员(这是我在这里的第一篇文章(，对浮点运算没有太多经验。如果我遗漏了一些显而易见的东西，我深表歉意。
我一直在尝试使用这里描述的通用方法，找到具有自定义权重函数的高斯求积的参数。当可以手动找到参数时，该方法在检查少量点时有效。
然而，对于大量的正交点，用数值计算参数是有意义的。矩可以通过超几何函数来表达，这是由快速收敛的级数给出的，我在这里使用它。
我的计算必要参数an和bn的算法包括明确地找到多项式的系数并使用参考文献中提供的公式。最后，我们有一个复杂的递归，它涉及到相当多的加法、减法、乘法和除法。
问题是：我非常确定，在我的情况下，所有的
an=0.5都是准确的。但我在R中所做的算法很快就丢失了数字，在第5步给出了0.4999999981034791707302。为了避免这个问题，我可以更改算法中的哪些内容？
这是代码：
#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中这样做，那里的精度可以设置得任意高。
在您引用的论文中：
然而，代数集的解系数aj和bj关于力矩µk的方程是极其病态的："即使在双倍精度下，在时间n=12时失去所有精度也并不罕见"[1]。
求解给定µk的aj和bj问题是极其病态的，并且随着点数的增加而呈指数级恶化。换句话说，µk的微小变化(由于浮点数的精度有限(会导致相应的aj和bj发生较大变化。
为了用这种方法获得准确的结果，有必要更准确地计算µk。例如，您所在的论文发现，对于n=64，有必要将µk计算到数千位数的精度。