Почему использование формулы ряда Тейлора для арксинуса так быстро сходится при x = 1 - PullRequest
Новая
       18

Почему использование формулы ряда Тейлора для арксинуса так быстро сходится при x = 1

0 голосов
/ 16 октября 2019

Я взял этот пример из видео на YouTube, где sen (x) разрешается с помощью этого подхода, я просто изменил формулу на arcsin (x).

Я пытаюсь выяснить, почемуэтот код работает не так, как ожидалось: 

 #include <stdio.h>
    //#include <float.h>
    #include <math.h>

    int factorial(int x)
    {
        int fac = 1;

        while (x!=0)
        {
            fac= fac*x;
            x--;
        }
        return fac;
    }

    int main()
    {
        int i;
        float x, sum, divisor = 1, dividendo = 1,  temp;

        printf("enter the value of X(in degrees): ");
        scanf("%f", &x);
        //x = x*3.141592/180;

        sum = 0;
        for(i=1; ; i++)
        {
            divisor = (factorial(2*i));
            dividendo = (pow(4,(i)) * pow(factorial(i),2)*(2*i+1));
            temp = (divisor/dividendo)*pow(x,(2*i+1));
            if (temp < FLT_EPSILON)
                          break;
            sum = sum + temp;
            printf("our result = %f\n", sum);
        }  
        printf("reference = %f\n", asin(x));
       // printf("our result = %f\n", sum);

    }

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

enter the value of X(in degrees): 1
our result = 0.166667
our result = 0.241667
our result = 0.286310
our result = 0.316691
our result = 0.339064
our result = 0.356416
our result = 0.356621
our result = 0.356622
our result = 0.356622
our result = 0.356622
**reference = 1.570796**

1 Ответ

4 голосов
/ 16 октября 2019

Ряд Тейлора для арсцина равен 

arcsin(x) = sum n = 0; inf; (2*n)! / (4**n * (n!)**2 * (2*n + 1)) * x**(2*n + 1)

, где ** - оператор степени, а ! обозначает факториал. Как уже отмечалось в комментариях, int может представлять факториалы только до 12!;64-битный long может представлять факториалы вплоть до 20!. Поскольку вы используете (2*n)!, эти пределы будут достигнуты быстро.

Лучшее решение этой проблемы - не вычислять каждый член отдельно, как указано в формуле, а развивать его, вычисляя его из предыдущеготермин: 

term(n + 1) = fact * term(n)

Поскольку каждый термин является одним большим фактором, вы можете сделать это для каждого подфактора: 

       2*(n + 1))!  =  (2*n)!         * (2*n + 1) * (2*n + 2)
        4**(n + 1)  =  4**n           * 4
     ((n + 1)!)**2  =  n!**2          * (n + 1)**2
     2*(n + 1) + 1  =  2*(n + 1)      * (2*n + 3) / (2*n + 1)
x**(2*(n + 1) + 1)  =  x**(2*n + 1)   * x**2

Собрав все это вместе: 

float res = x;        // first term ...
float fact = x;       // ...equals the first factor

for (int n = 0; n < nMax; n++) {
    float old = res;

    // calculate term(n + 1) as per the formulas above

    fact *= (2*n + 1) * (2*n + 2);
    fact /= 4.0 * (n + 1)*(n + 1) * (2*n + 3);
    fact *= x * x * (2*n + 1);

    res += fact;
    if (res == old) break;

    printf("[%d] %f\n", n, res);
}  

printf("ref %f\n", asin(x));

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

Этот ряд не сходится хорошо около ± 1, где наклон функции приближается к бесконечности. По приближению Стирлинга n-й член приблизительно равен 1/((2*n+1)*√(πn)), который сходится очень медленно. Чтобы n-й член был меньше FLT_EPSILON, n должен быть больше 26000. К сожалению, при таком количестве суммирований конечная точность чисел с плавающей запятой не позволила бы серии сходиться к правильному ответу.

...