Потеря значимости в сложном повторении в R, связанном с квадратурой Гаусса - PullRequest
Новая
       9

Потеря значимости в сложном повторении в R, связанном с квадратурой Гаусса

0 голосов
/ 21 ноября 2018

Я не программист (это мой первый пост здесь), и у меня нет большого опыта работы с арифметикой с плавающей запятой.Я извиняюсь, если пропустил что-то очевидное.

Я пытался найти параметры для квадратуры Гаусса с помощью пользовательской весовой функции, используя общий метод, описанный, например, здесь .Метод работает, как проверено для небольшого числа точек, когда параметры можно найти вручную.

Однако для большого количества квадратурных точек имеет смысл вычислять параметры численно.Моменты могут быть выражены через гипергеометрическую функцию, которая задается быстро сходящимся рядом, который я здесь использую.

Мой алгоритм вычисления необходимых параметров an и bn включает в себя нахождение коэффициентовполиномы явно и с использованием формул, приведенных в ссылке.В конце мы имеем сложное повторение, которое включает в себя довольно много сложений, вычитаний, умножений и делений.

Проблема в том, что я вполне уверен, что в моем случае все an=0.5 точно.Но алгоритм, который я сделал в R, быстро теряет цифры, давая 0.4999999981034791707302 вместо 5-го шага.Что я могу изменить в алгоритме, чтобы избежать этой проблемы?

Вот код: 

#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, где точность может быть установлена ​​произвольно высокой.

1 Ответ

0 голосов
/ 24 ноября 2018

В приведенной вами статье:

Однако решение системы алгебраических уравнений для коэффициентов a j и b j в терминахмоментов µ k крайне плохо обусловлен: «Даже с двойной точностью нет ничего необычного в том, чтобы потерять всю точность к моменту n = 12» [1].

Проблема решения j и b j при заданном µ k чрезвычайно плохо обусловлена ​​ и экспоненциально ухудшается с увеличением числаточки.Другими словами, незначительное изменение µ k (из-за ограниченной точности чисел с плавающей запятой) приводит к значительному изменению соответствующих a j и b j .

Чтобы получить точные результаты с помощью этого метода, необходимо значительно более точно вычислить µ k .Например, статья, которую вы разместили, считает необходимым вычислить µ k с точностью до тысяч цифр при n = 64.

...