I am trying to develop central finite difference scheme (lexicographical) coefficients. My code is,
import numpy as np from scipy.sparse import spdiags, identity x = np.linspace(0, 9, num=10) nx = len(x) ` hx = (x[-1] - x[0])/(nx - 1) one_x = np.ones((nx, 1)) x_coeffs = np.hstack((one_x, one_x)) x_coeffs = np.insert(x_coeffs, 1, -2, axis=1) # So far so good, then carry on x_diff = spdiags(x_coeffs.transpose(), np.asarray([0, 1, 2]), nx-2, nx) nx_id = identity(nx) Dx = np.kron(nx_id, x_diff.transpose()) Dx = 1/hx*DxLet us try it on a dummy example.
x = np.linspace(0, 4, num=5) f = np.power(x, 2) dfdx = Dx * fBut the result for this dummy example is (array([ -4., -6., -8., -10.]). Any idea about what could go wrong?