Ошибка в числовой производной 4-го порядка вблизи конечных точек интервала

Question

Я интегрирую уравнение в частных производных, в котором у меня есть частная производная четвертого порядка по $ x $, а численное интегрирование по времени дает мне абсурдные ошибки. Я полагаю, что причиной проблемы является то, что я получаю большие ошибки в производной четвертого порядка. Чтобы проиллюстрировать это, я беру численные производные функции $ y (x) = 1- \ cos (2 \ pi x) $. Ниже я строю производные $ y_ {xx} (x) $ и $ y_ {xxxx} (x) $ в области $ (0, 1.0) $. Смотрите рисунки ниже:

Как можно видеть, ошибки происходят в основном вблизи границ.

Численные производные были выполнены с использованием метода числового градиента. Код Python находится здесь:

import numpy as np
import matplotlib.pyplot as plt
N = 64
X = np.linspace(0.0, 1.0, N, endpoint = True)
dx = 1.0/(N-1)
Y= 1.0-np.cos(2*np.pi*X)
Y_x = np.gradient(Y, X, edge_order = 2)
Y_xx = np.gradient(Y_x, X, edge_order = 2)
plt.figure()
plt.title("$y(x)=1-\cos(2\pi x)$")
plt.xlabel("$x$")
plt.ylabel("$y_{xx}(x)$")
plt.plot(X, ((2*np.pi)**2)*np.cos(2*np.pi*X), 'r-', label="analytics")
plt.plot(X, Y_xx, 'b-', label="numerics")
plt.legend()
plt.grid()
plt.savefig("y_xx.png")
Y_xxx = np.gradient(Y_xx, X, edge_order = 2)
Y_xxxx = np.gradient(Y_xxx, X, edge_order = 2)
plt.figure()
plt.title("$y(x)=1-\cos(2\pi x)$")
plt.xlabel("$x$")
plt.ylabel("$y_{xxxx}(x)$")
plt.plot(X, -((2*np.pi)**4)*np.cos(2*np.pi*X), 'r-', label="analytics")
plt.plot(X, Y_xxxx, 'b-', label="numerics")
plt.legend()
plt.grid()
plt.savefig("y_xxxx.png")
plt.show()

Мой вопрос заключается в том, как алгоритмически уменьшить эту большую ошибку на границах за пределами очевидного увеличения N? Как это сходимость производной четвертого порядка не является равномерным. Можно ли сделать это равномерная? Возможно, использование экстраполяции на границах подойдет.

Answer 1

Здесь я использую экстраполяцию на границах методом NumPy Polyfit.

# Quadratic extrapolation of the left end-points of Y_xxxx
x = X[2:5]
y = Y_xxxx[2:5]
z = np.polyfit(x, y, 2)
p = np.poly1d(z)
for i in range(2):
    Y_xxxx[i] = p(X[i])

Кроме того, у меня есть меньшая сетка на границах, чтобы уменьшить погрешность числовой производной (Y_xxxx) без необходимости увеличивать число точек сетки равномерно по всей области интегрирования.

# grid with variable spacing
dx = 1.0/N
X1 = np.linspace(0.0, 4*dx, 16)
X2 = np.linspace(4*dx, 1.0-4*dx, N)
X3 = np.linspace(1.0-4*dx, 1.0, 16)
X= np.concatenate([X1, X2, X3])

Поскольку шаг интегрирования не является постоянным, я продолжаю использовать метод пустого градиента, поскольку он может учитывать это изменение при численном расчете производных.

Вот мой не очень питоническийкод для сравнения результатов с и без переменной сетки:

import numpy as np
import matplotlib.pyplot as plt
N = 512
X = np.linspace(0.0, 1.0, N, endpoint = True)
Y= 1.0-np.cos(2*np.pi*X)
Y_x = np.gradient(Y, X, edge_order = 2)
Y_xx = np.gradient(Y_x, X, edge_order = 2)
# Quadratic extrapolation of the left end-points of Y_xx
x = X[2:5]

y = Y_xx[2:5]
z = np.polyfit(x, y, 2)
print "z=", z
p = np.poly1d(z)
Y_xx[0] = p(X[0])
Y_xx[1] = p(X[1])

# Quadratic extrapolation of the right end-points of Y_xx
x = X[-5:-2]
y = Y_xx[-5:-2]
z = np.polyfit(x, y, 2)
p = np.poly1d(z)
Y_xx[-1] = p(X[-1])
Y_xx[-2] = p(X[-2])

Y_xxx = np.gradient(Y_xx, X, edge_order = 2)
Y_xxxx = np.gradient(Y_xxx, X, edge_order = 2)
# Quadratic extrapolation of the left end-points of Y_xxxx
x = X[2:5]
y = Y_xxxx[2:5]

z = np.polyfit(x, y, 2)
print z
p = np.poly1d(z)
for i in range(2):
    Y_xxxx[i] = p(X[i])
# Quadratic extrapolation of the right end-points of Y_xxxx
x = X[-5:-2]
y = Y_xxxx[-5:-2]

z = np.polyfit(x, y, 2)
#print z
p = np.poly1d(z)
for i in [-1, -2]:
    Y_xxxx[i] = p(X[i])
plt.figure()

plt.title("$y(x)=1-\cos(2\pi x)$")
plt.xlabel("$x$")
plt.ylabel("$y_{xxxx}(x)$")
plt.plot(X, -((2*np.pi)**4)*np.cos(2*np.pi*X), 'r-', label="analytics")
plt.plot(X, Y_xxxx, 'b-', label="numerics")
plt.legend()
plt.grid()

plt.savefig("y_xxxx.png")
plt.figure()
diff = Y_xxxx+((2*np.pi)**4)*np.cos(2*np.pi*X)
logErr = 0.5*np.log10(diff**2)
plt.plot(X, logErr, 'b-', label = "Fixed spacing")

# grid with variable spacing
dx = 1.0/N
X1 = np.linspace(0.0, 4*dx, 16)
X2 = np.linspace(4*dx, 1.0-4*dx, N)
X3 = np.linspace(1.0-4*dx, 1.0, 16)
X= np.concatenate([X1, X2, X3])

Y= 1.0-np.cos(2*np.pi*X)
Y_x = np.gradient(Y, X, edge_order = 2)
Y_xx = np.gradient(Y_x, X, edge_order = 2)
# Quadratic extrapolation of the left end-points of Y_xx
x = X[2:5]
y = Y_xx[2:5]

z = np.polyfit(x, y, 2)
print "z=", z
p = np.poly1d(z)
Y_xx[0] = p(X[0])
Y_xx[1] = p(X[1])
# Quadratic extrapolation of the right end-points of Y_xx

x = X[-5:-2]
y = Y_xx[-5:-2]
z = np.polyfit(x, y, 2)
p = np.poly1d(z)
Y_xx[-1] = p(X[-1])
Y_xx[-2] = p(X[-2])

Y_xxx = np.gradient(Y_xx, X, edge_order = 2)
Y_xxxx = np.gradient(Y_xxx, X, edge_order = 2)
# Quadratic extrapolation of the left end-points of Y_xxxx
x = X[2:5]
y = Y_xxxx[2:5]
z = np.polyfit(x, y, 2)
p = np.poly1d(z)
for i in range(2):
    Y_xxxx[i] = p(X[i])
# Quadratic extrapolation of the right end-points of Y_xxxx
x = X[-5:-2]
y = Y_xxxx[-5:-2]
z = np.polyfit(x, y, 2)
#print z
p = np.poly1d(z)
for i in [-1, -2]:
    Y_xxxx[i] = p(X[i])
diff = Y_xxxx+((2*np.pi)**4)*np.cos(2*np.pi*X)
logErr = 0.5*np.log10(diff**2)
plt.plot(X, logErr, 'r.-', label = "variable spacing")
plt.title("$log_{10}(|Error(x)|)$, N=%d" % (N))
plt.xlabel("$x$")
plt.xlim(0., 1.)
plt.legend()
plt.grid()
figure = "varStepLogErr.png"
print figure
plt.savefig(figure)
plt.show()

Вот сравнение с двумя методами (оба с экстраполяцией на границах), один из них с переменными шагами.

Здесь ошибка (x) = log_10 | Y_xxxx (x) + (2 * pi) ** 4 * cos (kx) |

Answer 2

Основная математическая проблема

Формула для 1-й производной, которую использует np.gradient, имеет ошибку O(h**2) (где h - ваш dx). Поскольку вы продолжаете принимать деривативы, это может быть плохо, поскольку изменение функции на величину z может изменить ее дериват на z/h. Обычно этого не происходит, потому что ошибка исходит от производных функции более высокого порядка, которые сами изменяются плавно; поэтому последующее дифференцирование «дифференцирует» большую часть ошибки предыдущего шага, а не увеличивает ее на 1/h.

Тем не менее, у границы другая история, где мы должны перейти от одной формулы конечных разностей (центрированная разница) к другой. Формула границы также имеет ошибку O(h**2), но она отличается O(h**2). И теперь у нас есть проблема с последующими этапами дифференциации, каждый из которых может дать коэффициент 1/h. Наихудший сценарий для 4-й производной - O(h**2) * (1/h**3), который, к счастью, здесь не реализуется, но вы все равно получаете довольно плохой граничный эффект. Я предлагаю два разных решения, которые работают примерно одинаково (второе несколько лучше, но дороже). Но сначала давайте подчеркнем значение того, что Уоррен Векессер сказал в комментарии:

Если функция периодическая , вы расширяете ее по периодичности, например, np.tile(Y, 3), вычисляете производную от по и берете среднюю часть, обрезая любые граничные эффекты ,

Прямое вычисление 4-й производной

Вместо того чтобы применять формулу конечных разностей для 1-й производной четыре раза, примените ее для 4-й производной один раз. Все еще будет проблема граничных значений, для которой моя лучшая идея - самая простая, постоянная экстраполяция. То есть:

Y_xxxx = (Y[4:] - 4*Y[3:-1] + 6*Y[2:-2] - 4*Y[1:-3] + Y[:-4])/(dx**4)
Y_xxxx = np.concatenate((2*[Y_xxxx[0]], Y_xxxx, 2*[Y_xxxx[-1]]))

с результатом

Если вам не нравятся маленькие плоские кусочки сбоку, возьмите меньшее dx (они имеют размер 2*dx).

Дифференцировать интерполирующий сплайн

Установите сплайн 5-й степени для данных, затем возьмите его 4-ю производную аналитически . Это требует SciPy и, вероятно, намного медленнее, чем прямой подход.

from scipy.interpolate import InterpolatedUnivariateSpline
spl = InterpolatedUnivariateSpline(X, Y, k=5)
Y_xxxx = spl.derivative(4)(X)

Результат:

Неоднородная сетка

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

def cheb_nodes(a, b, N):
    jj = 2.*np.arange(N) + 1
    x = np.cos(np.pi * jj / 2 / N)[::-1]
    x = (a + b + (b - a)*x)/2
    return x

N = 64
X = cheb_nodes(0, 1, N)

Остальное продолжается, как и раньше, используя InterpolatedUnivariateSpline степени 5. График самого 4-го производного теперь не показывает видимой разницы между числами и аналитикой, поэтому вместо этого приведен график Y_xxxx - analytics:

Ошибка не более 3, и она больше не максимальна на границе. Максимальная ошибка с равномерной сеткой была около 33.

Я также исследовал возможность наложения фиксированного состояния на интерполирующий сплайн для дальнейшего повышения его точности. Возможно, make_interp_spline может сделать это с

l, r = [(0, 0), (1, 0)], [(0, 0), (1, 0)]
spl = make_interp_spline(X, Y, k=5, bc_type=(l, r))

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