5

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.

Eugene B
  • 995
  • 2
  • 12
  • 27
  • 1
    Why are you bothering to move from C to Python if you're just going to write C and have Python execute it (inefficiently)? – ydaetskcoR Jul 22 '14 at 11:35
  • If you're dealing with PDEs, FEniCS (http://fenicsproject.org/documentation/tutorial/index.html) might be interesting for you – Dietrich Jul 22 '14 at 11:38
  • `X`, `Y`, `Z` mean? the shape of `x`? – emesday Jul 22 '14 at 11:43
  • That is it: I would like to move to pure Python, the trick is, the current Implementation I have is the fastest I can make, though I am sure, a pure numpy solution exists. Though, you are right: with the solution I have at the moment, there is no sense in moving to Python. – Eugene B Jul 22 '14 at 11:43
  • Yeah, just the dimensions. – Eugene B Jul 22 '14 at 11:44
  • Perhaps http://stackoverflow.com/questions/16078818/calculating-gradient-with-numpy might help. I found other interesting possibilities searching the web with `numpy partial derivative` – wwii Jul 22 '14 at 12:13
  • I've searched quite a lot, but I didn't find this link. I'll have a look! Thank you! – Eugene B Jul 22 '14 at 12:15

4 Answers4

8

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.

ali_m
  • 71,714
  • 23
  • 223
  • 298
2

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
Jaime
  • 65,696
  • 17
  • 124
  • 159
2
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
1

The numpy.gradient now supports evaluating derivative along a single direction as well.

a = np.array([[1,2,3],[2,3,5]])
np.gradient(a, axis=0)

which gives a single partial derivative along axis 0:

array([[1., 1., 2.],
       [1., 1., 2.]])

The argument axis specifies a set of directions to evaluate derivative.

ASJ
  • 11
  • 1