Разложение котангенса в ряд Тейлора и вычисление его приближенного значения

Question

Необходимо разложить котангенс в ряд Тейлора и вычислить его приближенное значение. Я написал программу, которая вычисляет эти значения в интервале от -2.4 до 2.4 , но меня это не устраивает, поскольку передо мной стоит задача - решить от -pi до pi. Помогите разобраться, что здесь не так.

#include <string>
#include <iomanip>
using namespace std;
constexpr auto PI = 3.14159265 ;
#include <cmath>
#include <iostream>
using namespace std;
 // Факториал x!
long double fact(int x)
{ 
    long double x0 = x / 1.0;
    if (x <= 1)
        return 1;
    else {
        long double fact = 1;
        for (; x != 0; x--) {
            fact *= x0;
            x0--;
        }
        return fact;
    }
}

// Биноминальный коэффициент Ньютона C(n/k)
long double comb(long long  n, long long  k)
{
    long double comb = 1;
    for (long i = 1; i <= k; i++) {
        comb *= (n + 1) - i;
        comb /= i;
    }
    return comb;
}
// Функция
long double function(long double x, int n, long double bn[]) {
    long double F = pow(-1, n)*pow(2.0, 2 * n)*bn[2 * n] * pow(x, 2 *n-1) / fact(2 * n);
    return F;
}
// Массив чисел Бернулли
void bernoulliNumbers(long double bn[],long long n){
    while (n != 225) {
        n++;
        bn[0] = 1;
            long double s = 0;
            for (long k = 1; k <= n; k++) {
                if ((n != 1) && (n % 2 != 0)) {
                    bn[n] = 0;
                    break;
                }
                s += comb(n + 1, k + 1)*bn[n - k];
            }
            bn[n] = (-1* s) / (n + 1);
        }
    }

int main() {
    long double xStart, xFinish, epsilon, dX, xIdeal, summa, x;
    cout << endl << "Input begin:";
    cin >> xStart;
    cout << "Input end:";
    cin >> xFinish;
    cout << endl << "Input dX";
    cin >> dX;
    cout << endl << "Input toch";
    cin >> epsilon;
    x = xStart;
    summa = 0;
    long long n = 0;
    long double bn[226];
    bn[0] = 1;
    bernoulliNumbers(bn, n);
    while (x <= xFinish) {  // Исключить зацикливание xStart = xFinish // x и эпсилон заменить
        for (int n = 0; (fabs(function( x, n, bn)) >= epsilon )  && (x != 0) &&  (fabs(x) < PI ) ; n++ ) {
            summa = summa + function(x,n, bn);
            cout << bn[n] << "\t" << function(x, n, bn) << "\t" << fact(2 * n) << " \t |" << pow(x, 2 * n - 1) << "\t" << pow(2.0, 2 * n) << endl;
        } 
        cout << endl;
        if (summa != 0) {
            long double MyFunction = summa;
            cout << MyFunction << " | " << 1 / tan(x) << endl;
        }
        x = x + dX; // В конце цикла 
        summa = 0;
        if (x == xFinish) {
            break;
        }
    }
    return 0;
}

Answer 1

Ваша проблема - в переполнении из-за способа вычисления: явно считать факториалы и степени, а делить потом. В результате вы вылетаете на пределы представления факториала используемым типом (конкретнее - у меня на VC++ 2019 при n == 86), и, соответственно получаете некорректный результат (вернее, NaN - не число).

Вам нужно изменить сам принцип расчета очередного члена, отталкиваясь, например, от отношения соседних членов ряда (посмотрите, на ruSO этот метод встречается очень часто).

Дополнение.

Признаю, сначала мне казалось, что проблема в факториалах, но нет. Проблема в самом числе Бернулли. Дело в том, что эти числа очень быстро растут,

так что формально ряд, который при этом при больших n превращается в

и сходится при

имеет огромные промежуточные члены, выпадающие за возможности представления типом double. Так что на значениях, близких к "краям", посчитать просто так, не прибегая к особым ухищрениям, не получится...

Вот код покороче, но он все равно слетает - не потому, что ряд не сходится, но в силу ограниченности представления математики в компьютерах :(

double bernoulli(int n)
{
    static vector<double> b{1.0, -0.5, 1.0/6.0, 0.0};
    if (n < b.size()) return b[n];
    if (n%2 && n > 2) return 0;
    double sum = 0;
    double a = -1;
    for(int k = 1; k <= n; ++k)
    {
        a = a*(n-k+1)/(k+1);
        sum += a*bernoulli(n-k);
    }
    b.resize(n+1);
    b[n] = sum;
    return sum;
}

double series(double x, double eps)
{
    double sum = 1, a = 1, z = 1;
    for(int n = 1; abs(z) > eps; ++n)
    {
        a = -a * 2*x*x/n/(2*n-1);
        z = a * bernoulli(2*n);
        sum += z;
    }
    return sum/x;
}

int main()
{
    for(double x = -3.1; x < 3.11; x += 0.2)
    {
        cout << setw(5) << x << setw(10) << 1/tan(x) << setw(10) << series(x,1e-7) << endl;
    }
}

Дополнение 2.

Проблема в вычислениях чисел Бернулли по рекурсивной формуле. Если взять функцию из boost - работает на ура... Но, откровенно говоря, разбираться в том, как они считают, не очень тянет.

Но применение аппроксимации для больших n уже существенно расширяет диапазон - можете попробовать:

#include <vector>
#include <iostream>
#include <iomanip>
#include <cmath>

using namespace std;

double bernoulli(int n)
{
    static vector<double> b{1.0, -0.5, 1.0/6.0, 0.0};
    if (n < b.size()) return b[n];
    if (n%2 && n > 2) return 0;
    if (n < 20)
    {
        double sum = 0;
        double a = -1;
        for(int k = 1; k <= n; ++k)
        {
            a = a*(n-k+1)/(k+1);
            sum += a*bernoulli(n-k);
        }
        b.resize(n+1);
        b[n] = sum;
        return sum;
    }
    else
    {
        const double pi = 3.14159265358, e = 2.718281828459045;
        return ((n/2%2)*2-1)*4*pi*pow(n/(2.*pi*e),0.5+n)*
            exp(0.5+1./(12.*n)-1./(360.*n*n*n)+1./(1260.*n*n*n*n*n));
    }
}

double series(double x, double eps)
{
    double sum = 1, a = 1, z = 1;
    for(int n = 1; abs(z) > eps; ++n)
    {
        a = -a * 2*x*x/n/(2*n-1);
        z = a * bernoulli(2*n);
        sum += z;
    }
    return sum/x;
}

int main()
{

    for(double x = -3.1; x < 3.11; x += 0.05)
    {
        double y = series(x,1e-9);
        cout
            << setw(5) << x << setw(10) << 1/tan(x)
            << setw(10) << y << endl;
    }
}

БЛОГ НА HUSL

Разложение котангенса в ряд Тейлора и вычисление его приближенного значения

Ответы (1 шт):