Return the second order central finite difference of a vector

second_order_fda.py

import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0, 50, 3001)
y = np.sin(3*x)
delta_x = x[1]
def second_order_central_fd(y, delta_x):
    """
    Return the second order central finite difference of a vector (i.e. d^2 y/ dx^2).
    Zero fills rollover. What it should probably do is a single order FDA
    or interpolate the endpoints.
    Parameters
    ----------
    y : np.array
        Vector of function values who's second derivative is to be calculated.
    delta_x : float
        Spacing in x between nodes.
    Returns
    -------
    np.array
    
    Notes
    -----
    Was surprised I couldn't find a vectorized implementation of this out there so figured I'd write one.
    
    """
    minus_1 = np.roll(y, 1)
    minus_1[0] = 0
    plus_1 = np.roll(y, -1)
    plus_1[-1] = 0
    ret = (minus_1 + plus_1 - 2*y) / delta_x**2
    return ret
y_guess = second_order_central_fd(y, delta_x)
plt.plot(x[1:-1], y_guess[1:-1])
# Exact
plt.plot(x, -9*np.sin(3*x))
plt.show()