What happens to the boundary conditions as a PDE is approximated by a lesser order PDE?

Consider a fourth order linear (biharmonic) PDE in two variables of the form

$\nabla^4u + c\nabla^2u-\lambda u = F(x,y)$; $(x,y) \in \Lambda$

To have uniqueness, we must specify two equations per point on $\partial \Lambda$. Now consider the limit where $c\rightarrow \infty$. The solution, $u$, is approximated by $\tilde u$, given by

$c\nabla^2\tilde u-\lambda \tilde u = F(x,y)$; $(x,y) \in \Lambda$

for which we must only specify one equation on the boundary.

The question is: what boundary conditions do we pick for $\tilde u$?

If I use a "clamped boundary" $u=u'=0$, the solution is approximated by using the boundary condition $\tilde u = 0$. If I use a "free boundary" $u''=u'''=0$, the appropriate approximation is $\tilde u'=0$.

-
Wouldn't it be more appropriate to consider the limit of $\epsilon\nabla^4u+\nabla^2u-\lambda u=F(x,y)$ as $\epsilon\to0$? The way you write it, as $c\to\infty$ I would expect to be left with just $\nabla^2\tilde u=0$. –  Harald Hanche-Olsen Jul 24 '10 at 8:10
also, does it even make sense to think of the limit? Since the derivative is generally an unbounded operator, I don't see how the limit can be justified to exist. (In other words, $\epsilon\nabla^4 u$ can still be potentially be much larger than $\nabla^2 u$ for arbitrarily small $\epsilon$, so one cannot really say that one can get an approximation by removing the term containing the "smallness parameter".) –  Willie Wong Jul 25 '10 at 23:29

The question seems to contain a misprint, as Harald Hanche-Olsen pointed out above. I am assuming that the real situation is as follows: $\epsilon\nabla^4u + c\nabla^2u-\lambda u=F(x,y)$, where $(x,y)\in\Gamma$, with a pair of boundary conditions, say $u=u'=0$, where $(x,y)\in\delta\Gamma$.
When $\epsilon$ is "small" it is natural to assume that solutions of this PDE should be "close" in some sense to solutions of the reduced PDE $c\nabla^2u'-\lambda u'=F(x,y)$. This is a typical singular perturbation problem (see Van Dyke's "Perturbation methods in fluid mechanics" or any other text on asymptotic expansions that contains a discussion of singular perturbation problems). It is clear that since the order of the original equation has been reduced, one can now only satisfy one boundary condition, not two of the original problem. This basically indicates the presence of a boundary layer. This boundary layer is a solution of another asymptotic limit of the problem, the one for which $\epsilon\nabla^4u$ can no longer be omitted (see the comment by Willie Wong). The problem can generally be non-trivial in 2D; let me just show the idea how this type of problems is solved in 1D.
Consider for $x\ge 0$ $$\epsilon^2u_{,xxxx}-c^2u_{,xx}-\lambda u=F(x) \quad\mbox{subject to}\quad u(0)=u_{,x}(0)=0$$ (I made some parameter changes to simplify the forthcoming algebra and changed the sign of $c^2$ as it actually makes more sense to pre-tension beam anyway). The "regular" solution $u_r$, as you pointed out, satisfies $c^2u_{r,xx}+\lambda u_r=-F(x)$. The boundary layer solution is typically extracted by "zooming" into the vicinity of the boundary by using a coordinate transformation $\xi=x/\epsilon^{\alpha}$ (so that $\partial/\partial x=\epsilon^{-\alpha}\partial/\partial\xi$). How do we find $\alpha$? Intuitively, we want to keep the previously vanishing term $\epsilon^2u_{,xxxx}$ in the equation and the non-trivial answer should balance it with at least one other term of the original equation. If $\alpha=1/2$ then first term is balanced by the third at $O(1)$, but the second term becomes $O(\epsilon^{-1})$. Hence, one needs to take $\alpha=1$ to balance the first two terms with the resulting equation $$\epsilon^{-2}u_{,\xi\xi\xi\xi}-c^2\epsilon^{-2}u_{,\xi\xi}-\lambda u=F(x).$$ Clearly, when $\epsilon$ is "small", the solution of this equation is approximated by the solution of boundary layer equation $u_{bl,\xi\xi}-c^2u_{bl}=0$. What happens next depends on the sign of $c^2$. I shall take $c^2$ to be strictly positive. Then the boundary layer solution is a rapidly decaying exponent $u_{bl}=U_{bl}\exp(-c\xi)=U_{bl}\exp(-cx/\epsilon)$.
Let us now assume that the general solution of our problem is given by sum $u=u_r+u_{bl}$ (this is so-called Vishik-Lyusternik method). This means that when you were posing boundary condition $u_{,x}(0)=0$, you effectively posed $$u_{r,x}(0)+u_{bl,x}(0)=0$$ Since $u_{bl}=U_{bl}\exp(-cx/\epsilon)$, then $u_{bl,x}=-cU_{bl}\exp(-cx/\epsilon)/\epsilon=-cu_{bl}/\epsilon$. Therefore, we can conclude that $$u_{r,x}(0)+u_{bl,x}(0)=u_{r,x}(0)-cu_{bl}(0)/\epsilon=0 \quad\mbox{so that}\quad u_{bl}(0)=\epsilon u_{r,x}(0)/c.$$ If we now turn our attention to the original first boundary condition, the same line of argument gives $$u(0)=u_r(0)+u_{bl}(0)=u_r(0)+\epsilon u_{r,x}(0)/c=0.$$
The conclusion then is: instead of solving the original fourth-order ODE with clamped boundary conditions, for sufficiently small $\epsilon$ one can simply solve the following second order boundary value problem $$c^2u_{r,xx}+\lambda u_r=-F(x) \quad\mbox{subject to the "effective" condition}\quad u_r(0)+\epsilon u_{r,x}(0)/c=0.$$ The resulting solution error will be of the order $O(\epsilon^2)$ (except in the close vicinity of the boundary where we made error $O(\epsilon)$). Your own suggested answer $u_r(0)=0$ is only accurate to within $O(\epsilon)$ everywhere. My explanation is quite loose here, but I hope it should give you the idea of how to tackle this kind of problems.