Generating a set of divergence-free basis

I'm trying to derive an algorithm that would generate a set of divergence-free vectors, which shall be used as basis vectors later on. Using a simple example, a 2D second-order divergence-free basis would look like this:

$$\begin{bmatrix} 1 \\ 0 \end{bmatrix}, \begin{bmatrix} 0 \\ 1 \end{bmatrix}, \begin{bmatrix} x \\ -y \end{bmatrix}, \begin{bmatrix} y \\ 0 \end{bmatrix}, \begin{bmatrix} 0 \\ x \end{bmatrix}, \begin{bmatrix} x^2 \\ -2xy \end{bmatrix}, \begin{bmatrix} xy \\ -\frac{1}{2} y^2 \end{bmatrix}, \begin{bmatrix} y^2 \\ 0 \end{bmatrix}, \begin{bmatrix} 0 \\ x^2 \end{bmatrix}.$$

So far, I have manually implemented these calculations, but it will be quite tedious to work it out for 3D cases and for higher order. As of now, my pseudocode for 2D case looks like this:

switch (order):
order=1:
// manually calculate vectors for first order
order=2:
// manually calculate vectors for second order
[...]

I would imagine things to be nastier in 3D... Currently, I am trying both to derive the whole basis manually and to come up with an algorithm to do it. Any suggestions?

• math.stackexchange.com or physics.stackexchange.com might be a better place to ask such a question. – Jakube Aug 1 at 20:44