2
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ 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). $\endgroup$ – AJK Dec 9 '13 at 1:47
  • $\begingroup$ I think this will actually work - so simple. Thanks! I'll post an answer tomorrow after some testing. $\endgroup$ – icurays1 Dec 9 '13 at 3:06
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ I just had this realization this morning - I'm trying both approaches now. The FFT time is negligible for me (GPU) so the accuracy gain might be worth it. $\endgroup$ – icurays1 Dec 9 '13 at 14:01
  • $\begingroup$ You could of course use higher order finite difference quotients instead of the first order one-sided one I've used above. $\endgroup$ – Wolfgang Bangerth Dec 9 '13 at 18:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.