Computing The Bilaplacian

The bilaplacian problem

\begin{array}{rcl} -\Delta^2 u & =& f \quad\mbox{on } \Omega=(0,1)^2,\\ u= g,\; \frac{\partial u}{\partial n}&=& h \quad \mbox{on } \partial \Omega. \end{array}

Take $uu=\Delta u$ and use $\partial_n\Delta u =\Delta \partial_n u$ to work with the equations for $u$ and $uu$.

$$\begin{array}{lr} \begin{array}{rcll} -\Delta uu & =& f, & \mbox{on } \Omega=(0,1)^2,\\ \frac{\partial uu}{\partial n}&=&\Delta h &\mbox{on } \partial \Omega. \end{array} & \begin{array}{rcl} -\Delta u +uu& =&0 \quad\mbox{on } \Omega=(0,1)^2,\\ u= g,\; \frac{\partial u}{\partial n}&=& h \quad \mbox{on } \partial \Omega. \end{array} \end{array}$$

Taking $g=h=0$ we end up with the variational formulation

$$\begin{array}{rcl} \int_\Omega \nabla u \cdot\nabla vv+\int_\Omega uu\cdot vv & =& 0\\ \int_\Omega \nabla uu\cdot\nabla v & =& \int_\Omega f \cdot v \end{array}$$

where $u=v=0$ on $\partial \Omega$ but $uu$ and $vv$ are not fixed on $\partial \Omega$.

int n=100,nn=n+10;
real[int] xx(nn),yy(nn);

mesh Th=square(40,40);  // mesh definition of $\Omega$
fespace Vh(Th,P1);      // finite element space
macro laplacien(u,v) (dx(u)*dx(v)+dy(u)*dy(v)) // fin macro
real f=1;
Vh u,uu,v,vv;

solve bilap([u,uu],[v,vv],solver=LU,eps=1.0e-6) =
int2d(Th)(  laplacien(u,vv)+uu*vv //  - Delta u + uu =0 (vv)
+ laplacien(uu,v)  )     //  - Delta uu = 1 (v)

- int2d(Th)(f*v)
+ on(1,2,3,4,u=0); // => v=0 also on 1,2,3,4

plot(u,wait=1);

for (int i=0;i<=n;i++)
{
xx[i]=real(i)/n;
yy[i]=u(0.5,real(i)/n); // value of uh at point (0.5, i/10.)
}
plot([xx(0:n),yy(0:n)],wait=1);