I am slowly moving from C to Python. This time I need to calculate partial derivatives numerically from a grid given. I know how to do it in C, so at the moment I just use inline adapter, i.e.
def dz(x,X,Y,Z,dx):
y = numpy.zeros((X,Y,Z), dtype='double');
code = """
int i, j, k;
for (i=0; i<X-1; i++){
for(k=0; k<Y; k++){
for (j=0; j<Z; j++){
y[i,k,j] = (x[i+1, k, j] - x[i, k, j])/dx;
}
}
}
for (j=0; j<Z; j++){
for(k=0; k<Y; k++){
y[X-1,k,j] = - x[X-1, k, j]/dx;
}
}
"""
weave.inline(code, ['x', 'y', 'dx', 'X', 'Y', 'Z'], \
type_converters=converters.blitz, compiler = 'gcc');
return y;
where x
and y
are 3D numpy arrays, as you can see, and the second loop stands for boundary conditions. Of course, I can implement the same logic in pure Python, but the code would be inefficient. I wonder, though, if it is possible to calculate a partial derivative using pure numpy? I would appreciate any help anyone can provide.
np.diff
might be the most idiomatic numpy way to do this:
y = np.empty_like(x)
y[:-1] = np.diff(x, axis=0) / dx
y[-1] = -x[-1] / dx
You may also be interested in np.gradient
, although this function takes the gradient over all dimensions of the input array rather than a single one.
If you are using numpy, this should do the same as your code above:
y = np.empty_like(x)
y[:-1] = (x[1:] - x[:-1]) / dx
y[-1] = -x[-1] / dx
To get the same result over the second axis, you would do:
y = np.empty_like(x)
y[:, :-1] = (x[:, 1:] - x[:, :-1]) / dx
y[:, -1] = -x[:, -1] / dx
def dz(x,dx):
y = numpy.zeros(x.shape, dtype='double')
y[:-1] = (x[1:] - x[:-1]) / dx
y[-1] = -x[-1] / dx
return y
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With