Ошибка интерполирования по Гауссу

Question

    from math import factorial

def Q(m, q):
    c = 1
    k = 1
    for i in range(m):
        if(i == 0):
            c *= q
        elif(i & 1):
          c *= q - (k - 1)
          k -= 1
        else:
            c *= q + (k - 1)
        k += 1
    return c

r = int(input('Укажите количество знаков после запятой: '))
X = [0.4, 0.7, 1.0, 1.3, 1.6, 1.9, 2.2]
Y1 = [round(i ** 5, r) for i in X] #Y1 = x ** 5
Y2 = [round(i ** -1, r) for i in X] #Y2 = x ** -1 

dY1 = [[round((Y1[i+1] - Y1[i]), r) for i in range(len(Y1)-1)]] 
# Дельты Y1
dY2 = [[round((Y2[i+1] - Y2[i]), r) for i in range(len(Y2)-1)]] 
# Дельты Y2
for i in range(1, len(X)-1):
    dY1.append([])
    dY2.append([])
    for j in range(len(X) - (i + 1)):
        dY1[i].append(round(dY1[i-1][j+1] - dY1[i-1][j], r))
        dY2[i].append(round(dY2[i-1][j+1] - dY2[i-1][j], r))

print('Xi', end='')
for i in range(len(X)): print('\t', X[i], end='')
print()

x = float(input('Введите x: '))
h = X[1] - X[0]
Px1 = Y1[len(Y1) // 2]
Px2 = Y2[len(Y2) // 2]
q = (x - X[len(X) // 2]) / h
k = 0
print(dY2)
for i in range(len(X)-1):
    Px1 += Q(i+1, q)*dY1[i][(len(X) // 2)-k] / factorial(i+1)
    Px2 += Q(i+1, q)*dY2[i][(len(X) // 2)-k] / factorial(i+1)
    if(i & 1): k -= 1
    k += 1
    
print('Для какой функции использовать первую формулу Гаусса:\n1. x^5\n2. x^-1')
a = int(input('Выбор: '))
if(a == 1): print('P(x) =', round(Px1, r))
elif(a == 2): print('P(x) =', round(Px2, r))

Почему-то для функции y = x^-1 интерполяционный многочлен считается неправильно. А чем может быть проблема, ведь для y = x^5 всё хорошо

Answer 1

у вас

def Q(m, q):
    c = 1
    k = 1
    for i in range(m):
        if(i == 0):
            c *= q
        elif(i & 1):
          c *= q - (k - 1)
          k -= 1
        else:
            c *= q + (k - 1)
        k += 1
    return c

не касаясь сущности задачи, точно не так должно быть:

def Q(m, q):
    c = 1
    k = 1

    for i in range(m):
        if i == 0:
            c *= q
        elif i & 1:
            c *= q - (k - 1)
            k -= 1
        else:
            c *= q + (k - 1)
            k += 1 # это находится внутри блока else
    
    return c

Answer 2

Так как при обратной функции 1/x в переменной q при больших значениях х, скапливается большое число, то я решил расширить диапазон интерполирования не превосходящий х+h, тем самым вычисляется всё очень хорошо

БЛОГ НА HUSL

Ошибка интерполирования по Гауссу

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