# Partial derivatives of a 3D array in Matlab

I'm interested in taking some partial derivatives of a 3 dimensional array in Matlab - say $A(i,j,k)$ approximates $f(x_i,y_j,z_k)$. I need to approximate things like $\partial_{xy}f$, $\partial_{yz}f$, etc. I've already built the finite difference matrices for a few of these derivatives, but ultimately this is going to be slow and inaccurate (I'd eventually like to operate on $256^3$ or even $512^3$ arrays).

Does anyone have a good suggestion for computing these partials with FFT? I've tried the following naive thing and it doesn't seem to work:

N=128;
L=2*pi;

x=L/N*(-N/2:N/2-1);y=x;z=x;
[X,Y,Z]=meshgrid(x,y,z);
f=sin(X).*sin(Y).*sin(Z);  % A 2-pi periodic function to test
fp_exact=cos(X).*sin(Y).*sin(Z); % Its exact x-partial

ik=1i*(2*pi/L)*[0:nx/2-1 0 -nx/2+1:-1]; % The spectral differentiation vector

fhat=fftn(Z);  % Compute the 3-D FFT

C=mat2cell(fhat,nx,ones(1,nx),ones(1,nx));  % So I can act on each column
result=cellfun(@(x) x.*ik',C,'UniformOutput',false);  % Multiply each column of Z by ik
fprimehat=cell2mat(result);  % Convert back to 3D array
fprime=ifftn(fprimehat);  % IFFT


The above code seems hopeless - I think I'm missing something big and haven't really thought about this hard enough. With the above definition of ik, I can compute a 1D spectral derivative with just

fprime=ifft(ik.*fft(f));


Thanks in advance for any tips.

• The usual approach I use is to generate [ikx,iky] = meshgrid(ik); (in 2D). Then ifft(ikx.*fft(f)) should get you the x derivative. The same approach should work in 3D. (I'm not sure what's wrong with your code, though).
– AJK
Dec 9, 2013 at 1:47
• I think this will actually work - so simple. Thanks! I'll post an answer tomorrow after some testing. Dec 9, 2013 at 3:06

It is far more complicated to compute derivatives with the FFT than necessary. If you just use the first order finite difference quotient, then you can approximate $\partial_x A(i,j,k)\approx \frac{A(i+i,j,k)-A(i,j,k)}{\Delta x}$ and similarly for the other derivatives. This can be done with a loop over all $i,j,k$ at a cost proportional to the number of elements, whereas the approach with the FFT is definitely more expensive.