I want to solve the biharmonic equation numerically, that is:

$$\Delta^2 u=f~~in~~\Omega$$ $$u=g_1~~on ~~\partial \Omega$$ $$\frac {\partial u}{\partial n}=g_2~~on ~~\partial \Omega$$

Using Green's formula, we have that if $u$ is smooth enough $$\int_{\partial \Omega}\frac {\partial \Delta u}{\partial n}v-\int_{\Omega}\nabla(\Delta u).\nabla v dx =\int_{\partial \Omega}\frac {\partial \Delta u}{\partial n}v-\int_{\partial \Omega}\frac {\partial v }{\partial n}\Delta u+\int_{\Omega}\Delta u \Delta v=\int_{\Omega}fv$$

the test space is $H^2(\Omega)$, How could I use the given boundary conditions$g_1,g_2$ .

how could I find the new terms $\frac {\partial v }{\partial n}$ and $\frac {\partial \Delta u}{\partial n}$

  • $\begingroup$ I would use Morley element here instead of the other suggestions. It's very easy to implement and the degrees of freedom directly correspond to the boundary conditions you want to give. $\endgroup$
    – knl
    Commented Oct 25, 2016 at 12:29

4 Answers 4


The biharmonic equation is the Euler-Lagrange equation of the Laplacian energy $\frac{1}{2} \langle \Delta u,\Delta u \rangle$. A systematic approach to discretize higher order problems is to convert the unconstrained problem to a constrained problem: Minimize $\frac{1}{2} \langle v,v\rangle$ s.t. $\Delta u=v$; that is, \begin{equation} \frac{1}{2} \langle v,v\rangle +\langle \lambda, \Delta u -v\rangle \\ \frac{1}{2} \langle v,v\rangle + \langle \lambda,\frac{\partial u}{\partial n}\rangle -\langle \lambda,v\rangle -\langle \nabla \lambda, \nabla v\rangle \end{equation} where we used Green's identity for the second equation and $\frac{\partial }{\partial n}$ is the normal derivative.

The new optimization problem can be discretized using piecewise linear elements as is commonly done for second order problems. This is the common thread in mixed finite-element discretizations.

  • $\begingroup$ But in this case where is $u=g_1$? $\endgroup$
    – rosa
    Commented Jun 11, 2014 at 4:15
  • $\begingroup$ @rosa: A Dirichlet boundary condition wouldn't be in the weak form. $\endgroup$ Commented Jun 11, 2014 at 4:47
  • $\begingroup$ To be more precise, a Dirichlet condition would enter as an additional constraint for the minimization. $\endgroup$ Commented Jun 11, 2014 at 6:58

I decided to expand my earlier comment into an answer. I'd suggest using the Morley element that uses $P_2$ basis in each element and the degrees of freedom are

  1. values at vertices
  2. normal derivatives at edge midpoints

It is probably the easiest element for biharmonic problem implementation wise.

Here is an example using sp.fem library (version 8019aa7):

from spfem.mesh import MeshTri                                                                                                                                                                                                           
from spfem.asm import AssemblerGlobal                                                                                                                                                                                                    
from spfem.element import ElementGlobalMorley                                                                                                                                                                                            
from spfem.utils import direct,const_cell                                                                                                                                                                                                
import numpy as np                                                                                                                                                                                                                       

# define mesh                                                                                                                                                                                                                            

# define element and assemble the stiffness matrix and load vector                                                                                                                                                                       

def bilinear(ddu,ddv):                                                                                                                                                                                                                   
    return ddu[0][0]*ddv[0][0]+\                                                                                                                                                                                                         

def F(x):                                                                                                                                                                                                                                
    return 50                                                                                                                                                                                                                            

f=a.iasm(lambda v,x: F(x)*v)                                                                                                                                                                                                             

# find boundary and interior DOFs                                                                                                                                                                                                        



# set boundary conditions                                                                                                                                                                                                                

def G1(x,y):                                                                                                                                                                                                                             
    return 0.05*np.sin(4*np.pi*y)                                                                                                                                                                                                        

def G2(x,y):                                                                                                                                                                                                                             
    return np.sin(4*np.pi*x)                                                                                                                                                                                                             

# check orientation of boundary normals from sign of jacobian determinant                                                                                                                                                                


# solve the system                                                                                                                                                                                                                       

# draw the solution                                                                                                                                                                                                                      

The solution looks like this:

enter image description here

Note that both $g_1$ and $g_2$ are easily defined as Dirichlet-type boundary conditions.

  • $\begingroup$ does this library allow for third order derivatives? $\endgroup$
    – uli.xu
    Commented Oct 26, 2016 at 21:18
  • $\begingroup$ I don't understand the line ddu[0][0]*ddv[0][0]+ddu[0][1]*ddv[0][1]+ddu[1][0]*ddv[1][0]+ddu[1][1]*ddv[1][1]. Is this computing $\Delta u \Delta v$ ? If yes, then I would have expect the code to read: ddu[0][0]*ddv[0][0]+ddu[0][0]*ddv[1][1]+ddu[1][1]*ddv[0][0]+ddu[1][1]*ddv[1][1]. I tried this, but then it seems the system is singular and direct crashes. What am I misunderstanding? $\endgroup$ Commented Oct 26, 2016 at 23:59
  • $\begingroup$ @uli.xu I support higher-order derivatives through performing projections into lower-order discontinuous finite element spaces. The usage is not trivial and probably cannot be used to implement e.g. $H^3$-conforming methods (haven't thought about it) but it works for evaluating functionals with arbitrary order derivatives. $\endgroup$
    – knl
    Commented Oct 27, 2016 at 11:33
  • $\begingroup$ @Alec Jacobson I answered your comment in GitHub. In short: the latter bilinear form is not elliptic. $\endgroup$
    – knl
    Commented Oct 27, 2016 at 11:35
  • $\begingroup$ @knl, is there an easier way to implement third order derivatives. I am trying to solve a Triharmonic PDE with an H^3-conforming method. Suggestions for any library especially where i can see how the third derivative is derived via the mapping? $\endgroup$
    – uli.xu
    Commented Oct 27, 2016 at 14:09

Using finite element methods, you also have the choice of discretizing the weak formulation directly, using two types of approaches:

  1. using H²-conforming elements like the Argyris element or the Hsieh-Clough-Tocher on triangles or the Bogner-Fox-Schmidt element on quadrilaterals.

  2. Using regular H¹-conforming elements in conjunction with the so called C0-interior penalty method.

If you take the weak formulation you derived through Green's formula, choosing test spaces with $v = \partial_n v = 0$ on the boundary makes the boundary terms go away.

Another approach very well suited to high order differential equations is isogeometric analysis.


The variational or weak formulation of the biharmonic problem reads as : find $u \in V_D$ such that $$ a(u,v) = \ell(v), \quad \forall v \in V_0,$$ where the bilinear and linear forms are given by $$a(u,v) = \int_\Omega \Delta u \Delta v \,d\Omega \quad \text{and} \quad \ell(v) = \int_\Omega f v\, d\Omega,$$ and the hyperplane and test spaces are given by $V_D :=\{ v \in H^2(\Omega), v = g_1, \frac{\partial v}{\partial n} = g_2 \quad on \quad \partial \Omega\}$ and $V_0 :=\{ v \in H^2(\Omega), v = 0, \frac{\partial v}{\partial n} = 0 \quad on \quad \partial \Omega\}.$

The higher order derivatives on the $\partial \Omega$ are actually the Neumann boundary conditions. These are not given hence your choice of test function is enough. This is a biharmonic problem with Dirichlet boundary conditions.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.