I would like to solve the biharmonic equation \Delta^2 u=f(x,y) on [-1,1]\times[-1,1] with u=\frac{\partial u}{\partial n}=0 on the boundary. I applied the weak form and then used Legendre-Galerkin spectral method. The Dirichlet condition u=0 can easily be incorporated but what about the other essential Neumann boundary condition? Can I use the method of Lagrange-multipliers for Neumann-BCs, that is append the system of equations with
\int_{\partial \Omega} \lambda \left( \frac{\partial u}{\partial n} - 0 \right) \mathrm{d} \Omega=0?

Since I use a nodal Galerkin method, I cannot choose the basis functions to fulfil the boundary conditions, opposed to the modal Galerkin method.

Thank you