Finding a scalar field in order to generate a solenoidal vector from a given vector

I am working on generating a (complex) solenoidal vector field $\mathbf{A}$ from a prescribed (complex) vector $\mathbf{a}$ and the gradient of a scalar, $b$, such that $$\mathbf{A} = \mathbf{a} + \nabla b.$$ The condition that I require is then $$\nabla\cdot\mathbf{A} = \nabla\cdot\mathbf{a} + \nabla^2 b = 0 \qquad \Longrightarrow \qquad \nabla^2 b = - B,$$ where $B = \nabla\cdot\mathbf{a}$.

To complicate this slightly, I want to decompose $\mathbf{A},\mathbf{a}, B$ and $b$ into $$\mathbf{A} = \hat{\mathbf{A}} \exp[i(\omega t - \mathbf{k}\cdot\mathbf{x})] = \hat{\mathbf{A}} \exp[i\phi],\\ \mathbf{a} = \hat{\mathbf{a}} \exp[i\phi],\\ B = \hat{B} \exp[i\phi],\\ b = \hat{b} \exp[i\phi].$$

With this decomposition, $B = (\nabla\cdot\hat{\mathbf{a}} - i\mathbf{k}\cdot\hat{\mathbf{a}}) \exp[i\phi]$, and the Poisson equation above becomes $$(\dagger)\qquad\nabla^2 b = (\nabla^2 \hat{b} - 2i\mathbf{k}\cdot\nabla b - \lvert\mathbf{k}\rvert^2 b) \exp[i\phi] = -\hat{B} \exp[i\phi].$$

This is then solved for $\hat{b}$ iteratively using the following solver scheme (given for 1-dimension, but it is easily extended to higher dimensions): $$(\ddagger)\qquad \hat{b}_j = \frac{\Delta x^2 \hat{B}_j + \hat{b}_{j+1} + \hat{b}_{j-1} - ik\Delta x(\hat{b}_{j+1} - \hat{b}_{j-1})}{2 + k^2 \Delta x^2},$$ where $j \in \{2,\ldots,N-1\}$ indicates the discrete grid location, the position $x_j = j \Delta x$, and $\Delta x$ is the grid spacing. The boundary $j=1,N$ is fixed at zero. $\hat{B}_j$ is evaluated as $\hat{B}_j = (\hat{a}_{j+1} - \hat{a}_{j-1})/2\Delta x -ik\hat{a}_j$.

The numerical solution to $\hat{b}$ produced satisfies $$\nabla^2 \hat{b} - 2i\mathbf{k}\cdot\nabla b - \lvert\mathbf{k}\rvert^2 b = -\hat{B}$$ very well, however if I construct $$\hat{A}_j = \hat{a}_j + \frac{\hat{b}_{j+1} - \hat{b}_{j-1}}{2\Delta{x}} - ik\hat{b}_j$$ and evaluate $$\nabla \cdot \mathbf{A} \simeq \Big(\frac{\hat{A}_{j+1} - \hat{A}_{j-1}}{2\Delta x} - ik \hat{A}_j\Big) \exp[i\phi],$$ then this is not only non-zero but also large. Putting the expression for $\hat{A}_j$ in terms of $\hat{a}$ and $\hat{b}$, one finds that the above equation reads $$(*)\qquad \nabla \cdot \mathbf{A} \simeq \Big( \frac{\hat{a}_{j+1} - \hat{a}_{j-1}}{2\Delta x} - ik \hat{a}_j + \frac{\hat{b}_{j+2} - 2\hat{b}_{j} + \hat{b}_{j-2}}{4\Delta x^2} - 2ik\frac{\hat{b}_{j+1} - \hat{b}_{j-1}}{2\Delta x} - k^2\hat{b}_j \Big) \exp[i\phi]\\ = \Big( \hat{B} + \nabla_2^2\hat{b} - 2ik \frac{\partial\hat{b}}{\partial x} - 4k^2 \hat{b} \Big) \exp[i\phi].$$ This is just the 1-dimensional version of $(\dagger)$, but with the Laplacian operator evaluated as $$\nabla_2^2 \hat{b} = \frac{\hat{b}_{j+2} - 2\hat{b}_j + \hat{b}_{j-2}}{4\Delta x^2},$$ suggesting that it is this difference which is responsible for $\nabla \cdot \mathbf{A} \neq 0$.

To overcome this, I attempted to invert ($*$) rather than use ($\ddagger$), this time with boundary conditions $\hat{b}_j = 0$ for $j = 1, 2, N-1, N$, but the (iterative) solution to the Poisson equation diverges.

I was wondering what is responsible for this? Could my setting the two boundary cells to zero cause such a problem? The solver scheme for point $j$ involves $j+2, j+1, j-1, j-2$. Given the non-solenoidal vector $\mathbf{a}$, is there a better way to find a scalar $b$ such that $\nabla\cdot(\mathbf{a} + \nabla b) = 0$?

Thank you for any help or suggestions!

• Is $k$ a fixed vector here? That doesn't seem to make any sense. Shouldn't you have a summation over $k$? – David Ketcheson Dec 13 '14 at 4:18
• Thanks for your comment, @DavidKetcheson. I should have said, I am decomposing the quanities into a carrier phase with angular frequency $\omega$ and wave vector $\mathbf{k}$, and a modulation or envelope ( given by the hatted fields). Hope that helps. – samwise Dec 16 '14 at 18:47