Hi all,

I am trying to solve a biharmonic equation, $\Delta^2 w = f$ via the SparseLU solver.

Unfortunately, the system is underdetermined (which is clear from the equation, I can add any constant to a solution) so in order to get a unique solution, I have an additional constraint, i.e. the integral over the entire domain of $w$ has to be zero.

The solver seems to be fine with solving the underdetermined system, but how can I interpret the result? Is this a special solution of the entire solution manifold, or just any of those?

When I add the integral constraint, I get a set of n+1 equations for n variables, and the solver dies with a complaint that it wants quadratic matrices. Is there a common trick to blow up the matrix to be (n+1)x(n+1) or how are similar problems solved usually?

thanks a lot
Hecke

To add hard constraints, you can use Lagrange multipliers. I have tutorial notes on that matter, see the http://giana.mmci.uni-saarland.de/EG14Course/, in particular the "Numerical Linear Algebra" slides address this issue from a practical point of view.

Hi Gael,

thanks a lot for your reply.
And sorry for the stupid thread title, I don;t know how this happened...

I read your slides, you seem to be really deep in the field, great!
And I was surprised that you work with Ivo, can you pass him my greetings? We don't know each other in person, I just worked really closely with his brother for some years

To the problem:
Is Lagrange multipliers the best (easiest, most robust...) way to introduce constraints?
I am in contact with one of the deal.ii guys in that matter, and he proposed to just fix some value and adjust the result as I already know the null space of the problem, but I have doubts that the methods derived from finite element shape functions are applicable in my case of finite differences and an integral constraint. Can you comment on that?

In your slides you say that I should not rely on the solvers methods to choose a solution from the multiple ones that exist. What methods does the sparseLU solver of Eigen use to determine the solution he returns?

thanks again
Hecke

small world

If some values are fixed, then re-order the unknowns and move the fixed ones to the right hand side. For an integral constraint you cannot do that, so Lagrange multiplier should do the job. FEM or finite-difference, no difference here.

If the problem is still undetermined, then SparseLU will pick an arbitrary solution, so better add constraints to get the one you want.

So, there is no way to say which solution sparseLU will take?

best
Hecke

Hi Gael,

just another short question:
If I know about the structure of my matrix, i.e. it is rank deficient by one, and I have the constraint equation, wouldn't it be a straightforward solution to append the constraint, get a (n+1)xn matrix and use the sparseQR solver?
As I have then a full rank of n, the least squares solution would be exact, right?

The approach with a Lagrange multiplicator seems well suited, but means a much higher coding effort for me...

thanks a lot
Hecke

Adding one row seems to be fine. BTW, is your initial problem symmetric?

cool. Than that is the way to go.

Currently my problem is symmetric, as I use a simple biharmonic equation with Neumann bc, but this is just the warm up phase.
Later I want to simulate the Karman-von Donnell Equations which contain the biharmonic as a subpart and which are no longer symmetric.

Another hint to me is to use Normal equations like described here on slide 7:
http://cims.nyu.edu/~donev/Teaching/Sci ... andout.pdf
As I know that my problem (even though not square) is full rank, I might end up with a Choleski decomposition, this would be even more attractive, right?

thanks and best wishes
Hecke

yes, I wanted to suggest you the normal equation as for sparse problems, Cholesky is really much faster, but on the other hand the matrix A^T * A will be much more dense, and since you start from a biharmonic equation which is already not extremely sparse, that might not be a good idea. Anyway, with Eigen it's just a matter of two lines of code to try it!

cool, then I'll give it a try and compare the methods for my case.

thanks a lot
Hecke