There are two issues that you're running into. One is mathematical and the other

is numerical.

Here's a version of the narrative that skips a lot of the Finite Element Method

behind it, but will ultimately lead you to a good solution.

The best way to understand polyharmonic weights is to consider the energies

they minimize. In the smooth setting the correspondence goes like this:

∫ ‖∇u‖² dA <--> ∆u = 0

∫ (∆u)² dA <--> ∆²u = 0

∫ ‖∇∆u‖² dA <--> ∆³u = 0

...

In the discrete setting, our primary building blocks will be two matrices: L and

M. These matrices discretize the following energies:

∫ uv dA ~~~> u' M v

∫ ∇u⋅∇v dA ~~~> u' L v

L is the "usual" symmetric cotangent matrix. And M is the mass matrix (often a

diagonal matrix of vertex areas). Often we think of L as the Laplacian, but

notice that here it's acting as integrating the product of the gradients of the

two functions u and y. We recognize it as the Laplacian because Green's identity

tells us this equality:

-∫ u∆v dA + boundary terms = ∫ ∇u⋅∇v dA ~~~> u' L v

So now we can think of L as computing *integrated* Laplacian of y times u.

Skipping the FEM behind it, we can "unintegrate" this by applying M⁻¹ so that:

∆y ~~~> M^{-1}L y

Notice the integral is gone and here we're dealing with pointwise quantities.

Now we have enough to go back to the original energies and discretize each of

them:

∫ ‖∇u‖² dA ~~~> u'Lu

∫ (∆u)² dA ~~~> (M^{-1}L u)' M (M^{-1}L) u = u'LM^{-1}L u

∫ ‖∇∆u‖² dA ~~~> (M^{-1}L u)' L (M^{-1}L) u = u' L M^{-1}L M^{-1}L u

Minimize these is as easy as stripping off the u' and setting the equation to 0

on the right hand side.

However... multipying all these L's and M^{-1}s will start to create a really

poorly conditioned matrix. I've also observed that for triharmonic and higher

the results of using a direct solver like Cholesky can be poor. So a better

solution is to solve the equivalent but better conditioned problems. For example

for the biharmonic equation we should minimize:

min u'LM^{-1}L u

u

we can introduce a auxiliary vector of variables v and solve the constrained

problem:

min y M y subject to Lu = My

u,y

Using the Lagrange multiplier method, introduce another vector of variables λ

and solve the unconstrained saddle problem:

saddle y M y + λ(Lu - My)

u,y,λ

taking derivatives to u,y, and λ and setting them to zero reveals the solution:

/ 0 0 L \ / u \ / 0 \

| 0 M -M | | y | = | 0 |

\ L -M 0 / \ λ / \ 0 /

Notice that the middle equation just says that y=λ, so this can be factored out

leaving:

/ 0 L \ / u \ / 0 \

\ L -M / \ y / = \ 0 /

Now, you could further factor out y = M^{-1}L u, and this would lead to the

LM^{-1}L system above. But the point is to keep the system matrix bigger (and

sparser) and better conditioned.

Solving this "unflattened" system will (in my experience) produce better

results, purely due to better numerics. Algebraically the solutions are the

same.

You can similarly work this out for triharmonic and so on.

This is skipping the FEM theory which you could find in my thesis among other

places. Hopefully this helps.

-Alec