Невозможно построить обратный полином Ньютона - PullRequest
Новая
       18

Невозможно построить обратный полином Ньютона

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

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


import numpy as np
import matplotlib.pyplot as plt

a = -1
b = 1
N = 201
x = np.linspace(a, b, N)
n = [3, 7, 10, 20]

plotsPfw = [None, None, None, None]
fig1, ((plotsPfw[0], plotsPfw[1]),(plotsPfw[2], plotsPfw[3])) = plt.subplots(2,2, gridspec_kw = {'hspace': 0.5, 'wspace': 0.5})
plotsPbw = [None, None, None, None]
fig2, ((plotsPbw[0], plotsPbw[1]),(plotsPbw[2], plotsPbw[3])) = plt.subplots(2,2, gridspec_kw = {'hspace': 0.5, 'wspace': 0.5})

def Runge(x):
    return 1/(1+25*np.power(x,2))

def findif_forw(func, order, i):
    if order == 1:
        return (func[i+1] - func[i])
    else:
        return (findif_forw(func, order-1, i+1) - findif_forw(func, order-1, i))

def findif_backw(func, order, i):
    if order == 1:
        return (func[-1-i] - func[-1-i+1])
    else:
        return (findif_backw(func, order-1, i) - findif_backw(func, order-1, i-1))

def Newton_forw(xi, func, n):
    h = xi[1] - xi[0]
    N = func[0]
    q = (x - xi[0])/h
    Q = 1
    for i in range(1, n):
        Q *= (q - i + 1)/i
        N += findif_forw(func, i, 0)*Q
    return N

def Newton_backw(xi, func, n):
    h = xi[1] - xi[0]
    N = func[-1]
    q = (x - xi[-1])/h
    Q = 1
    for i in range(1, n):
        Q *= (q + i - 1)/i
        N += findif_backw(func, i, i)*Q
    return N

yRunge = Runge(x)

for i in range(0, len(n)):
    xi = np.linspace(a,b,n[i])
    yRungei = Runge(xi)
    yNewton_forw = Newton_forw(xi, yRungei, n[i])
    yNewton_backw = Newton_backw(xi, yRungei, n[i])
    plotsPfw[i].plot(x, yRunge, ':r', x, yNewton_forw, '-g', linewidth = 0.25)
    plotsPfw[i].set_title('n = {}'.format(n[i]))
    plotsPfw[i].set_xlabel('x')
    plotsPfw[i].set_ylabel('Runge, fwd P_{}(x)'.format(n[i]))
    plotsPbw[i].plot(x, yRunge, ':r', x, yNewton_backw, '-b', linewidth = 0.25)
    plotsPbw[i].set_title('n = {}'.format(n[i]))
    plotsPbw[i].set_xlabel('x')
    plotsPbw[i].set_ylabel('Runge, bkwd P_{}(x)'.format(n[i]))

plt.show()

Конечные различия кажутся нормальными, после того, как я проверил их вручную, проблема, вероятно, возникает из части q полинома. Это наиболее очевидно для случая n = 3 - в прямой формуле после шага 1 Q - восходящий (в направлении интерполяции) наклон, а после шага 2 Q - нисходящая парабола, смещенная влево, поэтомув дальнем правом конце они, кажется, несколько компенсируют друг друга, приводя к тому, что полиномиальный граф интерполяции проходит через целевой граф n раз.

В обратной формуле Q почти всегда нисходящий (опять же, в направленииинтерполяции), поэтому крайний левый конец становится крайне несбалансированным, что видно на графике. Я пробовал следующее:

  • обратный знак значения, возвращенного конечной разницей порядка 1

  • обратный знак формулы для q

  • используйте q + i + 1, q + i, q-i + 1 в полиноме Q

Я считаю, что любой метод интерполяции предполагается, что станет неуравновешенным по отношению к противоположному концу (и прямой многочлен неправильный), или прямой многочлен ведет себя так, как он должен (а обратный многочлен неправильный)

1 Ответ

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

Ваша обратная разница очень странная. Почему -i в разнице порядка 1? Вы можете исправить это и еще больше упростить обе процедуры, завершив их в порядке 0 как 

def findif_forw(func, order, i):
    return  func[i] if order == 0 else findif_forw(func, order-1, i+1) - findif_forw(func, order-1, i);

def findif_backw(func, order, i):
    return  func[i] if order == 0 else findif_backw(func, order-1, i) - findif_backw(func, order-1, i-1)

. Это позволит вычислить обратные различия со значениями из исходного индекса, возвращаясь назад order слагаемых. Чтобы получить значения из конца массива, просто вызовите его с помощью 

        N += findif_backw(func, i, -1)*Q

Обратите внимание, что это не очень эффективно, поскольку вы пересчитываете несколько промежуточных разностей.

Результат одинаков вперед и назад, как и следовало ожидать от симметрии функции. enter image description here

...