在间隔的端点附近的4阶数值导数的误差

Question

我正在积分偏微分方程，其中我在$ x $中有一个四阶偏导数，时间上的数值积分给了我荒谬的误差。我相信问题的原因是我在四阶导数中得到了很大的误差。为了说明我采用函数$ y（x）= 1- \ cos（2 \ pi x）$的数值导数。下面我在域$（0,1.0）$中绘制衍生物$ y_ {xx}（x）$和$ y_ {xxxx}（x）$。见下图：

可以看出，错误主要发生在边界附近。

使用numpy梯度法进行数值导数。 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

Underlying mathematical issue

np.gradient使用的第一个导数的公式有误差O(h**2)（其中h是你的dx）。当你继续使用衍生品时，这有可能是坏的，因为用z量改变一个函数可以改变它的衍生物z/h。通常这不会发生，因为错误来自函数的高阶导数，它们本身变化平滑;因此，随后的区分“区分”前一步骤的大部分错误，而不是通过1/h放大它。

然而，在边界附近的故事是不同的，我们必须从一个有限差分公式（居中差异）切换到另一个。边界公式也有O(h**2)错误，但它是一个不同的O(h**2)。现在我们确实遇到了后续分化步骤的问题，每个分析步骤都可以贡献一个1/h因子。第四个导数的最坏情况是O(h**2) * (1/h**3)，幸运的是这里没有实现，但你仍然会得到一个非常糟糕的边界效应。我提供两种不同的解决方案，它们的表现大致相同（第二种方法略胜一筹但价格更贵）。但首先，让我们强调Warren Weckesser在评论中所说的含义：

如果函数是周期性的，则通过周期性扩展它，例如np.tile(Y, 3)，计算它的导数，并取中间部分，截断任何边界效应。

Direct computation of the 4th derivative

不是将一次导数的有限差分公式应用四次，而是将其应用于第四导数一次。仍然存在边界值的问题，我最好的想法是最简单，恒定的外推。那是：

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）。

Differentiate an interpolating spline

将5度样条拟合到数据中，然后分析得到其第四个导数。这需要SciPy并且可能比直接方法慢很多。

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

结果：

Nonuniform grid

在均匀网格上插值时，边界处的精度损失是典型的，因此我们应该期望使用切比雪夫节点进行改进。这里是：

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)

剩下的就像以前一样，使用5度的InterpolatedUnivariateSpline。绘制第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))

但是任何一种节点都存在错误：“搭配矩阵是单数”。我认为它对边界条件的处理是针对三次样条的。

Answer 2

在这里，我使用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])

因为积分步骤不是常数，所以我继续使用numpy梯度方法，因为在数值计算导数时可以考虑这种变化。

这是我用非常pythonic代码比较有和没有变量网格的结果：

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）|