Ошибка в функции вычисления интеграла методом Гаусса на Python

Считаю интеграл методом Гаусса. В функции применяются справочные значения (берутся из таблицы) коэффициентов А и t. У меня они приведены для порядка от 1 до 6 (количество узлов, на которые дробится интервал). Больше задать нельзя.

Если узлов для заданной точности недостаточно (т.е. err_gauss > delta), то количество интервалов дробления удваивается, и для каждого считается свой интеграл с заданным пользователем n. Для этого применяется аргумент j_els = n-1. По идее потом надо сложить отдельные интегралы в общий показатель. Это реализовано в пункте 6.5.1 кода. А вот формула:

введите сюда описание изображения

Однако при n=1 получаю нулевое значение интеграла, а при n=2 - значение 0,85 (хотя должно быть близко к 1,75). В чем ошибка?

Код:

import numpy as np

def func_for_integrate(x):
    """
    Исходная функция для последующего интергрирования
    """
    res = (1 + np.sqrt(x)) / (x ** 2)
    return res


def gauss_integrated():
    """
    Функция получения интеграла по методу Гаусса
    """
    # 1. запросим порядок формулы, точность и границы интервала
    n = int(input('Введите порядок формулы (количество узлов от 1 до 6): '))
    delta = float(input('Введите точность расчета (десятичная дробь): '))
    a = 1
    b = 4

    # 2. определим условие, при котором n не может быть меньше единицы и больше пяти
    assert n >= 1 and n <= 6, 'Определите n от 1 до 6'

    # 3. определим максимаьное количество итераций удвоения интервалов
    # при доведении результата до заданной точности расчета
    # (это позволит избежать бесконечный цикл while)
    max_iter = 100

    # 4. определим словарь с коэффициентами t:
    # здесь ключи верхнего порядка (1,2,3 и т.д.) - это n из таблицы 5.1.
    # ключи второго порядка - i из таблицы,
    # значения ключей второго порядка коэффициенты t
    dict_of_koeffs_t = {
        1:
            {1: [0]},
        2:
            {1: [0.57735027], 2: [-0.57735027]},
        3:
            {1: [0.77459667], 2: [0], 3: [-0.77459667]},
        4:
            {1: [0.86113631], 2: [0.33998104], 3: [-0.33998104], 4: [-0.86113631]},
        5:
            {1: [0.90617985], 2: [0.53846931], 3: [0.53846931], 4: [-0.53846931], 5: [-0.90617985]},
        6:
            {1: [0.932469515], 2: [0.66120939], 3: [0.23861919], 4: [-0.23861919], 5: [-0.66120939], 6: [-0.932469515]}
    }

    # 5. определим словарь с коэффициентами А
    dict_of_koeffs_A = {
        1:
            {1: [2]},
        2:
            {1: [1], 2: [1]},
        3:
            {1: [(5 / 9)], 2: [(8 / 9)], 3: [(5 / 9)]},
        4:
            {1: [0.34785484], 2: [0.65214516], 3: [0.65214516], 4: [0.34785484]},
        5:
            {1: [0.23692688], 2: [0.47862868], 3: [0.56888889], 4: [0.47862868], 5: [0.23692688]},
        6:
            {1: [0.17132450], 2: [0.36076158], 3: [0.46791394], 4: [0.46791394], 5: [0.36076158], 6: [0.17132450]}
    }

    # 6. определим функцию расчета интеграла
    def gauss(n, j_els=None):
        # 6.1. по заданному пользователем n извлекаем в список
        # значения коэффициента t словаря соответствующего порядка n
        t_list = []
        for t in dict_of_koeffs_t[n].keys():
            t_list.append(dict_of_koeffs_t[n][t])

        # 6.2. по заданному пользователем n извлекаем в список
        # значения коэффициента A словаря соответствующего порядка n
        a_list = []
        for a_els in dict_of_koeffs_A[n].keys():
            a_list.append(dict_of_koeffs_A[n][a_els])

        # 6.3. формируем список кортежей, в которых сложены соответствующие t и А:
        # первый t - первый А, второй t - второй А и т.д.
        t_a_total_list = list(zip(t_list, a_list))

        # 6.4. определим переменную для параметра j
        if j_els is None:
            j_els = n - 1

        # 6.5. Рассчитаем интеграл с параметром j
        # 6.5.1. определим цикл для получения интеграла с параметром j
        sum_outside = 0  # первая сумма
        # Ниже диапазон: по количеству коэф-ов в зависимости от n
        for k_elems in range(1, n):
            sum_inside = 0  # вторая сумма
            # i_elems - кол-во точек рассечения отрезка [a,b]
            for i_elems in range(0, j_els):
                x_k_j = ((b + (a * (2 * j_els + 1))) / (2 * (j_els + 1))) + (
                            ((b - a) / (2 * (j_els + 1))) * t_a_total_list[k_elems][0][0])
                integral_func = func_for_integrate(x_k_j + (((b - a) / (j_els + 1)) * i_elems))
                expression_under_the_sum_inside = t_a_total_list[k_elems][1][0] * integral_func
                sum_inside += expression_under_the_sum_inside
            sum_outside += sum_inside

        # 6.5.2. получим значение интеграла по формуле
        s_j = (((b - a) / (2 * (j_els + 1))) * sum_outside)

        return s_j, j_els

    # 7.1 Рассчитаем интеграл в первый раз с параметром j = n-1
    I_prev, j_els = gauss(n)

    # 8. Определим переменную для количества итераций цикла while
    iter_count = 0

    # 9. Начнем цикл while для доведения точности расчета до заданной
    while True:
        iter_count += 1

        # 9.1 Рассчитаем новое значение интеграла с параметром j = 2*j_els+1
        I_next, j_els = gauss(n, j_els=2 * j_els + 1)

        # 11. Рассчитаем погрешность расчета
        err_gauss = abs(I_next - I_prev)

        # 12. Если погрешность меньше заданной точности или достигнуто максимальное количество итераций,
        # то выходим из цикла while
        if err_gauss <= delta or iter_count >= max_iter:
            scientific_format_ = "{:e}".format(err_gauss)
            break

        # 13. Иначе присваиваем новое значение интеграла предыдущему и повторяем цикл while
        I_prev = I_next

    # 14. Выводим результат на экран
    print(
        f'Значение интеграла: {I_prev}, количество итераций: {iter_count}, итоговое количество интервалов: {j_els}')
    print('Ошибка интегрирования: ', scientific_format_)


gauss_integrated()

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