I would like to solve the biharmonic equation $\displaystyle \Delta^2 u=f(x,y)$ on $\displaystyle [-1,1]\times[-1,1]$ with $\displaystyle 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 $\displaystyle 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
$\displaystyle \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