How to calculate boundary points of numerical derivative in python?

Question

I'm trying to write a function to take the derivative of any general function / array of numbers. Specifically, I am using a Central difference formula. The issue is, I cannot compute the boundary points of the derivative as the central difference formula uses indices that would be out of bound. My code is below

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

When I plot y,y1,y2,y3 you can see it blows up at the end points

What I have tried to do is set the end points to their nearest neighbors in deriv as below. While this works for low order derivatives (1st and 2nd) it starts to break at higher order derivatives (3rd and greater).

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

The derivative in the middle, away from the boundaries, is calculating fine. The issue is with the boundaries.

How should I treat the boundary conditions here? Would a different numerical differentiation scheme work better than the central difference scheme?

EDIT: I am looking for a general method to solve this, not just a method that can be applied to the sine function or any other periodic function as I have used to illustrate the issue here.

Answer 1

This is more a numerical methods question than a programming question.

Anyways, if your function has periodic boundary conditions (it looks it is a sinusoidal wave so in this case you have periodicity) just create a new array with 2 additional elements: the new array start element will be your last element of the original array and the end element of the new array will be the start element of the original array. Here is a way to do it

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

You can now differentiate on f_periodic for which d[1] and d[-2] will be the correct derivative value on the boundaries (disregard d[0] and d[-1]) .

Edit after OP's new requirements...

For more general boundary conditions, say a specific value at the boundaries, there are different approaches one can follow:

1. Use ghost values:

Again extend the function and extrapolate values for the new boundaries. Depending on the order of numerical differentiation more ghost cells will be required. For the current scheme, a simple linear extrapolation will do (only 1 ghost value at each boundary are required):

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

Note that you have to also extend x. However, since you have a constant spacing just use h instead of spatial differences, e.g. x[-2]-x[-3]. You can now differentiate f_new and you will get an 1st-order approximation of the derivative on the boundaries (since you used a linear extrapolation to find the ghost value).

2. Use forwards and backwards schemes on the boundaries

I will not show code here, but basically you need to differentiate using an boundary value and the right (forwards) or left (backwards) value for the left and right boundaries respectively. This is a first-order approximation.