2
$\begingroup$

I am trying to extract the curl-free component of a discrete vector field. My plan is to take the Fourier transform of the vector field and then extract the radial component in Fourier space. The inverse transform of this radial component should give me the curl-free component.

As a test I tried a vector field that has only a curl-free component, $\vec{f}(x,y)=<x,y>$. This when taken to Fourier space,extracting radial component and inverting it should give me the original field back. But I am getting very big deviations from original field.

The code I used is

(*defining the continuous function, spacing and lattice length*)
fx[x_] := x;   (*x component of vector field*)
fy[y_] := y;   (*y component of vector field*)
a = 1;         (*lattice spacing*)   
m = 5;         (*there will be 2m+1 lattice sites*)

(*reordering function, to properly reorder for inbuilt MATHEMATICA Fourier 
command*)
shifta[n_] := Mod[n + m - 1, 2*m + 1] + 1;
shiftb[n_] := Mod[n + m, 2*m + 1] + 1;

(*discretizing the function*)
frealx = Table[fx[n*a], {n, -m, m}, {l, -m, m}]
frealy = Table[fy[l*a], {n, -m, m}, {l, -m, m}]

(*reordering freal for Fourier transform*)
grealx = Table[
  frealx[[shifta[n], shifta[l]]], {n, 1, 2*m + 1}, {l, 1, 2*m + 1}]
grealy = Table[
  frealy[[shifta[n], shifta[l]]], {n, 1, 2*m + 1}, {l, 1, 2*m + 1}]

(*taking Fourier transform. Chopping to get rid of small numerical values*)
ffourx = Chop[Fourier[grealx]]
ffoury = Chop[Fourier[grealy]]

(*taking radial component*)
ffourxrad = Chop[Table[If[(-m + n - 1 ) == 0 && (-m + l - 1 ) == 0, 0,
    ((ffourx[[n, l]]*(-m + n - 1 ) + 
       ffoury[[l, n]]*(-m + l - 1 ))*(-m + n - 1  ))/((-m + n - 
       1  )^2 + (-m + l - 1 )^2)], {n, 1, 2 m + 1}, {l, 1, 2 m + 1}]]
ffouryrad = Chop[Table[If[(-m + n - 1 ) == 0 && (-m + l - 1 ) == 0, 0,
    ((ffourx[[n, l]]*(-m + n - 1  ) + 
       ffoury[[l, n]]*(-m + l - 1 ))*(-m + n - 1  ))/((-m + n - 
       1  )^2 + (-m + l - 1 )^2)], {l, 1, 2 m + 1}, {n, 1, 2 m + 1}]]

(*inverting the radial component*)
realx = InverseFourier[ffourxrad]
realy = InverseFourier[ffouryrad];

(*reordering realx and realy*)
realxreord = 
 Table[realx[[shiftb[n], shiftb[l]]], {n, 1, 2*m + 1}, {l, 1, 2*m + 1}]
realyreord = 
  Table[realy[[shiftb[n], shiftb[l]]], {n, 1, 2*m + 1}, {l, 1, 
    2*m + 1}];

So, if all were correct, then frealx == realxreord and frealy == realyreord. But these equalities don't hold. Where am I going wrong?

$\endgroup$

Your Answer

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

Browse other questions tagged or ask your own question.