import numpy as np
from scipy.sparse import csr_matrix as sparse_matrix
# scipy.sparse.csr_matrix
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csr_matrix.html
def make_differential_ops(nx, dx):
# operators (matrix) that shift vector components (periodic boundary condition)
f0 = np.identity(nx, dtype=int) # f_{i}
f1 = np.roll(f0, 1, axis=1) # f_{i+1}
f2 = np.roll(f0, 2, axis=1) # f_{i+2}
f_1 = f1.transpose() # f_{i-1}
f_2 = f2.transpose() # f_{i-2}
# (f_{i+1} - f_{i-1}) / (2 dx)
deriv1 = sparse_matrix(f1 - f_1) / (2.0 * dx)
# (f_{i+1} - 2f_{i} + f_{i-1}) / (dx^2)
deriv2 = sparse_matrix(f1 - 2.0 * f0 + f_1) / dx**2
# (f_{i+2} - 2f_{i+1} + 2f_{i-1} - f_{i-2}) / (2 dx^3)
deriv3 = sparse_matrix(f2 - 2.0 * f1 + 2.0 * f_1 - f_2) / (2.0 * dx**3)
return deriv1, deriv2, deriv3