ifft2 not returning expected plots

2 views (last 30 days)
William Kett
William Kett on 3 Jun 2021
To start with, I am trying to solve the Poisson Equation using fft2 and ifft2. I am starting with the solution, with a=1, taking the Laplacian of u(x,y) to get and then trying to use the fft2 and ifft2 functions to work back to u(x,y). A more in-depth description of the problem and process can be found in the top answer found here. The problem comes at the 4th step in the linked answer when I take ifft2. The shape of the plot that I obtain for u(x,y) is wrong and the values are way too large. I will post my code below so that anyone can reproduce the problem as well as the plots of what I have obtained in the steps leading up to the problem.
syms u(x,y,a)
u(x,y,a) = 2^(4*a) * x^a * (1-x)^a * y^a * (1-y)^a;
u(x,y,a) = subs(u(x,y,a),a,1);
X = 1;Y=1;
Lx = 10; Ly = 10;
dx = X/Lx; dy = Y/Ly;
x_arr = meshgrid(0:dx:X,0:dx:X);
y_arr = transpose(meshgrid(0:dy:Y,0:dy:Y));
u_func(x,y,a) = u(x,y,a);
u(x,y,a) = subs(u(x,y,a),x,x_arr(1,:));
u(x,y,a) = subs(u(x,y,a),y,y_arr(:,1));
u_xy = double(u(x,y,a));
figure(1)
contour3(x_arr,y_arr,u(x,y,a));
surface(x_arr,y_arr,u(x,y,a));
title('u(x,y)');
xlabel('x = linspace(0,1,11)');ylabel('y = linspace(0,1,11)');
f_func = diff(u_func,x,2) + diff(u_func,y,2);
f_func = subs(f_func,x,x_arr(1,:));
f_func = subs(f_func,y,y_arr(:,1));
M = 10;
N = 10;
%m = meshgrid(linspace(0,M,M+1),linspace(0,N,N+1));
%n = transpose(meshgrid(linspace(0,N,N+1),linspace(0,N,N+1)));
P = M; Q = N;
Lp = P; Lq = Q;
dp = Lp/P;dq = Lq/Q;
p = meshgrid(0:dp:P,0:dp:P);
q = transpose(meshgrid(0:dq:Q,0:dq:Q));
f = double(f_func);
fhat = fft2(f,P+1,Q+1); %m = n = 10
km_axis = (-1j * 2 * pi/length(x_arr)).*p;
ln_axis = transpose((-1j * 2 * pi/length(y_arr)).*q);
figure(2)
contour3(real(km_axis) + imag(km_axis),real(ln_axis) + imag(ln_axis),abs(fftshift(fhat)));
surface(real(km_axis) + imag(km_axis),real(ln_axis) + imag(ln_axis),abs(fftshift(fhat)));
title('f(x,y) = $\nabla^2 u(x,y)$','Interpreter','LaTex')
xlabel('$k_{m}$ = $-j2 \pi p/Lx$','Interpreter','LaTex');ylabel('$l_{n}$ = $-j2 \pi q/Ly$','Interpreter','Latex');
uhat = -fhat./ (km_axis.^2 + ln_axis.^2);
mask = isinf(uhat);
uhat(mask) = 0;%assume that the values are 0 along the boundary
u = ifft2(uhat);
figure(3)
contour3(x_arr,y_arr,abs(u));
surface(x_arr,y_arr,abs(u));
title('$u(x,y) = ifft2\frac{\hat{f}}{(k_m^2 + l_n^2)}$','Interpreter','LaTex')
xlabel('x');ylabel('y');
I'm expecting my plot for figure 3 to match my plot for figure 1 in both shape and values. What I keep getting instead is a 3D parabola with values which are way too large and I am not certain as to why this is the case. I have been hung up on this for at least a month now and don't understand what I could be doing wrong whether it be in my methodology or if I am making a mistake somewhere in the code.
One thing I have noted is that if I change p and q to be centered at 0 such that:
p = meshgrid(-P/2:dp:P/2,-P/2:dp:P/2);
q = transpose(meshgrid(-Q/2:dq:Q/2,-Q/2:dq:Q/2));
I get a result for figure 3 that is much closer to figure 1 though I do not know why this is the case. Regardless, I cannot seem to figure out what is wrong with my approach whether it be in the methodology or in my usage of the functions.

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!