For a silly screen saver I'm trying to develop, I'd like to randomly generate a divergence-free 2D array of 2D vectors, and then use it to generate a line integral convolution plot. I've heard$^1$ that one way to do this is to generate random noise, and then project out the solenoidal component of its Helmholtz-Hodge decomposition. To do that, I tried using the following reasoning:
A function $\mathbf{f}:\mathbb{R}^2\rightarrow\mathbb{R}^2$ has Helmholtz-Hodge decomposition$^2$ $$\mathbf{f}=\mathbf{h}+\nabla\phi+J\nabla\psi$$ where $$J=\pmatrix{0&-1\\1&0}$$ and where $\phi,\psi$ are scalar functions. For now, assume the harmonic component $\mathbf{h}$ vanishes.
In Fourier space, this becomes $$\mathcal{F}\mathbf{f}=-i\mathbf{k}\hat{\phi}-iJ\mathbf{k}\hat{\psi}$$ and we can define a solenoidal projection operator on Fourier space as $$P=I-\frac{\mathbf{k}\otimes\mathbf{k}}{\mathbf{k}\cdot\mathbf{k}}$$ which projects a function on its solenoidal component, via $$\mathcal{F}^{-1}P\mathcal{F}\mathbf{f}=J\nabla\psi.$$
I then tried to implement this in Mathematica, applying it to a $21\times 21\times 2$ random array. First I generate the random array, and take apply the FFT to each of its two components:
arr = RandomReal[{-1, 1}, {2, 21, 21}];
fArr = Fourier /@ arr;
I then define $\mathbf{k}$ as a function of array index:
k[k1_, k2_] := Mod[{k1 - 1, k2 - 1}, 21, -10]/21;
Then I perform the projection on the Fourier components (the singularity at $\mathbf{k}=0$ is left alone using an If
statement):
dat = Transpose[
Table[If[k1 == 1 && k2 == 1, fArr[[;; , k1, k2]],
fArr[[;; , k1, k2]] -
k[k1, k2] (k[k1, k2].fArr[[;; , k1, k2]])/(k[k1, k2].k[k1,
k2])], {k1, 21}, {k2, 21}], {2, 3, 1}];
Then I iFFT the two components:
projArr = InverseFourier /@ dat;
This gives a purely real array, and I would naively expect the result to be an approximation of $J\nabla\psi$. My question is:
- In what sense does the result approximate $J\nabla\psi$?
Supposedly the Helmholtz-Hodge decomposition of 2D data is a nontrivial task, as Chris Beaumont's HH_DECOMP routine is supposed to use FFTs to perform Helmholtz-Hodge decomposition, but he also says (in the comments at the top of the code) that the method seems inaccurate. Likewise, there are more complicated variational methods to perform Helmholtz-Hodge decompositions$^3$ of 2D data, which would seem to suggest that the simpler FFT method is somehow inadequate. Why? What does the FFT method get wrong? And is it wrong to assume the harmonic component vanishes for my random noise?
(1): Stable Fluids, Jos Stam.
(2): Feature Detection in Vector Fields Using the Helmholtz-Hodge Decomposition, Alexander Wiebel, page 12.
(3): Discrete Multiscale Vector Field Decomposition, Yiying Tong.