# Boundary conditions in conforming Galerkin method for biharmonic equation

I am trying to solve simple scalar biharmonic equation using bubnov-galerkin finite element method. I am using $H^2$ conforming basis functions. I was wondering that if anyone can give me some pointers on how can I further debug my algorithm.

I proceeded by converting the strong form of equation to weak form: \begin{align} \Delta^2 u &= f \\ u &= g \;\;\; \text{on} \;\;\; \Gamma_D \\ \Gamma_N &= \emptyset \end{align}

As an example problem, I am trying to solve: \begin{align} u &= \cos(4 \pi x) \cos(4 \pi y) \end{align} Boundary conditions are implemented using penalty method. The weak form is follows:

\begin{align} a(u,v) &= L(v) \\ a(u,v) &= \sum_{\Omega_K \in \mathcal T} \int_{\Omega_K} \Delta \psi_{i} \Delta \psi_{j} + \sum_{E \in \mathcal E} \gamma \int_{E} \psi_{i} \psi_{h,j} \\ L(v) &= \sum_{\Omega_K \in \mathcal T} \int_{\Omega_K} \psi_{i} f + \sum_{E \in \mathcal E} \gamma \int_{E} \psi_{i} g\\ \end{align}

Problem: When I solve this equation on a square domain, my L2 error barely converges (convergence rate ~ 0.3). However, if I solve $u = \sin(4 \pi x) \sin(4 \pi y)$, I get correct convergence rates. I have tried the following things to debug my code:

• I solves the poisson equation, so replaced the first integrand with the stiffness integrand. I get correct convergence rates. I concluded from this that my penalty method implementation is right

• I am using Cartesian grid, so Jacobian lines up with expectations

• Since, my shape functions are defined in parametric ($\xi,\eta$) space, I had to work out the transformation for the laplacian. I tested this transformation on polar coordinates.

Many thanks.

Apologies. Editted.

• Welcome to Scicomp.SE! You haven't actually written what's going wrong... – Christian Clason Apr 20 '15 at 10:26
• What's the domain you are solving on? From what you describe, I suspect it's a problem with the boundary conditions. Note that the biharmonic equation is a fourth-order equation and hence you need two boundary conditions. The natural boundary conditions are $\Delta u=0$ on $\Gamma_D$, so that's what your weak formulation imposes (and $u=\sin(4\pi x)\sin(4\pi y)$ presumably satisfies). Depending on the domain, this isn't valid for your other $u$ and so you get the wrong solution. In this case, you need to add the corresponding nonhomogeneous Neumann condition to your weak form. – Christian Clason Apr 20 '15 at 10:51
• Sorry, that should have been $\partial_n u = 0$. – Christian Clason Apr 20 '15 at 11:17
• I don't understand what you mean by "I am trying to solve" when you pre-define the solution, u. Is u a so-called "manufactured solution" where you substitute u into your PDE and Dirichlet BC to compute f and g and then use that f and g to solve the FE equations? – Bill Greene Apr 20 '15 at 17:45
• @ChristianClason: My domain is $[0 \times 1]^2$. ahh may be that's where I am going wrong. However, I thought that since I am solving dirichlet problem, I won't need to add the neumann boundary integrand? – user21674 Apr 21 '15 at 10:55

I agree with Christian Clason that your problem is likely with the boundary conditions. $\partial_n u$ is another essential boundary condition (like $u$) so, in general, you will also need to use your penalty approach to satisfy it. Apparently, the prescribed $\partial_n u$ for your $u=\sin(4\pi x)\sin(4\pi y)$ case is consistent with the zero natural boundary conditions that you are getting by default from the weak form.