0

我正在尝试在 3 维中使用 Matlab 的 ifftn 来获得物理空间的解决方案。特别是我试图在 1/k^2 上使用 ifftn。物理空间的解析解是 1/(4*pi*r)。但是我没有恢复。请注意:$r = sqrt(x^2 + y^2 + z^2)$ 和 $k = sqrt(kx^2 + ky^2 + kz^2)$。

clc; clear;

n = 128; % no. of points for ifft
L = 2*pi; % size of the periodic domain
x = linspace(-L/2,L/2,n); y = x; z = x; % creating vectors for x-y-z direction

[X,Y,Z] = meshgrid(x,y,z); % creating meshgrid for physical space
R = sqrt(X.^2 + Y.^2 + Z.^2); % use for 1/(4*pi*r)


k = (2*pi/L)*[0:n/2 -n/2+1:-1]; % wave vector;
[k1,k2,k3] = meshgrid(k,k,k);
denom = (k1.^2 + k2.^2 + k3.^2); % This is k^2

F = 1./denom; F(1,1,1) = 0; % The first value is set to zero as it is infinite
phi = 1./(4*pi*R); % physical domain solution 

phys = fftshift(ifftn(F)); % Using ifftn
ph_abs = abs(phys); 

mid = ph_abs(n/2,:,:); % looking at the midplane of the output
mid = permute(mid,[3 2 1]); % permuting for contourplot.


PHI = phi(n/2,:,:); %looking at the midplane of the physical space.
PHI = permute(PHI,[3 2 1]);

figure(1)
surf(x,z,log(mid))
shading flat
colorbar();

figure(2)
surf(x,z,log10(abs(PHI)))
shading flat
colorbar();
4

0 回答 0