Error in numerical derivative of 4th order near the endpoints of an interval

Question

I am integrating a partial differential equation in which I have a fourth-order partial derivative in $x$ and the numerical integration in time is giving me absurd errors. The cause of the problem, I believe, is that I get large errors in the fourth-order derivative. To illustrate that I take numerical derivatives of the function $y(x)=1-\cos(2\pi x)$. Below I plot the derivatives $y_{xx}(x)$ and $y_{xxxx}(x)$ in the domain $(0, 1.0)$. See the figures below:

As one can see the errors occur mostly near the boundaries.

The numerical derivatives were performed using the numpy gradient method. The python code is here:

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()

My question is how do I reduce algorithmically this large error at boundaries beyond the obvious increasing of N? As it is the convergence of the fourth-order derivative is not uniform. Is it possible to make that uniform? Perhaps, using extrapolation at the boundaries will do.

Answer 1

Underlying mathematical issue

The formula for the 1st derivative that np.gradient uses has error O(h**2) (where h is your dx). As you keep on taking derivatives, this has potential to be bad because changing a function by an amount z can change its derivative by z/h. Usually this does not happen because the error comes from higher order derivatives of the function, which themselves vary smoothly; so, subsequent differentiation "differentiates away" most of the error of the preceding step, rather than magnifying it by 1/h.

However, the story is different near the boundary, where we must switch from one finite-difference formula (centered difference) to another. The boundary formula also has O(h**2) error but it's a different O(h**2). And now we do have a problem with subsequent differentiation steps, each of which can contribute a factor of 1/h. The worst-case scenario for the 4th derivative is O(h**2) * (1/h**3) which fortunately does not materialize here, but you still get a pretty bad boundary effect. I offer two different solutions which perform about the same (the second is marginally better but more expensive). But first, let's emphasize an implication of what Warren Weckesser said in a comment:

If the function is periodic, you extend its by periodicity, e.g., np.tile(Y, 3), compute the derivative of that, and take the middle part, truncating away any boundary effects.

Direct computation of the 4th derivative

Instead of applying a finite-difference formula for the 1st derivative four times, apply it for the 4th derivative once. There is still going to be the issue of the boundary values, for which my best idea is the simplest, constant extrapolation. That is:

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]]))

with the result

If you don't like the little flat bits on the side, take smaller dx (they are of size 2*dx).

Differentiate an interpolating spline

Fit a 5th degree spline to the data, then take its 4th derivative analytically. This requires SciPy and is probably a lot slower than the direct approach.

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

Result:

Nonuniform grid

A loss of precision at the boundaries is typical for interpolation on a uniform grid, so we should expect an improvement using Chebyshev nodes instead. Here it is:

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)

The rest proceeds as before, using InterpolatedUnivariateSpline of degree 5. Plotting the 4th derivative itself now shows no visible difference between numerics and analytics, so here is the plot of Y_xxxx - analytics instead:

The error is at most 3, and it's no longer maximal at the boundary. The maximal error with uniform grid was about 33.

I also explored the possibility of imposing the clamped condition on the interpolating spline to further improve its accuracy. Conceivably, make_interp_spline could do this with

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

but there are errors with either kind of nodes: "collocation matrix is singular". I think its handling of boundary conditions is targeted toward cubic splines.

Answer 2

Here I use extrapolation at the boundaries with numpy polyfit method.

# 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])

Furthermore, I have a smaller mesh at the boundaries so as to decrease the error of the numerical derivative (Y_xxxx) there, without having to increase the number meshpoints uniformly accross the domain of integration.

# 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])

Because the integration step is not constant, I keep using the numpy gradient method, since it can take that variation into account when calculating the derivatives numerically.

Here's my not very pythonic code to compare results with and without variable mesh:

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()

Here's the comparison with two methods (both with extrapolation at the boundaries), One of them with variable steps.

Here Error(x)= log_10|Y_xxxx(x)+(2*pi)**4*cos(kx)|