Ошибка в функции вычисления интеграла методом Гаусса на 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()