Как рассчитать граничные точки числовой производной в питоне?

Question

Я пытаюсь написать функцию для получения производной от любой общей функции / массива чисел.В частности, я использую Центральная разностная формула .Вопрос в том, что я не могу вычислить граничные точки производной, так как в формуле центральной разности используются индексы, которые не имеют границ.Мой код ниже

import numpy as np
n = 20000 # number of points in array
xs = np.linspace(start=-2*np.pi, stop=2*np.pi, num=n) # x values
y = np.array([np.sin(i) for i in xs]) # our function, sine

def deriv(f, h):
    """
    Calauclate the numerical derivative of any function
    :param f: numpy.array(float), the array of numbers we differentiate
    :param h: step size
    :rtype d: numpy.array(float)
    """
    d = np.zeros_like(f)
    # this loop misses the first and last points in f
    for i in range(1, f.shape[0]-1):
        # 2-point formula
        d[i] = (f[i+1] - f[i-1])/(2*h)

    return d

h = abs(xs[0] - xs[1]) # step size
y1 = deriv(y, h) # first derivative
y2 = deriv(y1, h) # second derivative
y3 = deriv(y2, h) # third derivative

Когда я строю график y,y1,y2,y3, вы можете видеть, как он взрывается в конечных точках

Я попытался установить конечные точки для ближайших соседей в deriv, как показано ниже.Хотя это работает для производных низкого порядка (1-го и 2-го), оно начинает разрушаться при производных более высокого порядка (3-го и более).

...
d = np.zeros_like(f)
for i in range(1, f.shape[0]-1):
    d[i] = (f[i+1] - f[i-1])/(2*h)

d[0] = d[1]
d[-1] = d[-2]
...

Производная в середине, вдали от границ, рассчитывается нормально.Проблема с границами.

Как мне относиться к граничным условиям здесь?Будет ли другая схема численного дифференцирования работать лучше, чем схема с центральным различием?

РЕДАКТИРОВАТЬ: Я ищу общий метод для решения этой проблемы, а не только метод, который может быть примененк функции синуса или любой другой периодической функции, как я использовал, чтобы проиллюстрировать эту проблему здесь.

Answer 1

Вы можете использовать прямую и обратную схему дифференцирования порядка 2 для граничных точек.По сути, мы знаем, что

(f(x+h)-f(x))/h = f'(x) + h/2*f''(x) + O(h²)                       (I)

и

(f(x+2h) - 2f(x+h) + f(x))/h² = f''(x+h) + O(h²) = f''(x) + O(h)   (II)

Используйте последний, чтобы заменить член первого порядка второй производной, то есть вычислите (I)-h/2*(II), чтобы получить

(-1/2*f(x+2h) + 2*f(x+h) -3/2*f(x))/h = f'(x) + O(h²)

Обратите внимание, что ошибка O (h²) в первой производной в общем случае приведет к ошибке O (h) во второй итерации разделенных разностей и O (1) в третьей.Кто-то может утверждать, что условия ошибки соответственно отменяются, но это произойдет только для внутренних точек, односторонние производные будут «портить» эту схему при увеличении расстояния от границы.

Answer 2

Это скорее вопрос численных методов, чем вопрос программирования.

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

f_periodic = np.zeros(f.size+2)
f_periodic[1:-1], f_periodic[0], f_periodic[-1] = f, f[-1], f[0]

Теперь вы можете дифференцировать f_periodic, для которого d[1] и d[-2] будут правильными производными значениями на границах (без учета d[0] и d[-1]).

Редактировать после новых требований OP ...

Для более общих граничных условий, скажем, конкретное значение на границах, существуют разные подходы, которым можно следовать:

1. Использовать значения-призраки:

Снова расширить функцию и экстраполировать значения для новых границ.В зависимости от порядка численной дифференциации потребуется больше ячеек-призраков.Для текущей схемы будет использоваться простая линейная экстраполяция (требуется только 1 призрачное значение на каждой границе):

f_new = np.zeros(f.size+2)
f_new[1:-1] = f
f_new[-1] = f[-2] + (f[-2]-f[-3])/(x[-2]-x[-3])*(x[-1]-x[-2])
f_new[0] = f[1] + (f[1]-f[2])/(x[1]-x[2])*(x[0]-x[1])

Обратите внимание, что вам также необходимо расширить x.Однако, поскольку у вас постоянный интервал, просто используйте h вместо пространственных различий, например, x[-2]-x[-3].Теперь вы можете дифференцировать f_new, и вы получите приближение 1-го порядка от производной на границах (поскольку вы использовали линейную экстраполяцию, чтобы найти значение-призрак).

Использование схем прямого и обратного на границах

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