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.

Improve this question
$\endgroup$
2
  • 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, 2013 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, 2013 at 3:06

1 Answer 1

Reset to default
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.

Improve this answer
$\endgroup$
2
  • $\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, 2013 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, 2013 at 18:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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