# How to calculate a numerical Fourier transform obtained from NDSolve?

We numerically solve for the electric potential on a box which obeys Ohm's law: $$$$\nabla\cdot(\overset{\scriptscriptstyle\leftrightarrow}{\sigma} \nabla \Phi) = 0, \label{Eqn:OhmsLaw}$$$$ where $$\overset{\scriptscriptstyle\leftrightarrow}{\sigma}$$ is the 2D conductivity tensor, and $$\Phi$$ are the electric potential and 2D current density respectively. The diagonal and off-diagonal elements of the conductivity tensor are given by the functions defined in the Mathematica code. They look like this:

The first plot is $$\sigma_{xx}$$ as a function of $$r$$ and second plot is $$\sigma_{xy}$$ and everything is radially symmetric. When $$\sigma_{xx}$$ is zero in that tiny interval the electric potential $$\hat{z}\cdot((\nabla\sigma_{xy})\times(\nabla\Phi))$$, contours of electric potential are parallel to contour of constant $$\sigma_{xy}$$ and since it is radially symmetric they should be circles.

Now here is the Mathematica code :

omegactimestau[B_] := 4.8 10^-10 B/(1.7 10^-29 3 10^10) 3 10^-4/10^8;

x1[r_] := -12.000000000000016 + 20.000000000000025 r;

(*Define conductivity tensor*)
sigmaxy2[x_, y_, B_] :=Piecewise[{{-1., Inequality[0.45, LessEqual, r, Less, 0.55]}, {0., r < 0},
{1., Inequality[0.6499999999999999, LessEqual, r, Less, 0.75]},
{-12.000000000000016 + 20.000000000000025*r, Inequality[0.55, LessEqual, r, Less,
0.6499999999999999]}, {18.849565896451544 - 32.1787625422147/Sqrt[1 + 4.*r^2], r >= 0.75}},
22.918285831105766 - 32.17876254221469/Sqrt[1 + 4.*r^2]]

sigmaxx2[x_, y_, B_] := 1/omegactimestau[B] (10^-3 + Piecewise[{
{Sqrt[x1[r]] (Sqrt[1 - x1[r]]), 0.6 <= r <= 0.65},
{-Sqrt[x1[r]] (Sqrt[-1 - x1[r]]), 0.55 <= r < 0.6},
{0, 0.45 <= r < 0.55},
{(37.69913179290309 - 64.3575250844294/Sqrt[1 + 4. r^2] - 2),
r > 0.75},
{Abs[
45.83657166221153 - 64.35752508442938/Sqrt[1 + 4. r^2]] -
2, 0 <= r < 0.45}
}] )/. {r -> Sqrt[(x - Lx/2)^2 + (y - Ly/2)^2]};

(*Solve Ohms law*)

gauss = 1 ;
tesla = 10^4 gauss ;
Lx = 15; Ly = 15;

rectangle = ImplicitRegion[0 <= x <= Lx && 0 <= y <= Ly, {x, y}];
frefine = Function[{vertices, area}, Block[{x1, y1}, {x1, y1} = Mean[vertices];
If[(x1 - Lx/2)^2 + (y1 - Ly/2)^2 < 100 rad^2 && area > 0.001, True,
False]]];
mesh = ToElementMesh[rectangle, MaxBoundaryCellMeasure -> 10^-2,
MeshRefinementFunction -> frefine];

<< NDSolveFEM

(*eqn[B_]=Inactive[Div][{{sigmaxx[x,y,B],sigmaxy[x,y,B]},{-sigmaxy[x,y,B],\
reg = mesh;

brange = Range[0.5, 2, 0.5]*tesla;
sol1 = NDSolveValue[{Inactive[
Div][{{sigmaxx2[x, y, #], sigmaxy2[x, y, #]}, {-sigmaxy2[x, y, #],
sigmaxx2[x, y, #]}} . Inactive[Grad][phi[x, y], {x, y}], {x, y}] ==
0, phi[Lx, y] == 5, phi[0, y] == -5}, phi, {x, y} \[Element] reg] & /@
brange;

(*Plots*)

ContourPlot[
sol1[x + Lx/2, y + Ly/2] /. sol1 -> sol1[[2]], {x, -1, 1}, {y, -1, 1},
AspectRatio -> Automatic, ColorFunction -> ColorData["SunsetColors"],
ImageSize -> 650, PlotRangePadding -> None, FrameStyle -> Black,
BaseStyle -> FontSize -> 22, PerformanceGoal -> "Quality", Contours -> 100,
PlotRange -> All, PlotPoints -> 100]


Having obtained a solution for phi[x, y] ($$\equiv \Phi(x,y)$$) I want to calculate $$\bar{\Phi}$$ defined as $$$$\bar{\Phi} ({r})= \int\,\frac{d^2 q}{2\pi^2}\tilde{\Phi}({q}) J_{0}\left(R_{c}|{q}|\right)e^{i {q}\cdot{r}}, \label{convphi}$$$$ where $$\tilde{\Phi}(q)$$ is the Fourier transform of $$\Phi(x,y)$$, $$J_{0}$$ is the zero order Bessel function of first kind, $$R_{c}$$ is some constant parameter (for the code above it is of the order 0.1). Essentially the steps are first take inverse Fourier transform of solution obtained from NDSolve and multiply by the Bessel function and take Fourier transform of result.

In Mathematica there is no 2d Fourier transform equivalent for NFourierTransform.

Is there some alternative for the 2d case?

Edit : This question is edited after being deleted by me because I realized that my previous question had no good answer or that my question was bad. I apologize to the stack exchange and the answer below for bad practice.

• Your code as posted returns an empty plot on my system (MM 13.1 on ARM Mac.) When I run it as is, it complains about how "{x,y} $\in$ GlobalToElementMesh[...] is not a valid region specification; I think this is because you invoke ToElementMesh before you load NDSolveFEM. Fixing that problem leads to another error "The PDE coefficient (big mess containing Piecewise functions) does not evaluate to a numeric matrix", which I'm not sure how to fix. Commented Jun 28, 2023 at 15:00
• That said: if you want to calculate the 2D Fourier transform then can't you just use NIntegrate to evaluate$$\frac{1}{2\pi} \iint \Phi(x,y) e^{i(q_x x + q_y y)}\, dx\, dy?$$ Commented Jun 28, 2023 at 15:03

To long for a comment:

Without Inactive Mathematica 12.2 evaluates

sol1 = NDSolveValue[{
Div [{{sigmaxx2[x, y, #],
sigmaxy2[x, y, #]}, {-sigmaxy2[x, y, #],
sigmaxx2[x, y, #]}} . Grad [phi[x, y], {x, y}], {x, y}] ==
0, phi[Lx, y] == 5, phi[0, y] == -5},
phi, {x, y} \[Element] reg] & /@ brange

ContourPlot[sol1[[2]][x + Lx/2, y + Ly/2] , {x, -1, 1}, {y, -1, 1},
AspectRatio -> Automatic, ColorFunction -> ColorData["SunsetColors"],
ImageSize -> 650, PlotRangePadding -> None, FrameStyle -> Black,
BaseStyle -> FontSize -> 22, PerformanceGoal -> "Quality",
Contours -> 100, PlotRange -> All, PlotPoints -> 100]
`